--- imach/src/imach.c	2016/08/26 09:20:19	1.237
+++ imach/src/imach.c	2016/08/26 14:23:35	1.238
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.237 2016/08/26 09:20:19 brouard Exp $
+/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.238  2016/08/26 14:23:35  brouard
+  Summary: Starting tests of 0.99
+
   Revision 1.237  2016/08/26 09:20:19  brouard
   Summary: to valgrind
 
@@ -911,12 +914,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.237 2016/08/26 09:20:19 brouard Exp $ */
+/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 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.237 $ $Date: 2016/08/26 09:20:19 $"; 
+char fullversion[]="$Revision: 1.238 $ $Date: 2016/08/26 14:23:35 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1145,6 +1148,8 @@ double *Tvalsel; /**< Selected modality
 int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
 int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ 
 int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ 
+int *DummyV; /** Dummy[v] 0=dummy (0 1), 1 quantitative */
+int *FixedV; /** FixedV[v] 0 fixed, 1 varying */
 int *Tage;
 int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */ 
 int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
@@ -1154,11 +1159,10 @@ int *Ndum; /** Freq of modality (tricode
 /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard;
 int *Tprod;/**< Gives the k position of the k1 product */
+/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3  */
 int *Tposprod; /**< Gives the k1 product from the k position */
-/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
-   if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2)
-   Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2 
-*/
+   /* if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2) */
+   /* Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5(V3*V2)]=2 (2nd product without age) */
 int cptcovprod, *Tvaraff, *invalidvarcomb;
 double *lsurv, *lpop, *tpop;
 
@@ -4918,7 +4922,7 @@ void  concatwav(int wav[], int **dh, int
   double ***p3mat;
   double eip;
 
-  pstamp(ficreseij);
+  /* pstamp(ficreseij); */
   fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
   fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++){
@@ -5283,6 +5287,14 @@ void  concatwav(int wav[], int **dh, int
    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
    pstamp(ficresprobmorprev);
    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
+   fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");
+   for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+     fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+   }
+   for(j=1;j<=cptcoveff;j++) 
+     fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,j)]);
+   fprintf(ficresprobmorprev,"\n");
+
    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
    for(j=nlstate+1; j<=(nlstate+ndeath);j++){
      fprintf(ficresprobmorprev," p.%-d SE",j);
@@ -6270,158 +6282,99 @@ void printinggnuplot(char fileresu[], ch
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
   /* 1eme*/
-  for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
-    for (k1=1; k1<= m ; k1 ++) /* For each valid combination of covariate */
+  for (cpt=1; cpt<= nlstate ; cpt ++){ /* For each live state */
+    for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-      if(TKresult[nres]!= k1)
-	continue;
-      /* We are interested in selected combination by the resultline */
-      printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
-      fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
-      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 */
-	printf(" V%d=%d ",Tvaraff[k],vlv);
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      }	
-      printf("\n#\n");
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-	continue;
-      }
+	/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+	if(TKresult[nres]!= k1)
+	  continue;
+	/* We are interested in selected combination by the resultline */
+	printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
+	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
+	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 */
+	  printf(" V%d=%d ",Tvaraff[k],vlv);
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	printf("\n#\n");
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
       
-      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\
+	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){ /* 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 */
-	if(cptcoveff ==0){
-	  fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",	 2+(cpt-1),  cpt );
-	}else{
-	  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' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
-		      4+(cpt-1),  cpt );  /* 4 or 6 ?*/
-	    }else{
-	      fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+	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 */
+	  if(cptcoveff ==0){
+	    fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",	 2+(cpt-1),  cpt );
+	  }else{
+	    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++;
-	    }
-	  } /* end covariate */
-	} /* end if no covariate */
-      } /* end if backcast */
-      fprintf(ficgp,"\nset out \n");
+	      /* 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' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+			4+(cpt-1),  cpt );  /* 4 or 6 ?*/
+	      }else{
+		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+		kl++;
+	      }
+	    } /* end covariate */
+	  } /* end if no covariate */
+	} /* end if backcast */
+	fprintf(ficgp,"\nset out \n");
+      } /* nres */
     } /* k1 */
   } /* cpt */
 
   
   /*2 eme*/
