GenerateurCruesMaxAnnuelles.py 16.05 KiB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
"""
CPoulard  ; mai 2021
côté "métier" :
                montrer le principe des formules Tchégodaiev / Hazen utilisées en hydrologie (TD 1ere année ENTPE)  et la relation fréquence/période de retour
                démo du caractère aléatoire des crues :  on  constitue progressivement une série d'  "observations"  issues d'une loi de Gumbel
coté "Python" : mettre en oeuvre deux widgets de matplotlib (Checkbutton et Slider), faciles d'emploi mais au rendu pas très abouti ; définir les fonctions associées à des événements sur ces boutons ; manipuler les fonctions en tant qu'objet (fonction_en_cours = une des  3 fonctions disponibles...)
"""

from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.widgets import Slider, CheckButtons, Button   # widgets du module matplotlib !
import numpy as np
from numpy.random import gumbel


# FONCTIONS
#
# trois fonctions correspondant chacune à une manière de définir un échantillon de nb valeurs
def echantillon_uniforme(nb):
    return sorted(range(1, nb+1), reverse=True)

def echantillon_constant(nb):
    # non encore testé
    return [1]*nb

def echantillon_gumbel(nb, x0=196, gr=102):
    # non encore testé

    return gumbel(loc=x0, scale=gr, size=nb)

class SerieGumbel:
    T_Gumbel_plot = [1.1, 2, 5, 10, 20, 50, 100]

    def __init__(self, gr, x0):
        self.gr = gr
        self.x0 = x0
        self.serie_chrono = []
        self.serie_triee = []
        self.serie_couleurs = []

        self.QT10 = self.QdeT(10)
        self.QT50 = self.QdeT(50)
        self.QT100 = self.QdeT(100)
        self.QGumbel_plot = [self.QdeT(T) for T in SerieGumbel.T_Gumbel_plot]
        self.gr_echantillon = None
        self.x0_echantillon = None
        print("quantiles ", self.QT10, self.QT50, self.QT100)

    def attribution_couleur(self, Q):

        if Q < self.QT10:
            couleur = "black"
        elif Q < self.QT50:
            couleur = "green"
        elif Q < self.QT100:
            couleur = "blue"
        else:
            couleur = 'red'
        return couleur

    def QdeT(self, T):
        """Donne la valeur estimée pour une  période de retour donnée selon une loi de Gumbel
           Args:
               T (float) :  période de retour
               gr, x0 (float) :  paramètres d'échelle et de position de la loi de Gumbel
           Returns :
               QdeT(float) : valeur estimée
        """
        Fnondepassement = 1 - (1 / T)  # fréquence au non dépassement = 1 - f
        u = -np.log(-np.log(Fnondepassement))  # log x s'obtient par np.log(x)
        Q = self.x0 + self.gr * u
        return Q

    def TdeQ(self, Q):
        """Donne la valeur estimée pour une  période de retour donnée selon une loi de Gumbel
           Args:
               Q (float) : valeur
               gr, x0 (float) :  paramètres d'échelle et de position de la loi de Gumbel
           Returns :
               TdeQ(float) : période de retour estimée
        """
        u = (Q - self.x0) / self.gr  # log x s'obtient par math.log(x)
        Fnondepassement = np.exp(-1 * np.exp(-u))

        T = 1 / (1 - Fnondepassement)  # fréquence au non dépassement = 1 - f

        return T

    def ajout_obs(self, nb):
        nouvelles_obs = gumbel(loc=self.x0, scale=self.gr, size=nb)
        self.serie_chrono.extend(nouvelles_obs)
        self.serie_couleurs.extend([self.attribution_couleur(Q) for Q in nouvelles_obs])
        self.serie_triee.extend(nouvelles_obs)
        self.serie_triee.sort(reverse=True)

        self.gr_echantillon = 0.78 * np.std(self.serie_chrono) # gradex = paramètre d'échelle
        self.x0_echantillon = np.mean(self.serie_chrono) - 0.577 * self.gr_echantillon  # "ordonnée à l'origine" = paramètre de position


