--- imach/src/imach.c	2022/08/21 09:06:25	1.332
+++ imach/src/imach.c	2022/08/31 08:23:16	1.335
@@ -1,6 +1,31 @@
-/* $Id: imach.c,v 1.332 2022/08/21 09:06:25 brouard Exp $
+/* $Id: imach.c,v 1.335 2022/08/31 08:23:16 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.335  2022/08/31 08:23:16  brouard
+  Summary: improvements...
+
+  Revision 1.334  2022/08/25 09:08:41  brouard
+  Summary: In progress for quantitative
+
+  Revision 1.333  2022/08/21 09:10:30  brouard
+  * src/imach.c (Module): Version 0.99r33 A lot of changes in
+  reassigning covariates: my first idea was that people will always
+  use the first covariate V1 into the model but in fact they are
+  producing data with many covariates and can use an equation model
+  with some of the covariate; it means that in a model V2+V3 instead
+  of codtabm(k,Tvaraff[j]) which calculates for combination k, for
+  three covariates (V1, V2, V3) the value of Tvaraff[j], but in fact
+  the equation model is restricted to two variables only (V2, V3)
+  and the combination for V2 should be codtabm(k,1) instead of
+  (codtabm(k,2), and the code should be
+  codtabm(k,TnsdVar[Tvaraff[j]]. Many many changes have been
+  made. All of these should be simplified once a day like we did in
+  hpxij() for example by using precov[nres] which is computed in
+  decoderesult for each nres of each resultline. Loop should be done
+  on the equation model globally by distinguishing only product with
+  age (which are changing with age) and no more on type of
+  covariates, single dummies, single covariates.
+
   Revision 1.332  2022/08/21 09:06:25  brouard
   Summary: Version 0.99r33
 
@@ -1259,25 +1284,25 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.332 2022/08/21 09:06:25 brouard Exp $ */
+/* $Id: imach.c,v 1.335 2022/08/31 08:23:16 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="August 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.332 $ $Date: 2022/08/21 09:06:25 $"; 
+char fullversion[]="$Revision: 1.335 $ $Date: 2022/08/31 08:23:16 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */
 /* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */
 /* Model(2)  V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */
-int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age */
+int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age but including products */
 int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
-int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 (dummy or quantit or time varying) */
-int cptcovsnq=0; /**< cptcovsnq number of simple covariates in the model but non quantitative V2+V1 =2 */
+int cptcovs=0; /**< cptcovs number of SIMPLE covariates in the model V2+V1 =2 (dummy or quantit or time varying) */
+int cptcovsnq=0; /**< cptcovsnq number of SIMPLE covariates in the model but non quantitative V2+V1 =2 */
 int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
 int cptcovprodnoage=0; /**< Number of covariate products without age */   
-int cptcoveff=0; /* Total number of covariates to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
+int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
 int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
 int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
 int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
@@ -1288,6 +1313,7 @@ int nqfveff=0; /**< nqfveff Number of Qu
 int ntveff=0; /**< ntveff number of effective time varying variables */
 int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */
 int cptcov=0; /* Working variable */
+int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs+1 declared globally ;*/
 int nobs=10;  /* Number of observations in the data lastobs-firstobs */
 int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
 int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */
@@ -1505,12 +1531,13 @@ int *TvarsQind;
 #define MAXRESULTLINESPONE 10+1
 int nresult=0;
 int parameterline=0; /* # of the parameter (type) line */
-int TKresult[MAXRESULTLINESPONE];
-int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model correspond to the k3 position in the resultline */
-int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */
+int TKresult[MAXRESULTLINESPONE]; /* TKresult[nres]=k for each resultline nres give the corresponding combination of dummies */
+int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model corresponds to the k3 position in the resultline */
+int modelresult[MAXRESULTLINESPONE][NCOVMAX];/* modelresult[k3]=k1: k1th position in the model corresponds to the k3 position in the resultline */
+int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline */
 int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line  */
 double TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */
-int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For dummy variable , variable # (output) */
+int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline */
 double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */
 double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */
 int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline */
@@ -4148,7 +4175,7 @@ double func( double *x)
 double funcone( double *x)
 {
   /* Same as func but slower because of a lot of printf and if */
-  int i, ii, j, k, mi, d, kk;
+  int i, ii, j, k, mi, d, kk, kf=0;
   int ioffset=0;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
@@ -4176,8 +4203,8 @@ double funcone( double *x)
     /* Fixed */
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
     /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
-    for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
-      cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
+    for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+      cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
 /*    cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i];  */
 /*    cov[2+6]=covar[Tvar[6]][i];  */
 /*    cov[2+6]=covar[2][i]; V2  */
@@ -4277,27 +4304,33 @@ double funcone( double *x)
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-      /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+      /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
       if(globpr){
 	fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \
 		num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
 		2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
+ /* 	printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
+ /* %11.6f %11.6f %11.6f ", \ */
+ /* 		num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
+ /* 		2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
 	for(k=1,llt=0.,l=0.; k<=nlstate; k++){
 	  llt +=ll[k]*gipmx/gsw;
 	  fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
+	  /* printf(" %10.6f",-ll[k]*gipmx/gsw); */
 	}
 	fprintf(ficresilk," %10.6f\n", -llt);
+	/* printf(" %10.6f\n", -llt); */
       }
