Modeling and Simulation in Python

Data604 Hw 2 Farhana Zahir

Chapter 4

Copyright 2017 Allen Downey

License: Creative Commons Attribution 4.0 International

In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='all'

# import functions from the modsim library
from modsim import *
import numpy as np

Returning values

Here's a simple function that returns a value:

In [2]:
def add_five(x):
    return x + 5

And here's how we call it.

In [3]:
y = add_five(3)

If you run a function on the last line of a cell, Jupyter displays the result:

In [4]:
add_five(5)
Out[4]:
10

But that can be a bad habit, because usually if you call a function and don't assign the result in a variable, the result gets discarded.

In the following example, Jupyter shows the second result, but the first result just disappears.

In [5]:
add_five(3)
add_five(5)
Out[5]:
8
Out[5]:
10

When you call a function that returns a variable, it is generally a good idea to assign the result to a variable.

In [6]:
y1 = add_five(3)
y2 = add_five(5)

print(y1, y2)

8 10

Exercise: Write a function called make_state that creates a State object with the state variables olin=10 and wellesley=2, and then returns the new State object.

Write a line of code that calls make_state and assigns the result to a variable named init.

In [7]:
#Function to make state object
def make_state():
    state = State(olin=10, wellesley=2) #object hasvariables olin and wellesley
    return state    #returns the object
In [8]:
init = make_state() #init variable calls the function make_state

Running simulations

Here's the code from the previous notebook.

In [32]:
def step(state, p1, p2):
    """Simulate one minute of time.
    
    state: bikeshare State object
    p1: probability of an Olin->Wellesley customer arrival
    p2: probability of a Wellesley->Olin customer arrival
    """
    if flip(p1):
        bike_to_wellesley(state)
    
    if flip(p2):
        bike_to_olin(state)
        
def bike_to_wellesley(state):
    """Move one bike from Olin to Wellesley.
    
    state: bikeshare State object
    """
    if state.olin == 0:
        state.olin_empty += 1
        return
    state.olin -= 1
    state.wellesley += 1
    
def bike_to_olin(state):
    """Move one bike from Wellesley to Olin.
    
    state: bikeshare State object
    """
    if state.wellesley == 0:
        state.wellesley_empty += 1
        return
    state.wellesley -= 1
    state.olin += 1
    
def decorate_bikeshare():
    """Add a title and label the axes."""
    decorate(title='Olin-Wellesley Bikeshare',
             xlabel='Time step (min)', 
             ylabel='Number of bikes')

Here's a modified version of run_simulation that creates a State object, runs the simulation, and returns the State object.

In [33]:
def run_simulation(p1, p2, num_steps):
    """Simulate the given number of time steps.
    
    p1: probability of an Olin->Wellesley customer arrival
    p2: probability of a Wellesley->Olin customer arrival
    num_steps: number of time steps
    """
    state = State(olin=10, wellesley=2, 
                  olin_empty=0, wellesley_empty=0)
                    
    for i in range(num_steps):
        step(state, p1, p2)
        
    return state

Now run_simulation doesn't plot anything:

In [34]:
state = run_simulation(0.4, 0.2, 60)

But after the simulation, we can read the metrics from the State object.

In [35]:
state.olin_empty
Out[35]:
4

Now we can run simulations with different values for the parameters. When p1 is small, we probably don't run out of bikes at Olin.

In [36]:
state = run_simulation(0.2, 0.2, 60)
state.olin_empty
Out[36]:
0

When p1 is large, we probably do.

In [37]:
state = run_simulation(0.6, 0.2, 60)
state.olin_empty
Out[37]:
13

More for loops

linspace creates a NumPy array of equally spaced numbers.

In [38]:
p1_array = linspace(0, 1, 5)

We can use an array in a for loop, like this:

In [39]:
for p1 in p1_array:
    print(p1)

0.0
0.25
0.5
0.75
1.0

This will come in handy in the next section.

linspace is defined in modsim.py. You can get the documentation using help.

In [17]:
help(linspace)

