--- imach/src/imach.c	2017/05/18 20:09:32	1.268
+++ imach/src/imach.c	2017/06/27 11:06:02	1.273
@@ -1,6 +1,21 @@
-/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $
+/* $Id: imach.c,v 1.273 2017/06/27 11:06:02 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.273  2017/06/27 11:06:02  brouard
+  Summary: More documentation on projections
+
+  Revision 1.272  2017/06/27 10:22:40  brouard
+  Summary: Color of backprojection changed from 6 to 5(yellow)
+
+  Revision 1.271  2017/06/27 10:17:50  brouard
+  Summary: Some bug with rint
+
+  Revision 1.270  2017/05/24 05:45:29  brouard
+  *** empty log message ***
+
+  Revision 1.269  2017/05/23 08:39:25  brouard
+  Summary: Code into subroutine, cleanings
+
   Revision 1.268  2017/05/18 20:09:32  brouard
   Summary: backprojection and confidence intervals of backprevalence
 
@@ -1007,12 +1022,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $ */
+/* $Id: imach.c,v 1.273 2017/06/27 11:06:02 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.268 $ $Date: 2017/05/18 20:09:32 $"; 
+char fullversion[]="$Revision: 1.273 $ $Date: 2017/06/27 11:06:02 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1084,10 +1099,7 @@ FILE *ficrescveij;
 char filerescve[FILENAMELENGTH];
 FILE  *ficresvij;
 char fileresv[FILENAMELENGTH];
-FILE  *ficresvpl;
-char fileresvpl[FILENAMELENGTH];
-FILE  *ficresvbl;
-char fileresvbl[FILENAMELENGTH];
+
 char title[MAXLINE];
 char model[MAXLINE]; /**< The model line */
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];
@@ -2932,14 +2944,15 @@ double **pmij(double **ps, double *cov,
   /* Diag(w_x) */
   /* Problem with prevacurrent which can be zero */
   sumnew=0.;
-  /* for (ii=1;ii<=nlstate+ndeath;ii++){ */
+  /*for (ii=1;ii<=nlstate+ndeath;ii++){*/
   for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
+    /* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]);  */
     sumnew+=prevacurrent[(int)agefin][ii][ij];
   }
   if(sumnew >0.01){  /* At least some value in the prevalence */
     for (ii=1;ii<=nlstate+ndeath;ii++){
       for (j=1;j<=nlstate+ndeath;j++)
-      doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
+	doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
     }
   }else{
     for (ii=1;ii<=nlstate+ndeath;ii++){
@@ -3239,7 +3252,7 @@ double ***hbxij(double ***po, int nhstep
       newm=savm;
       /* Covariates have to be included here again */
       cov[1]=1.;
-      agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */
+      agexact=age-( (h-1)*hstepm + (d)  )*stepm/YEARM; /* age just before transition, d or d-1? */
       /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
       cov[2]=agexact;
       if(nagesqr==1)
@@ -5463,7 +5476,7 @@ void  concatwav(int wav[], int **dh, int
   /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. 
      nhstepm is the number of hstepm from age to agelim 
      nstepm is the number of stepm from age to agelin. 
-     Look at hpijx to understand the reason of that which relies in memory size
+     Look at hpijx to understand the reason which relies in memory size consideration
      and note for a fixed period like estepm months */
   /* We decided (b) to get a life expectancy respecting the most precise curvature of the
      survival function given by stepm (the optimization length). Unfortunately it
@@ -5694,7 +5707,8 @@ void  concatwav(int wav[], int **dh, int
 	  /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
 					
 	}
-		
+
+    /* Standard deviation of expectancies ij */		
     fprintf(ficresstdeij,"%3.0f",age );
     for(i=1; i<=nlstate;i++){
       eip=0.;
@@ -5709,6 +5723,7 @@ void  concatwav(int wav[], int **dh, int
     }
     fprintf(ficresstdeij,"\n");
 		
+    /* Variance of expectancies ij */		
     fprintf(ficrescveij,"%3.0f",age );
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){
@@ -6062,7 +6077,7 @@ void  concatwav(int wav[], int **dh, int
  }  /* end varevsij */
 
 /************ Variance of prevlim ******************/
- void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres)
+ void varprevlim(char fileresvpl[], FILE *ficresvpl, double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres)
 {
   /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
@@ -6187,7 +6202,7 @@ void  concatwav(int wav[], int **dh, int
 
 
 /************ Variance of backprevalence limit ******************/
- void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)
+ void varbrevlim(char fileresvbl[], FILE  *ficresvbl, double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)
 {
   /* Variance of backward prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
@@ -6652,8 +6667,8 @@ void printinghtml(char fileresu[], char
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
 		  int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
-		  double jprev1, double mprev1,double anprev1, double dateprev1, \
-		  double jprev2, double mprev2,double anprev2, double dateprev2){
+		  double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \
+		  double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){
   int jj1, k1, i1, cpt, k4, nres;
 
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
@@ -6807,15 +6822,18 @@ divided by h: <sub>h</sub>P<sub>ij</sub>
      if(prevfcast==1){
        /* Projection of prevalence up to period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }
      }
      if(backcast==1){
       /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence  at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+	 fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
+ from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
+ account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+ <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
        }
      }
 	 
@@ -6921,7 +6939,7 @@ true period expectancies (those weighted
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
 
   char dirfileres[132],optfileres[132];
   char gplotcondition[132], gplotlabel[132];
@@ -6930,6 +6948,8 @@ void printinggnuplot(char fileresu[], ch
   int ng=0;
   int vpopbased;
   int ioffset; /* variable offset for columns */
+  int iyearc=1; /* variable column for year of projection  */
+  int iagec=1; /* variable column for age of projection  */
   int nres=0; /* Index of resultline */
   int istart=1; /* For starting graphs in projections */
 
@@ -7034,7 +7054,7 @@ void printinggnuplot(char fileresu[], ch
 	
 	fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_"));
         if(cptcoveff ==0){
-	  fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3",	 2+(cpt-1),  cpt );
+	  fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3",	 2+3*(cpt-1),  cpt );
 	}else{
 	  kl=0;
 	  for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
@@ -7092,17 +7112,17 @@ void printinggnuplot(char fileresu[], ch
 	      if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	      else        fprintf(ficgp," %%*lf (%%*lf)");
 	    }
-	    fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+	    fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 6 dt 3,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
 	    for (i=1; i<= nlstate ; i ++) {
 	      if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	      else fprintf(ficgp," %%*lf (%%*lf)");
 	    } 
-	    fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); 
+	    fprintf(ficgp,"\" t\"95%% CI\" w l lt 5,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); 
 	    for (i=1; i<= nlstate ; i ++) {
 	      if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	      else fprintf(ficgp," %%*lf (%%*lf)");
 	    } 
-	    fprintf(ficgp,"\" t\"\" w l lt 1");
+	    fprintf(ficgp,"\" t\"\" w l lt 5");
 	  } /* end if backprojcast */
 	} /* end if backcast */
 	fprintf(ficgp,"\nset out ;unset label;\n");
@@ -7503,24 +7523,22 @@ set ter svg size 640, 480\nunset log y\n
 	    /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
 	    fprintf(ficgp," u %d:(", ioffset); 
 	    if(i==nlstate+1){
-	      fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ",	\
+	      fprintf(ficgp," $%d/(1.-$%d)):1 t 'pw.%d' with line lc variable ",	\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
 	      fprintf(ficgp,",\\\n '' ");
 	      fprintf(ficgp," u %d:(",ioffset); 
-	      fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \
+	      fprintf(ficgp," (($1-$2) == %d ) ? $%d/(1.-$%d) : 1/0):1 with labels center not ", \
 		     offyear,				\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate );
 	    }else
 	      fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ",	\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
 	  }else{ /* more than 2 covariates */
-	    if(cptcoveff ==1){
-	      ioffset=4; /* Age is in 4 */
-	    }else{
-	      ioffset=6; /* Age is in 6 */
-	      /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
-	      /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
-	    }   
+	    ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/
+	    /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+	    /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+	    iyearc=ioffset-1;
+	    iagec=ioffset;
 	    fprintf(ficgp," u %d:(",ioffset); 
 	    kl=0;
 	    strcpy(gplotcondition,"(");
@@ -7542,13 +7560,13 @@ set ter svg size 640, 480\nunset log y\n
 	    /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
 	    /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
 	    if(i==nlstate+1){
-	      fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \
-		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+	      fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
+		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );
 	      fprintf(ficgp,",\\\n '' ");
-	      fprintf(ficgp," u %d:(",ioffset); 
-	      fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \
-		     offyear,				\
-		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate );
+	      fprintf(ficgp," u %d:(",iagec); 
+	      fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
+		      iyearc, iagec, offyear,				\
+		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );
 /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
 	    }else{
 	      fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
@@ -7618,24 +7636,22 @@ set ter svg size 640, 480\nunset log y\n
 	    /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
 	    fprintf(ficgp," u %d:(", ioffset); 
 	    if(i==nlstate+1){
-	      fprintf(ficgp," $%d/(1.-$%d)):5 t 'bw%d' with line lc variable ",	\
+	      fprintf(ficgp," $%d/(1.-$%d)):1 t 'bw%d' with line lc variable ",	\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
 	      fprintf(ficgp,",\\\n '' ");
 	      fprintf(ficgp," u %d:(",ioffset); 
-	      fprintf(ficgp," (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", \
+	      fprintf(ficgp," (($1-$2) == %d ) ? $%d : 1/0):1 with labels center not ", \
 		     offbyear,				\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
 	    }else
 	      fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ",	\
 		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i );
 	  }else{ /* more than 2 covariates */
-	    if(cptcoveff ==1){
-	      ioffset=4; /* Age is in 4 */
-	    }else{
-	      ioffset=6; /* Age is in 6 */
-	      /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
-	      /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
-	    }   
+	    ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/
+	    /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+	    /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+	    iyearc=ioffset-1;
+	    iagec=ioffset;
 	    fprintf(ficgp," u %d:(",ioffset); 
 	    kl=0;
 	    strcpy(gplotcondition,"(");
@@ -7657,14 +7673,14 @@ set ter svg size 640, 480\nunset log y\n
 	    /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
 	    /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
 	    if(i==nlstate+1){
-	      fprintf(ficgp,"%s ? $%d : 1/0):5 t 'bw%d' with line lc variable", gplotcondition, \
-		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt );
+	      fprintf(ficgp,"%s ? $%d : 1/0):%d t 'bw%d' with line lc variable", gplotcondition, \
+		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1),iyearc,cpt );
 	      fprintf(ficgp,",\\\n '' ");
-	      fprintf(ficgp," u %d:(",ioffset); 
+	      fprintf(ficgp," u %d:(",iagec); 
 	      /* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */
-	      fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \
-		     offbyear,				\
-		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
+	      fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d : 1/0):%d with labels center not ", gplotcondition, \
+		      iyearc,iagec,offbyear,				\
+		      ioffset+(cpt-1)*(nlstate+1)+1+(i-1), iyearc );
 /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
 	    }else{
 	      /* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */
@@ -8112,7 +8128,7 @@ set ter svg size 640, 480\nunset log y\n
  
 
 /************** Forecasting ******************/
-void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+ void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection 
      agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed
@@ -8151,7 +8167,12 @@ void prevforecast(char fileres[], double
   if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);
   }
-  else  hstepm=estepm;   
+  else{
+    hstepm=estepm;   
+  }
+  if(estepm > stepm){ /* Yes every two year */
+    stepsize=2;
+  }
 
   hstepm=hstepm/stepm; 
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
@@ -8196,7 +8217,8 @@ void prevforecast(char fileres[], double
     for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
       fprintf(ficresf,"\n");
       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
-      for (agec=fage; agec>=(ageminpar-1); agec--){ 
+      /* for (agec=fage; agec>=(ageminpar-1); agec--){  */
+      for (agec=fage; agec>=(bage); agec--){ 
 	nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
 	nhstepm = nhstepm/hstepm; 
 	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
@@ -8218,7 +8240,7 @@ void prevforecast(char fileres[], double
 	  ppij=0.;
 	  for(i=1; i<=nlstate;i++) {
 	    /* if (mobilav>=1)  */
-	    ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+	    ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k];
 	    /* else { */ /* even if mobilav==-1 we use mobaverage */
 	    /* 	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
 	    /* } */
@@ -8239,12 +8261,13 @@ void prevforecast(char fileres[], double
 
 }
 
-/* /\************** Back Forecasting ******************\/ */
-void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
+/************** Back Forecasting ******************/
+ void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
   /* back1, year, month, day of starting backection
      agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed
-     anback2 year of en of backection (same day and month as back1).
+     anback2 year of end of backprojection (same day and month as back1).
+     prevacurrent and prev are prevalences.
   */
   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
   double agec; /* generic age */
@@ -8283,7 +8306,12 @@ void prevbackforecast(char fileres[], do
   if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);
   }
-  else  hstepm=estepm;
+  else{
+    hstepm=estepm;   
+  }
+  if(estepm >= stepm){ /* Yes every two year */
+    stepsize=2;
+  }
   
   hstepm=hstepm/stepm;
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
@@ -8304,10 +8332,6 @@ void prevbackforecast(char fileres[], do
   
   fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
   
-  /* 	      if (h==(int)(YEARM*yearp)){ */
-  /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
-  /*   for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
-  /*     k=k+1; */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */
   for(k=1; k<=i1;k++){
     if(i1 != 1 && TKresult[nres]!= k)
@@ -8333,16 +8357,16 @@ void prevbackforecast(char fileres[], do
       /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  */
       fprintf(ficresfb,"\n");
       fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
-      printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
-      /* for (agec=fage; agec>=(ageminpar-1); agec--){ */
-      /* 	nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
-      for (agec=bage; agec<=agemax-1; agec++){  /* testing */
+      /* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
+      /* for (agec=bage; agec<=agemax-1; agec++){  /\* testing *\/ */
+      for (agec=bage; agec<=fage; agec++){  /* testing */
 	/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
-	nhstepm=(int) rint((agec-agelim)*YEARM/stepm);
+	nhstepm=(int) (agec-agelim) *YEARM/stepm;/*	nhstepm=(int) rint((agec-agelim)*YEARM/stepm);*/
 	nhstepm = nhstepm/hstepm;
 	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
 	oldm=oldms;savm=savms;
 	/* computes hbxij at age agec over 1 to nhstepm */
+	/* printf("####prevbackforecast debug  agec=%.2f nhstepm=%d\n",agec, nhstepm);fflush(stdout); */
 	hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
 	/* hpxij(p3mat,nhstepm,agec,hstepm,p,             nlstate,stepm,oldm,savm, k,nres); */
 	/* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
@@ -8360,8 +8384,10 @@ void prevbackforecast(char fileres[], do
 	  ppij=0.;ppi=0.;
 	  for(j=1; j<=nlstate;j++) {
 	    /* if (mobilav==1) */
-	    ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k];
-	    ppi=ppi+mobaverage[(int)agec][j][k];
+	    ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k];
+	    ppi=ppi+prevacurrent[(int)agec][j][k];
+	    /* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; */
+	    /* ppi=ppi+mobaverage[(int)agec][j][k]; */
 	      /* else { */
 	      /* 	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
 	      /* } */
@@ -8386,6 +8412,123 @@ void prevbackforecast(char fileres[], do
 	
 }
 
+/* Variance of prevalence limit: varprlim */
+ void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){
+    /*------- Variance of period (stable) prevalence------*/   
+ 
+   char fileresvpl[FILENAMELENGTH];  
+   FILE *ficresvpl;
+   double **oldm, **savm;
+   double **varpl; /* Variances of prevalence limits by age */   
+   int i1, k, nres, j ;
+   
+    strcpy(fileresvpl,"VPL_");
+    strcat(fileresvpl,fileresu);
+    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
+      printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
+      exit(0);
+    }
+    printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
+    fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
+    
+    /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
+    
+    i1=pow(2,cptcoveff);
+    if (cptcovn < 1){i1=1;}
+
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
+    for(k=1; k<=i1;k++){
+      if(i1 != 1 && TKresult[nres]!= k)
+	continue;
+      fprintf(ficresvpl,"\n#****** ");
+      printf("\n#****** ");
+      fprintf(ficlog,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }	
+      fprintf(ficresvpl,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
+      
+      varpl=matrix(1,nlstate,(int) bage, (int) fage);
+      oldm=oldms;savm=savms;
+      varprevlim(fileresvpl, ficresvpl, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyearp, k, strstart, nres);
+      free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
+      /*}*/
+    }
+    
+    fclose(ficresvpl);
+    printf("done variance-covariance of period prevalence\n");fflush(stdout);
+    fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
+
+ }
+/* Variance of back prevalence: varbprlim */
+ void varbprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **bprlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){
+      /*------- Variance of back (stable) prevalence------*/
+
+   char fileresvbl[FILENAMELENGTH];  
+   FILE  *ficresvbl;
+
+   double **oldm, **savm;
+   double **varbpl; /* Variances of back prevalence limits by age */   
+   int i1, k, nres, j ;
+
+   strcpy(fileresvbl,"VBL_");
+   strcat(fileresvbl,fileresu);
+   if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
+     printf("Problem with variance of back (stable) prevalence  resultfile: %s\n", fileresvbl);
+     exit(0);
+   }
+   printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
+   fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);
+   
+   
+   i1=pow(2,cptcoveff);
+   if (cptcovn < 1){i1=1;}
+   
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+     for(k=1; k<=i1;k++){
+       if(i1 != 1 && TKresult[nres]!= k)
+	 continue;
+       fprintf(ficresvbl,"\n#****** ");
+       printf("\n#****** ");
+       fprintf(ficlog,"\n#****** ");
+       for(j=1;j<=cptcoveff;j++) {
+	 fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	 fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	 printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       }
+       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	 printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	 fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	 fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+       }
+       fprintf(ficresvbl,"******\n");
+       printf("******\n");
+       fprintf(ficlog,"******\n");
+       
+       varbpl=matrix(1,nlstate,(int) bage, (int) fage);
+       oldm=oldms;savm=savms;
+       
+       varbrevlim(fileresvbl, ficresvbl, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, ncvyearp, k, strstart, nres);
+       free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
+       /*}*/
+     }
+   
+   fclose(ficresvbl);
+   printf("done variance-covariance of back prevalence\n");fflush(stdout);
+   fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);
+
+ } /* End of varbprlim */
+
 /************** Forecasting *****not tested NB*************/
 /* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
   
@@ -10494,11 +10637,12 @@ int main(int argc, char *argv[])
   double *delti; /* Scale */
   double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */
-  double **varbpl; /* Variances of back prevalence limits by age */
+
   double *epj, vepp;
 
-  double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
-  double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
+  double dateprev1, dateprev2;
+  double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0;
+  double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0;
 
   double **ximort;
   char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
@@ -11889,6 +12033,9 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
 	  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
 	  /* day and month of proj2 are not used but only year anproj2.*/
+	  dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.;
+	  dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.;
+
 	}
 	break;
       case 12:
@@ -11904,6 +12051,8 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
 	  fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
 	  /* day and month of proj2 are not used but only year anproj2.*/
+	  dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.;
+	  dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.;
 	}
 	break;
       case 13:
@@ -11954,11 +12103,12 @@ Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else{
-      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1);
+      /* printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1); */
+      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, backcast, pathc,p, (int)anproj1-bage, (int)anback1-fage);
     }
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
 		 model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
