--- imach/src/imach.c	2016/08/25 10:50:18	1.236
+++ imach/src/imach.c	2016/09/08 16:07:27	1.250
@@ -1,6 +1,50 @@
-/* $Id: imach.c,v 1.236 2016/08/25 10:50:18 brouard Exp $
+/* $Id: imach.c,v 1.250 2016/09/08 16:07:27 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.250  2016/09/08 16:07:27  brouard
+  Summary: continue
+
+  Revision 1.249  2016/09/07 17:14:18  brouard
+  Summary: Starting values from frequencies
+
+  Revision 1.248  2016/09/07 14:10:18  brouard
+  *** empty log message ***
+
+  Revision 1.247  2016/09/02 11:11:21  brouard
+  *** empty log message ***
+
+  Revision 1.246  2016/09/02 08:49:22  brouard
+  *** empty log message ***
+
+  Revision 1.245  2016/09/02 07:25:01  brouard
+  *** empty log message ***
+
+  Revision 1.244  2016/09/02 07:17:34  brouard
+  *** empty log message ***
+
+  Revision 1.243  2016/09/02 06:45:35  brouard
+  *** empty log message ***
+
+  Revision 1.242  2016/08/30 15:01:20  brouard
+  Summary: Fixing a lots
+
+  Revision 1.241  2016/08/29 17:17:25  brouard
+  Summary: gnuplot problem in Back projection to fix
+
+  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 ***
 
@@ -908,12 +952,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.236 2016/08/25 10:50:18 brouard Exp $ */
+/* $Id: imach.c,v 1.250 2016/09/08 16:07:27 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.236 $ $Date: 2016/08/25 10:50:18 $"; 
+char fullversion[]="$Revision: 1.250 $ $Date: 2016/09/08 16:07:27 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1026,7 +1070,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 
 
@@ -1114,9 +1159,11 @@ int *TvarsQind;
 #define MAXRESULTLINES 10
 int nresult=0;
 int TKresult[MAXRESULTLINES];
-double Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
 int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */
 double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
+double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
 int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */
 
 /* int *TDvar; /\**< TDvar[1]=4,  TDvarF[2]=3, TDvar[3]=6  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
@@ -1140,6 +1187,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*/
@@ -1149,11 +1198,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;
 
@@ -2052,10 +2100,10 @@ void powell(double p[], double **xi, int
  void linmin(double p[], double xi[], int n, double *fret, 
 	      double (*func)(double [])); 
 #else 
- void linmin(double p[], double xi[], int n, double *fret, 
-						 double (*func)(double []),int *flat); 
+ 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;
@@ -2087,31 +2135,64 @@ 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  ");
+	fprintf(ficlog,"  + age*age  ");
+    }
+    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);
-    if(*iter <=3){
+    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 && *iter >1){
       tml = *localtime(&rcurr_time);
       strcpy(strcurr,asctime(&tml));
       rforecast_time=rcurr_time; 
       itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
-				strcurr[itmp-1]='\0';
+	strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       for(niterf=10;niterf<=30;niterf+=10){
-				rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
-				forecast_time = *localtime(&rforecast_time);
-				strcpy(strfor,asctime(&forecast_time));
-				itmp = strlen(strfor);
-				if(strfor[itmp-1]=='\n')
-					strfor[itmp-1]='\0';
-				printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
-				fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+	rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
+	forecast_time = *localtime(&rforecast_time);
+	strcpy(strfor,asctime(&forecast_time));
+	itmp = strlen(strfor);
+	if(strfor[itmp-1]=='\n')
+	  strfor[itmp-1]='\0';
+	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
       }
     }
     for (i=1;i<=n;i++) { /* For each direction i */
@@ -2165,7 +2246,8 @@ void powell(double p[], double **xi, int
 				/* printf("\n"); */
 				/* fprintf(ficlog,"\n"); */
 			}
-    if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
+    /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
+    if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
       /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
       /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */
       /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */
@@ -2211,7 +2293,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]; 
@@ -2410,7 +2492,7 @@ void powell(double p[], double **xi, int
       cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
       /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
     }