Help on function linspace in module modsim.modsim:

linspace(start, stop, num=50, **options)
    Returns an array of evenly-spaced values in the interval [start, stop].
    
    start: first value
    stop: last value
    num: number of values
    
    Also accepts the same keyword arguments as np.linspace.  See
    https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html
    
    returns: array or Quantity

linspace is based on a NumPy function with the same name. Click here to read more about how to use it.

Exercise: Use linspace to make an array of 10 equally spaced numbers from 1 to 10 (including both).

In [40]:
np.linspace(1, 10, num=10, endpoint=True) #used the numpy function
Out[40]:
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.])

Exercise: The modsim library provides a related function called linrange. You can view the documentation by running the following cell:

In [19]:
help(linrange)

Help on function linrange in module modsim.modsim:

linrange(start=0, stop=None, step=1, endpoint=False, **options)
    Returns an array of evenly-spaced values in an interval.
    
    By default, the last value in the array is `stop-step`
    (at least approximately).
    If you provide the keyword argument `endpoint=True`,
    the last value in the array is `stop`.
    
    This function works best if the space between start and stop
    is divisible by step; otherwise the results might be surprising.
    
    start: first value
    stop: last value
    step: space between values
    
    returns: NumPy array

Use linrange to make an array of numbers from 1 to 11 with a step size of 2.

In [41]:
linrange(1, 11, 2, endpoint=True) #include endpoint 11 in the array with endpoint=True
Out[41]:
array([ 1,  3,  5,  7,  9, 11], dtype=int32)

Sweeping parameters

p1_array contains a range of values for p1.

In [42]:
p2 = 0.2
num_steps = 60
p1_array = linspace(0, 1, 11)

The following loop runs a simulation for each value of p1 in p1_array; after each simulation, it prints the number of unhappy customers at the Olin station:

In [43]:
for p1 in p1_array:
    state = run_simulation(p1, p2, num_steps)
    print(p1, state.olin_empty)

0.0 0
0.1 0
0.2 0
0.30000000000000004 0
0.4 6
0.5 13
0.6000000000000001 17
0.7000000000000001 29
0.8 29
0.9 29
1.0 28

Now we can do the same thing, but storing the results in a SweepSeries instead of printing them.

In [44]:
sweep = SweepSeries()

for p1 in p1_array:
    state = run_simulation(p1, p2, num_steps)
    sweep[p1] = state.olin_empty

And then we can plot the results.

In [45]:
plot(sweep, label='Olin')

decorate(title='Olin-Wellesley Bikeshare',
         xlabel='Arrival rate at Olin (p1 in customers/min)', 
         ylabel='Number of unhappy customers')
Out[45]:
[<matplotlib.lines.Line2D at 0x163f5956280>]

Exercises

Exercise: Wrap this code in a function named sweep_p1 that takes an array called p1_array as a parameter. It should create a new SweepSeries, run a simulation for each value of p1 in p1_array, store the results in the SweepSeries, and return the SweepSeries.

Use your function to plot the number of unhappy customers at Olin as a function of p1. Label the axes.

In [46]:
def sweep_p1(p1_array): #takes parameter p1_array
    p2 = 0.2 #p2 kept constant
    num_steps = 60 #iterations=60
    sweep = SweepSeries() #stores in array sweepseries
    
    for p1 in p1_array:
        state = run_simulation(p1, p2, num_steps)
        sweep[p1] = state.olin_empty
        
    return sweep
In [26]:
p1_array = linspace(0, 1, 11) 
np.random.seed(15061) #using seed to get the same results 
sweep = sweep_p1(p1_array)
plot(sweep, '-b', label='Olin') #'-b':# plot x and y using blue line
decorate(title='Olin-Wellesley Bikeshare',
         xlabel='Arrival rate at Olin (p1 in customers/min)',  #Add title, x axis label, y axis label
         ylabel='Number of unhappy customers')
Out[26]:
[<matplotlib.lines.Line2D at 0x163f478ad30>]