# Loi de Gumbel (max annuels)
#définition des lois d'ajustement
def AjustementParametresDeGumbel(moyenne, ecart_type):  # 2 arguments à calculer par ailleurs
    """Ajuste une loi de Gumbel, usuelle en hydrologie pour ajuster des valeurs max annuelles
    à partir "seulement" de la moyenne et l'écart-type de l'échantillon
    Elle renvoie ses 2 paramètres gr et x0 calculés par la méthode des moments"
    # ils sont donc moins sensibles aux points exceptionnels que la régression linéaire (distances au CARRE)
    Args:
        moyenne (float) :  moyenne de l'échantillon
        ecart_type (float) :  ecart-type de l'échantillon
    Returns :
        gr(float) : paramètre d'échelle, gradex
        x0(float) : paramètre de position

    """

    gr = 0.78 * ecart_type  # gradex = paramètre d'échelle
    x0 = moyenne - 0.577 * gr  # "ordonnée à l'origine" = paramètre de position

    return gr, x0


def QdeT_Gumbel(T, gr, x0):
    """Donne la valeur estimée pour une  période de retour donnée selon une loi de Gumbel
       Args:
           T (float) :  période de retour
           gr, x0 (float) :  paramètres d'échelle et de position de la loi de Gumbel
       Returns :
           QdeT(float) : valeur estimée
    """
    Fnondepassement = 1 - (1 / T)  # fréquence au non dépassement = 1 - f
    u = -np.log(-np.log(Fnondepassement))  # log x s'obtient par np.log(x)

    QdeT = x0 + gr * u

    return QdeT

def TdeQ_Gumbel(Q, gr, x0):
    """Donne la valeur estimée pour une  période de retour donnée selon une loi de Gumbel
       Args:
           Q (float) : valeur
           gr, x0 (float) :  paramètres d'échelle et de position de la loi de Gumbel
       Returns :
           TdeQ(float) : période de retour estimée
    """
    u = (Q - x0) / gr  # log x s'obtient par math.log(x)
    Fnondepassement = np.exp(-1 * np.exp(-u))

    TdeQ = 1 / (1 - Fnondepassement)   # fréquence au non dépassement = 1 - f

    return TdeQ


# fonction paramétrée calculant la fréquence par une fonction paramétrée par a et b
# puis 2 fonctions calculant la période de retour par la même fonction paramétrée et par une fonction avec valeurs fixes
# EVIDEMMENT IL Y A REDONDANCE pour l'exemple ; la première fonction suffirait !

def freq_emp_parametre(n_annees, a, b):
    # pour la formule de Chegodayev, a = 0.3 et b=0.4

    plotting_positions_freq = [ 1 - (((indice + 1) - a) / ((n_annees + b) )) for indice in range(n_annees)]
    # au lieu d'envoyer une liste brute pour l'affichage,
    # on crée une liste de freq formattées (:.2f = flottants 2 décimales) séparées par des "/"
    print("Fréquences empiriques : ", " / ".join([f"{freq:.2f}" for freq in plotting_positions_freq]))
    return plotting_positions_freq

def T_emp_parametre(n_annees, a , b):
    # T=1/freq de la fonction freq_emp_parametre, cas général de la précédente
    plotting_positions_param  = [((n_annees + b) / ((indice + 1) - a)) for indice in range(n_annees)]
    print(f"Périodes de retour empiriques, a={a:.1f}, b={b:.1f}, N={n_annees} : ", " / ".join([f"{freq:.2f}" for freq in plotting_positions_param]))
    return plotting_positions_param