-		 jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
+		 jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2);
 		
     /*------------ free_vector  -------------*/
     /*  chdir(path); */
@@ -11994,16 +12144,17 @@ Please run with mle=-1 to get a correct
     k=1;
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
     
-    /* Prevalence for each covariates in probs[age][status][cov] */
-    probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
-    for(i=1;i<=AGESUP;i++)
+    /* Prevalence for each covariate combination in probs[age][status][cov] */
+    probs= ma3x(AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+    for(i=AGEINF;i<=AGESUP;i++)
       for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */
 	for(k=1;k<=ncovcombmax;k++)
 	  probs[i][j][k]=0.;
-    prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+    prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, 
+	       ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     if (mobilav!=0 ||mobilavproj !=0 ) {
-      mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
-      for(i=1;i<=AGESUP;i++)
+      mobaverages= ma3x(AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+      for(i=AGEINF;i<=AGESUP;i++)
 	for(j=1;j<=nlstate+ndeath;j++)
 	  for(k=1;k<=ncovcombmax;k++)
 	    mobaverages[i][j][k]=0.;
@@ -12015,31 +12166,27 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
 	  printf(" Error in movingaverage mobilav=%d\n",mobilav);
 	}
-      }
-      /* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */
-      /* 	for(i=1;i<=AGESUP;i++) */
-      /* 	  for(j=1;j<=nlstate;j++) */
-      /* 	    for(k=1;k<=ncovcombmax;k++) */
-      /* 	      mobaverages[i][j][k]=probs[i][j][k]; */
-      /* 	/\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */
-      /* 	/\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */
-      /* } */
-      else if (mobilavproj !=0) {
+      } else if (mobilavproj !=0) {
 	printf("Movingaveraging projected observed prevalence\n");
 	fprintf(ficlog,"Movingaveraging projected observed prevalence\n");
 	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
 	  fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
 	  printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
 	}
+      }else{
+	printf("Internal error moving average\n");
+	fflush(stdout);
+	exit(1);
       }
     }/* end if moving average */
     
     /*---------- Forecasting ------------------*/