-  for (k1=1; k1<= m ; k1 ++)  
-  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k1)
-      continue;
-    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[k]][lv];
-      fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-    }
-    for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-      printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-    }
-    fprintf(ficgp,"\n#\n");
-    if(invalidvarcomb[k1]){
-      fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-      continue;
-    }
-			
-    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 end 2 eme*/
-	
-	
-  /*3eme*/
-  for (k1=1; k1<= m ; k1 ++) 
-  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k)
-      continue;
-
-    for (cpt=1; cpt<= nlstate ; cpt ++) {
-      fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate dummy combination and each value */
+  for (k1=1; k1<= m ; k1 ++){  
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
+	continue;
+      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 */
@@ -6431,139 +6384,197 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
       }
       /* for(k=1; k <= ncovds; k++){ */
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      }	
+      }
       fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){
 	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
 	continue;
       }
 			
-      /*       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\
+      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 */
+    } /* end nres */
+  } /* k1 end 2 eme*/
+	
+	
+  /*3eme*/
+  for (k1=1; k1<= m ; k1 ++){
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k)
+	continue;
+
+      for (cpt=1; cpt<= nlstate ; cpt ++) {
+	fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%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[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
+			
+	/*       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);
+	/*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 ; 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);*/
+	*/
+	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.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
-    }
-  }
+	} 
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
+      }
+    } /* end nres */
+  } /* end kl 3eme */
   
   /* 4eme */
   /* Survival functions (period) from state i in state j by initial state i */
-  for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
-  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k1)
-      continue;
-
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' 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[k]][lv];
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      }	
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+  for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
 	continue;
-      }
-			
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
-set ter svg size 640, 480\n																							\
-unset log y\n																														\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
-      k=3;
-      for (i=1; i<= nlstate ; i ++){
-	if(i==1){
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-	}else{
-	  fprintf(ficgp,", '' ");
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
+	fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' 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[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	}
-	l=(nlstate+ndeath)*(i-1)+1;
-	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
-	for (j=2; j<= nlstate+ndeath ; j ++)
-	  fprintf(ficgp,"+$%d",k+l+j-1);
-	fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
-      } /* nlstate */
-      fprintf(ficgp,"\nset out\n");
-    } /* end cpt state*/ 
-  } /* end covariate */  
-	
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
+      
+	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
+	k=3;
+	for (i=1; i<= nlstate ; i ++){
+	  if(i==1){
+	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	  }else{
+	    fprintf(ficgp,", '' ");
+	  }
+	  l=(nlstate+ndeath)*(i-1)+1;
+	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+	  for (j=2; j<= nlstate+ndeath ; j ++)
+	    fprintf(ficgp,"+$%d",k+l+j-1);
+	  fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+	} /* nlstate */
+	fprintf(ficgp,"\nset out\n");
+      } /* end cpt state*/ 
+    } /* end nres */
+  } /* end covariate k1 */  
+
 /* 5eme */
   /* Survival functions (period) from state i in state j by final state j */
-  for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
-  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k1)
-      continue;
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-			
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' 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[k]][lv];
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      }	
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+  for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
 	continue;
-      }
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+	fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' 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[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
       
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
-      k=3;
-      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-	if(j==1)
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-	else
-	  fprintf(ficgp,", '' ");
-	l=(nlstate+ndeath)*(cpt-1) +j;
-	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
-	/* for (i=2; i<= nlstate+ndeath ; i ++) */
-	/*   fprintf(ficgp,"+$%d",k+l+i-1); */
-	fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
-      } /* nlstate */
-      fprintf(ficgp,", '' ");
-      fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
-      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-	l=(nlstate+ndeath)*(cpt-1) +j;
-	if(j < nlstate)
-	  fprintf(ficgp,"$%d +",k+l);
-	else
-	  fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
-      }
-      fprintf(ficgp,"\nset out\n");
-    } /* end cpt state*/ 
-  } /* end covariate */  
+	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
+	k=3;
+	for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	  if(j==1)
+	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	  else
+	    fprintf(ficgp,", '' ");
+	  l=(nlstate+ndeath)*(cpt-1) +j;
+	  fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+	  /* for (i=2; i<= nlstate+ndeath ; i ++) */
+	  /*   fprintf(ficgp,"+$%d",k+l+i-1); */
+	  fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+	} /* nlstate */
+	fprintf(ficgp,", '' ");
+	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+	for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	  l=(nlstate+ndeath)*(cpt-1) +j;
+	  if(j < nlstate)
+	    fprintf(ficgp,"$%d +",k+l);
+	  else
+	    fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+	}
+	fprintf(ficgp,"\nset out\n");
+      } /* end cpt state*/ 
+    } /* end covariate */  
+  } /* end nres */
   
 /* 6eme */
   /* CV preval stable (period) for each covariate */
