--- imach/src/imach.c	2022/08/06 07:18:25	1.330
+++ imach/src/imach.c	2022/08/21 09:06:25	1.332
@@ -1,6 +1,30 @@
-/* $Id: imach.c,v 1.330 2022/08/06 07:18:25 brouard Exp $
+/* $Id: imach.c,v 1.332 2022/08/21 09:06:25 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.332  2022/08/21 09:06:25  brouard
+  Summary: Version 0.99r33
+
+  * 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.331  2022/08/07 05:40:09  brouard
+  *** empty log message ***
+
   Revision 1.330  2022/08/06 07:18:25  brouard
   Summary: last 0.99r31
 
@@ -1235,12 +1259,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.330 2022/08/06 07:18:25 brouard Exp $ */
+/* $Id: imach.c,v 1.332 2022/08/21 09:06:25 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
-char copyright[]="July 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.330 $ $Date: 2022/08/06 07:18:25 $"; 
+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 strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1403,6 +1427,7 @@ int *ncodemaxwundef;  /* ncodemax[j]= Nu
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs; /* Global pointer */
 double ***mobaverage, ***mobaverages; /* New global variable */
+double **precov; /* New global variable to store for each resultline, values of model covariates given by the resultlines (in order to speed up)  */
 double *ageexmed,*agecens;
 double dateintmean=0;
   double anprojd, mprojd, jprojd; /* For eventual projections */
@@ -1433,7 +1458,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv
        *       cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti
        *       cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1
        *       cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1
-       *       covar[k,i], value of kth fixed covariate dummy or quanti :
+       *       covar[Vk,i], value of the Vkth fixed covariate dummy or quanti for individual i:
        *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
        * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + V9 + V9*age + V10
        *   k=  1    2      3       4     5       6      7        8   9     10       11 
@@ -1483,12 +1508,12 @@ int parameterline=0; /* # of the paramet
 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 Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */
-int TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable or quanti value (output) */
+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) */
-double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */
+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]; /* For quantitative variable , variable # (output) */
+int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline */
 
 /* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1
   # States 1=Coresidence, 2 Living alone, 3 Institution
@@ -2784,7 +2809,7 @@ void powell(double p[], double **xi, int
   /*  0.51326036147820708, 0.48673963852179264} */
   /* If we start from prlim again, prlim tends to a constant matrix */
     
-  int i, ii,j,k;
+    int i, ii,j,k, k1;
   double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */
@@ -2815,48 +2840,87 @@ void powell(double p[], double **xi, int
      if(nagesqr==1){
       cov[3]= agefin*agefin;
      }
-    for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
- 			/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
-      cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])];
-      /* cov[++k1]=nbcode[TvarsD[k]][codtabm(ij,k)]; */
-      /* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
-    }
-    for (k=1; k<=nsq;k++) { /* For single varying covariates only */
- 			/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
-      cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
-      /* cov[++k1]=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 product with age */
-      if(Dummy[Tage[k]]==2){ /* dummy with age */
-	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,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 */
-	cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
-	/* cov[++k1]=Tqresult[nres][k];  */
-      }
-      /* 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 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]); */
-      if(Dummy[Tvard[k][1]]==0){
-	if(Dummy[Tvard[k][2]]==0){
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
-	  /* cov[++k1]=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,Tvard[k][1])] * Tqresult[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(ij,Tvard[k][2])] * Tqinvresult[nres][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][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
-	  /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */
-	}
-      }
-    }
+     /* Model(2)  V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */
+     /* total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age */
+     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{
+	 cov[2+nagesqr+k1]=precov[nres][k1];
+       }
+     }/* End of loop on model equation */
+     
+/* Start of old code (replaced by a loop on position in the model equation */
+    /* for (k=1; k<=nsd;k++) { /\* For single dummy covariates only of the model *\/ */
+    /* 			/\* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates *\/ */
+    /*   /\* cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])]; *\/ */
+    /*   cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TnsdVar[TvarsD[k]])]; */
+    /*   /\* model = 1 +age + V1*V3 + age*V1 + V2 + V1 + age*V2 + V3 + V3*age + V1*V2  */
+    /*    * k                  1        2      3    4      5      6     7        8 */
+    /*    *cov[]   1    2      3        4      5    6      7      8     9       10 */
+    /*    *TypeVar[k]          2        1      0    0      1      0     1        2 */
+    /*    *Dummy[k]            0        2      0    0      2      0     2        0 */
+    /*    *Tvar[k]             4        1      2    1      2      3     3        5 */
+    /*    *nsd=3                              (1)  (2)           (3) */
+    /*    *TvarsD[nsd]                      [1]=2    1             3 */
+    /*    *TnsdVar                          [2]=2 [1]=1         [3]=3 */
+    /*    *TvarsDind[nsd](=k)               [1]=3 [2]=4         [3]=6 */
+    /*    *Tage[]                  [1]=1                  [2]=2      [3]=3 */
+    /*    *Tvard[]       [1][1]=1                                           [2][1]=1 */
+    /*    *                   [1][2]=3                                           [2][2]=2 */
+    /*    *Tprod[](=k)     [1]=1                                              [2]=8 */
+    /*    *TvarsDp(=Tvar)   [1]=1            [2]=2             [3]=3          [4]=5 */
+    /*    *TvarD (=k)       [1]=1            [2]=3 [3]=4       [3]=6          [4]=6 */
+    /*    *TvarsDpType */
+    /*    *si model= 1 + age + V3 + V2*age + V2 + V3*age */
+    /*    * nsd=1              (1)           (2) */
+    /*    *TvarsD[nsd]          3             2 */
+    /*    *TnsdVar           (3)=1          (2)=2 */
+    /*    *TvarsDind[nsd](=k)  [1]=1        [2]=3 */
+    /*    *Tage[]                  [1]=2           [2]= 3    */
+    /*    *\/ */
+    /*   /\* cov[++k1]=nbcode[TvarsD[k]][codtabm(ij,k)]; *\/ */
+    /*   /\* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); *\/ */
+    /* } */
+    /* for (k=1; k<=nsq;k++) { /\* For single quantitative varying covariates only of the model *\/ */
+    /* 			/\* Here comes the value of quantitative after renumbering k with single quantitative covariates *\/ */
+    /*   /\* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline                                 *\/ */
+    /*   /\* cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; *\/ */
+    /*   cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][resultmodel[nres][k1]] */
+    /*   /\* cov[++k1]=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 product with age *\/ */
+    /*   if(Dummy[Tage[k]]==2){ /\* dummy with age *\/ */
+    /* 	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,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 *\/ */
+    /* 	cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; */
+    /* 	/\* cov[++k1]=Tqresult[nres][k];  *\/ */
+    /*   } */
+    /*   /\* 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 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]); *\/ */
+    /*   if(Dummy[Tvard[k][1]]==0){ */
+    /* 	if(Dummy[Tvard[k][2]]==0){ */
+    /* 	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+    /* 	  /\* cov[++k1]=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,Tvard[k][1])] * Tqresult[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(ij,Tvard[k][2])] * Tqinvresult[nres][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][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */
+    /* 	  /\* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; *\/ */
+    /* 	} */
+    /*   } */
+    /* } /\* End product without age *\/ */
+/* ENd of old code */
     /*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]);*/
     /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
