# SCEIR model

Using PyGOM, we will set up my simple SCEIR model ODE system
PyGOM – A Python Package for Simplifying Modelling with Systems of Ordinary Differential Equations https://arxiv.org/pdf/1803.06934.pdf

## Equations

\begin{equation}
\begin{split}
\dot{S} &= -(\beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3) S + c_1 S_c - c_2*S*I_3\\
\dot{S_c} &= - c_0(\beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3) S_c - c_1 S_c + c_2*S*I_3\\
\dot{E} &=(\beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3 ) (S + c_0 S_c) - a E \\
\dot{I_1} &= a E - \gamma_1 I_1 - p_1 I_1 \\
\dot{I_2} &= p_1 I_1 -\gamma_2 I_2 - p_2 I_2 \\
\dot{I_3} & = p_2 I_2 -\gamma_3 I_3 - \mu I_3 \\
\dot{R} & = \gamma_1 I_1 + \gamma_2 I_2 + \gamma_3 I_3 \\
\dot{D} & = \mu I_3
\end{split}
\end{equation}

The use of I_3 as a state variable triggering susceptibles to execute caution is just one choice. In contrast with deaths, it does not accumulate over time and so retains the property of an active threat to society, rather than an historical one.

Actually, the effect of caution may be quadratic, since both the individual doing the infection and individual potentially being infected may be executing caution. The current version assumes that infectives do not change their precautionary measures in response to I_3. To implement this we would need to further distinguish careful and non careful infectives at least up to the I_1 level. This is done in the SCCEIR model.

## Variables
* $S$: Susceptible individuals living as normal
* $S_c$: Susceptible individuals exercising pandemy precautions
* $E$: Exposed individuals - infected but not yet infectious or symptomatic
* $I_i$: Infected individuals in severity class $i$. Severity increaes with $i$ and we assume individuals must pass through all previous classes
  * $I_1$: Mild infection (hospitalization not required)
  * $I_2$: Severe infection (hospitalization required)
  * $I_3$: Critical infection (ICU required)
* $R$: individuals who have recovered from disease and are now immune
* $D$: Dead individuals
* $N=S+S_c+E+I_1+I_2+I_3+R+D$ Total population size (constant)

## Parameters
* $\beta_i$ rate at which infected individuals in class $I_i$ contact susceptibles and infect them
* $a$ rate of progression from the exposed to infected class
* $\gamma_i$ rate at which infected individuals in class $I_i$ recover from disease and become immune
* $p_i$ rate at which infected individuals in class $I_i$ progress to class $I_{I+1}$
* $\mu$ death rate for individuals in the most severe stage of disease

* $c_i$ three parameters characterizing cautionary response of population via class $S_c$

    - c_0 : reduction factor for exposure for cautioned susceptibles

    - c_1 : duration of caution (exponential decay time constant in days)

    - c_2 : rate constant for transition from uncautioned to cautioned susceptible




## Implementation


```python
# import required packages
import os 
from sympy import symbols, init_printing
import numpy as np
import matplotlib.pyplot as plt
import sympy
import itertools
import scipy
from pygom import DeterministicOde, Transition, SimulateOde, TransitionType, SquareLoss
from scipy.optimize import minimize

```


```python
def  print_ode2(self):
        '''
        Prints the ode in symbolic form onto the screen/console in actual
        symbols rather than the word of the symbol.
        
        Based on the PyGOM built-in but adapted for Jupyter
        Corrected by John McCaskill to avoid subscript format error
        '''
        A = self.get_ode_eqn()
        B = sympy.zeros(A.rows,2)
        for i in range(A.shape[0]):
            B[i,0] = sympy.symbols('d' + '{' + str(self._stateList[i]) + '}'+ '/dt=')
            B[i,1] = A[i]

        return B
```


```python
# set up the symbolic SCEIR model

state = ['S', 'E', 'I_1', 'I_2','I_3','R','D','S_1']
param_list = ['beta_1', 'beta_2','beta_3','alpha', 'gamma_1', 'gamma_2', 'gamma_3',
              'p_1','p_2','mu','c_0','c_1','c_2','N']

transition = [
    Transition(origin='S', destination='E', equation='(beta_1*I_1+beta_2*I_2+beta_3*I_3)*S',
               transition_type=TransitionType.T),
    Transition(origin='S', destination='S_1', equation='c_2*I_3*S',
               transition_type=TransitionType.T),
    Transition(origin='S_1', destination='S', equation='c_1*S_1',
               transition_type=TransitionType.T),
    Transition(origin='S_1', destination='E', equation='c_0*(beta_1*I_1+beta_2*I_2+beta_3*I_3)*S_1',
               transition_type=TransitionType.T),
    Transition(origin='E', destination='I_1', equation='alpha*E',
               transition_type=TransitionType.T),
    Transition(origin='I_1', destination='R', equation='gamma_1*I_1',
               transition_type=TransitionType.T),
    Transition(origin='I_2', destination='R', equation='gamma_2*I_2',
               transition_type=TransitionType.T),
    Transition(origin='I_3', destination='R', equation='gamma_3*I_3',
               transition_type=TransitionType.T),
    Transition(origin='I_1', destination='I_2', equation='p_1*I_1',
               transition_type=TransitionType.T),
    Transition(origin='I_2', destination='I_3', equation='p_2*I_2',
               transition_type=TransitionType.T),
    Transition(origin='I_3', destination='D', equation='mu*I_3',
               transition_type=TransitionType.T)    
    ]


SCEIR_model = DeterministicOde(state, param_list, transition=transition)
```


```python
# display equations
print_ode2(SCEIR_model)
```




$\displaystyle \left[\begin{matrix}d{S}/dt= & - I_{3} S c_{2} - S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) + S_{1} c_{1}\\d{E}/dt= & - E \alpha + S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) + S_{1} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right)\\d{I_1}/dt= & E \alpha - I_{1} \gamma_{1} - I_{1} p_{1}\\d{I_2}/dt= & I_{1} p_{1} - I_{2} \gamma_{2} - I_{2} p_{2}\\d{I_3}/dt= & I_{2} p_{2} - I_{3} \gamma_{3} - I_{3} \mu\\d{R}/dt= & I_{1} \gamma_{1} + I_{2} \gamma_{2} + I_{3} \gamma_{3}\\d{D}/dt= & I_{3} \mu\\d{S_1}/dt= & I_{3} S c_{2} - S_{1} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) - S_{1} c_{1}\end{matrix}\right]$




```python
# display graphical representation of the model
SCEIR_model.get_transition_graph()
```


```python
SCEIR_model.get_ode_eqn()
```




$\displaystyle \left[\begin{matrix}- I_{3} S c_{2} - S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) + S_{1} c_{1}\\- E \alpha + S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) + S_{1} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right)\\E \alpha - I_{1} \gamma_{1} - I_{1} p_{1}\\I_{1} p_{1} - I_{2} \gamma_{2} - I_{2} p_{2}\\I_{2} p_{2} - I_{3} \gamma_{3} - I_{3} \mu\\I_{1} \gamma_{1} + I_{2} \gamma_{2} + I_{3} \gamma_{3}\\I_{3} \mu\\I_{3} S c_{2} - S_{1} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) - S_{1} c_{1}\end{matrix}\right]$




```python
SCEIR_model.get_jacobian_eqn()
```




$\displaystyle \left[\begin{matrix}- I_{1} \beta_{1} - I_{2} \beta_{2} - I_{3} \beta_{3} - I_{3} c_{2} & 0 & - S \beta_{1} & - S \beta_{2} & - S \beta_{3} - S c_{2} & 0 & 0 & c_{1}\\I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} & - \alpha & S \beta_{1} + S_{1} \beta_{1} c_{0} & S \beta_{2} + S_{1} \beta_{2} c_{0} & S \beta_{3} + S_{1} \beta_{3} c_{0} & 0 & 0 & c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right)\\0 & \alpha & - \gamma_{1} - p_{1} & 0 & 0 & 0 & 0 & 0\\0 & 0 & p_{1} & - \gamma_{2} - p_{2} & 0 & 0 & 0 & 0\\0 & 0 & 0 & p_{2} & - \gamma_{3} - \mu & 0 & 0 & 0\\0 & 0 & \gamma_{1} & \gamma_{2} & \gamma_{3} & 0 & 0 & 0\\0 & 0 & 0 & 0 & \mu & 0 & 0 & 0\\I_{3} c_{2} & 0 & - S_{1} \beta_{1} c_{0} & - S_{1} \beta_{2} c_{0} & S c_{2} - S_{1} \beta_{3} c_{0} & 0 & 0 & - c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) - c_{1}\end{matrix}\right]$




```python
SCEIR_model.get_grad_eqn()
```




$\displaystyle \left[\begin{array}{cccccccccccccc}- I_{1} S & - I_{2} S & - I_{3} S & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & S_{1} & - I_{3} S & 0\\I_{1} S + I_{1} S_{1} c_{0} & I_{2} S + I_{2} S_{1} c_{0} & I_{3} S + I_{3} S_{1} c_{0} & - E & 0 & 0 & 0 & 0 & 0 & 0 & S_{1} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) & 0 & 0 & 0\\0 & 0 & 0 & E & - I_{1} & 0 & 0 & - I_{1} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & - I_{2} & 0 & I_{1} & - I_{2} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - I_{3} & 0 & I_{2} & - I_{3} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & I_{1} & I_{2} & I_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & I_{3} & 0 & 0 & 0 & 0\\- I_{1} S_{1} c_{0} & - I_{2} S_{1} c_{0} & - I_{3} S_{1} c_{0} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - S_{1} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3}\right) & - S_{1} & I_{3} S & 0\end{array}\right]$



# SCCEIR model with caution distinguished $E$ and  𝐼1

## Equations

\begin{equation}
\begin{split}
\dot{S} &= -(\beta_1 (I_1 + c_0 I_{1c}) +\beta_2 I_2 + \beta_3 I_3) S + c_1 S_c - c_2*S*I_3\\
\dot{S_c} &= - c_0(\beta_1 (I_1 + c_0 I_{1c}) +\beta_2 I_2 + \beta_3 I_3) S_c - c_1 S_c + c_2*S*I_3\\
\dot{E} &=(\beta_1 (I_1 + c_0 I_{1c}) +\beta_2 I_2 + \beta_3 I_3 ) S - a E \\
\dot{E_c} &=(\beta_1 (I_1 + c_0 I_{1c}) +\beta_2 I_2 + \beta_3 I_3 ) c_0 S_c - a E \\
\dot{I_1} &= a E - \gamma_1 I_1 - p_1 I_1 \\
\dot{I_{1c}} &= a E_c - \gamma_1 I_{1c} - p_1 I_{1c} \\
\dot{I_2} &= p_1 (I_1 + I_{1c}) -\gamma_2 I_2 - p_2 I_2 \\
\dot{I_3} & = p_2 I_2 -\gamma_3 I_3 - \mu I_3 \\
\dot{R} & = \gamma_1 (I_1 + I_{1c}) + \gamma_2 I_2 + \gamma_3 I_3 \\
\dot{D} & = \mu I_3
\end{split}
\end{equation}

The use of I_3 as a state variable triggering susceptibles to execute caution is just one choice. In contrast with deaths, it does not accumulate over time and so retains the property of an active threat to society, rather than an historical one.

