Upload New File : version de ProbaCruesMaxAn_SurNannees avec un autre type de...

Upload New File : version de ProbaCruesMaxAn_SurNannees avec un autre type de représentation possible, en escaliers

ProbaCruesMaxAn_SurNannees_stairs.py

+"""
+CPoulard  ; mai 2021
+côté "métier" :
+                appliquer la formule de calcul de proba du nombre de crues s en hydrologie (TD 1ere année ENTPE)
+
+coté "Python" : mettre en oeuvre deux widgets de matplotlib (Button et Slider), faciles d'emploi mais au rendu pas très abouti ;
+              définir les fonctions associées à des événements sur ces boutons
+"""
+
+from matplotlib import pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.widgets import Slider, Button   # widgets du module matplotlib !
+from matplotlib import ticker
+from scipy.special import comb
+import numpy as np
+
+# CONSTANTE
+
+# tracé par plot : on utilise "plot" qui est plus facile mais moins adapté à une fonction 'discrète"
+# sinon on utilise un tracé avec "hlines et fillbetween"
+# en attendant "stairs" (https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.stairs.html)
+TRACE_PAR_PLOT = False # False
+
+# FONCTIONS
+#
+# trois fonctions correspondant chacune à une manière de définir un échantillon de nb valeurs
+def proba_crue(n, p, f):
+    """  Probabilité d’avoir exactement p crues de fréquence f en n années : C(n,p)*f^p*(1-f)^(N-p)
+    """
+    return comb(n, p) * (f ** p) * ((1 - f) ** (n - p))
+
+
+def afficher_probas(nb, n, p, f):
+    nb = max(nb, p)
+    print(
+        f"proba d'avoir exactement p crue(s) supérieure(s) à la crue de période de retour {1 / f:.0f} (= fréquence {f:.02f] }en {n} années ")
+    for p in range(nb):
+        print(f"    exactement {p} crue(s) :  {proba_crue(n, p, f):.2f)*100:.2f}%")
+
+def calcul_courbes(n,f):
+    cumul_au_plus = 0
+    liste_probas_exactement=[]
+    liste_probas_auplus = []
+    p= 0
+    while cumul_au_plus < 0.999:
+        proba = proba_crue(n,p,f)
+        liste_probas_exactement.append(proba)
+        cumul_au_plus += proba
+        liste_probas_auplus.append(cumul_au_plus)
+        p+=1
+
+    return p-1, liste_probas_exactement, liste_probas_auplus
+
+def stair_plot(x,y, couleur):
+    x = [*x, x[-1]]  # duplicate last  to draw last line
+    ax.hlines(y, xmin=x[:-1], xmax=x[1:], linewidth=3, color=couleur)
+    # fill under line  :  y with duplicated last value to draw last step as well ; hide edge
+    ax.fill_between(x, [*y, y[-1]],
+                    step='post', facecolor=couleur, edgecolor=None, alpha=0.5)
+
+def update_graphique_plot(val):
+    global n, T
+    n = int(slider_n.val)
+    T = int(slider_T.val)
+
+    p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n,1/T)
+    fig.suptitle(f"Démo : probabilité d'avoir p crues > T={T}ans en {n} années")
+
+        courbe_exactement.set_data(range(p + 1), liste_probas_exactement)
+        courbe_cumul.set_data(range(p + 1), liste_probas_auplus)
+
+    ax.set_xlim(0, p + 1)
+
+    if p > 15:
+        ax.xaxis.set_major_locator(plt.MultipleLocator(5))
+        ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
+    else:
+        ax.xaxis.set_major_locator(plt.MultipleLocator(1))
+        ax.xaxis.set_minor_locator(plt.NullLocator())
+    ax.xaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
+    fig.canvas.draw_idle()
+
+def update_graphique_stairs(val):
+    global n, T
+    n = int(slider_n.val)
+    T = int(slider_T.val)
+
+    p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
+    fig.suptitle(f"Démo : probabilité d'avoir p crues > T={T}ans en {n} années")
+
+    ax.clear()
+    stair_plot(range(p + 1), liste_probas_exactement, couleur="blue")
+    stair_plot(range(p + 1), liste_probas_auplus, couleur="sienna")
+
+    ax.set_xlim(0, p + 1)
+    if p > 15:
+        ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 5)]))
+        ax.xaxis.set_minor_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 1)]))
+    else:
+        ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 1)]))
+        ax.xaxis.set_minor_locator(plt.NullLocator())
+    ax.xaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
+    fig.canvas.draw_idle()
+
+
+#CORPS DU PROGRAMME
+n = 100
+T = 100
+p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
+
+
+fig, (ax, ax_espace, ax_sn, ax_sT) = plt.subplots(nrows=4, ncols=1, gridspec_kw={'height_ratios':[6,1, 1,1]}, sharex=False)
+plt.subplots_adjust(left=0.2, bottom=None, right=None, top=None, wspace=None, hspace=0.1)
+fig.canvas.set_window_title("ScE - Hydrologie - Démo probas crue sur N années")
+fig.suptitle(f"Démo : probabilité d'avoir p crues > T={T}ans en {n} années")
+
+if TRACE_PAR_PLOT:
+    courbe_exactement, = ax.plot(range(p + 1), liste_probas_exactement, label="exactement p crues",ls='None', marker='o', markersize=10)
+    courbe_cumul, = ax.plot(range(p + 1), liste_probas_auplus, label="p crues ou moins de p", ls='None', marker='o', markersize=10)
+    update_graphique = update_graphique_plot
+else:
+    stair_plot(range(p + 1),liste_probas_exactement, couleur="blue")
+    stair_plot(range(p + 1),liste_probas_auplus, couleur="sienna")
+    update_graphique = update_graphique_stairs
+
+for ax_s in [ax_espace, ax_sn, ax_sT]:
+    ax_s.xaxis.set_visible(False)
+    ax_s.yaxis.set_visible(False)
+
+ax_espace.patch.set_alpha(0.01)  #sinon cet axe cache le titre de l'axe des x au-dessus !
+
+for pos in ['right', 'top', 'bottom', 'left']:
+            ax_espace.spines[pos].set_visible(False)
+
+ax.set_xlabel("nombre de crues p")
+ax.set_ylabel("probabilité")
+
+# nom connu même hors de la fonction pour éviter le GC ?
+# nombre d'années d'observation
+slider_n = Slider(
+    ax_sn, "nombre d'années d'observation ", valmin=1, valmax = 1000, valfmt='%0.0f', valinit=n, color="green")
+
+slider_n.on_changed(update_graphique)
+
+# T, période de retour
+# valeurs_possibles = np.array([2, 5, 10, 20, 25, 50, 100, 200])
+
+slider_T = Slider(
+    ax_sT, "période de retour T  ", 2, 1000, valinit=T, valstep=1, valfmt='%0.0f', color="blue")
+
+slider_T.on_changed(update_graphique)
+
+ax.legend()
+ax.set_xlim(0,p+1)
+if TRACE_PAR_PLOT :
+    if p > 15:
+        ax.xaxis.set_major_locator(plt.MultipleLocator(5))
+        ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
+    else:
+        ax.xaxis.set_major_locator(plt.MultipleLocator(1))
+        ax.xaxis.set_minor_locator(plt.NullLocator())
+else:
+    if p > 15:
+        ax.xaxis.set_major_locator(plt.FixedLocator([0.5+i for i in range(0,p,5)]))
+        ax.xaxis.set_minor_locator(plt.FixedLocator([0.5+i for i in range(0,p,1)]))
+    else:
+        ax.xaxis.set_major_locator(plt.FixedLocator([0.5+i for i in range(0,p,1)]))
+        ax.xaxis.set_minor_locator(plt.NullLocator())
+
+ax.xaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
+# plt.tight_layout() : cela raccourcit l'axe qui porte le graphique...
+plt.show()
+
+
+
+
+
+
+print("Premiers calculs, valeurs par défaut")
+
+