import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0., 2. * np.pi, 100)
x = np.cos(t) + np.cos(2. * t)
y = np.sin(t)

N = 100
rand = np.array([np.random.uniform(low=-3, high=3, size=N), np.random.uniform(low=-3, high=3, size=N)]).T

fig, ax = plt.subplots(1, 1, figsize=(7, 7))
ax.plot(rand[:,0], rand[:,1], '.k')
ax.plot(x, y, "-r", linewidth=2)
ax.plot([-3, -3, 3, 3, -3], [-3, 3, 3, -3, -3], "-b", linewidth=2)
ax.set_axis_off()
ax.set_xlim([-4, 4])
ax.set_ylim([-4, 4])

plt.show()

"""
Explanation: Calcule numérique avec la méthode Monte-Carlo (première partie)
TODO
- traduire en français certaines phrases restées en anglais
Approximation numérique d'une surface avec la méthode Monte-Carlo
$\newcommand{\ounk}{{\color{red}{O_1}}}$
$\newcommand{\aunk}{{\color{red}{\mathcal{A}_1}}}$
$\newcommand{\nunk}{{\color{red}{\mathcal{N}_1}}}$
$\newcommand{\okn}{{\color{blue}{O_2}}}$
$\newcommand{\akn}{{\color{blue}{\mathcal{A}_2}}}$
$\newcommand{\nkn}{{\color{blue}{\mathcal{N}_2}}}$
Conceptuellement, pour calculer la surface $\aunk$ d'un objet $\ounk$ avec la methode Monte-Carlo, il suffit:
1. de placer cet objet $\ounk$ entièrement dans une figure géométrique $\okn$ dont on connait la surface $\mathcal \akn$ (par exemple un carré ou un rectangle)
2. tirer aléatoirement un grand nombre de points dans cette figure  $\okn$ (tirage uniforme)
3. compter le nombre de points $\nunk$ tombés dans l'objet  $\ounk$ dont on veut calculer la surface
4. calculer le rapport $\frac{\nunk}{\nkn}$ où $\nkn$ est le nombre total de points tiré aléatoirement (en multipliant par 100 ce rapport on obtient le pourcentage de points tombés dans l'objet $\ounk$)
5. appliquer ce rapport à la surface $\mathcal \akn$ de la figure englobante $\okn$ (le carré, rectangle, ... dont on connait la surface) pour obtenir la surface $\aunk$ recherchée: $\aunk \simeq \frac{\nunk}{\nkn} \mathcal \akn$
End of explanation
"""


import sympy as sp
sp.init_printing()

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

"""
Explanation: Le même principe peut être appliqué pour calculer un volume.
Cette methode très simple est parfois très utile pour calculer la surface (ou le volume) de figures géometriques complexes. En revanche, elle suppose l'existance d'une procedure ou d'une fonction permettant de dire si un point est tombé dans l'objet $O_2$ ou non.
Application au calcul d'intégrales
Calculer l'intégrale d'une fonction, revient à calculer la surface entre la courbe décrivant cette fonction et l'axe des abscisses (les surfaces au dessus de l'axe des abscisses sont ajoutées et les surfaces au dessous sont retranchées).
Exemple written in Python
End of explanation
"""


def f(x):
    return -x**2 + 3.

"""
Explanation: The function to integrate
End of explanation
"""


N = 100000           # The number of random points

x_lower_bound = -4.0
x_upper_bound = 4.0
y_lower_bound = -16.0
y_upper_bound = 16.0

random_points = np.array([np.random.uniform(low=x_lower_bound, high=x_upper_bound, size=N),
                          np.random.uniform(low=y_lower_bound, high=y_upper_bound, size=N)]).T

"""
Explanation: Random points
End of explanation
"""


# Points between f and the abscissa
random_points_in_pos = np.array([p for p in random_points if 0 <= p[1] <= f(p[0])])
random_points_in_neg = np.array([p for p in random_points if 0 >  p[1] >= f(p[0])])

ratio_pos = float(len(random_points_in_pos)) / float(N)
ratio_neg = float(len(random_points_in_neg)) / float(N)
print('Percentage of "positive" points between f and the abscissa: {:.2f}%'.format(ratio_pos * 100))
print('Percentage of "negative" points between f and the abscissa: {:.2f}%'.format(ratio_neg * 100))

s2 = (x_upper_bound - x_lower_bound) * (y_upper_bound - y_lower_bound)
print("Box surface:", s2)

s1 = ratio_pos * s2 - ratio_neg * s2
print("Function integral (numerical computation using Monte-Carlo):", s1)

"""
Explanation: Numerical computation of the integral with Monte-Carlo
End of explanation
"""


x = sp.symbols("x")

integ = sp.Integral(f(x), (x, x_lower_bound, x_upper_bound))
sp.Eq(integ, integ.doit())

"""
Explanation: The actual integral value
End of explanation
"""


actual_s1 = float(integ.doit())

error = actual_s1 - s1
print("Error ratio = {:.6f}%".format(abs(error / actual_s1) * 100.))

"""
Explanation: The error ratio
End of explanation
"""


fig, ax = plt.subplots(1, 1, figsize=(12, 8))

x_array = np.arange(x_lower_bound, x_upper_bound, 0.01)
y_array = f(x_array)

plt.axis([x_lower_bound, x_upper_bound, y_lower_bound, y_upper_bound])
plt.plot(random_points[:,0], random_points[:,1], ',k')
plt.plot(random_points_in_pos[:,0], random_points_in_pos[:,1], ',r')
plt.plot(random_points_in_neg[:,0], random_points_in_neg[:,1], ',r')
plt.hlines(y=0, xmin=x_lower_bound, xmax=x_upper_bound)
plt.plot(x_array, y_array, '-r', linewidth=2)

plt.show()

"""
Explanation: Graphical illustration
End of explanation
"""