Here, the effect of caution is quadratic, since both the individual doing the infection and individual potentially being infected may be executing caution. To implement this we distinguish careful and non careful exposed and infectives up to the I_1 level. Once in hospital there is no difference, since all caution is executed wrt infected patients.

## Variables
* $S$: Susceptible individuals living as normal
* $S_c$: Susceptible individuals exercising pandemy precautions
* $E$: Exposed individuals living as normal - infected but not yet infectious or symptomatic
* $E_c$: Exposed individuals exercising pandemy precautions
* $I_i$: Infected individuals in severity class $i$. Severity increaes with $i$ and we assume individuals must pass through all previous classes. Split non hospital cases by caution.
  * $I_1$: Mild infection (hospitalization not required), living as normal
  * $I_{1c}$: Mild infection (hospitalization not required), exercising caution
  * $I_2$: Severe infection (hospitalization required)
  * $I_3$: Critical infection (ICU required)
* $R$: individuals who have recovered from disease and are now immune
* $D$: Dead individuals
* $N=S+S_c+E+E_c+I_{1c}+I_1+I_2+I_3+R+D$ Total population size (constant)

## Parameters
* $\beta_i$ rate at which infected individuals in class $I_i$ contact susceptibles and infect them
* $a$ rate of progression from the exposed to infected class
* $\gamma_i$ rate at which infected individuals in class $I_i$ recover from disease and become immune
* $p_i$ rate at which infected individuals in class $I_i$ progress to class $I_{I+1}$
* $\mu$ death rate for individuals in the most severe stage of disease

* $c_i$ three parameters characterizing cautionary response of population via class $S_c$

    - c_0 : reduction factor for exposure for cautioned susceptibles

    - c_1 : duration of caution (exponential decay time constant in days)

    - c_2 : rate constant for transition from uncautioned to cautioned susceptible

##  Implementation of SCCEIR model


```python
# import required packages
import os 
from sympy import symbols, init_printing
import numpy as np
import matplotlib.pyplot as plt
import sympy
import itertools
import scipy
from pygom import DeterministicOde, Transition, SimulateOde, TransitionType, SquareLoss
from scipy.optimize import minimize
from scipy.integrate import odeint
```


```python
# set up the symbolic SCCEIR model

state = ['S', 'E', 'I_1', 'I_2','I_3', 'R', 'D', 'S_c', 'E_c','I_c']
param_list = ['beta_1', 'beta_2','beta_3','alpha', 'gamma_1', 'gamma_2', 'gamma_3',
              'p_1','p_2','mu','c_0','c_1','c_2','N']

transition = [
    Transition(origin='S', destination='E', equation='(beta_1*I_1+beta_2*I_2+beta_3*I_3+c_0*beta_1*I_c)*S',
               transition_type=TransitionType.T),
    Transition(origin='S_c', destination='E_c', equation='c_0*(beta_1*I_1+beta_2*I_2+beta_3*I_3+c_0*beta_1*I_c)*S_c',
               transition_type=TransitionType.T),
    Transition(origin='S', destination='S_c', equation='c_2*I_3*S',
               transition_type=TransitionType.T),
    Transition(origin='S_c', destination='S', equation='c_1*S_c',
               transition_type=TransitionType.T),
    Transition(origin='E', destination='I_1', equation='alpha*E',
               transition_type=TransitionType.T),
    Transition(origin='E_c', destination='I_c', equation='alpha*E',
               transition_type=TransitionType.T),
    Transition(origin='I_1', destination='R', equation='gamma_1*I_1',
               transition_type=TransitionType.T),
    Transition(origin='I_c', destination='R', equation='gamma_1*I_c',
               transition_type=TransitionType.T),
    Transition(origin='I_2', destination='R', equation='gamma_2*I_2',
               transition_type=TransitionType.T),
    Transition(origin='I_3', destination='R', equation='gamma_3*I_3',
               transition_type=TransitionType.T),
    Transition(origin='I_1', destination='I_2', equation='p_1*I_1',
               transition_type=TransitionType.T),
    Transition(origin='I_c', destination='I_2', equation='p_1*I_c',
               transition_type=TransitionType.T),
    Transition(origin='I_2', destination='I_3', equation='p_2*I_2',
               transition_type=TransitionType.T),
    Transition(origin='I_3', destination='D', equation='mu*I_3',
               transition_type=TransitionType.T)    
    ]


SCCEIR_model = DeterministicOde(state, param_list, transition=transition)
```


```python
# display equations
print_ode2(SCCEIR_model)
```




$\displaystyle \left[\begin{matrix}d{S}/dt= & - I_{3} S c_{2} - S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) + S_{c} c_{1}\\d{E}/dt= & - E \alpha + S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right)\\d{I_1}/dt= & E \alpha - I_{1} \gamma_{1} - I_{1} p_{1}\\d{I_2}/dt= & I_{1} p_{1} - I_{2} \gamma_{2} - I_{2} p_{2} + I_{c} p_{1}\\d{I_3}/dt= & I_{2} p_{2} - I_{3} \gamma_{3} - I_{3} \mu\\d{R}/dt= & I_{1} \gamma_{1} + I_{2} \gamma_{2} + I_{3} \gamma_{3} + I_{c} \gamma_{1}\\d{D}/dt= & I_{3} \mu\\d{S_c}/dt= & I_{3} S c_{2} - S_{c} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) - S_{c} c_{1}\\d{E_c}/dt= & - E \alpha + S_{c} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right)\\d{I_c}/dt= & E \alpha - I_{c} \gamma_{1} - I_{c} p_{1}\end{matrix}\right]$




```python
# display graphical representation of the model
SCCEIR_model.get_transition_graph()
```


```python
SCCEIR_model.get_ode_eqn()
```




$\displaystyle \left[\begin{matrix}- I_{3} S c_{2} - S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) + S_{c} c_{1}\\- E \alpha + S \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right)\\E \alpha - I_{1} \gamma_{1} - I_{1} p_{1}\\I_{1} p_{1} - I_{2} \gamma_{2} - I_{2} p_{2} + I_{c} p_{1}\\I_{2} p_{2} - I_{3} \gamma_{3} - I_{3} \mu\\I_{1} \gamma_{1} + I_{2} \gamma_{2} + I_{3} \gamma_{3} + I_{c} \gamma_{1}\\I_{3} \mu\\I_{3} S c_{2} - S_{c} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) - S_{c} c_{1}\\- E \alpha + S_{c} c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right)\\E \alpha - I_{c} \gamma_{1} - I_{c} p_{1}\end{matrix}\right]$




```python
SCCEIR_model.get_jacobian_eqn()
```




$\displaystyle \left[\begin{matrix}- I_{1} \beta_{1} - I_{2} \beta_{2} - I_{3} \beta_{3} - I_{3} c_{2} - I_{c} \beta_{1} c_{0} & 0 & - S \beta_{1} & - S \beta_{2} & - S \beta_{3} - S c_{2} & 0 & 0 & c_{1} & 0 & - S \beta_{1} c_{0}\\I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0} & - \alpha & S \beta_{1} & S \beta_{2} & S \beta_{3} & 0 & 0 & 0 & 0 & S \beta_{1} c_{0}\\0 & \alpha & - \gamma_{1} - p_{1} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & p_{1} & - \gamma_{2} - p_{2} & 0 & 0 & 0 & 0 & 0 & p_{1}\\0 & 0 & 0 & p_{2} & - \gamma_{3} - \mu & 0 & 0 & 0 & 0 & 0\\0 & 0 & \gamma_{1} & \gamma_{2} & \gamma_{3} & 0 & 0 & 0 & 0 & \gamma_{1}\\0 & 0 & 0 & 0 & \mu & 0 & 0 & 0 & 0 & 0\\I_{3} c_{2} & 0 & - S_{c} \beta_{1} c_{0} & - S_{c} \beta_{2} c_{0} & S c_{2} - S_{c} \beta_{3} c_{0} & 0 & 0 & - c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) - c_{1} & 0 & - S_{c} \beta_{1} c_{0}^{2}\\0 & - \alpha & S_{c} \beta_{1} c_{0} & S_{c} \beta_{2} c_{0} & S_{c} \beta_{3} c_{0} & 0 & 0 & c_{0} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) & 0 & S_{c} \beta_{1} c_{0}^{2}\\0 & \alpha & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - \gamma_{1} - p_{1}\end{matrix}\right]$




```python
SCCEIR_model.get_grad_eqn()
```




