import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import scipy.integrate as scpi


def func_1(t):
    tmp = t*(1 - t)
    return tmp*np.sin(200*tmp)**2


def func(t, cMult):
    tmp = t*(1 - t)
    return tmp*np.sin(cMult*tmp)**2


def intMonteCarlo(f, xmin, xmax, n):
    w = xmax - xmin
    Y = [f(xmin + npr.rand()*w) for k in range(n)]
    S = sum(Y)
    I = (w/n)*S
    return I

def simpson(f, a, b, N):
    """
    Renvoie la valeur approchee de l'integrale de f sur [a, b]
    evaluee par la methode de Simpson.
    """
    h = (b - a)/N
    J = 0
    x0 = a
    for k in range(N):
        x1 = x0 + h/2
        x2 = x0 + h
        J += f(x0) + 4*f(x1) + f(x2)
        x0 = x2
    return h/6 * J


###
### tests
###

### graphe de la fonction
x = np.linspace(0, 1, 800)

fig0 = plt.figure()
ax0 = fig0.add_subplot(1, 1, 1, aspect="equal")
for k in [[20, "r"], [40, "g"], [80, "b"], [130, "k"]]:
    ax0.plot(x, func(x, k[0]), linewidth=.75, color=k[1], label=f"k = {k[0]}")
ax0.legend()
ax0.legend(bbox_to_anchor=(.6, -.17))
titre = "Les fonctions avec des constantes multiplicatives differentes"
ax0.set_title(titre)


### integrales
np.random.seed(0)
print("Les valeurs aleatoires sont fixees, i.e. toujours les memes "
      + "car on a utilise seed()\n")
n = 500
I = intMonteCarlo(func_1, 0, 1, n)
IE = scpi.quad(func_1, 0, 1)[0]
print(f"La valeur de l'integrale par la methode de Monte-Carlo avec n = {n} et"
      + f"\nla valeur obtenue avec quad() sont {I} et {IE}.")

print(simpson(func_1, 0, 1, 71))


N_list = [100, 200, 500, 1000, 2000, 3000, 5000, 10000, 20000, 50000]

e_list_MC = []
e_list_S = []
for n in N_list:
    I = intMonteCarlo(func_1, 0, 1, n)
    J = simpson(func_1, 0, 1, n)
    e_list_MC.append(abs(I - IE))
    e_list_S.append(abs(J - IE))


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)  # , aspect="equal")
ax.plot(N_list, e_list_MC, linewidth=1, marker="o", markersize=4, label='MC')
ax.plot(N_list, e_list_S, linewidth=1, marker="o", markersize=4, c='r',
        label='S')
ax.set_xscale('log')
ax.set_yscale('log')
titre = "L'erreur de la methode Monte-Carlo en coordonees logarithmiques"
ax.set_title(titre)
ax.legend()

plt.show()