#-*-coding:utf-8-*-
###########################################
########  Réalisation d'un graphique   ####
##  A partir de mesures expérimentales    #
###########################################

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress
from statistics import stdev, mean
plt.rcParams['font.size'] = 10

plt.close("all") # pour fermer les fenêtres graphiques


## Récupération du fichiers de points

# Mettre le fichier de données dans le même dossier

def recup(nomfichier, colonnes, sep = ";"):
    """
    nomfichier : fichier texte .txt ou .csv
    colonnes : type int, nbre de colonnes à récupérer
    sep : type str, séparateur de colonnes, ";" par défaut
    RETOUR : type array, tableau de données.
    """
    L = [ [] for _ in range(colonnes)] #initialisaition de la liste de données
    with open(nomfichier, "r") as contenu:
        for ligne in contenu:
            try : #pour éviter les entêtes
                ligne = ligne.replace(',','.') #  transforme les , en .
                ligne = ligne.split(sep) #sépare les valeurs
                for i in range(colonnes):
                    L[i].append(float(ligne[i]) )#remplissage et conversion en flottant
            except: ValueError
    return np.array(L) # conversion en tableau


#pour affichage


##Points expérimentaux factices
nomfichier ='TP.csv'
Donnees = recup(nomfichier, 2, sep = ";")

x =  Donnees[0] #abscisses
y = Donnees[1] #ordonnées


## Valeur de l'incertitudes de mesures : +/-0.2 en x, +/- 0.4 en y
u_x = [.001 for k in range(len(x))]
u_y = [.00001 for k in range(len(y))]

## Régression linéaire
(a,b,rho,_,_)=linregress(x,y)
print("Coefficient de correlation : ",rho)
print("Equation de la droite : ")
print("y =",a," x +",b)

## représentation graphique
plt.figure(figsize=(10,8))
plt.subplot(211)
plt.grid()
plt.plot(x,a*x+b,label="modèle",linewidth=3)
plt.errorbar(x,y,xerr=u_x,yerr=u_y,label="expérience",capsize = 5,fmt="go",linewidth=2)
plt.legend(loc="best",fontsize=12)
plt.xlabel(" X (unite)")
plt.ylabel(" Y (unite)")
plt.title(" Mon Titre ")



## Méthode de monté-Carlo
N =  1000 # nombre de tirages
liste_a,liste_b = [],[]
n = len(x)
for _ in range(N):
    n = len(x)
    x_alea = x + np.random.uniform(-u_x[0],u_x[0],n)
    y_alea = y + np.random.uniform(-u_y[0],u_y[0],n)
    (a_alea,b_alea,rho,_,_)=linregress(x_alea,y_alea)
    liste_a.append(a_alea)
    liste_b.append(b_alea)



plt.subplot(223)
_,bins,_ = plt.hist(liste_a, 50,
                            density = True,
                            color ='green',
                            alpha = 0.3)
# Répartition gaussienne
sigma, mu = stdev(liste_a), mean(liste_a)
y = ((1 / (np.sqrt(2 * np.pi) * sigma)) *
     np.exp(-0.5 * (1 / sigma * (bins - mu))**2))
plt.plot(bins, y, '--', color ='black')
# Affichage moyenne
plt.plot([mu,mu],[0,max(y)],'--r')
texte = "$\mu$ = {:.2e} ".format(mu)
plt.text(mu+sigma/10,max(y)/20,texte,color='r')
# Affichage écart type
plt.arrow(mu,max(y)*.55,sigma,0, head_width = max(y)/50,
            head_length = sigma/20,fc='k', ec='k')
texte = "$\sigma =  {:.0e} $".format(sigma)
plt.text(mu+sigma/10,max(y)*.56,texte)
# noms des axes et titre
plt.xlabel('a')
plt.ylabel('Fréquence')
plt.title('Coefficient a par MC')



plt.subplot(224)
_,bins,_ = plt.hist(liste_b, 50,
                            density = True,
                            color ='green',
                            alpha = 0.3)
# Répartition gaussienne
sigma, mu = stdev(liste_b), mean(liste_b)
y = ((1 / (np.sqrt(2 * np.pi) * sigma)) *
     np.exp(-0.5 * (1 / sigma * (bins - mu))**2))
plt.plot(bins, y, '--', color ='black')
# Affichage moyenne
plt.plot([mu,mu],[0,max(y)],'--r')
texte = "$\mu$ = {:.2e} ".format(mu)
plt.text(mu+sigma/10,max(y)/20,texte,color='r')
# Affichage écart type
plt.arrow(mu,max(y)*.55,sigma,0, head_width = max(y)/50,
            head_length = sigma/20,fc='k', ec='k')
texte = "$\sigma =  {:.0e} $".format(sigma)
plt.text(mu+sigma/10,max(y)*.56,texte)
# noms des axes et titre
plt.xlabel('b')
plt.ylabel('Fréquence')
plt.title('Coefficient b par MC')


plt.show()