$\displaystyle \left[\begin{array}{cccccccccccccc}- S \left(I_{1} + I_{c} c_{0}\right) & - I_{2} S & - I_{3} S & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - I_{c} S \beta_{1} & S_{c} & - I_{3} S & 0\\S \left(I_{1} + I_{c} c_{0}\right) & I_{2} S & I_{3} S & - E & 0 & 0 & 0 & 0 & 0 & 0 & I_{c} S \beta_{1} & 0 & 0 & 0\\0 & 0 & 0 & E & - I_{1} & 0 & 0 & - I_{1} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & - I_{2} & 0 & I_{1} + I_{c} & - I_{2} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - I_{3} & 0 & I_{2} & - I_{3} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & I_{1} + I_{c} & I_{2} & I_{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & I_{3} & 0 & 0 & 0 & 0\\- S_{c} c_{0} \left(I_{1} + I_{c} c_{0}\right) & - I_{2} S_{c} c_{0} & - I_{3} S_{c} c_{0} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - I_{c} S_{c} \beta_{1} c_{0} - S_{c} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) & - S_{c} & I_{3} S & 0\\S_{c} c_{0} \left(I_{1} + I_{c} c_{0}\right) & I_{2} S_{c} c_{0} & I_{3} S_{c} c_{0} & - E & 0 & 0 & 0 & 0 & 0 & 0 & I_{c} S_{c} \beta_{1} c_{0} + S_{c} \left(I_{1} \beta_{1} + I_{2} \beta_{2} + I_{3} \beta_{3} + I_{c} \beta_{1} c_{0}\right) & 0 & 0 & 0\\0 & 0 & 0 & E & - I_{c} & 0 & 0 & - I_{c} & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$



# Extract data from Johns Hopkins data base


```python
import numpy as np
import csv
import itertools
import matplotlib
%matplotlib inline
import seaborn as sb
from matplotlib import pyplot as plt
from cycler import cycler
import datetime
import matplotlib.dates as mdates
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
import pwlf
import sys
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
```


<style>.container { width:100% !important; }</style>



```python
def get_data(jhu_file):
    dat = []
    with open(jhu_file, newline='') as csvfile:
        myreader = csv.reader(csvfile, delimiter=',')
        popdat = []
        i = 0
        for row in myreader:
            if i != 0:
                poplist = []
                j = 0
                for elt in row:
                    if j >= 4:
                        poplist.append(int(elt))
                    elif j == 0:
                        poplist.append(elt)
                    elif j == 1:
                        poplist[0]=(elt,poplist[0])
                    j = j+1
                popdat.append(poplist)
            else:
                popdat.append(row)
            # print(popdat[i])
            i = i + 1;
    # dates
    popdat0=['dates']
    for elt in popdat[0][4:]:
        popdat0.append(elt)
    popdat[0] = [pop for pop in popdat0]
    # print('popdat[0]',popdat[0])
    # totals over all countries
    totals = np.zeros(len(popdat[0])-1,dtype=int)
    for row in popdat[1:]:
        totals = totals + np.array(row[1:])
    totals = list(np.asarray(totals))
    # print(totals)
    popkeyed = {poplist[0]: poplist[1:] for poplist in popdat}
    popkeyed.update({'dates':popdat[0][1:]})
    popkeyed.update({('World',''):totals})
    # del popkeyed[('d','a')]
    # assemble totals for countries with multiple regions
    total = np.zeros(len(popkeyed['dates']),dtype=int)
    poptotkeyed = {}
    for country,tseries in popkeyed.items():
        if country!='dates' and country[1] != '': # it seems that UK is single exception with both '' and non '' regions, UK total is then UK overseas
            countrytotal = (country[0],'Total')
            if countrytotal in poptotkeyed:
                # print(country,popkeyed[country],poptotkeyed[countrytotal])
                total = np.array(tseries)[:]+np.array(poptotkeyed[countrytotal])[:]
            else:
                total =  np.array(tseries)                        
            poptotkeyed.update({countrytotal:list(total)})
    for countrytotal,tseries in poptotkeyed.items():
        total = np.array(tseries)
        popkeyed.update({countrytotal:list(total)})
    return popkeyed
```

## JHU data


```python
base = '../../covid-19-JH/csse_covid_19_data/csse_covid_19_time_series/'
confirmed = get_data(base+'time_series_covid19_confirmed_global.csv')
deaths = get_data(base+'time_series_covid19_deaths_global.csv')
recovered = get_data(base+'time_series_covid19_recovered_global.csv')
covid_ts = {'confirmed':confirmed,'deaths':deaths,'recovered':recovered}
countries_jhu = [(row[0],row[1]) for row in confirmed][1:]
print("number of countries listed",len(countries_jhu))
i=0
for country in countries_jhu:
    print(i,country)
    i = i + 1
```

    number of countries listed 274
    0 ('Afghanistan', '')
    1 ('Albania', '')
    2 ('Algeria', '')
    3 ('Andorra', '')
    4 ('Angola', '')
    5 ('Antigua and Barbuda', '')
    6 ('Argentina', '')
    7 ('Armenia', '')
    8 ('Australia', 'Australian Capital Territory')
    9 ('Australia', 'New South Wales')
    10 ('Australia', 'Northern Territory')
    11 ('Australia', 'Queensland')
    12 ('Australia', 'South Australia')
    13 ('Australia', 'Tasmania')
    14 ('Australia', 'Victoria')
    15 ('Australia', 'Western Australia')
    16 ('Austria', '')
    17 ('Azerbaijan', '')
    18 ('Bahamas', '')
    19 ('Bahrain', '')
    20 ('Bangladesh', '')
    21 ('Barbados', '')
    22 ('Belarus', '')
    23 ('Belgium', '')
    24 ('Benin', '')
    25 ('Bhutan', '')
    26 ('Bolivia', '')
    27 ('Bosnia and Herzegovina', '')
    28 ('Brazil', '')
    29 ('Brunei', '')
    30 ('Bulgaria', '')
    31 ('Burkina Faso', '')
    32 ('Cabo Verde', '')
    33 ('Cambodia', '')
    34 ('Cameroon', '')
    35 ('Canada', 'Alberta')
    36 ('Canada', 'British Columbia')
    37 ('Canada', 'Grand Princess')
    38 ('Canada', 'Manitoba')
    39 ('Canada', 'New Brunswick')
    40 ('Canada', 'Newfoundland and Labrador')
    41 ('Canada', 'Nova Scotia')
    42 ('Canada', 'Ontario')
    43 ('Canada', 'Prince Edward Island')
    44 ('Canada', 'Quebec')
    45 ('Canada', 'Saskatchewan')
    46 ('Central African Republic', '')
    47 ('Chad', '')
    48 ('Chile', '')
    49 ('China', 'Anhui')
    50 ('China', 'Beijing')
    51 ('China', 'Chongqing')
    52 ('China', 'Fujian')
    53 ('China', 'Gansu')
    54 ('China', 'Guangdong')
    55 ('China', 'Guangxi')
    56 ('China', 'Guizhou')
    57 ('China', 'Hainan')
    58 ('China', 'Hebei')
    59 ('China', 'Heilongjiang')
    60 ('China', 'Henan')
    61 ('China', 'Hong Kong')
    62 ('China', 'Hubei')
    63 ('China', 'Hunan')
    64 ('China', 'Inner Mongolia')
    65 ('China', 'Jiangsu')
    66 ('China', 'Jiangxi')
    67 ('China', 'Jilin')
    68 ('China', 'Liaoning')
    69 ('China', 'Macau')
    70 ('China', 'Ningxia')
    71 ('China', 'Qinghai')
    72 ('China', 'Shaanxi')
    73 ('China', 'Shandong')
    74 ('China', 'Shanghai')
    75 ('China', 'Shanxi')
    76 ('China', 'Sichuan')
    77 ('China', 'Tianjin')
    78 ('China', 'Tibet')
    79 ('China', 'Xinjiang')
    80 ('China', 'Yunnan')
    81 ('China', 'Zhejiang')
    82 ('Colombia', '')
    83 ('Congo (Brazzaville)', '')
    84 ('Congo (Kinshasa)', '')
    85 ('Costa Rica', '')
    86 ("Cote d'Ivoire", '')
    87 ('Croatia', '')
    88 ('Diamond Princess', '')
    89 ('Cuba', '')
    90 ('Cyprus', '')
    91 ('Czechia', '')
    92 ('Denmark', 'Faroe Islands')
    93 ('Denmark', 'Greenland')
    94 ('Denmark', '')
    95 ('Djibouti', '')
    96 ('Dominican Republic', '')
    97 ('Ecuador', '')
    98 ('Egypt', '')
    99 ('El Salvador', '')
    100 ('Equatorial Guinea', '')
    101 ('Eritrea', '')
    102 ('Estonia', '')
    103 ('Eswatini', '')
    104 ('Ethiopia', '')
    105 ('Fiji', '')
    106 ('Finland', '')
    107 ('France', 'French Guiana')
    108 ('France', 'French Polynesia')
    109 ('France', 'Guadeloupe')
    110 ('France', 'Mayotte')
    111 ('France', 'New Caledonia')
    112 ('France', 'Reunion')
    113 ('France', 'Saint Barthelemy')
    114 ('France', 'St Martin')
    115 ('France', 'Martinique')
    116 ('France', '')
    117 ('Gabon', '')
    118 ('Gambia', '')
    119 ('Georgia', '')
    120 ('Germany', '')
    121 ('Ghana', '')
    122 ('Greece', '')
    123 ('Guatemala', '')
    124 ('Guinea', '')
    125 ('Guyana', '')
    126 ('Haiti', '')
    127 ('Holy See', '')
    128 ('Honduras', '')
    129 ('Hungary', '')
    130 ('Iceland', '')
    131 ('India', '')
    132 ('Indonesia', '')
    133 ('Iran', '')
    134 ('Iraq', '')
    135 ('Ireland', '')
    136 ('Israel', '')
    137 ('Italy', '')
    138 ('Jamaica', '')
    139 ('Japan', '')
    140 ('Jordan', '')
    141 ('Kazakhstan', '')
    142 ('Kenya', '')
    143 ('Korea, South', '')
    144 ('Kuwait', '')
    145 ('Kyrgyzstan', '')
    146 ('Latvia', '')
    147 ('Lebanon', '')
    148 ('Liberia', '')
    149 ('Liechtenstein', '')
    150 ('Lithuania', '')
    151 ('Luxembourg', '')
    152 ('Madagascar', '')
    153 ('Malaysia', '')
    154 ('Maldives', '')
    155 ('Malta', '')
    156 ('Mauritania', '')
    157 ('Mauritius', '')
    158 ('Mexico', '')
    159 ('Moldova', '')
    160 ('Monaco', '')
    161 ('Mongolia', '')
    162 ('Montenegro', '')
    163 ('Morocco', '')
    164 ('Namibia', '')
    165 ('Nepal', '')
    166 ('Netherlands', 'Aruba')
    167 ('Netherlands', 'Curacao')
    168 ('Netherlands', 'Sint Maarten')
    169 ('Netherlands', '')
    170 ('New Zealand', '')
    171 ('Nicaragua', '')
    172 ('Niger', '')
    173 ('Nigeria', '')
    174 ('North Macedonia', '')
    175 ('Norway', '')
    176 ('Oman', '')
    177 ('Pakistan', '')
    178 ('Panama', '')
    179 ('Papua New Guinea', '')
    180 ('Paraguay', '')
    181 ('Peru', '')
    182 ('Philippines', '')
    183 ('Poland', '')
    184 ('Portugal', '')
    185 ('Qatar', '')
    186 ('Romania', '')
    187 ('Russia', '')
    188 ('Rwanda', '')
    189 ('Saint Lucia', '')
    190 ('Saint Vincent and the Grenadines', '')
    191 ('San Marino', '')
    192 ('Saudi Arabia', '')
    193 ('Senegal', '')
    194 ('Serbia', '')
    195 ('Seychelles', '')
    196 ('Singapore', '')
    197 ('Slovakia', '')
    198 ('Slovenia', '')
    199 ('Somalia', '')
    200 ('South Africa', '')
    201 ('Spain', '')
    202 ('Sri Lanka', '')
    203 ('Sudan', '')
    204 ('Suriname', '')
    205 ('Sweden', '')
    206 ('Switzerland', '')
    207 ('Taiwan*', '')
    208 ('Tanzania', '')
    209 ('Thailand', '')
    210 ('Togo', '')
    211 ('Trinidad and Tobago', '')
    212 ('Tunisia', '')
    213 ('Turkey', '')
    214 ('Uganda', '')
    215 ('Ukraine', '')
    216 ('United Arab Emirates', '')
    217 ('United Kingdom', 'Bermuda')
    218 ('United Kingdom', 'Cayman Islands')
    219 ('United Kingdom', 'Channel Islands')
    220 ('United Kingdom', 'Gibraltar')
    221 ('United Kingdom', 'Isle of Man')
    222 ('United Kingdom', 'Montserrat')
    223 ('United Kingdom', '')
    224 ('Uruguay', '')
    225 ('US', '')
    226 ('Uzbekistan', '')
    227 ('Venezuela', '')
    228 ('Vietnam', '')
    229 ('Zambia', '')
    230 ('Zimbabwe', '')
    231 ('Canada', 'Diamond Princess')
    232 ('Dominica', '')
    233 ('Grenada', '')
    234 ('Mozambique', '')
    235 ('Syria', '')
    236 ('Timor-Leste', '')
    237 ('Belize', '')
    238 ('Laos', '')
    239 ('Libya', '')
    240 ('West Bank and Gaza', '')
    241 ('Guinea-Bissau', '')
    242 ('Mali', '')
    243 ('Saint Kitts and Nevis', '')
    244 ('Canada', 'Northwest Territories')
    245 ('Canada', 'Yukon')
    246 ('Kosovo', '')
    247 ('Burma', '')
    248 ('United Kingdom', 'Anguilla')
    249 ('United Kingdom', 'British Virgin Islands')
    250 ('United Kingdom', 'Turks and Caicos Islands')
    251 ('MS Zaandam', '')
    252 ('Botswana', '')
    253 ('Burundi', '')
    254 ('Sierra Leone', '')
    255 ('Netherlands', 'Bonaire, Sint Eustatius and Saba')
    256 ('Malawi', '')
    257 ('United Kingdom', 'Falkland Islands (Malvinas)')
    258 ('France', 'Saint Pierre and Miquelon')
    259 ('South Sudan', '')
    260 ('Western Sahara', '')
    261 ('Sao Tome and Principe', '')
    262 ('Yemen', '')
    263 ('Comoros', '')
    264 ('Tajikistan', '')
    265 ('Lesotho', '')
    266 ('World', '')
    267 ('Australia', 'Total')
    268 ('Canada', 'Total')
    269 ('China', 'Total')
    270 ('Denmark', 'Total')
    271 ('France', 'Total')
    272 ('Netherlands', 'Total')
    273 ('United Kingdom', 'Total')


## Get data for one country


```python
def Float(x):
    try:
        rtn = float(x)
    except:
        rtn = float('NaN')
    return rtn
```


```python
 def get_country_data(country_s='World', datatype='confirmed', firstdate=None, lastdate=None):
    if isinstance(country_s,str):
        country = (country_s,'')
    else:                               # single ('country','reg') entry
        country = country_s
    popkeyed = covid_ts[datatype]
    
    dates = popkeyed['dates']
    fmt = '%m/%d/%y'
    xx = [datetime.datetime.strptime(dd,fmt) for dd in dates ]
    if firstdate:
        firstdate_d = datetime.datetime.strptime(firstdate,fmt)
    else:
        firstdate_d = datetime.datetime.strptime(dates[0],fmt)
    if lastdate:
        lastdate_d = datetime.datetime.strptime(lastdate,fmt)
    else:
        lastdate_d = datetime.datetime.strptime(dates[-1],fmt)    
    daystart = (firstdate_d-xx[0]).days
    daystop = (lastdate_d-xx[-1]).days

    try:
        yy = popkeyed[country]
        # print(country)
    except:
            print('country data not found',country)
            return None,None      
    yyf = [Float(y) for y in yy]

    yy0 = []
    yy1 = []  
    if daystart>len(yyf):
        print('Error: start date does not overlap with available data')
        return None,None
    elif daystart>0:
        yyf = yyf[daystart:]
    elif daystart <0:
        yy0 = [0.]*(-daystart)
        
    if daystop < 0:
        yyf = yyf[:daystop]  
    elif daystop > 0:
        yy1 = [0.]*(daystop)
    yyf = yy0 + yyf + yy1
    xxf = [float(x) for x in range(len(yyf))]
    return xxf,yyf 
```


```python
# Italy assumed data starting on firstdate
N = 60000000
firstdate = '01/25/20'
lastdate = '06/07/20'
xxf,yy = get_country_data('Italy','confirmed',firstdate=firstdate,lastdate=lastdate)
# print(xxf)
xxf,yy1 = get_country_data('Italy','recovered',firstdate=firstdate,lastdate=lastdate)
xxf,yy2 = get_country_data('Italy','deaths',firstdate=firstdate,lastdate=lastdate)
y_jhu = np.array([[yy1[i],yy2[i]] for i in range(0,len(yy))])/N
data = np.array([[xxf[i],yy1[i],yy2[i]] for i in range(len(yy))])
# print(data)
lastday = len(y_jhu)
print(lastday)
```

    135



```python
plt.plot(data)
```

## Data using jhu module


```python
import sys
sys.path = ['..']+sys.path
import jhu
import matplotlib
%matplotlib inline
from matplotlib import pyplot as plt
```


```python
dat = jhu.get_ave_data('Italy')
```


```python
print(len(dat['Italy']))
```

    138



```python
plt.plot(range(len(dat['dates'])),dat['Italy']);
```


```python
dat = jhu.get_ave_data('Italy',datatype='deaths',dataaccum='daily',firstdate='02/10/20',lastdate='06/07/20')
```


```python
print(len(dat['Italy']))
```

    138



```python
plt.plot(range(len(dat['dates'])),dat['Italy']);
```

# Simulation of SCEIR model


```python
# setup time points for simulation, initial conditions and parameters
t = np.linspace(0, lastday -1, lastday)

# initial conditions assuming there is no natural immunity
E_0 =  0.00003
x0 = [1.0-E_0, E_0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

# Define parameters based on clinical observations Dr. Alison
Exposure=0.25 # Rate coefficient for exposure per individual in contact per day
IncubPeriod=5  #Incubation period, days
DurMildInf=10 #Duration of mild infections, days
FracMild=0.70  #Fraction of infections that are mild
FracSevere=0.20 #Fraction of infections that are severe
FracCritical=0.1 #Fraction of infections that are critical
CFR=0.05 #Case fatality rate (fraction of infections resulting in death)
TimeICUDeath=7 #Time from ICU admission to death, days
DurHosp=11 #Duration of hospitalization, days

# Model extension by John McCaskill to include caution 
CautionFactor= 0.4 # Fractional reduction of exposure rate for cautioned individuals
CautionRetention= 21. # Duration of cautionary state of susceptibles (2 weeks)
CautionICUFrac= 0.25 # Fraction of ICUs occupied leading to transition to caution @ 1/day
ICUFrac= 0.002 # Fraction of ICUs relative to population size N

# Model fitting extension to allow for incomplete detection
FracConfirmedDet=0.125  # Fraction of recovered individuals measured
FracRecoveredDet=FracConfirmedDet # Fraction of recovered individuals measured
FracDeathsDet=1.0

# convert above parameters to model parameters
params = {'beta_1' : Exposure/sum(x0),
          'beta_2' : 0.0,
          'beta_3' : 0.0,
          'alpha' : 1.0/IncubPeriod,
          'gamma_1': (1.0/DurMildInf)*FracMild,
          'gamma_2': (1.0/DurHosp)-(1/DurHosp)*(FracCritical/(FracCritical+FracSevere)),
          'gamma_3': (1.0/TimeICUDeath)-(1/TimeICUDeath)*(CFR/FracCritical),
          'p_1'    : (1.0/DurMildInf)-(1.0/DurMildInf)*FracMild,
          'p_2'    : (1.0/DurHosp)*(FracCritical/(FracCritical+FracSevere)),
          'mu'    : (1.0/TimeICUDeath)*(CFR/FracCritical),
          'c_0'    : CautionFactor,
          'c_1'    : 1.0/CautionRetention,
          'c_2'    : 1.0/(sum(x0)*ICUFrac*CautionICUFrac),
          'N'     : sum(x0)}
print(params)
# assign x0 and params to the model, integrate over t and plot
SCEIR_model.initial_values = (x0, t[0])
SCEIR_model.parameters = params
solution = SCEIR_model.integrate(t[1::])

SCEIR_model.plot()

# calculate time point when maximum number of people are infectious
peak_i = np.argmax(solution[:,2])
print('Peak infection (days)', t[peak_i])
```

## Compare data with SCEIR simulation


```python
# model with generating parameters 
print(params)
params_fit = params.copy()
print(params_fit['c_0'],params_fit['c_1'])
SCEIR_model.params = params_fit

x0_fit = x0.copy()
# x0_fit[2] = 0.00001
#t_fit = numpy.linspace(0, 150, 1000)
print(x0_fit)
t_fit = t
print(len(t))
SCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
#sol_fit =SCEIR_model.integrate(t_fit[1::])
# sol_fit = SCEIR_model.integrate(t_fit)
sol_fit = scipy.integrate.odeint(SCEIR_model.ode, x0_fit, t_fit[1::])
# print(len(sol_fit[0]))
#
plt.figure(figsize=(15,10))
#plt.plot(t,y_jhu[:,0], 'go',label='I_1') # infected observations
#plt.plot(t,y_jhu[:,1], 'go',label='I_2') # infected observations
#plt.plot(t,y_jhu[:,2], 'go',label='I_3') # infected observations
plt.plot(t,y_jhu[:,0]/FracRecoveredDet, 'bo',label='R') # recovered
plt.plot(t,y_jhu[:,1], 'ro',label='D') # died
plt.gca().set_prop_cycle(color=['grey','orange','green','green','green','blue','red', 'black'])
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.06])
plt.legend()
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

## Integration and plot using scipy and matplotlib directly


```python
# solution = scipy.integrate.odeint(SCEIR_model.ode, x0, t)
# print(len(t))
solution1 = scipy.integrate.odeint(SCEIR_model.ode, x0, t[1::])
ys = solution1.copy()
plt.figure(figsize=(15,7))
plt.subplot(1,2,1)
plt.xlabel("Time (days)")
plt.ylabel("Number per 1000000 People")
plt.ylim([0,1])
plt.legend(("S","E","I1","I2","I3","R","D","Sc"))
plt.plot(t[1::],ys)
plt.subplot(1,2,2)
plt.xlabel("Time (days)")
plt.ylabel("Number per 1000000 People")
plt.ylim([0.000001,1])
plt.legend(("S","E","I1","I2","I3","R","D","Sc"))
plt.semilogy()
plt.plot(t[1::],ys)
```

# Simulation of SCCEIR model


```python
# setup time points for simulation, initial conditions and parameters
t = np.linspace(0, lastday -1, lastday)

# initial conditions assuming there is no natural immunity
E_0 =  0.00003
x0_c = [1.0-E_0, E_0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

# Define parameters based on clinical observations Dr. Alison
Exposure=0.25 # Rate coefficient for exposure per individual in contact per day
IncubPeriod=5  #Incubation period, days
DurMildInf=10 #Duration of mild infections, days
FracMild=0.70  #Fraction of infections that are mild
FracSevere=0.20 #Fraction of infections that are severe
FracCritical=0.1 #Fraction of infections that are critical
CFR=0.07 #Case fatality rate (fraction of infections resulting in death)
TimeICUDeath=7 #Time from ICU admission to death, days
DurHosp=11 #Duration of hospitalization, days

# Model fitting extension to allow for incomplete detection
FracConfirmedDet=0.125  # Fraction of recovered individuals measured
FracRecoveredDet=FracConfirmedDet # Fraction of recovered individuals measured
FracDeathsDet=1.0

# Model extension by John McCaskill to include caution 
CautionFactor= 0.2 # Fractional reduction of exposure rate for cautioned individuals
CautionRetention= 21. # Duration of cautionary state of susceptibles (2 weeks)
CautionICUFrac= 0.3 # Fraction of ICUs occupied leading to transition to caution @ 1/day
ICUFrac= 0.002 # Fraction of ICUs relative to population size N


# convert above parameters to model parameters
params = {'beta_1' : Exposure/sum(x0),
          'beta_2' : 0.0,
          'beta_3' : 0.0,
          'alpha' : 1.0/IncubPeriod,
          'gamma_1': (1.0/DurMildInf)*FracMild,
          'gamma_2': (1.0/DurHosp)-(1/DurHosp)*(FracCritical/(FracCritical+FracSevere)),
          'gamma_3': (1.0/TimeICUDeath)-(1/TimeICUDeath)*(CFR/FracCritical),
          'p_1'    : (1.0/DurMildInf)-(1.0/DurMildInf)*FracMild,
          'p_2'    : (1.0/DurHosp)*(FracCritical/(FracCritical+FracSevere)),
          'mu'    : (1.0/TimeICUDeath)*(CFR/FracCritical),
          'c_0'    : CautionFactor,
          'c_1'    : 1.0/CautionRetention,
          'c_2'    : 1.0/(sum(x0)*ICUFrac*CautionICUFrac),
          'N'     : sum(x0)}
print(params)
# assign x0 and params to the model, integrate over t and plot
SCCEIR_model.initial_values = (x0_c, t[0])
SCCEIR_model.parameters = params
solution = SCCEIR_model.integrate(t[1::])

SCCEIR_model.plot()

# calculate time point when maximum number of people are infectious
peak_i = np.argmax(solution[:,2])
print('Peak infection (days)', t[peak_i])
```

## Compare data with SCCEIR simulation


```python
# model with generating parameters 
print(params)
params_fit = params.copy()
print(params_fit['c_0'],params_fit['c_1'])
SCCEIR_model.params = params_fit

x0_fit = x0_c.copy()
# x0_fit[2] = 0.00001
#t_fit = numpy.linspace(0, 150, 1000)
print(x0_fit)
t_fit = t
print(len(t))
SCCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
#sol_fit =SCCEIR_model.integrate(t_fit[1::])
# sol_fit = SCCEIR_model.integrate(t_fit)
sol_fit = scipy.integrate.odeint(SCCEIR_model.ode, x0_fit, t_fit[1::])
# print(len(sol_fit[0]))
#
plt.figure(figsize=(15,10))
plt.plot(t,y_jhu[:,0]/FracRecoveredDet, 'bo',label='R') # recovered
plt.plot(t,y_jhu[:,1]/FracDeathsDet, 'ro',label='D') # died
plt.gca().set_prop_cycle(color=['grey','orange','green','green','green','blue','red', 'black'])
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.06])
plt.legend()
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

