import random
import math

def stdDev(values):
    mean = float(sum(values))/len(values)
    diffsum = 0.0
    for item in values:
        diff = item - mean
        diffsum = diffsum + diff*diff
    return math.sqrt(diffsum/len(values))
 
def estimatePi(numRandomPts):
    inCircleCount = 0
    pointNumber = 0
    while pointNumber < numRandomPts:
        x = random.random()
        y = random.random()
        if x * x + y * y <= 1.0:
            inCircleCount = inCircleCount + 1
        pointNumber = pointNumber + 1
        
    return 4 * inCircleCount/float(numRandomPts)

def doPiTrials(numRandomPts, numTrials):
    estimates = []
    for trial in range(numTrials):
        piEst = estimatePi(numRandomPts)
        estimates.append(piEst)
    curEst = sum(estimates)/len(estimates)
    sDev = stdDev(estimates)
    print("Estimate: {}, SD: {}, num random pts: {}".format(curEst, sDev, numRandomPts))
    return (numRandomPts, curEst, sDev)

# Execute doPiTrials for different numbers of random points (and given number of trials),
# starting with minPoints and increasing by factor until number exceeds maxPoints.
# Return list of results returned by the calls to doPiTrials.
# Thus, results will be of form [(minPoints, ...std...), (factor*minPoints, ...std..), ...]
#
def findPi(minPoints, maxPoints, factor, numTrials=5):
    numPoints = minPoints
    results = []
    while numPoints <= maxPoints:
        results.append(doPiTrials(numPoints, numTrials ))
        numPoints = numPoints * factor
    return results

# Uncomment section below if you have pylab


import pylab

# results = findPi(10, 100000000, 10)
# plotResults(results)
# pylab.show()
#
#
# Create figure with two subplots, one showing avg. pi estimate for different numbers of
# random points, the other showing standard deviation of estimates.
#
def plotResults(results):
    numPoints, estimates, SDs = tuple(zip(*results))

    pylab.figure(1)
    pylab.clf()
    pylab.subplot(121)
    pylab.title("Pi estimates")
    pylab.xlabel("Number of random points")
    pylab.ylabel("Estimated pi value")
    pylab.semilogx()
    pylab.plot(numPoints, estimates, 'ro-') 
    
    pylab.subplot(122)
    pylab.title("Std dev. of pi estimates")
    pylab.xlabel("Number of random points")
    pylab.ylabel("Standard deviation")
    pylab.semilogx()
    pylab.plot(numPoints, SDs, 'ro-')