<jupyter_start><jupyter_text>Analytic Modeling of the Recurrent Outbreak of the Novel Coronavirus Disease 2019 Epidemics Case of The Republic of Croatia    <NAME>     <EMAIL>     Faculty of electrical engineering, mechanical engineering and naval architecture     University of Split     Split, Croatia **Description** $-$ This report is an attempt at analytical modeling of repeated outbreaks of the Coronavirus Disease 2019 (COVID-19) epidemic for the use case of The Republic of Croatia. The idea of the mentioned approach lies in observing the epidemic and the sudden jumps and falls of the total number of actively infected individuals as a series of *pulses* which represent the current presence of the virus in the population, conditioned by numerous factors. Results of the report show the direct dependence of the number of actively infected individuals on the effective reproduction number immediately before. Also, an exceptional fit of the dynamics of the number of active cases with the series of Heidler functions is achived through CoroPy, Python package for COVID-19 epidemics modeling.**Keywords** $-$ Coronavirus disease 2019, recurrent epidemic outbreaks modeling, Heidler function, multiwave SEIRD model, least-squares optimization**Software** $-$ Ubuntu 18.04 LTS, Python 3.7.6, CoroPy package<jupyter_code># IMPORT LIBS
import os
import datetime as dt

from coropy.growth_models import GrowthCOVIDModel
from coropy.compartmental_models import SEIRDModel
import coropy.reproduction_rate as R
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.metrics import r2_score

import plotting
# configure plotting
plotting.configure(grid=False)
hex_blue = plotting.default_colors('blue')
hex_red = plotting.default_colors('red')
hex_purple = plotting.default_colors('purple')
# load data
DATA_DIR = 'data'
FILE_NAMES = ['1st_wave_data_cro.csv',
              '2nd_wave_data_cro.csv']
FILE_PATHS = [os.path.join(DATA_DIR, FILE_NAMES[0]),
              os.path.join(DATA_DIR, FILE_NAMES[1])]
df_1 = pd.read_csv(FILE_PATHS[0])
df_2 = pd.read_csv(FILE_PATHS[1])
df = pd.concat([df_1, df_2[1:-1]], sort=False)

cum_positives = np.cumsum(df.confirmed_positives.values)
cum_recovered = df.estimate_recovered.values
cum_deceased = np.cumsum(df.confirmed_deceased.values)
active = cum_positives - cum_recovered - cum_deceased

# start of different epi waves
eff_dates=[dt.datetime(2020, 2, 25), dt.datetime(2020, 6, 1), dt.datetime(2020, 8, 8)]<jupyter_output><empty_output><jupyter_text>1 Heidler function<jupyter_code>def heidler(t, t1, t2, n, a, b):
    exp_1 = 1 - np.exp(-t/t1)
    exp_2 = np.exp(-t/t2)
    # if np.sign is not used NumPy will generate RuntimeWarning: invalid value encountered in power
    # https://stackoverflow.com/questions/45384602/numpy-runtimewarning-invalid-value-encountered-in-power#45384691
    return a/b * np.sign(exp_1) * (np.abs(exp_1))**np.abs(n) * exp_2  
X = []
Y = []

# past wave(s)
start_idx = 0
for start_date, end_date in zip(eff_dates[:-1], eff_dates[1:]):
    end_idx = start_idx+abs((end_date - start_date).days)
    
    model = GrowthCOVIDModel(heidler, normalize=True, calc_ci=False)
    model.fit(active[start_idx:end_idx])
    x, fitted = model.get_fitted
    X.extend(x)
    Y.extend(fitted)
    
    start_idx = end_idx  # update indexing

# current wave   
model = GrowthCOVIDModel(heidler, normalize=True, calc_ci=False)
model.fit(active[start_idx:])
x, fitted = model.get_fitted
X.extend(x)
Y.extend(fitted)
# goodness-of-fit
sre_fit = np.sqrt(np.mean((active - Y)**2))
r2_fit = r2_score(active, Y)
print(
    f'S = {sre_fit}',
    f'\nr2 = {r2_fit}'
)
# dates
t = np.arange(cum_positives.size)
start_dt = dt.datetime(2020, 2, 25)
end_dt = start_dt + dt.timedelta(int(t.max())+1)
dates = matplotlib.dates.drange(start_dt, end_dt, dt.timedelta(days=1))
fig, ax1 = plt.subplots()
ax1.plot(dates, active, marker='o', markersize=4, linestyle='none', c=hex_purple, markevery=1, label='$I(t)$')
ax1.plot(dates, Y, c=hex_blue, label=f'$\\hat I(t)$\n$S = {round(sre_fit, 2)}$, $R^2 = {round(r2_fit, 4)}$')
ax1.tick_params(axis='y', labelcolor=hex_blue)
ax1.set_ylabel('$N$', color=hex_blue)
ax1.legend(frameon=True)