# Parameter fitting

## Generate test data based on SCEIR simulation


```python
# Add noise
y = solution[:,2:7].copy()
#print('len(y)',len(y),'t',len(t),t[0],t[1],'...',t[-1])
np.random.seed(seed=6)
noise = np.random.normal(0,1.e-2,[len(t),5])
# ynoise = y *(1+noise)
ynoise = y *(1.0 + noise)
ynoise[ynoise<0] = 0
plt.figure(figsize=(15,10))
plt.plot(t,ynoise[:,0], 'go', label='I_1') 
plt.plot(t,ynoise[:,1], 'go', label='I_2') 
plt.plot(t,ynoise[:,2], 'go', label='I_3') 
plt.plot(t,ynoise[:,3], 'bo', label='R') 
plt.plot(t,ynoise[:,4], 'ro', label='D') 
plt.plot(t,y[:,0], 'g', label='I_1') 
plt.plot(t,y[:,1], 'g', label='I_2') 
plt.plot(t,y[:,2], 'g', label='I_3') 
plt.plot(t,y[:,3], 'b', label='R') 
plt.plot(t,y[:,4], 'r', label='D') 
plt.legend()
plt.ylim(0,0.003)
plt.show()
```


```python
# model with generating parameters 
print(params)
params_fit = params.copy()
print(params_fit['c_0'],params_fit['c_1'])
SCEIR_model.params = params_fit

x0_fit = x0.copy()
# x0_fit[2] = 0.00001
#t_fit = numpy.linspace(0, 150, 1000)
print(x0_fit)
t_fit = t
print(len(t))
SCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
#sol_fit =SCEIR_model.integrate(t_fit[1::])
# sol_fit = SCEIR_model.integrate(t_fit)
sol_fit = scipy.integrate.odeint(SCEIR_model.ode, x0_fit, t_fit[1::])
# print(len(sol_fit[0]))
#
plt.figure(figsize=(15,10))
plt.plot(t,ynoise[:,0], 'go',label='I_1') # infected observations
plt.plot(t,ynoise[:,1], 'go',label='I_2') # infected observations
plt.plot(t,ynoise[:,2], 'go',label='I_3') # infected observations
plt.plot(t,ynoise[:,3], 'bo',label='R') # recoverd
plt.plot(t,ynoise[:,4], 'ro',label='D') # died
plt.gca().set_prop_cycle(color=['grey','orange','green','green','green','blue','red', 'black'])
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.004])
plt.legend()
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```


