--- imach/src/imach.c	2016/02/12 11:29:23	1.218
+++ imach/src/imach.c	2016/02/15 00:48:12	1.219
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.218 2016/02/12 11:29:23 brouard Exp $
+/* $Id: imach.c,v 1.219 2016/02/15 00:48:12 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.219  2016/02/15 00:48:12  brouard
+  *** empty log message ***
+
   Revision 1.218  2016/02/12 11:29:23  brouard
   Summary: 0.99 Back projections
 
@@ -825,12 +828,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.218 2016/02/12 11:29:23 brouard Exp $ */
+/* $Id: imach.c,v 1.219 2016/02/15 00:48:12 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.218 $ $Date: 2016/02/12 11:29:23 $"; 
+char fullversion[]="$Revision: 1.219 $ $Date: 2016/02/15 00:48:12 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -976,7 +979,7 @@ int *ncodemaxwundef;  /* ncodemax[j]= Nu
 			     undefined. Usually 3: -1, 0 and 1. */
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs; /* Global pointer */
-double ***mobaverage; /* New global variable */
+double ***mobaverage, ***mobaverages; /* New global variable */
 double *ageexmed,*agecens;
 double dateintmean=0;
 
@@ -3969,7 +3972,7 @@ void prevalence(double ***probs, double
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
   first=1;
-  for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
+  for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
     for (i=1; i<=nlstate; i++)  
       for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
 				prop[i][iage]=0.0;
@@ -3977,11 +3980,11 @@ void prevalence(double ***probs, double
     for (i=1; i<=imx; i++) { /* Each individual */
       bool=1;
       if  (cptcovn>0) {  /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
-				for (z1=1; z1<=cptcoveff; z1++) 
+				for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
 					if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
 						bool=0;
       } 
-      if (bool==1) { 
+      if (bool==1) { /* For this combination of covariates values, this individual fits */
 				/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
 				for(mi=1; mi<wav[i];mi++){
 					m=mw[mi][i];
@@ -4242,20 +4245,20 @@ void tricode(int *Tvar, int **nbcode, in
   for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
     for (k=-1; k < maxncov; k++) Ndum[k]=0;
     for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the 
-			       modality of this covariate Vj*/ 
+																modality of this covariate Vj*/ 
       ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
-				    * If product of Vn*Vm, still boolean *:
-				    * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
-				    * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
+																		* If product of Vn*Vm, still boolean *:
+																		* If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+																		* 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
       /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
 				      modality of the nth covariate of individual i. */
       if (ij > modmaxcovj)
         modmaxcovj=ij; 
       else if (ij < modmincovj) 
-	modmincovj=ij; 
+				modmincovj=ij; 
       if ((ij < -1) && (ij > NCOVMAX)){
-	printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
-	exit(1);
+				printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+				exit(1);
       }else
       Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
       /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
@@ -4273,19 +4276,19 @@ void tricode(int *Tvar, int **nbcode, in
       printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
       fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
       if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
-	if( k != -1){
-	  ncodemax[j]++;  /* ncodemax[j]= Number of modalities of the j th
-			     covariate for which somebody answered excluding 
-			     undefined. Usually 2: 0 and 1. */
-	}
-	ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
-			     covariate for which somebody answered including 
-			     undefined. Usually 3: -1, 0 and 1. */
+				if( k != -1){
+					ncodemax[j]++;  /* ncodemax[j]= Number of modalities of the j th
+														 covariate for which somebody answered excluding 
+														 undefined. Usually 2: 0 and 1. */
+				}
+				ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
+																covariate for which somebody answered including 
+																undefined. Usually 3: -1, 0 and 1. */
       }
       /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for
-	 historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
+				 historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
     } /* Ndum[-1] number of undefined modalities */
-
+		
     /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
     /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. 
        If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
@@ -4302,8 +4305,8 @@ void tricode(int *Tvar, int **nbcode, in
     ij=0; /* ij is similar to i but can jump over null modalities */
     for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
     	if (Ndum[i] == 0) { /* If nobody responded to this modality k */
-	  break;
-	}
+				break;
+			}
 	ij++;
 	nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
 	cptcode = ij; /* New max modality for covar j */
@@ -4324,28 +4327,28 @@ void tricode(int *Tvar, int **nbcode, in
     /*   }  /\* end of loop on modality k *\/ */
   } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
   
- for (k=-1; k< maxncov; k++) Ndum[k]=0; 
+	for (k=-1; k< maxncov; k++) Ndum[k]=0; 
   
   for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ 
-   /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
-   ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ 
-   Ndum[ij]++; /* Might be supersed V1 + V1*age */
- } 
-
- ij=0;
- for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
-   /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
-   if((Ndum[i]!=0) && (i<=ncovcol)){
-     ij++;
-     /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
-     Tvaraff[ij]=i; /*For printing (unclear) */
-   }else{
-       /* Tvaraff[ij]=0; */
-   }
- }
- /* ij--; */
- cptcoveff=ij; /*Number of total covariates*/
-
+		/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
+		ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ 
+		Ndum[ij]++; /* Might be supersed V1 + V1*age */
+	} 
+	
+	ij=0;
+	for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+		/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
+		if((Ndum[i]!=0) && (i<=ncovcol)){
+			ij++;
+			/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+			Tvaraff[ij]=i; /*For printing (unclear) */
+		}else{
+			/* Tvaraff[ij]=0; */
+		}
+	}
+	/* ij--; */
+	cptcoveff=ij; /*Number of total covariates*/
+	
 }
 
 
@@ -5632,6 +5635,8 @@ true period expectancies (those weighted
   int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
+	int ioffset; /* variable offset for columns */
+
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
 /*     fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
@@ -5661,7 +5666,7 @@ true period expectancies (those weighted
       fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
       fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
       for (j=2; j<= nlstate+ndeath ; j ++) {
-	fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+				fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
       }
       fprintf(ficgp,";\nset out; unset ylabel;\n"); 
     }
@@ -5678,122 +5683,144 @@ true period expectancies (those weighted
     for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
       /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
       fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
-	fprintf(ficgp," V%d=%d ",k,vlv);
+      for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
+				lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+				vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
+			/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+				fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
 
-     fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
-     fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
-     fprintf(ficgp,"set xlabel \"Age\" \n\
-set ylabel \"Probability\" \n\
-set ter svg size 640, 480\n\
+			fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+			fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
+			fprintf(ficgp,"set xlabel \"Age\" \n\
+set ylabel \"Probability\" \n	\
+set ter svg size 640, 480\n	\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-
-     for (i=1; i<= nlstate ; i ++) {
-       if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-       else        fprintf(ficgp," %%*lf (%%*lf)");
-     }
-     fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-     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-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
-     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,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
-		 if(backcast==1){
-			 fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt);
-		 }
-     fprintf(ficgp,"\nset out \n");
+			
+			for (i=1; i<= nlstate ; i ++) {
+				if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+				else        fprintf(ficgp," %%*lf (%%*lf)");
+			}
+			fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+			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-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
+			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,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+			if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
+				/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
+				fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
+				kl=0;
+				for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
+					lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+					/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+					/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+					/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+					vlv= nbcode[Tvaraff[k]][lv];
+					kl++;
+					/* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+					/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
+					/*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(k==cptcoveff){
+							fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv], \
+										4+(cpt-1),  cpt );
+					}else{
+						fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv]);
+						kl++;
+					}
+				} /* end covariate */
+			}
+			fprintf(ficgp,"\nset out \n");
     } /* k1 */
   } /* cpt */
   /*2 eme*/
   for (k1=1; k1<= m ; k1 ++) { 
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
-	fprintf(ficgp," V%d=%d ",k,vlv);
+				lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+				vlv= nbcode[Tvaraff[k]][lv];
+				fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
-
-    fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
-    for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-      if(vpopbased==0)
-	fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-      else
-	fprintf(ficgp,"\nreplot ");
-      for (i=1; i<= nlstate+1 ; i ++) {
-	k=2*i;
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-	else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	fprintf(ficgp,"\" t\"\" w l lt 0,");
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-	else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
-      } /* state */
-    } /* vpopbased */
-    fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+			
+			fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+			for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+				if(vpopbased==0)
+					fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+				else
+					fprintf(ficgp,"\nreplot ");
+				for (i=1; i<= nlstate+1 ; i ++) {
+					k=2*i;
+					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+					for (j=1; j<= nlstate+1 ; j ++) {
+						if (j==i) fprintf(ficgp," %%lf (%%lf)");
+						else fprintf(ficgp," %%*lf (%%*lf)");
+					}   
+					if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+					else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+					for (j=1; j<= nlstate+1 ; j ++) {
+						if (j==i) fprintf(ficgp," %%lf (%%lf)");
+						else fprintf(ficgp," %%*lf (%%*lf)");
+					}   
+					fprintf(ficgp,"\" t\"\" w l lt 0,");
+					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+					for (j=1; j<= nlstate+1 ; j ++) {
+						if (j==i) fprintf(ficgp," %%lf (%%lf)");
+						else fprintf(ficgp," %%*lf (%%*lf)");
+					}   
+					if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+					else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+				} /* state */
+			} /* vpopbased */
+			fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
   } /* k1 */
