View source for EGR 103/Fall 2018/Lec 12

This page will go through the development of the program we wrote in class.  Each step along the way, I will describe what functionality we wanted to add and you will be able to see the code that accomplished that goal.

== Introduction ==
We will be writing a program that uses a [https://en.wikipedia.org/wiki/Monte_Carlo_method Monte Carlo method] to estimate <math>\pi</math>.   The program will simulate throwing some number of darts at a dart board of radius <math>r</math> mounted on a square back board that is <math>2r</math> on a side.  If the <math>x</math> and </math>y</math> coordinates are each uniformly distributed between <math>-r</math> and <math>r</math>, then the probability that a dart hits the dart board is the same as the ratio of the area of the dart board to the area of the back board - that is, <math>p_{hit}=\frac{\pi r^2}{(2r)^2}=\frac{\pi}{4}</math>.  We will thus estimate <math>\pi</math> by dividing the number of darts that hit the dart board by the total number of darts thrown and then multiplying that by 4.

== Detached figure windows ==
As a part of this code, we will be watching a plot as it evolves.  That is not something you want to do with inline graphics.  To make Spyder put figures in their own windows:
* For Windows, go to Tools->Preferences and with OSX, go to python->Preferences
* In the Preferences window, select IPython console from the list at left
* In the right window, select the Graphics tab
* In the Graphics tab, in the Graphics backend section, choose Automatic.  
* Click apply, then get out of the Preferences window.  You will need to restart Spyder for this to take effect.

== Tasks ==
We broke the tasks up into several parts and wrote to code to accomplish those tasks.  The main tasks were:
* Import modules
* Draw the board
* Throw a dart
* Throw lots of darts
* Plot the dart locations
* Determine how many darts land on the dart board
* Print how many darts land on the dart board and the estimates of <math>\pi</math>
* Refine how that is printed
* Plot the evolution of the estimates
* Have the value of a variable decide if we should plot while the simulation is running or only at the end
* Figure out what's next

=== Importing modules ===
As always, start by importing the necessary modules.  We will be doing math, using computational methods, and making graphs, so we will need:
<syntaxhighlight  lang='python' line highlight='1-3'>
import math as m
import numpy as np
import matplotlib.pyplot as plt
</syntaxhighlight>

=== Drawing the board ===
Assuming a dart board of radius <math>r</math>, the square that encloses it would be <math>2r</math> on a side.  Since we will later be using the distance from the center of the circle to determine if a dart is on the dart board, it is mathematically easiest to have the middle of the board at (0, 0).  Drawing the back board thus means drawing a square with vertices at (-r, -r), (r, -r), (r, r), and (-r, r).  The the circle, we will need a list of x and y coordinates that defines a circle; the easiest way to do that is to use polar coordinates with a constant radius and an angle that goes from 0 to <math>2\pi</math>.  We need to know how many values to use in making the circle, so that will be another input to the function.  We chose a default radius of 1 and a default number of points of 100:
<syntaxhighlight  lang='python' line highlight='5-10'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
</syntaxhighlight>

Test the code by calling it from the script:
<syntaxhighlight  lang='python' line highlight='12-15'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')

if __name__ == '__main__':
    plt.figure(1)
    plt.clf()
    draw_board()
</syntaxhighlight>

The problem with this is that Python will use all the figure window and depending on the shape of the window, the square and the circle might look non-square and non-circular.  To require that Python plot using the same scale along both axes, use the .axis('equal') method of plt:

<syntaxhighlight  lang='python' line highlight='11'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')

if __name__ == '__main__':
    plt.figure(1)
    plt.clf()
    draw_board()
</syntaxhighlight>

You will now have a square back board and a circular dart board.  As an aside, if you want to see something trippy, copy the following into the console after you have run your script:
<syntaxhighlight  lang='python'>
r = 1
for k in range(10):
     draw_board(r)
     r = r*m.sqrt(2)
</syntaxhighlight>

=== Throw a dart ===
Throwing a dart means figuring out the <math>x</math> and <math>y</math> coordinates for a dart.  To determine the valid location of darts, we need to know the radius of the dart board.  Given all that, adding a function to return the location of a thrown dart might look like:
<syntaxhighlight  lang='python' line highlight='13-16'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)

if __name__ == '__main__':
    plt.figure(1)
    plt.clf()
    draw_board()
</syntaxhighlight>
You can test this in the console by throwing darts at a board of a given size (say, 10):

<syntaxhighlight  lang='python'>
for k in range(10):
    print(throw_darts(10))
</syntaxhighlight>
Make sure both coordinates are in the proper range.

=== Throw lots of darts ===
Now that you can throw a single dart, you might want to throw several darts.  You should add variables to keep track of the total number of darts to throw and of the total number thrown: 
<syntaxhighlight  lang='python' line highlight='19-20, 26-28'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    
    plt.figure(1)
    plt.clf()
    draw_board()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts()
