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

@author: jbobowsk
"""

# This tutorial will use Monte Carlo methods to calculate the volume of an
# n-dimensional sphere.
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())

# This tutorial will attempt to numerically determine the volume of a sphere
# of arbitrary dimension using the Hit & Miss method.  This is equivalent to
# evaluating a d-dimensional integral.  In d-dimensions, a hypersphere is
# defined by x1^2+x^2+x3^3+...+xd^2 < R^2.

# To find the volume of a sphere we'll generate random numbers over some
# "square" volume and count how many land inside the sphere.  The
# probability of the randomly generated point landing within the sphere is
# simply p = Vsphere / Vsquare.  If we use spheres of radius 1, then the
# volume of the sqaure that contains the sphere is 2^d where d is the
# number of dimensions that we're working with. The Monte Carlo simulation
# is used to determine p = Zn/n where n is the number of trials and Zn is
# the number of "Hits".  Collecting all of these results we have
# p = Zn/n = Vsphere/2^d or Vsphere = 2^d*(Zn/n).

# We'll start with the trivial example of a 1-D sphere.  In 1-D, x1^2<R^2
# simply translates to -R < x1 < +R.  Throughout this tutorial will restrict
# ourselve to the positive quadrant.  That is, we'll generate random  values
# of x1 between [0, R] (where R = 1 in this tutorial).  The 1-D sphere is
# trivial because ALL of the randomly generated numbers will fall within
# the sphere.  The volume of the 1-D sphere is just the length of the line.
#  V1D = 2R.

d = 1 # The number of dimensions
n = int(1e5) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x1 = np.random.uniform()
    R = np.sqrt(x1**2) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V1D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V1Dexact = 2 # The known value of the volume of the sphere
ratio = V1D/V1Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V1D)
print('The exact volume is', V1Dexact, 'and the ratio of the two is', V1D/V1Dexact)

# Next a 2-D sphere (a circle).  In 2-D, x1^2+x2^2 < R^2.  Now we'll have
# to generate random values for x1 and x2 both being between [0, R]
# (where R = 1 in this tutorial).  The volume of the 2-D sphere is expected
# to be the area of a circle. V2D = pi*R^2.  One of the virtues of the
# Monte Carlo method is that going to higher dimensions requires only VERY
# minor changes to our simulation code.  Note that the example of the area
# of a circle is often used as a Monte Carlo calculation of the numerical
# value of pi.
d = 2 # The number of dimensions
n = int(1e5) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x1 = np.random.uniform()
    x2 = np.random.uniform()
    R = np.sqrt(x1**2 + x2**2) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V2D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V2Dexact = np.pi; # The known value of the volume of the sphere
ratio = V2D/V2Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V2D)
print('The exact volume is', V2Dexact, 'and the ratio of the two is', V2D/V2Dexact)

# A 3-D sphere is a ball. Here, x1^2+x2^2+x3^2 < R^2.  Now we'll have to
# generate random values for x1, x2, and x3 all being between [0, R]
# (where R = 1 in this tutorial).  The expected volume of the 3-D sphere is
# (4/3)*pi*R^3.
d = 3 # The number of dimensions
n = int(1e5) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x1 = np.random.uniform()
    x2 = np.random.uniform()
    x3 = np.random.uniform()
    R = np.sqrt(x1**2 + x2**2 + x3**2) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V3D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V3Dexact = (4/3)*np.pi # The known value of the volume of the sphere
ratio = V3D/V3Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V3D)
print('The exact volume is', V3Dexact, 'and the ratio of the two is', V3D/V3Dexact)

# Now we're onto something more interesting.  We cannot picture a 4-D sphere,
# but mathematically the definition is obvious: x1^2+x2^2+x3^2+x4^2 < R^2.
# Now we'll have to generate random values for x1, x2, x3, and x4 all being
# between [0, R] (where R = 1 in this tutorial).  So that we can compare our
# Monte Carlo simulation to expected results, we'll note that a 4-D sphere
# has a volume of (1/2)*pi^2*R^4.
d = 4 # The number of dimensions
n = int(1e5) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x1 = np.random.uniform()
    x2 = np.random.uniform()
    x3 = np.random.uniform()
    x4 = np.random.uniform()
    R = np.sqrt(x1**2 + x2**2 + x3**2 + x4**2) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V4D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V4Dexact = (1/2)*np.pi**2 # The known value of the volume of the sphere
ratio = V4D/V4Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V4D)
print('The exact volume is', V4Dexact, 'and the ratio of the two is', V4D/V4Dexact)

# We'll do two more examples.  First, a 5-D sphere: x1^2+x2^2+x3^2+x4^2+x5^2 < R^2.
# Now we'll have to generate random values for x1, x2, x3, x4, and x5 all
# being between [0, R] (where R = 1 in this tutorial).  So that we can
# compare our Monte Carlo simulation to expected results, we'll note that a
# 5-D sphere has a volume of (8/15)*pi^2*R^5.
d = 5 # The number of dimensions
n = int(1e5) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x1 = np.random.uniform()
    x2 = np.random.uniform()
    x3 = np.random.uniform()
    x4 = np.random.uniform()
    x5 = np.random.uniform()
    R = np.sqrt(x1**2 + x2**2 + x3**2 + x4**2 + x5**2) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V5D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V5Dexact = (8/15)*np.pi**2 # The known value of the volume of the sphere
ratio = V5D/V5Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V5D)
print('The exact volume is', V5Dexact, 'and the ratio of the two is', V5D/V5Dexact)

# By now you should have noticed how easy it is to extend the number of
# dimensions.  To prove a point, our last example will be used to find the
# volume of a 9-D sphere: x1^2+x2^2+x3^2+x4^2+x5^2+x6^2+x7^2+x8^2+x9^2 < R^2.
# Now we'll have to generate random values for x1, x2, x3, x4, x5, x6, x7, x8
# and x9 all being between [0, R] (where R = 1 in this tutorial).  Keep in
# mind that we are now effectively evaluating a 9-D integral with a tricky
# restriction on the integration limits.  So that we can compare our
# Monte Carlo simulation to expected results, we'll note that a 9-D sphere
# has a volume of (32/945)*pi^4*R^9. To save some typing, I'll define a
# vector of 9 random numbers...
d = 9 # The number of dimensions
n = int(1e6) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x = np.random.uniform(0, 1, d)
    Rsq = 0
    for j in range(d):
       Rsq = Rsq + x[j]**2
    R = np.sqrt(Rsq) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn
V9D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
V9Dexact = (32/945)*np.pi**4 # The known value of the volume of the sphere
ratio = V9D/V9Dexact # Compare the calculated volume to the known volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V9D)
print('The exact volume is', V9Dexact, 'and the ratio of the two is', V9D/V9Dexact)

# Below we import and display an image from Wikipedia that shows the
# surface areas and volumes of "n-spheres" up to n = 9
# (https://en.wikipedia.org/wiki/N-sphere).
import matplotlib.image as mpimg
img = mpimg.imread('hypersphere.png')
plt.figure(figsize = (12, 7))
imgplot = plt.imshow(img)

# Now let's do something more interesting...  Let's calculate the volume of
# a high-dimensional sphere -- one that we don't know ahead of time the
# correct answer.  Let's arbitrarily take d = 12.
d = 12 # The number of dimensions
n = int(1e6) # The number of Monte Carlo trials
Zn = 0
for i in range(n):
    x = np.random.uniform(0, 1, d)
    Rsq = 0
    for j in range(d):
       Rsq = Rsq + x[j]**2
    R = np.sqrt(Rsq) # Calculate distance from the orgin
    if R <= 1: # Check for a hit
        Zn += 1 # If there's a hit increment Zn

# Here's the result!  The volume of this sphere is this prefactor times the
# radius of the sphere to the 12th power.
V12D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume
print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V12D)

# A natural question is: What is the uncertainty in this value?  Well, we
# can run this d = 12 simulation many times and plot a distribution of the
# calculated prefactors.  Below, we run the same simulation 400 times.
d = 12
n = int(1e6)
maxk = 400
from datetime import datetime
print(datetime.now())
V12List = []
for k in range(maxk):
    Zn = 0
    for i in range(n):
        x = np.random.uniform(0, 1, d)
        Rsq = 0
        for j in range(d):
            Rsq = Rsq + x[j]**2
        R = np.sqrt(Rsq)
        if R <= 1:
            Zn += 1
    if k % 40 == 0: # This if statement is used to print the value of k
        # every 40th iteration.  It helps track the progress of the loop
        # during execution.
        print(int(k/4), '%', sep = '')
    V12List = V12List + [2**d*Zn/n]
print(datetime.now())
plt.figure()
nbins = 20
plt.hist(V12List, nbins, color = 'c', edgecolor = 'k')
plt.xlabel('Volume of 12-D Sphere with R = 1')
plt.ylabel('Counts')

# An estimate of the error in a Monte Carlo simulation of the volume of a
# sphere in 12-D is given by the standard deviation of the V12List
# distribution.
import statistics as stats
std = stats.stdev(V12List)

# However, we can do better!  A better estimate of 12-D sphere prefactor is
# given by the mean of V12List.
prefactor12D = stats.mean(V12List)
print('12-D sphere volume prefactor:', prefactor12D)
print('The standard deviation:', std)


# If you've taken PHYS 232, you know that the error in the mean of a
# distribution is given by the standard deviation (sigma) divided the
# square root of number of trials used to produce the distribution
# (N = maxk = 400 in our case).  This is called the standard error.
stdError = std/np.sqrt(maxk)
print('The standard error:', stdError)

# So we should be comfortable reporting that the prefactor used to
# calculate the volume of a 12-D sphere is:
print('The prefactor for the volume of a 12-D sphere is', prefactor12D,\
      '\u00B1', stdError)
# Thus, we've calculated this prefactor with a prefactor with a precision
# better than 1%.

# If you've looked at the "Hit & Miss" MATLAB tutorial, you would have seen
# that the error in a Hit & Miss integration decreases as 1/sqrt(N).  We
# first did 1e6 iterations of our calculation of the volume of a 12-D sphere.
# For this calculation the error is given by the standard deviation of the
# distribution that we calculated.  However, in calculating the
# distribution we actually ran a total of 400 million trials of our
# simulation (400*1e6).  Therefore, with a factor of 400 more trials, we
# should expect the error in our final volume (the mean of the
# distribution) to go down by a factor of sqrt(400) = 20.  This is exactly
# what happened when we calculated the standard error of our distribution.
# Everything is beautifully consistent!