-	} /* end of wave */
-} /* end of individual */
-for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
+    } /* end of wave */
+  } /* end of individual */
+  for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
 /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
-l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
-if(globpr==0){ /* First time we count the contributions and weights */
-	gipmx=ipmx;
-	gsw=sw;
-}
+  l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+  if(globpr==0){ /* First time we count the contributions and weights */
+    gipmx=ipmx;
+    gsw=sw;
+  }
 return -l;
 }
 
@@ -4956,7 +4989,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
   j1=0;
   
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
-  j=cptcoveff;  /* Only dummy covariates of the model */
+  j=cptcoveff;  /* Only simple dummy covariates used in the model */
   /* j=cptcovn;  /\* Only dummy covariates of the model *\/ */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
@@ -4977,7 +5010,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
 
   /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */
   /* Loop on nj=1 or 2 if dummy covariates j!=0
-   *   Loop on j1(1 to 2**cptcovn) covariate combination
+   *   Loop on j1(1 to 2**cptcoveff) covariate combination
    *     freq[s1][s2][iage] =0.
    *     Loop on iind
    *       ++freq[s1][s2][iage] weighted
@@ -5002,7 +5035,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
     if(nj==1)
       j=0;  /* First pass for the constant */
     else{
-      j=cptcovs; /* Other passes for the covariate values */
+      j=cptcoveff; /* Other passes for the covariate values number of simple covariates in the model V2+V1 =2 (simple dummy fixed or time varying) */
     }
     first=1;
     for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all dummy covariates combination of the model, ie excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
@@ -5038,13 +5071,16 @@ Title=%s 
Datafile=%s Firstpass=%d La
 	bool=1;
 	if(j !=0){
 	  if(anyvaryingduminmodel==0){ /* If All fixed covariates */
-	    if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
-	      for (z1=1; z1<=cptcovn; z1++) { /* loops on covariates in the model */
+	    if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+	      for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
 		/* if(Tvaraff[z1] ==-20){ */
 		/* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
 		/* }else  if(Tvaraff[z1] ==-10){ */
 		/* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
 		/* }else  */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/
+		/* if( iind >=imx-3) printf("Searching error iind=%d Tvaraff[z1]=%d covar[Tvaraff[z1]][iind]=%.f TnsdVar[Tvaraff[z1]]=%d, cptcoveff=%d, cptcovs=%d \n",iind, Tvaraff[z1], covar[Tvaraff[z1]][iind],TnsdVar[Tvaraff[z1]],cptcoveff, cptcovs); */
+		if(Tvaraff[z1]<1 || Tvaraff[z1]>=NCOVMAX)
+		  printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
 		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */
 		  /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */
 		  bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
@@ -5054,7 +5090,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
 		  /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
 		} /* Onlyf fixed */
 	      } /* end z1 */
-	    } /* cptcovn > 0 */
+	    } /* cptcoveff > 0 */
 	  } /* end any */
 	}/* end j==0 */
 	if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