-
-
+	
+	
   /*3eme*/
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  cov=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
-	fprintf(ficgp," V%d=%d ",k,vlv);
+				lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+				vlv= nbcode[Tvaraff[k]][lv];
+				fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
-
+			
       /*       k=2+nlstate*(2*cpt-2); */
       k=2+(nlstate+1)*(cpt-1);
       fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
       fprintf(ficgp,"set ter svg size 640, 480\n\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-	fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-	
+				for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+				fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+				fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+				for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+				fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+				
       */
       for (i=1; i< nlstate ; i ++) {
-	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
-	/*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
-	
+				fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+				/*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+				
       } 
       fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
     }
@@ -5808,7 +5835,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
+	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
@@ -5844,7 +5871,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
+	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
@@ -5888,7 +5915,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[lv]][lv];
+	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",k,vlv);
       }
       fprintf(ficgp,"\n#\n");
@@ -5923,7 +5950,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[lv]][lv];
+	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",k,vlv);
 	}
 	fprintf(ficgp,"\n#\n");
@@ -5959,81 +5986,100 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
     
     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
-	for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
-	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[lv]][lv];
-	  fprintf(ficgp," V%d=%d ",k,vlv);
-	}
-	fprintf(ficgp,"\n#\n");
-	
-	fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
-	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n\
-unset log y\n\
+				fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
+				for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
+					lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+					/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+					/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+					/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+					vlv= nbcode[Tvaraff[k]][lv];
+					fprintf(ficgp," V%d=%d ",k,vlv);
+				}
+				fprintf(ficgp,"\n#\n");
+				
+				fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+				fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+				fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\n	\
+unset log y\n	\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);
-	for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
-	  /*#  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 */   
-	  /*# 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 */   
-	  if(i==1){
-	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
-	  }else{
-	    fprintf(ficgp,",\\\n '' ");
-	  }
-	  if(cptcoveff ==0){ /* No covariate */
-	    fprintf(ficgp," u 2:("); /* Age is in 2 */
-	    /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
-	    /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
-	    if(i==nlstate+1)
-	      fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
-			2+(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,cpt );
-	    else
-	      fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
-		      2+(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
-	  }else{
-	    fprintf(ficgp,"u 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 */   
-	    kl=0;
-	    for (k=1; k<=cptcoveff; k++){    /* For each covariate  */
-	      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
-	      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	      vlv= nbcode[Tvaraff[lv]][lv];
-	      kl++;
-	      /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
-	      /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
-	      /*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(k==cptcoveff)
-		if(i==nlstate+1)
-		  fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
-			  6+(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,cpt );
-		else
-		  fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
-			  6+(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
-	      else{
-		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
-		kl++;
-	      }
-	    } /* end covariate */
-	  } /* end if covariate */
-	} /* nlstate */
-	fprintf(ficgp,"\nset out\n");
-      } /* end cpt state*/
-    } /* end covariate */
-  } /* End if prevfcast */
-
-
-  /* proba elementaires */
-  fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+				for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
+					/*#  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 */   
+					/*# 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 */   
+					if(i==1){
+						fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+					}else{
+						fprintf(ficgp,",\\\n '' ");
+					}
+					if(cptcoveff ==0){ /* No covariate */
+						ioffset=2; /* Age is in 2 */
+						/*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+						/*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
+						/*# V1  = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+						/*#  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)) t 'pw.%d' with line ",			\
+											ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+						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 */
+						}   
+						fprintf(ficgp," u %d:((",ioffset); 
+						kl=0;
+						for (k=1; k<=cptcoveff; k++){    /* For each covariate  */
+							lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+							/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+							/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+							/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+							vlv= nbcode[Tvaraff[k]][lv];
+							kl++;
+							/* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+							/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
+							/*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(k==cptcoveff){
+								if(i==nlstate+1){
+									if(cptcoveff ==1){
+									fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+													ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+									}else{
+									fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+													ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+									}
+								}else{
+									if(cptcoveff ==1){
+										fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+														ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+									}else{
+										fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+														ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+									}
+								}
+							}else{ /* k < cptcoveff */
+								fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[k]][lv]);
+								kl++;
+							}
+						} /* end covariate */
+					} /* end if covariate */
+				} /* nlstate */
+				fprintf(ficgp,"\nset out\n");
+			} /* end cpt state*/
+		} /* end covariate */
+	} /* End if prevfcast */
+	
+	
+	/* proba elementaires */
+	fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
   for(i=1,jk=1; i <=nlstate; i++){
     fprintf(ficgp,"# initial state %d\n",i);
     for(k=1; k <=(nlstate+ndeath); k++){
@@ -6175,25 +6221,34 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 
 
 /*************** Moving average **************/
-		/* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
+/* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
 int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
    
   int i, cpt, cptcod;
   int modcovmax =1;
   int mobilavrange, mob;
-  double age;
   int iage=0;
+
+  double sum=0.;
+  double age;
   double *sumnewp, *sumnewm;
   double *agemingood, *agemaxgood; /* Currently identical for all covariates */
   
+  
+  modcovmax=2*cptcoveff;/* Max number of modalities. We suppose 
+			   a covariate has 2 modalities, should be equal to ncovcombmax  */
+
   sumnewp = vector(1,modcovmax);
   sumnewm = vector(1,modcovmax);
   agemingood = vector(1,modcovmax);	
   agemaxgood = vector(1,modcovmax);
-  
-  
-  modcovmax=2*cptcoveff;/* Max number of modalities. We suppose 
-			   a covariate has 2 modalities, should be equal to ncovcombmax  */
+
+  for (cptcod=1;cptcod<=modcovmax;cptcod++){
+		sumnewm[cptcod]=0.;
+		sumnewp[cptcod]=0.;
+		agemingood[cptcod]=0;
+		agemaxgood[cptcod]=0;
+	}
   if (cptcovn<1) modcovmax=1; /* At least 1 pass */
   
   if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
@@ -6201,86 +6256,109 @@ int movingaverage(double ***probs, doubl
     else mobilavrange=mobilav;
     for (age=bage; age<=fage; age++)
       for (i=1; i<=nlstate;i++)
-	for (cptcod=1;cptcod<=modcovmax;cptcod++)
-	  mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+				for (cptcod=1;cptcod<=modcovmax;cptcod++)
+					mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
     /* We keep the original values on the extreme ages bage, fage and for 
        fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
        we use a 5 terms etc. until the borders are no more concerned. 
     */ 
     for (mob=3;mob <=mobilavrange;mob=mob+2){
       for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
-	for (i=1; i<=nlstate;i++){
-	  for (cptcod=1;cptcod<=modcovmax;cptcod++){
-	    mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
-	    for (cpt=1;cpt<=(mob-1)/2;cpt++){
-	      mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
-	      mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
-	    }
-	    mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
-	  }
-	}
+				for (i=1; i<=nlstate;i++){
+					for (cptcod=1;cptcod<=modcovmax;cptcod++){
+						mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
+						for (cpt=1;cpt<=(mob-1)/2;cpt++){
+							mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+							mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+						}
+						mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
+					}
+				}
       }/* end age */
     }/* end mob */
   }else
     return -1;
   for (cptcod=1;cptcod<=modcovmax;cptcod++){
     /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
-    agemingood[cptcod]=fage+(mob-1)/2;
+    agemingood[cptcod]=fage-(mob-1)/2;
     for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
       sumnewm[cptcod]=0.;
       for (i=1; i<=nlstate;i++){
         sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
       }
       if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
-	agemingood[cptcod]=age;
-      }else{ /* bad */
-	for (i=1; i<=nlstate;i++){
-	  mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
-	} /* i */
-      } /* end bad */
-    }/* age */
-    /* From youngest, finding the oldest wrong */
-    agemaxgood[cptcod]=bage+(mob-1)/2;
-    for (age=bage+(mob-1)/2; age<=fage; age++){
-      sumnewm[cptcod]=0.;
-      for (i=1; i<=nlstate;i++){
-        sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
-      }
-      if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
-	agemaxgood[cptcod]=age;
+				agemingood[cptcod]=age;
       }else{ /* bad */
-	for (i=1; i<=nlstate;i++){
-	  mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
-	} /* i */
+				for (i=1; i<=nlstate;i++){
+					mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+				} /* i */
       } /* end bad */
     }/* age */