-    for (k=1; k<=cptcovage;k++){
+    for (k=1; k<=cptcovage;k++){  /* For product with age */
       if(Dummy[Tvar[Tage[k]]]){
 	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
       } else{
@@ -2418,9 +2500,21 @@ void powell(double p[], double **xi, int
       }
       /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
     }
-    for (k=1; k<=cptcovprod;k++){ /*  */
+    for (k=1; k<=cptcovprod;k++){ /* For product without age */
       /* printf("prevalim Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+      if(Dummy[Tvard[k][1]==0]){
+	if(Dummy[Tvard[k][2]==0]){
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
+	}
+      }else{
+	if(Dummy[Tvard[k][2]==0]){
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
+	}else{
+	  cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
+	}
+      }
     }
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -2479,7 +2573,7 @@ Earliest age to start was %d-%d=%d, ncvl
 
  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
- double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij)
+  double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres)
 {
   /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
      matrix by transitions matrix until convergence is reached with precision ftolpl */
@@ -2502,6 +2596,7 @@ Earliest age to start was %d-%d=%d, ncvl
   /* If we start from prlim again, prlim tends to a constant matrix */
 
   int i, ii,j,k;
+  int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **bmij();
@@ -2538,15 +2633,49 @@ Earliest age to start was %d-%d=%d, ncvl
     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])]; */
-      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("bprevalim 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<=cptcovn;k++) { */
+    /*   /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
+    /*   cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
+    /*   /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\/ */
+    /* } */
+    for (k=1; k<=nsq;k++) { /* For single varying covariates only */
+ 			/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
+      cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
+      /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
+    }
+    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; */
+    /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */
+    /*   /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */
+    /*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
+    for (k=1; k<=cptcovage;k++){  /* For product with age */
+      if(Dummy[Tvar[Tage[k]]]){
+	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]];
+	}
+      }
     }
-    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]);*/
@@ -2593,7 +2722,12 @@ Earliest age to start was %d-%d=%d, ncvl
     }
   } /* age loop */
     /* After some age loop it doesn't converge */
-  printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+  if(first){
+    first=1;
+    printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\
+Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+  }
+  fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
 Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
   free_vector(min,1,nlstate);
@@ -2726,7 +2860,8 @@ double **pmij(double **ps, double *cov,
 	/* }else */
 	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
       }else{
-	printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin);
+	;
+	/* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */
       }
     } /*End ii */
   } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
@@ -3109,7 +3244,7 @@ double func( double *x)
 	 Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
 	 to be observed in j being in i according to the model.
       */
-      ioffset=2+nagesqr+cptcovage;
+      ioffset=2+nagesqr ;
    /* 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)*/
@@ -3128,7 +3263,8 @@ double func( double *x)
       */
       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];
+	  /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
+	  cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
 	}
 	for (ii=1;ii<=nlstate+ndeath;ii++)
 	  for (j=1;j<=nlstate+ndeath;j++){
@@ -3142,7 +3278,10 @@ double func( double *x)
 	  if(nagesqr==1)
 	    cov[3]= agexact*agexact;  /* Should be changed here */
 	  for (kk=1; kk<=cptcovage;kk++) {
+	  if(!FixedV[Tvar[Tage[kk]]])
 	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+	  else
+	    cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
 	  }
 	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 		       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
@@ -3424,7 +3563,8 @@ double funcone( double *x)
   for(k=1; k<=nlstate; k++) ll[k]=0.;
   ioffset=0;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){
-    ioffset=2+nagesqr+cptcovage;
+    /* ioffset=2+nagesqr+cptcovage; */
+    ioffset=2+nagesqr;
     /* Fixed */
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
     /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
@@ -3451,46 +3591,51 @@ double funcone( double *x)
     for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */
     /* Wave varying (but not age varying) */
       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];
-			}
+	/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
+	cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
+      }
       /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
-				/* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
-				/* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
-				/* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
-				/* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
-				/* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
+      /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
+      /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
+      /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
+      /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
+      /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
       /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */
-			/* 	iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
-			/* 	/\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
-			/* 	cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */
+      /* 	iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
+      /* 	/\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
+      /* 	cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][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 (j=1;j<=nlstate+ndeath;j++){
+	  oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	  savm[ii][j]=(ii==j ? 1.0 : 0.0);
+	}
       
       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
-				/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
-					and mw[mi+1][i]. dh depends on stepm.*/
-				newm=savm;
-				agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
-				cov[2]=agexact;
-				if(nagesqr==1)
-					cov[3]= agexact*agexact;
-				for (kk=1; kk<=cptcovage;kk++) {
-					cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
-				}
-				/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
-				/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
-				out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
-										 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
-				/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
-				/* 	     1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
-				savm=oldm;
-				oldm=newm;
+      /* for(d=0; d<=0; d++){  /\* Delay between two effective waves Only one matrix to speed up*\/ */
+	/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+	  and mw[mi+1][i]. dh depends on stepm.*/
+	newm=savm;
+	agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;  /* Here d is needed */
+	cov[2]=agexact;
+	if(nagesqr==1)
+	  cov[3]= agexact*agexact;
+	for (kk=1; kk<=cptcovage;kk++) {
+	  if(!FixedV[Tvar[Tage[kk]]])
+	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+	  else
+	    cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
+	}
+	/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
+	/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+		     1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+	/* 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 */
       
       s1=s[mw[mi][i]][i];