```python
params # use list(...) to convert to list
```

## Fit parameters to randomized simulation data based on square_loss

### Fit c_0 and c_1 only


```python
# Initial guess of parameters, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [value for value in cautionparams]
theta = [0.21,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.0,1.0],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'])
# perform optimization
res = minimize(fun=objSCEIR.cost,
                jac=objSCEIR.sensitivity,
                x0=theta,
                #bounds=boxBounds,
                method='BFGS',
                options={'disp':True,'maxiter':1000,'eps':0.00001})# ,'ftol':0.01})  #not BFGS
print(res)
```

    Optimization terminated successfully.
             Current function value: 0.000575
             Iterations: 28
             Function evaluations: 32
             Gradient evaluations: 32
          fun: 0.0005749380626160615
     hess_inv: array([[ 4.09498456e+10, -1.40156972e+10],
           [-1.40156972e+10,  4.79708205e+09]])
          jac: array([-2.65189056e-06, -7.75111654e-06])
      message: 'Optimization terminated successfully.'
         nfev: 32
          nit: 28
         njev: 32
       status: 0
      success: True
            x: array([-307.13613516,  105.2169699 ])


### Fit c_0 and c_1 as well as initial value of E

#### Fit c_0 and c_1 as well as initial value of E with 'SLSQP'
does not work well
note use of special methods IV for initial value fitting of target_state


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='SLSQP',
                options={'disp':True,'maxiter':1000,'eps':0.01,'ftol':0.01})
print(res)
```

#### Fit c_0 and c_1 as well as initial value of E with BFGS
works well:  no constraints and gtol not ftol


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='BFGS',
                options={'disp':True,'maxiter':1000,'eps':0.01,'gtol':0.01})
print(res)
```


```python
# model with fitted parameters 
print(params)
print(x0)
params_fit = params.copy()
#params_fit['c_0'] = res.x[0]
#params_fit['c_1'] = res.x[1]

SCEIR_model.params = params_fit
print(SCEIR_model.params)
#ode_fit = common_models.SEIR({'beta':res.x[0], 'gamma':res.x[1],'alpha':res.x[2]})
#x0_fit = [1-1.27e-6, 1.27e-6, 0]
x0_fit = x0.copy()
#x0_fit[2] = res.x[2]
#t_fit = numpy.linspace(0, 150, 1000)
t_fit = t
SCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
# sol_fit =SCEIR_model.integrate(t_fit[0::])
sol_fit = scipy.integrate.odeint(SCEIR_model.ode, x0_fit, t_fit[1::])
#
plt.figure(figsize=(15,10))
plt.plot(t,ynoise[:,0], 'go') # infected observations
plt.plot(t,ynoise[:,1], 'go') # infected observations
plt.plot(t,ynoise[:,2], 'go') # infected observations
plt.plot(t,ynoise[:,3], 'bo') # recoverd
plt.plot(t,ynoise[:,4], 'ro') # died
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.004])
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

#### Fit c_0 and c_1 as well as initial value of E using L-BFGS-B
this method does not work well


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                bounds=boxBounds+[(0.0000001,0.001)],
                method='L-BFGS-B',
                options={'disp':True,'maxiter':1000,'eps':0.0001,'ftol':0.001})
print(res)
```

          fun: 0.000906012800877765
     hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
          jac: array([ 3.96334113e-03,  9.90926438e-03, -1.71024414e+01])
      message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
         nfev: 3
          nit: 1
       status: 0
      success: True
            x: array([2.44936963e-01, 6.88053051e-02, 3.13794907e-05])



