<jupyter_start><jupyter_text>Covid-19: From model prediction to model predictive control A demo of a stochastic difference equation (SDE) SIR model*Code by <NAME> in consultation with the BIOMATH research unit headed by prof. <NAME> and prof. <NAME> of the KERMIT research unit.*Copyright (c) 2021 by <NAME>, BIOMATH, Ghent University. All Rights Reserved.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from covid19model.models import models
from covid19model.data import mobility, sciensano, model_parameters

# OPTIONAL: Load the "autoreload" extension so that package code can change
%load_ext autoreload
# OPTIONAL: always reload modules so that as you change code in src, it gets loaded
%autoreload 2
###############
## Load data ##
###############

age_stratification_size=3 #3, 9 or 10
# Initial conditions, time-integrated prepandemic contact matrices, model parameters of COVID19 SEIQRD
initN, Nc_dict, parameters = model_parameters.get_COVID19_SEIQRD_parameters(age_stratification_size)
# Google Mobility data
df_google = mobility.get_google_mobility_data(update=False)
# Sciensano public data
df_hosp, df_mort, df_cases, df_vacc = sciensano.get_sciensano_COVID19_data(update=False)
# Aggregate the data on the national level
df_hosp = df_hosp.groupby(level='date').sum()['H_in']
############################
## Define a draw function ##
############################

def draw_fcn(param_dict,samples_dict):
    return param_dict
###################################
## Define social policy function ##
###################################

from covid19model.models.time_dependant_parameter_fncs import ramp_fun
def compliance_func(t, states, param, l, effectivity):
    
    # Convert tau and l to dates
    l_days = pd.Timedelta(l, unit='D')
    # Define start of measures
    start_measures = pd.to_datetime('2020-03-15')
    # List of transitioning dates
    change_days = [start_measures,
                   start_measures + l_days]
    # List of transitioning values
    change_values = [param,
                     ramp_fun(param, effectivity*param, t, start_measures, l),
                     param * effectivity]
    # Return appropriate policy
    try:
        return change_values[[index for index,value in enumerate(change_days) if value > t][0]]
    except:
        return change_values[-1]
######################
## Initialize model ##
######################

# Define parameters
l = 11.5
effectivity = 0.13
params = {'beta': 0.020, 'gamma': 1/7, 'Nc': Nc_all['total'], 'l': l, 'effectivity': effectivity}
# Define initial states
initial_states = {"S": initN.values, "I": np.ones(initN.size)}
# Initialize model (SDEs are discrete time models and require the discrete=True flag upon initialization)
model = models.simple_stochastic_SIR(initial_states, params, time_dependent_parameters={'Nc': compliance_func}, discrete=True)
########################
## Perform simulation ##
########################

n_draws = 50
start_sim = '2020-03-15'
end_sim = '2020-07-15'
out = model.sim(end_sim,start_date=start_sim,warmup=4,N=n_draws,draw_fcn=draw_fcn,samples={})
#######################
## Visualize results ##
#######################

fig, ax = plt.subplots(figsize=(12,4))
for i in range(n_draws):
    ax.plot(out['time'],out['I'].sum(dim='Nc').isel(draws=i), color='red', alpha = 0.05, linewidth = 3)
ax.scatter(df_sciensano[start_sim:end_sim].index,df_sciensano['H_in'][start_sim:end_sim],color='black',alpha=0.8,linestyle='None',facecolors='none')
ax.xaxis.set_major_locator(plt.MaxNLocator(5))
ax.yaxis.set_major_locator(plt.MaxNLocator(4))
plt.show()
plt.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18127053669222343, 'fraction_numerical': 0.0208105147864184, 'mean_word_length': 5.666058394160584, 'pattern_counts': {'":': 2, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '161585', 'license': 'MIT', 'n_tokens_mistral': 1206, 'n_tokens_neox': 1132, 'n_words': 295, 'path': 'notebooks/templates/simple_stochastic_SIR_demo.ipynb', 'repo_name': 'UGentBiomath/COVID19-Model', 'size': '3676'}