--- imach/src/imach.c	2016/08/23 07:40:50	1.233
+++ imach/src/imach.c	2016/08/29 07:53:18	1.240
@@ -1,6 +1,29 @@
-/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 brouard Exp $
+/* $Id: imach.c,v 1.240 2016/08/29 07:53:18 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.240  2016/08/29 07:53:18  brouard
+  Summary: Better
+
+  Revision 1.239  2016/08/26 15:51:03  brouard
+  Summary: Improvement in Powell output in order to copy and paste
+
+  Author:
+
+  Revision 1.238  2016/08/26 14:23:35  brouard
+  Summary: Starting tests of 0.99
+
+  Revision 1.237  2016/08/26 09:20:19  brouard
+  Summary: to valgrind
+
+  Revision 1.236  2016/08/25 10:50:18  brouard
+  *** empty log message ***
+
+  Revision 1.235  2016/08/25 06:59:23  brouard
+  *** empty log message ***
+
+  Revision 1.234  2016/08/23 16:51:20  brouard
+  *** empty log message ***
+
   Revision 1.233  2016/08/23 07:40:50  brouard
   Summary: not working
 
@@ -899,12 +922,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.233 2016/08/23 07:40:50 brouard Exp $ */
+/* $Id: imach.c,v 1.240 2016/08/29 07:53:18 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.233 $ $Date: 2016/08/23 07:40:50 $"; 
+char fullversion[]="$Revision: 1.240 $ $Date: 2016/08/29 07:53:18 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -920,7 +943,8 @@ int cptcoveff=0; /* Total number of cova
 int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
 int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
 int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
-
+int nsd=0; /**< Total number of single dummy variables (output) */
+int nsq=0; /**< Total number of single quantitative variables (output) */
 int ncoveff=0; /* Total number of effective fixed dummy covariates in the model */
 int nqfveff=0; /**< nqfveff Number of Quantitative Fixed Variables Effective */
 int ntveff=0; /**< ntveff number of effective time varying variables */
@@ -978,6 +1002,7 @@ char fileresv[FILENAMELENGTH];
 FILE  *ficresvpl;
 char fileresvpl[FILENAMELENGTH];
 char title[MAXLINE];
+char model[MAXLINE]; /**< The model line */
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];
 char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
 char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH]; 
@@ -1015,7 +1040,8 @@ double dval;
 #define FTOL 1.0e-10
 
 #define NRANSI 
-#define ITMAX 200 
+#define ITMAX 200
+#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */ 
 
 #define TOL 2.0e-4 
 
@@ -1077,6 +1103,40 @@ double ***cotvar; /* Time varying covari
 double ***cotqvar; /* Time varying quantitative covariate itqv */
 double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+/*           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 */
+/* Tndvar[k]    1   2   3               4          5 */
+/*TDvar         4   3   6               7          1 */ /* For outputs only; combination of dummies fixed or varying */
+/* Tns[k]    1  2   2              4               5 */ /* Number of single cova */
+/* TvarsD[k]    1   2                              3 */ /* Number of single dummy cova */
+/* TvarsDind    2   3                              9 */ /* position K of single dummy cova */
+/* TvarsQ[k] 1                     2                 */ /* Number of single quantitative cova */
+/* TvarsQind 1                     6                 */ /* position K of single quantitative cova */
+/* Tprod[i]=k           4               7            */
+/* Tage[i]=k                  5               8      */
+/* */
+/* Type                    */
+/* V         1  2  3  4  5 */
+/*           F  F  V  V  V */
+/*           D  Q  D  D  Q */
+/*                         */
+int *TvarsD;
+int *TvarsDind;
+int *TvarsQ;
+int *TvarsQind;
+
+#define MAXRESULTLINES 10
+int nresult=0;
+int TKresult[MAXRESULTLINES];
+int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */
+double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
+double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
+int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */
+
+/* int *TDvar; /\**< TDvar[1]=4,  TDvarF[2]=3, TDvar[3]=6  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
 int *TvarF; /**< TvarF[1]=Tvar[6]=2,  TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
 int *TvarFind; /**< TvarFind[1]=6,  TvarFind[2]=7, Tvarind[3]=9  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
 int *TvarV; /**< TvarV[1]=Tvar[1]=5, TvarV[2]=Tvar[2]=4  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
@@ -1097,6 +1157,8 @@ double *Tvalsel; /**< Selected modality
 int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
 int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ 
 int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ 
+int *DummyV; /** Dummy[v] 0=dummy (0 1), 1 quantitative */
+int *FixedV; /** FixedV[v] 0 fixed, 1 varying */
 int *Tage;
 int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */ 
 int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
@@ -1106,11 +1168,10 @@ int *Ndum; /** Freq of modality (tricode
 /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard;
 int *Tprod;/**< Gives the k position of the k1 product */
+/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3  */
 int *Tposprod; /**< Gives the k1 product from the k position */
-/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
-   if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2)
-   Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2 
-*/
+   /* if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2) */
+   /* Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5(V3*V2)]=2 (2nd product without age) */
 int cptcovprod, *Tvaraff, *invalidvarcomb;
 double *lsurv, *lpop, *tpop;
 
@@ -1336,7 +1397,7 @@ int nbocc(char *s, char occ)
   i=0;
   lg=strlen(s);
   for(i=0; i<= lg; i++) {
-  if  (s[i] == occ ) j++;
+    if  (s[i] == occ ) j++;
   }
   return j;
 }
@@ -2012,7 +2073,7 @@ void powell(double p[], double **xi, int
  void linmin(double p[], double xi[], int n, double *fret, 
 						 double (*func)(double []),int *flat); 
 #endif
-  int i,ibig,j; 
+ int i,ibig,j,jk,k; 
   double del,t,*pt,*ptt,*xit;
   double directest;
   double fp,fptt;
@@ -2044,13 +2105,46 @@ void powell(double p[], double **xi, int
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
 /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
     for (i=1;i<=n;i++) {
-      printf(" %d %.12f",i, p[i]);
-      fprintf(ficlog," %d %.12lf",i, p[i]);
       fprintf(ficrespow," %.12lf", p[i]);
     }
+    fprintf(ficrespow,"\n");fflush(ficrespow);
+    printf("\n#model=  1      +     age ");
+    fprintf(ficlog,"\n#model=  1      +     age ");
+    if(nagesqr==1){
+	printf("  + age*age  ",Tvar[j]);
+	fprintf(ficlog,"  + age*age  ",Tvar[j]);
+    }
+    for(j=1;j <=ncovmodel-2;j++){
+      if(Typevar[j]==0) {
+	printf("  +      V%d  ",Tvar[j]);
+	fprintf(ficlog,"  +      V%d  ",Tvar[j]);
+      }else if(Typevar[j]==1) {
+	printf("  +    V%d*age ",Tvar[j]);
+	fprintf(ficlog,"  +    V%d*age ",Tvar[j]);
+      }else if(Typevar[j]==2) {
+	printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+	fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+      }
+    }
     printf("\n");
+/*     printf("12   47.0114589    0.0154322   33.2424412    0.3279905    2.3731903  */
+/* 13  -21.5392400    0.1118147    1.2680506    1.2973408   -1.0663662  */
     fprintf(ficlog,"\n");
-    fprintf(ficrespow,"\n");fflush(ficrespow);
+    for(i=1,jk=1; i <=nlstate; i++){
+      for(k=1; k <=(nlstate+ndeath); k++){
+	if (k != i) {
+	  printf("%d%d ",i,k);
+	  fprintf(ficlog,"%d%d ",i,k);
+	  for(j=1; j <=ncovmodel; j++){
+	    printf("%12.7f ",p[jk]);
+	    fprintf(ficlog,"%12.7f ",p[jk]);
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
+      }
+    }
     if(*iter <=3){
       tml = *localtime(&rcurr_time);
       strcpy(strcurr,asctime(&tml));
@@ -2168,7 +2262,7 @@ void powell(double p[], double **xi, int
       free_vector(pt,1,n); 
       return; 
     } /* enough precision */ 
-    if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
+    if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
       ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j]; 
@@ -2227,87 +2321,89 @@ void powell(double p[], double **xi, int
       if (directest < 0.0) { /* Then we use it for new direction */
 #endif
 #ifdef DEBUGLINMIN
-				printf("Before linmin in direction P%d-P0\n",n);
-				for (j=1;j<=n;j++) {
-					printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-					fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-					if(j % ncovmodel == 0){
-						printf("\n");
-						fprintf(ficlog,"\n");
-					}
-				}
+	printf("Before linmin in direction P%d-P0\n",n);
+	for (j=1;j<=n;j++) {
+	  printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  if(j % ncovmodel == 0){
+	    printf("\n");
+	    fprintf(ficlog,"\n");
+	  }
+	}
 #endif
 #ifdef LINMINORIGINAL
- 				linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+	linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
 #else
-				linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
-				flatdir[i]=flat; /* Function is vanishing in that direction i */
+	linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+	flatdir[i]=flat; /* Function is vanishing in that direction i */
 #endif
-
+	
 #ifdef DEBUGLINMIN
-				for (j=1;j<=n;j++) { 
-					printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-					fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-					if(j % ncovmodel == 0){
-						printf("\n");
-						fprintf(ficlog,"\n");
-					}
-				}
+	for (j=1;j<=n;j++) { 
+	  printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  if(j % ncovmodel == 0){
+	    printf("\n");
+	    fprintf(ficlog,"\n");
+	  }
+	}
 #endif
-				for (j=1;j<=n;j++) { 
-					xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
-					xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
-				}
+	for (j=1;j<=n;j++) { 
+	  xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
+	  xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
+	}
 #ifdef LINMINORIGINAL
 #else
-				for (j=1, flatd=0;j<=n;j++) {
-					if(flatdir[j]>0)
-						flatd++;
-				}
-				if(flatd >0){
-					printf("%d flat directions\n",flatd);
-					fprintf(ficlog,"%d flat directions\n",flatd);
-					for (j=1;j<=n;j++) { 
-						if(flatdir[j]>0){
-							printf("%d ",j);
-							fprintf(ficlog,"%d ",j);
-						}
-					}
-					printf("\n");
-					fprintf(ficlog,"\n");
-				}
+	for (j=1, flatd=0;j<=n;j++) {
+	  if(flatdir[j]>0)
+	    flatd++;
+	}
+	if(flatd >0){
+	  printf("%d flat directions\n",flatd);
+	  fprintf(ficlog,"%d flat directions\n",flatd);
+	  for (j=1;j<=n;j++) { 
+	    if(flatdir[j]>0){
+	      printf("%d ",j);
+	      fprintf(ficlog,"%d ",j);
+	    }
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
 #endif
-				printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
-				fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
-				
+	printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+	fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+	
 #ifdef DEBUG
-				printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
-				fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
-				for(j=1;j<=n;j++){
-					printf(" %lf",xit[j]);
-					fprintf(ficlog," %lf",xit[j]);
-				}
-				printf("\n");
-				fprintf(ficlog,"\n");
+	printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+	fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+	for(j=1;j<=n;j++){
+	  printf(" %lf",xit[j]);
+	  fprintf(ficlog," %lf",xit[j]);
+	}
+	printf("\n");
+	fprintf(ficlog,"\n");
 #endif
       } /* end of t or directest negative */
 #ifdef POWELLNOF3INFF1TEST
 #else
-    } /* end if (fptt < fp)  */
+      } /* end if (fptt < fp)  */
 #endif
 #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */
-		} /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */
+    } /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */
 #else
 #endif
-  } /* loop iteration */ 
+		} /* loop iteration */ 
 } 
-
+  
 /**** Prevalence limit (stable or period prevalence)  ****************/
-
-double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
-{
-  /* Computes the prevalence limit in each live state at age x and for covariate ij by left multiplying the unit
-     matrix by transitions matrix until convergence is reached with precision ftolpl */
+  
+  double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres)
+  {
+    /* Computes the prevalence limit in each live state at age x and for covariate combination ij 
+       (and selected quantitative values in nres)
+       by left multiplying the unit
+       matrix by transitions matrix until convergence is reached with precision ftolpl */
   /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
   /* Wx is row vector: population in state 1, population in state 2, population dead */
   /* or prevalence in state 1, prevalence in state 2, 0 */
@@ -2325,7 +2421,7 @@ double **prevalim(double **prlim, int nl
   /* {0.51571254859325999, 0.4842874514067399, */
   /*  0.51326036147820708, 0.48673963852179264} */
   /* If we start from prlim again, prlim tends to a constant matrix */
-
+    
   int i, ii,j,k;
   double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */
@@ -2355,19 +2451,40 @@ double **prevalim(double **prlim, int nl
     cov[2]=agefin;
     if(nagesqr==1)
       cov[3]= agefin*agefin;;
-    for (k=1; k<=cptcovn;k++) {
-      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
-			/* Here comes the value of the covariate 'ij' */
-      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<=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,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]; 
+      /* 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[Tvar[Tage[k]]]){
+	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+      } else{
+	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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
+	}
+      }else{
+	if(Dummy[Tvard[k][2]==0]){
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
+	}
+      }
     }
-    /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
-    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)];
-    
     /*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]);*/
@@ -2489,8 +2606,6 @@ Earliest age to start was %d-%d=%d, ncvl
       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])]); */
     }
-    /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
     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])]; */
@@ -2517,10 +2632,10 @@ Earliest age to start was %d-%d=%d, ncvl
     }
     for(j=1; j<=nlstate; j++){ 
       for(i=1;i<=nlstate;i++){
-				/* bprlim[i][j]= newm[i][j]/(1-sumnew); */
-				bprlim[i][j]= newm[i][j];
-				max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
-				min[i]=FMIN(min[i],bprlim[i][j]);
+	/* bprlim[i][j]= newm[i][j]/(1-sumnew); */
+	bprlim[i][j]= newm[i][j];
+	max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
+	min[i]=FMIN(min[i],bprlim[i][j]);
       }
     }
 		