@@ -6593,9 +6604,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
       
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){
 	if(i==1)
@@ -6641,9 +6650,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	
 	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3; /* Offset */
 	for (i=1; i<= nlstate ; i ++){
 	  if(i==1)
@@ -6695,9 +6702,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	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);
+set ter svg size 640, 480\nunset log y\nplot [%.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 */   
@@ -6764,8 +6769,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* End if prevfcast */
   
   
-  /* proba elementaires */
-  fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+  /* 9eme writing MLE parameters */
+  fprintf(ficgp,"\n##############\n#9eme 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++){
@@ -6782,7 +6787,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"##############\n#\n");
   
   /*goto avoid;*/
-  fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
+  /* 10eme Graphics of probabilities or incidences using written MLE parameters */
+  fprintf(ficgp,"\n##############\n#10eme Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
   fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
   fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -6797,9 +6803,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
   fprintf(ficgp,"#\n");
   for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
-    fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year \n");
+    fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
     fprintf(ficgp,"#model=%s \n",model);
-    fprintf(ficgp,"# ng=%d\n",ng);
+    fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
     for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
@@ -6856,7 +6862,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	      /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
 	      if(j==Tage[ij]) { /* Product by age */
 		if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
-		  if(Dummy[j]==0){
+		  if(DummyV[j]==0){
 		    fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
 		  }else{ /* quantitative */
 		    fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
@@ -6867,8 +6873,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	      }else if(j==Tprod[ijp]) { /* */ 
 		/* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
 		if(ijp <=cptcovprod) { /* Product */
-		  if(Dummy[Tvard[ijp][1]]==0){/* Vn is dummy */
-		    if(Dummy[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+		  if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+		    if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
 		      /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],nbcode[Tvard[ijp][2]][codtabm(jk,j)]); */
 		      fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
 		    }else{ /* Vn is dummy and Vm is quanti */
@@ -6876,12 +6882,13 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 		      fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
 		    }
 		  }else{ /* Vn*Vm Vn is quanti */
-		    if(Dummy[Tvard[ijp][2]]==0){
+		    if(DummyV[Tvard[ijp][2]]==0){
 		      fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
 		    }else{ /* Both quanti */
 		      fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
 		    }
 		  }
+		  ijp++;
   		}
 	      } else{  /* simple covariate */
 		/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
@@ -7175,8 +7182,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
     }
     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-      printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
     }
     fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){ 
@@ -7822,87 +7828,87 @@ int readdata(char datafile[], int firsto
     /* Loops on waves */
     for (j=maxwav;j>=1;j--){
       for (iv=nqtv;iv>=1;iv--){  /* Loop  on time varying quantitative variables */
-				cutv(stra, strb, line, ' '); 
-				if(strb[0]=='.') { /* Missing value */
-					lval=-1;
-					cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
-					cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
-					if(isalpha(strb[1])) { /* .m or .d Really Missing value */
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
-						return 1;
-					}
-				}else{
-					errno=0;
-					/* what_kind_of_number(strb); */
-					dval=strtod(strb,&endptr); 
-					/* if( strb[0]=='\0' || (*endptr != '\0')){ */
-					/* if(strb != endptr && *endptr == '\0') */
-					/*    dval=dlval; */
-					/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
-					if( strb[0]=='\0' || (*endptr != '\0')){
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
-						return 1;
-					}
-					cotqvar[j][iv][i]=dval; 
-					cotvar[j][ntv+iv][i]=dval; 
-				}
-				strcpy(line,stra);
+	cutv(stra, strb, line, ' '); 
+	if(strb[0]=='.') { /* Missing value */
+	  lval=-1;
+	  cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
+	  cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
+	  if(isalpha(strb[1])) { /* .m or .d Really Missing value */
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
+	    return 1;
+	  }
+	}else{
+	  errno=0;
+	  /* what_kind_of_number(strb); */
+	  dval=strtod(strb,&endptr); 
+	  /* if( strb[0]=='\0' || (*endptr != '\0')){ */
+	  /* if(strb != endptr && *endptr == '\0') */
+	  /*    dval=dlval; */
+	  /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
+	  if( strb[0]=='\0' || (*endptr != '\0')){
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
+	    return 1;
+	  }
+	  cotqvar[j][iv][i]=dval; 
+	  cotvar[j][ntv+iv][i]=dval; 
+	}
+	strcpy(line,stra);
       }/* end loop ntqv */
       
       for (iv=ntv;iv>=1;iv--){  /* Loop  on time varying dummies */
-				cutv(stra, strb, line, ' '); 
-				if(strb[0]=='.') { /* Missing value */
-					lval=-1;
-				}else{
-					errno=0;
-					lval=strtol(strb,&endptr,10); 
-					/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-					if( strb[0]=='\0' || (*endptr != '\0')){
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
-						return 1;
-					}
-				}
-				if(lval <-1 || lval >1){
-					printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+	cutv(stra, strb, line, ' '); 
+	if(strb[0]=='.') { /* Missing value */
+	  lval=-1;
+	}else{
+	  errno=0;
+	  lval=strtol(strb,&endptr,10); 
+	  /*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+	  if( strb[0]=='\0' || (*endptr != '\0')){
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
+	    return 1;
+	  }
+	}
+	if(lval <-1 || lval >1){
+	  printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
- For example, for multinomial values like 1, 2 and 3,\n									\
- build V1=0 V2=0 for the reference value (1),\n													\
-        V1=1 V2=0 for (2) \n																						\
+ For example, for multinomial values like 1, 2 and 3,\n			\
+ build V1=0 V2=0 for the reference value (1),\n				\
+        V1=1 V2=0 for (2) \n						\
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
- output of IMaCh is often meaningless.\n																\
+ output of IMaCh is often meaningless.\n				\
  Exiting.\n",lval,linei, i,line,j);
-					fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+	  fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
- For example, for multinomial values like 1, 2 and 3,\n									\
- build V1=0 V2=0 for the reference value (1),\n													\
-        V1=1 V2=0 for (2) \n																						\
+ For example, for multinomial values like 1, 2 and 3,\n			\
+ build V1=0 V2=0 for the reference value (1),\n				\
+        V1=1 V2=0 for (2) \n						\
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
- output of IMaCh is often meaningless.\n																\
+ output of IMaCh is often meaningless.\n				\
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
-					return 1;
-				}
-				cotvar[j][iv][i]=(double)(lval);
-				strcpy(line,stra);
+	  return 1;
+	}
+	cotvar[j][iv][i]=(double)(lval);
+	strcpy(line,stra);
       }/* end loop ntv */
       
       /* Statuses  at wave */
       cutv(stra, strb, line, ' '); 
       if(strb[0]=='.') { /* Missing value */
-				lval=-1;
+	lval=-1;
       }else{
-				errno=0;
-				lval=strtol(strb,&endptr,10); 
-				/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-				if( strb[0]=='\0' || (*endptr != '\0')){
-					printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
-					fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
-					return 1;
-				}
+	errno=0;
+	lval=strtol(strb,&endptr,10); 
+	/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+	if( strb[0]=='\0' || (*endptr != '\0')){
+	  printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+	  fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+	  return 1;
+	}
       }
       
       s[j][i]=lval;
@@ -8207,7 +8213,7 @@ int decodemodel( char model[], int lasto
 	* - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
 	*/
 {
-  int i, j, k, ks;
+  int i, j, k, ks, v;
   int  j1, k1, k2, k3, k4;
   char modelsav[80];
   char stra[80], strb[80], strc[80], strd[80],stre[80];
@@ -8414,6 +8420,26 @@ Typevar: 0 for simple covariate (dummy,
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
 
+  for(v=1; v <=ncovcol;v++){
+    DummyV[v]=0;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+1; v <=ncovcol+nqv;v++){
+    DummyV[v]=1;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
+    DummyV[v]=0;
+    FixedV[v]=1;
+  }
+  for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    DummyV[v]=1;
+    FixedV[v]=1;
+  }
+  for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+    fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+  }
   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
     if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;
@@ -9028,61 +9054,63 @@ int prevalence_limit(double *p, double *
   i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
   if (cptcovn < 1){i1=1;}
 
-  for(nres=1; nres <= nresult; nres++) /* For each resultline */
-  for(k=1; k<=i1;k++){
-    if(TKresult[nres]!= k)
-      continue;
+  for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k)
+	continue;
 
-  /* 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++) {/* all covariates */
-      fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
-      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)]);
-    }
-    for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-      printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-    }
-    fprintf(ficrespl,"******\n");
-    printf("******\n");
-    fprintf(ficlog,"******\n");
-    if(invalidvarcomb[k]){
-      printf("\nCombination (%d) ignored because no case \n",k); 
-      fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k); 
-      fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k); 
-						continue;
-    }
+      /* 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++) {/* all covariates */
+	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
+	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)]);
+      }
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficrespl," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
+      fprintf(ficrespl,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
+      if(invalidvarcomb[k]){
+	printf("\nCombination (%d) ignored because no case \n",k); 
+	fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k); 
+	fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k); 
+	continue;
+      }
 
