import numpy as np
import matplotlib.pyplot as plt
from matplotlib import interactive
interactive(True)

'''
  Pi by MteCarlo vs. various integration schemes of circle's area
  with the same number of regularly spaced points as those 
  randomly thrown onto a bounding square.
'''


def MteCarlo_4Int(N):
    N = N//2*2 +1
# in N experiments, sum three dice 
    xr = np.random.rand(N)      # rnd float array
    yr = np.random.rand(N)      # rnd float array
# plot a quarter circle  
    xx = np.linspace(0,1,N)
    yy = np.sqrt(1-xx*xx)
# compute pi by trapezoid rule (2nd order) & higher order methods
    dx = xx[1]-xx[0]
    print(' dx=',dx, ' x[-1]',xx[-1])
    subtr = yy[0]/2 + yy[-1]/2
    integ1_pi =  yy.sum() *dx*4 
    integ2_pi = (yy.sum()-subtr) *dx*4 
    integ3_pi = integ4_pi = 0.
    # Quadratic & cubic Simpson's rule
    for i in range(N):
        coeff2 = 2.*(1+i%2)     # 2,4,2,4,2,...
        coeff3 = 2+min(i%3,1)   # 2,3,3,2,3,3,... 
        integ3_pi += coeff2*yy[i]
        integ4_pi += coeff3*yy[i]
    integ3_pi = (integ3_pi - yy[0]+yy[-1]) * dx*4/3
    integ4_pi = (integ4_pi - yy[0]+yy[-1]) * dx*4*3/8  
# and now MteCarlo. 
# notice how we avoid any Python loops for efficiency.
    plt.figure(figsize=[7,6.5])
    p = plt.plot(xx,yy)
    plt.axis(limits=[-.3,1.3,0,1])
    inside = xr*xr+yr*yr < 1.0   # Boolean array
    co = inside.astype(int)      # (optionally) cast Boolean as int
    mypi = 4 * co.sum()/N        # summation of quarter-circle area 
    sig  = N**(-0.5)*mypi        # estimate of uncertainty  
     # LaTeX style text
    label = str(np.around(mypi,5))+' + '+ \
        str(np.around(np.pi-mypi,5))+' ('+ \
        str(np.around(sig,5))+')  Monte Carlo, N =  '+str(N) 
    logN = np.log10(N)  #  needed only for nicer plotting
	# this conditional cut speed up calculations by pruning points
    if (logN > 5.5): 
        xr = xr[0:200000]; yr = yr[0:200000]; co = co[0:200000]
    plt.scatter(xr,yr, s=7-logN//2, c=255-co*120, alpha=0.95-logN*.05)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title ('$\pi = $'+label)
    plt.grid(True)
    #print("MtC integration  "+label)
    print("MtC integr  ",mypi,     '+',np.pi-mypi, 'N =',N)
    print("integr. 1st ",integ1_pi,'+',np.pi-integ1_pi)
    print("integr. 2nd ",integ2_pi,'+',np.pi-integ2_pi) 
    print("integr. 3nd ",integ3_pi,'+',np.pi-integ3_pi) 
    print("integr. 4th ",integ4_pi,'+',np.pi-integ4_pi) 
# end of function


'''
 Main (driver) program

 Notice a well stuctured, functional programmig style: 
 a rather small driver program including interaction with user
 and self-sufficient function(s) that is called for different N.
 Here, the whole output is handled by the function. Thus it  does not 
 return any values, which can be made explicit by an  optional return 
 statement, returning None (nothing, void in C/C++). This and/or 
 the final non-indented comment line are needed to visually 
 separate the function from other functions or Main.
 The comment is sort of arrow-shaped by chance.
'''  
for n in range(1,9):
    N = 10**n
    MteCarlo_4Int(N)
    input(" next? ")
plt.close()