-    /*if((stepm == 1) && (strcmp(model,".")==0)){*/
     if(prevfcast==1){
       /*    if(stepm ==1){*/
-      prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+      prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
     }
+
+    /* Backcasting */
     if(backcast==1){
       ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
       ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
@@ -12048,71 +12195,23 @@ Please run with mle=-1 to get a correct
       /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
 
       bprlim=matrix(1,nlstate,1,nlstate);
+
       back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
       fclose(ficresplb);
 
       hBijx(p, bage, fage, mobaverage);
       fclose(ficrespijb);
-      /* free_matrix(bprlim,1,nlstate,1,nlstate); /\*here or after loop ? *\/ */
 
-      prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
-	 bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+      prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2,
+      		       mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+      varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
 
-      /*------- Variance of back (stable) prevalence------*/   
       
-      strcpy(fileresvbl,"VBL_");
-      strcat(fileresvbl,fileresu);
-      if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
-	printf("Problem with variance of back (stable) prevalence  resultfile: %s\n", fileresvbl);
-	exit(0);
-      }
-      printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
-      fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);
-      
-      /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-      
-      i1=pow(2,cptcoveff);
-      if (cptcovn < 1){i1=1;}
-      
-      for(nres=1; nres <= nresult; nres++) /* For each resultline */
-	for(k=1; k<=i1;k++){
-	  if(i1 != 1 && TKresult[nres]!= k)
-	    continue;
-	  fprintf(ficresvbl,"\n#****** ");
-	  printf("\n#****** ");
-	  fprintf(ficlog,"\n#****** ");
-	  for(j=1;j<=cptcoveff;j++) {
-	    fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	    fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	    printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	  }
-	  for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	    printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	    fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	    fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	  }	
-	  fprintf(ficresvbl,"******\n");
-	  printf("******\n");
-	  fprintf(ficlog,"******\n");
-	  
-	  varbpl=matrix(1,nlstate,(int) bage, (int) fage);
-	  oldm=oldms;savm=savms;
-	  
-	  varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres);
-	  free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
-	  /*}*/
-	}
-    
-      fclose(ficresvbl);
-      printf("done variance-covariance of back prevalence\n");fflush(stdout);
-      fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);
-
+      free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
       free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