@@ -2948,7 +3012,7 @@ void powell(double p[], double **xi, int
   /*  0.51326036147820708, 0.48673963852179264} */
   /* If we start from prlim again, prlim tends to a constant matrix */
 
-  int i, ii,j,k;
+  int i, ii,j,k, k1;
   int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */
@@ -2988,50 +3052,60 @@ void powell(double p[], double **xi, int
     if(nagesqr==1){
       cov[3]= agefin*agefin;;
     }
-    for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
- 			/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
-      cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])];
-      /* printf("bprevalim Dummy agefin=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agefin,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
-    }
-    /* for (k=1; k<=cptcovn;k++) { */
-    /*   /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
-    /*   cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
-    /*   /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\/ */
-    /* } */
-    for (k=1; k<=nsq;k++) { /* For single varying covariates only */
- 			/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
-      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++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; */
-    /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */
-    /*   /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */
-    /*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-    for (k=1; k<=cptcovage;k++){  /* For product with age */
-      /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ ERROR ???*/
-      if(Dummy[Tage[k]]== 2){ /* dummy with age */
-	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
-      } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
-	cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
-      }
-      /* 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 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]); */
-      if(Dummy[Tvard[k][1]]==0){
-	if(Dummy[Tvard[k][2]]==0){
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
-	}else{
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * Tqresult[nres][k];
-	}
+    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(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]];
-	}else{
-	  cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
-	}
+	cov[2+nagesqr+k1]=precov[nres][k1];
       }
-    }
+    }/* End of loop on model equation */
+
+/* Old code */ 
+
+    /* for (k=1; k<=nsd;k++) { /\* For single dummy covariates only *\/ */
+    /* 			/\* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates *\/ */
+    /*   cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])]; */
+    /*   /\* printf("bprevalim Dummy agefin=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agefin,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); *\/ */
+    /* } */
+    /* /\* for (k=1; k<=cptcovn;k++) { *\/ */
+    /* /\*   /\\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\\/ *\/ */
+    /* /\*   cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; *\/ */
+    /* /\*   /\\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\\/ *\/ */
+    /* /\* } *\/ */
+    /* for (k=1; k<=nsq;k++) { /\* For single varying covariates only *\/ */
+    /* 			/\* Here comes the value of quantitative after renumbering k with single quantitative covariates *\/ */
+    /*   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++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; *\/ */
+    /* /\* for (k=1; k<=cptcovprod;k++) /\\* Useless *\\/ *\/ */
+    /* /\*   /\\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\\/ *\/ */
+    /* /\*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
+    /* for (k=1; k<=cptcovage;k++){  /\* For product with age *\/ */
+    /*   /\* if(Dummy[Tvar[Tage[k]]]== 2){ /\\* dummy with age *\\/ ERROR ???*\/ */
+    /*   if(Dummy[Tage[k]]== 2){ /\* dummy with age *\/ */
+    /* 	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+    /*   } else if(Dummy[Tage[k]]== 3){ /\* quantitative with age *\/ */
+    /* 	cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; */
+    /*   } */
+    /*   /\* 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 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]); *\/ */
+    /*   if(Dummy[Tvard[k][1]]==0){ */
+    /* 	if(Dummy[Tvard[k][2]]==0){ */
+    /* 	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+    /* 	}else{ */
+    /* 	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * Tqresult[nres][k]; */
+    /* 	} */
+    /*   }else{ */
+    /* 	if(Dummy[Tvard[k][2]]==0){ */
+    /* 	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])] * 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]);*/
@@ -3407,7 +3481,7 @@ double **matprod2(double **out, double *
 
 double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres )
 {
-  /* Computes the transition matrix starting at age 'age' and combination of covariate values corresponding to ij over 
+  /* Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over 
      'nhstepm*hstepm*stepm' months (i.e. until
      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying 
      nhstepm*hstepm matrices. 
@@ -3445,89 +3519,103 @@ double ***hpxij(double ***po, int nhstep
       /* Model(2)  V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */
       /* total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age */
       for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */ 
-	if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy  */
-/*	   V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */
-/*       for (k=1; k<=nsd;k++) { /\* For single dummy covariates only *\/ */
-/* /\* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates *\/ */
-	/* codtabm(ij,k)  (1 & (ij-1) >> (k-1))+1 */
-/*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-/*    k        1  2   3   4     5    6    7     8    9 */
-/*Tvar[k]=     5  4   3   6     5    2    7     1    1 */
-/*    nsd         1   2                              3 */ /* Counting single dummies covar fixed or tv */
-/*TvarsD[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
-/*TvarsDind[k]    2   3                              9 */ /* position K of single dummy cova */
-	  /* cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];or [codtabm(ij,TnsdVar[TvarsD[k]] */
-	  cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]];
-	  /* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,TnsdVar[TvarsD[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,TnsdVar[TvarsD[k]])],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,TnsdVar[TvarsD[k]])); */
-	  printf("hpxij Dummy combi=%d k1=%d Tvar[%d]=V%d cov[2+%d+%d]=%lf resultmodel[nres][%d]=%d nres/nresult=%d/%d \n",ij,k1,k1, Tvar[k1],nagesqr,k1,cov[2+nagesqr+k1],k1,resultmodel[nres][k1],nres,nresult);
-	}else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative variables  */
-	  /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline */
-	  cov[2+nagesqr+k1]=Tqresult[nres][resultmodel[nres][k1]]; 
-	  /* for (k=1; k<=nsq;k++) { /\* For single varying covariates only *\/ */
-	  /* 	/\* Here comes the value of quantitative after renumbering k with single quantitative covariates *\/ */
-	  /* 	cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; */
-	  printf("hPxij Quantitative k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]);
-	}else if( Dummy[k1]==2 ){ /* For dummy with age product */
-	  /* Tvar[k1] Variable in the age product age*V1 is 1 */
-	  /* [Tinvresult[nres][V1] is its value in the resultline nres */
-	  cov[2+nagesqr+k1]=Tinvresult[nres][Tvar[k1]];
-	  printf("DhPxij Dummy with age k1=%d Tvar[%d]=%d Tinvresult[nres][%d]=%d,cov[2+%d+%d]=%.3f\n",k1,k1,Tvar[k1],Tinvresult[nres][Tvar[k1]],nagesqr,k1,cov[2+nagesqr+k1]);
-	  /* cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]];	 */
-	  /* for (k=1; k<=cptcovage;k++){ /\* For product with age V1+V1*age +V4 +age*V3 *\/ */
-	  /* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/
-	  /* */
-/*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-/*    k        1  2   3   4     5    6    7     8    9 */
-/*Tvar[k]=     5  4   3   6     5    2    7     1    1 */
-/*cptcovage=2                   1               2      */
-/*Tage[k]=                      5               8      */	
-	}else if( Dummy[k1]==2 ){ /* For quant with age product */
-	  cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]];	
-	  printf("QhPxij Quant with age k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]);
-	  /* if(Dummy[Tage[k]]== 2){ /\* dummy with age *\/ */
-	  /* /\* if(Dummy[Tvar[Tage[k]]]== 2){ /\\* dummy with age *\\/ *\/ */
-	  /*   /\* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */
-	  /*   /\* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; *\/ */
-	  /*   cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; */
-	  /*   printf("hPxij Age combi=%d k=%d cptcovage=%d Tage[%d]=%d Tvar[Tage[%d]]=V%d nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]]])]=%d nres=%d\n",ij,k,cptcovage,k,Tage[k],k,Tvar[Tage[k]], nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]])],nres); */
-	  /* } else if(Dummy[Tage[k]]== 3){ /\* quantitative with age *\/ */
-	  /*   cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; */
-	  /* } */
-	  /* printf("hPxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
-	}else if(Typevar[k1]==2 ){ /* For product (not with age) */
-/*       for (k=1; k<=cptcovprod;k++){ /\*  For product without age *\/ */
+	if(Typevar[k1]==1){ /* A product with age */
+	  cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
+	}else{
+	  cov[2+nagesqr+k1]=precov[nres][k1];
+	}
+      }/* End of loop on model equation */
+	/* Old code */ 
+/* 	if( Dummy[k1]==0 && Typevar[k1]==0 ){ /\* Single dummy  *\/ */
+/* /\*	   V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) *\/ */
+/* /\*       for (k=1; k<=nsd;k++) { /\\* For single dummy covariates only *\\/ *\/ */
+/* /\* /\\* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates *\\/ *\/ */
+/* 	/\* codtabm(ij,k)  (1 & (ij-1) >> (k-1))+1 *\/ */
 /* /\*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
 /* /\*    k        1  2   3   4     5    6    7     8    9 *\/ */
 /* /\*Tvar[k]=     5  4   3   6     5    2    7     1    1 *\/ */
