<jupyter_start><jupyter_text>**Calibration of the age-stratified deterministic model***Original code by <NAME>. Modified by <NAME> in consultation with the BIOMATH research unit headed by prof. <NAME>.*Copyright (c) 2020 by <NAME>, BIOMATH, Ghent University. All Rights Reserved.This notebook accompanies our preprint: "*A deterministic, age-stratified, extended SEIRD model for assessing the effect of non-pharmaceutical interventions on SARS-CoV-2 spread in Belgium*"(https://doi.org/10.1101/2020.07.17.20156034)  Load required packages<jupyter_code>import random
import os
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from IPython.display import Image
from ipywidgets import interact,fixed,FloatSlider,IntSlider,ToggleButtons
import pandas as pd
import datetime
import scipy
from scipy.integrate import odeint
import matplotlib.dates as mdates
import matplotlib
import scipy.stats as st

import math
import xarray as xr
import emcee
import json
import corner

from covid19model.optimization import objective_fcns
from covid19model.optimization import MCMC
from covid19model.models import models
from covid19model.data import google
from covid19model.data import sciensano
from covid19model.data import polymod
from covid19model.data import model_parameters
from covid19model.visualization.output import population_status, infected
from covid19model.visualization.optimization import plot_fit, traceplot
from covid19model.optimization.run_optimization import full_calibration

# 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
import dask

from dask.distributed import Client, progress
client = Client(threads_per_worker=16, n_workers=1)
client<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>initN, Nc_home, Nc_work, Nc_schools, Nc_transport, Nc_leisure, Nc_others, Nc_total = polymod.get_interaction_matrices()
levels = initN.size
df_sciensano = sciensano.get_sciensano_COVID19_data(update=False)
df_sciensano.tail()
nonpublic_ts = pd.read_csv('../data/interim/non_public/all_nonpublic_timeseries.csv', parse_dates=['DATE'])
nonpublic_ts.NIS.unique()<jupyter_output><empty_output><jupyter_text>Initialize the model<jupyter_code># Load the parameters using `get_COVID19_SEIRD_parameters()`.
params = model_parameters.get_COVID19_SEIRD_parameters()
# Add the delayed ramp parameters to the parameter dictionary.
params.update({'l': 1,
              'tau': 5})
# Define the initial condition: one exposed inidividual in every age category
initial_states = {'S': initN, 'E': np.ones(levels)}
# Load the compliance model
from covid19model.models.compliance import ramp_2
# Initialize the model
model = models.COVID19_SEIRD(initial_states, params, compliance=ramp_2)
start_date = '2020-03-15'
end_beta = '2020-03-22'
end_ramp = '2020-04-20'
fig_path = '../../results/calibrations_arrondissements/'
samples_path = '../../data/interim/model_parameters/arrondissements/'<jupyter_output><empty_output><jupyter_text>Run in parallel<jupyter_code>arr_list = list(nonpublic_ts.NIS.unique())
def run_me_parallel(arr):
    arrond_ts = nonpublic_ts.pivot(index='DATE', columns='NIS', values='hospitalised_IN')[arr]
    samples_dict = full_calibration(model, arrond_ts, arr, start_date, end_beta, end_ramp, 
                                fig_path, samples_path)
                                #maxiter=10, popsize=10, steps_mcmc=250)
    
    states = [['H_in']]
    end_date = '2020-09-20'
    data=[arrond_ts[start_date:end_ramp].values]
    fig,ax=plt.subplots()
    for i in range(200):
        idx,model.parameters['beta'] = random.choice(list(enumerate(samples_dict['beta'])))
        idx,model.parameters['l'] = random.choice(list(enumerate(samples_dict['l'])))
        model.parameters['tau'] = samples_dict['tau'][idx]
        prevention = samples_dict['prevention'][idx]
        # Create a dictionary of past policies
        chk = {'time':   [start_date], 
              'Nc':      [prevention*(Nc_home + 0.4*Nc_work + 0.3*Nc_transport + 0.7*Nc_others + 0.2*Nc_leisure)]
              }
        y_model = model.sim(time=end_date, excess_time=samples_dict['lag_time'],checkpoints=chk)
        ax = plot_fit(y_model,data,start_date,samples_dict['lag_time'],states,end_date=end_date,with_ints=False,ax=ax,plt_kwargs={'color':'blue','linewidth': 2,'alpha': 0.05})

    data_after_calib = arrond_ts[pd.to_datetime(end_ramp)+pd.to_timedelta('1d'):end_date]
    plt.scatter(data_after_calib.index, data_after_calib.values, marker='o',color='red',linestyle='None',facecolors='none')
    legend_text=['daily \nhospitalizations']
    ax.set_xlim('2020-03-10', '2020-08-03')
    fig.savefig('../../results/calibrations_arrondissements/'+arr+'_'+str(datetime.date.today())+'.pdf',
                bbox_inches='tight', dpi=600)
    return
tasks = []
for arr in arr_list:
    task = dask.delayed(run_me_parallel)(arr)
    tasks.append(task)
    
dask.compute(*tasks, scheduler='processes')
run_date = '2020-08-22'
with open(samples_path+'44000'+'_'+run_date+'.json', 'r') as fp:
    samples_dict = json.load(fp)<jupyter_output><empty_output><jupyter_text>Run one arrondissement<jupyter_code>run_me_parallel(58000)<jupyter_output>No constraints given.
New best for swarm at iteration 1: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 1: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 2: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 3: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 4: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 5: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 6: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 7: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 8: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 9: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610752739398591
Best after iteration 10: [1.00000000e+00 5.59021608e+01 2.00000000e-02] 7.610[...]
{'fraction_non_alphanumeric': 0.11382371591629677, 'fraction_numerical': 0.11033608116677235, 'mean_word_length': 5.47741273100616, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 22, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '402700', 'license': 'MIT', 'n_tokens_mistral': 2658, 'n_tokens_neox': 2161, 'n_words': 465, 'path': 'notebooks/JV-calibration-deterministic-arrondissements-parallel.ipynb', 'repo_name': 'JennaVergeynst/COVID19-Model', 'size': '6334'}