<jupyter_start><jupyter_text>Load Denmark age structure and contact matrix<jupyter_code>my_data = np.genfromtxt('../data/age_structures/Denmark-2019.csv', delimiter=',', skip_header=1)
aM, aF = my_data[:, 1], my_data[:, 2]

Ni0=aM+aF;

M=16  ## number of age classes

Ni = Ni0[:M]
N=np.sum(Ni)

print("Age groups are in brackets of 5 (i.e. 0-4, 5-9, 10-14, .. , 75-79).")
print("Number of individuals in each bracket:")
print(Ni.astype('int'))
print("Total number of individuals: {0}".format(np.sum(Ni.astype('int'))))
# Get individual contact matrices
CH, CW, CS, CO = pyross.contactMatrix.Denmark()

# By default, home, work, school, and others contribute to the contact matrix
C = CH + CW + CS + CO

# Illustrate the individual contact matrices:
fig,aCF =  plt.subplots(2,2);
aCF[0][0].pcolor(CH, cmap=plt.cm.get_cmap('GnBu', 10));
aCF[0][1].pcolor(CW, cmap=plt.cm.get_cmap('GnBu', 10));
aCF[1][0].pcolor(CS, cmap=plt.cm.get_cmap('GnBu', 10));
aCF[1][1].pcolor(CO, cmap=plt.cm.get_cmap('GnBu', 10));<jupyter_output><empty_output><jupyter_text>C is the sum of contributions from contacts at home, workplace, schools and all other public spheres. Using superscripts $H$, $W$, $S$ and $O$ for each of these, we write the contact matrix as$$C_{ij} = C^H_{ij} + C^W_{ij} + C^S_{ij} + C^O_{ij}$$We read in these contact matrices from the data sets provided in the paper *Projecting social contact matrices in 152 countries using contact surveys and demographic data* by Prem et al, sum them to obtain the total contact matrix. We also read in the age distribution of UK obtained from the *Population pyramid* website.  Fit with real data<jupyter_code># Get the latest data from Johns Hopkins University
!git clone https://github.com/CSSEGISandData/COVID-19
cases = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
cases.head()
deaths = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv')
deaths.shape
cols = cases.columns.tolist() 
case = cases.loc[94,][4:]
death = deaths.loc[94,][4:]
plt.scatter(np.arange(len(case)),case, label= 'case',lw=0.01);
plt.scatter(np.arange(len(death)),death, label='death', lw=0.1);
plt.title('Denmark case/death since January 22')
plt.legend();<jupyter_output><empty_output><jupyter_text>Using data from the repo (not exactly the same)<jupyter_code>denmark_age_structured_case = pd.read_csv('../data/covid-cases/denmark_cases_by_age.csv')
denmark_age_structured_case.tail(10)
case = denmark_age_structured_case.sum(axis=0)[1:]<jupyter_output><empty_output><jupyter_text>We will use the **infective number** for (52 days) for the initial conditions Is and Ia.<jupyter_code>## infective people
nbday = denmark_age_structured_case.shape[1]-1
nbgroup = denmark_age_structured_case.shape[0]
I = np.array(np.zeros([nbday, nbgroup-2]))
for i in range(1,nbday+1):
    I[i-1] = np.array([denmark_age_structured_case.iloc[:,i][:nbgroup-2]])
## we don't want a cum sum 
I[1:] -= I[:-1].copy()
I[3,0] = 1
##  duplicate for each subgroup of age 
M = 16
Is = np.array(np.zeros([nbday, M]))
for i in range(I.shape[0]):
    Is[i] = np.array([val for val in I[i]/2 for _ in (0, 1)])

## take only integer values, loosing cases !
for j in range(Is.shape[0]):
    Is[j] = np.array([int(i) for i in Is[j]])<jupyter_output><empty_output><jupyter_text>Define model parameters<jupyter_code>beta  = 0.036692            # infection rate 

gE    = 1/5.
gA    = 1/4.
gIa   = 1./3.                # recovery rate of asymptomatic infectives 
gIs   = 1./7                # recovery rate of symptomatic infectives 
alpha = 0.4                 # fraction of asymptomatic infectives 
fsa   = 0.2                 # the self-isolation parameter   
fh    = 0
gIh   = 1/7
gIc   = 1/14
 
sa  = 0*np.ones(M)      #  rate of additional/removal of population by birth etc
sa[0]  = 0            # birth
sa[3] = 0          # mortality


hh   = 0.1*np.ones(M)       # fraction which goes from Is to hospital
hh[:M] = 0.01
hh[M:2*M] = 0.025
cc   = 0.05*np.ones(M)      # fraction which goes from  hospital to ICU 
mm   = 0.4*np.ones(M)       # mortality from IC


