1"""
2Created on 20 dec. 2010
3
4This file illustrated the use of the workbench for a model
5specified in Python itself. The example is based on `Pruyt & Hamarat <https://www.systemdynamics.org/conferences/2010/proceed/papers/P1253.pdf>`_.
6For comparison, run both this model and the flu_vensim_no_policy_example.py and
7compare the results.
8
9
10.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
11 chamarat <c.hamarat (at) tudelft (dot) nl>
12
13"""
14
15import matplotlib.pyplot as plt
16import numpy as np
17from numpy import sin, min, exp
18
19from ema_workbench import Model, RealParameter, TimeSeriesOutcome, perform_experiments, ema_logging
20from ema_workbench import MultiprocessingEvaluator, SequentialEvaluator
21from ema_workbench.analysis import lines, Density
22
23# =============================================================================
24#
25# the model itself
26#
27# =============================================================================
28
29FINAL_TIME = 48
30INITIAL_TIME = 0
31TIME_STEP = 0.0078125
32
33switch_regions = 1.0
34switch_immunity = 1.0
35switch_deaths = 1.0
36switch_immunity_cap = 1.0
37
38
39def LookupFunctionX(variable, start, end, step, skew, growth, v=0.5):
40 return start + ((end - start) / ((1 + skew * exp(-growth * (variable - step))) ** (1 / v)))
41
42
43def flu_model(
44 x11=0,
45 x12=0,
46 x21=0,
47 x22=0,
48 x31=0,
49 x32=0,
50 x41=0,
51 x51=0,
52 x52=0,
53 x61=0,
54 x62=0,
55 x81=0,
56 x82=0,
57 x91=0,
58 x92=0,
59 x101=0,
60 x102=0,
61):
62 # Assigning initial values
63 additional_seasonal_immune_population_fraction_R1 = float(x11)
64 additional_seasonal_immune_population_fraction_R2 = float(x12)
65
66 fatality_rate_region_1 = float(x21)
67 fatality_rate_region_2 = float(x22)
68
69 initial_immune_fraction_of_the_population_of_region_1 = float(x31)
70 initial_immune_fraction_of_the_population_of_region_2 = float(x32)
71
72 normal_interregional_contact_rate = float(x41)
73 interregional_contact_rate = switch_regions * normal_interregional_contact_rate
74
75 permanent_immune_population_fraction_R1 = float(x51)
76 permanent_immune_population_fraction_R2 = float(x52)
77
78 recovery_time_region_1 = float(x61)
79 recovery_time_region_2 = float(x62)
80
81 susceptible_to_immune_population_delay_time_region_1 = 1
82 susceptible_to_immune_population_delay_time_region_2 = 1
83
84 root_contact_rate_region_1 = float(x81)
85 root_contact_rate_region_2 = float(x82)
86
87 infection_rate_region_1 = float(x91)
88 infection_rate_region_2 = float(x92)
89
90 normal_contact_rate_region_1 = float(x101)
91 normal_contact_rate_region_2 = float(x102)
92
93 ######
94 susceptible_to_immune_population_flow_region_1 = 0.0
95 susceptible_to_immune_population_flow_region_2 = 0.0
96 ######
97
98 initial_value_population_region_1 = 6.0 * 10**8
99 initial_value_population_region_2 = 3.0 * 10**9
100
101 initial_value_infected_population_region_1 = 10.0
102 initial_value_infected_population_region_2 = 10.0
103
104 initial_value_immune_population_region_1 = (
105 switch_immunity
106 * initial_immune_fraction_of_the_population_of_region_1
107 * initial_value_population_region_1
108 )
109 initial_value_immune_population_region_2 = (
110 switch_immunity
111 * initial_immune_fraction_of_the_population_of_region_2
112 * initial_value_population_region_2
113 )
114
115 initial_value_susceptible_population_region_1 = (
116 initial_value_population_region_1 - initial_value_immune_population_region_1
117 )
118 initial_value_susceptible_population_region_2 = (
119 initial_value_population_region_2 - initial_value_immune_population_region_2
120 )
121
122 recovered_population_region_1 = 0.0
123 recovered_population_region_2 = 0.0
124
125 infected_population_region_1 = initial_value_infected_population_region_1
126 infected_population_region_2 = initial_value_infected_population_region_2
127
128 susceptible_population_region_1 = initial_value_susceptible_population_region_1
129 susceptible_population_region_2 = initial_value_susceptible_population_region_2
130
131 immune_population_region_1 = initial_value_immune_population_region_1
132 immune_population_region_2 = initial_value_immune_population_region_2
133
134 deceased_population_region_1 = [0.0]
135 deceased_population_region_2 = [0.0]
136 runTime = [INITIAL_TIME]
137
138 # --End of Initialization--
139
140 Max_infected = 0.0
141
142 for time in range(int(INITIAL_TIME / TIME_STEP), int(FINAL_TIME / TIME_STEP)):
143 runTime.append(runTime[-1] + TIME_STEP)
144 total_population_region_1 = (
145 infected_population_region_1
146 + recovered_population_region_1
147 + susceptible_population_region_1
148 + immune_population_region_1
149 )
150 total_population_region_2 = (
151 infected_population_region_2
152 + recovered_population_region_2
153 + susceptible_population_region_2
154 + immune_population_region_2
155 )
156
157 infected_population_region_1 = max(0, infected_population_region_1)
158 infected_population_region_2 = max(0, infected_population_region_2)
159
160 infected_fraction_region_1 = infected_population_region_1 / total_population_region_1
161 infected_fraction_region_2 = infected_population_region_2 / total_population_region_2
162
163 impact_infected_population_on_contact_rate_region_1 = 1 - (
164 infected_fraction_region_1 ** (1 / root_contact_rate_region_1)
165 )
166 impact_infected_population_on_contact_rate_region_2 = 1 - (
167 infected_fraction_region_2 ** (1 / root_contact_rate_region_2)
168 )
169
170 # if ((time*TIME_STEP) >= 4) and ((time*TIME_STEP)<=10):
171 # normal_contact_rate_region_1 = float(x101)*(1 - 0.5)
172 # else:normal_contact_rate_region_1 = float(x101)
173
174 normal_contact_rate_region_1 = float(x101) * (
175 1 - LookupFunctionX(infected_fraction_region_1, 0, 1, 0.15, 0.75, 15)
176 )
177
178 contact_rate_region_1 = (
179 normal_contact_rate_region_1 * impact_infected_population_on_contact_rate_region_1
180 )
181 contact_rate_region_2 = (
182 normal_contact_rate_region_2 * impact_infected_population_on_contact_rate_region_2
183 )
184
185 recoveries_region_1 = (
186 (1 - (fatality_rate_region_1 * switch_deaths))
187 * infected_population_region_1
188 / recovery_time_region_1
189 )
190 recoveries_region_2 = (
191 (1 - (fatality_rate_region_2 * switch_deaths))
192 * infected_population_region_2
193 / recovery_time_region_2
194 )
195
196 flu_deaths_region_1 = (
197 fatality_rate_region_1
198 * switch_deaths
199 * infected_population_region_1
200 / recovery_time_region_1
201 )
202 flu_deaths_region_2 = (
203 fatality_rate_region_2
204 * switch_deaths
205 * infected_population_region_2
206 / recovery_time_region_2
207 )
208
209 infections_region_1 = (
210 susceptible_population_region_1
211 * contact_rate_region_1
212 * infection_rate_region_1
213 * infected_fraction_region_1
214 ) + (
215 susceptible_population_region_1
216 * interregional_contact_rate
217 * infection_rate_region_1
218 * infected_fraction_region_2
219 )
220 infections_region_2 = (
221 susceptible_population_region_2
222 * contact_rate_region_2
223 * infection_rate_region_2
224 * infected_fraction_region_2
225 ) + (
226 susceptible_population_region_2
227 * interregional_contact_rate
228 * infection_rate_region_2
229 * infected_fraction_region_1
230 )
231
232 infected_population_region_1_NEXT = infected_population_region_1 + (
233 TIME_STEP * (infections_region_1 - flu_deaths_region_1 - recoveries_region_1)
234 )
235 infected_population_region_2_NEXT = infected_population_region_2 + (
236 TIME_STEP * (infections_region_2 - flu_deaths_region_2 - recoveries_region_2)
237 )
238
239 if infected_population_region_1_NEXT < 0 or infected_population_region_2_NEXT < 0:
240 pass
241
242 recovered_population_region_1_NEXT = recovered_population_region_1 + (
243 TIME_STEP * recoveries_region_1
244 )
245 recovered_population_region_2_NEXT = recovered_population_region_2 + (
246 TIME_STEP * recoveries_region_2
247 )
248
249 if fatality_rate_region_1 >= 0.025:
250 qw = 1.0
251 elif fatality_rate_region_1 >= 0.01:
252 qw = 0.8
253 elif fatality_rate_region_1 >= 0.001:
254 qw = 0.6
255 elif fatality_rate_region_1 >= 0.0001:
256 qw = 0.4
257 else:
258 qw = 0.2
259
260 if (time * TIME_STEP) <= 10:
261 normal_immune_population_fraction_region_1 = (
262 additional_seasonal_immune_population_fraction_R1 / 2
263 ) * sin(4.5 + (time * TIME_STEP / 2)) + (
264 (
265 (2 * permanent_immune_population_fraction_R1)
266 + additional_seasonal_immune_population_fraction_R1
267 )
268 / 2
269 )
270 else:
271 normal_immune_population_fraction_region_1 = max(
272 (
273 float(qw),
274 (additional_seasonal_immune_population_fraction_R1 / 2)
275 * sin(4.5 + (time * TIME_STEP / 2))
276 + (
277 (
278 (2 * permanent_immune_population_fraction_R1)
279 + additional_seasonal_immune_population_fraction_R1
280 )
281 / 2
282 ),
283 )
284 )
285
286 normal_immune_population_fraction_region_2 = switch_immunity_cap * min(
287 (
288 (
289 sin((time * TIME_STEP / 2) + 1.5)
290 * additional_seasonal_immune_population_fraction_R2
291 / 2
292 )
293 + (
294 (
295 (2 * permanent_immune_population_fraction_R2)
296 + additional_seasonal_immune_population_fraction_R2
297 )
298 / 2
299 ),
300 (
301 permanent_immune_population_fraction_R1
302 + additional_seasonal_immune_population_fraction_R1
303 ),
304 ),
305 ) + (
306 (1 - switch_immunity_cap)
307 * (
308 (
309 sin((time * TIME_STEP / 2) + 1.5)
310 * additional_seasonal_immune_population_fraction_R2
311 / 2
312 )
313 + (
314 (
315 (2 * permanent_immune_population_fraction_R2)
316 + additional_seasonal_immune_population_fraction_R2
317 )
318 / 2
319 )
320 )
321 )
322
323 normal_immune_population_region_1 = (
324 normal_immune_population_fraction_region_1 * total_population_region_1
325 )
326 normal_immune_population_region_2 = (
327 normal_immune_population_fraction_region_2 * total_population_region_2
328 )
329
330 if switch_immunity == 1:
331 susminreg1_1 = (
332 normal_immune_population_region_1 - immune_population_region_1
333 ) / susceptible_to_immune_population_delay_time_region_1
334 susminreg1_2 = (
335 susceptible_population_region_1
336 / susceptible_to_immune_population_delay_time_region_1
337 )
338 susmaxreg1 = -(
339 immune_population_region_1 / susceptible_to_immune_population_delay_time_region_1
340 )
341 if (susmaxreg1 >= susminreg1_1) or (susmaxreg1 >= susminreg1_2):
342 susceptible_to_immune_population_flow_region_1 = susmaxreg1
343 elif (susminreg1_1 < susminreg1_2) and (susminreg1_1 > susmaxreg1):
344 susceptible_to_immune_population_flow_region_1 = susminreg1_1
345 elif (susminreg1_2 < susminreg1_1) and (susminreg1_2 > susmaxreg1):
346 susceptible_to_immune_population_flow_region_1 = susminreg1_2
347 else:
348 susceptible_to_immune_population_flow_region_1 = 0
349
350 if switch_immunity == 1:
351 susminreg2_1 = (
352 normal_immune_population_region_2 - immune_population_region_2
353 ) / susceptible_to_immune_population_delay_time_region_2
354 susminreg2_2 = (
355 susceptible_population_region_2
356 / susceptible_to_immune_population_delay_time_region_2
357 )
358 susmaxreg2 = -(
359 immune_population_region_2 / susceptible_to_immune_population_delay_time_region_2
360 )
361 if (susmaxreg2 >= susminreg2_1) or (susmaxreg2 >= susminreg2_2):
362 susceptible_to_immune_population_flow_region_2 = susmaxreg2
363 elif (susminreg2_1 < susminreg2_2) and (susminreg2_1 > susmaxreg2):
364 susceptible_to_immune_population_flow_region_2 = susminreg2_1
365 elif (susminreg2_2 < susminreg2_1) and (susminreg2_2 > susmaxreg2):
366 susceptible_to_immune_population_flow_region_2 = susminreg2_2
367 else:
368 susceptible_to_immune_population_flow_region_2 = 0
369
370 susceptible_population_region_1_NEXT = susceptible_population_region_1 - (
371 TIME_STEP * (infections_region_1 + susceptible_to_immune_population_flow_region_1)
372 )
373 susceptible_population_region_2_NEXT = susceptible_population_region_2 - (
374 TIME_STEP * (infections_region_2 + susceptible_to_immune_population_flow_region_2)
375 )
376
377 immune_population_region_1_NEXT = immune_population_region_1 + (
378 TIME_STEP * susceptible_to_immune_population_flow_region_1
379 )
380 immune_population_region_2_NEXT = immune_population_region_2 + (
381 TIME_STEP * susceptible_to_immune_population_flow_region_2
382 )
383
384 deceased_population_region_1_NEXT = deceased_population_region_1[-1] + (
385 TIME_STEP * flu_deaths_region_1
386 )
387 deceased_population_region_2_NEXT = deceased_population_region_2[-1] + (
388 TIME_STEP * flu_deaths_region_2
389 )
390
391 # Updating integral values
392 if Max_infected < (
393 infected_population_region_1_NEXT
394 / (
395 infected_population_region_1_NEXT
396 + recovered_population_region_1_NEXT
397 + susceptible_population_region_1_NEXT
398 + immune_population_region_1_NEXT
399 )
400 ):
401 Max_infected = infected_population_region_1_NEXT / (
402 infected_population_region_1_NEXT
403 + recovered_population_region_1_NEXT
404 + susceptible_population_region_1_NEXT
405 + immune_population_region_1_NEXT
406 )
407
408 recovered_population_region_1 = recovered_population_region_1_NEXT
409 recovered_population_region_2 = recovered_population_region_2_NEXT
410
411 infected_population_region_1 = infected_population_region_1_NEXT
412 infected_population_region_2 = infected_population_region_2_NEXT
413
414 susceptible_population_region_1 = susceptible_population_region_1_NEXT
415 susceptible_population_region_2 = susceptible_population_region_2_NEXT
416
417 immune_population_region_1 = immune_population_region_1_NEXT
418 immune_population_region_2 = immune_population_region_2_NEXT
419
420 deceased_population_region_1.append(deceased_population_region_1_NEXT)
421 deceased_population_region_2.append(deceased_population_region_2_NEXT)
422
423 # End of main code
424
425 return {"TIME": runTime, "deceased_population_region_1": deceased_population_region_1}
426
427
428if __name__ == "__main__":
429 ema_logging.log_to_stderr(ema_logging.INFO)
430
431 model = Model("mexicanFlu", function=flu_model)
432 model.uncertainties = [
433 RealParameter("x11", 0, 0.5),
434 RealParameter("x12", 0, 0.5),
435 RealParameter("x21", 0.0001, 0.1),
436 RealParameter("x22", 0.0001, 0.1),
437 RealParameter("x31", 0, 0.5),
438 RealParameter("x32", 0, 0.5),
439 RealParameter("x41", 0, 0.9),
440 RealParameter("x51", 0, 0.5),
441 RealParameter("x52", 0, 0.5),
442 RealParameter("x61", 0, 0.8),
443 RealParameter("x62", 0, 0.8),
444 RealParameter("x81", 1, 10),
445 RealParameter("x82", 1, 10),
446 RealParameter("x91", 0, 0.1),
447 RealParameter("x92", 0, 0.1),
448 RealParameter("x101", 0, 200),
449 RealParameter("x102", 0, 200),
450 ]
451
452 model.outcomes = [TimeSeriesOutcome("TIME"), TimeSeriesOutcome("deceased_population_region_1")]
453
454 nr_experiments = 500
455
456 with SequentialEvaluator(model) as evaluator:
457 results = perform_experiments(model, nr_experiments, evaluator=evaluator)
458
459 print("laat")
460
461 lines(
462 results,
463 outcomes_to_show="deceased_population_region_1",
464 show_envelope=True,
465 density=Density.KDE,
466 titles=None,
467 experiments_to_show=np.arange(0, nr_experiments, 10),
468 )
469 plt.show()