-    for (age=bage; age<=fage; age++){
-      printf("%d %d ", cptcod, (int)age);
-      sumnewp[cptcod]=0.;
-      sumnewm[cptcod]=0.;
-      for (i=1; i<=nlstate;i++){
-	sumnewp[cptcod]+=probs[(int)age][i][cptcod];
-        sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
-	printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]);
-      }
-      printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]);
+    sum=0.;
+    for (i=1; i<=nlstate;i++){
+      sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
     }
-    printf("\n");
+    if(fabs(sum - 1.) > 1.e-3) { /* bad */
+      printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod);
+      /* for (i=1; i<=nlstate;i++){ */
+      /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+      /* } /\* i *\/ */
+    } /* end bad */
+    /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
+		/* From youngest, finding the oldest wrong */
+		agemaxgood[cptcod]=bage+(mob-1)/2;
+		for (age=bage+(mob-1)/2; age<=fage; age++){
+			sumnewm[cptcod]=0.;
+			for (i=1; i<=nlstate;i++){
+				sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+			}
+			if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+				agemaxgood[cptcod]=age;
+			}else{ /* bad */
+				for (i=1; i<=nlstate;i++){
+					mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+				} /* i */
+			} /* end bad */
+		}/* age */
+		sum=0.;
+		for (i=1; i<=nlstate;i++){
+			sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+		}
+		if(fabs(sum - 1.) > 1.e-3) { /* bad */
+			printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod);
+			/* for (i=1; i<=nlstate;i++){ */
+			/*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+			/* } /\* i *\/ */
+		} /* end bad */
+		
+		for (age=bage; age<=fage; age++){
+			printf("%d %d ", cptcod, (int)age);
+			sumnewp[cptcod]=0.;
+			sumnewm[cptcod]=0.;
+			for (i=1; i<=nlstate;i++){
+				sumnewp[cptcod]+=probs[(int)age][i][cptcod];
+				sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+				/* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */
+			}
+			/* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */
+		}
+		/* printf("\n"); */
+    /* } */
     /* brutal averaging */
     for (i=1; i<=nlstate;i++){
       for (age=1; age<=bage; age++){
-	mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
-	printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+				mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+				/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
       }	
       for (age=fage; age<=AGESUP; age++){
-	mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
-	printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+				mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+				/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
       }
     } /* end i status */
     for (i=nlstate+1; i<=nlstate+ndeath;i++){
       for (age=1; age<=AGESUP; age++){
-	/*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
-	mobaverage[(int)age][i][cptcod]=0.;
+				/*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
+				mobaverage[(int)age][i][cptcod]=0.;
       }
     }
   }/* end cptcod */
@@ -6358,54 +6436,54 @@ void prevforecast(char fileres[], double
       k=k+1;
       fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+				fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficresf," yearproj age");
       for(j=1; j<=nlstate+ndeath;j++){ 
-	for(i=1; i<=nlstate;i++) 	      
+				for(i=1; i<=nlstate;i++) 	      
           fprintf(ficresf," p%d%d",i,j);
-	fprintf(ficresf," p.%d",j);
+				fprintf(ficresf," wp.%d",j);
       }
       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--){ 
-	  nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
-	  nhstepm = nhstepm/hstepm; 
-	  p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-	  oldm=oldms;savm=savms;
-	  hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
-	
-	  for (h=0; h<=nhstepm; h++){
-	    if (h*hstepm/YEARM*stepm ==yearp) {
+				fprintf(ficresf,"\n");
+				fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
+				for (agec=fage; agec>=(ageminpar-1); agec--){ 
+					nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
+					nhstepm = nhstepm/hstepm; 
+					p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+					oldm=oldms;savm=savms;
+					hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+					
+					for (h=0; h<=nhstepm; h++){
+						if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
                 fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	      fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
-	    } 
-	    for(j=1; j<=nlstate+ndeath;j++) {
-	      ppij=0.;
-	      for(i=1; i<=nlstate;i++) {
-		if (mobilav==1) 
-		  ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
-		else {
-		  ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
-		}
-		if (h*hstepm/YEARM*stepm== yearp) {
-		  fprintf(ficresf," %.3f", p3mat[i][j][h]);
-		}
-	      } /* end i */
-	      if (h*hstepm/YEARM*stepm==yearp) {
-		fprintf(ficresf," %.3f", ppij);
-	      }
-	    }/* end j */
-	  } /* end h */
-	  free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-	} /* end agec */
+							fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+						} 
+						for(j=1; j<=nlstate+ndeath;j++) {
+							ppij=0.;
+							for(i=1; i<=nlstate;i++) {
+								if (mobilav==1) 
+									ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
+								else {
+									ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
+								}
+								if (h*hstepm/YEARM*stepm== yearp) {
+									fprintf(ficresf," %.3f", p3mat[i][j][h]);
+								}
+							} /* end i */
+							if (h*hstepm/YEARM*stepm==yearp) {
+								fprintf(ficresf," %.3f", ppij);
+							}
+						}/* end j */
+					} /* end h */
+					free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+				} /* end agec */
       } /* end yearp */
     } /* end cptcod */
   } /* end  cptcov */