# initial conditions    
Is_0 = Is[0,]
# np.zeros((M)); 
Ia_0 = Is[0,]
# 10*np.ones((M)) # start with 10 asymptomatic infectives


R_0  = np.zeros((M))
E_0  = np.zeros((M))
A_0  = np.zeros((M))
Ih_0  = np.zeros((M))
Ic_0  = np.zeros((M))
Im_0  = np.zeros((M))

S_0  = Ni - (E_0 + A_0 + Ia_0 + Is_0 + Ih_0 + Ic_0 +Im_0 + R_0)<jupyter_output><empty_output><jupyter_text>Define events for protocol<jupyter_code># Dummy event for initial (standard) contact matrix
events = [lambda t: 1]
contactMatrices = [C]

# After 20 days, start lockdown
lockdown_threshold_0 = 20
def event0(t,rp):
    return t - lockdown_threshold_0
events.append(event0)
contactMatrices.append( CH + 0.1*CW + 0.4*CS )

# After 25 days, decrease contacts even further
lockdown_threshold_1 = 25
def event1(t,rp):
    return t- lockdown_threshold_1
events.append(event1)
contactMatrices.append( CH + 0.1*CW)

# After 70 days, add 50% of school contacts to contact matrix
lockdown_threshold_2 = 55
def event2(t,rp):
    return t - lockdown_threshold_2
events.append(event2)
contactMatrices.append( CH + 0.1*CW + 0.5*CS ) # everybody in lockdown<jupyter_output><empty_output><jupyter_text>Simulate protocol and analyse results Initialise pyross.control, run and plot a single test simulation<jupyter_code># duration of simulation
Tf=100;  Nf=Tf+1;  

# intantiate model
parameters = {'alpha':alpha,'beta':beta, 'gIa':gIa,'gIs':gIs,
              'gIh':gIh,'gIc':gIc,  'gE':gE, 'gA':gA,
              'fsa':fsa, 'fh':fh, 
              'sa':sa, 
              'hh':hh, 'cc':cc, 'mm':mm}
model = pyross.control.SEAI5R(parameters, M, Ni.copy())

# run model once
method='gillespie'
data=model.simulate(S_0, E_0, A_0, Ia_0, Is_0, Ih_0, Ic_0, Im_0,
                    events,contactMatrices, Tf, Nf,
                                method=method)
t = data['t']; 

y_plot = np.sum ( data['X'][:,4*M:7*M], axis = -1)


lw=2
fig,ax = plt.subplots(1,1,figsize=(8,5))
ax.plot(t,y_plot,lw=lw,)
ax.axvline(data['events_occured'][0][0],
           color='crimson',lw=lw,
          label='Beginning of lockdown',ls='--')
ax.axvline(data['events_occured'][-1][0],
 color='limegreen',lw=lw,
          label='Schools re-opened',ls='--')

ax.plot(case, 'o-', lw=4, color='#348ABD', ms=5, label='data', alpha=0.5)


ax.set_xlim(0,Tf)
fs=20
ax.legend(loc='best',fontsize=8,framealpha=1)
ax.set_xlabel('time since first infection [days]',fontsize=fs)
ax.set_ylabel('Number of known active cases',fontsize=fs)
plt.show(fig)


data_array = data['X']
t = data['t']
np.save('SEAI5R_latent_traj.npy', data_array)<jupyter_output><empty_output><jupyter_text>Infer latent data using first observations for Infected Ind<jupyter_code># load the data 
x = np.load('SEAI5R_latent_traj.npy').astype('float') # 9 equations * 16
x.shape<jupyter_output><empty_output><jupyter_text>We force the observed data into x.<jupyter_code>## Is non cumulative - 52 days from the data
dayobs = 52
# x[:dayobs,3*M:4*M].shape
x[:dayobs,3*M:4*M] = Is
x.shape
steps = 101 # number internal integration steps taken
# obs used for inference 
Tf_inference = 52
Nf_inference = Tf_inference + 1
x = x[:Nf_inference,:8*M]/N # remove the last equation only
x.shape 
# fraction of population in Age group
fi = Ni/sum(Ni)
Ni = N *fi
# initialise the estimator 
estimator = pyross.inference.SEAI5R(parameters, M, fi, int(N), steps)
## create filter to distinguish between observable classes
fltr=np.repeat(np.array([False, False, False, False, True, True, True, True]), M) # take only 4 equations
obs=x[:,fltr]
x0 = x[0]
obs.shape
## matrix with all events
def contactMatrix(t):
    return contactMatrices[3]<jupyter_output><empty_output><jupyter_text>import time more than 30 min ?start_time = time.time() logp = estimator.minus_logp_red(parameters, x0, obs[1:], fltr, Tf_inference, Nf_inference, contactMatrix)end_time = time.time()print(logp) print(end_time - start_time)<jupyter_code>Ia0_g = (Ia_0+0.2)/N 