@@ -2715,81 +2830,81 @@ double **bpmij(double **ps, double *cov,
   /*double t34;*/
   int i,j, nc, ii, jj;
 
-	for(i=1; i<= nlstate; i++){
-		for(j=1; j<i;j++){
-			for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-				/*lnpijopii += param[i][j][nc]*cov[nc];*/
-				lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
-				/* 	 printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-			}
-			ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-			/* 	printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-		}
-		for(j=i+1; j<=nlstate+ndeath;j++){
-			for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-				/*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
-				lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
-				/* 	  printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
-			}
-			ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-		}
-	}
-	
-	for(i=1; i<= nlstate; i++){
-		s1=0;
-		for(j=1; j<i; j++){
-			s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-			/*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-		}
-		for(j=i+1; j<=nlstate+ndeath; j++){
-			s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-			/*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-		}
-		/* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
-		ps[i][i]=1./(s1+1.);
-		/* Computing other pijs */
-		for(j=1; j<i; j++)
-			ps[i][j]= exp(ps[i][j])*ps[i][i];
-		for(j=i+1; j<=nlstate+ndeath; j++)
-			ps[i][j]= exp(ps[i][j])*ps[i][i];
-		/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
-	} /* end i */
-	
-	for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
-		for(jj=1; jj<= nlstate+ndeath; jj++){
-			ps[ii][jj]=0;
-			ps[ii][ii]=1;
-		}
-	}
-	/* Added for backcast */ /* Transposed matrix too */
-	for(jj=1; jj<= nlstate+ndeath; jj++){
-		s1=0.;
-		for(ii=1; ii<= nlstate+ndeath; ii++){
-			s1+=ps[ii][jj];
-		}
-		for(ii=1; ii<= nlstate; ii++){
-			ps[ii][jj]=ps[ii][jj]/s1;
-		}
-	}
-	/* Transposition */
-	for(jj=1; jj<= nlstate+ndeath; jj++){
-		for(ii=jj; ii<= nlstate+ndeath; ii++){
-			s1=ps[ii][jj];
-			ps[ii][jj]=ps[jj][ii];
-			ps[jj][ii]=s1;
-		}
-	}
-	/* for(ii=1; ii<= nlstate+ndeath; ii++){ */
-	/*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
-	/* 	printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
-	/*   } */
-	/*   printf("\n "); */
-	/* } */
-	/* printf("\n ");printf("%lf ",cov[2]);*/
-	/*
-		for(i=1; i<= npar; i++) printf("%f ",x[i]);
-		goto end;*/
-	return ps;
+  for(i=1; i<= nlstate; i++){
+    for(j=1; j<i;j++){
+      for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	/*lnpijopii += param[i][j][nc]*cov[nc];*/
+	lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+	/* 	 printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+      }
+      ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+      /* 	printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+    }
+    for(j=i+1; j<=nlstate+ndeath;j++){
+      for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	/*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+	lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+	/* 	  printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+      }
+      ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+    }
+  }
+  
+  for(i=1; i<= nlstate; i++){
+    s1=0;
+    for(j=1; j<i; j++){
+      s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+      /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+    }
+    for(j=i+1; j<=nlstate+ndeath; j++){
+      s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+      /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+    }
+    /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+    ps[i][i]=1./(s1+1.);
+    /* Computing other pijs */
+    for(j=1; j<i; j++)
+      ps[i][j]= exp(ps[i][j])*ps[i][i];
+    for(j=i+1; j<=nlstate+ndeath; j++)
+      ps[i][j]= exp(ps[i][j])*ps[i][i];
+    /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+  } /* end i */
+  
+  for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+    for(jj=1; jj<= nlstate+ndeath; jj++){
+      ps[ii][jj]=0;
+      ps[ii][ii]=1;
+    }
+  }
+  /* Added for backcast */ /* Transposed matrix too */
+  for(jj=1; jj<= nlstate+ndeath; jj++){
+    s1=0.;
+    for(ii=1; ii<= nlstate+ndeath; ii++){
+      s1+=ps[ii][jj];
+    }
+    for(ii=1; ii<= nlstate; ii++){
+      ps[ii][jj]=ps[ii][jj]/s1;
+    }
+  }
+  /* Transposition */
+  for(jj=1; jj<= nlstate+ndeath; jj++){
+    for(ii=jj; ii<= nlstate+ndeath; ii++){
+      s1=ps[ii][jj];
+      ps[ii][jj]=ps[jj][ii];
+      ps[jj][ii]=s1;
+    }
+  }
+  /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+  /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
+  /* 	printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+  /*   } */
+  /*   printf("\n "); */
+  /* } */
+  /* printf("\n ");printf("%lf ",cov[2]);*/
+  /*
+    for(i=1; i<= npar; i++) printf("%f ",x[i]);
+    goto end;*/
+  return ps;
 }
 
 
@@ -2815,7 +2930,7 @@ double **matprod2(double **out, double *
 
 /************* Higher Matrix Product ***************/
 
-double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
+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 
      'nhstepm*hstepm*stepm' months (i.e. until
@@ -2851,16 +2966,34 @@ double ***hpxij(double ***po, int nhstep
       cov[2]=agexact;
       if(nagesqr==1)
 	cov[3]= agexact*agexact;
-      for (k=1; k<=cptcovn;k++) 
-	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
-      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
-      for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
-	/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	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,Tvar[Tage[k]])]*cov[2]; */
-      for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
-      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+      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,k)];
+	/* printf("hpxij 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]; 
+	/* 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++){
+	if(Dummy[Tvar[Tage[k]]]){
+	  cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+	} else{
+	  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]); */
+      }
+      for (k=1; k<=cptcovprod;k++){ /*  */
+	/* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
+	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+      }
+      /* 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 *\/ */
+      /* 	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */
+      /* for (k=1; k<=cptcovprod;k++) /\* Useless because included in cptcovn *\/ */
+      /* 	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */
       
       
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -3041,9 +3174,9 @@ double func( double *x)
       */
       ioffset=2+nagesqr+cptcovage;
    /* Fixed */
-			for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
-				cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
-			}
+      for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
+	cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
+      }
       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
 	 is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] 
 	 has been calculated etc */
@@ -3057,85 +3190,85 @@ double func( double *x)
 	 meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
       */
       for(mi=1; mi<= wav[i]-1; mi++){
-				for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/
-					cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];
-				}
-				for (ii=1;ii<=nlstate+ndeath;ii++)
-					for (j=1;j<=nlstate+ndeath;j++){
-						oldm[ii][j]=(ii==j ? 1.0 : 0.0);
-						savm[ii][j]=(ii==j ? 1.0 : 0.0);
-					}
-				for(d=0; d<dh[mi][i]; d++){
-					newm=savm;
-					agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
-					cov[2]=agexact;
-					if(nagesqr==1)
-						cov[3]= agexact*agexact;  /* Should be changed here */
-					for (kk=1; kk<=cptcovage;kk++) {
-						cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
-					}
-					out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
-											 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
-					savm=oldm;
-					oldm=newm;
-				} /* end mult */
-				
-				/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
-				/* But now since version 0.9 we anticipate for bias at large stepm.
-				 * If stepm is larger than one month (smallest stepm) and if the exact delay 
-				 * (in months) between two waves is not a multiple of stepm, we rounded to 
-				 * the nearest (and in case of equal distance, to the lowest) interval but now
-				 * we keep into memory the bias bh[mi][i] and also the previous matrix product
-				 * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
-				 * probability in order to take into account the bias as a fraction of the way
+	for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/
+	  cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];
+	}
+	for (ii=1;ii<=nlstate+ndeath;ii++)
+	  for (j=1;j<=nlstate+ndeath;j++){
+	    oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	    savm[ii][j]=(ii==j ? 1.0 : 0.0);
+	  }
+	for(d=0; d<dh[mi][i]; d++){
+	  newm=savm;
+	  agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+	  cov[2]=agexact;
+	  if(nagesqr==1)
+	    cov[3]= agexact*agexact;  /* Should be changed here */
+	  for (kk=1; kk<=cptcovage;kk++) {
+	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+	  }
+	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+		       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+	  savm=oldm;
+	  oldm=newm;
+	} /* end mult */
+	
+	/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+	/* But now since version 0.9 we anticipate for bias at large stepm.
+	 * If stepm is larger than one month (smallest stepm) and if the exact delay 
+	 * (in months) between two waves is not a multiple of stepm, we rounded to 
+	 * the nearest (and in case of equal distance, to the lowest) interval but now
+	 * we keep into memory the bias bh[mi][i] and also the previous matrix product
+	 * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
+	 * probability in order to take into account the bias as a fraction of the way
 				 * from savm to out if bh is negative or even beyond if bh is positive. bh varies
 				 * -stepm/2 to stepm/2 .
 				 * For stepm=1 the results are the same as for previous versions of Imach.
 				 * For stepm > 1 the results are less biased than in previous versions. 
 				 */
