#!/usr/bin/env python

# ---------------------------------------------------------------------------------- #
# 
#  Python script for estimating the value of pi using a Monte Carlo technique.
#
#  The program picks random points in the area [-1,1] x [-1,1], and determines which
#  fraction of these are within the unit circle. This in turn gives a measure of pi
#  with associated statistical uncertainty. Performing such "experiments" many times
#  not only gives the value of pi, but also a feel for the use of Monte Carlo, and
#  experience in calculating averages, RMSs, and the error on these.
#
#  For more information see:
#    P. R. Bevington: page 75-78
# 
#   Author: Troels C. Petersen (NBI)
#   Email:  petersen@nbi.dk
#   Date:   3rd of December 2017
# 
# ---------------------------------------------------------------------------------- #


from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
plt.close("all")

SavePlots = False         # Determining if plots are saved or not
r = np.random
r.seed(4)

# Set parameters:
Nexperiments = 100        # Number of "experiments" determining pi
Npoints      = 1000       # Number of points per experiment in determining pi

pi_true = 3.14159265358979323846264338328 # np.pi

List_PiDist    = np.zeros(Nexperiments)
List_HitDist_x = np.zeros(Npoints)
List_HitDist_y = np.zeros(Npoints)




#----------------------------------------------------------------------------------- #
# Loop over process:
#----------------------------------------------------------------------------------- #

sum1 = 0.0
sum2 = 0.0

for iexp in range ( Nexperiments ) :
    Nhit = 0
    for j in range ( Npoints ) :
      
        # Get two random numbers uniformly distributed in [-1,1]:
        x1 = 2.0 * (r.uniform() - 0.5)
        x2 = 2.0 * (r.uniform() - 0.5)

        # Test if (x1,x2) is within unit circle:
        if (x1*x1 + x2*x2 < 1.0) : Nhit += 1

        # Plot 2D distribution of (x1,x2) for the first experiment:
        if (iexp == 0) :
            List_HitDist_x[j] = x1
            List_HitDist_y[j] = x2

    
    # Calculate the fraction of points within the circle and its error:
    f_hit = 0.0      # Calculate yourself!
    ef_hit = 0.0     # Calculate yourself!

    # From this we can get pi and its error:
    pi_estm = 0.0    # Calculate yourself!
    pi_error = 0.0   # Calculate yourself!

    # Plot these values and once outside the loop, calculate mean and RMS of this histogram:
    # Yes - also do this yourself!
    List_PiDist[iexp] = pi_estm

    # Print first couple of pi estimates:
    if (iexp < 5) :
        print("  {0:2d}. pi estimate:   {1:7.4f} +- {2:6.4f}".format(iexp+1, pi_estm, pi_error))


# ----------------------------------------------------------------------------------- #
# Calculate average and RMS, and print the result compared to the true value of pi:
if (Nexperiments > 1) :
    pi_ave = 0.0    # Fill this out yourself!!!
    pi_rms = 0.0    # ...and this!
    epi_ave = 0.0

# Print how many sigmas your final result is away from the true value of pi.
    if (epi_ave > 0.0) :
        print("  Average pi estimate: {0:6.4f} +- {1:6.4f}    ({2:4.2f} sigma)".format(pi_ave, epi_ave, (pi_ave - pi_true)/epi_ave))


# ----------------------------------------------------------------------------------- #
# Plot the histograms:
# ----------------------------------------------------------------------------------- #

# Distribution of points from one experiment:
fig1, ax1 = plt.subplots(figsize=(8, 6))     # binx,biny     x1 , x2    y1 , y2      or "Reds"
h = ax1.hist2d(List_HitDist_x, List_HitDist_y, (200,200), ((-1.0,1.0),(-1.0,1.0)), cmap="binary")
plt.colorbar(h[3]) # z-scale on the right of the figure

plt.tight_layout()
plt.show(block=False)

if (SavePlots) : fig1.savefig("HitDist.pdf")



# Finally, ensure that the program does not termine (and the plot disappears), before you press enter:
try:
    __IPYTHON__
except:
    raw_input('Press Enter to exit')



# ----------------------------------------------------------------------------------- #
# 
# First acquaint yourself with the program, and make sure that you understand what the
# parameters "Nexperiment" and "Npoints" refere to! Also, before running the program,
# calculate what precision you expect on pi in each experiment, when using the number
# of points chosen in the program (i.e. 1000 points).
# 
# Then, run the program, and then take a look at the result... which requires that you
# fill in the calculations yourself!
# 
# 
# Questions:
# ----------
#  1) Try to run 100 experiments with 1000 points in each. What is the approximate
#     uncertainty on pi in each experiment? Does that fit, what you calculated before
#     running the program? What is the uncertainty on the AVERAGE of all 100 experiments?
# 
#  2) How do you expect the values of pi to distribute themselves? And is this the
#     case in reality?
# 
#  3) Does it make any difference on the precision of the final pi value, whether
#     you make many experiments with few points, or one experiment with many points,
#     as long as the product of Nexperiment x Npoints remains constant?
# 
# The following problem is the REAL question in this exercise!
# 
#  4) Now try to use this method in three dimensions to estimate the constant in
#     front of the r^3 expression for the volume. Do you get 4/3*pi?
#     Increase the dimensionality (say up to 10), and see if you can figure out
#     the constants needed to calculate the hyper-volumes!
# 
# HINT: I'll reveal, that for Ndim of 4 and 5, the constant contains pi^2 and some
#       simple rational fraction, while for Ndim 6 and 7, it contains pi^3 and a
#       semi-simple rational fractions.
#
# ----------------------------------------------------------------------------------- #