```python
objSCEIR.residual()
```




    array([[-1.82936636e-07,  7.75738622e-08,  7.41002218e-10,
             1.90298452e-07,  2.61930903e-11],
           [-1.18714753e-07,  2.74592417e-07,  5.35455753e-09,
             7.02695193e-07,  3.73380566e-10],
           [-3.72984889e-07,  5.52824485e-07,  1.61628364e-08,
             1.46992842e-06,  1.73557012e-09],
           [-2.55156480e-07,  9.31197976e-07,  3.35222264e-08,
             2.43023007e-06,  5.01317221e-09],
           [-3.38990961e-07,  1.33814237e-06,  6.17848425e-08,
             3.66948629e-06,  1.11849951e-08],
           [ 4.43165958e-08,  1.69064251e-06,  9.74226592e-08,
             5.23887532e-06,  2.12464163e-08],
           [ 5.61800847e-07,  2.14049548e-06,  1.36856345e-07,
             6.76924192e-06,  3.65231556e-08],
           [ 1.44505990e-06,  2.74876078e-06,  1.88522242e-07,
             8.94973917e-06,  5.76158162e-08],
           [ 1.31880668e-06,  3.45150566e-06,  2.48565200e-07,
             1.08819831e-05,  8.49385764e-08],
           [ 2.60157903e-06,  4.10466029e-06,  3.21266928e-07,
             1.36379627e-05,  1.22436089e-07],
           [ 4.23480521e-06,  4.87106609e-06,  3.99530771e-07,
             1.65639519e-05,  1.72400885e-07],
           [ 4.63557060e-06,  5.43948086e-06,  4.93748009e-07,
             1.99050799e-05,  2.28784099e-07],
           [ 5.33624709e-06,  6.43738076e-06,  5.85779662e-07,
             2.23289600e-05,  2.97339293e-07],
           [ 7.61603426e-06,  7.39853035e-06,  7.31887474e-07,
             2.82904108e-05,  3.81306104e-07],
           [ 8.47031020e-06,  8.38645735e-06,  8.16229664e-07,
             3.19811038e-05,  4.76078665e-07],
           [ 1.08234655e-05,  9.74994616e-06,  9.42060075e-07,
             3.77058485e-05,  5.85754838e-07],
           [ 1.34121714e-05,  1.08878836e-05,  1.11708258e-06,
             4.36942113e-05,  7.14505329e-07],
           [ 1.57463055e-05,  1.21594019e-05,  1.30269821e-06,
             4.96253509e-05,  8.99964635e-07],
           [ 1.87300174e-05,  1.38702429e-05,  1.55419900e-06,
             5.68518451e-05,  1.04399249e-06],
           [ 2.05051355e-05,  1.56795264e-05,  1.72833848e-06,
             6.41939491e-05,  1.21117951e-06],
           [ 2.48025685e-05,  1.76191392e-05,  1.95063940e-06,
             7.28550461e-05,  1.44865342e-06],
           [ 2.84511788e-05,  2.01945427e-05,  2.23607797e-06,
             8.18774695e-05,  1.74364519e-06],
           [ 3.29605170e-05,  2.29912273e-05,  2.51115898e-06,
             9.55037117e-05,  2.02035738e-06],
           [ 3.89098895e-05,  2.48037620e-05,  2.85393660e-06,
             1.07823920e-04,  2.31898176e-06],
           [ 4.36400627e-05,  2.71782751e-05,  3.18462224e-06,
             1.20908517e-04,  2.68406730e-06],
           [ 4.97362573e-05,  3.08677234e-05,  3.64111750e-06,
             1.37169504e-04,  3.11234486e-06],
           [ 5.78335647e-05,  3.43295587e-05,  4.04737358e-06,
             1.56111775e-04,  3.47573280e-06],
           [ 6.62695159e-05,  3.81426211e-05,  4.27325433e-06,
             1.75924765e-04,  4.00142401e-06],
           [ 7.53053930e-05,  4.20932408e-05,  5.07101550e-06,
             1.97771440e-04,  4.58601108e-06],
           [ 8.04468932e-05,  4.74063681e-05,  5.81006475e-06,
             2.23031018e-04,  5.37573011e-06],
           [ 9.37568410e-05,  5.30838626e-05,  6.45123714e-06,
             2.47326470e-04,  5.90456070e-06],
           [ 1.05909359e-04,  5.79523344e-05,  6.97745486e-06,
             2.74308869e-04,  6.74530615e-06],
           [ 1.18264334e-04,  6.48780654e-05,  7.81251728e-06,
             3.04898892e-04,  7.55763747e-06],
           [ 1.27441639e-04,  7.27339771e-05,  8.79453244e-06,
             3.37796245e-04,  8.53916294e-06],
           [ 1.40130608e-04,  7.85460518e-05,  9.47647504e-06,
             3.81255143e-04,  9.60283997e-06],
           [ 1.59508710e-04,  8.79868292e-05,  1.07953368e-05,
             4.12681282e-04,  1.06296170e-05],
           [ 1.75768709e-04,  9.65783119e-05,  1.16408967e-05,
             4.57846725e-04,  1.22892962e-05],
           [ 1.88957184e-04,  1.07828939e-04,  1.32337049e-05,
             5.21151475e-04,  1.35201794e-05],
           [ 2.03405918e-04,  1.12365387e-04,  1.39034388e-05,
             5.77163744e-04,  1.52716203e-05],
           [ 2.20473490e-04,  1.27403946e-04,  1.56061474e-05,
             6.40618317e-04,  1.69449357e-05],
           [ 2.50579791e-04,  1.33640293e-04,  1.79026413e-05,
             6.98393715e-04,  1.88210264e-05],
           [ 2.53875322e-04,  1.49879205e-04,  1.93165795e-05,
             7.49976191e-04,  2.04322188e-05],
           [ 2.70234204e-04,  1.60763932e-04,  2.14243128e-05,
             8.44360943e-04,  2.36150188e-05],
           [ 2.98327842e-04,  1.71917897e-04,  2.27915910e-05,
             8.92681302e-04,  2.54917985e-05],
           [ 3.22127918e-04,  1.89206329e-04,  2.49759504e-05,
             9.99337666e-04,  2.85671999e-05],
           [ 3.35395024e-04,  2.01714229e-04,  2.70094774e-05,
             1.07336735e-03,  3.15503351e-05],
           [ 3.43704780e-04,  2.16635465e-04,  3.05116134e-05,
             1.16549341e-03,  3.43923415e-05],
           [ 3.53577014e-04,  2.32564784e-04,  3.21619875e-05,
             1.27532154e-03,  3.80529703e-05],
           [ 3.67235227e-04,  2.49338465e-04,  3.43097566e-05,
             1.40076981e-03,  4.15592285e-05],
           [ 3.81884745e-04,  2.69072579e-04,  3.89015042e-05,
             1.48615703e-03,  4.47407693e-05],
           [ 3.71290473e-04,  2.88380550e-04,  4.13543363e-05,
             1.64295841e-03,  5.05259643e-05],
           [ 3.90872601e-04,  2.95251928e-04,  4.27995434e-05,
             1.69743867e-03,  5.52098613e-05],
           [ 3.83331132e-04,  3.15449884e-04,  4.63774750e-05,
             1.83066777e-03,  5.95170374e-05],
           [ 3.60388395e-04,  3.37883058e-04,  4.94284640e-05,
             1.98290059e-03,  6.60606423e-05],
           [ 3.51088442e-04,  3.41351688e-04,  5.34061632e-05,
             2.09140247e-03,  7.08564691e-05],
           [ 3.21501026e-04,  3.61722888e-04,  5.49348658e-05,
             2.23677750e-03,  7.85744829e-05],
           [ 2.85701498e-04,  3.60725091e-04,  5.91493168e-05,
             2.36406027e-03,  8.39737485e-05],
           [ 2.70892905e-04,  3.79973978e-04,  6.11356978e-05,
             2.42160900e-03,  8.98514349e-05],
           [ 2.23787105e-04,  3.84664602e-04,  6.31265952e-05,
             2.60940148e-03,  9.73952174e-05],
           [ 1.75100508e-04,  4.08127802e-04,  6.71815163e-05,
             2.77583979e-03,  1.02005534e-04],
           [ 1.20964417e-04,  4.06972433e-04,  6.85920745e-05,
             2.87840443e-03,  1.11261662e-04],
           [ 5.85723397e-05,  4.09544162e-04,  7.14124333e-05,
             2.97813954e-03,  1.20036806e-04],
           [ 1.94499235e-05,  4.15352102e-04,  7.43461610e-05,
             3.13947904e-03,  1.28601932e-04],
           [-6.39338453e-05,  4.24528306e-04,  7.54744263e-05,
             3.26517538e-03,  1.35925687e-04],
           [-1.48793907e-04,  4.01292356e-04,  7.65437892e-05,
             3.27383330e-03,  1.47353761e-04],
           [-2.11258797e-04,  4.01961361e-04,  8.02038268e-05,
             3.41894822e-03,  1.53971954e-04],
           [-2.51883272e-04,  4.01515417e-04,  7.92776544e-05,
             3.50986623e-03,  1.66816676e-04],
           [-3.51480324e-04,  3.72044343e-04,  7.97068831e-05,
             3.63906859e-03,  1.74807650e-04],
           [-4.64082742e-04,  3.78033633e-04,  8.22837202e-05,
             3.77374278e-03,  1.80653631e-04],
           [-5.29767130e-04,  3.59506448e-04,  8.09727254e-05,
             3.59259644e-03,  1.94963649e-04],
           [-6.15957852e-04,  3.57378857e-04,  7.76915281e-05,
             3.82452480e-03,  1.99436553e-04],
           [-6.99145261e-04,  3.34005520e-04,  8.11605180e-05,
             3.94795448e-03,  2.05809354e-04],
           [-7.86365115e-04,  3.14261810e-04,  7.85870990e-05,
             3.98986859e-03,  2.19435182e-04],
           [-8.77515268e-04,  3.00728783e-04,  7.43073493e-05,
             3.98852021e-03,  2.28098793e-04],
           [-9.46485549e-04,  2.78584854e-04,  7.33947805e-05,
             3.93227804e-03,  2.39143614e-04],
           [-1.04156025e-03,  2.43999758e-04,  7.11799750e-05,
             3.90611675e-03,  2.45683275e-04],
           [-1.12837854e-03,  2.27522768e-04,  6.61722079e-05,
             4.03218973e-03,  2.62152609e-04],
           [-1.21014787e-03,  1.81872608e-04,  6.19449812e-05,
             4.05374941e-03,  2.69289773e-04],
           [-1.27359682e-03,  1.61320599e-04,  6.12399186e-05,
             3.87635526e-03,  2.78311161e-04],
           [-1.34303803e-03,  1.41033084e-04,  5.94682930e-05,
             3.89216762e-03,  2.80510664e-04],
           [-1.40816049e-03,  1.20698361e-04,  5.19004894e-05,
             3.86030302e-03,  2.92434579e-04],
           [-1.47224673e-03,  8.00658935e-05,  5.04144931e-05,
             4.08932646e-03,  3.00519592e-04],
           [-1.54276305e-03,  7.07091740e-05,  4.25244883e-05,
             3.95796464e-03,  3.12995295e-04],
           [-1.60822698e-03,  3.23977173e-05,  4.18092548e-05,
             3.76612104e-03,  3.04941946e-04],
           [-1.65778071e-03,  1.15120127e-05,  3.33096425e-05,
             3.82857041e-03,  3.16258663e-04],
           [-1.70717445e-03, -3.32477382e-05,  2.95769700e-05,
             3.71098512e-03,  3.28947468e-04],
           [-1.75756628e-03, -5.70049586e-05,  2.43475557e-05,
             3.55193590e-03,  3.31838359e-04],
           [-1.80094107e-03, -8.25704764e-05,  2.22579880e-05,
             3.50869704e-03,  3.41415007e-04],
           [-1.84348573e-03, -1.14471920e-04,  1.59229362e-05,
             3.37729656e-03,  3.47549736e-04],
           [-1.87164346e-03, -1.40985954e-04,  1.00518517e-05,
             3.38208721e-03,  3.53128946e-04],
           [-1.91670533e-03, -1.63051958e-04,  6.93562988e-07,
             3.21555925e-03,  3.54757705e-04],
           [-1.94653075e-03, -1.87285854e-04, -4.04193060e-06,
             3.02448182e-03,  3.49386281e-04],
           [-1.97328409e-03, -2.16152782e-04, -6.41565654e-06,
             2.89707270e-03,  3.63735739e-04],
           [-2.00070961e-03, -2.36317575e-04, -1.21579504e-05,
             3.00651312e-03,  3.71397140e-04],
           [-2.03121679e-03, -2.58374016e-04, -1.96222083e-05,
             2.75479033e-03,  3.60516433e-04],
           [-2.04191837e-03, -2.79375497e-04, -2.35756295e-05,
             2.48981971e-03,  3.57776124e-04],
           [-2.07179159e-03, -3.05586876e-04, -2.81375256e-05,
             2.55342205e-03,  3.75341489e-04],
           [-2.08654412e-03, -3.25114585e-04, -3.29315151e-05,
             2.37198991e-03,  3.69033499e-04],
           [-2.10267810e-03, -3.36952755e-04, -3.89255964e-05,
             1.99937598e-03,  3.72523657e-04],
           [-2.10922206e-03, -3.63722397e-04, -4.37759496e-05,
             1.89579182e-03,  3.79030801e-04],
           [-2.12145129e-03, -3.88288731e-04, -4.86582462e-05,
             1.59046775e-03,  3.71211734e-04],
           [-2.13301626e-03, -3.96416129e-04, -5.21594420e-05,
             1.72139139e-03,  3.85831034e-04],
           [-2.13654950e-03, -4.11068306e-04, -5.72034907e-05,
             1.55350802e-03,  3.63737061e-04],
           [-2.13879800e-03, -4.25663197e-04, -6.00365929e-05,
             1.36222186e-03,  3.62428442e-04],
           [-2.14773435e-03, -4.42252499e-04, -6.57073241e-05,
             1.06258038e-03,  3.64892410e-04],
           [-2.14915949e-03, -4.54384961e-04, -6.98368897e-05,
             1.10395651e-03,  3.71721416e-04],
           [-2.14748500e-03, -4.68150460e-04, -7.33358907e-05,
             9.48597038e-04,  3.67301347e-04],
           [-2.15133012e-03, -4.79653839e-04, -7.57810722e-05,
             7.16785373e-04,  3.50966149e-04],
           [-2.15146064e-03, -4.85740442e-04, -8.00086716e-05,
             4.42330871e-04,  3.61496342e-04],
           [-2.14873135e-03, -5.00603752e-04, -8.49666231e-05,
             3.71111267e-04,  3.47292101e-04],
           [-2.14428005e-03, -5.13460375e-04, -8.75992819e-05,
             1.07761285e-04,  3.53131223e-04],
           [-2.14417855e-03, -5.17663451e-04, -8.99823755e-05,
             1.76544121e-04,  3.37138446e-04],
           [-2.13687468e-03, -5.29290536e-04, -9.32575975e-05,
             5.20485825e-05,  3.29615199e-04],
           [-2.13283060e-03, -5.34335068e-04, -9.52050419e-05,
            -2.22397299e-04,  3.42362503e-04],
           [-2.12887053e-03, -5.38178865e-04, -9.69058906e-05,
            -3.67875831e-04,  3.29083226e-04],
           [-2.12348238e-03, -5.50399839e-04, -1.00712685e-04,
            -5.51171603e-04,  3.37077249e-04],
           [-2.11720938e-03, -5.51456294e-04, -1.02356073e-04,
            -7.02758561e-04,  3.20598302e-04],
           [-2.10701500e-03, -5.59527104e-04, -1.04471514e-04,
            -7.79741019e-04,  3.11578170e-04],
           [-2.09354278e-03, -5.60500310e-04, -1.06589813e-04,
            -1.04425851e-03,  3.18099420e-04],
           [-2.09099500e-03, -5.68863892e-04, -1.07735141e-04,
            -1.16661005e-03,  2.96375324e-04],
           [-2.08116164e-03, -5.69341738e-04, -1.09167365e-04,
            -1.49142973e-03,  2.88335910e-04],
           [-2.06646633e-03, -5.71680213e-04, -1.12100009e-04,
            -1.71471576e-03,  2.87211708e-04],
           [-2.05678775e-03, -5.74413129e-04, -1.12828801e-04,
            -2.00575668e-03,  3.03407422e-04],
           [-2.05056718e-03, -5.74620836e-04, -1.13277263e-04,
            -1.88282828e-03,  2.79876528e-04],
           [-2.03991350e-03, -5.81911032e-04, -1.15742866e-04,
            -2.29584031e-03,  2.68356570e-04],
           [-2.03028298e-03, -5.78198247e-04, -1.15060891e-04,
            -2.59150871e-03,  2.51064364e-04],
           [-2.01449887e-03, -5.81302783e-04, -1.17404776e-04,
            -2.48772388e-03,  2.59314760e-04],
           [-2.00566394e-03, -5.81075235e-04, -1.17456944e-04,
            -2.47704697e-03,  2.64846593e-04],
           [-1.99216455e-03, -5.81480714e-04, -1.19589651e-04,
            -2.83190967e-03,  2.33588534e-04],
           [-1.97826257e-03, -5.78750187e-04, -1.19596331e-04,
            -3.02344486e-03,  2.29829256e-04],
           [-1.96918761e-03, -5.83299406e-04, -1.19558847e-04,
            -3.17751884e-03,  2.41549818e-04],
           [-1.95496512e-03, -5.81623467e-04, -1.20678525e-04,
            -3.19645894e-03,  2.16495105e-04],
           [-1.94839952e-03, -5.80414837e-04, -1.21145379e-04,
            -3.21631019e-03,  2.12169845e-04],
           [-1.93342905e-03, -5.82545882e-04, -1.20553055e-04,
            -3.57470845e-03,  2.04961480e-04]])