Ia0_std = [5/N]*M
Ia0_b = [(0.1/N, 100/N)]*M
Is0_g = (Is_0+0.1)/N 
Is0_std = [5/N]*M  
Is0_b = [(0.1/N, 100/N)]*M
np.max(Is0_g)
max(S0_g)

[(1/N, f+0.8) for f in fi]
eps = 1e-4 

alpha_g = 0.25          
alpha_std = 0.1
alpha_b = (0.05, 0.5)

beta_g  = 0.06         
beta_std = 0.1
beta_b = (eps, 1)

gIa_g   = 0.13
gIa_std = 0.05 
gIa_b = (eps, 0.5)

gIs_g   = 0.13
gIs_std = 0.05 
gIs_b = (eps, 0.5) 

gE_g    = 0.21 
gE_std = 0.05
gE_b = (0.1, 0.3)

gA_g    = 0.35 
gA_std = 0.05 
gA_b = (0.2, 0.5) 
   

E0_g = (E_0+3)/N
E0_std = E0_g
E0_b = [(0.1/N, 100/N)]*M 

A0_g  = (A_0+3)/N
A0_std = A0_g
A0_b = [[0.1/N, 100/N]]*M

Ia0_g = (Ia_0+0.2)/N 
Ia0_std = [5/N]*M
Ia0_b = [(0.1/N, 200/N)]*M

Is0_g = (Is_0+0.1)/N 
Is0_std = [5/N]*M  
Is0_b = [(0.1/N, 200/N)]*M

Ih0_g = (Ih_0+0.1)/N 
Ih0_std = [5/N]*M 
Ih0_b = [(0.1/N, 100/N)]*M

Ic0_g = (Ic_0+0.1)/N 
Ic0_std = [5/N]*M 
Ic0_b = [(0.1/N, 100/N)]*M

Im0_g = (Im_0+0.1)/N 
Im0_std = [5/N]*M 
Im0_b = [(0.1/N, 100/N)]*M

S0_g = (S_0-10)/N 
S0_std = E0_std*4 
S0_b = [(1/N, f+0.8) for f in fi]

# optimisation parameters 
ftol = 1e-4 # the relative tol in (-logp)   
bounds = np.array([alpha_b, beta_b, gIa_b, gIs_b, gE_b, gA_b, 
                   *S0_b, *E0_b, *A0_b, *Ia0_b, *Is0_b, *Ih0_b, *Ic0_b, *Im0_b])
guess = np.array([alpha_g, beta_g, gIa_g, gIs_g, gE_g, gA_g, 
                   *S0_g, *E0_g, *A0_g, *Ia0_g, *Is0_g, *Ih0_g, *Ic0_g, *Im0_g])
stds = np.array([alpha_std, beta_std, gIa_std, gIs_std, gE_std, gA_std, 
                 *S0_std, *E0_std, *A0_std, *Ia0_std, *Is0_std, *Ih0_std, *Ic0_std, *Im0_std])


start_time = time.time()
params = estimator.latent_inference(guess, stds, obs, fltr, Tf_inference, Nf_inference, 
                                    contactMatrix, bounds, 
                                    global_max_iter=30, global_ftol_factor=1e3, verbose=True, ftol=ftol) 
end_time = time.time()
print(end_time - start_time)
# save the parameters 
np.save('SEAI5R_latent_param.npy', params)
print('True parameters')
print(parameters)

params = np.load('SEAI5R_latent_param.npy')
best_estimates = estimator.make_params_dict(params)
print('\n Best estimates')
print(best_estimates)

print('\n True initial conditions: ') 
print(S0, E0, A0, Ia0, Is0, Ih0, Ic0, Im0)
print('\n Inferred initial conditons: ') 
print((params[6:]*N).astype('int'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1573798791598734, 'fraction_numerical': 0.047089287426872545, 'mean_word_length': 4.015873015873016, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '484495', 'license': 'MIT', 'n_tokens_mistral': 4700, 'n_tokens_neox': 4392, 'n_words': 972, 'path': 'examples/inference/SEAI5R-Denmark.ipynb', 'repo_name': 'ineskris/pyross', 'size': '10427'}