covid-19 plot

main.py

#! /usr/bin/env python3
import math
import datetime
import csv
import matplotlib.pyplot as plt
import numpy as np

t_italy = None
italy = None

t_spain = None
spain = None

t_germany = None
germany = None

# https://systems.jhu.edu/research/public-health/ncov/
# https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv
with open('time_series_19-covid-Confirmed.csv') as csv_file:
    r = csv.reader(csv_file, delimiter=',')
    
    def parse_date(ts):
        return [datetime.datetime.strptime(t, '%m/%d/%y') for t in ts]
    def parse_value(vs):
        return [int(v) for v in vs]

    data_start = 4 # first entries are geographic location
    
    corona_start_in_italy = 33 # 2/20/20
    corona_start_in_spain = 33 # 2/20/20
    corona_start_in_germany = 39 # 2/26/20
    for row in r:
        # row[1] is the country

        # the header has the date
        if row[1] == 'Country/Region':
            t_italy = parse_date(row[corona_start_in_italy:])
            t_spain = parse_date(row[corona_start_in_spain:])
            t_germany = parse_date(row[corona_start_in_germany:])

        if row[1] == 'Italy':
            italy = parse_value(row[corona_start_in_italy:])
        if row[1] == 'Spain':
            spain = parse_value(row[corona_start_in_spain:])
        if row[1] == 'Germany':
            germany = parse_value(row[corona_start_in_germany:])

def date_relative(ts, offset=3):
    # Offset nimmt an, dass die infektionen schon seit 3 Tagen laufen.
    # Chosen quite arbitrarily to better fit the curve.
    start = ts[0]
    return [(d-start).days + offset for d in ts]

# Wenn ich das richtig gemacht habe, ist das eine Funktion die sich alle 2.8 Tage verdoppelt.
# Die Wachstumstrate im unteren graphen wird diese Annahme sanity-checken.
doubles_every = 2.8 # days
growth = math.e ** (math.log(2) / doubles_every)


initial_cases = 20 # chosen arbitrarily to fit the curve

def exp_model(ts):
    """Modell fuer expontntielles Wachstum."""

    # Formel von https://www.youtube.com/watch?v=Kas0tIxDvrg
    return [pow(growth, d)*initial_cases for d in ts]

print(f"exp_model mit {growth} growth rate")

extrapolation = 10 # days

def add_days(ts, days=extrapolation):
    # days=14 assumes we want to extrapolate 2 weeks.
    t2 = list(ts)
    for d in range(days):
        t2.append(t2[-1] + datetime.timedelta(days=1))
    return t2 

fig, (top, bot) = plt.subplots(nrows=2, sharex=True, gridspec_kw={'height_ratios': [3, 1]}, figsize=(18, 12))

fig.suptitle(f"""COVID-19 Modell für die nächsten {extrapolation} Tage (vom {datetime.datetime.now().date()}).
Keine Panik! Keine Hamsterkäufe (wir haben genug Nudeln mit Toilettenpapiersauce)! Hände Waschen!
Verbreitung durch Selbst-Quarantäne verlangsamen (auch wenn du keine Symptome hast)!
Keine Fakes verbreiten! Quellen prüfen (diese Grafik ist nur echt mit Quellcode, der Autor ist weder Virologe noch Mediziner; das Modell ist stark vereinfacht, Daten können falsch sein)""")


color_italy = 'lime'
color_spain = 'seagreen'
color_germany = 'blue'
color_model = 'orange'


top.set(xlabel='Zeit (tage)', ylabel='Infizierte')

top.plot(t_italy, italy, label='Italien', color=color_italy)
top.plot(t_italy, exp_model(date_relative(t_italy)), '--', label=f'Modell exponentielles Wachstum seit {t_italy[0].date()}', color=color_model)

top.plot(t_spain, spain, label='Spanien', color=color_spain)
#top.plot(t_spain, exp_model(date_relative(t_spain)), '--', label=f'Modell exponentielles Wachstum seit {t_spain[0].date()}', color=color_model)

top.plot(t_germany, germany, label='Deutschland', color=color_germany)
top.plot(add_days(t_germany), exp_model(date_relative(add_days(t_germany))), '--', label=f'Modell exponentielles Wachstum seit {t_germany[0].date()} (gleiche Funktion, anderes Startdatum)', color=color_model)



full_time = add_days(t_italy)