-				s1=s[mw[mi][i]][i];
-				s2=s[mw[mi+1][i]][i];
-				bbh=(double)bh[mi][i]/(double)stepm; 
-				/* bias bh is positive if real duration
-				 * is higher than the multiple of stepm and negative otherwise.
-				 */
-				/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
-				if( s2 > nlstate){ 
-					/* i.e. if s2 is a death state and if the date of death is known 
-						 then the contribution to the likelihood is the probability to 
-						 die between last step unit time and current  step unit time, 
-						 which is also equal to probability to die before dh 
-						 minus probability to die before dh-stepm . 
-						 In version up to 0.92 likelihood was computed
-						 as if date of death was unknown. Death was treated as any other
-						 health state: the date of the interview describes the actual state
-						 and not the date of a change in health state. The former idea was
-						 to consider that at each interview the state was recorded
-						 (healthy, disable or death) and IMaCh was corrected; but when we
-						 introduced the exact date of death then we should have modified
-						 the contribution of an exact death to the likelihood. This new
-						 contribution is smaller and very dependent of the step unit
-						 stepm. It is no more the probability to die between last interview
-						 and month of death but the probability to survive from last
-						 interview up to one month before death multiplied by the
-						 probability to die within a month. Thanks to Chris
-						 Jackson for correcting this bug.  Former versions increased
-						 mortality artificially. The bad side is that we add another loop
-						 which slows down the processing. The difference can be up to 10%
-						 lower mortality.
-					*/
-					/* If, at the beginning of the maximization mostly, the
-						 cumulative probability or probability to be dead is
-						 constant (ie = 1) over time d, the difference is equal to
-						 0.  out[s1][3] = savm[s1][3]: probability, being at state
-						 s1 at precedent wave, to be dead a month before current
-						 wave is equal to probability, being at state s1 at
-						 precedent wave, to be dead at mont of the current
-						 wave. Then the observed probability (that this person died)
-						 is null according to current estimated parameter. In fact,
-						 it should be very low but not zero otherwise the log go to
-						 infinity.
-					*/
+	s1=s[mw[mi][i]][i];
+	s2=s[mw[mi+1][i]][i];
+	bbh=(double)bh[mi][i]/(double)stepm; 
+	/* bias bh is positive if real duration
+	 * is higher than the multiple of stepm and negative otherwise.
+	 */
+	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+	if( s2 > nlstate){ 
+	  /* i.e. if s2 is a death state and if the date of death is known 
+	     then the contribution to the likelihood is the probability to 
+	     die between last step unit time and current  step unit time, 
+	     which is also equal to probability to die before dh 
+	     minus probability to die before dh-stepm . 
+	     In version up to 0.92 likelihood was computed
+	     as if date of death was unknown. Death was treated as any other
+	     health state: the date of the interview describes the actual state
+	     and not the date of a change in health state. The former idea was
+	     to consider that at each interview the state was recorded
+	     (healthy, disable or death) and IMaCh was corrected; but when we
+	     introduced the exact date of death then we should have modified
+	     the contribution of an exact death to the likelihood. This new
+	     contribution is smaller and very dependent of the step unit
+	     stepm. It is no more the probability to die between last interview
+	     and month of death but the probability to survive from last
+	     interview up to one month before death multiplied by the
+	     probability to die within a month. Thanks to Chris
+	     Jackson for correcting this bug.  Former versions increased
+	     mortality artificially. The bad side is that we add another loop
+	     which slows down the processing. The difference can be up to 10%
+	     lower mortality.
+	  */
+	  /* If, at the beginning of the maximization mostly, the
+	     cumulative probability or probability to be dead is
+	     constant (ie = 1) over time d, the difference is equal to
+	     0.  out[s1][3] = savm[s1][3]: probability, being at state
+	     s1 at precedent wave, to be dead a month before current
+	     wave is equal to probability, being at state s1 at
+	     precedent wave, to be dead at mont of the current
+	     wave. Then the observed probability (that this person died)
+	     is null according to current estimated parameter. In fact,
+	     it should be very low but not zero otherwise the log go to
+	     infinity.
+	  */
 /* #ifdef INFINITYORIGINAL */
 /* 	    lli=log(out[s1][s2] - savm[s1][s2]); */
 /* #else */
@@ -4003,7 +4136,7 @@ void  freqsummary(char fileres[], int ia
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);
   }
-
+  
   strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
   if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
     printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
@@ -4013,55 +4146,54 @@ void  freqsummary(char fileres[], int ia
   }
   else{
     fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
-<hr size=\"2\" color=\"#EC5E5E\"> \n\
+<hr size=\"2\" color=\"#EC5E5E\"> \n					\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
 	    fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
-  fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
-    
+  fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);
+  
   strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
     printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
     fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
     fflush(ficlog);
     exit(70); 
-  }
-  else{
+  } else{
     fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
-<hr size=\"2\" color=\"#EC5E5E\"> \n\
+<hr size=\"2\" color=\"#EC5E5E\"> \n					\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
 	    fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
-  fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
-
+  fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
+  
   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
   j1=0;
   
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
   j=cptcoveff;  /* Only dummy covariates of the model */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
-
+  
   first=1;
-
+  
   /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
      reference=low_education V1=0,V2=0
      med_educ                V1=1 V2=0, 
      high_educ               V1=0 V2=1
      Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff 
   */
-
+  
   for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
     posproptt=0.;
     /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
       scanf("%d", i);*/
     for (i=-5; i<=nlstate+ndeath; i++)  
       for (jk=-5; jk<=nlstate+ndeath; jk++)  
-				for(m=iagemin; m <= iagemax+3; m++)
-					freq[i][jk][m]=0;
-		
+	for(m=iagemin; m <= iagemax+3; m++)
+	  freq[i][jk][m]=0;
+    
     for (i=1; i<=nlstate; i++)  {
       for(m=iagemin; m <= iagemax+3; m++)
-				prop[i][m]=0;
+	prop[i][m]=0;
       posprop[i]=0;
       pospropt[i]=0;
     }
@@ -4071,233 +4203,239 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     /* 	meanqt[m][z1]=0.; */
     /*   } */
     /* } */
-		
+    
     dateintsum=0;
     k2cpt=0;
     /* For that combination of covariate j1, we count and print the frequencies in one pass */
     for (iind=1; iind<=imx; iind++) { /* For each individual iind */
       bool=1;
       if(anyvaryingduminmodel==0){ /* If All fixed covariates */
-				if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+	if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
 	  /* for (z1=1; z1<= nqfveff; z1++) {   */
 	  /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */
 	  /* } */
-					for (z1=1; z1<=cptcoveff; z1++) {  
-						/* if(Tvaraff[z1] ==-20){ */
-						/* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
-						/* }else  if(Tvaraff[z1] ==-10){ */
-						/* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
-						/* }else  */
-						if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
-							/* Tests if this individual iind responded to j1 (V4=1 V3=0) */
-							bool=0;
-							/* 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 */
-				} /* cptcovn > 0 */
+	  for (z1=1; z1<=cptcoveff; z1++) {  
+	    /* if(Tvaraff[z1] ==-20){ */
+	    /* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+	    /* }else  if(Tvaraff[z1] ==-10){ */
+	    /* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
+	    /* }else  */
+	    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
+	      /* Tests if this individual iind responded to j1 (V4=1 V3=0) */
+	      bool=0;
+	      /* 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 */
+	} /* cptcovn > 0 */
       } /* end any */
       if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
-				/* for(m=firstpass; m<=lastpass; m++){ */
-				for(mi=1; mi<wav[iind];mi++){ /* For that wave */
-					m=mw[mi][iind];
-					if(anyvaryingduminmodel==1){ /* Some are varying covariates */
-						for (z1=1; z1<=cptcoveff; z1++) {
-							if( Fixed[Tmodelind[z1]]==1){
-								iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
-								if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
-									bool=0;
-							}else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
-								if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
-									bool=0;
-								}
-							}
-						}
-					}/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
-					/* bool =0 we keep that guy which corresponds to the combination of dummy values */
-					if(bool==1){
-						/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
-							 and mw[mi+1][iind]. dh depends on stepm. */
-						agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
-						ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
-						if(m >=firstpass && m <=lastpass){
-							k2=anint[m][iind]+(mint[m][iind]/12.);
-							/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
-							if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
-							if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
-							if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
-								prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
-							if (m<lastpass) {
-								/* if(s[m][iind]==4 && s[m+1][iind]==4) */
-								/*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
-								if(s[m][iind]==-1)
-									printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
-								freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
-								/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
-								freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
-							}
-						} /* end if between passes */  
-						if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
-							dateintsum=dateintsum+k2;
-							k2cpt++;
-							/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
-						}
-					} /* end bool 2 */
-				} /* end m */
+	/* for(m=firstpass; m<=lastpass; m++){ */
+	for(mi=1; mi<wav[iind];mi++){ /* For that wave */
+	  m=mw[mi][iind];
+	  if(anyvaryingduminmodel==1){ /* Some are varying covariates */
+	    for (z1=1; z1<=cptcoveff; z1++) {
+	      if( Fixed[Tmodelind[z1]]==1){
+		iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+		if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
+		  bool=0;
+	      }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
+		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
+		  bool=0;
+		}
+	      }
+	    }
+	  }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
+	  /* bool =0 we keep that guy which corresponds to the combination of dummy values */
+	  if(bool==1){
+	    /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
+	       and mw[mi+1][iind]. dh depends on stepm. */
+	    agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
+	    ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
+	    if(m >=firstpass && m <=lastpass){
+	      k2=anint[m][iind]+(mint[m][iind]/12.);
+	      /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+	      if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
+	      if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
+	      if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
+		prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
+	      if (m<lastpass) {
+		/* if(s[m][iind]==4 && s[m+1][iind]==4) */
+		/*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
+		if(s[m][iind]==-1)
+		  printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
+		freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+		/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
+		freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
+	      }
+	    } /* end if between passes */  
+	    if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
+	      dateintsum=dateintsum+k2;
+	      k2cpt++;
+	      /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
+	    }
+	  } /* end bool 2 */
+	} /* end m */
       } /* end bool */
     } /* end iind = 1 to imx */
     /* prop[s][age] is feeded for any initial and valid live state as well as
        freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
-		
-		
+    
+    
     /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
     pstamp(ficresp);
-    /* if  (ncoveff>0) { */
-    if  (cptcoveff>0) {
+    if  (cptcoveff>0){
       fprintf(ficresp, "\n#********** Variable "); 
       fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
       fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
+      fprintf(ficlog, "\n#********** Variable "); 
       for (z1=1; z1<=cptcoveff; z1++){
-				fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-				fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-				fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	if(DummyV[z1]){
+	  fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	}else{
+	  fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	}
       }
       fprintf(ficresp, "**********\n#");
       fprintf(ficresphtm, "**********</h3>\n");
       fprintf(ficresphtmfr, "**********</h3>\n");
-      fprintf(ficlog, "\n#********** Variable "); 
-      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
       fprintf(ficlog, "**********\n");
     }
     fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
     for(i=1; i<=nlstate;i++) {
-      fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+      fprintf(ficresp, " Age Prev(%d)  N(%d)  N  ",i,i);
       fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
     }
     fprintf(ficresp, "\n");
     fprintf(ficresphtm, "\n");
-		
+    
     /* Header of frequency table by age */
     fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
     fprintf(ficresphtmfr,"<th>Age</th> ");
     for(jk=-1; jk <=nlstate+ndeath; jk++){
       for(m=-1; m <=nlstate+ndeath; m++){
-				if(jk!=0 && m!=0)
-					fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
+	if(jk!=0 && m!=0)
+	  fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
       }
     }
     fprintf(ficresphtmfr, "\n");
-		
+    
     /* For each age */
     for(iage=iagemin; iage <= iagemax+3; iage++){
       fprintf(ficresphtm,"<tr>");
       if(iage==iagemax+1){
-				fprintf(ficlog,"1");
-				fprintf(ficresphtmfr,"<tr><th>0</th> ");
+	fprintf(ficlog,"1");
+	fprintf(ficresphtmfr,"<tr><th>0</th> ");
       }else if(iage==iagemax+2){
-				fprintf(ficlog,"0");
-				fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+	fprintf(ficlog,"0");
+	fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
       }else if(iage==iagemax+3){
-				fprintf(ficlog,"Total");
-				fprintf(ficresphtmfr,"<tr><th>Total</th> ");
+	fprintf(ficlog,"Total");
+	fprintf(ficresphtmfr,"<tr><th>Total</th> ");
       }else{
-				if(first==1){
-					first=0;
-					printf("See log file for details...\n");
-				}
-				fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
-				fprintf(ficlog,"Age %d", iage);
+	if(first==1){
+	  first=0;
+	  printf("See log file for details...\n");
+	}
+	fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
+	fprintf(ficlog,"Age %d", iage);
       }
       for(jk=1; jk <=nlstate ; jk++){
-				for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-					pp[jk] += freq[jk][m][iage]; 
+	for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+	  pp[jk] += freq[jk][m][iage]; 
       }
       for(jk=1; jk <=nlstate ; jk++){
-				for(m=-1, pos=0; m <=0 ; m++)
-					pos += freq[jk][m][iage];
-				if(pp[jk]>=1.e-10){
-					if(first==1){
-						printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
-					}
-					fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
-				}else{
-					if(first==1)
-						printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
-					fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
-				}
+	for(m=-1, pos=0; m <=0 ; m++)
+	  pos += freq[jk][m][iage];
+	if(pp[jk]>=1.e-10){
+	  if(first==1){
+	    printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+	  }
+	  fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+	}else{
+	  if(first==1)
+	    printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+	  fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+	}
       }
-			
+      
       for(jk=1; jk <=nlstate ; jk++){ 
-				/* posprop[jk]=0; */
-				for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
-					pp[jk] += freq[jk][m][iage];
+	/* posprop[jk]=0; */
+	for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
+	  pp[jk] += freq[jk][m][iage];
       }	/* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
-			
+      
       for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
-				pos += pp[jk]; /* pos is the total number of transitions until this age */
-				posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
-																					from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
-				pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
-																				from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+	pos += pp[jk]; /* pos is the total number of transitions until this age */
+	posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
+					  from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+	pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
+					from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
       }
       for(jk=1; jk <=nlstate ; jk++){
-				if(pos>=1.e-5){
-					if(first==1)
-						printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-					fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-				}else{
-					if(first==1)
-						printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
-					fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
-				}
-				if( iage <= iagemax){
-					if(pos>=1.e-5){
-						fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-						fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-						/*probs[iage][jk][j1]= pp[jk]/pos;*/
-						/*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
-					}
-					else{
-						fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
-						fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
-					}
-				}
-				pospropt[jk] +=posprop[jk];
+	if(pos>=1.e-5){
+	  if(first==1)
+	    printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+	  fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+	}else{
+	  if(first==1)
+	    printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+	  fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+	}
+	if( iage <= iagemax){
+	  if(pos>=1.e-5){
+	    fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+	    fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+	    /*probs[iage][jk][j1]= pp[jk]/pos;*/
+	    /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
+	  }
+	  else{
+	    fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
+	    fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
+	  }
+	}
+	pospropt[jk] +=posprop[jk];
       } /* end loop jk */
       /* pospropt=0.; */
       for(jk=-1; jk <=nlstate+ndeath; jk++){
-				for(m=-1; m <=nlstate+ndeath; m++){
-					if(freq[jk][m][iage] !=0 ) { /* minimizing output */
-						if(first==1){
-							printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
-						}
-						fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
-					}
-					if(jk!=0 && m!=0)
-						fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
-				}
+	for(m=-1; m <=nlstate+ndeath; m++){
+	  if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+	    if(first==1){
+	      printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	    }
+	    fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	  }
+	  if(jk!=0 && m!=0)
+	    fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
+	}
       } /* end loop jk */
       posproptt=0.; 
       for(jk=1; jk <=nlstate; jk++){
-				posproptt += pospropt[jk];
+	posproptt += pospropt[jk];
       }
       fprintf(ficresphtmfr,"</tr>\n ");
       if(iage <= iagemax){
-				fprintf(ficresp,"\n");
-				fprintf(ficresphtm,"</tr>\n");
+	fprintf(ficresp,"\n");
+	fprintf(ficresphtm,"</tr>\n");
       }
       if(first==1)
-				printf("Others in log...\n");
+	printf("Others in log...\n");
       fprintf(ficlog,"\n");
     } /* end loop age iage */
     fprintf(ficresphtm,"<tr><th>Tot</th>");
     for(jk=1; jk <=nlstate ; jk++){
       if(posproptt < 1.e-5){
-				fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);	
+	fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);	
       }else{
-				fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	
+	fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	
       }
     }
     fprintf(ficresphtm,"</tr>\n");
@@ -4315,7 +4453,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     fprintf(ficresphtmfr,"</table>\n");
   } /* end selected combination of covariate j1 */
   dateintmean=dateintsum/k2cpt; 
-	
+  
   fclose(ficresp);
   fclose(ficresphtm);
   fclose(ficresphtmfr);
@@ -4821,7 +4959,7 @@ void  concatwav(int wav[], int **dh, int
 
 /*********** Health Expectancies ****************/
 
-void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
+ void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[], int nres )
 
 {
   /* Health expectancies, no variances */
@@ -4831,7 +4969,7 @@ void evsij(double ***eij, double x[], in
   double ***p3mat;
   double eip;
 
-  pstamp(ficreseij);
+  /* pstamp(ficreseij); */
   fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
   fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++){
@@ -4894,7 +5032,7 @@ void evsij(double ***eij, double x[], in
     /* Computed by stepm unit matrices, product of hstepma matrices, stored
        in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */
     
-    hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
+    hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij, nres);  
     
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
     
@@ -4929,7 +5067,7 @@ void evsij(double ***eij, double x[], in
   
 }
 
-void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
+ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[], int nres )
 
 {
   /* Covariances of health expectancies eij and of total life expectancies according
@@ -5042,8 +5180,8 @@ void cvevsij(double ***eij, double x[],
 	xp[i] = x[i] + (i==theta ?delti[theta]:0);
 	xm[i] = x[i] - (i==theta ?delti[theta]:0);
       }
-      hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);  
-      hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);  
+      hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij, nres);  
+      hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij, nres);  
 			
       for(j=1; j<= nlstate; j++){
 	for(i=1; i<=nlstate; i++){
@@ -5084,7 +5222,7 @@ void cvevsij(double ***eij, double x[],
     }
 		
     /* Computing expectancies */
-    hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
+    hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);  
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)
 	for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
@@ -5139,7 +5277,7 @@ void cvevsij(double ***eij, double x[],
 }
  
 /************ Variance ******************/
