--- imach/src/imach.c	2016/08/25 06:59:23	1.235
+++ imach/src/imach.c	2016/08/26 09:20:19	1.237
@@ -1,6 +1,12 @@
-/* $Id: imach.c,v 1.235 2016/08/25 06:59:23 brouard Exp $
+/* $Id: imach.c,v 1.237 2016/08/26 09:20:19 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.237  2016/08/26 09:20:19  brouard
+  Summary: to valgrind
+
+  Revision 1.236  2016/08/25 10:50:18  brouard
+  *** empty log message ***
+
   Revision 1.235  2016/08/25 06:59:23  brouard
   *** empty log message ***
 
@@ -905,12 +911,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.235 2016/08/25 06:59:23 brouard Exp $ */
+/* $Id: imach.c,v 1.237 2016/08/26 09:20:19 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.235 $ $Date: 2016/08/25 06:59:23 $"; 
+char fullversion[]="$Revision: 1.237 $ $Date: 2016/08/26 09:20:19 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1111,9 +1117,11 @@ int *TvarsQind;
 #define MAXRESULTLINES 10
 int nresult=0;
 int TKresult[MAXRESULTLINES];
-double Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
 int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */
 double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
+double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
 int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */
 
 /* int *TDvar; /\**< TDvar[1]=4,  TDvarF[2]=3, TDvar[3]=6  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
@@ -2407,7 +2415,7 @@ void powell(double p[], double **xi, int
       cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
       /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
     }
-    for (k=1; k<=cptcovage;k++){
+    for (k=1; k<=cptcovage;k++){  /* For product with age */
       if(Dummy[Tvar[Tage[k]]]){
 	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
       } else{
@@ -2415,9 +2423,21 @@ void powell(double p[], double **xi, int
       }
       /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
     }
-    for (k=1; k<=cptcovprod;k++){ /*  */
+    for (k=1; k<=cptcovprod;k++){ /* For product without age */
       /* printf("prevalim Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+      if(Dummy[Tvard[k][1]==0]){
+	if(Dummy[Tvard[k][2]==0]){
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
+	}
+      }else{
+	if(Dummy[Tvard[k][2]==0]){
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
+	}
+      }
     }
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -4084,7 +4104,7 @@ void  freqsummary(char fileres[], int ia
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
 	    fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
-  fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
+  fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);
     
   strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
@@ -5983,11 +6003,13 @@ void printinghtml(char fileresu[], char
 		  int popforecast, int prevfcast, int backcast, int estepm , \
 		  double jprev1, double mprev1,double anprev1, double dateprev1, \
 		  double jprev2, double mprev2,double anprev2, double dateprev2){
-  int jj1, k1, i1, cpt;
+  int jj1, k1, i1, cpt, k4, nres;
 
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
    <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
 </ul>");
+   fprintf(fichtm,"<ul><li> model=1+age+%s\n \
+</ul>", model);
    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
    fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
 	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
@@ -6022,16 +6044,27 @@ void printinghtml(char fileresu[], char
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
    jj1=0;
+
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){
+     if(TKresult[nres]!= k1)
+       continue;
 
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
-	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
+	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+	 printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout);
+	 /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+	 /* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
        }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);fflush(stdout);
+      }
+       
        /* if(nqfveff+nqtveff 0) */ /* Test to be done */
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        if(invalidvarcomb[k1]){
@@ -6143,13 +6176,22 @@ See page 'Matrix of variance-covariance
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
    jj1=0;
+
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){
+     if(TKresult[nres]!= k1)
+       continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++)  /**< cptcoveff number of variables */
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);
+	 /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
+
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
 
        if(invalidvarcomb[k1]){
@@ -6179,7 +6221,7 @@ void printinggnuplot(char fileresu[], ch
 
   char dirfileres[132],optfileres[132];
   char gplotcondition[132];
-  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0,l=0;
+  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
@@ -6229,8 +6271,8 @@ void printinggnuplot(char fileresu[], ch
   strcpy(optfileres,"vpl");
   /* 1eme*/
   for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
-    for(nres=1; nres <= nresult; nres++) /* For each resultline */
-    for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */
+    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;
@@ -6314,22 +6356,22 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 
   
   /*2 eme*/
-  for (k1=1; k1<= m ; k1 ++) { 
+  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(nres=1; nres <= nresult; nres++) /* For each resultline */
     for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-      if(TKresult[nres]!= k)
-	continue;
       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]);
-      }
+    }
+    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]){
@@ -6372,25 +6414,25 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	
 	
   /*3eme*/
-  for (k1=1; k1<= m ; k1 ++) { 
+  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:  cov=%d state=%d",k1, cpt);
-      for(nres=1; nres <= nresult; nres++) /* For each resultline */
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	if(TKresult[nres]!= k)
-	  continue;
+      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 */
 	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]);
-	}	
       }