@@ -3504,37 +3649,37 @@ double funcone( double *x)
        * is higher than the multiple of stepm and negative otherwise.
        */
       if( s2 > nlstate && (mle <5) ){  /* Jackson */
-				lli=log(out[s1][s2] - savm[s1][s2]);
+	lli=log(out[s1][s2] - savm[s1][s2]);
       } else if  ( s2==-1 ) { /* alive */
-				for (j=1,survp=0. ; j<=nlstate; j++) 
-					survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
-				lli= log(survp);
+	for (j=1,survp=0. ; j<=nlstate; j++) 
+	  survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+	lli= log(survp);
       }else if (mle==1){
-				lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+	lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
       } else if(mle==2){
-				lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
+	lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
       } else if(mle==3){  /* exponential inter-extrapolation */
-				lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
+	lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
       } else if (mle==4){  /* mle=4 no inter-extrapolation */
-				lli=log(out[s1][s2]); /* Original formula */
+	lli=log(out[s1][s2]); /* Original formula */
       } else{  /* mle=0 back to 1 */
-				lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
-				/*lli=log(out[s1][s2]); */ /* Original formula */
+	lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+	/*lli=log(out[s1][s2]); */ /* Original formula */
       } /* End of if */
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){
-				fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
+	fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \
-								num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
-								2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
-				for(k=1,llt=0.,l=0.; k<=nlstate; k++){
-					llt +=ll[k]*gipmx/gsw;
-					fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
-				}
-				fprintf(ficresilk," %10.6f\n", -llt);
+		num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
+		2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
+	for(k=1,llt=0.,l=0.; k<=nlstate; k++){
+	  llt +=ll[k]*gipmx/gsw;
+	  fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
+	}
+	fprintf(ficresilk," %10.6f\n", -llt);
       }
 	} /* end of wave */
 } /* end of individual */
@@ -3884,8 +4029,8 @@ double hessij( double x[], double **hess
       kmax=kmax+10;
     }
     if(kmax >=10 || firstime ==1){
-      printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
-      fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
+      printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
+      fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
       fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
     }
@@ -4041,12 +4186,12 @@ void pstamp(FILE *fichier)
 }
 
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
+void  freqsummary(char fileres[], double p[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
 		  int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],	\
 		  int firstpass,  int lastpass, int stepm, int weightopt, char model[])
-{  /* Some frequencies */
+{  /* Some frequencies as well as proposing some starting values */
   
-  int i, m, jk, j1, bool, z1,j, k, iv;
+  int i, m, jk, j1, bool, z1,j, k, iv, jj=0;
   int iind=0, iage=0;
   int mi; /* Effective wave */
   int first;
@@ -4073,7 +4218,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));
@@ -4083,55 +4228,57 @@ 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 
   */
+  dateintsum=0;
+  k2cpt=0;
 
+  for (j = 0; j <= cptcoveff; j+=cptcoveff){   
+  first=1;
   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;
     }
@@ -4141,26 +4288,28 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     /* 	meanqt[m][z1]=0.; */
     /*   } */
     /* } */
-		
-    dateintsum=0;
-    k2cpt=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(j !=0){
       if(anyvaryingduminmodel==0){ /* If All fixed covariates */
 	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++) {  
+	  for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
 	    /* if(Tvaraff[z1] ==-20){ */
 	    /* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
 	    /* }else  if(Tvaraff[z1] ==-10){ */
 	    /* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
 	    /* }else  */
-	    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;
+	    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */
+	      /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */
+	      bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
 	      /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
 		 bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
 		 j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
@@ -4169,16 +4318,20 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	  } /* end z1 */
 	} /* cptcovn > 0 */
       } /* end any */
+      }/* end j==0 */
       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(j!=0){
 	  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;
+		if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality. If covariate's 
+										  value is -1, we don't select. It differs from the 
+										  constant and age model which counts them. */
+		  bool=0; /* not selected */
 	      }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
 		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
 		  bool=0;
@@ -4186,6 +4339,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	      }
 	    }
 	  }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
+	  } /* end j==0 */
 	  /* 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]
@@ -4205,50 +4359,60 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 		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 */
+		/* if((int)agev[m][iind] == 55) */
+		/*   printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
 		/* 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;
+	    if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) {
+	      dateintsum=dateintsum+k2; /* on all covariates ?*/
 	      k2cpt++;
 	      /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
 	    }
-	  } /* end bool 2 */
+	  }else{
+	    bool=1;
+	  }/* 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 && j!=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> ");
@@ -4259,115 +4423,115 @@ Title=%s <br>Datafile=%s Firstpass=%d La
       }
     }
     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");
@@ -4384,8 +4548,47 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     }
     fprintf(ficresphtmfr,"</table>\n");
   } /* end selected combination of covariate j1 */
