# -*- coding: utf-8 -*-
"""
Created on Sat Oct 12 10:08:13 2024

@author: pbert

Mode d'emploi :
    1) Définir la fonction f avec la grandeur dont on cherche à évaluer l'incertitude,
    avec autant de paramètres que nécessaire.
    Exemple: def f(x,y): return sin(x*y) pour évaluer l'incertitude sur sin(x*y) connaissant
    x, y et leur incertitudes-types u(x), u(y)...)
    2) Définir la liste Lparams, sous la forme [[x,u(x)],[y,u(y)],...] 
    avec les paramètres x,y,... dans le même ordre que dans la déclaration de f.
    3) Exécuter le programme et relever les résultats en adaptant les CS 
    et sans oublier l'unité. Contrôler la cohérence de l'histogramme.
"""

from math import *
import numpy as np
import matplotlib.pyplot as plt


# =============================================================================
# Evaluation d'incertitudes par la méthode de Monte Carlo
# =============================================================================

def eval_incert(f,Lparams):
    n=10000 #nombre de points aléatoires pour les échantillons
    L_alea=np.empty((n,len(Lparams)))
    Lres=np.empty(n)
    for k in range(len(Lparams)): #Génération des chantillons sur les paramètres
        L_alea[:,k]=np.random.normal(Lparams[k][0],Lparams[k][1],n)
    for _ in range(n):
        Lres[_]=f(*L_alea[_])
    return np.mean(Lres),np.std(Lres),Lres

# =============================================================================
# Définition de la focntion f et des valeurs d'entrées
# =============================================================================
def f(x,y):
    return exp(x/y)

Lparams=[[0.2,2e-3],[1,9e-2]]

# =============================================================================
# Affichage résultats (ne pas modifier)
# =============================================================================

res,u_res,Lres=eval_incert(f,Lparams)

print("f a pour valeur {:.6} (unité) et son incertitude-type a pour valeur {:.2} (unité)".format(res,u_res))

plt.hist(Lres,bins=50)
plt.show()

    
    




