import numpy as np
import matplotlib.pyplot as plt
#This just needed for the Notebook to show plots inline. 
%matplotlib inline 

#Define the function
def f(x):
    fx = (np.sin(1/(x*(2-x))))**2
    return fx

#Integrate the function from x=0-2
#Note that you need to know the maximum value of the function
#over this range (which is y=1), and therefore the area of the box
#from which we draw random number is A=2. 
N=100000
k=0
for i in range(N):
    x=2*np.random.random()
    y=np.random.random()
    if y<f(x):
        k+=1
A=2.        
I=A*k/N
print("The integral is equal to I = ",I)

"""
Explanation: In Class Exercise: Monte Carlo Integration
Write a programme following the method of Monte Carlo Integration to calculate 
$$ I = \int_0^2 \sin^2 [\frac{1}{x(2-x)}] dx. $$
As you will need to calculate $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ many times please write a user defined function for this part of your programme.
End of explanation
"""


#Calculate the error: 

sigmaI = np.sqrt(I*(A-I))/np.sqrt(N)
print("The integral is equal to I = ",I)
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))

"""
Explanation: In Class Exercise
Returning to the programme you just wrote to calculate the integral of $f(x) = \sin^2 [\frac{1}{x(2-x)}]$, please now add an estimate of the accuracy of your integration. Using
$$\sigma_I = \frac{\sqrt{I(A-I)}}{\sqrt(N)} $$
End of explanation
"""


#Draw N values from the distribution, and calculate their mean to 
#use the mean method for integration. 
xvalues=[]
fvalues=[]
for i in range(N):
    x=2*np.random.random()
    y=np.random.random()
    if y<f(x):
        fvalues.append(f(x))
        xvalues.append(x)

fmean=np.mean(fvalues)
Imean = 2*fmean
Imeansigma = 2*np.sqrt(np.var(fvalues))/np.sqrt(len(fvalues))

print("Mean Method: ")
print("The integral is equal to I = {0:.4f}".format(Imean))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(Imeansigma))
print("The percent error is: {0:.2f} ".format(100*Imeansigma/Imean))
print("**********************")
print("Compare to the `hit or miss` Monte Carlo Method: ")
print("The integral is equal to I = {0:.4f}".format(I))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))
print("The percent error is: {0:.2f} ".format(100*sigmaI/I))

#Using mpl_style file.
import mpl_style
plt.style.use(mpl_style.style1)

plt.plot(xvalues,fvalues,'.')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()


np.mean?

"""
Explanation: In Class Exercise
Now calculate the integral of $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ using the mean value method, and compare the accuracy from 10000 sample to what you obtained previously.
End of explanation
"""