+  if(j==0){ /* We can estimate starting values from the occurences in each case */
+    printf("#Freqsummary\n");
+    fprintf(ficlog,"\n");
+    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(jj=1; jj <=ncovmodel; jj++){
+	    if(jj==1){
+	      printf("%12.7f ln(%12.1f/%12.1f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+	      fprintf(ficlog,"%12.7f ln(%12.1f/%12.1f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+	    }
+	    /* printf("%12.7f )", param[i][jj][k]); */
+	    /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
+      }
+    }
+    printf("#Freqsummary\n");
+    fprintf(ficlog,"\n");
+    for(jk=-1; jk <=nlstate+ndeath; jk++){
+      for(m=-1; m <=nlstate+ndeath; m++){
+	/* param[i]|j][k]= freq[jk][m][iagemax+3] */
+	  printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
+	  fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
+	/* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */
+	/*   printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
+	/*   fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
+	/* } */
+      }
+    } /* end loop jk */
+    printf("\n");
+    fprintf(ficlog,"\n");
+  } /* if j=0 */
+  } /* end j */
   dateintmean=dateintsum/k2cpt; 
-	
+  
   fclose(ficresp);
   fclose(ficresphtm);
   fclose(ficresphtmfr);
@@ -4720,173 +4923,173 @@ void  concatwav(int wav[], int **dh, int
 
 /*********** Tricode ****************************/
  void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum)