</syntaxhighlight>

=== Plot the dart locations ===
Now that you are throwing darts, you can plot the locations at which they land.  There are two fundamentally different ways to do this - you can plot each dart as it is thrown or you can wait until the end to plot all the darts.  We are going to first look at how to plot each dart and watch the board fill up with darts.  They key here will be a function that takes the dart coordinates and then updates the figure before moving on to the next command.  Having the figure update right now requires drawing the plot and then adding a pause to let the plot...plot:
<syntaxhighlight  lang='python' line highlight='18-22, 35'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    
    plt.figure(1)
    plt.clf()
    draw_board()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts()
        plot_darts(x_d, y_d)
</syntaxhighlight>

=== Determine how many darts land on the dart board ===
Now that you are throwing darts and see where they hit, you can add code to keep track of how many are landing on the dart board.  To do this, you will need to let your code know how big the dart board is and then calculate how far the throw is from the center.  You need to keep track of how many darts hit the dartboard.
<syntaxhighlight  lang='python' line highlight='27-28, 32, 36, 38-40'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    inside = 0

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        plot_darts(x_d, y_d)
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1 
</syntaxhighlight>
Note the changes on lines 32 and 36 to allow for board sizes other than 1.

=== Print how many darts land on the dart board and the estimates of <math>\pi</math> ===
While you can see the value of <code>inside</code> in the variable explorer, it would make more sense to print information out as the simulation is running.  This might include how many darts have been thrown, how many hit the dart board, and what the estimate for <math>\pi</math> is:
<syntaxhighlight  lang='python' line highlight='42'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    inside = 0

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        plot_darts(x_d, y_d)
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1

        print(darts_thrown, inside, 4*inside/darts_thrown)
</syntaxhighlight>

=== Refine how that is printed ===
That table is fine for a rough idea of how things are going, but it would be better to make a formatted table.  Assuming you throw fewer than 100,000 darts, you can format the number of darts thrown and the number of hits to be integers with 5 spaces reserved and the estimate of pi as a floating point number using scientific notation with three digits after the decimal point.  Note the spaces in the format strings - those make sure the three numbers are separated.
<syntaxhighlight  lang='python' line highlight='42-44'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    inside = 0

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        plot_darts(x_d, y_d)
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1

        print('{:5d}'.format(darts_thrown) + 
               ' {:5d}'.format(inside) +
               ' {:0.3e}'.format(4*inside/darts_thrown))
</syntaxhighlight>
=== Plot the evolution of the estimates ===
It may be interesting to see how the estimate for <math>\pi</math> changes as the simulation runs.  While the table now prints this information, making a graph might be useful.  We can once again either make a graph by adding each estimate or by graphing all the estimates at the end.  We will do the former here by creating a second figure window and plotting the estimates on it:
<syntaxhighlight  lang='python' line highlight='24-28, 40-41, 54'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)

def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    inside = 0

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    plt.figure(2)
    plt.clf()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        plot_darts(x_d, y_d)
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1

        print('{:5d}'.format(darts_thrown) + 
               ' {:5d}'.format(inside) +
               ' {:0.3e}'.format(4*inside/darts_thrown))
        plot_est(darts_thrown, 4*inside/darts_thrown)
</syntaxhighlight>
While this works, we are having to calculate the estimates twice.  It might make more sense to store all the estimates and then give the most recently calculated estimate to both the line that prints the estimate and the one that plots it:
<syntaxhighlight  lang='python' line highlight='24-28, 40-41, 54'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)

def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    inside = 0
    pi_est = []

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    plt.figure(2)
    plt.clf()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        plot_darts(x_d, y_d)
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1

        pi_est += [4*inside/darts_thrown]
        print('{:5d}'.format(darts_thrown) + 
               ' {:5d}'.format(inside) +
               ' {:0.3e}'.format(pi_est[-1]))
        plot_est(darts_thrown, pi_est[-1])
</syntaxhighlight>

=== Have the value of a variable decide if we should plot while the simulation is running or only at the end ===
Watching the animation is great, but very slow.  If I want to calculate the estimate after 50,000 trials, that will take a great deal of time.  Instead, I can add a variable that will indicate to the program whether it should plot as the simulation is running or only at the end.  For the latter, we will need to keep track of all the dart throws - not just the current one.  The do_plot variable will be the determining factor of whether I plot in the while loop or only when it is done.  Pay careful attention to the indentations and also the fact that the plot commands have been moved to the end of the loop:
<syntaxhighlight  lang='python' line highlight='34, 37-38, 50-51, 61-63, 65-67'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y):
    plt.figure(1)
    plt.plot(x, y, 'b.')
    plt.draw()
    plt.pause(0.01)

