"""A module to estimate the numerical value of pi using Monte Carlo"""

import numpy  # to generate arrays of random numbers
import pylab  # to make plots

def throwDarts(n):
    """returns an array of n uniformly random (x,y) pairs lying within the square that circumscribes the unit circle centered at the origin, i.e., the square with corners at (-1,-1), (-1,1), (1,1), (1,-1)"""
    darts = 2*numpy.random.random((n,2)) - 1
    return darts

def inUnitCircle(p):
    """return a boolean array, whose elements are True if the corresponding point in the array p is within the unit circle centered at the origin, and False otherwise -- hint: use numpy.linalg.norm to find the length of a vector"""
    #return numpy.linalg.norm(p)<=1.0
    return numpy.linalg.norm(p,axis=-1)<=1.0

def EstimatePi1(n):
    """returns an estimate of pi by drawing n random numbers in the square [[-1,1], [-1,1]] and calculating what fraction land within the unit circle"""
    number_in_circle = 0
    for i in range(n):
        dart = throwDarts(1)
        if inUnitCircle(dart):
            number_in_circle += 1
    return 4*number_in_circle/n

def EstimatePi2(n):
    """returns an estimate of pi by drawing n random numbers in the square [[-1,1], [-1,1]] and calculating what fraction land within the unit circle; in this version, draw all the random numbers at once and do a single array comparison to count what fraction lie within the unit circle"""
    darts = throwDarts(n)
    number_in_circle = numpy.sum(inUnitCircle(darts))
    return 4*number_in_circle/n

def EstimatePi3(n, block=10000):
    """returns an estimate of pi by drawing n random numbers in the square [[-1,1], [-1,1]] and calculating what fraction land within the unit circle; in this version, draw random numbers in blocks of the specified size, and keep a running total of the number of points within the unit circle"""
    total_number = 0
    i = 0
    while i<n:
        if n-i < block:
            block = n-i
        darts = throwDarts(block)
        number_in_circle = numpy.sum(inUnitCircle(darts))
        total_number += number_in_circle
        i += block
    return 4*total_number/n

def RunningEstimatePi1(n, block=10000):
    """returns a running estimate of pi in blocks by drawing n random numbers in the square [[-1,1], [-1,1]] and calculating what fraction land within the unit circle"""
    total_number = 0
    ntries = []
    estimates = []
    i = 0
    while i<n:
        if n-i < block:
            block = n-i
        darts = throwDarts(block)
        number_in_circle = numpy.sum(inUnitCircle(darts))
        total_number += number_in_circle
        i += block
        ntries.append(i)
        estimates.append(4*total_number/i)
    return numpy.array(ntries), numpy.array(estimates)

def PlotDarts(n):
    """Draw a scatter plot of n thrown darts, coloring those inside and outside the unit circle differently"""
    darts = throwDarts(n)
    inside = numpy.take(darts, numpy.where(inUnitCircle(darts))[0], axis=0)
    outside = numpy.take(darts, numpy.where(-inUnitCircle(darts))[0], axis=0)
    coll = pylab.scatter(inside[:,0], inside[:,1], c='r')
    pylab.scatter(outside[:,0], outside[:,1], c='b')
    coll.axes.set_aspect('equal')
    pylab.xlim(-1,1)
    pylab.ylim(-1,1)