--- imach/src/imach.c	2015/11/18 17:41:20	1.210
+++ imach/src/imach.c	2015/11/21 12:41:11	1.211
@@ -1,6 +1,11 @@
-/* $Id: imach.c,v 1.210 2015/11/18 17:41:20 brouard Exp $
+/* $Id: imach.c,v 1.211 2015/11/21 12:41:11 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.211  2015/11/21 12:41:11  brouard
+  Summary: 0.98r3 with some graph of projected cross-sectional
+
+  Author: Nicolas Brouard
+
   Revision 1.210  2015/11/18 17:41:20  brouard
   Summary: Start working on projected prevalences
 
@@ -757,6 +762,8 @@ typedef struct {
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
 #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
+/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
+#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 
 #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
@@ -773,12 +780,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.210 2015/11/18 17:41:20 brouard Exp $ */
+/* $Id: imach.c,v 1.211 2015/11/21 12:41:11 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.210 $ $Date: 2015/11/18 17:41:20 $"; 
+char fullversion[]="$Revision: 1.211 $ $Date: 2015/11/21 12:41:11 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2722,7 +2729,7 @@ void likelione(FILE *ficres,double p[],
     
       
     for (k=1; k<= nlstate ; k++) {
-      fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
+      fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
 <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
     }
     fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
@@ -4804,7 +4811,7 @@ To be simple, these graphs help to under
 void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
-		  int popforecast, int estepm ,\
+		  int popforecast, int prevfcast, int estepm ,		\
 		  double jprev1, double mprev1,double anprev1, \
 		  double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;
@@ -4822,12 +4829,14 @@ void printinghtml(char fileresu[], char
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
 	   subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
    fprintf(fichtm,"\
- - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
+ - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n",
 	   estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
-   fprintf(fichtm,"\
- - Population projections by age and states: \
+   if(prevfcast==1){
+     fprintf(fichtm,"\
+ - Prevalence projections by age and states:				\
    <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
+   }
 
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
@@ -4847,16 +4856,16 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      /* aij, bij */
-     fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
-<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+     fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
      /* Pij */
-     fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+     fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
 <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
      /* Quasi-incidences */
-     fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+     fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
  before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too,\
- incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
-divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \
+divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
 <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
      /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){
@@ -4874,6 +4883,14 @@ divided by h: hPij/h : <a href=\"%s_%d-3
        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
      }
+    if(prevfcast==1){
+      /* Projection of prevalence up to period (stable) prevalence in each health state */
+      for(cpt=1; cpt<=nlstate;cpt++){
+	fprintf(fichtm,"<br>\n- Projection of prevalece up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+      }
+    }
+
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
 <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
@@ -4961,10 +4978,11 @@ true period expectancies (those weighted
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+    void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+  int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
@@ -4977,11 +4995,6 @@ void printinggnuplot(char fileresu[], ch
     /*#endif */
   m=pow(2,cptcoveff);
 
-  /* Projected Prevalences */
-/* plot "NAGI0w_V1V2_monthlyb2b-proj/F_NAGI0w_V1V2_monthlyb2b-proj.txt" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $7/(1-$13):1/0) t 'p11' w line */
-/* replot ""  u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $8/(1-$14):1/0) t 'p21' w line */
-/* replot ""  u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0)&&($9!=0))? $9/(1-$15):1/0) t 'p.1' w line */
-
   /* Contribution to likelihood */
   /* Plot the probability implied in the likelihood */
     fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
@@ -5014,9 +5027,20 @@ void printinggnuplot(char fileresu[], ch
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
  /* 1eme*/
-  fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
-  for (cpt=1; cpt<= nlstate ; cpt ++) {
-    for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+  for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
+    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);
+      }
+      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\
@@ -5043,8 +5067,18 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
     } /* k1 */
   } /* cpt */
   /*2 eme*/
-  fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
   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);
+      }
+      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)
@@ -5077,10 +5111,22 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
     } /* 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);
+      }
+      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);
@@ -5106,13 +5152,23 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-      k=3;
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, 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[lv]][lv];
+	fprintf(ficgp," V%d=%d ",k,vlv);
+      }
+      fprintf(ficgp,"\n#\n");
+
       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_"));
@@ -5131,13 +5187,23 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-      k=3;
       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[lv]][lv];
+	fprintf(ficgp," V%d=%d ",k,vlv);
+      }
+      fprintf(ficgp,"\n#\n");
+
       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_"));
@@ -5163,15 +5229,25 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* end covariate */  
 
   /* CV preval stable (period) for each covariate */
-  for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+  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 */
-      k=3;
-      fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%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);
+      }
+      fprintf(ficgp,"\n#\n");
+
       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);