def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 100
    darts_thrown = 0
    board_rad = 1
    do_plot = 1
    inside = 0
    pi_est = []
    x_vec = []
    y_vec = []

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    plt.figure(2)
    plt.clf()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        x_vec += [x_d]
        y_vec += [y_d]
        r_d = m.sqrt(x_d**2 + y_d**2)
        if r_d < board_rad:
            inside += 1

        pi_est += [4*inside/darts_thrown]
        print('{:5d}'.format(darts_thrown) + 
               ' {:5d}'.format(inside) +
               ' {:0.3e}'.format(pi_est[-1]))
        
        if do_plot:
            plot_darts(x_d, y_d)
            plot_est(darts_thrown, pi_est[-1])

    if not do_plot:
        plot_darts(x_vec, y_vec)
        plot_est(np.linspace(1, darts_thrown, darts_thrown), pi_est)
</syntaxhighlight>

=== Figure out what's next ===
The program will nor run through a pre-determined number of simulations, printing a table of estimates as it goes and plotting the evolution of the dart board and the estimates either during or after the simulation.  The question is - what's next?  There were several recommendations during class, including:
* Re-print the dart board at the end so you can see it even when there are thousands of blue dots - this requires adding:
:<syntaxhighlight  lang='python' highlight='1-2'>
    plt.figure(1)
    draw_board(board_rad)
</syntaxhighlight>
to the end of the code
* Have the actual value of <math>\pi</math> show up on the estimate graph - this requires changing the <code>plot_est</code> function to:
:<syntaxhighlight  lang='python' highlight='4'>
def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.plot(nums, m.pi*np.ones(np.size(nums)), 'g.')
    plt.draw()
    plt.pause(0.01)
</syntaxhighlight>
: Note that the <code>np.ones(np.size(nums))</code> is there to make sure there are as many y values as there are x values in case nums is a list.
* Show the darts in different colors - several ways to do this but you have to carefully consider that the <code>plot_darts</code> command is currently set to work with both individual points and lists of points.  It would be more complicated to do this with the whole list of points since there is not currently a mechanism to decide which are in and which are out.  Some options include:
** Only choosing different colors when plotting one point at a time and plotting the whole list in a single color
** Keeping track of which points are inside and which are outside and plotting each list separately - here is an example of the code with changes needed to make that happen (and also re-print the board and show the actual value of pi on the estimate evolution):
<syntaxhighlight  lang='python' line highlight='18,20,27,38-41,54,57-62,70,74-75,78-79'>
import math as m
import numpy as np
import matplotlib.pyplot as plt

def draw_board(r=1, pts=100):
    plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
    theta = np.linspace(0, 2*m.pi, pts)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    plt.plot(x, y, 'r-')
    plt.axis('equal')
    
def throw_darts(r=1):
    x = np.random.uniform(-r, r)
    y = np.random.uniform(-r, r)
    return (x, y)
    
def plot_darts(x, y, color='b'):
    plt.figure(1)
    plt.plot(x, y, '.', mfc=color, mec=color)
    plt.draw()
    plt.pause(0.01)

def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.plot(nums, m.pi*np.ones(len(nums)), 'g.')
    plt.draw()
    plt.pause(0.01)
    
if __name__ == '__main__':
    darts_to_throw = 10000
    darts_thrown = 0
    board_rad = 1
    do_plot = 0
    inside = 0
    pi_est = []
    x_i_vec = []
    y_i_vec = []
    x_o_vec = []
    y_o_vec = []

    plt.figure(1)
    plt.clf()
    draw_board(r=board_rad)

    plt.figure(2)
    plt.clf()

    while darts_thrown < darts_to_throw:
        darts_thrown += 1
        x_d, y_d = throw_darts(board_rad)
        r_d = m.sqrt(x_d**2 + y_d**2)
        color='m'
        if r_d < board_rad:
            inside += 1
            x_i_vec += [x_d]
            y_i_vec += [y_d]
            color='g'
        else:
            x_o_vec += [x_d]
            y_o_vec += [y_d]

        pi_est += [4*inside/darts_thrown]
        print('{:5d}'.format(darts_thrown) + 
               ' {:5d}'.format(inside) +
               ' {:0.3e}'.format(pi_est[-1]))
        
        if do_plot:
            plot_darts(x_d, y_d, color)
            plot_est(darts_thrown, pi_est[-1])

    if not do_plot:
        plot_darts(x_i_vec, y_i_vec, 'g')
        plot_darts(x_o_vec, y_o_vec, 'm')
        plot_est(np.linspace(1, darts_thrown, darts_thrown), pi_est)

    plt.figure(1)
    draw_board(board_rad)
</syntaxhighlight>
* Finally, you may only want to print the table every once in a while; you could set a variable called <code>print_every</code> equal to some value and then add change the table printing code to something like:
<syntaxhighlight  lang='python' >
        if (darts_thrown in [1, darts_to_throw] or 
            darts_thrown % print_every == 0):
            print('{:5d}'.format(darts_thrown) + 
                   ' {:5d}'.format(inside) +
                   ' {:0.3e}'.format(pi_est[-1]))
</syntaxhighlight>