-/* /\*cptcovprod=1            1               2            *\/ */
-/* /\*Tprod[]=                4               7            *\/ */
-/* /\*Tvard[][1]             4               1             *\/ */
-/* /\*Tvard[][2]               3               2           *\/ */
+/* /\*    nsd         1   2                              3 *\/ /\* Counting single dummies covar fixed or tv *\/ */
+/* /\*TvarsD[nsd]     4   3                              1 *\/ /\* ID of single dummy cova fixed or timevary*\/ */
+/* /\*TvarsDind[k]    2   3                              9 *\/ /\* position K of single dummy cova *\/ */
+/* 	  /\* cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];or [codtabm(ij,TnsdVar[TvarsD[k]] *\/ */
+/* 	  cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]]; */
+/* 	  /\* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,TnsdVar[TvarsD[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,TnsdVar[TvarsD[k]])],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,TnsdVar[TvarsD[k]])); *\/ */
+/* 	  printf("hpxij Dummy combi=%d k1=%d Tvar[%d]=V%d cov[2+%d+%d]=%lf resultmodel[nres][%d]=%d nres/nresult=%d/%d \n",ij,k1,k1, Tvar[k1],nagesqr,k1,cov[2+nagesqr+k1],k1,resultmodel[nres][k1],nres,nresult); */
+/* 	  printf("hpxij new Dummy precov[nres=%d][k1=%d]=%.4f\n", nres, k1, precov[nres][k1]); */
+/* 	}else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /\* Single quantitative variables  *\/ */
+/* 	  /\* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline *\/ */
+/* 	  cov[2+nagesqr+k1]=Tqresult[nres][resultmodel[nres][k1]];  */
+/* 	  /\* for (k=1; k<=nsq;k++) { /\\* For single varying covariates only *\\/ *\/ */
+/* 	  /\* 	/\\* Here comes the value of quantitative after renumbering k with single quantitative covariates *\\/ *\/ */
+/* 	  /\* 	cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; *\/ */
+/* 	  printf("hPxij Quantitative k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]); */
+/* 	  printf("hpxij new Quanti precov[nres=%d][k1=%d]=%.4f\n", nres, k1, precov[nres][k1]); */
+/* 	}else if( Dummy[k1]==2 ){ /\* For dummy with age product *\/ */
+/* 	  /\* Tvar[k1] Variable in the age product age*V1 is 1 *\/ */
+/* 	  /\* [Tinvresult[nres][V1] is its value in the resultline nres *\/ */
+/* 	  cov[2+nagesqr+k1]=TinvDoQresult[nres][Tvar[k1]]*cov[2]; */
+/* 	  printf("DhPxij Dummy with age k1=%d Tvar[%d]=%d TinvDoQresult[nres=%d][%d]=%.f age=%.2f,cov[2+%d+%d]=%.3f\n",k1,k1,Tvar[k1],nres,TinvDoQresult[nres][Tvar[k1]],cov[2],nagesqr,k1,cov[2+nagesqr+k1]); */
+/* 	  printf("hpxij new Dummy with age product precov[nres=%d][k1=%d]=%.4f * age=%.2f\n", nres, k1, precov[nres][k1], cov[2]); */
+
+/* 	  /\* cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]];	 *\/ */
+/* 	  /\* for (k=1; k<=cptcovage;k++){ /\\* For product with age V1+V1*age +V4 +age*V3 *\\/ *\/ */
+/* 	  /\* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*\/ */
+/* 	  /\* *\/ */
+/* /\*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
+/* /\*    k        1  2   3   4     5    6    7     8    9 *\/ */
+/* /\*Tvar[k]=     5  4   3   6     5    2    7     1    1 *\/ */
+/* /\*cptcovage=2                   1               2      *\/ */
+/* /\*Tage[k]=                      5               8      *\/	 */
+/* 	}else if( Dummy[k1]==3 ){ /\* For quant with age product *\/ */
+/* 	  cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]];	 */
+/* 	  printf("QhPxij Quant with age k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]); */
+/* 	  printf("hpxij new Quanti with age product precov[nres=%d][k1=%d] * age=%.2f\n", nres, k1, precov[nres][k1], cov[2]); */
+/* 	  /\* if(Dummy[Tage[k]]== 2){ /\\* dummy with age *\\/ *\/ */
+/* 	  /\* /\\* if(Dummy[Tvar[Tage[k]]]== 2){ /\\\* dummy with age *\\\/ *\\/ *\/ */
+/* 	  /\*   /\\* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\\/ *\/ */
+/* 	  /\*   /\\* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; *\\/ *\/ */
+/* 	  /\*   cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; *\/ */
+/* 	  /\*   printf("hPxij Age combi=%d k=%d cptcovage=%d Tage[%d]=%d Tvar[Tage[%d]]=V%d nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]]])]=%d nres=%d\n",ij,k,cptcovage,k,Tage[k],k,Tvar[Tage[k]], nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]])],nres); *\/ */
+/* 	  /\* } else if(Dummy[Tage[k]]== 3){ /\\* quantitative with age *\\/ *\/ */
+/* 	  /\*   cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; *\/ */
+/* 	  /\* } *\/ */
+/* 	  /\* printf("hPxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); *\/ */
+/* 	}else if(Typevar[k1]==2 ){ /\* For product (not with age) *\/ */
+/* /\*       for (k=1; k<=cptcovprod;k++){ /\\*  For product without age *\\/ *\/ */
+/* /\* /\\*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\\/ *\/ */
+/* /\* /\\*    k        1  2   3   4     5    6    7     8    9 *\\/ *\/ */
+/* /\* /\\*Tvar[k]=     5  4   3   6     5    2    7     1    1 *\\/ *\/ */
+/* /\* /\\*cptcovprod=1            1               2            *\\/ *\/ */
+/* /\* /\\*Tprod[]=                4               7            *\\/ *\/ */
+/* /\* /\\*Tvard[][1]             4               1             *\\/ *\/ */
+/* /\* /\\*Tvard[][2]               3               2           *\\/ *\/ */
 	  