@@ -5063,7 +5099,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
 	    m=mw[mi][iind];
 	    if(j!=0){
 	      if(anyvaryingduminmodel==1){ /* Some are varying covariates */
-		for (z1=1; z1<=cptcovn; z1++) {
+		for (z1=1; z1<=cptcoveff; z1++) {
 		  if( Fixed[Tmodelind[z1]]==1){
 		    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
 		    if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's 
@@ -5071,7 +5107,12 @@ Title=%s 
Datafile=%s Firstpass=%d La
 										      constant and age model which counts them. */
 		      bool=0; /* not selected */
 		  }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
-		    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) {
+		    /* i1=Tvaraff[z1]; */
+		    /* i2=TnsdVar[i1]; */
+		    /* i3=nbcode[i1][i2]; */
+		    /* i4=covar[i1][iind]; */
+		    /* if(i4 != i3){ */
+		    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { /* Bug valgrind */
 		      bool=0;
 		    }
 		  }
@@ -5132,9 +5173,9 @@ Title=%s 
Datafile=%s Firstpass=%d La
       
       
       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
-      if(cptcovn==0 && nj==1) /* no covariate and first pass */
+      if(cptcoveff==0 && nj==1) /* no covariate and first pass */
         pstamp(ficresp);
-      if  (cptcovn>0 && j!=0){
+      if  (cptcoveff>0 && j!=0){
         pstamp(ficresp);
 	printf( "\n#********** Variable "); 
 	fprintf(ficresp, "\n#********** Variable "); 
@@ -5187,14 +5228,17 @@ Title=%s 
Datafile=%s Firstpass=%d La
       /* } */
 
       fprintf(ficresphtm,"
");
-      if((cptcovn==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
+      if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
         fprintf(ficresp, " Age");
-      if(nj==2) for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+      if(nj==2) for (z1=1; z1<=cptcoveff; z1++) {
+	  printf(" V%d=%d, z1=%d, Tvaraff[z1]=%d, j1=%d, TnsdVar[Tvaraff[%d]]=%d |",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])], z1, Tvaraff[z1], j1,z1,TnsdVar[Tvaraff[z1]]);
+	  fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+	}
       for(i=1; i<=nlstate;i++) {
-	if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d)  N(%d)  N  ",i,i);
+	if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d)  N(%d)  N  ",i,i);
 	fprintf(ficresphtm, "| Age | Prev(%d) | N(%d) | N",i,i);
       }
-      if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");
+      if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");
       fprintf(ficresphtm, "\n");
       
       /* Header of frequency table by age */
@@ -5262,14 +5306,14 @@ Title=%s | 
Datafile=%s Firstpass=%d La
 	}
 	
 	/* Writing ficresp */
-	if(cptcovn==0 && nj==1){ /* no covariate and first pass */
+	if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
           if( iage <= iagemax){
 	    fprintf(ficresp," %d",iage);
           }
         }else if( nj==2){
           if( iage <= iagemax){
 	    fprintf(ficresp," %d",iage);
-            for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+            for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
           }
 	}
 	for(s1=1; s1 <=nlstate ; s1++){
@@ -5284,7 +5328,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
 	  }
 	  if( iage <= iagemax){
 	    if(pos>=1.e-5){
-	      if(cptcovn==0 && nj==1){ /* no covariate and first pass */
+	      if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
 	        fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
               }else if( nj==2){
 	        fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
@@ -5293,7 +5337,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
 	      /*probs[iage][s1][j1]= pp[s1]/pos;*/
 	      /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/
 	    } else{
-	      if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
+	      if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
 	      fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f",iage, prop[s1][iage],pospropta);
 	    }
 	  }
@@ -5319,7 +5363,7 @@ Title=%s | 
Datafile=%s Firstpass=%d La
 	}
 	fprintf(ficresphtmfr,"\n ");
 	fprintf(ficresphtm,"\n");
-	if((cptcovn==0 && nj==1)|| nj==2 ) {
+	if((cptcoveff==0 && nj==1)|| nj==2 ) {
 	  if(iage <= iagemax)
  	    fprintf(ficresp,"\n");
 	}
@@ -5577,7 +5621,7 @@ void prevalence(double ***probs, double
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
   first=0;
-  for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
+  for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of simple dummy covariates */
     for (i=1; i<=nlstate; i++)  
       for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
 	prop[i][iage]=0.0;
@@ -5901,7 +5945,7 @@ void  concatwav(int wav[], int **dh, int
        nbcode[k][j]=0; /* Valgrind */
 
    /* Loop on covariates without age and products and no quantitative variable */
-   for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
+   for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
      for (j=-1; (j < maxncov); j++) Ndum[j]=0;
      if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
        switch(Fixed[k]) {
@@ -5999,8 +6043,12 @@ void  concatwav(int wav[], int **dh, int
 	 break;
        } /* end switch */
      } /* end dummy test */
-     if(Dummy[k]==1 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
+     if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */ 
        for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
+	 if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){
+	   printf("Error k=%d \n",k);
+	   exit(1);
+	 }
 	 if(isnan(covar[Tvar[k]][i])){
 	   printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
 	   fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
@@ -6008,7 +6056,7 @@ void  concatwav(int wav[], int **dh, int
 	   exit(1);
          }
        }
-     }
+     } /* end Quanti */
    } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/  
   
    for (k=-1; k< maxncov; k++) Ndum[k]=0; 
@@ -6022,13 +6070,22 @@ void  concatwav(int wav[], int **dh, int
   
    ij=0;
    /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
-   for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+   for (k=1; k<=  cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
+     /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
      /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
      /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */
-     if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */
+     if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy simple and non empty in the model */
+       /* Typevar[k] =0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */
+       /* Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product*/
        /* If product not in single variable we don't print results */
        /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
-       ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */
+       ++ij;/*    V5 + V4 + V3 + V4*V3 + V5*age + V2 +  V1*V2 + V1*age + V1, *//* V5 quanti, V2 quanti, V4, V3, V1 dummies */
+       /* k=       1    2   3     4       5       6      7       8        9  */
+       /* Tvar[k]= 5    4    3    6       5       2      7       1        1  */
+       /* ij            1    2                                            3  */  
+       /* Tvaraff[ij]=  4    3                                            1  */
+       /* Tmodelind[ij]=2    3                                            9  */
+       /* TmodelInvind[ij]=2 1                                            1  */
        Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
        Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
        TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
@@ -6044,7 +6101,7 @@ void  concatwav(int wav[], int **dh, int
    } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
    /* ij--; */
    /* cptcoveff=ij; /\*Number of total covariates*\/ */
-   *cptcov=ij; /* cptcov= Number of total real effective covariates: effective (used as cptcoveff in other functions)
+   *cptcov=ij; /* cptcov= Number of total real effective simple dummies (fixed or time  arying) effective (used as cptcoveff in other functions)
 		* because they can be excluded from the model and real
 		* if in the model but excluded because missing values, but how to get k from ij?*/
    for(j=ij+1; j<= cptcovt; j++){
@@ -6445,11 +6502,11 @@ void  concatwav(int wav[], int **dh, int
    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 */
+   for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ /* To be done*/
      fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
    }
    for(j=1;j<=cptcoveff;j++) 
-     fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]);
+     fprintf(ficresprobmorprev," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]);
    fprintf(ficresprobmorprev,"\n");
 
    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
@@ -7080,7 +7137,7 @@ To be simple, these graphs help to under
 
    for(nres=1;nres <=nresult; nres++){ /* For each resultline */
    for(j1=1; j1<=tj;j1++){ /* For any combination of dummy covariates, fixed and varying */
-     printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d \n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj);
+     printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs);
      if(tj != 1 && TKresult[nres]!= j1)
        continue;
 
@@ -7088,25 +7145,63 @@ To be simple, these graphs help to under
      /* for(nres=1;nres <=1; nres++){ /\* For each resultline *\/ */
      /* /\* for(nres=1;nres <=nresult; nres++){ /\\* For each resultline *\\/ *\/ */
      if  (cptcovn>0) {
-       fprintf(ficresprob, "\n#********** Variable "); 
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
-       fprintf(ficresprob, "**********\n#\n");
+       fprintf(ficresprob, "\n#********** Variable ");
        fprintf(ficresprobcov, "\n#********** Variable "); 
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+       fprintf(ficgp, "\n#********** Variable ");
+       fprintf(fichtmcov, "\n
********** Variable "); 
+       fprintf(ficresprobcor, "\n#********** Variable ");    
+
+       /* Including quantitative variables of the resultline to be done */
+       for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */
+	 printf("Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
+	 fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
+	 /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */
+	 if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline  */
+	   if(Fixed[modelresult[nres][z1]]==0){ /* Fixed referenced to model equation */
+	     fprintf(ficresprob,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	     fprintf(ficresprobcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	     fprintf(ficgp,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	     fprintf(fichtmcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	     fprintf(ficresprobcor,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	     fprintf(ficresprob,"fixed ");
+	     fprintf(ficresprobcov,"fixed ");
+	     fprintf(ficgp,"fixed ");
+	     fprintf(fichtmcov,"fixed ");
+	     fprintf(ficresprobcor,"fixed ");
+	   }else{
+	     fprintf(ficresprob,"varyi ");
+	     fprintf(ficresprobcov,"varyi ");
+	     fprintf(ficgp,"varyi ");
+	     fprintf(fichtmcov,"varyi ");
+	     fprintf(ficresprobcor,"varyi ");
+	   }
+	 }else if(Dummy[modelresult[nres][z1]]==1){ /* Quanti variable */
+	   /* For each selected (single) quantitative value */
+	   fprintf(ficresprob," V%d=%f ",Tvqresult[nres][z1],Tqresult[nres][z1]);
+	   if(Fixed[modelresult[nres][z1]]==0){ /* Fixed */
+	     fprintf(ficresprob,"fixed ");
+	     fprintf(ficresprobcov,"fixed ");
+	     fprintf(ficgp,"fixed ");
+	     fprintf(fichtmcov,"fixed ");
+	     fprintf(ficresprobcor,"fixed ");
+	   }else{
+	     fprintf(ficresprob,"varyi ");
+	     fprintf(ficresprobcov,"varyi ");
+	     fprintf(ficgp,"varyi ");
+	     fprintf(fichtmcov,"varyi ");
+	     fprintf(ficresprobcor,"varyi ");
+	   }
+	 }else{
+	   printf("Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff);  /* end if dummy  or quanti */
+	   fprintf(ficlog,"Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff);  /* end if dummy  or quanti */
+	   exit(1);
+	 }
+       } /* End loop on variable of this resultline */
+       /* for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); */
+       fprintf(ficresprob, "**********\n#\n");
        fprintf(ficresprobcov, "**********\n#\n");
-			
-       fprintf(ficgp, "\n#********** Variable "); 
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
        fprintf(ficgp, "**********\n#\n");
-			
-			
-       fprintf(fichtmcov, "\n
********** Variable "); 
-       /* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
        fprintf(fichtmcov, "**********\n
");
-			
-       fprintf(ficresprobcor, "\n#********** Variable ");    
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
        fprintf(ficresprobcor, "**********\n#");    
        if(invalidvarcomb[j1]){
 	 fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); 
@@ -7118,57 +7213,66 @@ To be simple, these graphs help to under
      trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
      gp=vector(1,(nlstate)*(nlstate+ndeath));
      gm=vector(1,(nlstate)*(nlstate+ndeath));
-     for (age=bage; age<=fage; age ++){ 
+     for (age=bage; age<=fage; age ++){ /* Fo each age we feed the model equation with covariates, using precov as in hpxij() ? */
        cov[2]=age;
        if(nagesqr==1)
 	 cov[3]= age*age;
-       /* for (k=1; k<=cptcovn;k++) { */
-       /* 	 cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; */
-       for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
-	 /* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */
-	 cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])];
-	 /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
-								    * 1  1 1 1 1
-								    * 2  2 1 1 1
-								    * 3  1 2 1 1
-								    */
-	 /* nbcode[1][1]=0 nbcode[1][2]=1;*/
-       }
-       /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
-       /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */
-       /*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-       for (k=1; k<=cptcovage;k++){  /* For product with age */
-	 if(Dummy[Tage[k]]==2){ /* dummy with age */
-	   cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2];
-	   /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */
-	 } else if(Dummy[Tage[k]]==3){ /* quantitative with age */
-	   printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]);
-	   /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */
-	   /* exit(1); */
-	   /* cov[++k1]=Tqresult[nres][k];  */
-	 }
-	 /* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */
-       }
-       for (k=1; k<=cptcovprod;k++){/* For product without age */
-	 if(Dummy[Tvard[k][1]]==0){
-	   if(Dummy[Tvard[k][2]]==0){
-	     cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])];
-	     /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-	   }else{ /* Should we use the mean of the quantitative variables? */
-	     cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[nres][k]];
-	     /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
-	   }
+       /* New code end of combination but for each resultline */
+       for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */ 
+	 if(Typevar[k1]==1){ /* A product with age */
+	   cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
 	 }else{
-	   if(Dummy[Tvard[k][2]]==0){
-	     cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[Tvard[k][1]]];
-	     /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
-	   }else{
-	     cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][TnsdVar[Tvard[k][1]]]*  Tqinvresult[nres][TnsdVar[Tvard[k][2]]];
-	     /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */
-	   }
+	   cov[2+nagesqr+k1]=precov[nres][k1];
 	 }
-	 /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-       }			
+       }/* End of loop on model equation */
+/* Old code */
+       /* /\* for (k=1; k<=cptcovn;k++) { *\/ */
+       /* /\* 	 cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; *\/ */
+       /* for (k=1; k<=nsd;k++) { /\* For single dummy covariates only *\/ */
+       /* 	 /\* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates *\/ */
+       /* 	 cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])]; */
+       /* 	 /\*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*\//\* j1 1 2 3 4 */
+       /* 								    * 1  1 1 1 1 */
+       /* 								    * 2  2 1 1 1 */
+       /* 								    * 3  1 2 1 1 */
+       /* 								    *\/ */
+       /* 	 /\* nbcode[1][1]=0 nbcode[1][2]=1;*\/ */
+       /* } */
+       /* /\* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 *\/ */
+       /* /\* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] *\/ */
+       /* /\*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; *\/ */
+       /* for (k=1; k<=cptcovage;k++){  /\* For product with age *\/ */
+       /* 	 if(Dummy[Tage[k]]==2){ /\* dummy with age *\/ */
+       /* 	   cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2]; */
+       /* 	   /\* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */
+       /* 	 } else if(Dummy[Tage[k]]==3){ /\* quantitative with age *\/ */
+       /* 	   printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); */
+       /* 	   /\* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\\* Using the mean of quantitative variable Tvar[Tage[k]] /\\* Tqresult[nres][k]; *\\/ *\/ */
+       /* 	   /\* exit(1); *\/ */
+       /* 	   /\* cov[++k1]=Tqresult[nres][k];  *\/ */
+       /* 	 } */
+       /* 	 /\* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */
+       /* } */
+       /* for (k=1; k<=cptcovprod;k++){/\* For product without age *\/ */
+       /* 	 if(Dummy[Tvard[k][1]]==0){ */
+       /* 	   if(Dummy[Tvard[k][2]]==0){ */
+       /* 	     cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])]; */
+       /* 	     /\* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
+       /* 	   }else{ /\* Should we use the mean of the quantitative variables? *\/ */
+       /* 	     cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[nres][k]]; */
+       /* 	     /\* cov[++k1]=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(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[Tvard[k][1]]]; */
+       /* 	     /\* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; *\/ */
+       /* 	   }else{ */
+       /* 	     cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][TnsdVar[Tvard[k][1]]]*  Tqinvresult[nres][TnsdVar[Tvard[k][2]]]; */
+       /* 	     /\* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; *\/ */
+       /* 	   } */
+       /* 	 } */
+       /* 	 /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
+       /* } */			
 /* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/			
        for(theta=1; theta <=npar; theta++){
 	 for(i=1; i<=npar; i++)
@@ -10246,7 +10350,7 @@ int decoderesult( char resultline[], int
   int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0;
   char resultsav[MAXLINE];
   /* int resultmodel[MAXLINE]; */