-       
+	
   fclose(ficresf);
   printf("End of Computing forecasting \n");
   fprintf(ficlog,"End of Computing forecasting\n");
@@ -7656,9 +7734,9 @@ void syscompilerinfo(int logged)
 #endif
    
 
- }
+}
 
- int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
+int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
   /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
   int i, j, k, i1 ;
   /* double ftolpl = 1.e-10; */
@@ -7679,55 +7757,55 @@ void syscompilerinfo(int logged)
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");
   
-    /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+  /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
 
-    agebase=ageminpar;
-    agelim=agemaxpar;
+  agebase=ageminpar;
+  agelim=agemaxpar;
 
-    i1=pow(2,cptcoveff);
-    if (cptcovn < 1){i1=1;}
+  i1=pow(2,cptcoveff);
+  if (cptcovn < 1){i1=1;}
 
-    for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+  for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
-      //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-	k=k+1;
-	/* to clean */
-	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-	fprintf(ficrespl,"#******");
-	printf("#******");
-	fprintf(ficlog,"#******");
-	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	}
-	fprintf(ficrespl,"******\n");
-	printf("******\n");
-	fprintf(ficlog,"******\n");
-
-	fprintf(ficrespl,"#Age ");
-	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	}
-	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
-	fprintf(ficrespl,"Total Years_to_converge\n");
-	
-	for (age=agebase; age<=agelim; age++){
-	/* for (age=agebase; age<=agebase; age++){ */
-	  prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
-	  fprintf(ficrespl,"%.0f ",age );
-	  for(j=1;j<=cptcoveff;j++)
-	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-	  tot=0.;
-	  for(i=1; i<=nlstate;i++){
-	    tot +=  prlim[i][i];
-	    fprintf(ficrespl," %.5f", prlim[i][i]);
-	  }
-	  fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
-	} /* Age */
-	/* was end of cptcod */
-    } /* cptcov */
-	return 0;
+    //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+    k=k+1;
+    /* to clean */
+    //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+    fprintf(ficrespl,"#******");
+    printf("#******");
+    fprintf(ficlog,"#******");
+    for(j=1;j<=cptcoveff;j++) {
+      fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+    }
+    fprintf(ficrespl,"******\n");
+    printf("******\n");
+    fprintf(ficlog,"******\n");
+
+    fprintf(ficrespl,"#Age ");
+    for(j=1;j<=cptcoveff;j++) {
+      fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+    }
+    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
+    fprintf(ficrespl,"Total Years_to_converge\n");
+	
+    for (age=agebase; age<=agelim; age++){
+      /* for (age=agebase; age<=agebase; age++){ */
+      prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
+      fprintf(ficrespl,"%.0f ",age );
+      for(j=1;j<=cptcoveff;j++)
+	fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      tot=0.;
+      for(i=1; i<=nlstate;i++){
+	tot +=  prlim[i][i];
+	fprintf(ficrespl," %.5f", prlim[i][i]);
+      }
+      fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
+    } /* Age */
+    /* was end of cptcod */
+  } /* cptcov */
+  return 0;
 }
 
 int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){
@@ -7798,22 +7876,22 @@ int back_prevalence_limit(double *p, dou
       if(mobilavproj > 0){
 	/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
 	/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-	bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+				bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
       }else if (mobilavproj == 0){
-	printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-	fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-	exit(1);
+				printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+				fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+				exit(1);
       }else{
-	/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-	bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+				/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+				bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
       }
       fprintf(ficresplb,"%.0f ",age );
       for(j=1;j<=cptcoveff;j++)
-	fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+				fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       tot=0.;
       for(i=1; i<=nlstate;i++){
-	tot +=  bprlim[i][i];
-	fprintf(ficresplb," %.5f", bprlim[i][i]);
+				tot +=  bprlim[i][i];
+				fprintf(ficresplb," %.5f", bprlim[i][i]);
       }
       fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
     } /* Age */
@@ -8004,7 +8082,7 @@ int main(int argc, char *argv[])
   double agedeb=0.;
 
   double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
-	double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
+  double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
 
   double fret;
   double dum=0.; /* Dummy variable */
@@ -8702,11 +8780,11 @@ Please run with mle=-1 to get a correct
      *              bbbbbbbb
      *              76543210     
      *   h-1        00000101 (6-1=5)
-     *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+     *(h-1)>>(k-1)= 00000010 >> (2-1) = 1 right shift
      *           &
      *     1        00000001 (1)
-     *              00000001        = 1 & ((h-1) >> (k-1))
-     *          +1= 00000010 =2 
+     *              00000000        = 1 & ((h-1) >> (k-1))
+     *          +1= 00000001 =1 
      *
      * h=14, k=3 => h'=h-1=13, k'=k-1=2
      *          h'      1101 =2^3+2^2+0x2^1+2^0
@@ -9284,8 +9362,8 @@ Please run with mle=-1 to get a correct
     if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
 
     if (num_filled != 6) {
-      printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
-      printf("but line=%s\n",line);
+      printf("Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
+      fprintf(ficlog,"Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
       goto end;
     }
     printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
@@ -9360,9 +9438,9 @@ Please run with mle=-1 to get a correct
     ungetc(c,ficpar);
     
     fscanf(ficpar,"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);
-    fscanf(ficparo,"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);
-    fscanf(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);
-    fscanf(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);
+    fprintf(ficparo,"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(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.*/
     
     
@@ -9413,36 +9491,41 @@ Please run with mle=-1 to get a correct
     hPijx(p, bage, fage);
     fclose(ficrespij);
 
-		ncovcombmax=  pow(2,cptcoveff);
-  /*-------------- Variance of one-step probabilities---*/
+    ncovcombmax=  pow(2,cptcoveff);
+    /*-------------- Variance of one-step probabilities---*/
     k=1;
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
 
-		/* Prevalence for each covariates in probs[age][status][cov] */
+    /* Prevalence for each covariates in probs[age][status][cov] */
     probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
     for(i=1;i<=AGESUP;i++)
-      for(j=1;j<=nlstate;j++)
+      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);
-		if (mobilav!=0 ||mobilavproj !=0 ) {
-			mobaverage= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
-			if (mobilav!=0) {
+    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++)
+				for(j=1;j<=nlstate;j++)
+					for(k=1;k<=ncovcombmax;k++)
+						mobaverages[i][j][k]=0.;
+      mobaverage=mobaverages;
+      if (mobilav!=0) {
 				if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
 					fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
 					printf(" Error in movingaverage mobilav=%d\n",mobilav);
 				}
-			}
-			/* /\* 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) {
+      }
+      /* /\* 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) {
 				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);
 				}
-			}
-		}/* end if moving average */
-
+      }
+    }/* end if moving average */
+		
     /*---------- Forecasting ------------------*/
     /*if((stepm == 1) && (strcmp(model,".")==0)){*/
     if(prevfcast==1){
@@ -9450,60 +9533,36 @@ Please run with mle=-1 to get a correct
       prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
     }
     if(backcast==1){
-		ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
-		ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
-		ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
-
-			/*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
-			/*#include "prevlim.h"*/  /* Use ficresplb, ficlog */
-			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, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
-			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);
-		}
-  /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-  /* if (mobilavproj!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-
-    /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
-    /*      }  */
-    /*      else{ */
-    /*        erreur=108; */
-    /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
-    /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
-    /*      } */
+      ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
+      ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); 	
+      ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+
+      /*--------------- 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, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
+	 bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+      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);
+    }
     
  
     /* ------ Other prevalence ratios------------ */
 
-    /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
-
-    /* prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-    /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
-	ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
-    */
     free_ivector(wav,1,imx);
     free_imatrix(dh,1,lastpass-firstpass+2,1,imx);
     free_imatrix(bh,1,lastpass-firstpass+2,1,imx);
     free_imatrix(mw,1,lastpass-firstpass+2,1,imx);   
 		
 		
-		/*   if (mobilav!=0) { */
-		/*     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-		/*     if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
-		/* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
-		/* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
-		/*     } */
-		/*   } */
-		
-		
     /*---------- Health expectancies, no variances ------------*/
 		
     strcpy(filerese,"E_");
@@ -9514,22 +9573,19 @@ Please run with mle=-1 to get a correct
     }
     printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