-	  /* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]=%d nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])]=%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2],nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])],nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]); */
-	  /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-	  cov[2+nagesqr+k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];	
-	  printf("hPxij Prod ij=%d k1=%d  cov[2+%d+%d]=%.5f Tvard[%d][1]=V%d * Tvard[%d][2]=V%d ; TinvDoQresult[nres][Tvardk[k1][1]]=%.4f * TinvDoQresult[nres][Tvardk[k1][1]]=%.4f\n",ij,k1,nagesqr,k1,cov[2+nagesqr+k1],k1,Tvard[k1][1], k1,Tvard[k1][2], TinvDoQresult[nres][Tvardk[k1][1]], TinvDoQresult[nres][Tvardk[k1][2]]);
-	  /* if(Dummy[Tvardk[k1][1]]==0){ */
-	  /*   if(Dummy[Tvardk[k1][2]]==0){ /\* Product of dummies *\/ */
-	      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-	      /* cov[2+nagesqr+k1]=Tinvresult[nres][Tvardk[k1][1]] * Tinvresult[nres][Tvardk[k1][2]];	 */
-	      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])]; */
-	    /* }else{ /\* Product of dummy by quantitative *\/ */
-	      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * Tqresult[nres][k]; */
-	      /* cov[2+nagesqr+k1]=Tresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]]; */
-	  /*   } */
-	  /* }else{ /\* Product of quantitative by...*\/ */
-	  /*   if(Dummy[Tvard[k][2]]==0){  /\* quant by dummy *\/ */
-	  /*     /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][Tvard[k][1]]; *\/ */
-	  /*     cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tresult[nres][Tinvresult[nres][Tvardk[k1][2]]]  ; */
-	  /*   }else{ /\* Product of two quant *\/ */
-	  /*     /\* cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; *\/ */
-	  /*     cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]]  ; */
-	  /*   } */
-	  /* }/\*end of products quantitative *\/ */
-	}/*end of products */
-      } /* End of loop on model equation */
+/* 	  /\* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]=%d nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])]=%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2],nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])],nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]); *\/ */
+/* 	  /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
+/* 	  cov[2+nagesqr+k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];	 */
+/* 	  printf("hPxij Prod ij=%d k1=%d  cov[2+%d+%d]=%.5f Tvard[%d][1]=V%d * Tvard[%d][2]=V%d ; TinvDoQresult[nres][Tvardk[k1][1]]=%.4f * TinvDoQresult[nres][Tvardk[k1][1]]=%.4f\n",ij,k1,nagesqr,k1,cov[2+nagesqr+k1],k1,Tvardk[k1][1], k1,Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][1]], TinvDoQresult[nres][Tvardk[k1][2]]); */
+/* 	  printf("hpxij new Product no age product precov[nres=%d][k1=%d]=%.4f\n", nres, k1, precov[nres][k1]); */
+
+/* 	  /\* if(Dummy[Tvardk[k1][1]]==0){ *\/ */
+/* 	  /\*   if(Dummy[Tvardk[k1][2]]==0){ /\\* Product of dummies *\\/ *\/ */
+/* 	      /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
+/* 	      /\* cov[2+nagesqr+k1]=Tinvresult[nres][Tvardk[k1][1]] * Tinvresult[nres][Tvardk[k1][2]];	 *\/ */
+/* 	      /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])]; *\/ */
+/* 	    /\* }else{ /\\* Product of dummy by quantitative *\\/ *\/ */
+/* 	      /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * Tqresult[nres][k]; *\/ */
+/* 	      /\* cov[2+nagesqr+k1]=Tresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]]; *\/ */
+/* 	  /\*   } *\/ */
+/* 	  /\* }else{ /\\* Product of quantitative by...*\\/ *\/ */
+/* 	  /\*   if(Dummy[Tvard[k][2]]==0){  /\\* quant by dummy *\\/ *\/ */
+/* 	  /\*     /\\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][Tvard[k][1]]; *\\/ *\/ */
+/* 	  /\*     cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tresult[nres][Tinvresult[nres][Tvardk[k1][2]]]  ; *\/ */
+/* 	  /\*   }else{ /\\* Product of two quant *\\/ *\/ */
+/* 	  /\*     /\\* cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; *\\/ *\/ */
+/* 	  /\*     cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]]  ; *\/ */
+/* 	  /\*   } *\/ */
+/* 	  /\* }/\\*end of products quantitative *\\/ *\/ */
+/* 	}/\*end of products *\/ */
+      /* } /\* End of loop on model equation *\/ */
       /* for (k=1; k<=cptcovn;k++)  */
       /* 	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
       /* for (k=1; k<=cptcovage;k++) /\* Should start at cptcovn+1 *\/ */
@@ -3573,7 +3661,8 @@ double ***hpxij(double ***po, int nhstep
 /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
 double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij, int nres )
 {
-  /* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over
+  /* For dummy covariates given in each resultline (for historical, computes the corresponding combination ij),
+     computes the transition matrix starting at age 'age' over
      'nhstepm*hstepm*stepm' months (i.e. until
      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying
      nhstepm*hstepm matrices.
@@ -3585,7 +3674,7 @@ double ***hbxij(double ***po, int nhstep
      The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output
   */
 
-  int i, j, d, h, k;
+  int i, j, d, h, k, k1;
   double **out, cov[NCOVMAX+1], **bmij();
   double **newm, ***newmm;
   double agexact;
@@ -3611,47 +3700,59 @@ double ***hbxij(double ***po, int nhstep
         /* Debug */
       /* printf("hBxij age=%lf, agexact=%lf\n", age, agexact); */
       cov[2]=agexact;
-      if(nagesqr==1)
+      if(nagesqr==1){
 	cov[3]= agexact*agexact;
-      for (k=1; k<=nsd;k++){ /* For single dummy covariates only *//* cptcovn error */
-      /* 	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
-      /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
-	cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])];/* Bug valgrind */
-        /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
-      }
-      for (k=1; k<=nsq;k++) { /* For single varying covariates only */
-	/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
-	cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
-	/* printf("hPxij 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++){ /* Should start at cptcovn+1 *//* For product with age */
-	/* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age error!!!*\/ */
-	if(Dummy[Tage[k]]== 2){ /* dummy with age */
-	  cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
-	} else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
-	  cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; 
-	}
-	/* printf("hBxij 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++){ /* Useless because included in cptcovn */
-	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
-	if(Dummy[Tvard[k][1]]==0){
-	  if(Dummy[Tvard[k][2]]==0){
-	    cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])];
-	  }else{
-	    cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * Tqresult[nres][k];
-	  }
+      }
+      /** New code */
+      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(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]];
-	  }else{
-	    cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
-	  }
+	  cov[2+nagesqr+k1]=precov[nres][k1];
 	}
