Introduction
Ce tutoriel montre comment estimer des paramètres de diffusion et de réaction sur les données images provenant du réseau Sentinelles. Les deux paramètres (p1,p2) représentent la diffusion et le taux de réaction. La période de simulation s'étale sur les semaines 2 à 5 du mois de janvier 2019. La réaction et la diffusion sont ajustés afin de minimiser l'erreur quadratique entre la simulation et les données des semaines 3 à 5.
L'intérêt de ce tutoriel est de montrer qu'il est facile d'ajuster un modèle sur des données images. Le modèle mis en oeuvre n'a que peu d'intérêt épidémiologique.
La suite de cette page décrit comment mener cette étude dans MSE, d'une part via l'interface, d'autre part via un script python
Description de l'étude
Il s'agit donc d'une étude portant sur 1 espèce avec 2 paramètres :
import mse
MY_OUTPUT_DIR=os.getenv("MSE_OUTPUT_DIR", "/home/usermse/GRIPPE/")
mse.setStringParam(mse.MSE_KEY_OUTPUT_DIR,MY_OUTPUT_DIR)
mse.setIntParam(mse.MSE_KEY_MODEL_NBSPECIES,1)
mse.setIntParam(mse.MSE_KEY_MODEL_NBPARAMS,2)
mse.initGradEdp()
mse.initModel()
Description de la zone d'étude
La zone d'étude est une partie de la France :
mse.setStringParam(mse.MSE_KEY_GEOM_DIR,MY_GEOM_DIR)
mse.setRealParam(mse.MSE_KEY_GEOM_NPERLENGTH,1e-4)
mse.setIntParam(mse.MSE_KEY_GEOM_SYM,0)
mse.setIntParam(mse.MSE_KEY_GEOM_REVEDGE,1)
mse.buildGeom()
mse.buildDirectories()
os.system("cd "+MY_OUTPUT_DIR+";/opt/freefem++/src/nw/FreeFem++-nw "+MY_OUTPUT_DIR+"/buildXY.edp")
Description du modèle de réaction diffusion
Le modèle étudié est le suivant :
Le terme de réaction peut se définir de la façon suivante :
La définition de la diffusion se fait avec le paramètre p1 :
L'état initial permet de référencer une covariable D02 :
D02 est créée à partir de l'image de la semaine 02. Le tutoriel 'import d'une image' permet de voir comment la créer.
mse.setfxExp(0,"a*(1-a)")
mse.setdyfxExp(0,0,"-2*a + 1")
mse.setdPiFj(0,0,"0")
mse.setdPiFj(1,0,"0")
mse.buildModel()
mse.buildModelParams()
mse.buildWVarf()
mse.printTSOneStepW()
mse.setDiffusionxExp(0,0,"p1")
mse.setdpiDiffusionxExp(0,0,0,"1")
mse.setdpiDiffusionxExp(1,0,0,"0")
mse.setDiffusionyExp(0,0,"p1")
mse.setdpiDiffusionyExp(0,0,0,"1")
mse.setdpiDiffusionyExp(1,0,0,"0")
mse.buildDiffusion()
mse.initInitialState()
mse.setIntParam(mse.MSE_KEY_FORCE_INITIAL_STATE,1)
mse.setInitialState(0,0,"D02")
mse.setdPICiFj(0,0,"0")
mse.setdPICiFj(1,0,"0")
mse.printInitialState()
Définition des paramètres du simulateur
La durée de simulation est de 3 semaines :
mse.setRealParam(mse.MSE_KEY_SIM_MINSTEP,0.1)
mse.setIntParam(mse.MSE_KEY_SIM_MAXMULTSTEP,1)
mse.setRealParam(mse.MSE_KEY_SIM_NEWTON_CRIT,0.01)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MIN_IT,1)
mse.setIntParam(mse.MSE_KEY_SIM_NEWTON_MAX_IT,5)
mse.setIntParam(mse.MSE_KEY_SIM_FORCE_POSITIV,1)
mse.buildTS()
Définition de la fonction objectif à minimiser
Il s'agit ici de minimiser l'intégrale de la différence quadratique entre la simulation et les données pour les pas de temps spécifiés.
La fonction de coût associé au paramètre (p1,p2) est:
Avec uref la simulation servant de référence.
Cela se spécifie de la façon suivante dans l'interface :
mse.buildContrast()
Description des données
Les données proviennent du réseau Sentinelles. Il s'agit d'images représentant des cartes d'interpolation des incidences des semaines 2, 3, 4 et 5 du mois de janvier 2019. Ces images sont importées comme décrit dans le tutoriel 'import d'images'. Les fichiers obtenus sont ensuite copiés dans le répertoire "/home/mseuser/GRIPPE/PROCESSUS/DATA" sont les noms "0s02D0Ref 0s12D0Ref 0s22D0Ref".
Les dates d'observations sont précisées dans un fichier dataDates.txt qu'il faudra créer dans le répertoire de l'étude initialement spécifié "/home/mseuser/GRIPPE/PROCESSUS/DATA".
Description du problème d'optimisation
Il s'agit d'un problème d'optimisation sous contraintes. Pour le résoudre, un algorithme de Quasi-Newton BFGS est couplé au simulateur. Il s'agit d'un algorithme itératif pour lequel le point de départ est précisé ainsi que les conditions de convergence :
mse.setBfgsLowerBound(0,0.1)
mse.setBfgsLowerBound(1,1e3)
mse.setBfgsUpperBound(0,8)
mse.setBfgsUpperBound(1,1e10)
mse.setBfgsInitialState(0,0.2)
mse.setBfgsInitialState(1,1e6)
mse.setBfgsTol(1e-3)
mse.buildBfgs()
Résultat de l'optimisation
13 appels au simulateur suffisent pour voir converger l'algorithme d'optimisation vers un minimum local (p1,p2)=(0.8,2.2e9). On peut remarquer que ces paramètres correspondent à une vitesse de propagation de 90km/semaine.
La figure suivante présente les trajectoires des itérations de l'algorithme d'optimisation pour 3 points de départ.
Image visualisant les données et la simulation optimale :
