Le modèle SEIRD de la Suisse

In [2]:
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
import pandas as pd
import time
In [3]:
class Data():
    def __init__(self):
        self.data = None
        self.dataByCountry = {}
        self.lengthDataByCountry = {}

    def extractDataCountry(self, countryCode):
        """
        Exemple de codes :  SWE (Suède), FRA (France), CHE (Suisse)
        Retourne :          dataframe Pandas pour le pays demandé
        Noms colonnes :     dateRep, day, cases, deaths, countryterritoryCode
        Méthode publique
        """
        if self.data is None:
            self._getDataFromECDC()

        self.dataByCountry[countryCode]       = self.data[self.data["countryterritoryCode"] == countryCode]
        self.lengthDataByCountry[countryCode] = len(self.dataByCountry[countryCode])

    def getDataFrameForCountry(self, countryCode):
        self.extractDataCountry(countryCode)
        return self.dataByCountry[countryCode]

    def getDataFrameLengthForCountry(self, countryCode):
        # TODO rendre ceci utile
        return self.lengthDataByCountry[countryCode]

    def _getDataFromECDC(self):
        # Chercher les données à jour
        self.data = pd.read_excel("https://www.ecdc.europa.eu/sites/default/files/documents/COVID-19-geographic-disbtribution-worldwide-{}.xlsx".format(time.strftime("%Y-%m-%d")))

        #Local pour les tests
        #self.data = pd.read_excel('./COVID-19-geographic-disbtribution-worldwide-2020-05-04.xlsx')
        #self.data = pd.read_excel('./COVID-19-geographic-disbtribution-worldwide-2020-05-25.xlsx')
        self.data = self.data[["dateRep", "day", "cases", "deaths", "countryterritoryCode"]]
        self.data = self.data.sort_values(by='dateRep', ascending=1)   # inversion de l'ordre des données (défaut : ordre anti-chronologique)
In [4]:
class Seird():

    def __init__(self, t, country, covid, i=1, exp=0, recov=0, deaths=0, seird=None):
        # Pays
        self.population = country['pop']
        self.beta       = country['beta']
        # Covid
        self.sigma      = 1/covid['dureeIncubationJour']
        self.gamma      = 1/covid['dureeContagieuxJour']
        self.mu         = country['tauxLetalite']/(100*covid['dureeContagieuxJour'])
        # modele
        self.t          = t
        self.seird      = seird
        self.inf0       = i
        self.exp0       = exp
        self.recov0     = recov
        self.death0     = deaths
        self.suscept0   = self.population - self.exp0 - self.inf0 - self.recov0 - self.death0
        self.y0         = [self.suscept0, self.exp0, self.inf0, self.recov0, self.death0]

    def computeSEIRD(self):
        """
        Retourne un dataframe Pandas  [s, e, i, r, d]
        """

        def deriv2(y, t):
            S, E, I, R, D = y
            dSdt = -(self.beta/self.population)*I*S
            dEdt = (self.beta/self.population)*I*S - self.sigma*E
            dIdt = self.sigma*E - self.gamma*I - self.mu*I
            dRdt = self.gamma*I
            dDdt = self.mu*I
            return dSdt, dEdt, dIdt, dRdt, dDdt

        resol      = odeint(deriv2, self.y0, self.t)
        self.seird = pd.DataFrame(resol, columns=('s','e','i','r','d'))


    def getSEIRD(self, para = None):
        if para is None:
            if self.seird is None:
                self.computeSEIRD()
            return self.seird
        elif para == 's':
            return self.seird['s']
        elif para == 'e':
            return self.seird['e']
        elif para == 'i':
            return self.seird['i']
        elif para == 'r':
            return self.seird['r']
        elif para == 'd':
            return self.seird['d']

    def getGamma(self):
        return self.gamma
    def getSigma(self):
        return self.sigma
    def getMu(self):
        return self.mu
    def getdDdt(self):
        return self.seird['i']*self.mu
    def getY0(self):
        return self.y0
    def getPopulation(self):
        return self.population
In [5]:
def calclMoyMobile(T, k=2):                            #Entier k : Nombre de valeurs prises
    l = len(T)                               #de chaque côté des éléments de T
    L = np.zeros(l)

    if ((2*k+1)>l):
        pass

    else :
        for i in range(l):
            if ((i-k)<0):                           #Soit k < i
                m = i
                L[i] = np.mean( T[(i-m):(i+m+1)] )  #i-m = 0

            elif ((i+k)>(l-1)):                     #Soit k > l-1-i
                m = l-1-i
                L[i] = np.mean( T[(i-m):(i+m+1)] )  #i+m+1 = l

            else :
                L[i] = np.mean( T[(i-k):(i+k+1)] )

    return L