-    /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-		
+ 		
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
-			fprintf(ficreseij,"\n#****** ");
-			for(j=1;j<=cptcoveff;j++) {
-				fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-			}
-			fprintf(ficreseij,"******\n");
-			
-			eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
-			oldm=oldms;savm=savms;
-			evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
+      fprintf(ficreseij,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      fprintf(ficreseij,"******\n");
       
-			free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-      /*}*/
+      eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+      oldm=oldms;savm=savms;
+      evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
+      
+      free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
     }
     fclose(ficreseij);
     printf("done evsij\n");fflush(stdout);
@@ -9615,57 +9671,57 @@ Please run with mle=-1 to get a correct
       
       
       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 */
-				printf("varevsij %d \n",vpopbased);
-				fprintf(ficlog, "varevsij %d \n",vpopbased);
-				varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
-				fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
-				if(vpopbased==1)
-					fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
-				else
-					fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
-				fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
-				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++){
-					prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
-					if (vpopbased==1) {
-						if(mobilav ==0){
-							for(i=1; i<=nlstate;i++)
-								prlim[i][i]=probs[(int)age][i][k];
-						}else{ /* mobilav */ 
-							for(i=1; i<=nlstate;i++)
-								prlim[i][i]=mobaverage[(int)age][i][k];
-						}
-					}
-					
-					fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
-					/* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
-					/* printf(" age %4.0f ",age); */
-					for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
-						for(i=1, epj[j]=0.;i <=nlstate;i++) {
-							epj[j] += prlim[i][i]*eij[i][j][(int)age];
-							/*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
-							/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
-						}
-						epj[nlstate+1] +=epj[j];
-					}
-					/* printf(" age %4.0f \n",age); */
-					
-					for(i=1, vepp=0.;i <=nlstate;i++)
-						for(j=1;j <=nlstate;j++)
-							vepp += vareij[i][j][(int)age];
-					fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
-					for(j=1;j <=nlstate;j++){
-						fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
-					}
-					fprintf(ficrest,"\n");
-				}
+	oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+	cptcod= 0; /* To be deleted */
+	printf("varevsij %d \n",vpopbased);
+	fprintf(ficlog, "varevsij %d \n",vpopbased);
+	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
+	if(vpopbased==1)
+	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+	else
+	  fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+	fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+	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++){
+	  prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+	  if (vpopbased==1) {
+	    if(mobilav ==0){
+	      for(i=1; i<=nlstate;i++)
+		prlim[i][i]=probs[(int)age][i][k];
+	    }else{ /* mobilav */ 
+	      for(i=1; i<=nlstate;i++)
+		prlim[i][i]=mobaverage[(int)age][i][k];
+	    }
+	  }
+	  
+	  fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+	  /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+	  /* printf(" age %4.0f ",age); */
+	  for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+	    for(i=1, epj[j]=0.;i <=nlstate;i++) {
+	      epj[j] += prlim[i][i]*eij[i][j][(int)age];
+	      /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+	      /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
+	    }
+	    epj[nlstate+1] +=epj[j];
+	  }
+	  /* printf(" age %4.0f \n",age); */
+	  
+	  for(i=1, vepp=0.;i <=nlstate;i++)
+	    for(j=1;j <=nlstate;j++)
+	      vepp += vareij[i][j][(int)age];
+	  fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+	  for(j=1;j <=nlstate;j++){
+	    fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+	  }
+	  fprintf(ficrest,"\n");
+	}
       } /* End vpopbased */
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
@@ -9722,7 +9778,8 @@ Please run with mle=-1 to get a correct
     fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
 
     /*---------- End : free ----------------*/
-    if (mobilav!=0 ||mobilavproj !=0) free_ma3x(mobaverage,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
+    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);
   }  /* mle==-3 arrives here for freeing */
  /* endfree:*/