- void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
+ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
  {
    /* Variance of health expectancies */
    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
@@ -5196,6 +5334,14 @@ void cvevsij(double ***eij, double x[],
    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
    pstamp(ficresprobmorprev);
    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
+   fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");
+   for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+     fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+   }
+   for(j=1;j<=cptcoveff;j++) 
+     fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,j)]);
+   fprintf(ficresprobmorprev,"\n");
+
    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
    for(j=nlstate+1; j<=(nlstate+ndeath);j++){
      fprintf(ficresprobmorprev," p.%-d SE",j);
@@ -5265,7 +5411,7 @@ void cvevsij(double ***eij, double x[],
 	 xp[i] = x[i] + (i==theta ?delti[theta]:0);
        }
 			
-       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);
 			
        if (popbased==1) {
 	 if(mobilav ==0){
@@ -5277,7 +5423,7 @@ void cvevsij(double ***eij, double x[],
 	 }
        }
 			
-       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
+       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
        for(j=1; j<= nlstate; j++){
 	 for(h=0; h<=nhstepm; h++){
 	   for(i=1, gp[h][j]=0.;i<=nlstate;i++)
@@ -5297,7 +5443,7 @@ void cvevsij(double ***eij, double x[],
        for(i=1; i<=npar; i++) /* Computes gradient x - delta */
 	 xp[i] = x[i] - (i==theta ?delti[theta]:0);
 			
-       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
+       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult);
 			
        if (popbased==1) {
 	 if(mobilav ==0){
@@ -5309,7 +5455,7 @@ void cvevsij(double ***eij, double x[],
 	 }
        }
 			
-       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
+       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  
 			
        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
 	 for(h=0; h<=nhstepm; h++){
@@ -5374,7 +5520,7 @@ void cvevsij(double ***eij, double x[],
      /* end ppptj */
      /*  x centered again */
 		
-     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
+     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);
 		
      if (popbased==1) {
        if(mobilav ==0){
@@ -5390,7 +5536,7 @@ void cvevsij(double ***eij, double x[],
 	computed over hstepm (estepm) matrices product = hstepm*stepm months) 
 	as a weighted average of prlim.
      */
-     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);  
+     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);  
      for(j=nlstate+1;j<=nlstate+ndeath;j++){
        for(i=1,gmp[j]=0.;i<= nlstate; i++) 
 	 gmp[j] += prlim[i][i]*p3mat[i][j][1]; 
@@ -5453,7 +5599,7 @@ void cvevsij(double ***eij, double x[],
  }  /* end varevsij */
 
 /************ Variance of prevlim ******************/
- void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
+ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres)
 {
   /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
@@ -5496,9 +5642,9 @@ void cvevsij(double ***eij, double x[],
 	xp[i] = x[i] + (i==theta ?delti[theta]:0);
       }
       if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       else
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       for(i=1;i<=nlstate;i++){
 	gp[i] = prlim[i][i];
 	mgp[theta][i] = prlim[i][i];
@@ -5506,9 +5652,9 @@ void cvevsij(double ***eij, double x[],
       for(i=1; i<=npar; i++) /* Computes gradient */
 	xp[i] = x[i] - (i==theta ?delti[theta]:0);
       if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       else
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       for(i=1;i<=nlstate;i++){
 	gm[i] = prlim[i][i];
 	mgm[theta][i] = prlim[i][i];
@@ -5916,11 +6062,13 @@ void printinghtml(char fileresu[], char
 		  int popforecast, int prevfcast, int backcast, int estepm , \
 		  double jprev1, double mprev1,double anprev1, double dateprev1, \
 		  double jprev2, double mprev2,double anprev2, double dateprev2){
-  int jj1, k1, i1, cpt;
+  int jj1, k1, i1, cpt, k4, nres;
 
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
    <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
 </ul>");
+   fprintf(fichtm,"<ul><li> model=1+age+%s\n \
+</ul>", model);
    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
    fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
 	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
@@ -5955,16 +6103,27 @@ void printinghtml(char fileresu[], char
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
    jj1=0;
+
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){
+     if(TKresult[nres]!= k1)
+       continue;
 
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
-	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
+	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+	 printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout);
+	 /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+	 /* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
        }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);fflush(stdout);
+      }
+       
        /* if(nqfveff+nqtveff 0) */ /* Test to be done */
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        if(invalidvarcomb[k1]){
@@ -6076,13 +6235,22 @@ See page 'Matrix of variance-covariance
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
    jj1=0;
+
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){
+     if(TKresult[nres]!= k1)
+       continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++)  /**< cptcoveff number of variables */
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);
+	 /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
+
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
 
        if(invalidvarcomb[k1]){
@@ -6112,11 +6280,12 @@ void printinggnuplot(char fileresu[], ch
 
   char dirfileres[132],optfileres[132];
   char gplotcondition[132];
-  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
   int ioffset; /* variable offset for columns */
+  int nres=0; /* Index of resultline */
 
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
@@ -6160,135 +6329,95 @@ void printinggnuplot(char fileresu[], ch
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
   /* 1eme*/
-  for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
-    for (k1=1; k1<= m && selected(k1) ; k1 ++) { /* For each valid combination of covariate */
-      /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-      fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
-	/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-	continue;
-      }
-
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
-      fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \n\
-set ylabel \"Probability\" \n	\
-set ter svg size 640, 480\n	\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-			
-      for (i=1; i<= nlstate ; i ++) {
-	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-	else        fprintf(ficgp," %%*lf (%%*lf)");
-      }
-      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-      for (i=1; i<= nlstate ; i ++) {
-	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      } 
-      fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
-      for (i=1; i<= nlstate ; i ++) {
-	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }  
-      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
-      if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
-	/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
-	fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
-	if(cptcoveff ==0){
-	  fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",	 2+(cpt-1),  cpt );
-	}else{
-	  kl=0;
-	  for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
-	    lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
-	    /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	    /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	    /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	    vlv= nbcode[Tvaraff[k]][lv];
-	    kl++;
-	    /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
-	    /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
-	    /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
-	    /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
-	    if(k==cptcoveff){
-	      fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
-		      4+(cpt-1),  cpt );  /* 4 or 6 ?*/
-	    }else{
-	      fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+  for (cpt=1; cpt<= nlstate ; cpt ++){ /* For each live state */
+    for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
+      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+	/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+	if(TKresult[nres]!= k1)
+	  continue;
+	/* We are interested in selected combination by the resultline */
+	printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
+	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
+	for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
+	  lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	  vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
+	  /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+	  printf(" V%d=%d ",Tvaraff[k],vlv);
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	printf("\n#\n");
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
+      
+	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+	fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+      
+	for (i=1; i<= nlstate ; i ++) {
+	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+	  else        fprintf(ficgp," %%*lf (%%*lf)");
+	}
+	fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+	for (i=1; i<= nlstate ; i ++) {
+	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	} 
+	fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
+	for (i=1; i<= nlstate ; i ++) {
+	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}  
+	fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+	if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
+	  /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
+	  fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
+	  if(cptcoveff ==0){
+	    fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",	 2+(cpt-1),  cpt );
+	  }else{
+	    kl=0;
+	    for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
+	      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	      vlv= nbcode[Tvaraff[k]][lv];
 	      kl++;
-	    }
-	  } /* end covariate */
-	} /* end if no covariate */
-      } /* end if backcast */
-      fprintf(ficgp,"\nset out \n");
+	      /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+	      /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
+	      /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
+	      /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+	      if(k==cptcoveff){
+		fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+			4+(cpt-1),  cpt );  /* 4 or 6 ?*/
+	      }else{
+		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+		kl++;
+	      }
+	    } /* end covariate */
+	  } /* end if no covariate */
+	} /* end if backcast */
+	fprintf(ficgp,"\nset out \n");
+      } /* nres */
     } /* k1 */
   } /* cpt */
-  /*2 eme*/
-  for (k1=1; k1<= m ; k1 ++) { 
-
-    fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
-    for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-      vlv= nbcode[Tvaraff[k]][lv];
-      fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-    }
-    fprintf(ficgp,"\n#\n");
-    if(invalidvarcomb[k1]){
-      fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-      continue;
-    }
-			
-    fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
-    for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-      if(vpopbased==0)
-	fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-      else
-	fprintf(ficgp,"\nreplot ");
-      for (i=1; i<= nlstate+1 ; i ++) {
-	k=2*i;
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-	else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	fprintf(ficgp,"\" t\"\" w l lt 0,");
-	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-	for (j=1; j<= nlstate+1 ; j ++) {
-	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	  else fprintf(ficgp," %%*lf (%%*lf)");
-	}   
-	if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-	else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
-      } /* state */
-    } /* vpopbased */
-    fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
-  } /* k1 */
-	
-	
-  /*3eme*/
-  for (k1=1; k1<= m ; k1 ++) { 
 
-    for (cpt=1; cpt<= nlstate ; cpt ++) {
-      fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  cov=%d state=%d",k1, cpt);
+  
+  /*2 eme*/
+  for (k1=1; k1<= m ; k1 ++){  
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
+	continue;
+      fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -6297,129 +6426,206 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
+      /* for(k=1; k <= ncovds; k++){ */
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
       fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){
 	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
 	continue;
       }
 			
-      /*       k=2+nlstate*(2*cpt-2); */
-      k=2+(nlstate+1)*(cpt-1);
-      fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
-      fprintf(ficgp,"set ter svg size 640, 480\n\
+      fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+      for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+	if(vpopbased==0)
+	  fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+	else
+	  fprintf(ficgp,"\nreplot ");
+	for (i=1; i<= nlstate+1 ; i ++) {
+	  k=2*i;
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+	  for (j=1; j<= nlstate+1 ; j ++) {
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	    else fprintf(ficgp," %%*lf (%%*lf)");
+	  }   
+	  if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+	  else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	  for (j=1; j<= nlstate+1 ; j ++) {
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	    else fprintf(ficgp," %%*lf (%%*lf)");
+	  }   
+	  fprintf(ficgp,"\" t\"\" w l lt 0,");
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	  for (j=1; j<= nlstate+1 ; j ++) {
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	    else fprintf(ficgp," %%*lf (%%*lf)");
+	  }   
+	  if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+	  else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+	} /* state */
+      } /* vpopbased */
+      fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+    } /* end nres */
+  } /* k1 end 2 eme*/
+	
+	
+  /*3eme*/
+  for (k1=1; k1<= m ; k1 ++){
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
+	continue;
+
+      for (cpt=1; cpt<= nlstate ; cpt ++) {
+	fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
+	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	  vlv= nbcode[Tvaraff[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
+			
+	/*       k=2+nlstate*(2*cpt-2); */
+	k=2+(nlstate+1)*(cpt-1);
+	fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+	fprintf(ficgp,"set ter svg size 640, 480\n\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
-      /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-	fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+	/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+	  for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+	  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+	  fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+	  for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+	  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
 				
-      */
-      for (i=1; i< nlstate ; i ++) {
-	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
-	/*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+	*/
+	for (i=1; i< nlstate ; i ++) {
+	  fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+	  /*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
 				
-      } 
-      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
-    }
-  }
+	} 
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
+      }
+    } /* end nres */
+  } /* end kl 3eme */
   
   /* 4eme */
   /* Survival functions (period) from state i in state j by initial state i */
-  for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
-
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv];
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+  for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
 	continue;
-      }
-			
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
-set ter svg size 640, 480\n																							\
-unset log y\n																														\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
-      k=3;
-      for (i=1; i<= nlstate ; i ++){
-	if(i==1){
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-	}else{
-	  fprintf(ficgp,", '' ");
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
+	fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
+	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	  vlv= nbcode[Tvaraff[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	}
-	l=(nlstate+ndeath)*(i-1)+1;
-	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
-	for (j=2; j<= nlstate+ndeath ; j ++)
-	  fprintf(ficgp,"+$%d",k+l+j-1);
-	fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
-      } /* nlstate */
-      fprintf(ficgp,"\nset out\n");
-    } /* end cpt state*/ 
-  } /* end covariate */  
-	
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
+      
+	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
+	k=3;
+	for (i=1; i<= nlstate ; i ++){
+	  if(i==1){
+	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	  }else{
+	    fprintf(ficgp,", '' ");
+	  }
+	  l=(nlstate+ndeath)*(i-1)+1;
+	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+	  for (j=2; j<= nlstate+ndeath ; j ++)
+	    fprintf(ficgp,"+$%d",k+l+j-1);
+	  fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+	} /* nlstate */
+	fprintf(ficgp,"\nset out\n");
+      } /* end cpt state*/ 
+    } /* end nres */
+  } /* end covariate k1 */  
+
 /* 5eme */
   /* Survival functions (period) from state i in state j by final state j */
-  for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-			
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-	vlv= nbcode[Tvaraff[k]][lv];
-	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      fprintf(ficgp,"\n#\n");
-      if(invalidvarcomb[k1]){
-	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+  for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
 	continue;
-      }
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+	fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	  /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	  vlv= nbcode[Tvaraff[k]][lv];
+	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
+	fprintf(ficgp,"\n#\n");
+	if(invalidvarcomb[k1]){
+	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	  continue;
+	}
       
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
-      k=3;
-      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-	if(j==1)
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-	else
-	  fprintf(ficgp,", '' ");
-	l=(nlstate+ndeath)*(cpt-1) +j;
-	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
-	/* for (i=2; i<= nlstate+ndeath ; i ++) */
-	/*   fprintf(ficgp,"+$%d",k+l+i-1); */
-	fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
-      } /* nlstate */
-      fprintf(ficgp,", '' ");
-      fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
-      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-	l=(nlstate+ndeath)*(cpt-1) +j;
-	if(j < nlstate)
-	  fprintf(ficgp,"$%d +",k+l);
-	else
-	  fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
-      }
-      fprintf(ficgp,"\nset out\n");
-    } /* end cpt state*/ 
-  } /* end covariate */  
+	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
+	k=3;
+	for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	  if(j==1)
+	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	  else
+	    fprintf(ficgp,", '' ");
+	  l=(nlstate+ndeath)*(cpt-1) +j;
+	  fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+	  /* for (i=2; i<= nlstate+ndeath ; i ++) */
+	  /*   fprintf(ficgp,"+$%d",k+l+i-1); */
+	  fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+	} /* nlstate */
+	fprintf(ficgp,", '' ");
+	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+	for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	  l=(nlstate+ndeath)*(cpt-1) +j;
+	  if(j < nlstate)
+	    fprintf(ficgp,"$%d +",k+l);
+	  else
+	    fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+	}
+	fprintf(ficgp,"\nset out\n");
+      } /* end cpt state*/ 
+    } /* end covariate */  
+  } /* end nres */
   
 /* 6eme */
   /* CV preval stable (period) for each covariate */
