fix bug about w1 edo

AEDESINTER_AUTOROUTES_ADAPTIV/POSTSIMU/computeBS.edp

-real BSScore=0;
+include "defUsedDep.edp";
 
-int indexYear=CURYEAR-2003;
+ BSScore=0;
 
+ indexYear=LASTYEAR-2003;
+ cmpData=0;
 for (int i=0;i<114;i++){
+  if (useDep[i]==0)
+    continue;
+  cmpData++;
   real lBS=0;
   //premier terme
   lBS= (1-vdepP[i])^2;
-  if (Tobs[indexYear*114+i] == 0)
+  if (Tobs[indexYear*114+i] == 1)
     lBS+=-1;
   BSScore+=lBS*lBS;
-  
+  cout<<lBS*lBS<<endl;
   //second terme
-  lBS= (1-vdepP[i])*vdepP[i];
-  if (Tobs[indexYear*114+i] == 1)
+  lBS= 2*(1-vdepP[i])*vdepP[i];
+  if (Tobs[indexYear*114+i] == 2)
     lBS+= -1;
   BSScore+=lBS*lBS;
-  
+  cout<<lBS*lBS<<endl;
   //troisieme terme
   lBS= (vdepP[i])^2;
-  if (Tobs[indexYear*114+i] == 2)
+  if (Tobs[indexYear*114+i] == 3)
     lBS+= -1;
   BSScore+=lBS*lBS;
-  
+  cout<<lBS*lBS<<endl;
+  cout<<"Tobs[indexYear*114+i]="<<Tobs[indexYear*114+i]<<" vdepP[i]="<<vdepP[i]<<" BSSCORE="<<BSScore<<endl<<endl;
 }
+BSScore=BSScore/(2.0*cmpData);
+cout<<"n datas="<<cmpData<<endl;

AEDESINTER_AUTOROUTES_ADAPTIV/POSTSIMU/defP.edp

 real[int] vdepP(114);
-real coefP=1;
+
 for (int i=0;i<114;i++){
   vdepP[i]=1;
   if (coefP*vdepMax[i]<1)
     vdepP[i]=coefP*vdepMax[i];
-
+  cout<<" vdepP[i]="<<vdepP[i]<<endl;
 }

AEDESINTER_AUTOROUTES_ADAPTIV/POSTSIMU/loopOnSimu.edp

@@ -5,8 +5,9 @@ include "../defGeomGen.edp"
 include "../defDofsGen.edp"
 include "../buildTobs.edp"
 
-//string MODE="DOVISU";
-string MODE="COMPUTEBS";
+string SRCDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/";
+string MODE="DOVISU";
+//string MODE="COMPUTEBS";
 
 include "defVdepMax.edp"
 include "initVdepMax.edp"
@@ -15,71 +16,11 @@ include "../postTraitementLoad.edp"
 int curStep=0;
 int numRun=0;
 int size;
