Newer
Older
Poulard Christine
committed
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
"""
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)
if p < 2:
texte = "crue supérieure"
else:
texte = "crues supérieures"
print(
f"proba d'avoir exactement p {texte} à 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} {texte} : {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, 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)
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", label="exactement p crues")
stair_plot(range(p + 1), liste_probas_auplus, couleur="sienna", label="p crues ou moins de p")
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", label="exactement p crues")
stair_plot(range(p + 1),liste_probas_auplus, couleur="sienna", label="p crues ou moins de p")
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")