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
from matplotlib.widgets import Slider, Button   # widgets du module matplotlib !
from matplotlib import ticker
from scipy.special import comb

# 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, k, f):
    """  Probabilité d’avoir exactement p crues de fréquence f en n années : C(n,k)*f^k*(1-f)^(N-k)
    """
    return comb(n, k) * (f ** k) * ((1 - f) ** (n - k))


def afficher_probas(kmax, n, f):
    kmax = max(n, kmax)
    print(f" pour k variant de 0 à {kmax}"
        f"proba d'avoir exactement 'k crue(s)' supérieures à la crue de période de retour {1 / f:.0f} (= fréquence {f:.02f})  en {n} années ")

    for k in range(kmax+1):
        if k < 2:
            texte = "crue supérieure"
        else:
            texte = "crues supérieures"

        print(f"    exactement {k} {texte} :  {proba_crue(n, k, f):.2f)*100:.2f}%")

def calcul_courbes(n,f):
    cumul_au_plus = 0
    liste_probas_exactement=[]
    liste_probas_auplus = []
    k= 0
    while cumul_au_plus < 0.999:
        proba = proba_crue(n,k,f)
        liste_probas_exactement.append(proba)
        cumul_au_plus += proba
        liste_probas_auplus.append(cumul_au_plus)
        k+=1
    # on retourne k-1 car on a effectué le dernier calcul pour p-1
    return k-1, liste_probas_exactement, liste_probas_auplus

def stair_plot(x,y, couleur, label=None):
    x = [*x, x[-1]]  # duplicate last  to draw last line
    ax.hlines(y, xmin=x[:-1], xmax=x[1:], linewidth=3, color=couleur, label=label)
    # 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)

    k_dernier, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n,1/T)
    fig.suptitle(f"Démo : probabilité d'avoir k crues > T={T}ans en {n} années")

    courbe_exactement.set_data(range(k_dernier + 1), liste_probas_exactement)
    courbe_cumul.set_data(range(k_dernier + 1), liste_probas_auplus)

    ax.set_xlim(0, k_dernier + 1)

    if k_dernier > 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)

    k_dernier, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
    fig.suptitle(f"Démo : probabilité d'avoir k crues > T={T}ans en {n} années")

    ax.clear()
    stair_plot(range(k_dernier + 1), liste_probas_exactement, couleur="blue", label="exactement k crues")
    stair_plot(range(k_dernier + 1), liste_probas_auplus, couleur="sienna", label="k crues ou moins de k")

    ax.axhline(y=0, color='grey', ls=':')
    ax.axhline(y=1, c='grey', ls=':')
    ax.set_xlim(0, k_dernier + 1)
    ax.set_xlabel("nombre de crues k")
    ax.set_ylabel("probabilité")

    if k_dernier > 15:
        ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, k_dernier, 5)]))
        ax.xaxis.set_minor_locator(plt.FixedLocator([0.5 + i for i in range(0, k_dernier, 1)]))
    else:
        ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, k_dernier, 1)]))
        ax.xaxis.set_minor_locator(plt.NullLocator())
    ax.xaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
    ax.legend(loc='center right')
    fig.canvas.draw_idle()


#CORPS DU PROGRAMME
n = 100
T = 100
k_dernier, 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 k_dernier crues > T={T}ans en {n} années")

if TRACE_PAR_PLOT:
    courbe_exactement, = ax.plot(range(k_dernier + 1), liste_probas_exactement, label="exactement k crues",ls='None', marker='o', markersize=10)
    courbe_cumul, = ax.plot(range(k_dernier + 1), liste_probas_auplus, label="k crues ou moins de k", ls='None', marker='o', markersize=10)
    update_graphique = update_graphique_plot
else:
    stair_plot(range(k_dernier + 1),liste_probas_exactement, couleur="blue", label="exactement k crues")
    stair_plot(range(k_dernier + 1),liste_probas_auplus, couleur="sienna",  label="k crues ou moins de k")
    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 k")
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)

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(loc='center right')
ax.set_xlim(0,k_dernier+1)
ax.axhline(y=0,color='grey',ls=':')
ax.axhline(y=1, c='grey',ls=':')

if TRACE_PAR_PLOT :
    if k_dernier > 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,k_dernier,5)]))
        ax.xaxis.set_minor_locator(plt.FixedLocator([0.5+i for i in range(0,k_dernier,1)]))
    else:
        ax.xaxis.set_major_locator(plt.FixedLocator([0.5+i for i in range(0,k_dernier,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")