-    }
-    free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+    }    /* end  Backcasting */
  
  
     /* ------ Other prevalence ratios------------ */
@@ -12165,10 +12264,10 @@ Please run with mle=-1 to get a correct
     fclose(ficreseij);
     printf("done evsij\n");fflush(stdout);
     fprintf(ficlog,"done evsij\n");fflush(ficlog);
+
 		
     /*---------- State-specific expectancies and variances ------------*/
 		
-		
     strcpy(filerest,"T_");
     strcat(filerest,fileresu);
     if((ficrest=fopen(filerest,"w"))==NULL) {
@@ -12177,8 +12276,6 @@ Please run with mle=-1 to get a correct
     }
     printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
-		
-
     strcpy(fileresstde,"STDE_");
     strcat(fileresstde,fileresu);
     if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
@@ -12206,9 +12303,6 @@ Please run with mle=-1 to get a correct
     printf("      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);
     fprintf(ficlog,"      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);
 
-    /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-          
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}
     
@@ -12270,7 +12364,7 @@ Please run with mle=-1 to get a correct
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       pstamp(ficrest);
       
-      
+      epj=vector(1,nlstate+1);
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	oldm=oldms;savm=savms; /* ZZ Segmentation fault */
 	cptcod= 0; /* To be deleted */
@@ -12286,7 +12380,6 @@ Please run with mle=-1 to get a correct
 	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	fprintf(ficrest,"\n");
 	/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
-	epj=vector(1,nlstate+1);
 	printf("Computing age specific period (stable) prevalences in each health state \n");
 	fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
 	for(age=bage; age <=fage ;age++){
@@ -12324,66 +12417,19 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficrest,"\n");
 	}
       } /* End vpopbased */