+      k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){
 	if(i==1)
 	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
@@ -5187,6 +5263,84 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
     } /* end cpt state*/ 
   } /* end covariate */  
 
+  if(prevfcast==1){
+  /* Projection from cross-sectional to stable (period) for each covariate */
+
+    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\
+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,jk=1; i <=nlstate; i++){
@@ -5378,6 +5532,10 @@ void prevforecast(char fileres[], double
   char fileresf[FILENAMELENGTH];
 
   agelim=AGESUP;
+  /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+     in each health status at the date of interview (if between dateprev1 and dateprev2).
+     We still use firstpass and lastpass as another selection.
+  */
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
  
   strcpy(fileresf,"F_"); 
@@ -5428,12 +5586,11 @@ void prevforecast(char fileres[], double
   for(cptcov=1, k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;
-      fprintf(ficresf,"\n#******");
+      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, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
-      fprintf(ficresf,"******\n");
-      fprintf(ficresf,"# Covariate valuofcovar yearproj age");
+      fprintf(ficresf," yearproj age");
       for(j=1; j<=nlstate+ndeath;j++){ 
 	for(i=1; i<=nlstate;i++) 	      
           fprintf(ficresf," p%d%d",i,j);
@@ -7382,13 +7539,20 @@ Please run with mle=-1 to get a correct
   Ndum =ivector(-1,NCOVMAX);  
   if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
     tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
-  /* Nbcode gives the value of the lth modality of jth covariate, in
+  /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
      V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
-  /* 1 to ncodemax[j] is the maximum value of this jth covariate */
+  /* 1 to ncodemax[j] which is the maximum value of this jth covariate */
 
   /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
   /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
+  /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j, 
+   * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded 
+   * (currently 0 or 1) in the data.
+   * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of 
+   * corresponding modality (h,j).
+   */
+
   h=0;
 
 
@@ -7398,8 +7562,9 @@ Please run with mle=-1 to get a correct
   m=pow(2,cptcoveff);
  
 	  /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
-	   * For k=4 covariates, h goes from 1 to 2**k
-	   * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
+	   * For k=4 covariates, h goes from 1 to m=2**k
+	   * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;
+           * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 	   *     h\k   1     2     3     4
 	   *______________________________  
 	   *     1 i=1 1 i=1 1 i=1 1 i=1 1
@@ -7419,6 +7584,49 @@ Please run with mle=-1 to get a correct
 	   *    15 i=8 1     2     2     2
 	   *    16     2     2     2     2
 	   */
+  /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates?
+     /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
+     * and the value of each covariate?
+     * V1=1, V2=1, V3=2, V4=1 ?
+     * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok.
+     * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st.
+     * In order to get the real value in the data, we use nbcode
+     * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0
+     * We are keeping this crazy system in order to be able (in the future?) 
+     * to have more than 2 values (0 or 1) for a covariate.
+     * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
+     * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1
+     *              bbbbbbbb
+     *              76543210     
+     *   h-1        00000101 (6-1=5)
+     *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+     *           &
+     *     1        00000001 (1)
+     *              00000001        = 1 & ((h-1) >> (k-1))
+     *          +1= 00000010 =2 
+     *
+     * h=14, k=3 => h'=h-1=13, k'=k-1=2
+     *          h'      1101 =2^3+2^2+0x2^1+2^0
+     *    >>k'            11
+     *          &   00000001
+     *            = 00000001
+     *      +1    = 00000010=2    =  codtabm(14,3)   
+     * Reverse h=6 and m=16?
+     * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1.
+     * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff)
+     * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1 
+     * decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1)
+     * V3=decodtabm(14,3,2**4)=2
+     *          h'=13   1101 =2^3+2^2+0x2^1+2^0
+     *(h-1) >> (j-1)    0011 =13 >> 2
+     *          &1 000000001
+     *           = 000000001
+     *         +1= 000000010 =2
+     *                  2211
+     *                  V1=1+1, V2=0+1, V3=1+1, V4=1+1
+     *                  V3=2
+     */
+
   /* /\* for(h=1; h <=100 ;h++){  *\/ */
   /*   /\* printf("h=%2d ", h); *\/ */
   /*    /\* for(k=1; k <=10; k++){ *\/ */
@@ -7990,8 +8198,8 @@ Please run with mle=-1 to get a correct
     }
     
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
-    fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
-    fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+    fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
+    fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
 
     /* Other stuffs, more or less useful */    
     while((c=getc(ficpar))=='#' && c!= EOF){
@@ -8054,10 +8262,10 @@ Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else
-      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);
     
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
-		 model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
+		 model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \
 		 jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
       
    /*------------ free_vector  -------------*/