-{
-  /**< Uses cptcovn+2*cptcovprod as the number of covariates */
-  /*	  Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
-   * Boring subroutine which should only output nbcode[Tvar[j]][k]
-   * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable
-   * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);
-  */
+ {
+   /**< Uses cptcovn+2*cptcovprod as the number of covariates */
+   /*	  Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
+    * Boring subroutine which should only output nbcode[Tvar[j]][k]
+    * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable
+    * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);
+    */
 
-  int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
-  int modmaxcovj=0; /* Modality max of covariates j */
-  int cptcode=0; /* Modality max of covariates j */
-  int modmincovj=0; /* Modality min of covariates j */
+   int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
+   int modmaxcovj=0; /* Modality max of covariates j */
+   int cptcode=0; /* Modality max of covariates j */
+   int modmincovj=0; /* Modality min of covariates j */
 
 
-  /* cptcoveff=0;  */
-	/* *cptcov=0; */
+   /* cptcoveff=0;  */
+   /* *cptcov=0; */
  
-  for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
+   for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
 
-  /* Loop on covariates without age and products and no quantitative variable */
-  /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
-  for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
-    for (j=-1; (j < maxncov); j++) Ndum[j]=0;
-    if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
-      switch(Fixed[k]) {
-      case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
-				for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
-					ij=(int)(covar[Tvar[k]][i]);
-					/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
-					 * If product of Vn*Vm, still boolean *:
-					 * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
-					 * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
-					/* Finds for covariate j, n=Tvar[j] of Vn . ij is the
-						 modality of the nth covariate of individual i. */
-					if (ij > modmaxcovj)
-						modmaxcovj=ij; 
-					else if (ij < modmincovj) 
-						modmincovj=ij; 
-					if ((ij < -1) && (ij > NCOVMAX)){
-						printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
-						exit(1);
-					}else
-						Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
-					/*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
-					/*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
-					/* getting the maximum value of the modality of the covariate
-						 (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
-						 female ies 1, then modmaxcovj=1.
-					*/
-				} /* end for loop on individuals i */
-				printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
-				fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
-				cptcode=modmaxcovj;
-				/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
-				/*for (i=0; i<=cptcode; i++) {*/
-				for (j=modmincovj;  j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */
-					printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
-					fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
-					if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */
-						if( j != -1){
-							ncodemax[k]++;  /* ncodemax[k]= Number of modalities of the k th
-																 covariate for which somebody answered excluding 
-																 undefined. Usually 2: 0 and 1. */
-						}
-						ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th
-																		covariate for which somebody answered including 
-																		undefined. Usually 3: -1, 0 and 1. */
-					}	/* In fact  ncodemax[k]=2 (dichotom. variables only) but it could be more for
-						 * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
-				} /* Ndum[-1] number of undefined modalities */
+   /* Loop on covariates without age and products and no quantitative variable */
+   /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
+   for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
+     for (j=-1; (j < maxncov); j++) Ndum[j]=0;
+     if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
+       switch(Fixed[k]) {
+       case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
+	 for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
+	   ij=(int)(covar[Tvar[k]][i]);
+	   /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
+	    * If product of Vn*Vm, still boolean *:
+	    * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+	    * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
+	   /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
+	      modality of the nth covariate of individual i. */
+	   if (ij > modmaxcovj)
+	     modmaxcovj=ij; 
+	   else if (ij < modmincovj) 
+	     modmincovj=ij; 
+	   if ((ij < -1) && (ij > NCOVMAX)){
+	     printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+	     exit(1);
+	   }else
+	     Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
+	   /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
+	   /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
+	   /* getting the maximum value of the modality of the covariate
+	      (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
+	      female ies 1, then modmaxcovj=1.
+	   */
+	 } /* end for loop on individuals i */
+	 printf(" Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
+	 fprintf(ficlog," Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
+	 cptcode=modmaxcovj;
+	 /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
+	 /*for (i=0; i<=cptcode; i++) {*/
+	 for (j=modmincovj;  j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */
+	   printf("Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
+	   fprintf(ficlog, "Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
+	   if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */
+	     if( j != -1){
+	       ncodemax[k]++;  /* ncodemax[k]= Number of modalities of the k th
+				  covariate for which somebody answered excluding 
+				  undefined. Usually 2: 0 and 1. */
+	     }
+	     ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th
+				     covariate for which somebody answered including 
+				     undefined. Usually 3: -1, 0 and 1. */
+	   }	/* In fact  ncodemax[k]=2 (dichotom. variables only) but it could be more for
+		 * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
+	 } /* Ndum[-1] number of undefined modalities */
 			
-				/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
-				/* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */
-				/* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */
-				/* modmincovj=3; modmaxcovj = 7; */
-				/* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */
-				/* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */
-			  /*		 defining two dummy variables: variables V1_1 and V1_2.*/
- 	      /* nbcode[Tvar[j]][ij]=k; */
-	      /* nbcode[Tvar[j]][1]=0; */
-	      /* nbcode[Tvar[j]][2]=1; */
-	      /* nbcode[Tvar[j]][3]=2; */
-	      /* To be continued (not working yet). */
-	      ij=0; /* ij is similar to i but can jump over null modalities */
-				for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
-      	  if (Ndum[i] == 0) { /* If nobody responded to this modality k */
-	          break;
-	        }
-					ij++;
-					nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
-					cptcode = ij; /* New max modality for covar j */
-				} /* end of loop on modality i=-1 to 1 or more */
-				break;
-      case 1: /* Testing on varying covariate, could be simple and
-	       * should look at waves or product of fixed *
-	       * varying. No time to test -1, assuming 0 and 1 only */
-				ij=0;
-				for(i=0; i<=1;i++){
-					nbcode[Tvar[k]][++ij]=i;
-				}
-				break;
-      default:
-				break;
-      } /* end switch */
-    } /* end dummy test */
-    
-    /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
-    /* 	/\*recode from 0 *\/ */
-    /* 				     k is a modality. If we have model=V1+V1*sex  */
-    /* 				     then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
-    /* 				  But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
-    /* 	} */
-    /* 	/\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
-    /* 	if (ij > ncodemax[j]) { */
-    /* 	  printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */
-    /* 	  fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
-    /* 	  break; */
-    /* 	} */
-    /*   }  /\* end of loop on modality k *\/ */
-  } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
-  
-  for (k=-1; k< maxncov; k++) Ndum[k]=0; 
-  /* Look at fixed dummy (single or product) covariates to check empty modalities */
-  for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ 
-    /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
-    ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ 
-    Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */
-    /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1,  {2, 1, 1, 1, 2, 1, 1, 0, 0} */
-  } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */
-  
-  ij=0;
-  /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
-  for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
-    /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
-    /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */
-    if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */
-      /* If product not in single variable we don't print results */
-      /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
-      ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */
-      Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
-      Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
-      TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
-      if(Fixed[k]!=0)
-        anyvaryingduminmodel=1;
-			/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */
-			/*   Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */
-			/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */
-			/*   Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */
-			/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */
-			/*   Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */
-    } 
-  } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
-  /* ij--; */
-  /* cptcoveff=ij; /\*Number of total covariates*\/ */
-  *cptcov=ij; /*Number of total real effective covariates: effective
-							 * because they can be excluded from the model and real
-							 * if in the model but excluded because missing values, but how to get k from ij?*/
-  for(j=ij+1; j<= cptcovt; j++){
-    Tvaraff[j]=0;
-    Tmodelind[j]=0;
-  }
-  for(j=ntveff+1; j<= cptcovt; j++){
-    TmodelInvind[j]=0;
-  }
-  /* To be sorted */
-  ;
-}
+	 /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
+	 /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */
+	 /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */
+	 /* modmincovj=3; modmaxcovj = 7; */
+	 /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */
+	 /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */
+	 /*		 defining two dummy variables: variables V1_1 and V1_2.*/
+	 /* nbcode[Tvar[j]][ij]=k; */
+	 /* nbcode[Tvar[j]][1]=0; */
+	 /* nbcode[Tvar[j]][2]=1; */
+	 /* nbcode[Tvar[j]][3]=2; */
+	 /* To be continued (not working yet). */
+	 ij=0; /* ij is similar to i but can jump over null modalities */
+	 for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+	   if (Ndum[i] == 0) { /* If nobody responded to this modality k */
+	     break;
+	   }
+	   ij++;
+	   nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
+	   cptcode = ij; /* New max modality for covar j */
+	 } /* end of loop on modality i=-1 to 1 or more */
+	 break;
+       case 1: /* Testing on varying covariate, could be simple and
+		* should look at waves or product of fixed *
+		* varying. No time to test -1, assuming 0 and 1 only */
+	 ij=0;
+	 for(i=0; i<=1;i++){
+	   nbcode[Tvar[k]][++ij]=i;
+	 }
+	 break;
+       default:
+	 break;
+       } /* end switch */
+     } /* end dummy test */
+    
+     /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
+     /* 	/\*recode from 0 *\/ */
+     /* 				     k is a modality. If we have model=V1+V1*sex  */
+     /* 				     then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
+     /* 				  But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
+     /* 	} */
+     /* 	/\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
+     /* 	if (ij > ncodemax[j]) { */
+     /* 	  printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */
+     /* 	  fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
+     /* 	  break; */
+     /* 	} */
+     /*   }  /\* end of loop on modality k *\/ */
+   } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
+  
+   for (k=-1; k< maxncov; k++) Ndum[k]=0; 
+   /* Look at fixed dummy (single or product) covariates to check empty modalities */
+   for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ 
+     /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
+     ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ 
+     Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */
+     /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1,  {2, 1, 1, 1, 2, 1, 1, 0, 0} */
+   } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */
+  
+   ij=0;
+   /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
+   for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+     /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
+     /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */
+     if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */
+       /* If product not in single variable we don't print results */
+       /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+       ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */
+       Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
+       Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
+       TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
+       if(Fixed[k]!=0)
+	 anyvaryingduminmodel=1;
+       /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */
+       /*   Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */
+       /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */
+       /*   Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */
+       /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */
+       /*   Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */
+     } 
+   } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
+   /* ij--; */
+   /* cptcoveff=ij; /\*Number of total covariates*\/ */
+   *cptcov=ij; /*Number of total real effective covariates: effective
+		* because they can be excluded from the model and real
+		* if in the model but excluded because missing values, but how to get k from ij?*/
+   for(j=ij+1; j<= cptcovt; j++){
+     Tvaraff[j]=0;
+     Tmodelind[j]=0;
+   }
+   for(j=ntveff+1; j<= cptcovt; j++){
+     TmodelInvind[j]=0;
+   }
+   /* To be sorted */
+   ;
+ }
 
 
 /*********** Health Expectancies ****************/