#### Fit c_0 and c_1 as well as initial value of E with Nelder-Mead
no use of Jacobian and no constraints


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                #jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='Nelder-Mead',
                options={'disp':True,'maxiter':1000}) #,'eps':0.0001,'ftol':0.01}) #not NM
print(res)
```

    Optimization terminated successfully.
             Current function value: 0.000058
             Iterations: 106
             Function evaluations: 191
     final_simplex: (array([[-8.68812924e-01,  3.58400644e-01,  1.41863206e-04],
           [-8.68738149e-01,  3.58376078e-01,  1.41867431e-04],
           [-8.68791267e-01,  3.58397481e-01,  1.41863367e-04],
           [-8.68827485e-01,  3.58402559e-01,  1.41869459e-04]]), array([5.80545054e-05, 5.80545055e-05, 5.80545055e-05, 5.80545056e-05]))
               fun: 5.805450540921994e-05
           message: 'Optimization terminated successfully.'
              nfev: 191
               nit: 106
            status: 0
           success: True
                 x: array([-8.68812924e-01,  3.58400644e-01,  1.41863206e-04])



```python
# model with fitted parameters 
print(params)
print(x0)
params_fit = params.copy()
#params_fit['c_0'] = res.x[0]
#params_fit['c_1'] = res.x[1]

SCEIR_model.params = params_fit
print(SCEIR_model.params)
#ode_fit = common_models.SEIR({'beta':res.x[0], 'gamma':res.x[1],'alpha':res.x[2]})
#x0_fit = [1-1.27e-6, 1.27e-6, 0]
x0_fit = x0.copy()
#x0_fit[2] = res.x[2]
#t_fit = numpy.linspace(0, 150, 1000)
t_fit = t
SCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
# sol_fit =SCEIR_model.integrate(t_fit[0::])
sol_fit = scipy.integrate.odeint(SCEIR_model.ode, x0_fit, t_fit[1::])
#
plt.figure(figsize=(15,10))
plt.plot(t,ynoise[:,0], 'go') # infected observations
plt.plot(t,ynoise[:,1], 'go') # infected observations
plt.plot(t,ynoise[:,2], 'go') # infected observations
plt.plot(t,ynoise[:,3], 'bo') # recoverd
plt.plot(t,ynoise[:,4], 'ro') # died
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.004])
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

## Fit SCCEIR parameters to jhu data based on square_loss

### Fit c_0 and c_1 only


```python
# Initial guess of parameters, and bounding constraints
cautionparams = list(params.values())[-4:-3]
theta = [value for value in cautionparams]
theta = [0.3,0.08,2500.]
# boxBounds = [(0.2,0.8),(0.05,0.15),(100.,10000.)]
objSCCEIR = SquareLoss(theta=theta, ode=SCCEIR_model, x0=x0_c, t0=t[0], t=t[1::], y=y_jhu[1::,:],
                      state_weight=[1.,1.],state_name=['R','D'],
                      target_param=['c_0','c_1','c_2'])
# perform optimization
res = minimize(fun=objSCCEIR.cost,
                #jac=objSCCEIR.sensitivity,
                x0=theta,
                #bounds=boxBounds,
                method='Nelder-Mead',
                options={'disp':True,'maxiter':1000,'eps':0.00001})# ,'ftol':0.01})  #not BFGS
print(res)
```

    /usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:15: OptimizeWarning: Unknown solver options: eps
      from ipykernel import kernelapp as app


    Optimization terminated successfully.
             Current function value: 0.000012
             Iterations: 544
             Function evaluations: 931
     final_simplex: (array([[-0.40574359, -0.14005259,  0.05948037],
           [-0.40574432, -0.14006047,  0.05943691],
           [-0.40574314, -0.14003868,  0.05954975],
           [-0.40574384, -0.14004428,  0.05952122]]), array([1.15955941e-05, 1.15955941e-05, 1.15955941e-05, 1.15955941e-05]))
               fun: 1.1595594101394368e-05
           message: 'Optimization terminated successfully.'
              nfev: 931
               nit: 544
            status: 0
           success: True
                 x: array([-0.40574359, -0.14005259,  0.05948037])


### Fit c_0 and c_1 as well as initial value of E