-    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");
+      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, nres);
-      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 */
+      for (age=agebase; age<=agelim; age++){
+	/* for (age=agebase; age<=agebase; age++){ */
+	prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres);
+	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 */
+  } /* nres */
   return 0;
 }
 
@@ -9123,69 +9151,69 @@ int back_prevalence_limit(double *p, dou
   i1=pow(2,cptcoveff);
   if (cptcovn < 1){i1=1;}
   
-  for(nres=1; nres <= nresult; nres++) /* For each resultline */
-  for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-    if(TKresult[nres]!= k)
-      continue;
-    //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-    fprintf(ficresplb,"#******");
-    printf("#******");
-    fprintf(ficlog,"#******");
-    for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-      fprintf(ficresplb," 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)]);
-    }
-    for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-      printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-      fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-      fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-    }
-    fprintf(ficresplb,"******\n");
-    printf("******\n");
-    fprintf(ficlog,"******\n");
-    if(invalidvarcomb[k]){
-      printf("\nCombination (%d) ignored because no cases \n",k); 
-      fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k); 
-      fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k); 
-      continue;
-    }
+  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
+      //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+      fprintf(ficresplb,"#******");
+      printf("#******");
+      fprintf(ficlog,"#******");
+      for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+	fprintf(ficresplb," 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)]);
+      }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }
+      fprintf(ficresplb,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
+      if(invalidvarcomb[k]){
+	printf("\nCombination (%d) ignored because no cases \n",k); 
+	fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k); 
+	fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k); 
+	continue;
+      }
     