-  int modelresult[MAXLINE];
+  /* int modelresult[MAXLINE]; */
   char stra[80], strb[80], strc[80], strd[80],stre[80];
 
   removefirstspace(&resultline);
@@ -10255,18 +10359,18 @@ int decoderesult( char resultline[], int
   strcpy(resultsav,resultline);
   printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline);
   if (strlen(resultsav) >1){
-    j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
+    j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
   }
   if(j == 0){ /* Resultline but no = */
     TKresult[nres]=0; /* Combination for the nresult and the model */
     return (0);
   }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
-    printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model);
-    fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model);
+    printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
+    fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
     /* return 1;*/
   }
-  for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */
+  for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
     if(nbocc(resultsav,'=') >1){
       cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//*     resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */
       /* If resultsav= "V4= 1 V5=25.1 V3=0" with a blank then strb="V4=" and stra="1 V5=25.1 V3=0" */
@@ -10290,13 +10394,13 @@ int decoderesult( char resultline[], int
   }
   /* Checking for missing or useless values in comparison of current model needs */
   /* Feeds resultmodel[nres][k1]=k2 for k1th product covariate with age in the model equation fed by the index k2 of the resutline*/
-  for(k1=1; k1<= cptcovt ;k1++){ /* Loop on model. model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+  for(k1=1; k1<= cptcovt ;k1++){ /* Loop on MODEL LINE V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
     if(Typevar[k1]==0){ /* Single covariate in model */
       /* 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
       match=0;
       for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
 	if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
-	  modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
+	  modelresult[nres][k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
 	  match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
 	  break;
 	}
@@ -10311,7 +10415,7 @@ int decoderesult( char resultline[], int
       match=0;
       for(k2=1; k2 <=j;k2++){/* Loop on resultline.  jth occurence of = signs in the result line. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
 	if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
-	  modelresult[k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
+	  modelresult[nres][k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
 	  resultmodel[nres][k1]=k2; /* Added here */
 	  printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]);
 	  match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