def retracer_parametree(echantillon, a, b):
    # si switch :  au lieu de "updater", on retrace complètement en attendant se pouvoir updater des stems..
    global stemlines_emp, markerline_emp , baseline_emp

    if en_periode_de_retour:
        x_emp =  T_emp_parametre(len(echantillon), a, b)
    else:
        # freq = [ 1 - (1 / T) for T in T_emp]
        x_emp = freq_emp_parametre(len(echantillon), a, b)

    for item in [stemlines_emp, markerline_emp, baseline_emp]:
        item.remove()

    markerline_emp, stemlines_emp, baseline_emp = ax_T.stem(x_emp, serie_Gumbel.serie_triee, label=f"T ou f empirique, fonction de a et b",
                                        use_line_collection=True)

    plt.setp(stemlines_emp, color='purple')
    plt.setp(markerline_emp, color='cyan', markersize=2, markeredgecolor='blue', markeredgewidth=3)
    plt.setp(baseline_emp, visible=False)

    # on retrace
    fig_pp.canvas.draw_idle()

def retracer_Gumbel(serie_Gumbel):

    if en_periode_de_retour:
        courbe_Gumbel.set_xdata(SerieGumbel.T_Gumbel_plot)
    else:
        #freq = [ 1 - (1 / T) for T in T_emp]
        freq = [1 - (1/T) for T in SerieGumbel.T_Gumbel_plot]
        courbe_Gumbel.set_xdata(freq)

    # on retrace
    fig_pp.canvas.draw_idle()

def retracer_estim_Gumbel(serie_Gumbel):
    gr, x0 = AjustementParametresDeGumbel(np.mean(serie_Gumbel), np.std(serie_Gumbel))
    T_estim_Gumbel = [TdeQ_Gumbel(Q, gr, x0) for Q in serie_Gumbel]
    if en_periode_de_retour:

        courbe_estim_Gumbel.set_data(T_estim_Gumbel, serie_Gumbel)
    else:
        # freq = [ 1 - (1 / T) for T in T_emp]
        freq = [ 1 - (1 / T) for T in T_estim_Gumbel]
        courbe_estim_Gumbel.set_data(freq, serie_Gumbel)

    # on retrace
    fig_pp.canvas.draw_idle()


# https://matplotlib.org/stable/gallery/widgets/slider_demo.html
# https://stackoverflow.com/questions/13656387/can-i-make-matplotlib-sliders-more-discrete
def update_a(val):
    # val = valeur renvoyée par le widget
    global a
    a = val   # comme on a précisé "global a", changer a dans la fonction = changer a du prog principal
    retracer_parametree(serie_Gumbel.serie_triee, a, b)

def update_b(val):
    # val = valeur renvoyée par le widget
    global b
    b = val # comme on a précisé "global b", changer b dans la fonction = changer b du prog principal
    retracer_parametree(serie_Gumbel.serie_triee, a, b)