-    fprintf(ficresplb,"#Age ");
-    for(j=1;j<=cptcoveff;j++) {
-      fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    }
-    for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
-    fprintf(ficresplb,"Total Years_to_converge\n");
+      fprintf(ficresplb,"#Age ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
+      fprintf(ficresplb,"Total Years_to_converge\n");
     
     
-    for (age=agebase; age<=agelim; age++){
-      /* for (age=agebase; age<=agebase; age++){ */
-      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);
-      }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);
-      }else{
-	/* 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)]);
-      tot=0.;
-      for(i=1; i<=nlstate;i++){
-	tot +=  bprlim[i][i];
-	fprintf(ficresplb," %.5f", bprlim[i][i]);
-      }
-      fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
-    } /* Age */
-    /* was end of cptcod */
-  } /* cptcov */
-  
+      for (age=agebase; age<=agelim; age++){
+	/* for (age=agebase; age<=agebase; age++){ */
+	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);
+	}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);
+	}else{
+	  /* 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)]);
+	tot=0.;
+	for(i=1; i<=nlstate;i++){
+	  tot +=  bprlim[i][i];
+	  fprintf(ficresplb," %.5f", bprlim[i][i]);
+	}
+	fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+      } /* Age */
+      /* was end of cptcod */
+    } /* end of any combination */
+  } /* end of nres */  
   /* hBijx(p, bage, fage); */
   /* fclose(ficrespijb); */
   
@@ -9277,7 +9305,7 @@ int hPijx(double *p, int bage, int fage)
 	int ageminl;
   int hstepm;
   int nhstepm;
-  int h, i, i1, j, k;
+  int h, i, i1, j, k, nres;
 	
   double agedeb;
   double ***p3mat;
@@ -9305,48 +9333,54 @@ int hPijx(double *p, int bage, int fage)
   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
   /*  	k=k+1;  */
