# -*- coding: utf-8 -*-
"""
Created on Fri Oct  2 12:48:44 2020

@author: jbobowsk
"""

# This tutorial will numerically evaluate a definite integral using the
# Monte Carlo Hit & Miss method.
import numpy as np
import matplotlib.pyplot as plt

# In all Monte Carlo simulations it is necessary to generate random or
# pseudo-random numbers.  The Python command 'np.random.uniform()'
# from the NumPy module generates random numbers uniformly distributed 
# between zero and one. 
print(np.random.uniform())

# To confirm that the numbers generated are uniformly distribted on
# the interval [0, 1], we can generate a long list of these numbers and
# plot a histogram. 'np.random.uniform(a, b, n)' returns n random number
# uniformly distributed between a and b.
nos = np.random.uniform(0, 1, int(1e4))

# Here's the histogram.
nbins = 50
plt.hist(nos, nbins, color = 'c', edgecolor = 'k')
plt.xlabel('x');
plt.ylabel('counts');

# This tutorial will first attempt to evaluate an integral for which the
# solution is easily obtained.  This approach has been taken purposely so
# that we can confirm that our numerical techniques are reliable.  The
# function that we will integrate is a simply polynomial.  Below, the
# function is plotted on the interval x = [0, 1] and the exact value of the
# integral is evaluated over the same interval.  We are going to use this
# polynomial over and over again in this tutorial, so we define it as a
# function.
fcn = lambda x : (-65536*x**8 + 262144*x**7 - 409600*x**6 + 311296*x**5 -\
                  114688*x**4 + 16384*x**3)/27

xx = np.arange(0, 1, 0.001)
plt.figure()
plt.plot(xx, fcn(xx), 'k-', linewidth = 1.2)
plt.xlabel('x')
plt.ylabel('f(x)')

# Here's how we can use the SymPy module to do symbolic integration in 
# Python.  This is an indefinite integral.
import sympy as sym
x = sym.Symbol('x')
fint = sym.integrate(fcn(x), x)
print(fint)

# Here is the definite integral of our function over 0 < x < 1.
fint = sym.integrate(fcn(x), (x, 0, 1))
print(fint)
print('The exact integral of f(x) on 0 < x < 1 is', sym.N(fint))

# We now implement is the Hit & Miss numerical integration technique.  In
# this method pairs of random numbers (x_i, y_i) will be generated.  For
# the x-coordinate the x_i values will be between 0 and 1 (the integration
# interval).  Notice that our function of interest always lies between
# y = 0..1.  Therefore, we will also restrict our y_i values to be between
# 0 and 1.  The randomly generated points (x_i, y_i) have equal probability
# of landing anywhere in the box which has area 1.  The probability of a
# point landing beneath the function is equal to the area A beneath the
# curve (which is just the integral of the function of interest) divided
# by the area of the box.  Therefore, if we can determine the probability
# of a point landing beneath the curve, we can easily approximate the value
# of the definite integral.  The probability will be estimated simply by
# spraying n points into the box and counting how many land beneath the
# curve Zn (i.e. counting the Hits).  Then the probability p is simply
# p = Zn/n.
n = int(1e4)
hits = 0
misses = 0
xH = []
yH = []
xM = []
yM = []
for i in range(n):
    xx = np.random.uniform()
    yy = np.random.uniform()
    if yy <= fcn(xx):
        hits += 1
        xH = xH + [xx]
        yH = yH + [yy]
    else:
        misses += 1
        xM = xM + [xx]
        yM = yM + [yy]
intEst = hits/n
print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', intEst)

# Plot the function as  a line and the add the hits and misses as scatter
# plots.
plt.figure()
plt.plot(xH, yH, 'ro', fillstyle = 'none', markersize = 1.5)
plt.plot(xM, yM, 'bo', fillstyle = 'none', markersize = 1.5)
xx = np.arange(0, 1, 0.001)
plt.plot(xx, fcn(xx), 'k-', linewidth = 3)
plt.xlabel('x')
plt.ylabel('f(x)')

# If you don't need to keep track of the actual coordinates of each hit and
# miss, then you can get away with a more compact loop that just counts the
# number of hits.
n = int(1e4)
hits = 0
for i in range(n):
    xx = np.random.uniform()
    yy = np.random.uniform()
    if yy <= fcn(xx):
        hits += 1
print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', hits/n)

# Now we will numberically approximate the integral using n = 10,000 five 
# thousand times and plot the resulting distribution of our determination of
# the integral.  The distribution is expected to be Gaussian.

# To keep track of how long this block of codes takes to complete, I'll use the
# "datetime.now()" function to display the time before and after the loop.
from datetime import datetime
print(datetime.now())
n = int(1e4)
m = int(1e3)
intList = []
for j in range(m):
    hits = 0
    if j % int(m/10) == 0:
        print(int(j/10), '%', sep = '')
    for i in range(n):
        xx = np.random.uniform()
        yy = np.random.uniform()
        if yy <= fcn(xx):
            hits += 1
    intList = intList + [hits/n]
print(datetime.now())
plt.figure()
nbins = 20
plt.hist(intList, nbins, color = 'c', edgecolor = 'k')
plt.xlabel('integral')
plt.ylabel('counts')

# The width of the distribution can be calculated from the standard
# deviation of the list of 1000 approximations of the integral and is an
# estimate in the uncertainty of our determination of the definite integral.

# The standard deviation can be calculated easily using the statistics module.
import statistics
statistics.stdev(intList)

# Finally, the last thing we'll attempt to do is to understand how our
# uncertainly (standard deviation) depends on the n.  So far, all of our
# calculations have used n = 1e4.  Now we'll determine the standard
# deviation for values of n that range from 100 to 1e5.  This
# block of code took my laptop about 30 minutes to complete...
# If you've got a lot time to kill, trying running the same loop in Maple!
print(datetime.now())
nList = [100, 200, 500, 1000, 2000, 5000, 1e4, 2e4, 5e4, 1e5]
sigmaList = []
for k in range(len(nList)):
    intList = []
    print(k+1, 'of', len(nList))    
    n = int(nList[k])
    for j in range(m):
        hits = 0
        for i in range(n):
            xx = np.random.uniform()
            yy = np.random.uniform()
            if yy <= fcn(xx):
                hits += 1
        intList = intList + [hits/n]
    sigmaList = sigmaList + [statistics.stdev(intList)]
print(datetime.now())

# Below we plot the uncertainty in the numerical integral estimation as
# a function of n (the number of trials in the Monte Carlo simulation).
# As expected, the uncertainty decreases as the number of trials increases.
# The sigma values are proportional to 1/sqrt(n).
plt.figure()
plt.plot(nList, sigmaList, 'bo', fillstyle = 'none')
plt.xlabel('n')
plt.ylabel(r'$\sigma$ (std)')
plt.axis((-5e3, 1.2e5, 0, 0.055))

# To get a better appreciation of the dependence of sigma on n, below we
# plot sigma as a function of 1/sqrt(n) and observe the linear relationship
# between the two.  Beautiful!  All of this generated from uniformly
# distributed random numbers!  Take a moment to reflect on what we've
# accomplished.  We've used Monte Carlo simulations to study the behaviour
# of Monte Carlo simulations!  The objective of a Monte Carlo calculation
# is always to study the characteristics of some system (often a physical
# system) by simulating data using random numbers.
plt.figure()
plt.plot(1/np.sqrt(nList), sigmaList, 'bo', fillstyle = 'none')
plt.xlabel(r'$1/\sqrt{n}$')
plt.ylabel(r'$\sigma$ (std)')
plt.axis((0, .12, 0, 0.055))