## MÉTHODE MONTE-CARLO - Relation Z = f(X,Y)
## Armel MARTIN - Septembre 2021
## distribué sous licence CC-by-nc-sa
import numpy as np
import matplotlib.pyplot as plt
import random

#### MESURE D'UNE DISTANCE FOCALE PAR LA RELATION DE DESCARTES
## Fabrication des fausses mesures
##fp = 10. # focale cm
##ux = 0.3 # précision pointé viseur
##OA = -np.array([15.,20.,25.,30.,35.,40.,45.,50.])
##OAp = 1/( 1/fp + 1/OA )
##uOA = np.sqrt(2/3)*ux
##uOAp = uOA
##X = 1/OA
##Y = 1/OAp
##uX = uOA/OA**2
##uY = uOAp/OAp**2
##Ym = np.copy(Y)
##for i in range(len(OA)):
##    Ym[i] = Y[i] + random.uniform(-np.sqrt(3)*uY[i],np.sqrt(3)*uY[i])
##plt.plot(X,Ym,'ok-')
##plt.plot(X,Y,'r-')
##plt.xlabel('X')
##plt.ylabel('Y')
##plt.grid()
##titre = "Vérification loi Y = f(X) = aX + b"
##plt.title(titre)
##plt.show()

# Mesures
OA = -np.array([15.,20.,25.,30.,35.,40.,45.,50.])
OAp = np.array([29.9, 19.9, 16.7, 15.3, 14.3, 13.0 , 12.8, 12.9])
X = 1/OA
Y = 1/OAp

a,b = np.polyfit(X,Y,1)
# Valeurs calculées par le modèle
Y_reg = b + a*X
# Analyse des composantes de variance
SST = np.mean((Y - np.mean(Y))**2) # variance totale
SSE = np.mean((Y_reg - np.mean(Y))**2) # variance expliquée (modèle)
# Coeff de détermination (R²)
R2 = SSE / SST
# Affichage des résultats en console
print('b=',b,' ; a=',a,' ; R2=',R2)
## Fabriquer l'équation de la courbe (type string)
equation = 'Y = '+format(a,'.3e')+' X + '+format(b,'.3e')

plt.figure(figsize=[5,3])
plt.plot(X,Y,'s',color='black')
plt.plot(X,Y_reg,'-k')
plt.xlabel('X')
plt.ylabel('Y')
plt.title("Vérification loi Y = f(X) = X + 1/fp")
plt.annotate(equation,(-0.065,0.07))
plt.grid()
plt.tight_layout()
plt.savefig('reglin.png')
plt.show()
#plt.clf()

## Incertitudes
Delta_x = 0.3 # précision pointé viseur
uOA = np.sqrt(2/3)*Delta_x
uOAp = uOA
uX = uOA/OA**2
uY = uOAp/OAp**2

## Génération aléatoire d'une série de séries de mesures
N = 100000 # nombre de simulations à exécuter.

## En utilisant directement des vecteurs Numpy
veca,vecb = np.zeros(N),np.zeros(N)
Xmod,Ymod = np.copy(X),np.copy(Y)
for i in range(0,N):
    for k in range(len(X)):
        Xmod[k] = random.uniform(X[k]-uX[k],X[k]+uX[k])
        Ymod[k] = random.uniform(Y[k]-uY[k],Y[k]+uY[k])
    veca[i],vecb[i] = np.polyfit(Xmod,Ymod,1)

#### Statistique sur Z
am = np.mean(veca)
bm = np.mean(vecb)
ua = np.std(veca,ddof=1) # ddof=1 => non biaisé
ub = np.std(vecb,ddof=1) # ddof=1 => non biaisé
print('am =',am,' ; u(a) =',ua)
print('bm =',bm,' ; u(b) =',ub)

plt.figure(1,figsize=[5,3])
plt.hist(veca,bins='rice') #,color='grey')
#plt.xlim(0,1.1*max(10**6*vecY))
plt.xlabel('a')# (10^-6 s)')
plt.ylabel('frequence')
plt.grid()
titre = "Coeff dir a ="+ format(am,'.4f') + ' +- '+ format(ua,'.4f') 
plt.title(titre)
plt.tight_layout()
plt.savefig('reglin_hist_a.png')
plt.show()

plt.figure(1,figsize=[5,3])
plt.hist(vecb,bins='rice') #,color='grey')
#plt.xlim(0,1.1*max(10**6*vecY))
plt.xlabel('b')# (10^-6 s)')
plt.ylabel('frequence')
plt.grid()
titre = "Ordonnée à l'origine b ="+ format(bm,'.4f') + ' +- '+ format(ub,'.4f') 
plt.title(titre)
plt.tight_layout()
plt.savefig('reglin_hist_b.png')
plt.show()


#plt.clf()
##
##plt.figure(2,figsize=[5,3])
##plt.hist(1e6*vecZ,bins='rice')
###plt.hist(vecZ,range=(0,100000),bins=100)
##plt.xlim(0,1.1*max(1e6*vecZ))
##plt.xlabel('10^6*Z')
##plt.ylabel('frequence')
##plt.grid()
##titre = "Distribution de Z=f(X,Y)=X/Y"
##plt.title(titre)
##resultat_MC = 'MC: Moyenne ='+format(1e6*Zm,'.3f')+'; Incert. relat. ='+format(100*uZ/Zm,'.1f')+'%'
##resultat_diff = 'Diff: Moyenne ='+format(1e6*f(xm,ym),'.3f')+'; Incert. relat. ='+ \
##               format(100*np.sqrt(uX**2/xm**2+uY**2/ym**2),'.1f')+'%'
##plt.annotate(resultat_MC,(0,1300))#(40850,700))
##plt.annotate(resultat_diff,(0,900))#(40850,500))
##plt.tight_layout()
##plt.savefig('hist_f(X,Y)_1.png')
##plt.show()
###plt.clf()
##