ax2 = ax1.twinx()
ax2.plot(dates[:-15-3+1], R._estimate(cum_positives, 15, 3), 'kx', markevery=3, markersize=4, label='$R_t$(t)')
ax2.axhline(y=1, linestyle='--', linewidth=1, c='k', label='$R_t^{critical}$')
ax2.tick_params(axis='y', labelcolor='k')
ax2.set_ylabel('$R_t$', color='k')
ax2.legend(loc='lower right', frameon=True)

fig.tight_layout()
_ = fig.gca().xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y-%m-%d'))
_ = fig.gca().xaxis.set_major_locator(matplotlib.dates.DayLocator(interval=19))
_ = plt.gcf().autofmt_xdate()
plt.show()
fname = f'multiwave_Heidler'
fig.savefig(fname=os.path.join('figs', fname + '.pdf'), format='pdf', bbox_inches='tight')
fig.savefig(fname=os.path.join('figs', fname + '.png'), format='png', bbox_inches='tight', dpi=300)<jupyter_output><empty_output><jupyter_text>2 SEIRD model<jupyter_code>eff_dates=[dt.datetime(2020, 2, 25), dt.datetime(2020, 6, 6), dt.datetime(2020, 8, 8)]
eff_population_scaler = 1
first_wave_eff_population = 2200
S0 = first_wave_eff_population * eff_population_scaler
E0 = 3 * cum_positives[0]
I0 = cum_positives[0]
R0 = cum_recovered[0]
D0 = cum_deceased[0]
IC = (S0, E0, I0, R0, D0)

S_tot, E_tot, I_tot, R_tot, D_tot = [], [], [], [], []

# past wave(s)
start_idx = 0
for start_date, end_date in zip(eff_dates[:-1], eff_dates[1:]):
    end_idx = start_idx+abs((end_date - start_date).days)
    model = SEIRDModel()
    _, _ = model.fit(cum_positives[start_idx:end_idx], cum_recovered[start_idx:end_idx], cum_deceased[start_idx:end_idx], IC)
    (S, E, I, R, D) = model.simulate()
    S_tot.extend(S.tolist())
    E_tot.extend(E.tolist())
    I_tot.extend(I.tolist())
    R_tot.extend(R.tolist())
    D_tot.extend(D.tolist())
    
    eff_population_scaler += 1
    S0 = S0 * eff_population_scaler
    IC = (S0, 3 * I[-1], I[-1], R[-1], D[-1])  # update initial conditions
    start_idx = end_idx  # update indexing

# current wave   
model = SEIRDModel()
_, _ = model.fit(cum_positives[start_idx:], cum_recovered[start_idx:], cum_deceased[start_idx:], IC)
(S, E, I, R, D) = model.simulate()
S_tot.extend(S.tolist())
E_tot.extend(E.tolist())
I_tot.extend(I.tolist())
R_tot.extend(R.tolist())
D_tot.extend(D.tolist())
# goodness-of-fit
sre_fit = np.sqrt(np.mean((active - I_tot)**2))
r2_fit = r2_score(active, I_tot)
print(
    f'S = {sre_fit}',
    f'\nr2 = {r2_fit}'
)
# visualize
fig, ax = plt.subplots()
ax.plot(dates, active, marker='o', markersize=4, linestyle='none', c=hex_purple, markevery=3, label='$I(t)$')
ax.plot(dates, I_tot, c=hex_blue, label=f'$\\hat I(t)$\n$S = {round(sre_fit, 2)}$, $R^2 = {round(r2_fit, 4)}$')
ax.plot(dates, cum_recovered, marker='^', markersize=4, linestyle='none', c=hex_purple, markevery=3, label='$R(t)$')
ax.plot(dates, R_tot, '--', c=hex_blue, label=f'$\\hat R(t)$')
ax.plot(dates, cum_deceased, marker='x', markersize=5, linestyle='none', c=hex_purple, markevery=3, label='$D(t)$')
ax.plot(dates, D_tot, '-.', c=hex_blue, label=f'$\\hat D(t)$')
ax.tick_params(axis='y', labelcolor=hex_blue)
ax.set_ylabel('$N$', color=hex_blue)
ax.legend(frameon=True)

fig.tight_layout()
_ = fig.gca().xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y-%m-%d'))
_ = fig.gca().xaxis.set_major_locator(matplotlib.dates.DayLocator(interval=19))
_ = plt.gcf().autofmt_xdate()
plt.show()
fname = f'multiwave_SEIRD'
fig.savefig(fname=os.path.join('figs', fname + '.pdf'), format='pdf', bbox_inches='tight')
fig.savefig(fname=os.path.join('figs', fname + '.png'), format='png', bbox_inches='tight', dpi=300)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15690565066289097, 'fraction_numerical': 0.026515639078388466, 'mean_word_length': 5.573604060913706, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '274711', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3146, 'n_tokens_neox': 2949, 'n_words': 627, 'path': '02_multiwave_modeling.ipynb', 'repo_name': 'antelk/covid-19-modelling-croatia', 'size': '7794'}