In [6]:
def graphe(df, seird, nomVariable, xNom, yNom, titreGraphe, moyMob=2):

    def _plot(data, couleur, labelName, typePlot="ligneD"):
        if typePlot == "bar":
            plt.bar(t, data, color=couleur, label=labelName)
        elif typePlot == "ligneD":
            plt.plot(t, data, couleur, alpha=1, lw=2, label=labelName)

    n = len(df)
    t = np.linspace(0, n, n)

    if nomVariable == 'morts':
        _plot(df['deaths'],    '#a6cee3', 'Nombre réel de décès', 'bar')
        _plot(calclMoyMobile(df['deaths'], moyMob), '#1f78b4', 'MoyenneMobile ({} jours)'.format(moyMob))
        #if moyMob is not None:
        _plot(seird.getdDdt(), '#ba4f83', 'Modèle')

    elif nomVariable == 'cas':
        # TODO ne fonctionne pas, le modèle donne un nombre de cas bien trop grand.
        _plot(df['cases'] , 'r', 'Nombre réel de cas', 'bar')
        _plot(seird.getSEIRD('i'), 'g', 'Modèle')


    plt.xlabel(xNom)
    plt.ylabel(yNom)
    plt.title(titreGraphe)
    plt.grid(b=True, which='major', c='grey', lw=1, ls='-', alpha=0.2)
    legend = plt.legend()
    legend.get_frame().set_alpha(0.5)

    plt.show()
In [8]:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
Le code python est caché par défault.
Pour l'activer ou le désactiver cliquez <a href="javascript:code_toggle()">ici</a>.''')
Out[8]:
Le code python est caché par défault. Pour l'activer ou le désactiver cliquez ici.

Suisse

Mesures prises par le gouvernement

  • 28 février : interdiction des rassemblements de plus de 1000 personnes
  • 13 mars : fermeture des écoles, contrôle des frontières, et interdiction des rassemblements de plus de 100 personnes.
  • 16 mars : le gouvernement parle de "situation extraordinaire". Fermeture des commerces non essentiels, déploiement de l'armée en soutient au système médical, fermeture partielle des frontières.
  • A partir du 20 mars, interdiction des rassemblements de plus de 5 personnes. "Pas de confinement" d'après le gouvernement.

Le retour à la normale est réalisé en trois étapes :

  • 27 avril : Autorisation des métiers avec contacts en nombre restreint (dentistes, coiffeurs, etc.)
  • 11 mai : Commerces, écoles primaires ouvertes.
  • 8 juin : Musées, zoos, bibliothèques, etc. ouverts
  • La "situation extraordinaire" sera levée le 19 juin.
  • Les évènements de plus de 1000 personnes restent interdits jusqu'à la fin août.

Modélisation

Modèle

Le modèle utilisé est un modèle SEIRD simple sans phases. On obtient une bonne approximation de la moyenne mobile (sur 3 jours) des décès quotidiens avec un beta = 0.52 et un taux de létalité de seulement 0.02%. Cette valeur est très surprenante étant donné que les estimations sont plutôt de l'ordre de [0,4, 0,9]%.

Remarques

On remarque d'une part que les mesures prises par le gouvernement semblent avoir cassé la dynamique de la courbe des décès. On observe un plateau dès le 90ième jour environ. Une seconde remarque est qu'il y a eu un rebond dans les environs du 135-ième jour. Ceci coïncide plus ou moins à la seconde pĥase des diminutions des restrictions (réouverture des écoles et des commerces).

In [7]:
def main():
    ##############################################################################################################
    #           Il suffit de modifier LISTE_PAYS, COVID, et COUNTRY_CODE pour adapter le modèle.                 #
    ##############################################################################################################
    #                                                                                                            #
    #                                                                                                            #
    #  LISTE_PAYS permet d'éviter de créer plusieurs fichiers pour plusieurs pays ou de se trimbaler             #
    #       une tonne de variables. C'est un dictionnaire dont les valeurs sont des dictionnaires.               #
    #       On accède au dictionnaire de CHE par LISTE_PAYS['CHE'], et on accède à la population                 #
    #       la Suisse en tapant LISTE_PAYS['CHE']['pop'].                                                        #
    #       Attention à la syntax si on veut rajouter des pays ou en enlever : virgules, doubles points...       #
    #                                                                                                            #
    #  COVID est aussi un dictionnaire.                                                                          #
    #                                                                                                            #
    #  COUNTRY_CODE est le pays d'étude. Il faut le changer pour passer de la Suisse à la France ou la Suède.    #
    #
    #  Les données sont une dataframe Pandas avec pour variables dateRep, day, cases, deaths,
    #    countryterritoryCode
    #                                                                                                            #
    ##############################################################################################################

    global seird, Tn, dDdtMoyen, dDdtTheo


    LISTE_PAYS = {"CHE":
                    {'pop':8.6*10**6,
                     'beta': 0.52,
                     'tauxLetalite':0.02},

                  "FRA":{   # Je n'ai pas modelise la France.
                      'pop':6.7*10**7,
                      'beta':1/3.00,
                      'tauxLetalite': 0.8},

                  "SWE":{
                      'pop':1*10**6,
                      'beta':1/3.10,
                      'tauxLetalite': 0.47}}

    COVID = {'dureeIncubationJour': 5.5,  'dureeContagieuxJour': 10}       # Suisse
    

    COUNTRY_CODE = "CHE"   # le pays que l'on choisi d'etudier



    data = Data()
    donneePays = data.getDataFrameForCountry(COUNTRY_CODE)
    n = len(donneePays)
    t = np.linspace(0, n, n)
    division = 7


    seird = Seird(t, LISTE_PAYS[COUNTRY_CODE], COVID)
    seird.computeSEIRD()
    dDdt = seird.getdDdt()



    #Suisse
    graphe(donneePays, seird,'morts', 'jours', 'Décès quotidiens', 'Suisse - nombre de décès')


main()