# -*- coding: utf-8 -*-
"""
Created on Mon Apr  6 14:43:35 2020

@author: neil
"""

import numpy as np
# Monte Carlo techniques can be used for a wide variety of problems where introducing randomness
# can improve the solution in some way

# just to get our feet wet lets go back to simple integration.  We integrated a function using the trapiziodal rule.
# lets integrate the area of a diameter 1 circle using the trapiziod and using basic Monte Carlo
# we only need to integrate a 1/4 circle (x,y > 0) and multiply by 4, that way we stay away from sign problems

# Analytic soln
print('Analytic solution of area of a unit circle,  {:.4}'.format( np.pi*(.5)**2 ))

# numerical integration
# trapizoid rule, 10 points, for a circle, y = sqrt(r**2-x^2)
area   = 0      # place to keep the area calcs
slices = 10     # number of trapizoids
X      = 0.5    # distance from 0 to max radius
delx   = X/(slices) #width of trapizoids


for x in np.arange(slices):
    x1 = x*delx
    x2 = x1 + delx
    area += delx*(np.sqrt(.25 - x1**2) + np.sqrt(.25 - x2**2))/2
        
print('Numerical trapizoidal integration area,      {:.4}'.format(4*area))
    
# there are several approaches to MC integration, here we throw darts (random guesses)
# at a circle inscribed in a square and calculate the number of hits vs the total throws, which will approx
# the ratio of the size of the circle to the size of the square.

import numpy as np
import matplotlib.pyplot as plt

# Define number of samples
N = 2000

# Define rectangle boundaries, in this case a unit square, centered on 0,0
# although we only use the quadrant x,y > 0 (area of .25 instead of 1)
rec_x = [0, .5]
rec_y = [0, .5]
# Sample from uniform distribution, this illustrated how to get a list of uniform random numbers
mc_x = np.random.uniform(min(rec_x), max(rec_x), N)
mc_y = np.random.uniform(min(rec_y), max(rec_y), N)

# Define the points under the curve
points_under = [True if mc_y[i] <= np.sqrt(.25-mc_x[i]**2) else False for i in range(len(mc_x))]
# this is called a LIST COMPREHENSION, it is a short hand way of writing a simple for loop!!
# note 'points_under' is a list of TRUE and FALSE results, this is called a BOOLEAN array
# Calculate integral by taking the ratio of points in and out to the area of the quadrant of the 
# unit square
area_rect  = (max(rec_x) - min(rec_x)) * (max(rec_y) - min(rec_y))
integral = area_rect* sum(points_under) / N
print("Monte Carlo integration of unit circle area: {:.4}".format(integral*4))    

# Plot the points to illustrate the result (some nice plotting ideas here)
x = np.linspace(0,.5,100)
fig=plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, np.sqrt(.25 - x**2), linewidth=3, c='k')
ax.scatter(mc_x[points_under], mc_y[points_under],
           c='r', s=15)
# note using the TRUE/FALSE list, the same length as the mc lists, you can  use
# mc_x[points_under] to only plot the TRUE points (this is advanced but very cute)
ax.scatter(mc_x[np.logical_not(points_under)], 
            mc_y[np.logical_not(points_under)], 
            c='b', s=15)
ax.set_xlim(rec_x)
ax.set_ylim(rec_y)
ax.axis('equal')
ax.set_title('Dart-throwing Monte Carlo Method (only 1/4 of the circle)')
# place a pretty box around the results
ax.plot([-.001,.501,.501,-.001,-.001],[-.001,-.001,.501,.501,-.001] , linewidth=3, c='g')
plt.show()