def create_series(event):
    # bouton "recréer une série de N obs" activé

    # comme on a précisé "global N_total, echantillon", sa modifications sera propagée jusque dans le prog principal
    global N_total, serie_Gumbel,  gr_est, x0_est, courbe_estim_Gumbel, stemlines_emp, markerline_emp, baseline_emp

    serie_Gumbel = SerieGumbel(x0=196, gr=102)
    serie_Gumbel.ajout_obs(N_serie)

    N_total = N_serie

    print("Première série d'obs : ", serie_Gumbel.serie_chrono)

    ax_chrono.clear()
    ax_T.clear()
    # LES OBS
    derniere_annee = 1 + len(serie_Gumbel.serie_chrono)
    ax_chrono.vlines(x=range(1, derniere_annee), ymin=0, ymax=serie_Gumbel.serie_chrono)
    ax_chrono.scatter(range(1, derniere_annee), serie_Gumbel.serie_chrono,
                                      c=serie_Gumbel.serie_couleurs, alpha=0.7,
                                      label=f"Q(t) observé sur {N_total} années, ", marker="*", s=60)

    gr_est, x0_est = serie_Gumbel.gr_echantillon, serie_Gumbel.x0_echantillon

    # LA COURBE DE GUMBEL, inchangée (sauf si tracé en fréquence...)
    if en_periode_de_retour:
        x_Gumbel = SerieGumbel.T_Gumbel_plot
        x_emp = T_emp_parametre(N_total, a, b)
        x_estimGumbel = [TdeQ_Gumbel(Q, gr_est, x0_est) for Q in serie_Gumbel.serie_triee]
    else:
        x_Gumbel = [1- (1/T) for T in SerieGumbel.T_Gumbel_plot]
        x_emp = [1- (1/T) for T in T_emp_parametre(N_total, a, b)]
        x_estimGumbel = [1 - (1 / TdeQ_Gumbel(Q, gr_est, x0_est)) for Q in serie_Gumbel.serie_triee]

    courbe_Gumbel, = ax_T.plot(x_Gumbel, serie_Gumbel.QGumbel_plot, alpha=0.7, color='orangered',
                               linewidth=5, solid_capstyle='round',
                               zorder=10, label="Q(T) du bassin versant", marker="None", ls='--', markersize=7)
    # courbe_estim_T_stem = markerline, stemlines, baseline (objet retournés par la fonction)
    markerline_emp, stemlines_emp, baseline_emp = ax_T.stem(x_emp, serie_Gumbel.serie_triee, label=f"T ou f empirique, fonction de a et b", use_line_collection=True)

    plt.setp(stemlines_emp, color='purple')
    plt.setp(markerline_emp, color ='cyan', markersize=5, markeredgecolor ='blue', markeredgewidth = 3)
    plt.setp(baseline_emp, visible=False)

    courbe_estim_Gumbel, = ax_T.plot(x_estimGumbel, serie_Gumbel.serie_triee, color='blue', alpha=0.7,
                                     label=" ajustement Gumbel sur obs",
                                     zorder=10, marker="x", markersize=7)
    ax_T.legend()
    fig_pp.canvas.draw_idle()

    return serie_Gumbel, courbe_Gumbel, gr_est, x0_est, courbe_estim_Gumbel, stemlines_emp, markerline_emp, baseline_emp


def update_series(event):
    # bouton "ajouter des obs" activé

    # comme on a précisé "global N_total, echantillon", sa modifications sera propagée jusque dans le prog principal
    global N_total
    N_total += N_serie
    serie_Gumbel.ajout_obs(N_serie)

    ax_chrono.clear()
    derniere_annee = 1 + len(serie_Gumbel.serie_chrono)
    ax_chrono.set_title(f"Qmax annuels observés sur {N_total} années")
    ax_chrono.vlines(x=range(1, derniere_annee), ymin=0, ymax=serie_Gumbel.serie_chrono, colors=serie_Gumbel.serie_couleurs) #, ls="--')
    ax_chrono.scatter(range(1, derniere_annee), serie_Gumbel.serie_chrono,
                                      c=serie_Gumbel.serie_couleurs, alpha=0.7, marker="*", s=20)
    retracer_estim_Gumbel(serie_Gumbel.serie_triee)
    retracer_parametree(serie_Gumbel.serie_triee, a, b)
    fig_pp.canvas.draw_idle()

def switch_freq(label):
    global en_periode_de_retour

    en_periode_de_retour = not en_periode_de_retour

    retracer_parametree(serie_Gumbel.serie_triee, a, b)
    retracer_Gumbel(serie_Gumbel.serie_triee)
    retracer_estim_Gumbel(serie_Gumbel.serie_triee)
    if en_periode_de_retour:
        ax_T.set_xlim(left=0, right= (N_total + b)/(1 - a) * 1.1)
        ax_T.set_xlabel("Période de retour, années")
        bouton_T2F.label.set_text('Période de retour => fréquence')
    else:
        # freq = [ 1 - (1 / T) for T in T_emp]
        ax_T.set_xlim(left=0, right=1.1)
        ax_T.set_xlabel("fréquence de non-dépassement")
        bouton_T2F.label.set_text('fréquence => Période de retour')

    ax_T.set_ylim(bottom=0, top=max(serie_Gumbel.serie_triee) * 1.1)

    fig_pp.canvas.draw_idle()

