#! /usr/bin/env python
"""
 Mcint.py is a program  which performs Monte Carlo
 integration of sin(x *y**2) from x,y=0,0 to x,y = pi,pi

 The integration is performed many times for various numbers of random
 samples.  Graphs for the statistical average and SDOM are 
 displayed which exhibits the mean-value integration error decreasing
 linearly with 1/sqrt(N). 


 Paul Eugenio
 PHZ4151C
 20 Feb 2019
"""

from __future__ import division, print_function
import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt



def integrateMC(func, dim, limit, N=100):
    """
    Monte carlo mean-value integration of any dimension.
    Func is the user defined integrand function 
    which has a list argument x containing dim values. 

    for example: sin(x*y) is f(x)= sin(x[0]*x[1])

    Dim is the number of dimensions, limit is a list of [a,b] 
    values containg the intgration limits for  all dimensions, and 
    N is the number of Monte carlo sampling points.
    """
    I = 1/N
    sum = 0
    for n in range(dim):
        I *= (limit[n][1] - limit[n][0])
        
    for k in range(N):
        x = []
        for n in range(dim):
            x += [limit[n][0] + (limit[n][1] - limit[n][0])*rnd.random()]
        sum += func(x)
    return I*sum

    
def f(x):
    """ integrand function sin(xy**2)"""
    return np.sin(x[0] * x[1]**2)


#
# main
dim, limit = 2, [[0,np.pi],[0,np.pi]]
nSamples = 10000         # starting value for number of MC samples

nPoints,avgPoints,sdomPoints = [],[],[] # lists for plotting 

n = 20            # the number of times each integration is performed to
                  # define average and SDOM 

while nSamples < 2e6:
    # loop over values of N and average results over n-number of repeated calculations
    vSum = 0.0
    vSumSq = 0.0
    for _ in range(n):
        v = integrateMC(f, dim, limit, nSamples)
        vSum += v
        vSumSq += v*v
    
    avg = vSum/n
    SDOM = np.sqrt(1/n) * np.sqrt( 1/(n - 1) * (vSumSq - n*avg*avg) )
    print(nSamples, avg, SDOM)
    nSamples *= 2
    nPoints += [nSamples]
    avgPoints += [avg]
    sdomPoints += [SDOM]

# 
# generate plots
fig1 = plt.figure(1)
plt.plot(nPoints,avgPoints,"-o")
plt.xlabel("N")
plt.ylabel("AVG")
fig1.savefig("avg"+str(n)+".png")
fig1.show()

fig2 = plt.figure(2)
plt.plot(1/np.sqrt(nPoints),sdomPoints,"-o")
plt.xlabel(r"$1/\sqrt{N}$")
plt.ylabel("SDOM")
fig2.savefig("sdom"+str(n)+".png")
fig2.show()

# pause before program exist
raw_input("Enter [return] to exit")