-      }			
-      /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
-      /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
-
+      }/* End of loop on model equation */
+      /** End of new code */
+  /** This was old code */
+      /* for (k=1; k<=nsd;k++){ /\* For single dummy covariates only *\//\* cptcovn error *\/ */
+      /* /\* 	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; *\/ */
+      /* /\* /\\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\\/ *\/ */
+      /* 	cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])];/\* Bug valgrind *\/ */
+      /*   /\* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); *\/ */
+      /* } */
+      /* for (k=1; k<=nsq;k++) { /\* For single varying covariates only *\/ */
+      /* 	/\* Here comes the value of quantitative after renumbering k with single quantitative covariates *\/ */
+      /* 	cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];  */
+      /* 	/\* printf("hPxij 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++){ /\* Should start at cptcovn+1 *\//\* For product with age *\/ */
+      /* 	/\* if(Dummy[Tvar[Tage[k]]]== 2){ /\\* dummy with age error!!!*\\/ *\/ */
+      /* 	if(Dummy[Tage[k]]== 2){ /\* dummy with age *\/ */
+      /* 	  cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+      /* 	} else if(Dummy[Tage[k]]== 3){ /\* quantitative with age *\/ */
+      /* 	  cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];  */
+      /* 	} */
+      /* 	/\* printf("hBxij 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++){ /\* Useless because included in cptcovn *\/ */
+      /* 	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+      /* 	if(Dummy[Tvard[k][1]]==0){ */
+      /* 	  if(Dummy[Tvard[k][2]]==0){ */
+      /* 	    cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])]; */
+      /* 	  }else{ */
+      /* 	    cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * Tqresult[nres][k]; */
+      /* 	  } */
+      /* 	}else{ */
+      /* 	  if(Dummy[Tvard[k][2]]==0){ */
+      /* 	    cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]]; */
+      /* 	  }else{ */
+      /* 	    cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */
+      /* 	  } */
+      /* 	} */
+      /* }			 */
+      /* /\*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*\/ */
+      /* /\*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*\/ */
+/** End of old code */
+      
       /* Careful transposed matrix */
       /* age is in cov[2], prevacurrent at beginning of transition. */
       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
@@ -4787,7 +4888,7 @@ void  freqsummary(char fileres[], double
 		  int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],	\
 		  int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies as well as proposing some starting values */
-  
+  /* Frequencies of any combination of dummy covariate used in the model equation */ 
   int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1;
   int iind=0, iage=0;
   int mi; /* Effective wave */
@@ -4904,7 +5005,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
       j=cptcovs; /* Other passes for the covariate values */
     }
     first=1;
-    for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all covariates combination of the model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
+    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 */
       posproptt=0.;
       /*printf("cptcovn=%d Tvaraff=%d", cptcovn,Tvaraff[1]);
 	scanf("%d", i);*/
@@ -4944,12 +5045,12 @@ Title=%s 
Datafile=%s Firstpass=%d La
 		/* }else  if(Tvaraff[z1] ==-10){ */
 		/* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
 		/* }else  */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/
-		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]){ /* for combination j1 of covariates */
+		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 */
-		  /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
-		     bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
-		     j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+		  /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", */
+		  /*   bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),*/
+                  /*   j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
 		  /* 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 */
@@ -4965,12 +5066,12 @@ Title=%s 
Datafile=%s Firstpass=%d La
 		for (z1=1; z1<=cptcovn; z1++) {
 		  if( Fixed[Tmodelind[z1]]==1){
 		    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
-		    if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) /* iv=1 to ntv, right modality. If covariate's 
+		    if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's 
 										      value is -1, we don't select. It differs from the 
 										      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,Tvaraff[z1])]) {
+		    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) {
 		      bool=0;
 		    }
 		  }
@@ -4998,10 +5099,10 @@ Title=%s 
Datafile=%s Firstpass=%d La
 		  freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
 		  for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean on known values only */
 		    if(!isnan(covar[ncovcol+z1][iind])){
-			idq[z1]=idq[z1]+weight[iind];
-			meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind];  /* Computes mean of quantitative with selected filter */
-			/* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; *//*error*/
-			stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]; /* *weight[iind];*/  /* Computes mean of quantitative with selected filter */
+		      idq[z1]=idq[z1]+weight[iind];
+		      meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind];  /* Computes mean of quantitative with selected filter */
+		      /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; *//*error*/
+		      stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]; /* *weight[iind];*/  /* Computes mean of quantitative with selected filter */
 		    }
 		  }
 		  /* if((int)agev[m][iind] == 55) */
@@ -5088,7 +5189,7 @@ Title=%s 
Datafile=%s Firstpass=%d La
       fprintf(ficresphtm,"
");
       if((cptcovn==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,Tvaraff[z1])]);
+      if(nj==2) for (z1=1; z1<=cptcovn; 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);
 	fprintf(ficresphtm, "| Age | Prev(%d) | N(%d) | N",i,i);
@@ -5168,7 +5269,7 @@ Title=%s | 
Datafile=%s Firstpass=%d La
         }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,Tvaraff[z1])]);
+            for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
           }
 	}
 	for(s1=1; s1 <=nlstate ; s1++){
@@ -5494,10 +5595,10 @@ void prevalence(double ***probs, double
 	for (z1=1; z1<=cptcoveff; z1++){
 	  if( Fixed[Tmodelind[z1]]==1){
 	    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
-	    if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) /* iv=1 to ntv, right modality */
+	    if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */
 	      bool=0;
 	  }else if( Fixed[Tmodelind[z1]]== 0)  /* fixed */
-	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) {
+	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) {
 	      bool=0;
 	    }
 	}
@@ -6345,10 +6446,10 @@ void  concatwav(int wav[], int **dh, int
    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
    fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");
    for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-     fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+     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,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);
@@ -6976,29 +7077,36 @@ To be simple, these graphs help to under
    tj = (int) pow(2,cptcoveff);
    if (cptcovn<1) {tj=1;ncodemax[1]=1;}
    j1=0;
-   for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates or only once*/
-     for(nres=1;nres <=1; nres++){ /* For each resultline */
-     /* for(nres=1;nres <=nresult; nres++){ /\* For each resultline *\/ */
+
+   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);
+     if(tj != 1 && TKresult[nres]!= j1)
+       continue;
+
+   /* for(j1=1; j1<=tj;j1++){  /\* For each valid combination of covariates or only once*\/ */
+     /* 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,Tvaraff[z1])]);
+       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#********** Variable "); 
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]);
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
        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,Tvaraff[z1])]);
+       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,Tvaraff[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,Tvaraff[z1])]);
+       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); 
@@ -7018,7 +7126,7 @@ To be simple, these graphs help to under
        /* 	 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,TvarsD[k])];
+	 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
@@ -7031,12 +7139,12 @@ To be simple, these graphs help to under
        /*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,Tvar[Tage[k]])]*cov[2];
+	   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]);
-	   exit(1);
-	     /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */
+	   /* 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]; */
@@ -7044,18 +7152,18 @@ To be simple, these graphs help to under
        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,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(j1,Tvard[k][2])];
+	     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,Tvard[k][1])] * Tqresult[nres][k];
+	     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,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]];
+	     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][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
+	     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]]; */
 	   }
 	 }
@@ -7246,8 +7354,8 @@ To be simple, these graphs help to under
 	 } /* k12 */
        } /*l1 */
      }/* k1 */