-  for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+  for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+    if(TKresult[nres]!= k1)
+      continue;
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
@@ -6431,6 +6637,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }	
       fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){
 	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6439,9 +6648,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
       
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){
 	if(i==1)
@@ -6462,7 +6669,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 /* 7eme */
   if(backcast == 1){
     /* CV back preval stable (period) for each covariate */
-    for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+    for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
+	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
 	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
@@ -6473,6 +6683,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6481,9 +6694,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	
 	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3; /* Offset */
 	for (i=1; i<= nlstate ; i ++){
 	  if(i==1)
@@ -6509,7 +6720,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   if(prevfcast==1){
     /* Projection from cross-sectional to stable (period) for each covariate */
     
-    for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+    for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k1)
+	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
 	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
@@ -6520,6 +6734,9 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	}
+	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}	
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6529,9 +6746,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
 	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n						\
-unset log y\n								\
-plot [%.f:%.f]  ", ageminpar, agemaxpar);
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
 	  /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
 	  /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */   
@@ -6598,8 +6813,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* End if prevfcast */
   
   
-  /* proba elementaires */
-  fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+  /* 9eme writing MLE parameters */
+  fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n");
   for(i=1,jk=1; i <=nlstate; i++){
     fprintf(ficgp,"# initial state %d\n",i);
     for(k=1; k <=(nlstate+ndeath); k++){
@@ -6616,7 +6831,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"##############\n#\n");
   
   /*goto avoid;*/
-  fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
+  /* 10eme Graphics of probabilities or incidences using written MLE parameters */
+  fprintf(ficgp,"\n##############\n#10eme Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
   fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
   fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -6631,10 +6847,19 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
   fprintf(ficgp,"#\n");
   for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
-    fprintf(ficgp,"# ng=%d\n",ng);
-    fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);
-    for(jk=1; jk <=m; jk++) {
-      fprintf(ficgp,"#    jk=%d\n",jk);
+    fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
+    fprintf(ficgp,"#model=%s \n",model);
+    fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
+    fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
+    for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= jk)
+	continue;
+      fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }	
+      fprintf(ficgp,"\n#\n");
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
       fprintf(ficgp,"\nset ter svg size 640, 480 ");
       if (ng==1){
@@ -6675,18 +6900,50 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	      break;
 	    }
 	    ij=1;/* To be checked else nbcode[0][0] wrong */
-	    for(j=3; j <=ncovmodel-nagesqr; j++) {
+	    ijp=1; /* product no age */
+	    /* for(j=3; j <=ncovmodel-nagesqr; j++) { */
+	    for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
 	      /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
-	      if(ij <=cptcovage) { /* Bug valgrind */
-		if((j-2)==Tage[ij]) { /* Bug valgrind */
-		  fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
-		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+	      if(j==Tage[ij]) { /* Product by age */
+		if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+		  if(DummyV[j]==0){
+		    fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
+		  }else{ /* quantitative */
+		    fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
+		    /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		  }
 		  ij++;
 		}
-	      }
-	      else
-		fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */
-	    }
+	      }else if(j==Tprod[ijp]) { /* */ 
+		/* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+		if(ijp <=cptcovprod) { /* Product */
+		  if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+		    if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+		      /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],nbcode[Tvard[ijp][2]][codtabm(jk,j)]); */
+		      fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
+		    }else{ /* Vn is dummy and Vm is quanti */
+		      /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+		      fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+		    }
+		  }else{ /* Vn*Vm Vn is quanti */
+		    if(DummyV[Tvard[ijp][2]]==0){
+		      fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
+		    }else{ /* Both quanti */
+		      fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+		    }
+		  }
+		  ijp++;
+  		}
+	      } else{  /* simple covariate */
+		/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
+		if(Dummy[j]==0){
+		  fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
+		}else{ /* quantitative */
+		  fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* */
+		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		}
+	      } /* end simple */
+	    } /* end j */
 	  }else{
 	    i=i-ncovmodel;
 	    if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
@@ -6704,8 +6961,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	       
 	      ij=1;
 	      for(j=3; j <=ncovmodel-nagesqr; j++){
-		if(ij <=cptcovage) { /* Bug valgrind */
-		  if((j-2)==Tage[ij]) { /* Bug valgrind */
+		if((j-2)==Tage[ij]) { /* Bug valgrind */
+		  if(ij <=cptcovage) { /* Bug valgrind */
 		    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 		    /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
 		    ij++;
@@ -6857,7 +7114,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
      } /* end bad */
 		
      for (age=bage; age<=fage; age++){
-       printf("%d %d ", cptcod, (int)age);
+       /* printf("%d %d ", cptcod, (int)age); */
        sumnewp[cptcod]=0.;
        sumnewm[cptcod]=0.;
        for (i=1; i<=nlstate;i++){
@@ -6896,13 +7153,13 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
  
 
 /************** Forecasting ******************/
-void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+ void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection 
      agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed
      anproj2 year of en of projection (same day and month as proj1).
   */
-  int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
+   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
   double agec; /* generic age */
   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;
@@ -6924,8 +7181,8 @@ void prevforecast(char fileres[], double
     printf("Problem with forecast resultfile: %s\n", fileresf);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
   }
-  printf("Computing forecasting: result on file '%s', please wait... \n", fileresf);
-  fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf);
+  printf("\nComputing forecasting: result on file '%s', please wait... \n", fileresf);
+  fprintf(ficlog,"\nComputing forecasting: result on file '%s', please wait... \n", fileresf);
 
   if (cptcoveff==0) ncodemax[cptcoveff]=1;
 
@@ -6956,7 +7213,10 @@ void prevforecast(char fileres[], double
   fprintf(ficresf,"#****** Routine prevforecast **\n");
   
 /* 	      if (h==(int)(YEARM*yearp)){ */
-  for(k=1;k<=i1;k++){
+  for(nres=1; nres <= nresult; nres++) /* For each resultline */
+  for(k=1; k<=i1;k++){
+    if(TKresult[nres]!= k)
+      continue;
     if(invalidvarcomb[k]){
       printf("\nCombination (%d) projection ignored because no cases \n",k); 
       continue;
@@ -6965,6 +7225,9 @@ void prevforecast(char fileres[], double
     for(j=1;j<=cptcoveff;j++) {
       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
     }
+    for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+      fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+    }
     fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){ 
       for(i=1; i<=nlstate;i++) 	      
@@ -6979,7 +7242,7 @@ void prevforecast(char fileres[], double
 	nhstepm = nhstepm/hstepm; 
 	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
 	oldm=oldms;savm=savms;
-	hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+	hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k,nres);
 	
 	for (h=0; h<=nhstepm; h++){
 	  if (h*hstepm/YEARM*stepm ==yearp) {
@@ -7572,14 +7835,36 @@ int readdata(char datafile[], int firsto
   /*-------- data file ----------*/
   FILE *fic;
   char dummy[]="                         ";
-  int i=0, j=0, n=0, iv=0;
+  int i=0, j=0, n=0, iv=0, v;
   int lstra;
   int linei, month, year,iout;
   char line[MAXLINE], linetmp[MAXLINE];
   char stra[MAXLINE], strb[MAXLINE];
   char *stratrunc;
 
+  DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
+  FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
 
+  for(v=1; v <=ncovcol;v++){
+    DummyV[v]=0;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+1; v <=ncovcol+nqv;v++){
+    DummyV[v]=1;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
+    DummyV[v]=0;
+    FixedV[v]=1;
+  }
+  for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    DummyV[v]=1;
+    FixedV[v]=1;
+  }
+  for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+    fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+  }
 
   if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
@@ -7609,87 +7894,87 @@ int readdata(char datafile[], int firsto
     /* Loops on waves */
     for (j=maxwav;j>=1;j--){
       for (iv=nqtv;iv>=1;iv--){  /* Loop  on time varying quantitative variables */
-				cutv(stra, strb, line, ' '); 
-				if(strb[0]=='.') { /* Missing value */
-					lval=-1;
-					cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
-					cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
-					if(isalpha(strb[1])) { /* .m or .d Really Missing value */
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
-						return 1;
-					}
-				}else{
-					errno=0;
-					/* what_kind_of_number(strb); */
-					dval=strtod(strb,&endptr); 
-					/* if( strb[0]=='\0' || (*endptr != '\0')){ */
-					/* if(strb != endptr && *endptr == '\0') */
-					/*    dval=dlval; */
-					/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
-					if( strb[0]=='\0' || (*endptr != '\0')){
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
-						return 1;
-					}
-					cotqvar[j][iv][i]=dval; 
-					cotvar[j][ntv+iv][i]=dval; 
-				}
-				strcpy(line,stra);
+	cutv(stra, strb, line, ' '); 
+	if(strb[0]=='.') { /* Missing value */
+	  lval=-1;
+	  cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
+	  cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
+	  if(isalpha(strb[1])) { /* .m or .d Really Missing value */
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
+	    return 1;
+	  }
+	}else{
+	  errno=0;
+	  /* what_kind_of_number(strb); */
+	  dval=strtod(strb,&endptr); 
+	  /* if( strb[0]=='\0' || (*endptr != '\0')){ */
+	  /* if(strb != endptr && *endptr == '\0') */
+	  /*    dval=dlval; */
+	  /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
+	  if( strb[0]=='\0' || (*endptr != '\0')){
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
+	    return 1;
+	  }
+	  cotqvar[j][iv][i]=dval; 
+	  cotvar[j][ntv+iv][i]=dval; 
+	}
+	strcpy(line,stra);
       }/* end loop ntqv */
       
       for (iv=ntv;iv>=1;iv--){  /* Loop  on time varying dummies */
-				cutv(stra, strb, line, ' '); 
-				if(strb[0]=='.') { /* Missing value */
-					lval=-1;
-				}else{
-					errno=0;
-					lval=strtol(strb,&endptr,10); 
-					/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-					if( strb[0]=='\0' || (*endptr != '\0')){
-						printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
-						fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
-						return 1;
-					}
-				}
-				if(lval <-1 || lval >1){
-					printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+	cutv(stra, strb, line, ' '); 
+	if(strb[0]=='.') { /* Missing value */
+	  lval=-1;
+	}else{
+	  errno=0;
+	  lval=strtol(strb,&endptr,10); 
+	  /*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+	  if( strb[0]=='\0' || (*endptr != '\0')){
+	    printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
+	    fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
+	    return 1;
+	  }
+	}
+	if(lval <-1 || lval >1){
+	  printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
- For example, for multinomial values like 1, 2 and 3,\n									\
- build V1=0 V2=0 for the reference value (1),\n													\
-        V1=1 V2=0 for (2) \n																						\
+ For example, for multinomial values like 1, 2 and 3,\n			\
+ build V1=0 V2=0 for the reference value (1),\n				\
+        V1=1 V2=0 for (2) \n						\
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
- output of IMaCh is often meaningless.\n																\
+ output of IMaCh is often meaningless.\n				\
  Exiting.\n",lval,linei, i,line,j);
-					fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+	  fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
- For example, for multinomial values like 1, 2 and 3,\n									\
- build V1=0 V2=0 for the reference value (1),\n													\
-        V1=1 V2=0 for (2) \n																						\
+ For example, for multinomial values like 1, 2 and 3,\n			\
+ build V1=0 V2=0 for the reference value (1),\n				\
+        V1=1 V2=0 for (2) \n						\
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
- output of IMaCh is often meaningless.\n																\
+ output of IMaCh is often meaningless.\n				\
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
-					return 1;
-				}
-				cotvar[j][iv][i]=(double)(lval);
-				strcpy(line,stra);
+	  return 1;
+	}
+	cotvar[j][iv][i]=(double)(lval);
+	strcpy(line,stra);
       }/* end loop ntv */
       
       /* Statuses  at wave */
       cutv(stra, strb, line, ' '); 
       if(strb[0]=='.') { /* Missing value */
-				lval=-1;
+	lval=-1;
       }else{
-				errno=0;
-				lval=strtol(strb,&endptr,10); 
-				/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-				if( strb[0]=='\0' || (*endptr != '\0')){
-					printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
-					fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
-					return 1;
-				}
+	errno=0;
+	lval=strtol(strb,&endptr,10); 
+	/*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+	if( strb[0]=='\0' || (*endptr != '\0')){
+	  printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+	  fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+	  return 1;
+	}
       }
       
       s[j][i]=lval;
@@ -7847,23 +8132,29 @@ int readdata(char datafile[], int firsto
   return (1);
 }
 
-void removespace(char **stri){/*, char stro[]) {*/
+void removefirstspace(char **stri){/*, char stro[]) {*/
   char *p1 = *stri, *p2 = *stri;
-  do
-    while (*p2 == ' ')
-      p2++;
-  while (*p1++ == *p2++);
-  *stri=p1; 
+  while (*p2 == ' ')
+    p2++; 
+  /* while ((*p1++ = *p2++) !=0) */
+  /*   ; */
+  /* do */
+  /*   while (*p2 == ' ') */
+  /*     p2++; */
+  /* while (*p1++ == *p2++); */
+  *stri=p2; 
 }
 
-int decoderesult ( char resultline[])
+int decoderesult ( char resultline[], int nres)
 /**< This routine decode one result line and returns the combination # of dummy covariates only **/
 {
-  int j=0, k=0;
+  int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0;
   char resultsav[MAXLINE];
+  int resultmodel[MAXLINE];
+  int modelresult[MAXLINE];
   char stra[80], strb[80], strc[80], strd[80],stre[80];
 
-  removespace(&resultline);
+  removefirstspace(&resultline);
   printf("decoderesult:%s\n",resultline);
 
   if (strstr(resultline,"v") !=0){
@@ -7875,13 +8166,19 @@ int decoderesult ( char resultline[])
   if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
   }
-
-  for(k=1; k<=j;k++){ /* Loop on total covariates of the model */
-    cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' 
-				     resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */
-    cutl(strc,strd,strb,'=');  /* strb:V4=1 strc=1 strd=V4 */
+  if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
+    printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
+    fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
+  }
+  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 ' ' 
+				      resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */
+       cutl(strc,strd,strb,'=');  /* strb:V4=1 strc=1 strd=V4 */
+    }else
+      cutl(strc,strd,resultsav,'=');
     Tvalsel[k]=atof(strc); /* 1 */
-
+    
     cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */;
     Tvarsel[k]=atoi(strc);
     /* Typevarsel[k]=1;  /\* 1 for age product *\/ */
@@ -7889,14 +8186,83 @@ int decoderesult ( char resultline[])
     if (nbocc(stra,'=') >0)
       strcpy(resultsav,stra); /* and analyzes it */
   }
+  /* Checking for missing or useless values in comparison of current model needs */
+  for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+    if(Typevar[k1]==0){ /* Single covariate in model */
+      match=0;
+      for(k2=1; k2 <=j;k2++){/* result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+	if(Tvar[k1]==Tvarsel[k2]) {/* Tvar[1]=5 == Tvarsel[2]=5   */
+	  modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
+	  match=1;
+	  break;
+	}
+      }
+      if(match == 0){
+	printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
+      }
+    }
+  }
+  /* Checking for missing or useless values in comparison of current model needs */
+  for(k2=1; k2 <=j;k2++){ /* result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+    match=0;
+    for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+      if(Typevar[k1]==0){ /* Single */
+	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
+	  resultmodel[k1]=k2;  /* resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
+	  ++match;
+	}
+      }
+    }
+    if(match == 0){
+      printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model);
+    }else if(match > 1){
+      printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
+    }
+  }
+      
+  /* We need to deduce which combination number is chosen and save quantitative values */
+  /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+  /* result line V4=1 V5=25.1 V3=0  V2=8 V1=1 */
+  /* should give a combination of dummy V4=1, V3=0, V1=1 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 5 + (1offset) = 6*/
+  /* result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
+  /* should give a combination of dummy V4=1, V3=1, V1=0 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 3 + (1offset) = 4*/
+  /*    1 0 0 0 */
+  /*    2 1 0 0 */
+  /*    3 0 1 0 */ 
+  /*    4 1 1 0 */ /* V4=1, V3=1, V1=0 */
+  /*    5 0 0 1 */
+  /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */
+  /*    7 0 1 1 */
+  /*    8 1 1 1 */
+  /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */
+  /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */
+  /* V5*age V5 known which value for nres?  */
+  /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */
+  for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */
+    if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */
+      k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */
+      k2=(int)Tvarsel[k3]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
+      k+=Tvalsel[k3]*pow(2,k4);  /*  Tvalsel[1]=1  */
+      Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1)  Tresult[nres][2]=0(V3=0) */
+      Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+      Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
+      printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4);
+      k4++;;
+    }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */
+      k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
+      k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
+      Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
+      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
+      Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
+      printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
+      k4q++;;
+    }
+  }
+  
+  TKresult[nres]=++k; /* Combination for the nresult and the model */
   return (0);
 }
-int selected( int kvar){ /* Selected combination of covariates */
-  if(Tvarsel[kvar])
-    return (0);
-  else
-    return(1);
-}
+
 int decodemodel( char model[], int lastobs)
  /**< This routine decodes the model and returns:
 	* Model  V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age
@@ -7913,7 +8279,7 @@ int decodemodel( char model[], int lasto
 	* - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
 	*/
 {
-  int i, j, k, ks;
+  int i, j, k, ks, v;
   int  j1, k1, k2, k3, k4;
   char modelsav[80];
   char stra[80], strb[80], strc[80], strd[80],stre[80];
@@ -7936,21 +8302,21 @@ int decodemodel( char model[], int lasto
     if ((strpt=strstr(model,"age*age")) !=0){
       printf(" strpt=%s, model=%s\n",strpt, model);
       if(strpt != model){
-				printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+	printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model);
-				fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+	fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model); fflush(ficlog);
-				return 1;
+	return 1;
       }
       nagesqr=1;
       if (strstr(model,"+age*age") !=0)
-				substrchaine(modelsav, model, "+age*age");
+	substrchaine(modelsav, model, "+age*age");
       else if (strstr(model,"age*age+") !=0)
-				substrchaine(modelsav, model, "age*age+");
+	substrchaine(modelsav, model, "age*age+");
       else 
-				substrchaine(modelsav, model, "age*age");
+	substrchaine(modelsav, model, "age*age");
     }else
       nagesqr=0;
     if (strlen(modelsav) >1){
@@ -8020,67 +8386,67 @@ int decodemodel( char model[], int lasto
       }
       cptcovage=0;
       for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
-				cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
+	cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
 					 modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ 
-				if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
-				/*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
-				/*scanf("%d",i);*/
-				if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
-					cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
-					if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
-						/* covar is not filled and then is empty */
-						cptcovprod--;
-						cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
-						Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
-						Typevar[k]=1;  /* 1 for age product */
-						cptcovage++; /* Sums the number of covariates which include age as a product */
-						Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
-						/*printf("stre=%s ", stre);*/
-					} else if (strcmp(strd,"age")==0) { /* or age*Vn */
-						cptcovprod--;
-						cutl(stre,strb,strc,'V');
-						Tvar[k]=atoi(stre);
-						Typevar[k]=1;  /* 1 for age product */
-						cptcovage++;
-						Tage[cptcovage]=k;
-					} else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
-						/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
-						cptcovn++;
-						cptcovprodnoage++;k1++;
-						cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-						Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
-																								because this model-covariate is a construction we invent a new column
-																								which is after existing variables ncovcol+nqv+ntv+nqtv + k1
-																								If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
-																								Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
-						Typevar[k]=2;  /* 2 for double fixed dummy covariates */
-						cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
-						Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
-						Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */
-						Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
-						Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
-						k2=k2+2;  /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
-						/* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */
-						/* Tvar[cptcovt+k2+1]=Tvard[k1][2];  /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
+	if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+	/*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+	/*scanf("%d",i);*/
+	if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
+	  cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+	  if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+	    /* covar is not filled and then is empty */
+	    cptcovprod--;
+	    cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+	    Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+	    Typevar[k]=1;  /* 1 for age product */
+	    cptcovage++; /* Sums the number of covariates which include age as a product */
+	    Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+	    /*printf("stre=%s ", stre);*/
+	  } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+	    cptcovprod--;
+	    cutl(stre,strb,strc,'V');
+	    Tvar[k]=atoi(stre);
+	    Typevar[k]=1;  /* 1 for age product */
+	    cptcovage++;
+	    Tage[cptcovage]=k;
+	  } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
+	    /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+	    cptcovn++;
+	    cptcovprodnoage++;k1++;
+	    cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+	    Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
+						because this model-covariate is a construction we invent a new column
+						which is after existing variables ncovcol+nqv+ntv+nqtv + k1
+						If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
+						Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
+	    Typevar[k]=2;  /* 2 for double fixed dummy covariates */
+	    cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
+	    Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
+	    Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */
+	    Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+	    Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+	    k2=k2+2;  /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
+	    /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */
+	    /* Tvar[cptcovt+k2+1]=Tvard[k1][2];  /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
             /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
-						/*                     1  2   3      4     5 | Tvar[5+1)=1, Tvar[7]=2   */
-						for (i=1; i<=lastobs;i++){
-							/* Computes the new covariate which is a product of
-								 covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
-							covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
-						}
-					} /* End age is not in the model */
-				} /* End if model includes a product */
-				else { /* no more sum */
-					/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
-					/*  scanf("%d",i);*/
-					cutl(strd,strc,strb,'V');
-					ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */
-					cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */
-					Tvar[k]=atoi(strd);
-					Typevar[k]=0;  /* 0 for simple covariates */
-				}
-				strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
+	    /*                     1  2   3      4     5 | Tvar[5+1)=1, Tvar[7]=2   */
+	    for (i=1; i<=lastobs;i++){
+	      /* Computes the new covariate which is a product of
+		 covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+	      covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+	    }
+	  } /* End age is not in the model */
+	} /* End if model includes a product */
+	else { /* no more sum */
+	  /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+	  /*  scanf("%d",i);*/
+	  cutl(strd,strc,strb,'V');
+	  ks++; /**< Number of simple covariates dummy or quantitative, fixe or varying */
+	  cptcovn++; /** V4+V3+V5: V4 and V3 timevarying dummy covariates, V5 timevarying quantitative */
+	  Tvar[k]=atoi(strd);
+	  Typevar[k]=0;  /* 0 for simple covariates */
+	}
+	strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
 				/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
 				  scanf("%d",i);*/
       } /* end of loop + on total covariates */
@@ -8119,223 +8485,245 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
 Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
-
-  for(k=1, ncovf=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
-    if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy (<=ncovcol) covariates */
+  for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
+  for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
+    if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;
       Dummy[k]= 0;
       ncoveff++;
       ncovf++;
-			modell[k].maintype= FTYPE;
-			TvarF[ncovf]=Tvar[k];
-			TvarFind[ncovf]=k;
+      nsd++;
+      modell[k].maintype= FTYPE;
+      TvarsD[nsd]=Tvar[k];
+      TvarsDind[nsd]=k;
+      TvarF[ncovf]=Tvar[k];
+      TvarFind[ncovf]=k;
       TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-    }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */
+    }else if( Tvar[k] <=ncovcol &&  Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */
+      Fixed[k]= 0;
+      Dummy[k]= 0;
+      ncoveff++;
+      ncovf++;
+      modell[k].maintype= FTYPE;
+      TvarF[ncovf]=Tvar[k];
+      TvarFind[ncovf]=k;
+      TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+      TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
+    }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */
       Fixed[k]= 0;
       Dummy[k]= 1;
       nqfveff++;
-			modell[k].maintype= FTYPE;
-			modell[k].subtype= FQ;
+      modell[k].maintype= FTYPE;
+      modell[k].subtype= FQ;
+      nsq++;
+      TvarsQ[nsq]=Tvar[k];
+      TvarsQind[nsq]=k;
       ncovf++;
-			TvarF[ncovf]=Tvar[k];
-			TvarFind[ncovf]=k;
+      TvarF[ncovf]=Tvar[k];
+      TvarFind[ncovf]=k;
       TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
       TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
-    }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){
+    }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */
       Fixed[k]= 1;
       Dummy[k]= 0;
       ntveff++; /* Only simple time varying dummy variable */
-			modell[k].maintype= VTYPE;
-			modell[k].subtype= VD;
-			ncovv++; /* Only simple time varying variables */
-			TvarV[ncovv]=Tvar[k];
-			TvarVind[ncovv]=k;
+      modell[k].maintype= VTYPE;
+      modell[k].subtype= VD;
+      nsd++;
+      TvarsD[nsd]=Tvar[k];
+      TvarsDind[nsd]=k;
+      ncovv++; /* Only simple time varying variables */
+      TvarV[ncovv]=Tvar[k];
+      TvarVind[ncovv]=k;
       TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4  TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
       TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
       printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
       printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv);
     }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv  && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/
-			Fixed[k]= 1;
-			Dummy[k]= 1;
-			nqtveff++;
-			modell[k].maintype= VTYPE;
-			modell[k].subtype= VQ;
-			ncovv++; /* Only simple time varying variables */
-			TvarV[ncovv]=Tvar[k];
-			TvarVind[ncovv]=k;
+      Fixed[k]= 1;
+      Dummy[k]= 1;
+      nqtveff++;
+      modell[k].maintype= VTYPE;
+      modell[k].subtype= VQ;
+      ncovv++; /* Only simple time varying variables */
+      nsq++;
+      TvarsQ[nsq]=Tvar[k];
+      TvarsQind[nsq]=k;
+      TvarV[ncovv]=Tvar[k];
+      TvarVind[ncovv]=k;
       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 */
       TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
-			TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
-			/* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
-			printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);
+      TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
+      /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
+      printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);
       printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv);
     }else if (Typevar[k] == 1) {  /* product with age */
-			ncova++;
-			TvarA[ncova]=Tvar[k];
-			TvarAind[ncova]=k;
+      ncova++;
+      TvarA[ncova]=Tvar[k];
+      TvarAind[ncova]=k;
       if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
-				Fixed[k]= 2;
-				Dummy[k]= 2;
-				modell[k].maintype= ATYPE;
-				modell[k].subtype= APFD;
-				/* ncoveff++; */
+      	Fixed[k]= 2;
+      	Dummy[k]= 2;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APFD;
+      	/* ncoveff++; */
       }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