@@ -4901,7 +5104,7 @@ void  concatwav(int wav[], int **dh, int
   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++){
@@ -5266,6 +5469,14 @@ void  concatwav(int wav[], int **dh, int
    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);
@@ -5335,7 +5546,7 @@ void  concatwav(int wav[], int **dh, int
 	 xp[i] = x[i] + (i==theta ?delti[theta]:0);
        }
 			
-       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);
+       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
 			
        if (popbased==1) {
 	 if(mobilav ==0){
@@ -5367,7 +5578,7 @@ void  concatwav(int wav[], int **dh, int
        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, nresult);
+       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres);
 			
        if (popbased==1) {
 	 if(mobilav ==0){
@@ -5444,7 +5655,7 @@ void  concatwav(int wav[], int **dh, int
      /* end ppptj */
      /*  x centered again */
 		
-     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);
+     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
 		
      if (popbased==1) {
        if(mobilav ==0){
@@ -5539,7 +5750,9 @@ void  concatwav(int wav[], int **dh, int
   
   pstamp(ficresvpl);
   fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
-  fprintf(ficresvpl,"# Age");
+  fprintf(ficresvpl,"# Age ");
+  if(nresult >=1)
+    fprintf(ficresvpl," Result# ");
   for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);
   fprintf(ficresvpl,"\n");
@@ -5625,6 +5838,8 @@ void  concatwav(int wav[], int **dh, int
       varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
 
     fprintf(ficresvpl,"%.0f ",age );
+    if(nresult >=1)
+      fprintf(ficresvpl,"%d ",nres );
     for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
     fprintf(ficresvpl,"\n");
@@ -5986,11 +6201,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"));
@@ -6025,16 +6242,27 @@ void printinghtml(char fileresu[], char
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
    jj1=0;
-   for(k1=1; k1<=m;k1++){
+
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
+     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]){
@@ -6044,51 +6272,51 @@ void printinghtml(char fileresu[], char
        }
      }
      /* aij, bij */
-     fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
-<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+     fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
+<img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
      /* Pij */
-     fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
-<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
+     fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \
+<img src=\"%s_%d-2-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);     
      /* Quasi-incidences */
      fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
  before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
  incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \
-divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
-<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
+divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-%d.svg</a><br> \
+<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }
      /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
  Or probability to survive in various states (1 to %d) being in state %d at different ages.	\
- <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
+ <a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }
      /* Period (stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }
      if(backcast==1){
        /* Period (stable) back prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
+	 fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
        }
      }
      if(prevfcast==1){
        /* Projection of prevalence up to period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }
      }
 	 
      for(cpt=1; cpt<=nlstate;cpt++) {
-       fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
-<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
+       fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \
+<img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
      }
      /* } /\* end i1 *\/ */
    }/* End k1 */
@@ -6146,13 +6374,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]){
@@ -6162,17 +6399,18 @@ See page 'Matrix of variance-covariance
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\
-<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);  
+prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\
+<img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\
-<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
+ observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\
+<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
      /* } /\* end i1 *\/ */
    }/* End k1 */