@@ -10320,7 +10424,7 @@ int decoderesult( char resultline[], int
       }
       if(match == 0){
 	printf("Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
-	fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=%s\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
+	fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
       return 1;
       }
     }else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/
@@ -10336,7 +10440,7 @@ int decoderesult( char resultline[], int
       }
       if(match == 0){
 	printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
-	fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",k1,Tvardk[k1][1], resultline, model);
+	fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
 	return 1;
       }
       match=0;
@@ -10350,19 +10454,21 @@ int decoderesult( char resultline[], int
       }
       if(match == 0){
 	printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][2], resultline, model);
-	fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=%s\n",k1,Tvardk[k1][2], resultline, model);
+	fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=%s\n",Tvardk[k1][2], resultline, model);
 	return 1;
       }
     }/* End of testing */
-  }/* End loop cptcovt
+  }/* End loop cptcovt */
   /* Checking for missing or useless values in comparison of current model needs */
   /* Feeds resultmodel[nres][k1]=k2 for single covariate (k1) in the model equation */
-  for(k2=1; k2 <=j;k2++){ /* Loop on resultline variables: result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+  for(k2=1; k2 <=j;k2++){ /* j or cptcovs is the number of single covariates used either in the model line as well as in the result line (dummy or quantitative)
+			   * Loop on resultline variables: result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
     match=0;
     for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       if(Typevar[k1]==0){ /* Single only */
 	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
 	  resultmodel[nres][k1]=k2;  /* k1th position in the model equation corresponds to k2th position in the result line. resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
+	  modelresult[nres][k2]=k1; /* k1th position in the model equation corresponds to k2th position in the result line. modelresult[1]=2 modelresult[2]=1  modelresult[3]=3  remodelresult[4]=6 modelresult[5]=9 */
 	  ++match;
 	}
       }