-				Fixed[k]= 2;
-				Dummy[k]= 3;
-				modell[k].maintype= ATYPE;
-				modell[k].subtype= APFQ;		/*	Product age * fixed quantitative */
-				/* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
+      	Fixed[k]= 2;
+      	Dummy[k]= 3;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APFQ;		/*	Product age * fixed quantitative */
+      	/* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv ){
-				Fixed[k]= 3;
-				Dummy[k]= 2;
-				modell[k].maintype= ATYPE;
-				modell[k].subtype= APVD;		/*	Product age * varying dummy */
-				/* ntveff++; /\* Only simple time varying dummy variable *\/ */
+      	Fixed[k]= 3;
+      	Dummy[k]= 2;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APVD;		/*	Product age * varying dummy */
+      	/* ntveff++; /\* Only simple time varying dummy variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
-				Fixed[k]= 3;
-				Dummy[k]= 3;
-				modell[k].maintype= ATYPE;
-				modell[k].subtype= APVQ;		/*	Product age * varying quantitative */
-				/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
+      	Fixed[k]= 3;
+      	Dummy[k]= 3;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APVQ;		/*	Product age * varying quantitative */
+      	/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
       }
     }else if (Typevar[k] == 2) {  /* product without age */
       k1=Tposprod[k];
       if(Tvard[k1][1] <=ncovcol){
-				if(Tvard[k1][2] <=ncovcol){
-					Fixed[k]= 1;
-					Dummy[k]= 0;
-					modell[k].maintype= FTYPE;
-					modell[k].subtype= FPDD;		/*	Product fixed dummy * fixed dummy */
-					ncovf++; /* Fixed variables without age */
-					TvarF[ncovf]=Tvar[k];
-					TvarFind[ncovf]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv){
-					Fixed[k]= 0;  /* or 2 ?*/
-					Dummy[k]= 1;
-					modell[k].maintype= FTYPE;
-					modell[k].subtype= FPDQ;		/*	Product fixed dummy * fixed quantitative */
-					ncovf++; /* Varying variables without age */
-					TvarF[ncovf]=Tvar[k];
-					TvarFind[ncovf]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-					Fixed[k]= 1;
-					Dummy[k]= 0;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDD;		/*	Product fixed dummy * varying dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product fixed dummy * varying quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDD;		/*	Product fixed dummy * fixed dummy */
+    	  ncovf++; /* Fixed variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 0;  /* or 2 ?*/
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDQ;		/*	Product fixed dummy * fixed quantitative */
+    	  ncovf++; /* Varying variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product fixed dummy * varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product fixed dummy * varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv){
-				if(Tvard[k1][2] <=ncovcol){
-					Fixed[k]= 0;  /* or 2 ?*/
-					Dummy[k]= 1;
-					modell[k].maintype= FTYPE;
-					modell[k].subtype= FPDQ;		/*	Product fixed quantitative * fixed dummy */
-					ncovf++; /* Fixed variables without age */
-					TvarF[ncovf]=Tvar[k];
-					TvarFind[ncovf]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product fixed quantitative * varying dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPQQ;		/*	Product fixed quantitative * varying quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 0;  /* or 2 ?*/
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDQ;		/*	Product fixed quantitative * fixed dummy */
+    	  ncovf++; /* Fixed variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product fixed quantitative * varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product fixed quantitative * varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
-				if(Tvard[k1][2] <=ncovcol){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDD;		/*	Product time varying dummy * fixed dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product time varying dummy * fixed quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-					Fixed[k]= 1;
-					Dummy[k]= 0;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDD;		/*	Product time varying dummy * time varying dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product time varying dummy * time varying quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product time varying dummy * fixed dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * fixed quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product time varying dummy * time varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * time varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
-				if(Tvard[k1][2] <=ncovcol){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product time varying quantitative * fixed dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPQQ;		/*	Product time varying quantitative * fixed quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPDQ;		/*	Product time varying quantitative * time varying dummy */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-					Fixed[k]= 1;
-					Dummy[k]= 1;
-					modell[k].maintype= VTYPE;
-					modell[k].subtype= VPQQ;		/*	Product time varying quantitative * time varying quantitative */
-					ncovv++; /* Varying variables without age */
-					TvarV[ncovv]=Tvar[k];
-					TvarVind[ncovv]=k;
-				} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * fixed dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * fixed quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * time varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * time varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else{
-				printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
-				fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
-      } /* end k1 */
+    	printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+    	fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+      } /*end k1*/
     }else{
       printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
       fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
@@ -8348,28 +8736,28 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
   for(k1=1; k1<= cptcovt;k1++){
     for(k2=1; k2 <k1;k2++){
       if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){
-				if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
-					if(Tvar[k1]==Tvar[k2]){
-						printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
-						fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
-						return(1);
-					}
-				}else if (Typevar[k1] ==2){
-					k3=Tposprod[k1];
-					k4=Tposprod[k2];
-					if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
-						printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
-						fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
-						return(1);
-					}
-				}
+	if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
+	  if(Tvar[k1]==Tvar[k2]){
+	    printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+	    return(1);
+	  }
+	}else if (Typevar[k1] ==2){
+	  k3=Tposprod[k1];
+	  k4=Tposprod[k2];
+	  if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
+	    printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+	    return(1);
+	  }
+	}
       }
     }
   }
   printf("ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
   fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
-  printf("ncovf=%d, ncovv=%d, ncova=%d\n",ncovf,ncovv,ncova);
-  fprintf(ficlog,"ncovf=%d, ncovv=%d, ncova=%d\n",ncovf,ncovv,ncova);
+  printf("ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd,nsq);
+  fprintf(ficlog,"ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd, nsq);
   return (0); /* with covar[new additional covariate if product] and Tage if age */ 
   /*endread:*/
   printf("Exiting decodemodel: ");
@@ -8684,7 +9072,7 @@ void syscompilerinfo(int logged)
 
 int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
   /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
-  int i, j, k, i1 ;
+  int i, j, k, i1, k4=0, nres=0 ;
   /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;
   double tot;
@@ -8709,56 +9097,66 @@ int prevalence_limit(double *p, double *
   agelim=agemaxpar;
 
   /* i1=pow(2,ncoveff); */
-  i1=pow(2,cptcoveff); /* Number of dummy covariates */
+  i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
   if (cptcovn < 1){i1=1;}
 
-  for(k=1; k<=i1;k++){
-  /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
-    /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
-    //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-    /* k=k+1; */
-    /* to clean */
-    //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-    fprintf(ficrespl,"#******");
-    printf("#******");
-    fprintf(ficlog,"#******");
-    for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-      fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
-      printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    }
-    fprintf(ficrespl,"******\n");
-    printf("******\n");
-    fprintf(ficlog,"******\n");
-    if(invalidvarcomb[k]){
-      printf("\nCombination (%d) ignored because no case \n",k); 
-      fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k); 
-      fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k); 
-						continue;
-    }
+  for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+      if(TKresult[nres]!= k)
+	continue;
 
-    fprintf(ficrespl,"#Age ");
-    for(j=1;j<=cptcoveff;j++) {
-      fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    }
-    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
-    fprintf(ficrespl,"Total Years_to_converge\n");
+      /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+      /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
+      //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+      /* k=k+1; */
+      /* to clean */
+      //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+      fprintf(ficrespl,"#******");
+      printf("#******");
+      fprintf(ficlog,"#******");
+      for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
+	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficrespl," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
+      fprintf(ficrespl,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
+      if(invalidvarcomb[k]){
+	printf("\nCombination (%d) ignored because no case \n",k); 
+	fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k); 
+	fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k); 
+	continue;
+      }
+
+      fprintf(ficrespl,"#Age ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
+      fprintf(ficrespl,"Total Years_to_converge\n");
     
-    for (age=agebase; age<=agelim; age++){
-      /* for (age=agebase; age<=agebase; age++){ */
-      prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
-      fprintf(ficrespl,"%.0f ",age );
-      for(j=1;j<=cptcoveff;j++)
-	fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      tot=0.;
-      for(i=1; i<=nlstate;i++){
-	tot +=  prlim[i][i];
-	fprintf(ficrespl," %.5f", prlim[i][i]);
-      }
-      fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
-    } /* Age */
-    /* was end of cptcod */
-  } /* cptcov */
+      for (age=agebase; age<=agelim; age++){
+	/* for (age=agebase; age<=agebase; age++){ */
+	prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres);
+	fprintf(ficrespl,"%.0f ",age );
+	for(j=1;j<=cptcoveff;j++)
+	  fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	tot=0.;
+	for(i=1; i<=nlstate;i++){
+	  tot +=  prlim[i][i];
+	  fprintf(ficrespl," %.5f", prlim[i][i]);
+	}
+	fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
+      } /* Age */
+      /* was end of cptcod */
+    } /* cptcov */
+  } /* nres */
   return 0;
 }
 