+      free_vector(epj,1,nlstate+1);
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-      free_vector(epj,1,nlstate+1);
       printf("done selection\n");fflush(stdout);
       fprintf(ficlog,"done selection\n");fflush(ficlog);
       
-      /*}*/
     } /* End k selection */
 
     printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
 
-    /*------- Variance of period (stable) prevalence------*/   
-    
-    strcpy(fileresvpl,"VPL_");
-    strcat(fileresvpl,fileresu);
-    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
-      printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
-      exit(0);
-    }
-    printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
-    fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
-    
-    /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-    
-    i1=pow(2,cptcoveff);
-    if (cptcovn < 1){i1=1;}
-
-    for(nres=1; nres <= nresult; nres++) /* For each resultline */
-    for(k=1; k<=i1;k++){
-      if(i1 != 1 && TKresult[nres]!= k)
-	continue;
-      fprintf(ficresvpl,"\n#****** ");
-      printf("\n#****** ");
-      fprintf(ficlog,"\n#****** ");
-      for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      }
-      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-      }	
-      fprintf(ficresvpl,"******\n");
-      printf("******\n");
-      fprintf(ficlog,"******\n");
-      
-      varpl=matrix(1,nlstate,(int) bage, (int) fage);
-      oldm=oldms;savm=savms;
-      varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres);
-      free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
-      /*}*/
-    }
-    
-    fclose(ficresvpl);
-    printf("done variance-covariance of period prevalence\n");fflush(stdout);
-    fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
+    /* variance-covariance of period prevalence*/
+    varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
 
     
     free_vector(weight,1,n);
@@ -12402,8 +12448,8 @@ Please run with mle=-1 to get a correct
     
     /*---------- End : free ----------------*/
     if (mobilav!=0 ||mobilavproj !=0)
-      free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
-    free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+      free_ma3x(mobaverages,AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
+    free_ma3x(probs,AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
     free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
   }  /* mle==-3 arrives here for freeing */
@@ -12420,6 +12466,7 @@ Please run with mle=-1 to get a correct
   /*free_vector(delti,1,npar);*/
   free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
   free_matrix(agev,1,maxwav,1,imx);
+  free_ma3x(paramstart,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
   free_ivector(ncodemax,1,NCOVMAX);