#https://www.tagesschau.de/investigativ/ndr-wdr/coronavirus-unikliniken-101.html
"""Bundesgesundheitsminister Jens Spahn wies am Mittwoch in der Bundespressekonferenz darauf hin, dass hierzulande rund 28.000 Intensivbetten und 25.000 Beatmungsgeräte zur Verfügung stehen."""
intensivbetten = 28000
# beatmungsgeraete = 25000

"""Das Problem ist nur: Sie sind größtenteils belegt mit anderen schwerkranken Patienten. Im Durchschnitt sind sie zu etwa 80 Prozent ausgelastet."""
aktuelle_krankenhausauslastung = 0.8

# https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Steckbrief.html
# Annahme: Da 80% leicht bis milde verlaufen muessen maximal 20% ins Krankenhaus. Oder eventuell sogar weniger.
infizierte_die_ins_krankenhaus_muessen = 0.2


kapazitaetsgrenze_deutschland = intensivbetten / infizierte_die_ins_krankenhaus_muessen
#top.plot(full_time[-5:], [kapazitaetsgrenze_deutschland]*5, '-x', color='red', label='Kapazitätsgrenze deutsche Krankenhäuser')


kapazitaetsgrenze2_deutschland = (intensivbetten * (1-aktuelle_krankenhausauslastung)) / infizierte_die_ins_krankenhaus_muessen
top.plot(full_time[-5:], [kapazitaetsgrenze2_deutschland]*5, '-^', color='darkred', label=f'Kapazitätsgrenze deutsche Krankenhäuser wenn bereits {aktuelle_krankenhausauslastung*100}% der Betten belegt sind')

top.annotate(f"""Annahmen:
Alle {intensivbetten} Intensivbetten haben Beatmungsgeräte,
Es sind bereits {aktuelle_krankenhausauslastung*100}% der Betten belegt,
Nur {infizierte_die_ins_krankenhaus_muessen*100}% der Infizierten müssen tatsächlich behandelt werden""", xy=(full_time[16],25000))


top.annotate(f"""Exponentielles Modell: $f(d) = wachstumsrate^d \\cdot anfangsfälle$
wobei ich annehme
  $wachstumsrate={growth}$
  $anfangsfälle={initial_cases}$
Bedeutet, wir fangen mit {initial_cases} Infizierten an und
  jeder Infizierte steckt durchschnittlich pro Tag
  {growth:.2f} Leute an (einschl. sich selbst).
Dies entspricht einer Verdopplung alle {doubles_every} Tage.
Die angenommene Wachstumrate ist im unteren
   Graph mit den realen Zahlen verglichen.""", xy=(full_time[0],5000))


top.legend()


def growth_factor(vs):
    res = []
    for i, v in enumerate(vs):
        if i == 0:
            v_prev = vs[0]
        else:
            v_prev = vs[i-1]
        res.append(v/v_prev)
    return res


print("Median Wachstumsrate Italien:", np.median(growth_factor(italy)))
print("Median Wachstumsrate Deutschland:", np.median(growth_factor(germany)))
print("Median Wachstumsrate Modell:",np.median(growth_factor(exp_model(date_relative(full_time)))))

def annotated_median_growth_factor(vs):
    median = np.median(growth_factor(vs))
    langsamer_schneller = 'langsamer' if median < np.median(growth_factor(exp_model(date_relative(full_time)))) else 'schneller'
    return f"""{median:.4f} (wächst {langsamer_schneller} als im Modell)"""

bot.annotate(f"""Median Wachstumsrate Italien: {annotated_median_growth_factor(italy)}
Median Wachstumsrate Spanien: {annotated_median_growth_factor(spain)}
Median Wachstumsrate Deutschland: {annotated_median_growth_factor(germany)}
Median Wachstumsrate Modell: {np.median(growth_factor(exp_model(date_relative(full_time)))):.4f}""", xy=(full_time[-15],1.5))

bot.set(ylabel='Wachstumsrate\n(fälle_heute / fälle_gestern)')
bot.set_ylim(0.8, 2.3)
bot.plot(t_italy, growth_factor(italy), label='Italien', color=color_italy)
bot.plot(t_spain, growth_factor(spain), label='Spanien', color=color_spain)
bot.plot(t_germany, growth_factor(germany), label='Deutschland', color=color_germany)
bot.plot(full_time, growth_factor(exp_model(date_relative(full_time))), '--', label='Exponentielles Modell', color=color_model)
#bot.legend()

fig.autofmt_xdate()

plt.show()