.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_zombies.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_zombies.py: Zombie Apocalypse ================= author: Barron H. Henderson date: 2020-02-11 Inspired by numerous mathematics examples using Zombies (e.g., [1,2,3]), this is a test of the pykpp system to solve for a zombie outbreak. The basic model is that zombies (Z) both kill and infect humans (H), and humans kill zombies. I have also added a natural death and birth rate for humans, though it makes little difference on the time scales here. I have made additional assumptions as follows: 1. uniform distribution of humans and zombies, 2. interactions are proportional to density, 3. like infectous diseases, the zombie infection rates increase in winter 4. human militarization against zombies increases with the total zombies. To run this exmaple, you must have installed pykpp python -m pip install https://github.com/barronh/pykpp/archive/main.zip [1] https://people.maths.ox.ac.uk/maini/PKM%20publications/384.pdf [2] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4798798/ [3] https://www.amazon.com/dp/0691173206/ref=cm_sw_em_r_mt_dp_688KKQ5HMHJW9G6Z6654 .. GENERATED FROM PYTHON SOURCE LINES 31-39 .. code-block:: default import io import pykpp import pandas as pd import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 40-56 Define Reactions ----------------- * R1 : Zombies infect humans with a seasonal variation * k_inf is set later as a constant * A multiplier ranging from 0.25 in the summer to 1.75 in the winter * R2 : zombies kill humans using a constant coefficient * R3 : humans kill zombies using a coefficient made of two parts: * k_hkillz : humans naturally defend themselves killing zombies * Z * AREA : the military or police gets involved proportional to the total zombie population. * Right now, human militarization increase with zombie population, but it decreases immediately in response. This is unrealistic. It would be good to add a "memory" so that the zombies are removed more effectively. * R4, R5 : humans naturally procreate and die .. GENERATED FROM PYTHON SOURCE LINES 56-67 .. code-block:: default mechanism = """ #EQUATIONS {R1} Zombie + Human = 2 Zombie + Infected : k_inf * (np.cos(np.radians(t / 365 * 360)) * 0.5 + 1); {infections increase in winter} {R2} Zombie + Human = Zombie + KilledHuman : k_zkillh; {R3} Zombie + Human = Human + KilledZombie : k_hkillz * (Zombie * AREA); {humans militarize against zombies proportional to total zombie population} {R4} Human = 2 Human : k_hbirth; {R5} Human = DeadHuman : k_hdeath; """ .. GENERATED FROM PYTHON SOURCE LINES 68-84 Initialize the model ----------------- * time (t) is initialized at the start * the end is set to 10 years hences * Latitude/Longitude/Temperature/M are unused, but added to satisfy exectations of pykpp * AREA is set to USA area in km2 * constant rate coefficients are set * monitoring the Human and Zombie population * using the scipy.ode integrator * Finally setting the initial values * CFACTOR is used to convert populations to population/km2 * Humans (H) is set to an approximate USA population * Zombies (Z) is set to a very small initial population .. GENERATED FROM PYTHON SOURCE LINES 84-118 .. code-block:: default setup = """ #INLINE PY_INIT t=TSTART=0 # days TEND=365.25*10+TSTART TEMP = 288.; # unused M = 0.; # unused DT = 3/24. StartDate = 'datetime(2010, 7, 14)' Latitude_Degrees = 40. # Unused Longitude_Degrees = 0.00E+00 # Unused AREA = 9.834e6 {USA area in km2} k_inf = 0.5/365 # 0.5 km2/person/year k_zkillh = 0.2/365 # 0.2 km2/person/year k_hkillz = 0.1/365 # 0.1 km2/person/year k_hbirth = 0.01/365 # 1% per year k_hdeath = 0.01/365 # 1% per year #ENDINLINE """ config = """ #MONITOR Human; Zombie; #INTEGRATOR odeint; """ init = """ #INITVALUES CFACTOR = 1./AREA {people-to-people/km2 based on USA area} ALL_SPEC=1e-32*CFACTOR; Human = 300e6 {approximate US population} Zombie = 10 """ mechdef = setup + '\n' + config + '\n' + init + '\n' + mechanism .. GENERATED FROM PYTHON SOURCE LINES 119-127 Load the model -------------- * Using IO to create an in-memory file * Creating a mechanism: * called 'zombies', so the output will be 'zombies.pykpp.tsv' * solving at a 1 day increment * monitoring as it is solved per year .. GENERATED FROM PYTHON SOURCE LINES 127-134 .. code-block:: default mech = pykpp.mech.Mech( io.StringIO(mechdef), mechname='zombies', incr=1, monitor_incr=365, verbose=0 ) .. rst-class:: sphx-glr-script-out .. code-block:: none Species: 6 Reactions: 5 .. GENERATED FROM PYTHON SOURCE LINES 135-139 Reset and Run ------------- * Run the model and print the total time to solve .. GENERATED FROM PYTHON SOURCE LINES 139-143 .. code-block:: default runtime = mech.run(verbose=1, debug=False, solver='odeint'); # vode or lsoda print('Runtime (s)', runtime) .. rst-class:: sphx-glr-script-out .. code-block:: none tstart: 0 tend: 3652.5 dt: 0.125 solver_keywords: {} {t:0,Human:3E+08,Zombie:10} odeint {'atol': 0.001, 'rtol': 0.0001, 'mxstep': 1000, 'hmax': 0.125, 'mxords': 2, 'mxordn': 2} {t:365,Human:3E+08,Zombie:7.3} {t:730,Human:3E+08,Zombie:7.3} {t:1095,Human:3E+08,Zombie:7.3} {t:1460,Human:3E+08,Zombie:7.3} {t:1826,Human:3E+08,Zombie:7.3} {t:2191,Human:3E+08,Zombie:7.3} {t:2556,Human:3E+08,Zombie:7.3} {t:2921,Human:3E+08,Zombie:7.3} {t:3286,Human:3E+08,Zombie:7.3} {t:3651,Human:3E+08,Zombie:7.3} Runtime (s) 8.203604698181152 .. GENERATED FROM PYTHON SOURCE LINES 144-150 Plot the Zombie Statistics -------------------------- * Using pandas to read the file * Using matplotlib for control over plotting. * note that the seasonality of infection rate causes the wobble. .. GENERATED FROM PYTHON SOURCE LINES 150-184 .. code-block:: default def zombieplot(tsvpath, timelabel='Years'): data = pd.read_csv(tsvpath, delimiter='\t', index_col='t') if timelabel.lower() == 'years': tfac = 1 / 365.25 elif timelabel.lower() == 'days': tfac = 1 else: raise KeyError(f'timelabel must be Years or Days; got {timelabel}') t = data.index.values * tfac kH = data.Human / 1e6 Z = data.Zombie dZ = data.KilledZombie kI = data.Infected gskw = dict(wspace=0.2, left=.2) fig, axx = plt.subplots(1, 3, figsize=(14, 4), gridspec_kw=gskw) plt.setp(axx, xlabel=timelabel) plt.rc('lines', marker='o', linestyle='none', markersize=1) axx[0].plot(t, kH, label='Human Population (1e6)') mfmt = plt.matplotlib.ticker.FuncFormatter('{0:.5f}'.format) axx[0].yaxis.set_major_formatter(mfmt) axx[1].plot(t, Z, label='Zombies Population') axx[2].plot(t[1:], np.diff(dZ), label='Zombies Killed Rate (per day)') axx[2].plot(t[1:], np.diff(kI), label='Infection Rate (per day)') for ai, ax in enumerate(axx): ax.legend() axx[1].set_ylim(0, None); axx[2].set_ylim(0, None); return fig fig = zombieplot('zombies.pykpp.tsv') fig.savefig('zombies.pykpp.png') .. image-sg:: /auto_examples/images/sphx_glr_plot_zombies_001.png :alt: plot zombies :srcset: /auto_examples/images/sphx_glr_plot_zombies_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 185-188 Re run with adjustments parameters ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 188-203 .. code-block:: default mech = pykpp.mech.Mech( io.StringIO(mechdef), mechname='zombies2', incr=1, monitor_incr=365, verbose=0 ) # Increase initial Zombies by 2 and the infection rate by 20 mech.world['Zombie'] *= 2 mech.world['k_inf'] *= 20 runtime = mech.run(verbose=1, debug=False, solver='odeint'); # vode or lsoda fig = zombieplot('zombies2.pykpp.tsv'); fig.savefig('zombies2.pykpp.png') .. image-sg:: /auto_examples/images/sphx_glr_plot_zombies_002.png :alt: plot zombies :srcset: /auto_examples/images/sphx_glr_plot_zombies_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Species: 6 Reactions: 5 tstart: 0 tend: 3652.5 dt: 0.125 solver_keywords: {} {t:0,Human:3E+08,Zombie:20} odeint {'atol': 0.001, 'rtol': 0.0001, 'mxstep': 1000, 'hmax': 0.125, 'mxords': 2, 'mxordn': 2} {t:365,Human:3E+08,Zombie:1.5E+02} {t:730,Human:3E+08,Zombie:1.5E+02} {t:1095,Human:3E+08,Zombie:1.5E+02} {t:1460,Human:3E+08,Zombie:1.5E+02} {t:1826,Human:3E+08,Zombie:1.5E+02} {t:2191,Human:3E+08,Zombie:1.5E+02} {t:2556,Human:3E+08,Zombie:1.5E+02} {t:2921,Human:3E+08,Zombie:1.5E+02} {t:3286,Human:3E+08,Zombie:1.5E+02} {t:3651,Human:3E+08,Zombie:1.5E+02} .. GENERATED FROM PYTHON SOURCE LINES 204-218 Equations of the End: Teaching Mathematical Modeling Using the Zombie Apocalypse -------------------------------------------------------------------------------- Lofgren et al. 2016[1] published a teaching example and an on-line playground[2] with solutions. This offers a really cool opportunity to check our results. We get the exact same results -- hooray! Try changing the parameters. For example, adjust the parameters to represent COVID. Consider adding effects of temperature and humidity in transmission. [1] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4798798/#s1-jmbe-17-134 [2] http://cartwrig.ht/apps/whitezed/ .. GENERATED FROM PYTHON SOURCE LINES 218-271 .. code-block:: default import pykpp import pandas as pd zeddef = """ #EQUATIONS {R1} Zombie + Susceptible = Zombie + Infected : k_inf; {R2} Zombie + Susceptible = Zombie + KilledHuman : k_zkillh; {R3} Zombie + Susceptible = Susceptible + KilledZombie : k_hkillz; {R5} Zombie = DecayedZombie : k_zdeath; {R6} Susceptible = Vaccinated : k_vaccination; {R7} Infected = Zombie : k_incubation; #INLINE PY_INIT t=TSTART=0 # days TEND=TSTART + 60 TEMP = 288.; # unused, but represents average surface temperature P = 101325.; # unused, but represents average surface pressure M = P / 8.314 / TEMP * 6.022e23; # unused DT = 1/24. StartDate = 'datetime(2010, 7, 14)' Latitude_Degrees = 40. # Unused Longitude_Degrees = 0.00E+00 # Unused k_inf = 1. {infection rate} k_zkillh = 0.0 {zombie kill human rate} k_hkillz = 0.0 {human kill zombie rate} k_zdeath = 0.0 {zombie decay} k_vaccination = 0.0 {rate of susceptible vaccination} incubation_time_days = 1e-3 {must be greater than zero} k_incubation = 1 / (incubation_time_days) { tau = 1 / k; so tau = 1 seconds} #ENDINLINE #MONITOR Susceptible; Vaccinated; Zombie; #INTEGRATOR odeint; #INITVALUES CFACTOR = 1./10001 { denominator should match sum of people = Susceptible + Vaccinated + Zombie} ALL_SPEC=1e-32*CFACTOR; Susceptible = 10000 Vaccinated = 0 Zombie = 1 """ mech = pykpp.mech.Mech( io.StringIO(zeddef), mechname='zed', incr=1/24, monitor_incr=5, verbose=0 ) runtime = mech.run(verbose=1, debug=False, solver='odeint'); # vode or lsoda data = pd.read_csv('zed.pykpp.tsv', delimiter='\t', index_col='t') ax = data.drop(['TEMP', 'P', 'CFACTOR'], axis=1).plot() ax.figure.savefig('zed.pykpp.png') .. image-sg:: /auto_examples/images/sphx_glr_plot_zombies_003.png :alt: plot zombies :srcset: /auto_examples/images/sphx_glr_plot_zombies_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Species: 7 Reactions: 6 tstart: 0 tend: 60 dt: 0.041666666666666664 solver_keywords: {} {t:0,Susceptible:1E+04,Vaccinated:0,Zombie:1} odeint {'atol': 0.001, 'rtol': 0.0001, 'mxstep': 1000, 'hmax': 0.041666666666666664, 'mxords': 2, 'mxordn': 2} {t:5,Susceptible:9.9E+03,Vaccinated:0,Zombie:1.5E+02} {t:10,Susceptible:3.1E+03,Vaccinated:0,Zombie:6.9E+03} {t:15,Susceptible:30,Vaccinated:0,Zombie:1E+04} {t:20,Susceptible:0.22,Vaccinated:0,Zombie:1E+04} {t:25,Susceptible:0.0016,Vaccinated:0,Zombie:1E+04} {t:30,Susceptible:1.1E-05,Vaccinated:0,Zombie:1E+04} {t:35,Susceptible:8.3E-08,Vaccinated:0,Zombie:1E+04} {t:40,Susceptible:5.8E-10,Vaccinated:0,Zombie:1E+04} {t:45,Susceptible:4E-12,Vaccinated:0,Zombie:1E+04} {t:50,Susceptible:2.9E-14,Vaccinated:0,Zombie:1E+04} {t:55,Susceptible:2E-16,Vaccinated:0,Zombie:1E+04} .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 18.513 seconds) .. _sphx_glr_download_auto_examples_plot_zombies.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_zombies.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_zombies.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_