# Source code for aftercovid.models.covid_sird_cst

# coding: utf-8
"""
Implementation of a model for epidemics propagation.
"""
import numpy.random
from sklearn.linear_model import ElasticNet
from ._base_sir import BaseSIR
from .covid_sird import CovidSIRD

[docs]class CovidSIRDc(BaseSIR):
"""
Inspiration Modelling some COVID-19 data
<http://webpopix.org/covidix19.html>_.
This model considers that observed data are not the
true ones.

.. math::

\\begin{array}{rcl}
S &=& S_{obs} (1 + a)\\\\
I &=& I_{obs} (1 + b)\\\\
R &=& R_{obs} (1 + c)\\\\
D &=& D_{obs}
\\end{array}

Where :math:S, :math:I, :math:R are
:math:R_{obs} are observed.
hidden, only :math:S_{obs}, :math:I_{obs},
The following equality should be verified:
:math:S + I + R + D = N = S_{obs} + I_{obs} + R_{obs} + D_{obs}.
We also get from the previous equations:

.. math::

\\begin{array}{rcl}
dS &=& dS_{obs} (1 + a) = - \\beta \\frac{IS}{N} =
- \\beta \\frac{I_{obs}S_{obs}}{N}(1+a)(1+b) \\\\
\\Longrightarrow dS_{obs} &=& - \\beta \\frac{I_{obs}S_{obs}}{N}(1+b)
\\end{array}

And also:

.. math::

\\begin{array}{rcl}
dD &=& dD_{obs} = \\nu I = \\nu I_{obs} (1+b)
\\end{array}

And as well:

.. math::

\\begin{array}{rcl}
dR &=& dR_{obs} (1 + c) = \\mu I = \\mu (1 + b) I_{obs}  \\\\
\\Longrightarrow dR_{obs} &=& - \\nu I_{obs} \\frac{1+b}{1+c}
\\end{array}

And finally:

.. math::

\\begin{array}{rcl}
dI &=& dI_{obs} (1 + b) = -dR - dS - dD =
- \\mu \\frac{1 + b}{1+ c} I_{obs} - \\nu (1+b) I_{obs} -
- \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a)(1 + b)
\\\\
\\Longrightarrow dI_{obs} &=& - \\nu I_{obs} - \\mu I_{obs}
- \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a)
\\end{array}

This model should still verify:

.. math::

\\begin{array}{rcl}
S_{obs} + I_{obs} + R_{obs} + D_{obs} &=& N = S + I + R + D \\\\
&=& S_{obs}(1+a) + I_{obs}(1+b) + R_{obs}(1+c) + D_{obs}
\\end{array}

That gives :math:a S_{obs} + b I_{obs} + c R_{obs} = 0.

.. runpython::
:showcode:
:rst:

from aftercovid.models import CovidSIRDc

model = CovidSIRDc()
print(model.to_rst())

.. exref::
:title: SIRDc simulation and plotting

.. plot::

from pandas import DataFrame
import matplotlib.pyplot as plt
from aftercovid.models import CovidSIRDc

model = CovidSIRDc()
sims = list(model.iterate(60))
df = DataFrame(sims)
ax = df.plot(y=['S', 'I', 'R', 'D'], kind='line')
ax.set_xlabel("jours")
ax.set_ylabel("population")
r0 = model.R0()
ax.set_title("Simulation SIRDc\\nR0=%f" % r0)

plt.show()

Visual representation:

.. gdot::
:script:

from aftercovid.models import CovidSIRDc
model = CovidSIRDc()
print(model.to_dot())

See :ref:l-base-model-sir to get the methods
common to SIRx models. This model is not really working
better than :class:CovidSIRD <aftercovid.covid_sird.CovidSIRD>.
"""

P0 = [
('beta', 0.5, 'taux de transmission dans la population'),
('mu', 1 / 14., '1/. : durée moyenne jusque la guérison'),
('nu', 1 / 21., '1/. : durée moyenne jusqu\'au décès'),
('a', -1.5025135094805093e-08,
'paramètre gérant les informations cachées (S)'),
('b', 1e-5, 'paramètre gérant les informations cachées (I)'),
('c', 1e-5, 'paramètre gérant les informations cachées (R)'),
]

Q0 = [
('S', 9990., 'personnes non contaminées'),
('I', 10., 'nombre de personnes malades ou contaminantes'),
('R', 0., 'personnes guéries (recovered)'),
('D', 0., 'personnes décédées'),
]

C0 = [
('N', 10000., 'population'),
]

# c = b (nu + mu) / (mu - b * nu)
# (1 + b) / (1 + c) = 1 - nu * b / mu
eq = {
'S': '-beta * (1 + b) * I * S / N',
'I': ('beta * (1 + a) * I * S / N'
'- nu * I - mu * I'),
'R': 'mu * (1 + b) * I / (1 + c)',
'D': 'nu * (1 + b) * I'}

def __init__(self):
BaseSIR.__init__(
self,
p=CovidSIRDc.P0.copy(),
q=CovidSIRDc.Q0.copy(),
c=CovidSIRDc.C0.copy(),
eq=CovidSIRDc.eq.copy())

[docs]    def R0(self, t=0):
'''Returns R0 coefficient.
See :meth:CovidSIRD.R0 <aftercovid.models.CovidSIRD.R0>'''
return self['beta'] / (self['nu'] + self['mu'])

[docs]    def correctness(self, X=None):
'''
Returns :math:a S_{obs} + b I_{obs} + c R_{obs} = 0.

:param X: None to use inner quantities
:return: a number
'''
if X is None:
X = self.vect().reshape((1, -1))
return (X[:, 0] * self['a'] + X[:, 1] * self['b'] +
X[:, 2] * self['c']) / self['N']

[docs]    def update_abc(self, X=None, update=True, alpha=1.0, l1_ratio=0.5):
'''
Updates coefficients *a*, *b*, *c* so that method
:meth:correctness <aftercovid.models.CovidSIRDc.correctness>
returns 0. It uses ElasticNet
<https://scikit-learn.org/stable/modules/generated/
sklearn.linear_model.ElasticNet.html>_.

:param X: None to use inner quantities
:param update: True to update to the coefficients
or False to just return the results
:param alpha: see ElasticNet
:param l1_ratio: see ElasticNet
:return: dictionary
'''
if X is None:
X = self.vect().reshape((1, -1))
X = X / self['N']
cst = - self.correctness(X)
model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, fit_intercept=False)
model.fit(X, cst)
coef = model.coef_
res = dict(a=self['a'] + coef[0], b=self['b'] + coef[1],
c=self['c'] + coef[2])
if update:
self.update(**res)
return res

[docs]    def rnd(self):
'''
Draws random parameters.
Not perfect.
'''
self['beta'] = numpy.random.randn(1) * 0.1 + 0.5
self['mu'] = numpy.random.randn(1) * 0.1 + 1. / 14
self['nu'] = numpy.random.randn(1) * 0.1 + 1. / 21
self['a'] = numpy.random.rand() * 1e-5
self['b'] = numpy.random.rand() * 1e-5
self['c'] = numpy.random.rand() * 1e-5

[docs]    @staticmethod