# -*- 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 "f-average" 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))

# The f-average method requires that we repeatedly generate x_i random
# numbers (n times) uniformly distributed between the integration interval
# (0 to 1 in this example).  For each x_i we calculate the corresponding
# f(x_i) and sum all of these to find f_total.  f_average is obtained
# from f_average = f_total/n.  An estimate of the integral is given by
# the integration interval (b-a)*f_average.  In the plot below, we show
# each of the f(x_i) and f_average which is our estimate of the integral.
# The square area below the red line is approximately equal to the area
# beneath the blue curve which is the integral that we're trying to evaluate.

n = int(1e4)
ftot = 0
xList = []
fList = []
for i in range(n):
    xx = np.random.uniform()
    xList = xList + [xx]
    fList = fList + [fcn(xx)]
    ftot = ftot + fcn(xx)
favg = ftot/n
print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', favg)

# Plot the function as  a line and the add the hits and misses as scatter
# plots.
plt.figure()
plt.plot(xList, fList, 'bo', fillstyle = 'none')
plt.plot(xList, favg*np.ones(n), 'r-', linewidth = 3)
plt.xlabel('x')
plt.ylabel('f(x)')

# If you don't need to keep track of the actual coordinates of each
# (x, f(x)), then you can get away with a more compact loop keeps a running
# total of all the f(x) values.
n = int(1e4)
ftot = 0
for i in range(n):
    xx = np.random.uniform()
    ftot = ftot + fcn(xx)
favg = ftot/n
print('The Hit & Miss estimate of the integral of f(x) on 0 < x < 1 is', favg)

# Now we will numberically approximate the integral using n = 10,000 one
# 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):
    ftot = 0
    if j % int(m/10) == 0:
        print(int(j/10), '%', sep = '')
    for i in range(n):
        xx = np.random.uniform()
        ftot = ftot + fcn(xx)
    intList = intList + [ftot/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 18 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):
        ftot = 0
        for i in range(n):
            xx = np.random.uniform()
            ftot = ftot + fcn(xx)
        intList = intList + [ftot/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.04))

# 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, .11, 0, 0.04))