+  }/* End nres */
    fprintf(fichtm,"</ul>");
    fflush(fichtm);
 }
@@ -6182,7 +6420,7 @@ void printinggnuplot(char fileresu[], ch
 
   char dirfileres[132],optfileres[132];
   char gplotcondition[132];
-  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0,l=0;
+  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
@@ -6231,302 +6469,303 @@ 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 ; k1 ++) /* For each valid combination of covariate */
+  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 */
+	/* 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-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
+	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
+      
+	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==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
+	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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); 
+	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, nres in 2 to be fixed */
+	  if(cptcoveff ==0){
+	    fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3",	 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' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+			2+cptcoveff*2+(cpt-1),  cpt );  /* 4 or 6 ?*/
+	      }else{
+		fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+		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 ++){  
+    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       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 */
+      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]; /* 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);
+	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]);
-      }	
-      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\" \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]);
-	      kl++;
-	    }
-	  } /* end covariate */
-	} /* end if no covariate */
-      } /* end if backcast */
-      fprintf(ficgp,"\nset out \n");
-    } /* 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 */
-      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-      vlv= nbcode[Tvaraff[k]][lv];
-      fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-    }
-    for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-      printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-    }
-    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 end 2 eme*/
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
+      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-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* 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]!= k)
-      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); 
+  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\
+	/*       k=2+nlstate*(2*cpt-2); */
+	k=2+(nlstate+1)*(cpt-1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
+	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 (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-  for(nres=1; nres <= nresult; nres++) /* For each resultline */
-    if(TKresult[nres]!= k)
-      continue;
-
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
-	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); 
+  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-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+	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 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); 
+  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-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+	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);
@@ -6538,17 +6777,18 @@ 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); 
 	continue;
       }
       
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
+      fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
       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)
@@ -6569,7 +6809,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 */
@@ -6580,17 +6823,18 @@ 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); 
 	  continue;
 	}
 	
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
 	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)
@@ -6616,7 +6860,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  */
@@ -6627,6 +6874,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); 
@@ -6634,11 +6884,9 @@ 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,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
 	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 */   
@@ -6705,8 +6953,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++){
@@ -6723,7 +6971,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");
@@ -6738,11 +6987,20 @@ 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,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
+    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-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres);
       fprintf(ficgp,"\nset ter svg size 640, 480 ");
       if (ng==1){
 	fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
@@ -6782,18 +7040,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 */
@@ -6811,8 +7101,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++;
@@ -7076,8 +7366,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
       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 */
-      printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-      fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+      fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
     }
     fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){ 
@@ -7686,14 +7975,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);
@@ -7723,87 +8034,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;
@@ -8020,7 +8331,7 @@ int decoderesult ( char resultline[], in
     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[2]=5 == Tvarsel[1]=4   */
+	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;
@@ -8036,7 +8347,7 @@ int decoderesult ( char resultline[], in
     match=0;
     for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       if(Typevar[k1]==0){ /* Single */
-	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=5 == Tvarsel[1]=4   */
+	if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
 	  resultmodel[k1]=k2;  /* resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
 	  ++match;
 	}
@@ -8063,20 +8374,26 @@ int decoderesult ( char resultline[], in
   /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */
   /*    7 0 1 1 */
   /*    8 1 1 1 */
+  /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */
+  /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */
+  /* V5*age V5 known which value for nres?  */
+  /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */
   for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */
     if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */
-      k3= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
+      k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */
       k2=(int)Tvarsel[k3]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
       k+=Tvalsel[k3]*pow(2,k4);  /*  Tvalsel[1]=1  */
-      Tresult[nres][k4+1]=Tvalsel[k3];
-      Tvresult[nres][k4+1]=(int)Tvarsel[k3];
+      Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1)  Tresult[nres][2]=0(V3=0) */
+      Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+      Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
       printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4);
       k4++;;
     }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */
       k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
-      Tqresult[nres][k4q+1]=Tvalsel[k3q];
-      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q];
+      Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
+      Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
+      Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
       k4q++;;
     }