@@ -8768,7 +9166,7 @@ int back_prevalence_limit(double *p, dou
 	/* Computes the back prevalence limit  for any combination	of covariate values 
    * at any age between ageminpar and agemaxpar
 	 */
-  int i, j, k, i1 ;
+  int i, j, k, i1, nres=0 ;
   /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;
   double tot;
@@ -8799,61 +9197,69 @@ int back_prevalence_limit(double *p, dou
   i1=pow(2,cptcoveff);
   if (cptcovn < 1){i1=1;}
   
-  for(k=1; k<=i1;k++){ 
-    //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-    fprintf(ficresplb,"#******");
-    printf("#******");
-    fprintf(ficlog,"#******");
-    for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-      fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    }
-    fprintf(ficresplb,"******\n");
-    printf("******\n");
-    fprintf(ficlog,"******\n");
-    if(invalidvarcomb[k]){
-      printf("\nCombination (%d) ignored because no cases \n",k); 
-      fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k); 
-      fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k); 
-      continue;
-    }
+  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
+      //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+      fprintf(ficresplb,"#******");
+      printf("#******");
+      fprintf(ficlog,"#******");
+      for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+	fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }
+      fprintf(ficresplb,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
+      if(invalidvarcomb[k]){
+	printf("\nCombination (%d) ignored because no cases \n",k); 
+	fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k); 
+	fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k); 
+	continue;
+      }
     
-    fprintf(ficresplb,"#Age ");
-    for(j=1;j<=cptcoveff;j++) {
-      fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    }
-    for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
-    fprintf(ficresplb,"Total Years_to_converge\n");
+      fprintf(ficresplb,"#Age ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
+      fprintf(ficresplb,"Total Years_to_converge\n");
     
     
-    for (age=agebase; age<=agelim; age++){
-      /* for (age=agebase; age<=agebase; age++){ */
-      if(mobilavproj > 0){
-	/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
-	/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-	bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
-      }else if (mobilavproj == 0){
-	printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-	fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-	exit(1);
-      }else{
-	/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-	bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
-      }
-      fprintf(ficresplb,"%.0f ",age );
-      for(j=1;j<=cptcoveff;j++)
-	fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-      tot=0.;
-      for(i=1; i<=nlstate;i++){
-	tot +=  bprlim[i][i];
-	fprintf(ficresplb," %.5f", bprlim[i][i]);
-      }
-      fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
-    } /* Age */
-    /* was end of cptcod */
-  } /* cptcov */
-  
+      for (age=agebase; age<=agelim; age++){
+	/* for (age=agebase; age<=agebase; age++){ */
+	if(mobilavproj > 0){
+	  /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
+	  /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+	  bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+	}else if (mobilavproj == 0){
+	  printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+	  fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+	  exit(1);
+	}else{
+	  /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+	  bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+	}
+	fprintf(ficresplb,"%.0f ",age );
+	for(j=1;j<=cptcoveff;j++)
+	  fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	tot=0.;
+	for(i=1; i<=nlstate;i++){
+	  tot +=  bprlim[i][i];
+	  fprintf(ficresplb," %.5f", bprlim[i][i]);
+	}
+	fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+      } /* Age */
+      /* was end of cptcod */
+    } /* end of any combination */
+  } /* end of nres */  
   /* hBijx(p, bage, fage); */
   /* fclose(ficrespijb); */
   
@@ -8867,7 +9273,7 @@ int hPijx(double *p, int bage, int fage)
   int agelim;
   int hstepm;
   int nhstepm;
-  int h, i, i1, j, k;
+  int h, i, i1, j, k, k4, nres=0;
 
   double agedeb;
   double ***p3mat;
@@ -8894,10 +9300,17 @@ int hPijx(double *p, int bage, int fage)
 		/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
 		/*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
 		/*  	k=k+1;  */
-    for (k=1; k <= (int) pow(2,cptcoveff); k++){
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
+    for(k=1; k<=i1;k++){
+      if(TKresult[nres]!= k)
+	continue;
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
 	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	fprintf(ficrespij," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      }
       fprintf(ficrespij,"******\n");
       
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -8908,7 +9321,7 @@ int hPijx(double *p, int bage, int fage)
 	
 	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
 	oldm=oldms;savm=savms;
-	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
+	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k, nres);  
 	fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
 	for(i=1; i<=nlstate;i++)
 	  for(j=1; j<=nlstate+ndeath;j++)
@@ -8938,7 +9351,7 @@ int hPijx(double *p, int bage, int fage)
 	int ageminl;
   int hstepm;
   int nhstepm;
-  int h, i, i1, j, k;
+  int h, i, i1, j, k, nres;
 	
   double agedeb;
   double ***p3mat;
@@ -8966,48 +9379,54 @@ int hPijx(double *p, int bage, int fage)
   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
   /*  	k=k+1;  */
-  for (k=1; k <= (int) pow(2,cptcoveff); k++){
-    fprintf(ficrespijb,"\n#****** ");
-    for(j=1;j<=cptcoveff;j++)
-      fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-    fprintf(ficrespijb,"******\n");
-    if(invalidvarcomb[k]){
-      fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); 
-      continue;
-    }
-    
-    /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
-    for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
-      /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
-      nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-      nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
-      
-      /*	  nhstepm=nhstepm*YEARM; aff par mois*/
-      
-      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-      /* oldm=oldms;savm=savms; */
-      /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
-      hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
-      /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
-      fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
-      for(i=1; i<=nlstate;i++)
-	for(j=1; j<=nlstate+ndeath;j++)
-	  fprintf(ficrespijb," %1d-%1d",i,j);
-      fprintf(ficrespijb,"\n");
-      for (h=0; h<=nhstepm; h++){
-	/*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
-	fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
-	/* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+  for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
+      fprintf(ficrespijb,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++)
+	fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }
+      fprintf(ficrespijb,"******\n");
+      if(invalidvarcomb[k]){
+	fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); 
+	continue;
+      }
+      
+      /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+      for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
+	/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
+	nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+	nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
+	
+	/*	  nhstepm=nhstepm*YEARM; aff par mois*/
+	
+	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	/* oldm=oldms;savm=savms; */
+	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
+	hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
+	fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
 	for(i=1; i<=nlstate;i++)
 	  for(j=1; j<=nlstate+ndeath;j++)
-	    fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+	    fprintf(ficrespijb," %1d-%1d",i,j);
 	fprintf(ficrespijb,"\n");
-      }
-      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-      fprintf(ficrespijb,"\n");
-    }
-    /*}*/
-  }
+	for (h=0; h<=nhstepm; h++){
+	  /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+	  fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
+	  /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+	  for(i=1; i<=nlstate;i++)
+	    for(j=1; j<=nlstate+ndeath;j++)
+	      fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+	  fprintf(ficrespijb,"\n");
+	}
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	fprintf(ficrespijb,"\n");
+      } /* end age deb */
+    } /* end combination */
+  } /* end nres */
   return 0;
  } /*  hBijx */
 
@@ -9034,6 +9453,7 @@ int main(int argc, char *argv[])
   int itimes;
   int NDIM=2;
   int vpopbased=0;
+  int nres=0;
 
   char ca[32], cb[32];
   /*  FILE *fichtm; *//* Html File */