#plt.ion()

a_Tchego, b_Tchego = 0.3 , 0.4  # affectation de plusieurs valeurs simultanées par tuples
# initialisation des paramètres de départ (modifiables par curseur)
a = a_Tchego
b = b_Tchego

# on va constituer une série d'  "observations"  par ajouts de N_series obs à la fois
N_serie = 10

#rangs = echantillon_uniforme(N_serie)

en_periode_de_retour = True

fig_pp = plt.figure()
gs = gridspec.GridSpec(17, 2)
ax_chrono = plt.subplot(gs[0:6, :])
ax_bouton = plt.subplot(gs[7, 0])
ax_bouton_RAZ = plt.subplot(gs[7, 1])
ax_T = plt.subplot(gs[8:13, :])
ax_slide_a = plt.subplot(gs[16, 0])
ax_slide_b = plt.subplot(gs[16, 1])
ax_T2f = plt.subplot(gs[15, :])

# ax_chb = fig.add_subplot()
plt.subplots_adjust(wspace=1, hspace=0.5,left=0.1,top=0.90,right=0.9,bottom=0.1)
fig_pp.canvas.set_window_title("ScE - Hydrologie -Démo Gumbel")

fig_pp.suptitle(f"Démo Gumbel : tirages et leurs périodes de retour empiriques")

ax_T.set_xlabel("Période de retour (années)")


for ax_s in [ax_slide_a, ax_bouton_RAZ, ax_T2f, ax_slide_b, ax_bouton]:
    ax_s.xaxis.set_visible(False)
    ax_s.yaxis.set_visible(False)

"""
for pos in ['right', 'top', 'bottom', 'left']:
    ax_T2f.spines[pos].set_visible(False)
"""

# note : les variables après serie_Gumbel pourraient devenenir des attributs...
serie_Gumbel, courbe_Gumbel, gr_est, x0_est, courbe_estim_Gumbel, stemlines_emp, markerline_emp, baseline_emp = create_series('button_press_event')
ax_T.legend(bbox_to_anchor=(1.05, 0.05), loc='lower right')
print("Première série de max annuels créée")

# nom connu même hors de la fonction pour éviter le GC ?
# a
slider_a = Slider(
    ax_slide_a, "paramètre a ", valmin=0.01, valmax = 0.99, valfmt='%0.2f', valinit=a, color="green")

slider_a.on_changed(update_a)

# b
slider_b = Slider(
    ax_slide_b, "paramètre b ", valmin=0.01, valmax = 0.99, valfmt='%0.2f', valinit=b, color="blue")

slider_b.on_changed(update_b)

# bouton pour ajouter N_serie valeurs
couleur_bouton = 'linen'
bouton_obs = Button(ax_bouton, f" Ajouter {N_serie} obs ",  color=couleur_bouton)
bouton_obs.on_clicked(update_series)

# bouton pour ajouter N_serie valeurs
bouton_obs_RAZ = Button(ax_bouton_RAZ, f"Nouvelle série, {N_serie} obs ",  color=couleur_bouton)
bouton_obs_RAZ.on_clicked(create_series)

# switch T / freq
bouton_T2F = Button(ax_T2f, "axe des x en T => fréquence (de 0 à 1) ",  color=couleur_bouton)
bouton_T2F.on_clicked(switch_freq)
"""
for rect in chb_enT.rectangles:
    rect.set_width(0.1)
    rect.set_height(1)
    rect.xy = (0.05,0)
"""


#plt.tight_layout()
#fig_pp.savefig("ma_figure.png")
#fig_pp.canvas.draw()
plt.show()
#plt.ioff()