-CURYEAR=2019;
+CURYEAR=2015;
 {
-  cout<<"../SAVE/u2d0_"+string(CURYEAR)+"_rep_"+string(numRun)<<endl;
-  GetSizeVec(size,"../SAVE/u2d0_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u0(size);
-  //cout<<"size="<<size<<" "<<Vh0.ndof<<endl;
-  int nbDays=size/Vh0.ndof;
+  int nbDays;
+  include "loadState.edp";
   
-  LoadVec(u0,"../SAVE/u2d0_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d1_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u1(size);
-  LoadVec(u1,"../SAVE/u2d1_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d2_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u2(size);
-  LoadVec(u2,"../SAVE/u2d2_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d3_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u3(size);
-  LoadVec(u3,"../SAVE/u2d3_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d4_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u4(size);
-  LoadVec(u4,"../SAVE/u2d4_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d5_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u5(size);
-  LoadVec(u5,"../SAVE/u2d5_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d6_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u6(size);
-  LoadVec(u6,"../SAVE/u2d6_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d7_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u7(size);
-  LoadVec(u7,"../SAVE/u2d7_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d8_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u8(size);
-  LoadVec(u8,"../SAVE/u2d8_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d9_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u9(size);
-  LoadVec(u9,"../SAVE/u2d9_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d10_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u10(size);
-  LoadVec(u10,"../SAVE/u2d10_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d11_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u11(size);
-  LoadVec(u11,"../SAVE/u2d11_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d12_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u12(size);
-  LoadVec(u12,"../SAVE/u2d12_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d13_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u13(size);
-  LoadVec(u13,"../SAVE/u2d13_"+string(CURYEAR)+"_rep_"+string(numRun));
-  
-  GetSizeVec(size,"../SAVE/u2d14_"+string(CURYEAR)+"_rep_"+string(numRun));
-  real[int] u14(size);
-  LoadVec(u14,"../SAVE/u2d14_"+string(CURYEAR)+"_rep_"+string(numRun));
   int nbCallSave=0;
   int firstStep=364-nbDays;
   for ( curStep=firstStep;curStep<364;curStep++){

AEDESINTER_AUTOROUTES_ADAPTIV/buildTobs.edp

+../DATA_AEDES/SRC/buildTobs.edp
\ No newline at end of file

AEDESINTER_AUTOROUTES_ADAPTIV/cas1.edp

@@ -12,8 +12,11 @@ include "buildTobs.edp";
 include "defDataAssimilation.edp";
 
 int numRun=0;
-int NnumRun=10;
-string SAVEDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE/WORK_VERSION/";
+int NnumRun=1;
+string SAVEDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE/";
+string LOADDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE/";
+string SRCDIR="./";
+/*
 if (ARGV.n>2){
   numRun=atoi(ARGV[2]);
   NnumRun=10;
@@ -21,10 +24,89 @@ if (ARGV.n>2){
 if (ARGV.n>3){
   NnumRun=atoi(ARGV[3]);
 }
+*/
+/*
+  u: pop adulte
+  w1: oeuf diaposant en début d'hier
+  w0: oeuf ayant passe l'hiver
+  
+  Simulation avec assimilation de 2008 a 1019
+  ---------------------------------------------
+  * 2008 --> 2009
+     T0 = 1 fev 2008
+     IC=0
+     dc u=0 jusqu'au 1 septembre
+     1 septembre assimilation des donnees 2008 avec thetaAssim=1
+     Sauvegarde ds SAVE/u2d0_2008_rep_0
+     en fin d'annee t= 31 dec : interAnnualStage.edp
+     calcul de w0 a partir de w1
+     sauvegarde de w1 dans SAVE/w1s2d0_0_2008
+     sauvegarde de w0 dans SAVE/w0s2d0_0_2008
+  * 2009 --> 2010
+     T0 = 1 fev 2008
+     u=0 w1=0 ou touche pas a w0
+     simulation Sauvegarde ds SAVE/u2d0_2009_rep_0
+     sauvegarde de w1 dans SAVE/w1s2d0_0_2009
+     sauvegarde de w0 dans SAVE/w0s2d0_0_2009
+
+   ect
+
+   Simulation sans assimilation de FIRSTYEAR a FIRSTYEAR+k
+   --------------------------------------------------------
+   
+   Il s'agit de simuler la periode du 1 fev FIRSTYEAR au 31 dec FIRSTYEAR+(k-1) sans assimilation de donnee.
+   Principe: on simule seleument la derinere annee FIRSTYEAR+(k-1) a partir des de la simu FIRSTYEAR+((k-1)-1)
+
+   FIRSTYEAR ds {2009,..,2019}
+   k ds {1,..,2020-FIRSTYEAR} (note NBPREDYEAR=k dans le code)
+   ie, couple (FIRSTYEAR,k)possible :
+     * pour FIRSTYEAR=2019
+        (2019,1)
+     * pour FIRSTYEAR=2018
+        (2018,1), (2018,2)
+     *
+     *
+     * pour FIRSTEYEAR=2009
+        (2009,1),(2009,2),...,(2009,11) 
+
+   Sauvegarde des trajectoires des etat de pop
+   
+   si k=1: ie simulation de l'annee FIRSTYEAR sans assimilation
+   ------------------------------------------------------------
+     On utilise la simulation de 2008 a 2019 avec assimilation de donnees:
+     -> Etat initial a partir des oeufs diaposant de l'annee FIRSTYEAR ayant passe l'hiver,
+         c'est a dire w0s2d5_0_$(FIRSTYEAR-1) (ici 0 est le num de repete)
+   
+   si k>1: ie simulation de l'annee FIRSTYEAR+(k-1) sans assimilation
+   ------------------------------------------------------------------
+    On utilise la simulation FIRSTYEAR+(k-1):
+    -> Etat initial a partir des oeufs diaposant de l'annee FIRSTYEAR+(k-1) ayant passe l'hiver,
+       c'est a dire w0s2d4_0_$(k-1)_$FIRSTYEAR (ici 0 est le num de repete)
+
+   pour tout k:
+   -> sauvegarde de la population d'adulte au cours de l'annee FIRSTYEAR+(k-1) dans
+        u2d5_$k_$FIRSTYEAR_rep_0 (ici 0 est le num de repete) (code: saveState.edp)
+   -> sauvegarde de la population d'oeufs diaposant ayant passe l'hiver (code: saveIntraAnnuel.edp)
+        de l'annee FIRSTYEAR+k dans w0s2d4_0_$(k)_$FIRSTYEAR (ici 0 est le num de repete)
+
+
+dans le repertoire POSTSIMU, Puis calcul des Brier Score pour tout les couples (FIRSTYEAR,k) a partir des sauvegarde u2d5_$k_$FIRSTYEAR_rep_0
+ */
+if (ARGV.n>2){
+  FIRSTYEAR=atoi(ARGV[2]);
+}
+if (ARGV.n>3){
+  NBPREDYEAR=atoi(ARGV[3]);
+  if (NBPREDYEAR<=0){
+    cout<<"usage, Freefem++ cas1.edp FIRSTYEAR>=2008 NBPREDYEAR>=0\n";
+    exit(0);
+  }
+  LASTYEAR=FIRSTYEAR+NBPREDYEAR-1;
+}
 int scenarioLow=0;
 int scenarioMedian=1;
 int scenarioHight=2;
-int scenario=scenarioLow;
+int scenario=scenarioMedian;
 int seed = 1234+numRun;
 int callSolver=1;
 randinit ( seed );
@@ -52,20 +134,23 @@ real deathrate=12;
 real mupowm1=0;
 int freqVtk=100000;
 real Tw=1.0;
-real Rstar=0.14*365.0;
+//real Rstar=0.14*365.0;
 real timeScal=365;
 string GDDTHRESHOLD="GDD_269/";
 if (scenario==scenarioLow){
-  Rstar=0.12*365;
+  //Rstar=0.12*365;
   Tw=2;
   GDDTHRESHOLD="GDD_294/";
-  SAVEDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE/LOW/";
+}
+if (scenario==scenarioMedian){
+  //Rstar=0.12*365;
+  Tw=1;
+  GDDTHRESHOLD="GDD_269/";
 }
 if (scenario==scenarioHight){
-  Rstar=0.15*365;
+  //Rstar=0.15*365;
   Tw=0;
   GDDTHRESHOLD="GDD_246/";
-  SAVEDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE/HIGHT/";
 }
 
 real K=1.0;
@@ -141,22 +226,27 @@ func real J(real[int] & u)
   //include "buildMatGenW01.edp";
       
   cout<<"end ic\n";
-  ofstream infoFileNull("infosimu.txt");
+//  ofstream infoFileNull("infosimu.txt");
   //nbCallSave=0;
   cout<<"begin loop\n";
-  for (CURYEAR=FIRSTYEAR; CURYEAR<=LASTYEAR; CURYEAR++){
+  if (NBPREDYEAR>0)
+    CURYEAR=FIRSTYEAR+NBPREDYEAR-1;
+  else
+    CURYEAR=FIRSTYEAR;
+  
+  for (CURYEAR=CURYEAR; CURYEAR<=LASTYEAR; CURYEAR++){
     nbCallSave=0;
     if (CURYEAR>2019){
       real p;
       p=randreal1();
       int rg=p/0.2+1;
-      CURYEARRAND=2011+rg;
+      CURYEARRAND=2008+rg;
       if (CURYEARRAND>2019)
         CURYEARRAND=2019;
       cout<<"CURYEAR "<<CURYEAR<<" use data "<<CURYEARRAND<<endl;
     }
-    ofstream infoFile("infosimu.txt",append);
-    if (CURYEAR==2011){
+//    ofstream infoFile("infosimu.txt",append);
+    if (CURYEAR==2008){
       //nbSteps=nbSteps2011;
       T0=1-122.0/365.0;
     }else{
@@ -175,13 +265,15 @@ func real J(real[int] & u)
     w1max=-100000;
     RRmin=100000;
     RRmax=-100000;
-    if (CURYEAR!=FIRSTYEAR){
+//    if (CURYEAR!=FIRSTYEAR){
       include "zeroAdultAndW1.edp";
-    }
-   
+//    }
+   plot(w0s2d0,w0s2d1,w0s2d2,w0s2d3,w0s2d4,w0s2d5,w0s2d6,w0s2d7,w0s2d8,w0s2d9,w0s2d10,w0s2d11,w0s2d12,w0s2d13,w0s2d14,cmm=string(CURYEAR),value=true,fill=true,wait=true);
+   if (CURYEAR>=2010)
+   plot(u2d0,u2d1,u2d2,u2d3,u2d4,u2d5,u2d6,u2d7,u2d8,u2d9,u2d10,u2d11,u2d12,u2d13,uk2d14,cmm="uk k="+"curT="+string(curT),value=true,fill=true,wait=true);
     cout<<"doing adaptiv timestepping "<<endl;//with nbSteps="<<nbSteps<<endl;
     include "timeSteppingGenAdaptiv.edp";
-    printInfoYear(infoFile);
+//    printInfoYear(infoFile);
     printInfoYear(cout);
         
     /*
@@ -221,7 +313,7 @@ func real J(real[int] & u)
 
 
 real resTllh;
-for (int ii=0; ii<NnumRun;ii++)
+//for (int ii=0; ii<NnumRun;ii++)
        resTllh=J(u6);
 //grad6=DJ(u6);
 //cout<<"PostTrait "<<grad6(0)<<" "<<grad6(1)<<" "<<grad6(2)<<" "<<grad6(3)<<" "<<grad6(4)<<endl;;

AEDESINTER_AUTOROUTES_ADAPTIV/casNoAssim.edp

+load "iovtk"
+load "resizeCSR"
+load "BinaryIO";
+include "myMacro.edp";
+include "defGeomGen.edp"
+include "defDofsGen.edp"
+
+include "defTigerVar.edp"
+include "defSimuParamsGen.edp";
+include "paramsTS.edp";
+include "buildTobs.edp";
+include "defDataAssimilation.edp";
+thetaAssim=0;
+int numRun=0;
+int NnumRun=1;
+string SAVEDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE_BUFF/";
+string LOADDIR="/home/olivierb/aedes/AEDESINTER_AUTOROUTES_ADAPTIV/SAVE_NOASSIM/";
+string SRCDIR="./";
+/*
+if (ARGV.n>2){
+  numRun=atoi(ARGV[2]);
+  NnumRun=10;
+}
+if (ARGV.n>3){
+  NnumRun=atoi(ARGV[3]);
+}
+*/
+/*
+  u: pop adulte
+  w1: oeuf diaposant en début d'hier
+  w0: oeuf ayant passe l'hiver
+  
+  Simulation avec assimilation de 2008 a 1019
+  ---------------------------------------------
+  * 2008 --> 2009
+     T0 = 1 fev 2008
+     IC=0
+     dc u=0 jusqu'au 1 septembre
+     1 septembre assimilation des donnees 2008 avec thetaAssim=1
+     Sauvegarde ds SAVE/u2d0_2008_rep_0
+     en fin d'annee t= 31 dec : interAnnualStage.edp
+     calcul de w0 a partir de w1
+     sauvegarde de w1 dans SAVE/w1s2d0_0_2008
+     sauvegarde de w0 dans SAVE/w0s2d0_0_2008
+  * 2009 --> 2010
+     T0 = 1 fev 2008
+     u=0 w1=0 ou touche pas a w0
+     simulation Sauvegarde ds SAVE/u2d0_2009_rep_0
+     sauvegarde de w1 dans SAVE/w1s2d0_0_2009
+     sauvegarde de w0 dans SAVE/w0s2d0_0_2009
+
+   ect
+
+   Simulation sans assimilation de FIRSTYEAR a FIRSTYEAR+k
+   --------------------------------------------------------
+   
+   Il s'agit de simuler la periode du 1 fev FIRSTYEAR au 31 dec FIRSTYEAR+(k-1) sans assimilation de donnee.
+   Principe: on simule seleument la derinere annee FIRSTYEAR+(k-1) a partir des de la simu FIRSTYEAR+((k-1)-1)
+
+   FIRSTYEAR ds {2009,..,2019}
+   k ds {1,..,2020-FIRSTYEAR} (note NBPREDYEAR=k dans le code)
+   ie, couple (FIRSTYEAR,k)possible :
+     * pour FIRSTYEAR=2019
+        (2019,1)
+     * pour FIRSTYEAR=2018
+        (2018,1), (2018,2)
+     *
+     *
+     * pour FIRSTEYEAR=2009
+        (2009,1),(2009,2),...,(2009,11) 
+
+   Sauvegarde des trajectoires des etat de pop
+   
+   si k=1: ie simulation de l'annee FIRSTYEAR sans assimilation
+   ------------------------------------------------------------
+     On utilise la simulation de 2008 a 2019 avec assimilation de donnees:
+     -> Etat initial a partir des oeufs diaposant de l'annee FIRSTYEAR ayant passe l'hiver,
+         c'est a dire w0s2d5_0_$(FIRSTYEAR-1) (ici 0 est le num de repete)
+   
+   si k>1: ie simulation de l'annee FIRSTYEAR+(k-1) sans assimilation
+   ------------------------------------------------------------------
+    On utilise la simulation FIRSTYEAR+(k-1):
+    -> Etat initial a partir des oeufs diaposant de l'annee FIRSTYEAR+(k-1) ayant passe l'hiver,
+       c'est a dire w0s2d4_0_$(k-1)_$FIRSTYEAR (ici 0 est le num de repete)
+
+   pour tout k:
+   -> sauvegarde de la population d'adulte au cours de l'annee FIRSTYEAR+(k-1) dans
+        u2d5_$k_$FIRSTYEAR_rep_0 (ici 0 est le num de repete) (code: saveState.edp)
+   -> sauvegarde de la population d'oeufs diaposant ayant passe l'hiver (code: saveIntraAnnuel.edp)
+        de l'annee FIRSTYEAR+k dans w0s2d4_0_$(k)_$FIRSTYEAR (ici 0 est le num de repete)
+
+
+dans le repertoire POSTSIMU, Puis calcul des Brier Score pour tout les couples (FIRSTYEAR,k) a partir des sauvegarde u2d5_$k_$FIRSTYEAR_rep_0
+ */
+if (ARGV.n>2){
+  FIRSTYEAR=atoi(ARGV[2]);
+}
+if (ARGV.n>3){
+  NBPREDYEAR=atoi(ARGV[3]);
+  if (NBPREDYEAR<=0){
+    cout<<"usage, Freefem++ cas1.edp FIRSTYEAR>=2008 NBPREDYEAR>=0\n";
+    exit(0);
+  }
+  LASTYEAR=FIRSTYEAR+NBPREDYEAR-1;
+}
+int scenarioLow=0;
+int scenarioMedian=1;
+int scenarioHight=2;
+int scenario=scenarioMedian;
+int seed = 1234+numRun;
+int callSolver=1;
+randinit ( seed );
+//include "defPluginVaraibleGen.edp"
+//real nu=0.0;//1.0;//2D to 1D
+//real mu=0.0;//0.5;//1D to 2D /ok
+//real alpha1D=0.0;
+int freqAff=1000;
+real B0=0;
+real betamodel=0;
+real alphaw0,alphaw1;
+//exp(-2*h0/365)=0.8;
+//h0=-365*ln(0.8)/2=40.7;
+alphaw0=40.7;
+alphaw1=1;
+//h0 tel que exp(-(2 jour)*h0) = 0.8
+//exp(-2*h0/364)=0.8
+//h0/364~0.1
+//h0=36.4
+//real lifeexpectancy=1.0/12.0;
+
+real h0=36.4;
+//deathrate=1/lifeexpected
+real deathrate=12;
+real mupowm1=0;
+int freqVtk=100000;
+real Tw=1.0;
+//real Rstar=0.14*365.0;
+real timeScal=365;
+string GDDTHRESHOLD="GDD_269/";
+if (scenario==scenarioLow){
+  //Rstar=0.12*365;
+  Tw=2;
+  GDDTHRESHOLD="GDD_294/";
+}
+if (scenario==scenarioMedian){
+  //Rstar=0.12*365;
+  Tw=1;
+  GDDTHRESHOLD="GDD_269/";
+}
+if (scenario==scenarioHight){
+  //Rstar=0.15*365;
+  Tw=0;
+  GDDTHRESHOLD="GDD_246/";
+}
+
+real K=1.0;
+//real Ru2=(Rstar-deathrate)/K;
+real maxGt=0.11*365.0;
+real Ru2=maxGt + deathrate;
+include "defOutputGen.edp";
+real[int] u6(6);
+//real[int] grad6(6);
+include "inputParam.edp";
+cout<<"u6="<<u6<<endl;
+real dirichlet0=1;
+include "defParamsGen.edp";
+include "mat1D2Duseful.edp";
+include "postTraitementLoad.edp";
+//system("export OUTPUTDIR="+OUTPUTDIR+";./outputScript.sh");
+
+real gamma0min=100000;
+real gamma0max=-100000;
+real gamma1min=100000;
+real gamma1max=-100000;
+real w0min=100000;
+real w0max=-100000;
+real w1min=100000;
+real w1max=-100000;
+
+real RRmin=100000;
+real RRmax=-100000;
+
+macro  printInfoYear(XX){
+        XX<<"simulation year "<<CURYEAR<<" ended"<<endl;
+        XX<<"uminDensity="<<uminDensity<<endl;
+        XX<<"umaxDensity="<<umaxDensity<<endl;
+        XX<<"w0min="<<w0min<<endl;
+        XX<<"w0max="<<w0max<<endl;
+        XX<<"w1min="<<w1min<<endl;
+        XX<<"w1max="<<w1max<<endl;
+        XX<<"gamma0min="<<h0*gamma0min<<endl;
+        XX<<"gamma0max="<<h0*gamma0max<<endl;
+        XX<<"gamma1min="<<gamma1min<<endl;
+        XX<<"gamma1max="<<gamma1max<<endl;
+        XX<<"RRmin="<<RRmin<<endl;
+        XX<<"RRmax="<<RRmax<<endl;
+}//
+
+
+
+
+int nbCallSave;
+int numDay=0;
+int lastRRPlot=0;
+func real J(real[int] & u)
+{
+  real coefLength=969000/12;
+  //real coefExpo=2*log(coefLength)/log(10.0);
+  D2DXu0=pow(10,u[0])*coefLength*coefLength;
+  D2DYu0=D2DXu0;
+  D1Du1=pow(10,u[1])*coefLength*coefLength;
+  D1Du0=D1Du1;
+  real R=pow(10,u[2]);
+  muu1=u[3]*u[3];
+  nuu1=coefLength*u[4]*u[4];
+  dirichlet0=0.1*u[5]*u[5]/coefLength;
+  betamodel=coefLength*coefLength;
+  RR0=R;RR1=R;RR2=R;RR3=R;RR4=R;RR5=R;RR6=R;RR7=R;RR8=R;RR9=R;RR10=R;RR11=R;RR12=R;RR13=R;RR14=R;
+  cout<<"PostTrait u("<<u[0]<<","<<u[1]<<","<<u[2]<<","<<u[3]<<","<<u[4]<<","<<u[5]<<")"<<endl;
+  cout<<"PostTrait J("<<D2DXu0<<","<<D2DYu0<<","<<D1Du1<<","<<R<<","<<muu1<<","<<nuu1<<","<<dirichlet0<<")"<<endl;
+  include "defVarfGen.edp";
+  include "defIC.edp";
+  include "defIC1dNull.edp";
+  include "defMatGen.edp";
+
+  //include "buildMatGenW01.edp";
+      
+  cout<<"end ic\n";
+//  ofstream infoFileNull("infosimu.txt");
+  //nbCallSave=0;
+  cout<<"begin loop\n";
+  if (NBPREDYEAR>0)
+    CURYEAR=FIRSTYEAR+NBPREDYEAR-1;
+  else
+    CURYEAR=FIRSTYEAR;
+  
+  for (CURYEAR=CURYEAR; CURYEAR<=LASTYEAR; CURYEAR++){
+    nbCallSave=0;
+    if (CURYEAR>2019){
+      real p;
+      p=randreal1();
+      int rg=p/0.2+1;
+      CURYEARRAND=2008+rg;
+      if (CURYEARRAND>2019)
+        CURYEARRAND=2019;
+      cout<<"CURYEAR "<<CURYEAR<<" use data "<<CURYEARRAND<<endl;
+    }
+//    ofstream infoFile("infosimu.txt",append);
+    if (CURYEAR==2008){
+      //nbSteps=nbSteps2011;
+      T0=1-122.0/365.0;
+    }else{
+      //nbSteps=nbStepsByYear;