Source code for aftercovid.models.covid_sird

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


[docs]class CovidSIRD(BaseSIR): """ Inspiration `Modelling some COVID-19 data <http://webpopix.org/covidix19.html>`_. .. runpython:: :showcode: :rst: from aftercovid.models import CovidSIRD model = CovidSIRD() print(model.to_rst()) .. exref:: :title: SIRD simulation and plotting .. plot:: from pandas import DataFrame import matplotlib.pyplot as plt from aftercovid.models import CovidSIRD model = CovidSIRD() sims = list(model.iterate(60)) df = DataFrame(sims) print(df.head()) 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 SIRD\\nR0=%f" % r0) plt.show() Visual representation: .. gdot:: :script: from aftercovid.models import CovidSIRD model = CovidSIRD() print(model.to_dot()) See :ref:`l-base-model-sir` to get the methods common to SIRx models. """ 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'), ] 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'), ] eq = { 'S': '- beta * S / N * I', 'I': 'beta * S / N * I - mu * I - nu * I', 'R': 'mu * I', 'D': 'nu * I' } def __init__(self): BaseSIR.__init__( self, p=CovidSIRD.P0.copy(), q=CovidSIRD.Q0.copy(), c=CovidSIRD.C0.copy(), eq=CovidSIRD.eq.copy())
[docs] def R0(self, t=0): ''' Returns R0 coefficient. :param t: unused This coefficients tells how many people one contagious infects in average. Let's consider the contagious population :math:`I_t`. A time *t+1*, there are :math:`\\beta I_t` new cases and :math:`(\\nu + \\mu)I_t` disappear. So at time *t+2*, there :math:`\\beta I_t(1 - \\nu - \\mu)` new cases, at time *t+d+1*, it is :math:`\\beta I_t(1 - \\nu - \\mu)^d`. We finally find that: :math:`R_0=\\beta \\sum_d (1 - \\nu - \\mu)^d` and :math:`R_0=\\frac{\\beta}{\\nu + \\mu}`. ''' return self['beta'] / (self['nu'] + self['mu'])
[docs] def correctness(self, X=None): """ Unused. """ if X is None: X = self.vect().reshape((1, -1)) return numpy.zeros(X.shape)
[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
[docs] @staticmethod def add_noise(X, epsilon=1.): """ Tries to add reasonable noise to the quantities stored in *X*. :param epsilon: amplitude :return: new X """ rnd = numpy.random.randn(*X.shape) * epsilon + 1. rnd = numpy.maximum(rnd, 0) X2 = X.copy().astype(numpy.float64) grad = (X2[1:] - X2[:-1]) grad[:, 3] *= rnd[1:, 3] grad[:, 2] *= (rnd[1:, 2] - 1) / 5 + 1. X2[1:] = X2[0, :] + numpy.cumsum(grad, axis=0) fact = numpy.multiply(numpy.sum(X2, axis=1), 1. / numpy.sum(X, axis=1)) X2 = numpy.multiply(X2, fact.reshape(X.shape[0], 1)) delta = numpy.sum(X, axis=1) - numpy.sum(X2, axis=1) delta /= X.shape[1] X2 = numpy.add(X2, delta.reshape(-1, 1)) return X2