@@ -9052,7 +9472,7 @@ int main(int argc, char *argv[])
   char line[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
 
-  char model[MAXLINE], modeltemp[MAXLINE];
+  char  modeltemp[MAXLINE];
   char resultline[MAXLINE];
   
   char pathr[MAXLINE], pathimach[MAXLINE]; 
@@ -9419,41 +9839,41 @@ int main(int argc, char *argv[])
     
     param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     for(i=1; i <=nlstate; i++){
-			j=0;
+      j=0;
       for(jj=1; jj <=nlstate+ndeath; jj++){
-				if(jj==i) continue;
-				j++;
-				fscanf(ficpar,"%1d%1d",&i1,&j1);
-				if ((i1 != i) || (j1 != jj)){
-					printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
+	if(jj==i) continue;
+	j++;
+	fscanf(ficpar,"%1d%1d",&i1,&j1);
+	if ((i1 != i) || (j1 != jj)){
+	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
 It might be a problem of design; if ncovcol and the model are correct\n \
 run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
-					exit(1);
-				}
-				fprintf(ficparo,"%1d%1d",i1,j1);
-				if(mle==1)
-					printf("%1d%1d",i,jj);
-				fprintf(ficlog,"%1d%1d",i,jj);
-				for(k=1; k<=ncovmodel;k++){
-					fscanf(ficpar," %lf",&param[i][j][k]);
-					if(mle==1){
-						printf(" %lf",param[i][j][k]);
-						fprintf(ficlog," %lf",param[i][j][k]);
-					}
-					else
-						fprintf(ficlog," %lf",param[i][j][k]);
-					fprintf(ficparo," %lf",param[i][j][k]);
-				}
-				fscanf(ficpar,"\n");
-				numlinepar++;
-				if(mle==1)
-					printf("\n");
-				fprintf(ficlog,"\n");
-				fprintf(ficparo,"\n");
+	  exit(1);
+	}
+	fprintf(ficparo,"%1d%1d",i1,j1);
+	if(mle==1)
+	  printf("%1d%1d",i,jj);
+	fprintf(ficlog,"%1d%1d",i,jj);
+	for(k=1; k<=ncovmodel;k++){
+	  fscanf(ficpar," %lf",&param[i][j][k]);
+	  if(mle==1){
+	    printf(" %lf",param[i][j][k]);
+	    fprintf(ficlog," %lf",param[i][j][k]);
+	  }
+	  else
+	    fprintf(ficlog," %lf",param[i][j][k]);
+	  fprintf(ficparo," %lf",param[i][j][k]);
+	}
+	fscanf(ficpar,"\n");
+	numlinepar++;
+	if(mle==1)
+	  printf("\n");
+	fprintf(ficlog,"\n");
+	fprintf(ficparo,"\n");
       }
     }  
     fflush(ficlog);
-
+    
     /* Reads scales values */
     p=param[1][1];
     
@@ -9470,29 +9890,29 @@ run imach with mle=-1 to get a correct t
 
     for(i=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath-1; j++){
-				fscanf(ficpar,"%1d%1d",&i1,&j1);
-				if ( (i1-i) * (j1-j) != 0){
-					printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
-					exit(1);
-				}
-				printf("%1d%1d",i,j);
-				fprintf(ficparo,"%1d%1d",i1,j1);
-				fprintf(ficlog,"%1d%1d",i1,j1);
-				for(k=1; k<=ncovmodel;k++){
-					fscanf(ficpar,"%le",&delti3[i][j][k]);
-					printf(" %le",delti3[i][j][k]);
-					fprintf(ficparo," %le",delti3[i][j][k]);
-					fprintf(ficlog," %le",delti3[i][j][k]);
-				}
-				fscanf(ficpar,"\n");
-				numlinepar++;
-				printf("\n");
-				fprintf(ficparo,"\n");
-				fprintf(ficlog,"\n");
+	fscanf(ficpar,"%1d%1d",&i1,&j1);
+	if ( (i1-i) * (j1-j) != 0){
+	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
+	  exit(1);
+	}
+	printf("%1d%1d",i,j);
+	fprintf(ficparo,"%1d%1d",i1,j1);
+	fprintf(ficlog,"%1d%1d",i1,j1);
+	for(k=1; k<=ncovmodel;k++){
+	  fscanf(ficpar,"%le",&delti3[i][j][k]);
+	  printf(" %le",delti3[i][j][k]);
+	  fprintf(ficparo," %le",delti3[i][j][k]);
+	  fprintf(ficlog," %le",delti3[i][j][k]);
+	}
+	fscanf(ficpar,"\n");
+	numlinepar++;
+	printf("\n");
+	fprintf(ficparo,"\n");
+	fprintf(ficlog,"\n");
       }
     }
     fflush(ficlog);
-		
+    
     /* Reads covariance matrix */
     delti=delti3[1][1];
 		
@@ -9582,15 +10002,15 @@ Please run with mle=-1 to get a correct
   agedc=vector(1,n);
   cod=ivector(1,n);
   for(i=1;i<=n;i++){
-		num[i]=0;
-		moisnais[i]=0;
-		annais[i]=0;
-		moisdc[i]=0;
-		andc[i]=0;
-		agedc[i]=0;
-		cod[i]=0;
-		weight[i]=1.0; /* Equal weights, 1 by default */
-	}
+    num[i]=0;
+    moisnais[i]=0;
+    annais[i]=0;
+    moisdc[i]=0;
+    andc[i]=0;
+    agedc[i]=0;
+    cod[i]=0;
+    weight[i]=1.0; /* Equal weights, 1 by default */
+  }
   mint=matrix(1,maxwav,1,n);
   anint=matrix(1,maxwav,1,n);
   s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ 
@@ -9603,14 +10023,18 @@ Please run with mle=-1 to get a correct
     goto end;
 
   /* Calculation of the number of parameters from char model */
-    /*    modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 
+  /*    modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 
 	k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tag[cptcovage=1]=4
 	k=3 V4 Tvar[k=3]= 4 (from V4)
 	k=2 V1 Tvar[k=2]= 1 (from V1)
 	k=1 Tvar[1]=2 (from V2)
-    */
-
-	Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
+  */
+  
+  Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
+  TvarsDind=ivector(1,NCOVMAX); /*  */
+  TvarsD=ivector(1,NCOVMAX); /*  */
+  TvarsQind=ivector(1,NCOVMAX); /*  */
+  TvarsQ=ivector(1,NCOVMAX); /*  */
   TvarF=ivector(1,NCOVMAX); /*  */
   TvarFind=ivector(1,NCOVMAX); /*  */
   TvarV=ivector(1,NCOVMAX); /*  */
@@ -10407,6 +10831,7 @@ Please run with mle=-1 to get a correct
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);
+      fputs(line,ficres);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
@@ -10423,6 +10848,7 @@ Please run with mle=-1 to get a correct
       fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);
       fputs(line,ficparo);
+      fputs(line,ficres);
     }
     ungetc(c,ficpar);
     
@@ -10433,6 +10859,7 @@ Please run with mle=-1 to get a correct
     /* day and month of proj2 are not used but only year anproj2.*/
     
     /* Results */
+    nresult=0;
     while(fgets(line, MAXLINE, ficpar)) {
       /* If line starts with a # it is a comment */
       if (line[0] == '#') {
@@ -10440,24 +10867,37 @@ Please run with mle=-1 to get a correct
 	fputs(line,stdout);
 	fputs(line,ficparo);
 	fputs(line,ficlog);
+	fputs(line,ficres);
 	continue;
       }else
 	break;
     }
+    if (!feof(ficpar))
     while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
-      if (num_filled == 0)
+      if (num_filled == 0){
 	resultline[0]='\0';
-      else if (num_filled != 1){
+      break;
+      } else if (num_filled != 1){
 	printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
       }
-      printf("Result %d: result line should be at minimum 'line=' %s, result=%s\n",num_filled, line, resultline);
-      decoderesult(resultline);
+      nresult++; /* Sum of resultlines */
+      printf("Result %d: result=%s\n",nresult, resultline);
+      if(nresult > MAXRESULTLINES){
+	printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	goto end;
+      }
+      decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+      fprintf(ficparo,"result: %s\n",resultline);
+      fprintf(ficres,"result: %s\n",resultline);
+      fprintf(ficlog,"result: %s\n",resultline);
       while(fgets(line, MAXLINE, ficpar)) {
 	/* If line starts with a # it is a comment */
 	if (line[0] == '#') {
 	  numlinepar++;
 	  fputs(line,stdout);
 	  fputs(line,ficparo);
+	  fputs(line,ficres);
 	  fputs(line,ficlog);
 	  continue;
 	}else
@@ -10467,7 +10907,7 @@ Please run with mle=-1 to get a correct
 	break;
       else{ /* Processess output results for this combination of covariate values */
       }				   
-    }
+    } /* end while */
 
 
     
@@ -10538,6 +10978,7 @@ Please run with mle=-1 to get a correct
 	    mobaverages[i][j][k]=0.;
       mobaverage=mobaverages;
       if (mobilav!=0) {
+	printf("Movingaveraging observed prevalence\n");
 	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
 	  fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
 	  printf(" Error in movingaverage mobilav=%d\n",mobilav);
@@ -10546,6 +10987,7 @@ Please run with mle=-1 to get a correct
       /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */
       /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
       else if (mobilavproj !=0) {
+	printf("Movingaveraging projected observed prevalence\n");
 	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
 	  fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
 	  printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
@@ -10600,17 +11042,32 @@ Please run with mle=-1 to get a correct
     }
     printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
+
+    pstamp(ficreseij);
  		
-    for (k=1; k <= (int) pow(2,cptcoveff); k++){ /* For any combination of dummy covariates, fixed and varying */
+    i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
+    if (cptcovn < 1){i1=1;}
+    
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
       fprintf(ficreseij,"\n#****** ");
+      printf("\n#****** ");
       for(j=1;j<=cptcoveff;j++) {
 	fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficreseij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
       }
       fprintf(ficreseij,"******\n");
+      printf("******\n");
       
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;
-      evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
+      evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart, nres);  
       
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
     }
@@ -10661,15 +11118,26 @@ Please run with mle=-1 to get a correct
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
           
-    for (k=1; k <= (int) pow(2,cptcoveff); k++){
-      printf("\n#****** ");
-      fprintf(ficrest,"\n#****** ");
-      fprintf(ficlog,"\n#****** ");
+    i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
+    if (cptcovn < 1){i1=1;}
+    
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
+    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+      if(TKresult[nres]!= k)
+	continue;
+      printf("\n#****** Selected:");
+      fprintf(ficrest,"\n#****** Selected:");
+      fprintf(ficlog,"\n#****** Selected:");
       for(j=1;j<=cptcoveff;j++){ 
 	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }	
       fprintf(ficrest,"******\n");
       fprintf(ficlog,"******\n");
       printf("******\n");
@@ -10680,19 +11148,27 @@ Please run with mle=-1 to get a correct
 	fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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,"******\n");
       fprintf(ficrescveij,"******\n");
       
       fprintf(ficresvij,"\n#****** ");
+      /* pstamp(ficresvij); */
       for(j=1;j<=cptcoveff;j++) 
 	fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }	
       fprintf(ficresvij,"******\n");
       
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;
-      printf(" cvevsij combination#=%d, ",k);
-      fprintf(ficlog, " cvevsij combination#=%d, ",k);
-      cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+      printf(" cvevsij ");
+      fprintf(ficlog, " cvevsij ");
+      cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart, nres);
       printf(" end cvevsij \n ");
       fprintf(ficlog, " end cvevsij \n ");
       
@@ -10709,7 +11185,7 @@ Please run with mle=-1 to get a correct
 	cptcod= 0; /* To be deleted */
 	printf("varevsij vpopbased=%d \n",vpopbased);
 	fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);
-	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */
 	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
 	if(vpopbased==1)
 	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
@@ -10723,7 +11199,7 @@ Please run with mle=-1 to get a correct
 	printf("Computing age specific period (stable) prevalences in each health state \n");
 	fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
 	for(age=bage; age <=fage ;age++){
-	  prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+	  prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */
 	  if (vpopbased==1) {
 	    if(mobilav ==0){
 	      for(i=1; i<=nlstate;i++)
@@ -10760,11 +11236,11 @@ Please run with mle=-1 to get a correct
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_vector(epj,1,nlstate+1);
-      printf("done \n");fflush(stdout);
-      fprintf(ficlog,"done\n");fflush(ficlog);
+      printf("done selection\n");fflush(stdout);
+      fprintf(ficlog,"done selection\n");fflush(ficlog);
       
       /*}*/
-    } /* End k */
+    } /* End k selection */
 
     printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
@@ -10783,7 +11259,13 @@ Please run with mle=-1 to get a correct
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
     
-    for (k=1; k <= (int) pow(2,cptcoveff); k++){
+    i1=pow(2,cptcoveff);
+    if (cptcovn < 1){i1=1;}
+
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
+    for(k=1; k<=i1;k++){
+      if(TKresult[nres]!= k)
+	continue;
       fprintf(ficresvpl,"\n#****** ");
       printf("\n#****** ");
       fprintf(ficlog,"\n#****** ");
@@ -10792,13 +11274,18 @@ Please run with mle=-1 to get a correct
 	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
+      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+      }	
       fprintf(ficresvpl,"******\n");
       printf("******\n");
       fprintf(ficlog,"******\n");
       
       varpl=matrix(1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;
-      varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
+      varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres);
       free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
       /*}*/
     }
@@ -10847,8 +11334,14 @@ Please run with mle=-1 to get a correct
   free_ivector(ncodemaxwundef,1,NCOVMAX);
   free_ivector(Dummy,-1,NCOVMAX);
   free_ivector(Fixed,-1,NCOVMAX);
+  free_ivector(DummyV,1,NCOVMAX);
+  free_ivector(FixedV,1,NCOVMAX);
   free_ivector(Typevar,-1,NCOVMAX);
   free_ivector(Tvar,1,NCOVMAX);
+  free_ivector(TvarsQ,1,NCOVMAX);
+  free_ivector(TvarsQind,1,NCOVMAX);
+  free_ivector(TvarsD,1,NCOVMAX);
+  free_ivector(TvarsDind,1,NCOVMAX);
   free_ivector(TvarFD,1,NCOVMAX);
   free_ivector(TvarFDind,1,NCOVMAX);
   free_ivector(TvarF,1,NCOVMAX);