-  for (k=1; k <= (int) pow(2,cptcoveff); k++){
-    fprintf(ficrespijb,"\n#****** ");
-    for(j=1;j<=cptcoveff;j++)
-      fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    fprintf(ficrespijb,"******\n");
-    if(invalidvarcomb[k]){
-      fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); 
-      continue;
-    }
-    
-    /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
-    for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
-      /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
-      nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-      nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
-      
-      /*	  nhstepm=nhstepm*YEARM; aff par mois*/
-      
-      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-      /* oldm=oldms;savm=savms; */
-      /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
-      hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
-      /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
-      fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
-      for(i=1; i<=nlstate;i++)
-	for(j=1; j<=nlstate+ndeath;j++)
-	  fprintf(ficrespijb," %1d-%1d",i,j);
-      fprintf(ficrespijb,"\n");
-      for (h=0; h<=nhstepm; h++){
-	/*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
-	fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
-	/* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
+      fprintf(ficrespijb,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++)
+	fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }
+      fprintf(ficrespijb,"******\n");
+      if(invalidvarcomb[k]){
+	fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); 
+	continue;
+      }
+      
+      /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+      for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
+	/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
+	nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+	nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
+	
+	/*	  nhstepm=nhstepm*YEARM; aff par mois*/
+	
+	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	/* oldm=oldms;savm=savms; */
+	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
+	hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
+	fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
 	for(i=1; i<=nlstate;i++)
 	  for(j=1; j<=nlstate+ndeath;j++)
-	    fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+	    fprintf(ficrespijb," %1d-%1d",i,j);
 	fprintf(ficrespijb,"\n");
-      }
-      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-      fprintf(ficrespijb,"\n");
-    }
-    /*}*/
-  }
+	for (h=0; h<=nhstepm; h++){
+	  /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+	  fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
+	  /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+	  for(i=1; i<=nlstate;i++)
+	    for(j=1; j<=nlstate+ndeath;j++)
+	      fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+	  fprintf(ficrespijb,"\n");
+	}
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	fprintf(ficrespijb,"\n");
+      } /* end age deb */
+    } /* end combination */
+  } /* end nres */
   return 0;
  } /*  hBijx */
 
@@ -9975,6 +10009,8 @@ Please run with mle=-1 to get a correct
   Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
   Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
   Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
+  DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
+  FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
   /*  V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). 
       For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, 
       Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
@@ -10751,6 +10787,7 @@ Please run with mle=-1 to get a correct
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);
+      fputs(line,ficres);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
@@ -10767,6 +10804,7 @@ Please run with mle=-1 to get a correct
       fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);
       fputs(line,ficparo);
+      fputs(line,ficres);
     }
     ungetc(c,ficpar);
     
@@ -10785,6 +10823,7 @@ Please run with mle=-1 to get a correct
 	fputs(line,stdout);
 	fputs(line,ficparo);
 	fputs(line,ficlog);
+	fputs(line,ficres);
 	continue;
       }else
 	break;
@@ -10803,12 +10842,16 @@ Please run with mle=-1 to get a correct
 	goto end;
       }
       decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+      fprintf(ficparo,"result: %s\n",resultline);
+      fprintf(ficres,"result: %s\n",resultline);
+      fprintf(ficlog,"result: %s\n",resultline);
       while(fgets(line, MAXLINE, ficpar)) {
 	/* If line starts with a # it is a comment */
 	if (line[0] == '#') {
 	  numlinepar++;
 	  fputs(line,stdout);
 	  fputs(line,ficparo);
+	  fputs(line,ficres);
 	  fputs(line,ficlog);
 	  continue;
 	}else
@@ -10953,6 +10996,8 @@ 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);
+
+    pstamp(ficreseij);
  		
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}
@@ -11065,6 +11110,7 @@ Please run with mle=-1 to get a correct
       fprintf(ficrescveij,"******\n");
       
       fprintf(ficresvij,"\n#****** ");
+      /* pstamp(ficresvij); */
       for(j=1;j<=cptcoveff;j++) 
 	fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
@@ -11242,6 +11288,8 @@ Please run with mle=-1 to get a correct
   free_ivector(ncodemaxwundef,1,NCOVMAX);
   free_ivector(Dummy,-1,NCOVMAX);
   free_ivector(Fixed,-1,NCOVMAX);
+  free_ivector(DummyV,1,NCOVMAX);
+  free_ivector(FixedV,1,NCOVMAX);
   free_ivector(Typevar,-1,NCOVMAX);
   free_ivector(Tvar,1,NCOVMAX);
   free_ivector(TvarsQ,1,NCOVMAX);