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

#setting the random seed
random.seed(1234)

#defining the integrating function
def f(x): return (np.sin(x)**4*np.cos(x/3)**2+np.sin(x/4)**6)*np.exp(-x/8)*1.4

#number of extraction
nMax = 300 

#domain of definition of the function
xmin = 0. 
xmax = 13.
#range of the function values 
ymin = 0.
ymax = 1.
#area of the rectangle used for the random points extraction
A = (xmax-xmin)*(ymax-ymin)

#initialization of variables and arrays for the output
n = 0
N = 0
extractedX = []
extractedY = []
goodX = []
goodY = []
values = []

#loop for the extractions
for i in range(0,nMax):
    #increasing the counter with total extractions performed
    N +=1
    #extract a point in the range
    thisX = random.uniform(xmin,xmax)
    thisY = random.uniform(ymin,ymax)
    #add the point to the list with all the extracted points
    extractedX.append(thisX)
    extractedY.append(thisY)
    #check if the point is beyond the function
    if thisY < f(thisX):
        #if so increasing the counter with the number of points beyond the function
        n +=1
        #and adding the point to the relative list
        goodX.append(thisX)
        goodY.append(thisY)
    #list with the extimation of the integral value as a function of the number of extraction
    values.append(n/N*A)

print("n     =",n)
print("N     =",N)
print("n/N*A =",n/N*A)

#making the plots
w, h = plt.figaspect(1.)
plt.figure(figsize=(w,h), dpi=160)

plt.subplot(211)
x1 = np.arange(xmin,xmax,0.01)
plt.grid(True)
plt.xlabel('x',labelpad=0.5)
plt.ylabel('y',labelpad=0.5)
plt.scatter(extractedX, extractedY)
plt.scatter(goodX, goodY)
plt.plot(x1, f(x1), 'r--')

plt.subplot(212)
plt.grid(True)
plt.xlabel('n',labelpad=0.5)
plt.ylabel(u'n/N*A',labelpad=0.5)
plt.plot(range(0,nMax), values)

plt.show()