#!/usr/bin/env python
# coding: utf-8

# In[1]:


from math import pi,sqrt
import random
import matplotlib.pyplot as plt
import numpy as np


# ## The number _pi_ ($\pi$)
# The number $\pi$ is a mathematical constant. Originally defined as the ratio of a circle's circumference to its diameter, it now has various equivalent definitions and appears in many formulas in all areas of mathematics and physics. It is approximately equal to 3.14159. It has been represented by the Greek letter __"$\pi$"__ since the mid-18th century, though it is also sometimes spelled out as "pi". It is also called Archimedes' constant.
# 
# Being an irrational number, $\pi$ cannot be expressed as a common fraction (equivalently, its decimal representation never ends and never settles into a permanently repeating pattern). 

# ### What is the value of $\pi$

# In[2]:


print("Value of pi: ",pi)


# ## What is the logic behind computing pi by throwing dart randomly?
# Imagine a square dartboard.
# 
# Then, the dartboard with a circle drawn inside it touching all its sides.
# 
# And then, you throw darts at it. Randomly. That means some fall inside the circle, some outside. But assume that no dart falls outside the board.
# 
# ![boards](https://raw.githubusercontent.com/tirthajyoti/Stats-Maths-with-Python/master/images/boards.png)
# 
# At the end of your dart throwing session, you count the fraction of darts that fell inside the circle of the total number of darts thrown. Multiply that number by 4.
# 
# The resulting number should be pi. Or, a close approximation if you had thrown a lot of darts.
# 
# The idea is extremely simple. If you throw a large number of darts, then the **probability of a dart falling inside the circle is just the ratio of the area of the circle to that of the area of the square board**. With the help of basic mathematics, you can show that this ratio turns out to be pi/4. So, to get pi, you just multiply that number by 4.
# 
# The key here is to simulate the throwing of a lot of darts so as to make the fraction equal to the probability, an assertion valid only in the limit of a large number of trials of this random event. This comes from the [law of large number](https://en.wikipedia.org/wiki/Law_of_large_numbers) or the [frequentist definition of probability](https://en.wikipedia.org/wiki/Frequentist_probability).
# 
# See also the concept of [Buffon's Needle](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem)

# ### Center point and the side of the square

# In[3]:


# Center point
x,y = 0,0
# Side of the square
a = 2


# ### Function to simulate a random throw of a dart aiming at the square

# In[4]:


def throw_dart():
    """
    Simulates the randon throw of a dirt. It can land anywhere in the square (uniformly randomly)
    """
    # Random final landing position of the dirt between -a/2 and +a/2 around the center point
    position_x = x+a/2*(-1+2*random.random())
    position_y = y+a/2*(-1+2*random.random())
    
    return (position_x,position_y)


# In[5]:


throw_dart()


# ### Function to determine if the dart landed inside the circle

# In[6]:


def is_within_circle(x,y):
    """
    Given the landing coordinate of a dirt, determines if it fell inside the circle
    """
    # Side of the square
    a = 2
    
    distance_from_center = sqrt(x**2+y**2)
    
    if distance_from_center < a/2:
        return True
    else:
        return False


# In[7]:


is_within_circle(1.9,1.9)


# In[8]:


is_within_circle(1.2,1.9)


# In[9]:


is_within_circle(0.4,-0.74)


# In[10]:


r1,r2=throw_dart()
print(r1,r2)
if is_within_circle(r1,r2):
    print("This one landed inside the circle!")
else:
    print("This one did not land inside the circle!")


# ### Now throw a few darts!

# In[11]:


n_throws = 10


# In[12]:


count_inside_circle=0


# In[13]:


for i in range(n_throws):
    r1,r2=throw_dart()
    if is_within_circle(r1,r2):
        count_inside_circle+=1


# ### Compute the ratio of `count_inside_circle` and `n_throws`

# In[14]:


ratio = count_inside_circle/n_throws


# ### Is it approximately equal to $\pi$?

# In[15]:


print(4*ratio)


# ### Not exactly. Let's try with a lot more darts!

# In[16]:


n_throws = 10000
count_inside_circle=0


# In[17]:


for i in range(n_throws):
    r1,r2=throw_dart()
    if is_within_circle(r1,r2):
        count_inside_circle+=1


# In[18]:


ratio = count_inside_circle/n_throws


# In[19]:


print(4*ratio)


# ### Now we are approaching $\pi$ :-). Let's functionalize this process and run a number of times

# In[20]:


def compute_pi_throwing_dart(n_throws):
    """
    Computes pi by throwing a bunch of darts at the square
    """
    n_throws = n_throws
    count_inside_circle=0
    for i in range(n_throws):
        r1,r2=throw_dart()
        if is_within_circle(r1,r2):
            count_inside_circle+=1
            
    result = 4*(count_inside_circle/n_throws)
    
    return result


# ### Now let us run this experiment a few times and see what happens.

# In[22]:


n_exp=[]
pi_exp=[]
n = [int(10**(0.5*i)) for i in range(1,15)]
for i in n:
    p = compute_pi_throwing_dart(i)
    pi_exp.append(p)
    n_exp.append(i)
    print("Computed value of pi by throwing {} darts is: {}".format(i,p))


# In[43]:


plt.figure(figsize=(8,5))
plt.title("Computing pi with \nincreasing number of random throws",fontsize=20)
plt.semilogx(n_exp, pi_exp,c='k',marker='o',lw=3)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.xlabel("Number of random throws",fontsize=15)
plt.ylabel("Computed value of pi",fontsize=15)
plt.hlines(y=3.14159,xmin=1,xmax=1e7,linestyle='--')
plt.text(x=10,y=3.05,s="Value of pi",fontsize=17)
plt.grid(True)
plt.show()


# ### So, effectively we can average a few experiments and take that value

# In[31]:


n = 5000000
sum=0
for i in range(20):
    p=compute_pi_throwing_dart(n)
    sum+=p
    print("Experiment number {} done. Computed value: {}".format(i+1,p))
print("-"*75)
pi_computed = round(sum/20,4)
print("Average value from 20 experiments:",pi_computed)


# ### Error percentage can be easily computed

# In[32]:


error_pct = 100*(pi - pi_computed)/pi
print("Error percentage: ", error_pct)