#### Fit c_0 and c_1 as well as initial value of E with 'SLSQP'
does not work well
note use of special methods IV for initial value fitting of target_state


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-3]
theta = [value for value in cautionparams]
theta = [0.21,0.08,2500.]
objSCCEIR = SquareLoss(theta=theta, ode=SCCEIR_model, x0=x0_c, t0=t[0], t=t[1::], y=y_jhu[1::,:],
                      state_weight=[1.,1.],state_name=['R','D'],
                      target_param=['c_0','c_1','c_2'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCCEIR.costIV,
                jac=objSCCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='BFGS',
                options={'disp':True,'maxiter':1000,'eps':0.01,'gtol':0.01})
print(res)
```

#### Fit c_0 and c_1 as well as initial value of E with BFGS
works well:  no constraints and gtol not ftol


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='BFGS',
                options={'disp':True,'maxiter':1000,'eps':0.01,'gtol':0.01})
print(res)
```


```python
# model with fitted parameters 
print(params)
print(x0)
params_fit = params.copy()
#params_fit['c_0'] = res.x[0]
#params_fit['c_1'] = res.x[1]

SCEIR_model.params = params_fit
print(SCEIR_model.params)
#ode_fit = common_models.SEIR({'beta':res.x[0], 'gamma':res.x[1],'alpha':res.x[2]})
#x0_fit = [1-1.27e-6, 1.27e-6, 0]
x0_fit = x0.copy()
#x0_fit[2] = res.x[2]
#t_fit = numpy.linspace(0, 150, 1000)
t_fit = t
SCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
# sol_fit =SCEIR_model.integrate(t_fit[0::])
sol_fit = scipy.integrate.odeint(SCEIR_model.ode, x0_fit, t_fit[1::])
#
plt.figure(figsize=(15,10))
plt.plot(t,ynoise[:,0], 'go') # infected observations
plt.plot(t,ynoise[:,1], 'go') # infected observations
plt.plot(t,ynoise[:,2], 'go') # infected observations
plt.plot(t,ynoise[:,3], 'bo') # recoverd
plt.plot(t,ynoise[:,4], 'ro') # died
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.004])
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

#### Fit c_0 and c_1 as well as initial value of E using L-BFGS-B
this method does not work well


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                bounds=boxBounds+[(0.0000001,0.001)],
                method='L-BFGS-B',
                options={'disp':True,'maxiter':1000,'eps':0.0001,'ftol':0.001})
print(res)
```


```python
objSCEIR.residual()
```

#### Fit c_0 and c_1 as well as initial value of E with Nelder-Mead
no use of Jacobian and no constraints


```python
# Initial guess of parameters and initial condition, and bounding constraints
cautionparams = list(params.values())[-4:-2]
theta = [0.25,0.08]
boxBounds = [(0.2,0.4),(0.05,0.15)]
objSCEIR = SquareLoss(theta=theta, ode=SCEIR_model, x0=x0, t0=t[0], t=t[1::], y=ynoise[1::,:],
                      state_weight=[1.,1.,1.,1.,1.],state_name=['I_1','I_2','I_3','R','D'],
                      target_param=['c_0','c_1'],target_state=['E'])
# perform optimization
res = minimize(fun=objSCEIR.costIV,
                #jac=objSCEIR.sensitivityIV,
                x0=theta+[0.00005],
                #bounds=boxBounds+[(0.0000001,0.001)],
                method='Nelder-Mead',
                options={'disp':True,'maxiter':1000}) #,'eps':0.0001,'ftol':0.01}) #not NM
print(res)
```


```python
# model with fitted parameters 
print(params)
print(x0_c)
params_fit = params.copy()
params_fit['c_0'] = res.x[0]
params_fit['c_1'] = res.x[1]

SCCEIR_model.params = params_fit
print(SCCEIR_model.params)
x0_fit = x0_c.copy()
#x0_fit[2] = res.x[2]
#t_fit = numpy.linspace(0, 150, 1000)
t_fit = t
SCCEIR_model.initial_values = (x0_fit, t_fit[0])
# %timeit sol_fit =SCCEIR_model.integrate(t_fit[1::]) # use magic %timeit to time
# sol_fit =SCCEIR_model.integrate(t_fit[0::])
sol_fit = scipy.integrate.odeint(SCCEIR_model.ode, x0_fit, t_fit[1::])
#
plt.figure(figsize=(15,10))
plt.plot(t,y_jhu[:,0], 'bo') # recoverd
plt.plot(t,y_jhu[:,1], 'ro') # died
plt.plot(t_fit[1::], sol_fit)
plt.ylim([0,0.004])
#plt.show(())
#ode_fit.plot()

peak_i = np.argmax(sol_fit[:,2])
print('Peak infection (days)', t_fit[peak_i])
```

### Information on method options


```python
scipy.optimize.show_options(solver='minimize', method='SLSQP', disp=True)
print(' ')
scipy.optimize.show_options(solver='minimize', method='L-BFGS-B', disp=True)
```

## Plot using full control


```python
def plotmodel(solns,t,scale='linear',species='no_susc',plottitle= '',label='',
              newplot = True,models=['SEIR','SCEIR','SCCEIR']):
    """
    plot solns over 
    times t interpreted as models indicated in models parameter
    scale: alternative 'linear' or 'log'
    species alternatives 'all', 'confirmed', 'deaths', 'daily confirmed', 'daily deaths'
    plottitle : title for plot
    label : label for curve when called as part of multicurve plot
    newplot : whether to open new plot True/False
    models : list of models to include, default all three of those possible
    """
   
    nmodels = len(models)
    if len(solns) != len(models):
        print("Error: number of models must match number of solutions")
        return None
    nm = 0
    
    if newplot == True:
        plt.figure(figsize=(nmodels*8,6))
 
    for nm in range(nmodels):
        soln = solns[nm]
        if models[nm] == 'SEIR':   #SEIR
            plt.subplot(1,nmodels,nm+1)
            if scale == 'log': #Plot on log scale
                plt.semilogy()
                plt.ylim([1,10000])
            elif species != 'daily confirmed': # Plot on normal linear scale
                #plt.ylim([0,10000])
                pass
            if species == 'no_susc':
                plt.plot(tvec,soln[:,1:5],label=label)
                plt.legend(("E","I1","I2","I3"))
            elif species == 'confirmed' or species == 'daily confirmed':
                suma = np.sum(soln[:,2:7],axis=1)
                # print('length=',len(suma))
                if species == 'daily confirmed':
                    sumd = np.zeros(len(suma))
                    for i in range(1,len(suma)):
                        sumd[i] = suma[i]-suma[i-1]
                    #plt.ylim([0,1000])
                    plt.plot(tvec,sumd,label=label)
                else:
                    #plt.ylim([0,200000])
                     plt.plot(t,suma,label=label)                   
            elif species == 'all':
                plt.plot(tvec,soln,label=label)
                plt.legend(("S","E","I1","I2","I3","R","D"))
            plt.xlabel("Time (days)")
            plt.ylabel("Portion of population N")
            plt.title('SEIR %s' % plottitle)

        elif models[nm] == 'SCEIR':   #SCEIR
            #Plot
            plt.subplot(1,nmodels,nm+1)
            if scale == 'log': #Plot on log scale
                plt.semilogy()
                plt.ylim([1,10000])
            elif species != 'daily confirmed': # Plot on normal linear scale
                #plt.ylim([0,10000])
                pass
            if species == 'no_susc':
                plt.plot(t,soln[:,1:5],label=label)
                plt.legend(("E","I1","I2","I3"))
            elif species == 'confirmed' or species == 'daily confirmed':
                suma = np.sum(soln[:,2:7],axis=1)
                # print('length=',len(suma))
                if species == 'daily confirmed':
                    sumd = np.zeros(len(suma))
                    for i in range(1,len(suma)):
                        sumd[i] = suma[i]-suma[i-1]
                    #plt.ylim([0,1000])
                    plt.plot(t,sumd,label=label)
                else:
                    #plt.ylim([0,200000])
                    plt.plot(t,suma,label=label)
            elif species == 'all':
                plt.plot(t,soln,label=label)
                plt.legend(("S","E","I1","I2","I3","R","D","Sc"))
            plt.xlabel("Time (days)")
            plt.ylabel("Portion of population N")
            plt.title('SCEIR %s' % plottitle)
        elif models[nm] == 'SCCEIR':   #SCCEIR
            plt.subplot(1,nmodels,nm+1)
            if scale == 'log': #Plot on log scale
                plt.semilogy()
                plt.ylim([1,10000])
            elif species != 'daily confirmed': # Plot on normal linear scale
                #plt.ylim([0,10000])
                pass
            if species == 'no_susc':
                plt.plot(t,sol[:,1:5])
                plt.legend(("E","I1","I2","I3"))
            elif species == 'confirmed' or species == 'daily confirmed':
                suma = np.sum(soln[:,2:7],axis=1) + soln[:,9]
                if species == 'daily confirmed':
                    sumd = np.zeros(len(suma))
                    for i in range(1,len(suma)):
                        sumd[i] = suma[i]-suma[i-1]
                    # plt.ylim([0,1000])
                    plt.plot(t,sumd,label=label)
                else:
                    # plt.ylim([0,200000])
                    plt.plot(t,suma,label=label)
            elif species == 'all':
                plt.plot(t,soln,label=label)
                plt.legend(("S","E","I1","I2","I3","R","D","Sc","Ec","I1c"))
            plt.xlabel("Time (days)")
            plt.ylabel("Portion of population N")
            plt.title('SCCEIR %s' % plottitle)
    return True
```


```python
plotmodel([sol_fit],t_fit[1:],scale='linear',species='no_susc',plottitle= 'test',label='',
              newplot = True,models=['SCEIR'])
```

# Not yet implemented


```python
# R0 from arXiv paper
# does not work yet properly ! reason unclear
from pygom.model.epi_analysis import R0
ode1 = SCEIR_model.ode
ode2 = SimulateOde(state, param_list, transition=transition)
# print(ode2)
# ode = ode2.get_unrolled_obj() # Exception: Object was not initialized using a set of ode
# so need to define model via odes for this to work?
# R0( ode2, ['I_1','I_2','I_3'])
```


```python
# function for altering parameters
print(param_list) 
model = DeterministicOde(state, param_list, transition=transition)
def parameterize_model(t=np.linspace(0,365,366), c_0=0.3, c_1=0.07,
                       c_2=0.004,I1_init=0.00001, default_params = param_list, model=model):
    ic=[1-I1_init, I1_init, 0, 0, 0, 0, 0, 0]
    cautionparams = {'c_0': c_0,
              'c_1': c_1, 
              'c_2': c_2,
              'N': sum(ic)}
    
    model.initial_values = (ic, t[0])
    model.parameters = params
    sol = model.integrate(t[1::])
    model.plot()
    peak_i = np.argmax(sol[:,2])
    print('Peak infection (days)', t[peak_i] )
    
```

In this simple framework, reducing $\beta$ results in a smaller epidemic:
- the peak infection time is delayed
- the magnitude of peak infection is reduced.

Reducing beta may crudely represent giving out anti-virals, which make a person less infectious.


```python
parameterize_model(beta=1.2, t=np.linspace(0,500,5001))
```

Vaccinating 5% of the population (assuming instantaneous rollout) or natural immunity, delays the peak period, and reduces its magnitude.


```python
parameterize_model(ic=[9490,5, 5, 500], beta=0.5, gamma=0.3, t=np.linspace(0,150,10))
```

# gompy test example for FitzHugh model


```python
from unittest import main, TestCase

import numpy as np
from scipy.optimize import minimize

from pygom import SquareLoss
from pygom.model import common_models


class TestFHEstimate(TestCase):

    def setUp(self):
        # initial values
        x0 = [-1.0, 1.0]
        # params
        param_eval = [('a', 0.2), ('b', 0.2), ('c', 3.0)]
        self.target = np.array([0.2, 0.2, 3.0])
        # the time points for our observations
        t = np.linspace(0, 20, 30).astype('float64')
        ode = common_models.FitzHugh(param_eval)
        ode.initial_values = (x0, t[0])
        solution = ode.integrate(t[1::])
        self.theta = np.array([0.5, 0.5, 0.5])

        self.obj = SquareLoss(self.theta, ode, x0, t[0],
                              t[1::], solution[1::, :], ['V', 'R'])

        g = self.obj.gradient()
        assert np.linalg.norm(g) > 0

        EPSILON = np.sqrt(np.finfo(np.float).eps)

        self.box_bounds = [(EPSILON, 5.0)]*len(self.theta)

    def test_FH_sensitivity(self):
        res = minimize(fun=self.obj.cost,
                       jac=self.obj.sensitivity,
                       x0=self.theta,
                       bounds=self.box_bounds,
                       method='L-BFGS-B')

        self.assertTrue(np.allclose(self.target, res['x'], 1e-2, 1e-2))

    def test_FH_adjoint(self):
        res = minimize(fun=self.obj.cost,
                       jac=self.obj.adjoint,
                       x0=self.theta,
                       bounds=self.box_bounds,
                       method='L-BFGS-B')

        self.assertTrue(np.allclose(self.target, res['x'], 1e-2, 1e-2))

    def test_FH_IV(self):
        box_bounds = self.box_bounds + [(None, None)]*2

        res = minimize(fun=self.obj.costIV,
                       jac=self.obj.sensitivityIV,
                       x0=self.theta.tolist() + [-0.5, 0.5],
                       bounds=box_bounds,
                       method='L-BFGS-B')

        target = np.array([0.2, 0.2, 3.0, -1.0, 1.0])
        self.assertTrue(np.allclose(res['x'], target, 1e-2, 1e-2))
```


```python
test=TestFHEstimate()
```


```python
test.setUp()
```


```python
test.test_FH_IV()
```


```python

```

a198186cb19253dcb6e6448f956404fd918a4c3b
ipynb
Jupyter Notebook
Notebooks/covid-19-caution/SCEIR_pygom_fitting_N.ipynb
ProtoLife/covid-recovery
b3989f110c6961cc51da673fc33c6384cf7d6de6
MIT
Notebooks/covid-19-caution/SCEIR_pygom_fitting_N.ipynb
ProtoLife/covid-recovery
b3989f110c6961cc51da673fc33c6384cf7d6de6
MIT
2020-05-24T00:06:54.000Z
2020-05-25T19:59:44.000Z
Notebooks/covid-19-caution/SCEIR_pygom_fitting_N.ipynb
ProtoLife/covid-recovery
b3989f110c6961cc51da673fc33c6384cf7d6de6
MIT
Qwen/Qwen-72B
1. YES
2. YES