pi.py

from graphics import GraphWin, Point, Circle, Rectangle, Text
from random import randrange
from math import pi, sqrt

height, width = 1044, 1044

# The module graphics.py can be obtained from, http://mcsp.wartburg.edu/zelle/python/graphics.py.
# It may also be installable via pip.

def intro():
    s = """
    Calculating π using a Monte Carlo integration method

    Source: https://en.wikipedia.org/wiki/Monte_Carlo_method#Introduction

    For example, consider a circle inscribed in a unit square. Given that the
    circle and the square have a ratio of areas that is π/4, the value of π can
    be approximated using a Monte Carlo method:

    Draw a square on the ground, then inscribe a circle within it.  Uniformly
    scatter some objects of uniform size (grains of rice or sand) over the
    square.  Count the number of objects inside the circle and the total number
    of objects.  The ratio of the two counts is an estimate of the ratio of the
    two areas, which is π/4.  Multiply the result by 4 to estimate π.

    In this procedure the domain of inputs is the square that circumscribes our
    circle. We generate random inputs by scattering grains over the square then
    perform a computation on each input (test whether it falls within the
    circle).  Finally, we aggregate the results to obtain our final result, the
    approximation of π.

    There are two important points to consider here: Firstly, if the grains are
    not uniformly distributed, then our approximation will be poor. Secondly,
    there should be a large  number of inputs. The approximation is generally
    poor if only a few grains are randomly dropped into the whole square. On
    average, the approximation improves as more grains are dropped.


    """

    print(s)


def prepare_window():
    global height, width

    # Set our window dimensions and display it.
    win = GraphWin("Monte Carlo Calculation of π", height, height)

    # Prepare our circle for drawing.
    center = Point(height/2, width/2)
    unit_circle = Circle(center,
                         (height - 20)/2 if height < width else (width - 20)/2)
    unit_circle.setFill("red1")
    unit_circle.setOutline("red1")

    p1, p2 = Point((height - 10), 10), Point(10, (width - 10))
    unit_square = Rectangle(p1, p2)
    unit_square.setFill("blue")
    unit_square.setOutline("blue")

    # The order of drawing matters, the square first, then the circle.
    # Otherwise, the circle may be obscured.
    unit_square.draw(win)
    unit_circle.draw(win)

    return unit_circle, unit_square, win


def get_input():
    num_points = int(float(input("How many points to draw? ")))
    return num_points


def in_circle(point, circ):

    # get the distance between pt1 and circ using the
    # distance formula
    dx = point.getX() - circ.getCenter().getX()
    dy = point.getY() - circ.getCenter().getY()
    dist = sqrt(dx*dx + dy*dy)

    # check whether the distance is less than the radius
    return dist <= circ.getRadius()


def calc_pi(points, circle):
    global height, width
    total = len(points)
    total_in_circle = 0

    for p in points:
        if in_circle(p, circle):
            total_in_circle += 1

    return (total_in_circle/total) * 4


def scatter(win, circle, n):
    global height, width
    points = []

    # Create Text object to show us current stats.
    status = Text(Point(80, 50), "")
    status.setSize(16)
    status.draw(win)

    for i in range(n):
        # Random position in side square
        x = randrange(10, width - 10 + 1)
        y = randrange(10, height - 10 + 1)

        # Create and save point in a list
        point = Point(x, y)
        points.append(point)

        # Set the points properties, including drawing.
        point.setFill("white")
        point.setOutline("white")
        point.draw(win)

        # Update the stats text label with new info.
        pi_sim = calc_pi(points, circle)
        status.setText("Points drawn: {0}\n\
                       π     = {1:0.6f}\n\
                       π_sim = {2:0.6f}\n\
                       diff = {3:0.6f}%".format(len(points),
                                                pi,
                                                pi_sim,
                                                (1 - (pi_sim/pi))*100))


def main():
    # Print an introduction
    intro()

    # Get input from user on how many points to draw.
    n = get_input()

    # Prepare our window for simulation
    u_circle, u_square, window = prepare_window()

    # Scatter points all over the unit square.
    scatter(window, u_circle, n)

    # Get mouse click to exit (useful when run from a command line terminal)
    window.getMouse()

main()