-   } /* loop on nres */
    }  /* loop on combination of covariates j1 */
+   } /* loop on nres */
    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
    free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
@@ -7278,7 +7386,7 @@ void printinghtml(char fileresu[], char
    fprintf(fichtm,"- \n");
    fprintf(fichtm,"
- - Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file)
 \n",
 	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
-   fprintf(fichtm,"
-  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file) ",
+   fprintf(fichtm,"
-  - Observed prevalence (cross-sectional prevalence) in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s (html file) ",
 	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm"));
    fprintf(fichtm,",  %s (text file) 
 \n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
    fprintf(fichtm,"\
@@ -7676,7 +7784,8 @@ void printinggnuplot(char fileresu[], ch
 	fprintf(ficgp,"\n# 1st: Forward (stable period) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
 	strcpy(gplotlabel,"(");
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+	  /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the value of the covariate corresponding to k1 combination *\/ */
+	  lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
@@ -7729,7 +7838,8 @@ void printinggnuplot(char fileresu[], ch
 	}else{
 	  kl=0;
 	  for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
-	    lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	    /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to k1 combination and kth covariate *\/ */
+	    lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	    /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	    /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	    /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
@@ -7757,11 +7867,13 @@ void printinggnuplot(char fileresu[], ch
 	  }else{
 	    kl=0;
 	    for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
-	      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	      /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to k1 combination and kth covariate *\/ */
+	      lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	      vlv= nbcode[Tvaraff[k]][lv];
+	      /* vlv= nbcode[Tvaraff[k]][lv]; */
+	      vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	      kl++;
 	      /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
 	      /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
@@ -7771,7 +7883,7 @@ void printinggnuplot(char fileresu[], ch
 		fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
 			2+cptcoveff*2+(cpt-1),  cpt );  /* 4 or 6 ?*/
 	      }else{
-		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]);
 		kl++;
 	      }
 	    } /* end covariate */
@@ -7811,11 +7923,13 @@ void printinggnuplot(char fileresu[], ch
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       strcpy(gplotlabel,"(");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	/* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate number corresponding to k1 combination *\/ */
+	lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv];
+	/* vlv= nbcode[Tvaraff[k]][lv]; */
+	vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
       }
@@ -7874,21 +7988,23 @@ void printinggnuplot(char fileresu[], ch
       if(m != 1 && TKresult[nres]!= k1)
 	continue;
 
-      for (cpt=1; cpt<= nlstate ; cpt ++) {
+      for (cpt=1; cpt<= nlstate ; cpt ++) { /* Fragile no verification of covariate values */
 	fprintf(ficgp,"\n\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
 	strcpy(gplotlabel,"(");
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate number corresponding to k1 combination *\/ */
+	  lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /* Should be the covariate value corresponding to combination k1 and covariate k */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[k]][lv];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -7932,17 +8048,19 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	strcpy(gplotlabel,"(");
 	fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
+	  /* 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];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -7983,17 +8101,19 @@ set ter svg size 640, 480\nunset log y\n
 	strcpy(gplotlabel,"(");
 	fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
+	  /* 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];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -8042,17 +8162,19 @@ set ter svg size 640, 480\nunset log y\n
       strcpy(gplotlabel,"(");      
       fprintf(ficgp,"\n#\n#\n#CV preval stable (forward): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	/* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate number corresponding to k1 combination *\/ */
+	lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv];
+	/* vlv= nbcode[Tvaraff[k]][lv]; */
+	vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
       }	
       strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
@@ -8093,17 +8215,19 @@ set ter svg size 640, 480\nunset log y\n
 	strcpy(gplotlabel,"(");      
  	fprintf(ficgp,"\n#\n#\n#CV Backward stable prevalence: 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate number corresponding to k1 combination *\/ */
+	  lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[k]][lv];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -8149,17 +8273,19 @@ set ter svg size 640, 480\nunset log y\n
 	strcpy(gplotlabel,"(");      
 	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to forward stable prevalence (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	  /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to k1 combination and kth covariate *\/ */
+	  lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[k]][lv];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -8215,11 +8341,13 @@ set ter svg size 640, 480\nunset log y\n
 	    kl=0;
 	    strcpy(gplotcondition,"(");
 	    for (k=1; k<=cptcoveff; k++){    /* For each covariate writing the chain of conditions */
-	      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+	      /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
+	      lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
 	      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	      vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
+	      /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
+	      vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	      kl++;
 	      sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
 	      kl++;
@@ -8262,17 +8390,19 @@ set ter svg size 640, 480\nunset log y\n
 	strcpy(gplotlabel,"(");      
 	fprintf(ficgp,"\n#\n#\n#Back projection of prevalence to stable (mixed) back prevalence: 'BPROJ_' files, covariatecombination#=%d originstate=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
-	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	  /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to k1 combination and kth covariate *\/ */
+	  lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /* Should be the covariate value corresponding to combination k1 and covariate k */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	  vlv= nbcode[Tvaraff[k]][lv];
+	  /* vlv= nbcode[Tvaraff[k]][lv]; */
+	  vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	  sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
 	}	
 	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
@@ -8328,11 +8458,13 @@ set ter svg size 640, 480\nunset log y\n
 	    kl=0;
 	    strcpy(gplotcondition,"(");
 	    for (k=1; k<=cptcoveff; k++){    /* For each covariate writing the chain of conditions */
-	      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+	      /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
+	      lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /* Should be the covariate value corresponding to combination k1 and covariate k */
 	      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	      vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
+	      /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
+	      vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	      kl++;
 	      sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
 	      kl++;
@@ -8413,17 +8545,19 @@ set ter svg size 640, 480\nunset log y\n
       strcpy(gplotlabel,"(");
       /*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/
       for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	/* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to k1 combination and kth covariate *\/ */
+	lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /* Should be the covariate value corresponding to combination k1 and covariate k */
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv];
+	/* vlv= nbcode[Tvaraff[k]][lv]; */
+	vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	sprintf(gplotlabel+strlen(gplotlabel)," 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]);
-	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficgp," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][resultmodel[nres][k4]],Tqresult[nres][resultmodel[nres][k4]]);
       }	
       strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
@@ -8963,7 +9097,7 @@ void prevforecast(char fileres[], double
   
 /* 	      if (h==(int)(YEARM*yearp)){ */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */
-  for(k=1; k<=i1;k++){
+    for(k=1; k<=i1;k++){ /* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) */
     if(i1 != 1 && TKresult[nres]!= k)
       continue;
     if(invalidvarcomb[k]){
@@ -8972,7 +9106,8 @@ void prevforecast(char fileres[], double
     }
     fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
     for(j=1;j<=cptcoveff;j++) {
-      fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+      /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */
+      fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
     }
     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
       fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
@@ -9002,7 +9137,8 @@ void prevforecast(char fileres[], double
 	}
 	fprintf(ficresf,"\n");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	  /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Tvaraff not correct *\/ */
+	  fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* TnsdVar[Tvaraff]  correct */
 	fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
 	
 	for(j=1; j<=nlstate+ndeath;j++) {
@@ -9112,7 +9248,7 @@ void prevforecast(char fileres[], double
     }
     fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
     for(j=1;j<=cptcoveff;j++) {
-      fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+      fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
     }
     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
       fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
@@ -9148,7 +9284,7 @@ void prevforecast(char fileres[], double
 	}
 	fprintf(ficresfb,"\n");
 	for(j=1;j<=cptcoveff;j++)
-	  fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	  fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
 	fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
 	for(i=1; i<=nlstate+ndeath;i++) {
 	  ppij=0.;ppi=0.;
@@ -9208,21 +9344,21 @@ void prevforecast(char fileres[], double
     if (cptcovn < 1){i1=1;}
 
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
-    for(k=1; k<=i1;k++){
+      for(k=1; k<=i1;k++){ /* We find the combination equivalent to result line values of dummies */
       if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficresvpl,"\n#****** ");
       printf("\n#****** ");
       fprintf(ficlog,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficresvpl,"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]])]);
+	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][j]);
-	fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+	fprintf(ficresvpl," 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(ficresvpl,"******\n");
       printf("******\n");
@@ -9272,14 +9408,14 @@ void prevforecast(char fileres[], double
        printf("\n#****** ");
        fprintf(ficlog,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) {
-	 fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	 fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	 printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	 fprintf(ficresvbl,"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]])]);
+	 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][j]);
-	 fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	 fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	 printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+	 fprintf(ficresvbl," 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(ficresvbl,"******\n");
        printf("******\n");
@@ -10114,13 +10250,10 @@ int decoderesult( char resultline[], int
   char stra[80], strb[80], strc[80], strd[80],stre[80];
 
   removefirstspace(&resultline);
+  printf("decoderesult:%s\n",resultline);
 
-  if (strstr(resultline,"v") !=0){
-    printf("Error. 'v' must be in upper case 'V' result: %s ",resultline);
-    fprintf(ficlog,"Error. 'v' must be in upper case result: %s ",resultline);fflush(ficlog);
-    return 1;
-  }
-  trimbb(resultsav, resultline);
+  strcpy(resultsav,resultline);
+  printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline);
   if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
   }
@@ -10129,13 +10262,21 @@ int decoderesult( char resultline[], int
     return (0);
   }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
-    printf("ERROR: the number of variables in this result line, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs);
-    fprintf(ficlog,"ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs);
+    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);
+    /* return 1;*/
   }
   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" */
       cutl(strc,strd,strb,'=');  /* strb:"V4=1" strc="1" strd="V4" */
+      /* If a blank, then strc="V4=" and strd='\0' */
+      if(strc[0]=='\0'){
+      printf("Error in resultline, probably a blank after the \"%s\", \"result:%s\", stra=\"%s\" resultsav=\"%s\"\n",strb,resultline, stra, resultsav);
+	fprintf(ficlog,"Error in resultline, probably a blank after the \"V%s=\", resultline=%s\n",strb,resultline);
+	return 1;
+      }
     }else
       cutl(strc,strd,resultsav,'=');
     Tvalsel[k]=atof(strc); /* 1 */ /* Tvalsel of k is the float value of the kth covariate appearing in this result line */
@@ -10148,8 +10289,10 @@ int decoderesult( char resultline[], int
       strcpy(resultsav,stra); /* and analyzes it */
   }
   /* 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 */
-    if(Typevar[k1]==0){ /* Single covariate in model *//*0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
+    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   */
@@ -10159,17 +10302,65 @@ int decoderesult( char resultline[], int
 	}
       }
       if(match == 0){
-	printf("Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model);
-	fprintf(ficlog,"Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model);
+	printf("Error in result line (Dummy single): V%d is missing in result: %s according to model=%s. Tvar[k1=%d]=%d is different from Tvarsel[k2=%d]=%d.\n",Tvar[k1], resultline, model,k1, Tvar[k1], k2, Tvarsel[k2]);
+	fprintf(ficlog,"Error in result line (Dummy single): V%d is missing in result: %s according to model=%s\n",Tvar[k1], resultline, model);
 	return 1;
       }
-    }
-  }
+    }else if(Typevar[k1]==1){ /* Product with age We want to get the position k2 in the resultline of the product k1 in the model line*/
+      /* We feed resultmodel[k1]=k2; */
+      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 */
+	  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 */
+	  break;
+	}
+      }
+      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]);
+      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*/
+      /* resultmodel[nres][of such a Vn * Vm product k1] is not unique, so can't exist, we feed Tvard[k1][1] and [2] */ 
+      match=0;
+      printf("Decoderesult very first Product Tvardk[k1=%d][1]=%d Tvardk[k1=%d][2]=%d V%d * V%d\n",k1,Tvardk[k1][1],k1,Tvardk[k1][2],Tvardk[k1][1],Tvardk[k1][2]);
+      for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+	if(Tvardk[k1][1]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
+	  /* modelresult[k2]=k1; */
+	  printf("Decoderesult first Product modelresult[k2=%d]=%d (k1) V%d * \n",k2,k1,Tvarsel[k2]);
+	  match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
+	}
+      }
+      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);
+	return 1;
+      }
+      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(Tvardk[k1][2]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
+	  /* modelresult[k2]=k1;*/
+	  printf("Decoderesult second Product modelresult[k2=%d]=%d (k1) * V%d \n ",k2,k1,Tvarsel[k2]);
+	  match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
+	  break;
+	}
+      }
+      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);
+	return 1;
+      }
+    }/* End of testing */
+  }/* 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 */
     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 */