@@ -10377,7 +10483,7 @@ int decoderesult( char resultline[], int
       return 1;
     }
   }
-      
+  /* cptcovres=j /\* Number of variables in the resultline is equal to cptcovs and thus useless *\/     */
   /* We need to deduce which combination number is chosen and save quantitative values */
   /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
   /* nres=1st result line: V4=1 V5=25.1 V3=0  V2=8 V1=1 */
@@ -10396,12 +10502,13 @@ int decoderesult( char resultline[], int
   /* 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++){ /* loop k1 on position in the model line (excluding product) */
+  for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* cptcovt number of covariates (excluding 1 and age or age*age) in the MODEL equation.
+						   * loop on position k1 in the MODEL LINE */
     /* k counting number of combination of single dummies in the equation model */
     /* k4 counting single dummies in the equation model */
     /* k4q counting single quantitatives in the equation model */
-    if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single */
-       /* k4+1= position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */
+    if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, k1 is sorting according to MODEL, but k3 to resultline */
+       /* k4+1= (not always if quant in model) position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */
       /* modelresult[k3]=k1: k3th position in the result line corresponds to the k1 position in the model line (doesn't work with products)*/
       /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */
       /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline                        */
@@ -10409,19 +10516,21 @@ int decoderesult( char resultline[], int
       /* Tvarsel[k3]: Name of the variable at the k3th position in the result line.                                                  */
       /* Tvalsel[k3]: Value of the variable at the k3th position in the result line.                                                 */
       /* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline                   */
-      /* Tvresult[nres][result_position]= id of the dummy variable at the result_position in the nres resultline                     */
+      /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline                     */
       /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line                                        */
       /* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line                                                      */
       k3= resultmodel[nres][k1]; /* From position k1 in model get position k3 in result line */
       /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
       k2=(int)Tvarsel[k3]; /* from position k3 in resultline get name k2: nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/
       k+=Tvalsel[k3]*pow(2,k4);  /* nres=1 k1=2 Tvalsel[1]=1 (V4=1); k1=3 k3=2 Tvalsel[2]=0 (V3=0) */
-      TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Stores the value into the name of the variable. */
+      TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][Name]=Value; stores the value into the name of the variable. */
       /* Tinvresult[nres][4]=1 */
-      Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) */
-      Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+      /* Tresult[nres][k4+1]=Tvalsel[k3];/\* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) *\/ */
+      Tresult[nres][k3]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) */
+      /* Tvresult[nres][k4+1]=(int)Tvarsel[k3];/\* Tvresult[nres][1]=4 Tvresult[nres][3]=1 *\/ */
+      Tvresult[nres][k3]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
       Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
-      precov[nres][k1]=Tvalsel[k3];
+      precov[nres][k1]=Tvalsel[k3]; /* Value from resultline of the variable at the k1 position in the model */
       printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4);
       k4++;;
     }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */
@@ -10431,8 +10540,12 @@ int decoderesult( char resultline[], int
       k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 5 =k3q */
       k2q=(int)Tvarsel[k3q]; /*  Name of variable at k3q th position in the resultline */
       /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
-      Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
-      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
+      /* Tqresult[nres][k4q+1]=Tvalsel[k3q]; /\* Tqresult[nres][1]=25.1 *\/ */
+      /* Tvresult[nres][k4q+1]=(int)Tvarsel[k3q];/\* Tvresult[nres][1]=4 Tvresult[nres][3]=1 *\/ */
+      /* Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /\* Tvqresult[nres][1]=5 *\/ */
+      Tqresult[nres][k3q]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
+      Tvresult[nres][k3q]=(int)Tvarsel[k3q];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+      Tvqresult[nres][k3q]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
       Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       precov[nres][k1]=Tvalsel[k3q];
@@ -10444,15 +10557,15 @@ int decoderesult( char resultline[], int
       
       k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
       k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/
-      TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
+      TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */
       precov[nres][k1]=Tvalsel[k3];
       printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]);
     }else if( Dummy[k1]==3 ){ /* For quant with age product */
       k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
-      TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
+      TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */
       precov[nres][k1]=Tvalsel[k3q];
-      printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
+      printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
     }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */
       precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];      
       printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]);
@@ -10462,7 +10575,7 @@ int decoderesult( char resultline[], int
     }
   }
   
-  TKresult[nres]=++k; /* Combination for the nresult and the model */
+  TKresult[nres]=++k; /* Number of combinations of dummies for the nresult and the model =Tvalsel[k3]*pow(2,k4) + 1*/
   return (0);
 }
 
@@ -10622,9 +10735,17 @@ int decodemodel( char model[], int lasto
 	    Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
 						because this model-covariate is a construction we invent a new column
 						which is after existing variables ncovcol+nqv+ntv+nqtv + k1
-						If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2
+						If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
 						thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
 						Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
+	    /* Please remark that the new variables are model dependent */
+	    /* If we have 4 variable but the model uses only 3, like in
+	     * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
+	     *  k=     1     2       3   4     5        6        7       8
+	     * Tvar[k]=1     1       2   3     2        3       (5       6) (and not 4 5 because of V4 missing)
+	     * Tage[kk]    [1]= 2           [2]=5      [3]=6                  kk=1 to cptcovage=3
+	     * Tvar[Tage[kk]][1]=2          [2]=2      [3]=3
+	     */
 	    Typevar[k]=2;  /* 2 for double fixed dummy covariates */
 	    cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
 	    Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
@@ -10727,8 +10848,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       modell[k].maintype= FTYPE;
       modell[k].subtype= FQ;
       nsq++;
-      TvarsQ[nsq]=Tvar[k];
-      TvarsQind[nsq]=k;
+      TvarsQ[nsq]=Tvar[k]; /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary see below */
+      TvarsQind[nsq]=k;    /* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */
       ncovf++;
       TvarF[ncovf]=Tvar[k];
       TvarFind[ncovf]=k;
@@ -10759,8 +10880,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       modell[k].subtype= VQ;
       ncovv++; /* Only simple time varying variables */
       nsq++;
-      TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */
-      TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate */
+      TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary here) */
+      TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate *//* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */
       TvarV[ncovv]=Tvar[k];
       TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
@@ -11703,7 +11824,7 @@ int main(int argc, char *argv[])
   
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
-  int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/
+  /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */
   int c,  h , cpt, c2;
   int jl=0;
   int i1, j1, jk, stepsize=0;
@@ -12346,6 +12467,7 @@ Please run with mle=-1 to get a correct
   Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
   TvarsDind=ivector(1,NCOVMAX); /*  */
   TnsdVar=ivector(1,NCOVMAX); /*  */
+    /* for(i=1; i<=NCOVMAX;i++) TnsdVar[i]=3; */
   TvarsD=ivector(1,NCOVMAX); /*  */
   TvarsQind=ivector(1,NCOVMAX); /*  */
   TvarsQ=ivector(1,NCOVMAX); /*  */
@@ -12468,7 +12590,7 @@ Please run with mle=-1 to get a correct
   Ndum =ivector(-1,NCOVMAX);  
   cptcoveff=0;
   if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */
-    tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
+    tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; as well as calculate cptcoveff or number of total effective dummy covariates*/
   }
   
   ncovcombmax=pow(2,cptcoveff);
@@ -12607,11 +12729,18 @@ Title=%s 
Datafile=%s Firstpass=%d La
 	  optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
 
-  fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain 
\nSponsored by Copyright (C)  2002-2015 INED-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
  \
-
 \n\