+      /* for(k=1; k <= ncovds; k++){ */
+      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); 
@@ -6421,7 +6463,10 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   
   /* 4eme */
   /* 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 (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);
@@ -6433,6 +6478,9 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	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); 
@@ -6463,7 +6511,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	
 /* 5eme */
   /* 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 (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);
@@ -6475,6 +6526,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	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); 
@@ -6513,7 +6567,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   
 /* 6eme */
   /* CV preval 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 (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#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
@@ -6525,6 +6582,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	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); 
@@ -6556,7 +6616,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 /* 7eme */
   if(backcast == 1){
     /* CV back preval 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 (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#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
@@ -6567,6 +6630,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	  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); 
@@ -6603,7 +6669,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   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 (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#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  */
@@ -6614,6 +6683,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	  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); 
@@ -6725,10 +6797,19 @@ 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,"#model=%s \n",model);
     fprintf(ficgp,"# ng=%d\n",ng);
-    fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);
-    for(jk=1; jk <=m; jk++) {
-      fprintf(ficgp,"#    jk=%d\n",jk);
+    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 */
+      if(TKresult[nres]!= jk)
+	continue;
+      fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
+      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");
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
       fprintf(ficgp,"\nset ter svg size 640, 480 ");
       if (ng==1){
@@ -6769,18 +6850,49 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	      break;
 	    }
 	    ij=1;/* To be checked else nbcode[0][0] wrong */
-	    for(j=3; j <=ncovmodel-nagesqr; j++) {
+	    ijp=1; /* product no age */
+	    /* for(j=3; j <=ncovmodel-nagesqr; j++) { */
+	    for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
 	      /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
-	      if(ij <=cptcovage) { /* Bug valgrind */
-		if((j-2)==Tage[ij]) { /* Bug valgrind */
-		  fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
-		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+	      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){
+		    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 */
+		    /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		  }
 		  ij++;
 		}
-	      }
-	      else
-		fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */
-	    }
+	      }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 */
+		      /* 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 */
+		      /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+		      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){
+		      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]]);
+		    }
+		  }
+  		}
+	      } else{  /* simple covariate */
+		/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
+		if(Dummy[j]==0){
+		  fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
+		}else{ /* quantitative */
+		  fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* */
+		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		}
+	      } /* end simple */
+	    } /* end j */
 	  }else{
 	    i=i-ncovmodel;
 	    if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
@@ -6798,8 +6910,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	       
 	      ij=1;
 	      for(j=3; j <=ncovmodel-nagesqr; j++){
-		if(ij <=cptcovage) { /* Bug valgrind */
-		  if((j-2)==Tage[ij]) { /* Bug valgrind */
+		if((j-2)==Tage[ij]) { /* Bug valgrind */
+		  if(ij <=cptcovage) { /* Bug valgrind */
 		    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 		    /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
 		    ij++;
@@ -8003,12 +8115,12 @@ int decoderesult ( char resultline[], in
       strcpy(resultsav,stra); /* and analyzes it */
   }
   /* Checking for missing or useless values in comparison of current model needs */
-  for(k1=1; k1<= cptcovt ;k1++){ /* model line */
-    if(Typevar[k1]==0){
+  for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+    if(Typevar[k1]==0){ /* Single covariate in model */
       match=0;
-      for(k2=1; k2 <=j;k2++){
-	if(Tvar[k1]==Tvarsel[k2]) {
-	  modelresult[k2]=k1;
+      for(k2=1; k2 <=j;k2++){/* result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+	if(Tvar[k1]==Tvarsel[k2]) {/* Tvar[1]=5 == Tvarsel[2]=5   */
+	  modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
 	  match=1;
 	  break;
 	}
@@ -8023,7 +8135,7 @@ int decoderesult ( char resultline[], in
     match=0;
     for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       if(Typevar[k1]==0){ /* Single */
-	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=5 == Tvarsel[1]=4   */
+	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
 	  resultmodel[k1]=k2;  /* resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
 	  ++match;
 	}
@@ -8050,20 +8162,26 @@ int decoderesult ( char resultline[], in
   /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */
   /*    7 0 1 1 */
   /*    8 1 1 1 */
+  /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */
+  /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */
+  /* V5*age V5 known which value for nres?  */
+  /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */
   for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */
     if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */
-      k3= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
+      k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */
       k2=(int)Tvarsel[k3]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
       k+=Tvalsel[k3]*pow(2,k4);  /*  Tvalsel[1]=1  */
-      Tresult[nres][k4+1]=Tvalsel[k3];
-      Tvresult[nres][k4+1]=(int)Tvarsel[k3];
+      Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1)  Tresult[nres][2]=0(V3=0) */
+      Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+      Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
       printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4);
       k4++;;
     }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */
       k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
-      Tqresult[nres][k4q+1]=Tvalsel[k3q];
-      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q];
+      Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
+      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
+      Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
       k4q++;;
     }