+      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 */
 	  ++match;
@@ -10177,8 +10368,8 @@ int decoderesult( char resultline[], int
       }
     }
     if(match == 0){
-      printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
-      fprintf(ficlog,"Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
+      printf("Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
+      fprintf(ficlog,"Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
       return 1;
     }else if(match > 1){
       printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
@@ -10206,42 +10397,68 @@ int decoderesult( char resultline[], int
   /* 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) */
-    if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */
-      /* k4+1= 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 correspond to the k1 position 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) */
+      /* 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 */
-      /*      k3 is the position in the nres result line of the k1th variable of the model equation                                          */
-      /* Tvarsel: Name of the variable at the k3th position in the result line Tvarsel[k3].                                                  */
-      /* Tvalsel: Value of the variable at the k3th position in the result line Tvarsel[k3].                                                 */
-      /* 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                                   */
-      /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line                                                      */
+      /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline                        */
+      /*      k3 is the position in the nres result line of the k1th variable of the model equation                                  */
+      /* 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                     */
+      /* 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]; /* 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)*/
+      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. */
+      /* 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 */
       Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
-      TinvDoQresult[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);
+      precov[nres][k1]=Tvalsel[k3];
+      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 ){ /* Single quantitative */
+    }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */
       /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline                                 */
-      /* Tqvresult[nres][result_position]= id of the 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                                 */
       /* Tqinvresult[nres][Name of a quantitative variable]= value of the variable in the result line                                                      */
-      k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
-      k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
+      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 */
       Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       TinvDoQresult[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]);
+      precov[nres][k1]=Tvalsel[k3q];
+      printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
       k4q++;;
+    }else if( Dummy[k1]==2 ){ /* For dummy with age product */
+      /* Tvar[k1]; */ /* Age variable */
+      /* Wrong we want the value of variable name Tvar[k1] */
+      
+      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 */
+      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 */
+      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]);
+    }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]]);
     }else{
-      printf("Decodemodel probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
-      fprintf(ficlog,"Decodemodel probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
+      printf("Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
+      fprintf(ficlog,"Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
     }
   }
   
@@ -10543,7 +10760,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       ncovv++; /* Only simple time varying variables */
       nsq++;
       TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */
-      TvarsQind[nsq]=k;
+      TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate */
       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 */
@@ -11069,6 +11286,7 @@ void syscompilerinfo(int logged)
 
 int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
   /*--------------- Prevalence limit  (forward period or forward stable prevalence) --------------*/
+  /* Computes the prevalence limit for each combination of the dummy covariates */
   int i, j, k, i1, k4=0, nres=0 ;
   /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;
@@ -11107,14 +11325,15 @@ int prevalence_limit(double *p, double *
       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       /* k=k+1; */
       /* to clean */
-      //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+      /*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*/
       fprintf(ficrespl,"#******");
       printf("#******");
       fprintf(ficlog,"#******");
       for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /* Here problem for varying dummy*/
-	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	/* fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Here problem for varying dummy*\/ */
+	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* Here problem for varying dummy*/
+	printf(" 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 (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
@@ -11133,7 +11352,7 @@ int prevalence_limit(double *p, double *
 
       fprintf(ficrespl,"#Age ");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
       }
       for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
       fprintf(ficrespl,"Total Years_to_converge\n");
@@ -11143,7 +11362,7 @@ int prevalence_limit(double *p, double *
 	prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres);
 	fprintf(ficrespl,"%.0f ",age );
 	for(j=1;j<=cptcoveff;j++)
-	  fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	  fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
 	tot=0.;
 	for(i=1; i<=nlstate;i++){
 	  tot +=  prlim[i][i];
@@ -11198,19 +11417,19 @@ int back_prevalence_limit(double *p, dou
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
-      //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+     /*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*/
       fprintf(ficresplb,"#******");
       printf("#******");
       fprintf(ficlog,"#******");
       for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-	fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+	printf(" 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][j]);
-	fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+	fprintf(ficresplb," 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(ficresplb,"******\n");
       printf("******\n");
@@ -11224,7 +11443,7 @@ int back_prevalence_limit(double *p, dou
     
       fprintf(ficresplb,"#Age ");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
       }
       for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
       fprintf(ficresplb,"Total Years_to_converge\n");
@@ -11248,7 +11467,7 @@ int back_prevalence_limit(double *p, dou
 	}
 	fprintf(ficresplb,"%.0f ",age );
 	for(j=1;j<=cptcoveff;j++)
-	  fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	  fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
 	tot=0.;
 	for(i=1; i<=nlstate;i++){
 	  tot +=  bprlim[i][i];
@@ -11306,7 +11525,7 @@ int hPijx(double *p, int bage, int fage)
 	continue;
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
-	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	fprintf(ficrespij," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
@@ -11385,9 +11604,9 @@ int hPijx(double *p, int bage, int fage)
 	continue;
       fprintf(ficrespijb,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)
-	fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficrespijb,"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(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
       }
       fprintf(ficrespijb,"******\n");
       if(invalidvarcomb[k]){  /* Is it necessary here? */
@@ -11480,7 +11699,7 @@ int main(int argc, char *argv[])
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
 
   char  modeltemp[MAXLINE];
-  char resultline[MAXLINE];
+  char resultline[MAXLINE], resultlineori[MAXLINE];
   
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
@@ -12388,7 +12607,7 @@ 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-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
 \
+  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\
 IMaCh-%s
 %s \
 \n\
@@ -13059,6 +13278,9 @@ Please run with mle=-1 to get a correct
     }
      
     /* Results */
+    /* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */
+    /* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */  
+    precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
     endishere=0;
     nresult=0;
     parameterline=0;
@@ -13132,14 +13354,24 @@ Please run with mle=-1 to get a correct
 	}
 	break;
       case 13:
-	num_filled=sscanf(line,"result:%[^\n]\n",resultline);
+	num_filled=sscanf(line,"result:%[^\n]\n",resultlineori);
 	nresult++; /* Sum of resultlines */
-	printf("Result %d: result:%s\n",nresult, resultline);
+	printf("Result %d: result:%s\n",nresult, resultlineori);
+	/* removefirstspace(&resultlineori); */
+	
+	if(strstr(resultlineori,"v") !=0){
+	  printf("Error. 'v' must be in upper case 'V' result: %s ",resultlineori);
+	  fprintf(ficlog,"Error. 'v' must be in upper case result: %s ",resultlineori);fflush(ficlog);
+	  return 1;
+	}
+	trimbb(resultline, resultlineori); /* Suppressing double blank in the resultline */
+	printf("Decoderesult resultline=\"%s\" resultlineori=\"%s\"\n", resultline, resultlineori);
 	if(nresult > MAXRESULTLINESPONE-1){
 	  printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
 	  fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
 	  goto end;
 	}
+	
 	if(!decoderesult(resultline, nresult)){ /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
 	  fprintf(ficparo,"result: %s\n",resultline);
 	  fprintf(ficres,"result: %s\n",resultline);
@@ -13243,16 +13475,19 @@ Please run with mle=-1 to get a correct
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
     prlim=matrix(1,nlstate,1,nlstate);
+    /* Computes the prevalence limit for each combination k of the dummy covariates by calling prevalim(k) */
     prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);
     fclose(ficrespl);
 
     /*------------- h Pij x at various ages ------------*/
     /*#include "hpijx.h"*/
+    /** h Pij x Probability to be in state j at age x+h being in i at x, for each combination k of dummies in the model line or to nres?*/
+    /* calls hpxij with combination k */
     hPijx(p, bage, fage);
     fclose(ficrespij);
     
     /* ncovcombmax=  pow(2,cptcoveff); */
-    /*-------------- Variance of one-step probabilities---*/
+    /*-------------- Variance of one-step probabilities for a combination ij or for nres ?---*/
     k=1;
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
     
@@ -13373,12 +13608,12 @@ Please run with mle=-1 to get a correct
       fprintf(ficreseij,"\n#****** ");
       printf("\n#****** ");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
+	fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+	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][j]);
-	fprintf(ficreseij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	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]]);
       }
       fprintf(ficreseij,"******\n");
       printf("******\n");
@@ -13443,14 +13678,14 @@ Please run with mle=-1 to get a correct
       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,Tvaraff[j])]);
-	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]);
-	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[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][j]);
-	fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	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]]);
       }	
       fprintf(ficrest,"******\n");
       fprintf(ficlog,"******\n");
@@ -13459,12 +13694,12 @@ 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,Tvaraff[j])]);
-	fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[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][j]);
-	fprintf(ficrescveij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	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,"******\n");
       fprintf(ficrescveij,"******\n");
@@ -13472,9 +13707,10 @@ 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,Tvaraff[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]);
+	/* 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 */
       }	
       fprintf(ficresvij,"******\n");
       
@@ -13636,7 +13872,9 @@ Please run with mle=-1 to get a correct
   free_ivector(Tmodelind,1,NCOVMAX);
   free_ivector(TmodelInvind,1,NCOVMAX);
   free_ivector(TmodelInvQind,1,NCOVMAX);
-  
+
+  free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
+
   free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
   /* free_imatrix(codtab,1,100,1,10); */
   fflush(fichtm);