Exercise: Write a function called sweep_p2 that runs simulations with p1=0.5 and a range of values for p2. It should store the results in a SweepSeries and return the SweepSeries.

In [47]:
def sweep_p2(p2_array): #This function keeps p1 constant but takes an array of values for p2
    p1 = 0.5
    num_steps = 60
    sweep = SweepSeries() #values stored in sweep
    
    for p2 in p2_array:
        state = run_simulation(p1, p2, num_steps)
        sweep[p2] = state.olin_empty
        
    return sweep
In [28]:
p2_array = linspace(0, 1, 11)
np.random.seed(15062)
sweep = sweep_p2(p2_array)
plot(sweep, '-b', label='Olin')

decorate(title='Olin-Wellesley Bikeshare',
         xlabel='Arrival rate at Wellesley (p2 in customers/min)', 
         ylabel='Number of unhappy customers')
Out[28]:
[<matplotlib.lines.Line2D at 0x163f47fd5b0>]

Optional Exercises

The following two exercises are a little more challenging. If you are comfortable with what you have learned so far, you should give them a try. If you feel like you have your hands full, you might want to skip them for now.

Exercise: Because our simulations are random, the results vary from one run to another, and the results of a parameter sweep tend to be noisy. We can get a clearer picture of the relationship between a parameter and a metric by running multiple simulations with the same parameter and taking the average of the results.

Write a function called run_multiple_simulations that takes as parameters p1, p2, num_steps, and num_runs.

num_runs specifies how many times it should call run_simulation.

After each run, it should store the total number of unhappy customers (at Olin or Wellesley) in a TimeSeries. At the end, it should return the TimeSeries.

Test your function with parameters

p1 = 0.3
p2 = 0.3
num_steps = 60
num_runs = 10

Display the resulting TimeSeries and use the mean function provided by the TimeSeries object to compute the average number of unhappy customers (see Section 2.7).

In [48]:
def run_multiple_simulations(p1, p2, num_steps, num_runs):#takes as parameter no of runs for simulation
    results = TimeSeries()
    
    for i in range(num_runs):
        state = run_simulation(p1, p2, num_steps) #runs the simulation for num_runs
        results[i] = state.olin_empty + state.wellesley_empty #total unhappy customers for both
       
    return results
In [53]:
p1 = 0.3
p2 = 0.3
num_steps = 60
num_runs = 10
np.random.seed(15063)
run_multiple_simulations(p1, p2, num_steps, num_runs)
Out[53]:
values
0 0
1 4
2 1
3 0
4 0
5 1
6 0
7 0
8 1
9 0
In [54]:
plot(results, label='total', color='green')
    
decorate(title='Olin-Wellesley Bikeshare',
         xlabel='Arrival rate at Olin (p1 in customers/min)', 
         ylabel='Average total unhappy customers')
Out[54]:
[<matplotlib.lines.Line2D at 0x163f59ead60>]

Exercise: Continuting the previous exercise, use run_multiple_simulations to run simulations with a range of values for p1 and

p2 = 0.3
num_steps = 60
num_runs = 20

Store the results in a SweepSeries, then plot the average number of unhappy customers as a function of p1. Label the axes.

What value of p1 minimizes the average number of unhappy customers?

In [55]:
p1_array = linspace(0, 1, 11)
p2 = 0.3
num_steps = 60
num_runs = 20

sweep = SweepSeries()
for p1 in p1_array:
    results = run_multiple_simulations(p1, p2, num_steps, num_runs)
    sweep[p1] = results.mean()
In [60]:
plot(sweep, label='total', color='green')
    
decorate(title='Olin-Wellesley Bikeshare',
         xlabel='Arrival rate at Olin (p1 in customers/min)', 
         ylabel='Average total unhappy customers')

print('Value of p1 that minimizes the average number of unhappy customers:', 0.4)
Out[60]:
[<matplotlib.lines.Line2D at 0x163f5b6ca30>]
Value of p1 that minimizes the average number of unhappy customers: 0.4