+  fprintf(fichtm,"\n\n\
+IMaCh %s\n\
+ IMaCh for Interpolated Markov Chain 
\n\
+Sponsored by Copyright (C)  2002-2015 INED\
+-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \
+(Grant-in-Aid for Scientific Research 25293121) - \
+Intel Software 2015-2018
 \n", optionfilehtm);
+  
+  fprintf(fichtm,"
 \n\
 IMaCh-%s 
 %s \
 
 \n\
-Title=%s 
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\
+This file: %sTitle=%s 
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\
 \n\
 
\
  - Parameter files\n\
@@ -12620,7 +12749,7 @@ Title=%s
 Datafile=%s Firstpass=%d La
  - Log file of the run: %s
 \n\
  - Gnuplot file name: %s
 \n\
  - Date and time at start: %s
\n",\
-	  optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+	  version,fullversion,optionfilehtm,optionfilehtm,title,datafile,datafile,firstpass,lastpass,stepm, weightopt, model, \
 	  optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
 	  fileres,fileres,\
 	  filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
@@ -12934,6 +13063,7 @@ Please run with mle=-1 to get a correct
     globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
+          /* exit(0); */
     for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);
     printf("\n");
@@ -13612,8 +13742,8 @@ Please run with mle=-1 to get a correct
 	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
       }
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficreseij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+	printf(" V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */
+	fprintf(ficreseij,"V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]);
       }
       fprintf(ficreseij,"******\n");
       printf("******\n");
@@ -13670,23 +13800,70 @@ Please run with mle=-1 to get a correct
     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++){ /* For any combination of dummy covariates, fixed and varying */
-      if(i1 != 1 && TKresult[nres]!= k)
+    for(nres=1; nres <= nresult; nres++) /* For each resultline, find the combination and output results according to the values of dummies and then quanti.  */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying. For each nres and each value at position k
+			  * we know Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline
+			  * Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline 
+			  * and Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */
+      /* */
+      if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
 	continue;
       printf("\n# model %s \n#****** Result for:", model);
       fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
       fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
-      for(j=1;j<=cptcoveff;j++){ 
-	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-      }
-      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-      }	
+      /* It might not be a good idea to mix dummies and quantitative */
+      /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */
+      for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
+	/* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* Output by variables in the resultline *\/ */
+	/* Tvaraff[j] is the name of the dummy variable in position j in the equation model:
+	 * Tvaraff[1]@9={4, 3, 0, 0, 0, 0, 0, 0, 0}, in model=V5+V4+V3+V4*V3+V5*age
+	 * (V5 is quanti) V4 and V3 are dummies
+	 * TnsdVar[4] is the position 1 and TnsdVar[3]=2 in codtabm(k,l)(V4  V3)=V4  V3
+	 *                                                              l=1 l=2
+	 *                                                           k=1  1   1   0   0
+	 *                                                           k=2  2   1   1   0
+	 *                                                           k=3 [1] [2]  0   1
+	 *                                                           k=4  2   2   1   1
+	 * If nres=1 result: V3=1 V4=0 then k=3 and outputs
+	 * If nres=2 result: V4=1 V3=0 then k=2 and outputs
+	 * nres=1 =>k=3 j=1 V4= nbcode[4][codtabm(3,1)=1)=0; j=2  V3= nbcode[3][codtabm(3,2)=2]=1
+	 * nres=2 =>k=2 j=1 V4= nbcode[4][codtabm(2,1)=2)=1; j=2  V3= nbcode[3][codtabm(2,2)=1]=0
+	 */
+	/* Tvresult[nres][j] Name of the variable at position j in this resultline */
+	/* Tresult[nres][j] Value of this variable at position j could be a float if quantitative  */
+/* We give up with the combinations!! */
+	printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]);  /* end if dummy  or quanti */
+
+	if(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline  */
+	  printf("V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	  fprintf(ficlog,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	  fprintf(ficrest,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+	  if(Fixed[modelresult[nres][j]]==0){ /* Fixed */
+	    printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
+	  }else{
+	    printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi ");
+	  }
+	  /* fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+	  /* fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+	}else if(Dummy[modelresult[nres][j]]==1){ /* Quanti variable */
+	  /* For each selected (single) quantitative value */
+	  printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	  if(Fixed[modelresult[nres][j]]==0){ /* Fixed */
+	    printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
+	  }else{
+	    printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi ");
+	  }
+	}else{
+	  printf("Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff);  /* end if dummy  or quanti */
+	  fprintf(ficlog,"Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff);  /* end if dummy  or quanti */
+	  exit(1);
+	}
+      } /* End loop for each variable in the resultline */
+      /* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */
+      /* 	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /\* Wrong j is not in the equation model *\/ */
+      /* 	fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
+      /* 	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
+      /* }	 */
       fprintf(ficrest,"******\n");
       fprintf(ficlog,"******\n");
       printf("******\n");
@@ -13694,12 +13871,13 @@ Please run with mle=-1 to get a correct
       fprintf(ficresstdeij,"\n#****** ");
       fprintf(ficrescveij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-	fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-      }
-      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	fprintf(ficresstdeij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficrescveij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+	fprintf(ficresstdeij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]);
+	/* fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+	/* fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+      }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value, TvarsQind gives the position of a quantitative in model equation  */
+	fprintf(ficresstdeij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]);
+	fprintf(ficrescveij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]);
       }	
       fprintf(ficresstdeij,"******\n");
       fprintf(ficrescveij,"******\n");
@@ -13707,7 +13885,8 @@ Please run with mle=-1 to get a correct
       fprintf(ficresvij,"\n#****** ");
       /* pstamp(ficresvij); */
       for(j=1;j<=cptcoveff;j++) 
-	fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]);
+	fprintf(ficresvij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]);
+	/* fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]); */
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
 	/* fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); /\* To solve *\/ */
 	fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */
@@ -13741,7 +13920,7 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
 	else
 	  fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
-	fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+	fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
 	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	fprintf(ficrest,"\n");
 	/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
@@ -13788,12 +13967,12 @@ Please run with mle=-1 to get a correct
       printf("done selection\n");fflush(stdout);
       fprintf(ficlog,"done selection\n");fflush(ficlog);
       
-    } /* End k selection */
+    } /* End k selection or end covariate selection for nres */
 
     printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
 
-    /* variance-covariance of forward period prevalence*/
+    /* variance-covariance of forward period prevalence */
     varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);