@@ -8102,7 +8419,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];
@@ -8308,7 +8625,7 @@ 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;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;
@@ -8333,7 +8650,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       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+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++;
@@ -8347,7 +8664,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       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){/* Only simple time varying variables */
+    }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */
       Fixed[k]= 1;
       Dummy[k]= 0;
       ntveff++; /* Only simple time varying dummy variable */
@@ -8358,7 +8675,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarsDind[nsd]=k;
       ncovv++; /* Only simple time varying variables */
       TvarV[ncovv]=Tvar[k];
-      TvarVind[ncovv]=k;
+      TvarVind[ncovv]=k; /* TvarVind[2]=2  TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       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);
@@ -8374,7 +8691,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarsQ[nsq]=Tvar[k];
       TvarsQind[nsq]=k;
       TvarV[ncovv]=Tvar[k];
-      TvarVind[ncovv]=k;
+      TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
       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 */
@@ -8386,167 +8703,167 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       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]);
@@ -8588,7 +8905,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
 }
 
 int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
-{
+{/* Check ages at death */
   int i, m;
   int firstone=0;
   
@@ -8923,61 +9240,63 @@ int prevalence_limit(double *p, double *
   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++){
-    if(TKresult[nres]!= 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;
 
-  /* 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(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;
-    }
+      /* 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");
+      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, 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 */
+      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;
 }
 
@@ -9018,69 +9337,69 @@ int back_prevalence_limit(double *p, dou
   i1=pow(2,cptcoveff);
   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("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;
-    }
+  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, nres);
+	}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,nres);
+	}
+	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); */
   
@@ -9172,7 +9491,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;
@@ -9200,48 +9519,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 */
 
@@ -9624,6 +9949,12 @@ int main(int argc, char *argv[])
   delti=delti3[1][1];
   /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */
+/* We could also provide initial parameters values giving by simple logistic regression 
+ * only one way, that is without matrix product. We will have nlstate maximizations */
+      /* for(i=1;i<nlstate;i++){ */
+      /* 	/\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
+      /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
+      /* } */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
@@ -9632,6 +9963,11 @@ int main(int argc, char *argv[])
     fclose (ficlog);
     goto end;
     exit(0);
+  } else if(mle==-2) { /* Guessing from means */
+    prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
+    printf(" You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
+    fprintf(ficlog," You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
+   
   }  else if(mle==-5) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
@@ -9955,7 +10291,7 @@ Please run with mle=-1 to get a correct
   */
 
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
-  /* */
+  /* Concatenates waves */
  
   free_vector(moisdc,1,n);
   free_vector(andc,1,n);
@@ -10136,7 +10472,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   /* Calculates basic frequencies. Computes observed prevalence at single age 
 		 and for any valid combination of covariates
      and prints on file fileres'p'. */
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
+  freqsummary(fileres, p, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
 	      firstpass, lastpass,  stepm,  weightopt, model);
 
   fprintf(fichtm,"\n");
@@ -10158,9 +10494,9 @@ Interval (in months) between two waves:
   /* For mortality only */
   if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM); 
-		for(i=1;i<=NDIM;i++)
-			for(j=1;j<=NDIM;j++)
-				ximort[i][j]=0.;
+    for(i=1;i<=NDIM;i++)
+      for(j=1;j<=NDIM;j++)
+	ximort[i][j]=0.;
     /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
     cens=ivector(1,n);
     ageexmed=vector(1,n);
@@ -10396,6 +10732,10 @@ Please run with mle=-1 to get a correct
     printf("\n");
     if(mle>=1){ /* Could be 1 or 2, Real Maximization */
       /* mlikeli uses func not funcone */
+      /* for(i=1;i<nlstate;i++){ */
+      /* 	/\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
+      /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
+      /* } */
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
     }
     if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
@@ -10646,6 +10986,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);
@@ -10662,6 +11003,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);
     
@@ -10680,14 +11022,17 @@ 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);
       }
       nresult++; /* Sum of resultlines */
@@ -10698,12 +11043,16 @@ Please run with mle=-1 to get a correct
 	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
@@ -10713,7 +11062,7 @@ Please run with mle=-1 to get a correct
 	break;
       else{ /* Processess output results for this combination of covariate values */
       }				   
-    }
+    } /* end while */
 
 
     
@@ -10848,6 +11197,8 @@ 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);
  		
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}
@@ -10929,9 +11280,9 @@ Please run with mle=-1 to get a correct
     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:");
+      printf("\n#****** Result for:");
+      fprintf(ficrest,"\n#****** Result for:");
+      fprintf(ficlog,"\n#****** Result for:");
       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)]);
@@ -10960,6 +11311,7 @@ Please run with mle=-1 to get a correct
       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 */
@@ -11137,6 +11489,8 @@ 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);