--- imach/src/imach.c	2016/02/15 23:35:36	1.221
+++ imach/src/imach.c	2016/02/19 09:23:35	1.223
@@ -1,6 +1,12 @@
-/* $Id: imach.c,v 1.221 2016/02/15 23:35:36 brouard Exp $
+/* $Id: imach.c,v 1.223 2016/02/19 09:23:35 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.223  2016/02/19 09:23:35  brouard
+  Summary: temporary
+
+  Revision 1.222  2016/02/17 08:14:50  brouard
+  Summary: Probably last 0.98 stable version 0.98r6
+
   Revision 1.221  2016/02/15 23:35:36  brouard
   Summary: minor bug
 
@@ -739,7 +745,7 @@ Back prevalence and projections:
 /* #define DEBUGLINMIN */
 /* #define DEBUGHESS */
 #define DEBUGHESSIJ
-#define LINMINORIGINAL  /* Don't use loop on scale in linmin (accepting nan)*/
+/* #define LINMINORIGINAL  /\* Don't use loop on scale in linmin (accepting nan)*\/ */
 #define POWELL /* Instead of NLOPT */
 #define POWELLF1F3 /* Skip test */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
@@ -831,12 +837,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.221 2016/02/15 23:35:36 brouard Exp $ */
+/* $Id: imach.c,v 1.223 2016/02/19 09:23:35 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.221 $ $Date: 2016/02/15 23:35:36 $"; 
+char fullversion[]="$Revision: 1.223 $ $Date: 2016/02/19 09:23:35 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -854,6 +860,7 @@ int npar=NPARMAX;
 int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */
 int ncovmodel=0, ncovcol=0;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int  nqv=0, ntv=0, nqtv=0;    /* Total number of quantitative variables, time variable (dummy), quantitative and time variable */ 
 int popbased=0;
 
 int *wav; /* Number of waves for this individuual 0 is possible */
@@ -992,6 +999,9 @@ double *agedc;
 double  **covar; /**< covar[j,i], value of jth covariate for individual i,
 		  * covar=matrix(0,NCOVMAX,1,n); 
 		  * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
+double ***cotvar; /* Time varying covariate */
+double ***cotqvar; /* Time varying quantitative covariate */
+double **coqvar; /* Fixed quantitative covariate */
 double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
 int *Tage;
@@ -2331,65 +2341,65 @@ double **pmij(double **ps, double *cov,
   /*double t34;*/
   int i,j, nc, ii, jj;
 
-	for(i=1; i<= nlstate; i++){
-		for(j=1; j<i;j++){
-			for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-				/*lnpijopii += param[i][j][nc]*cov[nc];*/
-				lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
-				/* 	 printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-			}
-			ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-			/* 	printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-		}
-		for(j=i+1; j<=nlstate+ndeath;j++){
-			for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-				/*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
-				lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
-				/* 	  printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
-			}
-			ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-		}
-	}
+  for(i=1; i<= nlstate; i++){
+    for(j=1; j<i;j++){
+      for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	/*lnpijopii += param[i][j][nc]*cov[nc];*/
+	lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+	/* 	 printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+      }
+      ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+      /* 	printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+    }
+    for(j=i+1; j<=nlstate+ndeath;j++){
+      for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	/*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+	lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+	/* 	  printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+      }
+      ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+    }
+  }
   
-	for(i=1; i<= nlstate; i++){
-		s1=0;
-		for(j=1; j<i; j++){
-			s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-			/*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-		}
-		for(j=i+1; j<=nlstate+ndeath; j++){
-			s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-			/*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-		}
-		/* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
-		ps[i][i]=1./(s1+1.);
-		/* Computing other pijs */
-		for(j=1; j<i; j++)
-			ps[i][j]= exp(ps[i][j])*ps[i][i];
-		for(j=i+1; j<=nlstate+ndeath; j++)
-			ps[i][j]= exp(ps[i][j])*ps[i][i];
-		/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
-	} /* end i */
+  for(i=1; i<= nlstate; i++){
+    s1=0;
+    for(j=1; j<i; j++){
+      s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+      /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+    }
+    for(j=i+1; j<=nlstate+ndeath; j++){
+      s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+      /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+    }
+    /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+    ps[i][i]=1./(s1+1.);
+    /* Computing other pijs */
+    for(j=1; j<i; j++)
+      ps[i][j]= exp(ps[i][j])*ps[i][i];
+    for(j=i+1; j<=nlstate+ndeath; j++)
+      ps[i][j]= exp(ps[i][j])*ps[i][i];
+    /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+  } /* end i */
   
-	for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
-		for(jj=1; jj<= nlstate+ndeath; jj++){
-			ps[ii][jj]=0;
-			ps[ii][ii]=1;
-		}
-	}
+  for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+    for(jj=1; jj<= nlstate+ndeath; jj++){
+      ps[ii][jj]=0;
+      ps[ii][ii]=1;
+    }
+  }
   
   
-	/* for(ii=1; ii<= nlstate+ndeath; ii++){ */
-	/*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
-	/* 	printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
-	/*   } */
-	/*   printf("\n "); */
-	/* } */
-	/* printf("\n ");printf("%lf ",cov[2]);*/
-	/*
-		for(i=1; i<= npar; i++) printf("%f ",x[i]);
+  /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+  /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
+  /* 	printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+  /*   } */
+  /*   printf("\n "); */
+  /* } */
+  /* printf("\n ");printf("%lf ",cov[2]);*/
+  /*
+    for(i=1; i<= npar; i++) printf("%f ",x[i]);
 		goto end;*/
-	return ps;
+  return ps;
 }
 
 /*************** backward transition probabilities ***************/ 
@@ -2398,60 +2408,60 @@ double **pmij(double **ps, double *cov,
 /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
  double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, int ij )
 {
-	/* Computes the backward probability at age agefin and covariate ij
-	 * and returns in **ps as well as **bmij.
-	 */
+  /* Computes the backward probability at age agefin and covariate ij
+   * and returns in **ps as well as **bmij.
+   */
   int i, ii, j,k;
-
-	double **out, **pmij();
-	double sumnew=0.;
+  
+  double **out, **pmij();
+  double sumnew=0.;
   double agefin;
-
-	double **dnewm, **dsavm, **doldm;
-	double **bbmij;
-
+  
+  double **dnewm, **dsavm, **doldm;
+  double **bbmij;
+  
   doldm=ddoldms; /* global pointers */
-	dnewm=ddnewms;
-	dsavm=ddsavms;
-
-	agefin=cov[2];
-	/* bmij *//* age is cov[2], ij is included in cov, but we need for
-		 the observed prevalence (with this covariate ij) */
-	dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate);
-	/* We do have the matrix Px in savm  and we need pij */
-	for (j=1;j<=nlstate+ndeath;j++){
-		sumnew=0.; /* w1 p11 + w2 p21 only on live states */
-		for (ii=1;ii<=nlstate;ii++){
-			sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij];
-		} /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
-		for (ii=1;ii<=nlstate+ndeath;ii++){
-			if(sumnew >= 1.e-10){
-				/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
-				/* 	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
-				/* }else if(agefin >= agemaxpar+stepm/YEARM){ */
-				/* 	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
-				/* }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);
-			}
-		} /*End ii */
-	} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
-		/* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */
-	bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */
-	/* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */
-	/* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
-	/* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
-	/* left Product of this matrix by diag matrix of prevalences (savm) */
-	for (j=1;j<=nlstate+ndeath;j++){
-		for (ii=1;ii<=nlstate+ndeath;ii++){
-			dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
-		}
-	} /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */
-	ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */
-	/* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
-	/* end bmij */
-	return ps; 
+  dnewm=ddnewms;
+  dsavm=ddsavms;
+  
+  agefin=cov[2];
+  /* bmij *//* age is cov[2], ij is included in cov, but we need for
+     the observed prevalence (with this covariate ij) */
+  dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate);
+  /* We do have the matrix Px in savm  and we need pij */
+  for (j=1;j<=nlstate+ndeath;j++){
+    sumnew=0.; /* w1 p11 + w2 p21 only on live states */
+    for (ii=1;ii<=nlstate;ii++){
+      sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij];
+    } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
+    for (ii=1;ii<=nlstate+ndeath;ii++){
+      if(sumnew >= 1.e-10){
+	/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
+	/* 	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+	/* }else if(agefin >= agemaxpar+stepm/YEARM){ */
+	/* 	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+	/* }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);
+      }
+    } /*End ii */
+  } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
+  /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */
+  bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */
+  /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */
+  /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
+  /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
+  /* left Product of this matrix by diag matrix of prevalences (savm) */
+  for (j=1;j<=nlstate+ndeath;j++){
+    for (ii=1;ii<=nlstate+ndeath;ii++){
+      dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
+    }
+  } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */
+  ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */
+  /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
+  /* end bmij */
+  return ps; 
 }
 /*************** transition probabilities ***************/ 
 
@@ -2657,7 +2667,7 @@ double ***hpxij(double ***po, int nhstep
 
 /************* Higher Back Matrix Product ***************/
 /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
- double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
+double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
 {
   /* Computes the transition matrix starting at age 'age' over
      'nhstepm*hstepm*stepm' months (i.e. until
@@ -2669,16 +2679,16 @@ double ***hpxij(double ***po, int nhstep
      Model is determined by parameters x and covariates have to be
      included manually here.
 
-     */
+  */
 
   int i, j, d, h, k;
   double **out, cov[NCOVMAX+1];
   double **newm;
   double agexact;
   double agebegin, ageend;
-	double **oldm, **savm;
+  double **oldm, **savm;
 
-	oldm=oldms;savm=savms;
+  oldm=oldms;savm=savms;
   /* Hstepm could be zero and should return the unit matrix */
   for (i=1;i<=nlstate+ndeath;i++)
     for (j=1;j<=nlstate+ndeath;j++){
@@ -2695,27 +2705,27 @@ double ***hpxij(double ***po, int nhstep
       /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
       cov[2]=agexact;
       if(nagesqr==1)
-				cov[3]= agexact*agexact;
+	cov[3]= agexact*agexact;
       for (k=1; k<=cptcovn;k++)
-				cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
-			/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
-				/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-				cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
-			/* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+	/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+      /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-				cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
-			/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
 			
 			
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
       /* Careful transposed matrix */
-			/* age is in cov[2] */
+      /* age is in cov[2] */
       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
-			/* 						 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
+      /* 						 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
       out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\
-									 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+		   1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
       /* if((int)age == 70){ */
       /* 	printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
       /* 	for(i=1; i<=nlstate+ndeath; i++) { */
@@ -2735,12 +2745,12 @@ double ***hpxij(double ***po, int nhstep
     }
     for(i=1; i<=nlstate+ndeath; i++)
       for(j=1;j<=nlstate+ndeath;j++) {
-				po[i][j][h]=newm[i][j];
-				/*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+	po[i][j][h]=newm[i][j];
+	/*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
       }
     /*printf("h=%d ",h);*/
   } /* end h */
-	/*     printf("\n H=%d \n",h); */
+  /*     printf("\n H=%d \n",h); */
   return po;
 }
 
@@ -2769,11 +2779,13 @@ double ***hpxij(double ***po, int nhstep
 double func( double *x)
 {
   int i, ii, j, k, mi, d, kk;
+	int ioffset;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
   int s1, s2;
+  int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */
   double bbh, survp;
   long ipmx;
   double agexact;
@@ -2793,17 +2805,38 @@ double func( double *x)
   if(mle==1){
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       /* Computes the values of the ncovmodel covariates of the model
-	 depending if the covariates are fixed or variying (age dependent) and stores them in cov[]
+	 depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
 	 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;
       for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */
-	  cov[2+nagesqr+k]=covar[Tvar[k]][i];
+	 cov[++ioffset]=covar[Tvar[k]][i];
       }
+      for(iqv=1; iqv <= nqv; iqv++){ /* Varying quantitatives covariates */
+	/* cov[2+nagesqr+cptcovn+iqv]=varq[mw[mi+1][i]][iqv][i]; */
+      }
+
       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
 	 is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] 
 	 has been calculated etc */
+      /* For an individual i, wav[i] gives the number of effective waves */
+      /* We compute the contribution to Likelihood of each effective transition
+	 mw[mi][i] is real wave of the mi th effectve wave */
+      /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
+	s2=s[mw[mi+1][i]][i];
+	And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
+	But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
+	meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
+      */
       for(mi=1; mi<= wav[i]-1; mi++){
+	for(itv=1; itv <= ntv; itv++){ /* Varying dummy covariates */
+	  cov[++ioffset]=cotvar[mw[mi+1][i]][itv][i];
+	}
+	for(iqtv=1; iqtv <= nqtv; iqtv++){ /* Varying quantitatives covariates */
+	  /* cov[2+nagesqr+cptcovn+nqv+ntv+iqtv]=varq[mw[mi+1][i]][iqtv][i]; */
+	}
+	ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv;
 	for (ii=1;ii<=nlstate+ndeath;ii++)
 	  for (j=1;j<=nlstate+ndeath;j++){
 	    oldm[ii][j]=(ii==j ? 1.0 : 0.0);
@@ -2814,7 +2847,7 @@ double func( double *x)
 	  agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
 	  cov[2]=agexact;
 	  if(nagesqr==1)
-	    cov[3]= agexact*agexact;
+	    cov[3]= agexact*agexact;  /* Should be changed here */
 	  for (kk=1; kk<=cptcovage;kk++) {
 	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
 	  }
@@ -3995,97 +4028,97 @@ Title=%s <br>Datafile=%s Firstpass=%d La
  }
 
 /************ Prevalence ********************/
-void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
-{  
-  /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
-     in each health status at the date of interview (if between dateprev1 and dateprev2).
-     We still use firstpass and lastpass as another selection.
-  */
+ void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
+ {  
+   /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+      in each health status at the date of interview (if between dateprev1 and dateprev2).
+      We still use firstpass and lastpass as another selection.
+   */
  
-  int i, m, jk, j1, bool, z1,j;
-  int mi; /* Effective wave */
-  int iage;
-  double agebegin, ageend;
-
-  double **prop;
-  double posprop; 
-  double  y2; /* in fractional years */
-  int iagemin, iagemax;
-  int first; /** to stop verbosity which is redirected to log file */
-
-  iagemin= (int) agemin;
-  iagemax= (int) agemax;
-  /*pp=vector(1,nlstate);*/
-  prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); 
-  /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
-  j1=0;
-  
-  /*j=cptcoveff;*/
-  if (cptcovn<1) {j=1;ncodemax[1]=1;}
-  
-  first=1;
-  for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
-    for (i=1; i<=nlstate; i++)  
-      for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
-				prop[i][iage]=0.0;
-    
-    for (i=1; i<=imx; i++) { /* Each individual */
-      bool=1;
-      if  (cptcovn>0) {  /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
-				for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
-					if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
-						bool=0;
-      } 
-      if (bool==1) { /* For this combination of covariates values, this individual fits */
-				/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
-				for(mi=1; mi<wav[i];mi++){
-					m=mw[mi][i];
-					agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
-					/* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
-					if(m >=firstpass && m <=lastpass){
-						y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
-						if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
-							if(agev[m][i]==0) agev[m][i]=iagemax+1;
-							if(agev[m][i]==1) agev[m][i]=iagemax+2;
-							if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
-								printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); 
-								exit(1);
-							}
-							if (s[m][i]>0 && s[m][i]<=nlstate) { 
-								/*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
-								prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
-								prop[s[m][i]][iagemax+3] += weight[i]; 
-							} /* end valid statuses */ 
-						} /* end selection of dates */
-					} /* end selection of waves */
-				} /* end effective waves */
-      } /* end bool */
-    }
-    for(i=iagemin; i <= iagemax+3; i++){  
-      for(jk=1,posprop=0; jk <=nlstate ; jk++) { 
-				posprop += prop[jk][i]; 
-      } 
-      
-      for(jk=1; jk <=nlstate ; jk++){	    
-				if( i <=  iagemax){ 
-					if(posprop>=1.e-5){ 
-						probs[i][jk][j1]= prop[jk][i]/posprop;
-					} else{
-						if(first==1){
-							first=0;
-							printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
-						}
-					}
-				} 
-      }/* end jk */ 
-    }/* end i */ 
-    /*} *//* end i1 */
-  } /* end j1 */
-  
-  /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
-  /*free_vector(pp,1,nlstate);*/
-  free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
-}  /* End of prevalence */
+   int i, m, jk, j1, bool, z1,j;
+   int mi; /* Effective wave */
+   int iage;
+   double agebegin, ageend;
+
+   double **prop;
+   double posprop; 
+   double  y2; /* in fractional years */
+   int iagemin, iagemax;
+   int first; /** to stop verbosity which is redirected to log file */
+
+   iagemin= (int) agemin;
+   iagemax= (int) agemax;
+   /*pp=vector(1,nlstate);*/
+   prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); 
+   /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
+   j1=0;
+  
+   /*j=cptcoveff;*/
+   if (cptcovn<1) {j=1;ncodemax[1]=1;}
+  
+   first=1;
+   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
+     for (i=1; i<=nlstate; i++)  
+       for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
+	 prop[i][iage]=0.0;
+    
+     for (i=1; i<=imx; i++) { /* Each individual */
+       bool=1;
+       if  (cptcovn>0) {  /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+	 for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
+	   if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
+	     bool=0;
+       } 
+       if (bool==1) { /* For this combination of covariates values, this individual fits */
+	 /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
+	 for(mi=1; mi<wav[i];mi++){
+	   m=mw[mi][i];
+	   agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+	   /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
+	   if(m >=firstpass && m <=lastpass){
+	     y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
+	     if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
+	       if(agev[m][i]==0) agev[m][i]=iagemax+1;
+	       if(agev[m][i]==1) agev[m][i]=iagemax+2;
+	       if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
+		 printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); 
+		 exit(1);
+	       }
+	       if (s[m][i]>0 && s[m][i]<=nlstate) { 
+		 /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
+		 prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
+		 prop[s[m][i]][iagemax+3] += weight[i]; 
+	       } /* end valid statuses */ 
+	     } /* end selection of dates */
+	   } /* end selection of waves */
+	 } /* end effective waves */
+       } /* end bool */
+     }
+     for(i=iagemin; i <= iagemax+3; i++){  
+       for(jk=1,posprop=0; jk <=nlstate ; jk++) { 
+	 posprop += prop[jk][i]; 
+       } 
+      
+       for(jk=1; jk <=nlstate ; jk++){	    
+	 if( i <=  iagemax){ 
+	   if(posprop>=1.e-5){ 
+	     probs[i][jk][j1]= prop[jk][i]/posprop;
+	   } else{
+	     if(first==1){
+	       first=0;
+	       printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
+	     }
+	   }
+	 } 
+       }/* end jk */ 
+     }/* end i */ 
+     /*} *//* end i1 */
+   } /* end j1 */
+  
+   /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
+   /*free_vector(pp,1,nlstate);*/
+   free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+ }  /* End of prevalence */
 
 /************* Waves Concatenation ***************/
 
@@ -4116,36 +4149,36 @@ void  concatwav(int wav[], int **dh, int
     m=firstpass;
     while(s[m][i] <= nlstate){  /* a live state */
       if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
-	mw[++mi][i]=m;
+				mw[++mi][i]=m;
       }
       if(m >=lastpass){
-	if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
-	  if(firsthree == 0){
-	    printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
-	    firsthree=1;
-	  }
-	  fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
-	  mw[++mi][i]=m;
-	}
-	if(s[m][i]==-2){ /* Vital status is really unknown */
-	  nbwarn++;
-	  if((int)anint[m][i] == 9999){  /*  Has the vital status really been verified? */
-	    printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
-	    fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
-	  }
-	  break;
-	}
-	break;
+				if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
+					if(firsthree == 0){
+						printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+						firsthree=1;
+					}
+					fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+					mw[++mi][i]=m;
+				}
+				if(s[m][i]==-2){ /* Vital status is really unknown */
+					nbwarn++;
+					if((int)anint[m][i] == 9999){  /*  Has the vital status really been verified? */
+						printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+						fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+					}
+					break;
+				}
+				break;
       }
       else
-	m++;
+				m++;
     }/* end while */
     
     /* After last pass */
     if (s[m][i] > nlstate){  /* In a death state */
       mi++;	/* Death is another wave */
       /* if(mi==0)  never been interviewed correctly before death */
-	 /* Only death is a correct wave */
+			/* Only death is a correct wave */
       mw[mi][i]=m;
     }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */
       /* m++; */
@@ -4154,117 +4187,117 @@ void  concatwav(int wav[], int **dh, int
       /* mw[mi][i]=m; */
       nberr++;
       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
-	if(firstwo==0){
-	  printf("Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
-	  firstwo=1;
-	}
-	fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+				if(firstwo==0){
+					printf("Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+					firstwo=1;
+				}
+				fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
       }else{ /* end date of interview is known */
-	/* death is known but not confirmed by death status at any wave */
-	if(firstfour==0){
-	  printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
-	  firstfour=1;
-	}
-	fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+				/* death is known but not confirmed by death status at any wave */
+				if(firstfour==0){
+					printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+					firstfour=1;
+				}
+				fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
       }
     }
     wav[i]=mi;
     if(mi==0){
       nbwarn++;
       if(first==0){
-	printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
-	first=1;
+				printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
+				first=1;
       }
       if(first==1){
-	fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
+				fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
       }
     } /* end mi==0 */
   } /* End individuals */
   /* wav and mw are no more changed */
-
+	
   
   for(i=1; i<=imx; i++){
     for(mi=1; mi<wav[i];mi++){
       if (stepm <=0)
-	dh[mi][i]=1;
+				dh[mi][i]=1;
       else{
-	if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
-	  if (agedc[i] < 2*AGESUP) {
-	    j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
-	    if(j==0) j=1;  /* Survives at least one month after exam */
-	    else if(j<0){
-	      nberr++;
-	      printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-	      j=1; /* Temporary Dangerous patch */
-	      printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
-	      fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-	      fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
-	    }
-	    k=k+1;
-	    if (j >= jmax){
-	      jmax=j;
-	      ijmax=i;
-	    }
-	    if (j <= jmin){
-	      jmin=j;
-	      ijmin=i;
-	    }
-	    sum=sum+j;
-	    /*if (j<0) printf("j=%d num=%d \n",j,i);*/
-	    /*	  printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
-	  }
-	}
-	else{
-	  j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
+				if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
+					if (agedc[i] < 2*AGESUP) {
+						j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
+						if(j==0) j=1;  /* Survives at least one month after exam */
+						else if(j<0){
+							nberr++;
+							printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+							j=1; /* Temporary Dangerous patch */
+							printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+							fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+							fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
+						}
+						k=k+1;
+						if (j >= jmax){
+							jmax=j;
+							ijmax=i;
+						}
+						if (j <= jmin){
+							jmin=j;
+							ijmin=i;
+						}
+						sum=sum+j;
+						/*if (j<0) printf("j=%d num=%d \n",j,i);*/
+						/*	  printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
+					}
+				}
+				else{
+					j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
 /* 	  if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
-
-	  k=k+1;
-	  if (j >= jmax) {
-	    jmax=j;
-	    ijmax=i;
-	  }
-	  else if (j <= jmin){
-	    jmin=j;
-	    ijmin=i;
-	  }
-	  /*	    if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
-	  /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
-	  if(j<0){
-	    nberr++;
-	    printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-	    fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-	  }
-	  sum=sum+j;
-	}
-	jk= j/stepm;
-	jl= j -jk*stepm;
-	ju= j -(jk+1)*stepm;
-	if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
-	  if(jl==0){
-	    dh[mi][i]=jk;
-	    bh[mi][i]=0;
-	  }else{ /* We want a negative bias in order to only have interpolation ie
-		  * to avoid the price of an extra matrix product in likelihood */
-	    dh[mi][i]=jk+1;
-	    bh[mi][i]=ju;
-	  }
-	}else{
-	  if(jl <= -ju){
-	    dh[mi][i]=jk;
-	    bh[mi][i]=jl;	/* bias is positive if real duration
-				 * is higher than the multiple of stepm and negative otherwise.
-				 */
-	  }
-	  else{
-	    dh[mi][i]=jk+1;
-	    bh[mi][i]=ju;
-	  }
-	  if(dh[mi][i]==0){
-	    dh[mi][i]=1; /* At least one step */
-	    bh[mi][i]=ju; /* At least one step */
-	    /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
-	  }
-	} /* end if mle */
+					
+					k=k+1;
+					if (j >= jmax) {
+						jmax=j;
+						ijmax=i;
+					}
+					else if (j <= jmin){
+						jmin=j;
+						ijmin=i;
+					}
+					/*	    if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
+					/*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
+					if(j<0){
+						nberr++;
+						printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+						fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+					}
+					sum=sum+j;
+				}
+				jk= j/stepm;
+				jl= j -jk*stepm;
+				ju= j -(jk+1)*stepm;
+				if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
+					if(jl==0){
+						dh[mi][i]=jk;
+						bh[mi][i]=0;
+					}else{ /* We want a negative bias in order to only have interpolation ie
+									* to avoid the price of an extra matrix product in likelihood */
+						dh[mi][i]=jk+1;
+						bh[mi][i]=ju;
+					}
+				}else{
+					if(jl <= -ju){
+						dh[mi][i]=jk;
+						bh[mi][i]=jl;	/* bias is positive if real duration
+													 * is higher than the multiple of stepm and negative otherwise.
+													 */
+					}
+					else{
+						dh[mi][i]=jk+1;
+						bh[mi][i]=ju;
+					}
+					if(dh[mi][i]==0){
+						dh[mi][i]=1; /* At least one step */
+						bh[mi][i]=ju; /* At least one step */
+						/*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
+					}
+				} /* end if mle */
       }
     } /* end wave */
   }
@@ -4520,7 +4553,7 @@ void cvevsij(double ***eij, double x[],
 
 {
   /* Covariances of health expectancies eij and of total life expectancies according
-   to initial status i, ei. .
+     to initial status i, ei. .
   */
   int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
   int nhstepma, nstepma; /* Decreasing with age */
@@ -4626,47 +4659,47 @@ void cvevsij(double ***eij, double x[],
        decrease memory allocation */
     for(theta=1; theta <=npar; theta++){
       for(i=1; i<=npar; i++){ 
-				xp[i] = x[i] + (i==theta ?delti[theta]:0);
-				xm[i] = x[i] - (i==theta ?delti[theta]:0);
+	xp[i] = x[i] + (i==theta ?delti[theta]:0);
+	xm[i] = x[i] - (i==theta ?delti[theta]:0);
       }
       hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);  
       hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);  
 			
       for(j=1; j<= nlstate; j++){
-				for(i=1; i<=nlstate; i++){
-					for(h=0; h<=nhstepm-1; h++){
-						gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
-						gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
-					}
-				}
+	for(i=1; i<=nlstate; i++){
+	  for(h=0; h<=nhstepm-1; h++){
+	    gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
+	    gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
+	  }
+	}
       }
 			
       for(ij=1; ij<= nlstate*nlstate; ij++)
-				for(h=0; h<=nhstepm-1; h++){
-					gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
-				}
+	for(h=0; h<=nhstepm-1; h++){
+	  gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
+	}
     }/* End theta */
     
     
     for(h=0; h<=nhstepm-1; h++)
       for(j=1; j<=nlstate*nlstate;j++)
-				for(theta=1; theta <=npar; theta++)
-					trgradg[h][j][theta]=gradg[h][theta][j];
+	for(theta=1; theta <=npar; theta++)
+	  trgradg[h][j][theta]=gradg[h][theta][j];
     
 		
-		for(ij=1;ij<=nlstate*nlstate;ij++)
+    for(ij=1;ij<=nlstate*nlstate;ij++)
       for(ji=1;ji<=nlstate*nlstate;ji++)
-				varhe[ij][ji][(int)age] =0.;
+	varhe[ij][ji][(int)age] =0.;
 		
-		printf("%d|",(int)age);fflush(stdout);
-		fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
-		for(h=0;h<=nhstepm-1;h++){
+    printf("%d|",(int)age);fflush(stdout);
+    fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
+    for(h=0;h<=nhstepm-1;h++){
       for(k=0;k<=nhstepm-1;k++){
-				matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
-				matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
-				for(ij=1;ij<=nlstate*nlstate;ij++)
-					for(ji=1;ji<=nlstate*nlstate;ji++)
-						varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
+	matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
+	matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
+	for(ij=1;ij<=nlstate*nlstate;ij++)
+	  for(ji=1;ji<=nlstate*nlstate;ji++)
+	    varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }
     }
 		
@@ -4674,22 +4707,22 @@ void cvevsij(double ***eij, double x[],
     hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)
-				for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
-					eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
+	for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
+	  eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
 					
-					/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
+	  /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
 					
-				}
+	}
 		
     fprintf(ficresstdeij,"%3.0f",age );
     for(i=1; i<=nlstate;i++){
       eip=0.;
       vip=0.;
       for(j=1; j<=nlstate;j++){
-				eip += eij[i][j][(int)age];
-				for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
-					vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
-				fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
+	eip += eij[i][j][(int)age];
+	for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
+	  vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
+	fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
       }
       fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
     }
@@ -4698,13 +4731,13 @@ void cvevsij(double ***eij, double x[],
     fprintf(ficrescveij,"%3.0f",age );
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){
-				cptj= (j-1)*nlstate+i;
-				for(i2=1; i2<=nlstate;i2++)
-					for(j2=1; j2<=nlstate;j2++){
-						cptj2= (j2-1)*nlstate+i2;
-						if(cptj2 <= cptj)
-							fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
-					}
+	cptj= (j-1)*nlstate+i;
+	for(i2=1; i2<=nlstate;i2++)
+	  for(j2=1; j2<=nlstate;j2++){
+	    cptj2= (j2-1)*nlstate+i2;
+	    if(cptj2 <= cptj)
+	      fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
+	  }
       }
     fprintf(ficrescveij,"\n");
 		
@@ -5161,86 +5194,86 @@ void cvevsij(double ***eij, double x[],
 
 /************ Variance of one-step probabilities  ******************/
 void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[])
-{
-  int i, j=0,  k1, l1, tj;
-  int k2, l2, j1,  z1;
-  int k=0, l;
-  int first=1, first1, first2;
-  double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
-  double **dnewm,**doldm;
-  double *xp;
-  double *gp, *gm;
-  double **gradg, **trgradg;
-  double **mu;
-  double age, cov[NCOVMAX+1];
-  double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
-  int theta;
-  char fileresprob[FILENAMELENGTH];
-  char fileresprobcov[FILENAMELENGTH];
-  char fileresprobcor[FILENAMELENGTH];
-  double ***varpij;
-
-  strcpy(fileresprob,"PROB_"); 
-  strcat(fileresprob,fileres);
-  if((ficresprob=fopen(fileresprob,"w"))==NULL) {
-    printf("Problem with resultfile: %s\n", fileresprob);
-    fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
-  }
-  strcpy(fileresprobcov,"PROBCOV_"); 
-  strcat(fileresprobcov,fileresu);
-  if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
-    printf("Problem with resultfile: %s\n", fileresprobcov);
-    fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
-  }
-  strcpy(fileresprobcor,"PROBCOR_"); 
-  strcat(fileresprobcor,fileresu);
-  if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
-    printf("Problem with resultfile: %s\n", fileresprobcor);
-    fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
-  }
-  printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
-  fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
-  printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
-  fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
-  printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
-  fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
-  pstamp(ficresprob);
-  fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
-  fprintf(ficresprob,"# Age");
-  pstamp(ficresprobcov);
-  fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
-  fprintf(ficresprobcov,"# Age");
-  pstamp(ficresprobcor);
-  fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
-  fprintf(ficresprobcor,"# Age");
-
+ {
+   int i, j=0,  k1, l1, tj;
+   int k2, l2, j1,  z1;
+   int k=0, l;
+   int first=1, first1, first2;
+   double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
+   double **dnewm,**doldm;
+   double *xp;
+   double *gp, *gm;
+   double **gradg, **trgradg;
+   double **mu;
+   double age, cov[NCOVMAX+1];
+   double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
+   int theta;
+   char fileresprob[FILENAMELENGTH];
+   char fileresprobcov[FILENAMELENGTH];
+   char fileresprobcor[FILENAMELENGTH];
+   double ***varpij;
+
+   strcpy(fileresprob,"PROB_"); 
+   strcat(fileresprob,fileres);
+   if((ficresprob=fopen(fileresprob,"w"))==NULL) {
+     printf("Problem with resultfile: %s\n", fileresprob);
+     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
+   }
+   strcpy(fileresprobcov,"PROBCOV_"); 
+   strcat(fileresprobcov,fileresu);
+   if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
+     printf("Problem with resultfile: %s\n", fileresprobcov);
+     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
+   }
+   strcpy(fileresprobcor,"PROBCOR_"); 
+   strcat(fileresprobcor,fileresu);
+   if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
+     printf("Problem with resultfile: %s\n", fileresprobcor);
+     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
+   }
+   printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+   fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+   printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+   fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+   printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
+   fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
+   pstamp(ficresprob);
+   fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
+   fprintf(ficresprob,"# Age");
+   pstamp(ficresprobcov);
+   fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
+   fprintf(ficresprobcov,"# Age");
+   pstamp(ficresprobcor);
+   fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
+   fprintf(ficresprobcor,"# Age");
 
-  for(i=1; i<=nlstate;i++)
-    for(j=1; j<=(nlstate+ndeath);j++){
-      fprintf(ficresprob," p%1d-%1d (SE)",i,j);
-      fprintf(ficresprobcov," p%1d-%1d ",i,j);
-      fprintf(ficresprobcor," p%1d-%1d ",i,j);
-    }  
- /* fprintf(ficresprob,"\n");
-  fprintf(ficresprobcov,"\n");
-  fprintf(ficresprobcor,"\n");
- */
-  xp=vector(1,npar);
-  dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
-  doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
-  mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
-  varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage);
-  first=1;
-  fprintf(ficgp,"\n# Routine varprob");
-  fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
-  fprintf(fichtm,"\n");
 
-  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
-  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
+   for(i=1; i<=nlstate;i++)
+     for(j=1; j<=(nlstate+ndeath);j++){
+       fprintf(ficresprob," p%1d-%1d (SE)",i,j);
+       fprintf(ficresprobcov," p%1d-%1d ",i,j);
+       fprintf(ficresprobcor," p%1d-%1d ",i,j);
+     }  
+   /* fprintf(ficresprob,"\n");
+      fprintf(ficresprobcov,"\n");
+      fprintf(ficresprobcor,"\n");
+   */
+   xp=vector(1,npar);
+   dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+   doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+   mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
+   varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage);
+   first=1;
+   fprintf(ficgp,"\n# Routine varprob");
+   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
+   fprintf(fichtm,"\n");
+
+   fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
+   fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+   fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
  They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
-  fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
+   fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
 It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \
 would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \
 standard deviations wide on each axis. <br>\
@@ -5248,252 +5281,252 @@ standard deviations wide on each axis. <
  and made the appropriate rotation to look at the uncorrelated principal directions.<br>\
 To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n");
 
-  cov[1]=1;
-  /* tj=cptcoveff; */
-  tj = (int) pow(2,cptcoveff);
-  if (cptcovn<1) {tj=1;ncodemax[1]=1;}
-  j1=0;
-  for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates */
-		if  (cptcovn>0) {
-			fprintf(ficresprob, "\n#********** Variable "); 
-			for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-			fprintf(ficresprob, "**********\n#\n");
-			fprintf(ficresprobcov, "\n#********** Variable "); 
-			for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-			fprintf(ficresprobcov, "**********\n#\n");
+   cov[1]=1;
+   /* tj=cptcoveff; */
+   tj = (int) pow(2,cptcoveff);
+   if (cptcovn<1) {tj=1;ncodemax[1]=1;}
+   j1=0;
+   for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates */
+     if  (cptcovn>0) {
+       fprintf(ficresprob, "\n#********** Variable "); 
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+       fprintf(ficresprob, "**********\n#\n");
+       fprintf(ficresprobcov, "\n#********** Variable "); 
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+       fprintf(ficresprobcov, "**********\n#\n");
 			
-			fprintf(ficgp, "\n#********** Variable "); 
-			for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-			fprintf(ficgp, "**********\n#\n");
+       fprintf(ficgp, "\n#********** Variable "); 
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+       fprintf(ficgp, "**********\n#\n");
 			
 			
-			fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
-			for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-			fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
+       fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+       fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
 			
-			fprintf(ficresprobcor, "\n#********** Variable ");    
-			for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-			fprintf(ficresprobcor, "**********\n#");    
-			if(invalidvarcomb[j1]){
-			  fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); 
-			  fprintf(fichtmcov,"\n<h3>Combination (%d) ignored because no cases </h3>\n",j1); 
-			  continue;
-			}
-		}
-		gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
-		trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
-		gp=vector(1,(nlstate)*(nlstate+ndeath));
-		gm=vector(1,(nlstate)*(nlstate+ndeath));
-		for (age=bage; age<=fage; age ++){ 
-			cov[2]=age;
-			if(nagesqr==1)
-				cov[3]= age*age;
-			for (k=1; k<=cptcovn;k++) {
-				cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
-				/*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
-																																	 * 1  1 1 1 1
-																																	 * 2  2 1 1 1
-																																	 * 3  1 2 1 1
-																																	 */
-				/* nbcode[1][1]=0 nbcode[1][2]=1;*/
-			}
-			/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-			for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
-			for (k=1; k<=cptcovprod;k++)
-				cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+       fprintf(ficresprobcor, "\n#********** Variable ");    
+       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+       fprintf(ficresprobcor, "**********\n#");    
+       if(invalidvarcomb[j1]){
+	 fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); 
+	 fprintf(fichtmcov,"\n<h3>Combination (%d) ignored because no cases </h3>\n",j1); 
+	 continue;
+       }
+     }
+     gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
+     trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+     gp=vector(1,(nlstate)*(nlstate+ndeath));
+     gm=vector(1,(nlstate)*(nlstate+ndeath));
+     for (age=bage; age<=fage; age ++){ 
+       cov[2]=age;
+       if(nagesqr==1)
+	 cov[3]= age*age;
+       for (k=1; k<=cptcovn;k++) {
+	 cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
+	 /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
+								    * 1  1 1 1 1
+								    * 2  2 1 1 1
+								    * 3  1 2 1 1
+								    */
+	 /* nbcode[1][1]=0 nbcode[1][2]=1;*/
+       }
+       /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+       for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+       for (k=1; k<=cptcovprod;k++)
+	 cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
 			
 			
-			for(theta=1; theta <=npar; theta++){
-				for(i=1; i<=npar; i++)
-					xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);
+       for(theta=1; theta <=npar; theta++){
+	 for(i=1; i<=npar; i++)
+	   xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);
 				
-				pmij(pmmij,cov,ncovmodel,xp,nlstate);
+	 pmij(pmmij,cov,ncovmodel,xp,nlstate);
 				
-				k=0;
-				for(i=1; i<= (nlstate); i++){
-					for(j=1; j<=(nlstate+ndeath);j++){
-						k=k+1;
-						gp[k]=pmmij[i][j];
-					}
-				}
+	 k=0;
+	 for(i=1; i<= (nlstate); i++){
+	   for(j=1; j<=(nlstate+ndeath);j++){
+	     k=k+1;
+	     gp[k]=pmmij[i][j];
+	   }
+	 }
 				
-				for(i=1; i<=npar; i++)
-					xp[i] = x[i] - (i==theta ?delti[theta]:(double)0);
+	 for(i=1; i<=npar; i++)
+	   xp[i] = x[i] - (i==theta ?delti[theta]:(double)0);
 				
-				pmij(pmmij,cov,ncovmodel,xp,nlstate);
-				k=0;
-				for(i=1; i<=(nlstate); i++){
-					for(j=1; j<=(nlstate+ndeath);j++){
-						k=k+1;
-						gm[k]=pmmij[i][j];
-					}
-				}
+	 pmij(pmmij,cov,ncovmodel,xp,nlstate);
+	 k=0;
+	 for(i=1; i<=(nlstate); i++){
+	   for(j=1; j<=(nlstate+ndeath);j++){
+	     k=k+1;
+	     gm[k]=pmmij[i][j];
+	   }
+	 }
 				
-				for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) 
-					gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta];  
-			}
+	 for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) 
+	   gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta];  
+       }
 
-			for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
-				for(theta=1; theta <=npar; theta++)
-					trgradg[j][theta]=gradg[theta][j];
+       for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
+	 for(theta=1; theta <=npar; theta++)
+	   trgradg[j][theta]=gradg[theta][j];
 			
-			matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); 
-			matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
+       matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); 
+       matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
 			
-			pmij(pmmij,cov,ncovmodel,x,nlstate);
+       pmij(pmmij,cov,ncovmodel,x,nlstate);
 			
-			k=0;
-			for(i=1; i<=(nlstate); i++){
-				for(j=1; j<=(nlstate+ndeath);j++){
-					k=k+1;
-					mu[k][(int) age]=pmmij[i][j];
-				}
-			}
-     	for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
-				for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
-					varpij[i][j][(int)age] = doldm[i][j];
+       k=0;
+       for(i=1; i<=(nlstate); i++){
+	 for(j=1; j<=(nlstate+ndeath);j++){
+	   k=k+1;
+	   mu[k][(int) age]=pmmij[i][j];
+	 }
+       }
+       for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
+	 for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
+	   varpij[i][j][(int)age] = doldm[i][j];
 			
-			/*printf("\n%d ",(int)age);
-				for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
-				printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
-				fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
-				}*/
+       /*printf("\n%d ",(int)age);
+	 for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+	 printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+	 fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+	 }*/
 			
-			fprintf(ficresprob,"\n%d ",(int)age);
-			fprintf(ficresprobcov,"\n%d ",(int)age);
-			fprintf(ficresprobcor,"\n%d ",(int)age);
+       fprintf(ficresprob,"\n%d ",(int)age);
+       fprintf(ficresprobcov,"\n%d ",(int)age);
+       fprintf(ficresprobcor,"\n%d ",(int)age);
 			
-			for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
-				fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
-			for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
-				fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
-				fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
-			}
-			i=0;
-			for (k=1; k<=(nlstate);k++){
-				for (l=1; l<=(nlstate+ndeath);l++){ 
-					i++;
-					fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
-					fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
-					for (j=1; j<=i;j++){
-						/* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */
-						fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
-						fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
-					}
-				}
-			}/* end of loop for state */
-		} /* end of loop for age */
-		free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
-		free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
-		free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
-		free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
-    
-		/* Confidence intervalle of pij  */
-		/*
-			fprintf(ficgp,"\nunset parametric;unset label");
-			fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
-			fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
-			fprintf(fichtm,"\n<br>Probability with  confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
-			fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
-			fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
-			fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
-		*/
-		
-		/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
-		first1=1;first2=2;
-		for (k2=1; k2<=(nlstate);k2++){
-			for (l2=1; l2<=(nlstate+ndeath);l2++){ 
-				if(l2==k2) continue;
-				j=(k2-1)*(nlstate+ndeath)+l2;
-				for (k1=1; k1<=(nlstate);k1++){
-					for (l1=1; l1<=(nlstate+ndeath);l1++){ 
-						if(l1==k1) continue;
-						i=(k1-1)*(nlstate+ndeath)+l1;
-						if(i<=j) continue;
-						for (age=bage; age<=fage; age ++){ 
-							if ((int)age %5==0){
-								v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
-								v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM;
-								cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
-								mu1=mu[i][(int) age]/stepm*YEARM ;
-								mu2=mu[j][(int) age]/stepm*YEARM;
-								c12=cv12/sqrt(v1*v2);
-								/* Computing eigen value of matrix of covariance */
-								lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
-								lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
-								if ((lc2 <0) || (lc1 <0) ){
-									if(first2==1){
-										first1=0;
-										printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
-									}
-									fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog);
-									/* lc1=fabs(lc1); */ /* If we want to have them positive */
-									/* lc2=fabs(lc2); */
-								}
+       for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
+	 fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
+       for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+	 fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
+	 fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
+       }
+       i=0;
+       for (k=1; k<=(nlstate);k++){
+	 for (l=1; l<=(nlstate+ndeath);l++){ 
+	   i++;
+	   fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
+	   fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
+	   for (j=1; j<=i;j++){
+	     /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */
+	     fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
+	     fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
+	   }
+	 }
+       }/* end of loop for state */
+     } /* end of loop for age */
+     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+     free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+     free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+     free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+    
+     /* Confidence intervalle of pij  */
+     /*
+       fprintf(ficgp,"\nunset parametric;unset label");
+       fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
+       fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+       fprintf(fichtm,"\n<br>Probability with  confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
+       fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
+       fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
+       fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
+     */
+		
+     /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
+     first1=1;first2=2;
+     for (k2=1; k2<=(nlstate);k2++){
+       for (l2=1; l2<=(nlstate+ndeath);l2++){ 
+	 if(l2==k2) continue;
+	 j=(k2-1)*(nlstate+ndeath)+l2;
+	 for (k1=1; k1<=(nlstate);k1++){
+	   for (l1=1; l1<=(nlstate+ndeath);l1++){ 
+	     if(l1==k1) continue;
+	     i=(k1-1)*(nlstate+ndeath)+l1;
+	     if(i<=j) continue;
+	     for (age=bage; age<=fage; age ++){ 
+	       if ((int)age %5==0){
+		 v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
+		 v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM;
+		 cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
+		 mu1=mu[i][(int) age]/stepm*YEARM ;
+		 mu2=mu[j][(int) age]/stepm*YEARM;
+		 c12=cv12/sqrt(v1*v2);
+		 /* Computing eigen value of matrix of covariance */
+		 lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+		 lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+		 if ((lc2 <0) || (lc1 <0) ){
+		   if(first2==1){
+		     first1=0;
+		     printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
+		   }
+		   fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog);
+		   /* lc1=fabs(lc1); */ /* If we want to have them positive */
+		   /* lc2=fabs(lc2); */
+		 }
 								
-								/* Eigen vectors */
-								v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
-								/*v21=sqrt(1.-v11*v11); *//* error */
-								v21=(lc1-v1)/cv12*v11;
-								v12=-v21;
-								v22=v11;
-								tnalp=v21/v11;
-								if(first1==1){
-									first1=0;
-									printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
-								}
-								fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
-								/*printf(fignu*/
-								/* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
-								/* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
-								if(first==1){
-									first=0;
-									fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
-									fprintf(ficgp,"\nset parametric;unset label");
-									fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
-									fprintf(ficgp,"\nset ter svg size 640, 480");
-									fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
+		 /* Eigen vectors */
+		 v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
+		 /*v21=sqrt(1.-v11*v11); *//* error */
+		 v21=(lc1-v1)/cv12*v11;
+		 v12=-v21;
+		 v22=v11;
+		 tnalp=v21/v11;
+		 if(first1==1){
+		   first1=0;
+		   printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
+		 }
+		 fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);
+		 /*printf(fignu*/
+		 /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
+		 /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
+		 if(first==1){
+		   first=0;
+		   fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
+		   fprintf(ficgp,"\nset parametric;unset label");
+		   fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
+		   fprintf(ficgp,"\nset ter svg size 640, 480");
+		   fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
  :<a href=\"%s_%d%1d%1d-%1d%1d.svg\">																		\
 %s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
-													subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,	\
-													subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
-									fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
-									fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
-									fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
-		    fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
-		    fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
-		    fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",	\
-								mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),										\
-								mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
-								}else{
-									first=0;
-									fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
-									fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
-									fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
-									fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \
-													mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),					\
-													mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
-								}/* if first */
-							} /* age mod 5 */
-						} /* end loop age */
-						fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
-						first=1;
-					} /*l12 */
-				} /* k12 */
-			} /*l1 */
-		}/* k1 */
-	}  /* loop on combination of covariates j1 */
-  free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
-  free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
-  free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
-  free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
-  free_vector(xp,1,npar);
-  fclose(ficresprob);
-  fclose(ficresprobcov);
-  fclose(ficresprobcor);
-  fflush(ficgp);
-  fflush(fichtmcov);
-}
+			   subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,	\
+			   subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+		   fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+		   fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
+		   fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+		   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+		   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+		   fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",	\
+			   mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),										\
+			   mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+		 }else{
+		   first=0;
+		   fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
+		   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
+		   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
+		   fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \
+			   mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),					\
+			   mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+		 }/* if first */
+	       } /* age mod 5 */
+	     } /* end loop age */
+	     fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+	     first=1;
+	   } /*l12 */
+	 } /* k12 */
+       } /*l1 */
+     }/* k1 */
+   }  /* loop on combination of covariates j1 */
+   free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
+   free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
+   free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+   free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
+   free_vector(xp,1,npar);
+   fclose(ficresprob);
+   fclose(ficresprobcov);
+   fclose(ficresprobcor);
+   fflush(ficgp);
+   fflush(fichtmcov);
+ }
 
 
 /******************* Printing html file ***********/
@@ -5536,81 +5569,81 @@ void printinghtml(char fileresu[], char
    <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
    }
 
-fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
+   fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
- m=pow(2,cptcoveff);
- if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+   m=pow(2,cptcoveff);
+   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
- jj1=0;
- for(k1=1; k1<=m;k1++){
+   jj1=0;
+   for(k1=1; k1<=m;k1++){
 
-   /* 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," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
-		 if(invalidvarcomb[k1]){
-			 fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); 
-			 printf("\nCombination (%d) ignored because no cases \n",k1); 
-			 continue;
-		 }
-	 }
-	 /* 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> \
+     /* 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," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
+       if(invalidvarcomb[k1]){
+	 fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); 
+	 printf("\nCombination (%d) ignored because no cases \n",k1); 
+	 continue;
+       }
+     }
+     /* 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);
-	 /* 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> \
+     /* 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);     
-	 /* 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\
+     /* 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); 
-	 /* 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> \
+     /* 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);
-	 }
-	 /* 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.\
+     }
+     /* 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);
-	 }
-	 /* 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);
-	 }
-	 if(backcast==1){
-     /* Period (stable) back prevalence in each health state */
+     }
+     /* Period (stable) 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> \
+       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);
+     }
+     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);
+       }
      }
-	 }
-	 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> \
+     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);
-		 }
-	 }
+       }
+     }
 	 
-	 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> \
+     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);
-	 }
-   /* } /\* end i1 *\/ */
- }/* End k1 */
- fprintf(fichtm,"</ul>");
+     }
+     /* } /\* end i1 *\/ */
+   }/* End k1 */
+   fprintf(fichtm,"</ul>");
 
- fprintf(fichtm,"\
+   fprintf(fichtm,"\
 \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
  - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).<br> \
@@ -5622,32 +5655,32 @@ variances but at the covariance matrix.
 covariance matrix of the one-step probabilities. \
 See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
 
- fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
- fprintf(fichtm,"\
+   fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
+	   subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
+   fprintf(fichtm,"\
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
+	   subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
 
- fprintf(fichtm,"\
+   fprintf(fichtm,"\
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_"));
- fprintf(fichtm,"\
+	   subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_"));
+   fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",
 	   estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_"));
- fprintf(fichtm,"\
+   fprintf(fichtm,"\
  - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",
 	   estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_"));
- fprintf(fichtm,"\
+   fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
-	 estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
- fprintf(fichtm,"\
+	   estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
+   fprintf(fichtm,"\
  - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n",
-	 estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
- fprintf(fichtm,"\
+	   estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
+   fprintf(fichtm,"\
  - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
-	 subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
+	   subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
 
 /*  if(popforecast==1) fprintf(fichtm,"\n */
 /*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
@@ -5655,26 +5688,26 @@ See page 'Matrix of variance-covariance
 /* 	<br>",fileres,fileres,fileres,fileres); */
 /*  else  */
 /*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
- fflush(fichtm);
- fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
+   fflush(fichtm);
+   fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
 
- m=pow(2,cptcoveff);
- if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+   m=pow(2,cptcoveff);
+   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
- jj1=0;
- for(k1=1; k1<=m;k1++){
-  /* for(i1=1; i1<=ncodemax[k1];i1++){ */
-		 jj1++;
+   jj1=0;
+   for(k1=1; k1<=m;k1++){
+     /* 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)]);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
 
-			 if(invalidvarcomb[k1]){
-				 fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); 
-				 continue;
-			 }
+       if(invalidvarcomb[k1]){
+	 fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); 
+	 continue;
+       }
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
@@ -5687,22 +5720,22 @@ true period expectancies (those weighted
  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);
-   /* } /\* end i1 *\/ */
- }/* End k1 */
- fprintf(fichtm,"</ul>");
- fflush(fichtm);
+     /* } /\* end i1 *\/ */
+   }/* End k1 */
+   fprintf(fichtm,"</ul>");
+   fflush(fichtm);
 }
 
 /******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
-	char gplotcondition[132];
+  char gplotcondition[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
   int lv=0, vlv=0, kl=0;
   int ng=0;
   int vpopbased;
-	int ioffset; /* variable offset for columns */
+  int ioffset; /* variable offset for columns */
 
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
@@ -5711,158 +5744,162 @@ true period expectancies (those weighted
 
   /*#ifdef windows */
   fprintf(ficgp,"cd \"%s\" \n",pathc);
-    /*#endif */
+  /*#endif */
   m=pow(2,cptcoveff);
 
   /* Contribution to likelihood */
   /* Plot the probability implied in the likelihood */
-    fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
-    fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
-    /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
-    fprintf(ficgp,"\nset ter pngcairo size 640, 480");
+  fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
+  fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
+  /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
+  fprintf(ficgp,"\nset ter pngcairo size 640, 480");
 /* nice for mle=4 plot by number of matrix products.
    replot  "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
 /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */
-    /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
-    fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
-    fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
-    fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
-    fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
-    for (i=1; i<= nlstate ; i ++) {
-      fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
-      fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
-      fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
-      for (j=2; j<= nlstate+ndeath ; j ++) {
-				fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
-      }
-      fprintf(ficgp,";\nset out; unset ylabel;\n"); 
-    }
-    /* unset log; plot  "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u  2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */		 
-    /* fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
-    /* fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
-    fprintf(ficgp,"\nset out;unset log\n");
-    /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+  /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
+  fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
+  fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
+  fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
+  fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
+  for (i=1; i<= nlstate ; i ++) {
+    fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
+    fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
+    fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
+    for (j=2; j<= nlstate+ndeath ; j ++) {
+      fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+    }
+    fprintf(ficgp,";\nset out; unset ylabel;\n"); 
+  }
+  /* unset log; plot  "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u  2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */		 
+  /* fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
+  /* fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
+  fprintf(ficgp,"\nset out;unset log\n");
+  /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
 
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
- /* 1eme*/
+  /* 1eme*/
   for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
     for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */
       /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
       fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
-				lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
-				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-				vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
-			/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
-				fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
+	/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
       fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-						fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-						continue;
-			}
+      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\
+      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 */
-				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' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
-											6+(cpt-1),  cpt );
-					}else{
-						fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
-						kl++;
-					}
-				} /* end covariate */
-			}
-			fprintf(ficgp,"\nset out \n");
+      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' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+		      6+(cpt-1),  cpt );
+	    }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 ++) { 
 
-      fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
-      for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-				lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-				vlv= nbcode[Tvaraff[k]][lv];
-				fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-      }
-      fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-						fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-						continue;
-			}
+    fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+    for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+      lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+      /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+      /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+      /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+      vlv= nbcode[Tvaraff[k]][lv];
+      fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+    }
+    fprintf(ficgp,"\n#\n");
+    if(invalidvarcomb[k1]){
+      fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+      continue;
+    }
 			
-			fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
-			for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-				if(vpopbased==0)
-					fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-				else
-					fprintf(ficgp,"\nreplot ");
-				for (i=1; i<= nlstate+1 ; i ++) {
-					k=2*i;
-					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
-					for (j=1; j<= nlstate+1 ; j ++) {
-						if (j==i) fprintf(ficgp," %%lf (%%lf)");
-						else fprintf(ficgp," %%*lf (%%*lf)");
-					}   
-					if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-					else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
-					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-					for (j=1; j<= nlstate+1 ; j ++) {
-						if (j==i) fprintf(ficgp," %%lf (%%lf)");
-						else fprintf(ficgp," %%*lf (%%*lf)");
-					}   
-					fprintf(ficgp,"\" t\"\" w l lt 0,");
-					fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-					for (j=1; j<= nlstate+1 ; j ++) {
-						if (j==i) fprintf(ficgp," %%lf (%%lf)");
-						else fprintf(ficgp," %%*lf (%%*lf)");
-					}   
-					if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-					else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
-				} /* state */
-			} /* vpopbased */
-			fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+    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 */
 	
 	
@@ -5872,18 +5909,18 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ 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);
+	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	vlv= nbcode[Tvaraff[k]][lv];
+	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
       fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-						fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-						continue;
-			}
+      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);
@@ -5891,41 +5928,41 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
       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);
+	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);*/
+	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);
     }
   }
   
-	/* 4eme */
+  /* 4eme */
   /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
 
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-				lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-				/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-				/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-				/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-				vlv= nbcode[Tvaraff[k]][lv];
-				fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	vlv= nbcode[Tvaraff[k]][lv];
+	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
       fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-							fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-							continue;
-			}
+      if(invalidvarcomb[k1]){
+	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+	continue;
+      }
 			
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
@@ -5934,16 +5971,16 @@ 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,", '' ");
-				}
-				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);
+	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*/ 
@@ -5953,7 +5990,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-
+			
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 				lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
@@ -5964,10 +6001,10 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 				fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
       fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-						fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-						continue;
-			}
+      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\
@@ -6003,7 +6040,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* 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 (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);
       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 */
@@ -6014,15 +6051,15 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 				fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
       }
       fprintf(ficgp,"\n#\n");
-			if(invalidvarcomb[k1]){
-						fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-						continue;
-			}
-
+      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,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n\
-unset log y\n\
+set ter svg size 640, 480\n																						 \
+unset log y\n																													 \
 plot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){
@@ -6039,8 +6076,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
       fprintf(ficgp,"\nset out\n");
     } /* end cpt state*/ 
   } /* end covariate */  
-
-
+	
+	
 /* 7eme */
   if(backcast == 1){
     /* CV back preval stable (period) for each covariate */
@@ -6051,14 +6088,14 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 					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 */
+	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 					vlv= nbcode[Tvaraff[k]][lv];
 					fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 				}
 				fprintf(ficgp,"\n#\n");
 				if(invalidvarcomb[k1]){
-								fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-								continue;
+					fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+					continue;
 				}
 				
 				fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
@@ -6087,7 +6124,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
     } /* end covariate */  
   } /* End if backcast */
   
-	/* 8eme */
+  /* 8eme */
   if(prevfcast==1){
     /* Projection from cross-sectional to stable (period) for each covariate */
     
@@ -6104,15 +6141,15 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 				}
 				fprintf(ficgp,"\n#\n");
 				if(invalidvarcomb[k1]){
-								fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
-								continue;
+					fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
+					continue;
 				}
 				
 				fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
 				fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
 				fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n	\
-unset log y\n	\
+set ter svg size 640, 480\n																							\
+unset log y\n																														\
 plot [%.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*/
@@ -6142,8 +6179,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 							ioffset=4; /* Age is in 4 */
 						}else{
 							ioffset=6; /* Age is in 6 */
-						/*#  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 */
+							/*#  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 */
 						}   
 						fprintf(ficgp," u %d:(",ioffset); 
 						kl=0;
@@ -6169,34 +6206,34 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 							fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \
 											ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
 						}else{
-								fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
-												ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+							fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
+											ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
 						}
 					} /* end if covariate */
 				} /* nlstate */
 				fprintf(ficgp,"\nset out\n");
-			} /* end cpt state*/
-		} /* end covariate */
-	} /* End if prevfcast */
+      } /* end cpt state*/
+    } /* end covariate */
+  } /* End if prevfcast */
 	
 	
-	/* proba elementaires */
-	fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+  /* proba elementaires */
+  fprintf(ficgp,"\n##############\n#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++){
       if (k != i) {
-	fprintf(ficgp,"#   current state %d\n",k);
-	for(j=1; j <=ncovmodel; j++){
-	  fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
-	  jk++; 
-	}
-	fprintf(ficgp,"\n");
+				fprintf(ficgp,"#   current state %d\n",k);
+				for(j=1; j <=ncovmodel; j++){
+					fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
+					jk++; 
+				}
+				fprintf(ficgp,"\n");
       }
     }
-   }
+  }
   fprintf(ficgp,"##############\n#\n");
-
+	
   /*goto avoid;*/
   fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
@@ -6212,264 +6249,269 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
   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,"\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 */
-	 fprintf(ficgp,"\nunset log y");
-       }else if (ng==2){
-	 fprintf(ficgp,"\nset ylabel \"Probability\"\n");
-	 fprintf(ficgp,"\nset log y");
-       }else if (ng==3){
-	 fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
-	 fprintf(ficgp,"\nset log y");
-       }else
-	 fprintf(ficgp,"\nunset title ");
-       fprintf(ficgp,"\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
-       i=1;
-       for(k2=1; k2<=nlstate; k2++) {
-	 k3=i;
-	 for(k=1; k<=(nlstate+ndeath); k++) {
-	   if (k != k2){
-	     switch( ng) {
-	     case 1:
-	       if(nagesqr==0)
-		 fprintf(ficgp," p%d+p%d*x",i,i+1);
-	       else /* nagesqr =1 */
-		 fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
-	       break;
-	     case 2: /* ng=2 */
-	       if(nagesqr==0)
-		 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-	       else /* nagesqr =1 */
-		   fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
-	       break;
-	     case 3:
-	       if(nagesqr==0)
-		 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
-	       else /* nagesqr =1 */
-		 fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
-	       break;
-	     }
-	     ij=1;/* To be checked else nbcode[0][0] wrong */
-	     for(j=3; j <=ncovmodel-nagesqr; j++) {
-	       /* 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])]); */
-		   ij++;
-		 }
-	       }
-	       else
-		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
-	     }
-	   }else{
-	     i=i-ncovmodel;
-	     if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
-	       fprintf(ficgp," (1.");
-	   }
+  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,"\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 */
+	fprintf(ficgp,"\nunset log y");
+      }else if (ng==2){
+	fprintf(ficgp,"\nset ylabel \"Probability\"\n");
+	fprintf(ficgp,"\nset log y");
+      }else if (ng==3){
+	fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
+	fprintf(ficgp,"\nset log y");
+      }else
+	fprintf(ficgp,"\nunset title ");
+      fprintf(ficgp,"\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+      i=1;
+      for(k2=1; k2<=nlstate; k2++) {
+	k3=i;
+	for(k=1; k<=(nlstate+ndeath); k++) {
+	  if (k != k2){
+	    switch( ng) {
+	    case 1:
+	      if(nagesqr==0)
+		fprintf(ficgp," p%d+p%d*x",i,i+1);
+	      else /* nagesqr =1 */
+		fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	      break;
+	    case 2: /* ng=2 */
+	      if(nagesqr==0)
+		fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+	      else /* nagesqr =1 */
+		fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	      break;
+	    case 3:
+	      if(nagesqr==0)
+		fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+	      else /* nagesqr =1 */
+		fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
+	      break;
+	    }
+	    ij=1;/* To be checked else nbcode[0][0] wrong */
+	    for(j=3; j <=ncovmodel-nagesqr; j++) {
+	      /* 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])]); */
+		  ij++;
+		}
+	      }
+	      else
+		fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+	    }
+	  }else{
+	    i=i-ncovmodel;
+	    if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
+	      fprintf(ficgp," (1.");
+	  }
 	   
-	   if(ng != 1){
-	     fprintf(ficgp,")/(1");
+	  if(ng != 1){
+	    fprintf(ficgp,")/(1");
 	     
-	     for(k1=1; k1 <=nlstate; k1++){ 
-	       if(nagesqr==0)
-		 fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
-	       else /* nagesqr =1 */
-		 fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+	    for(k1=1; k1 <=nlstate; k1++){ 
+	      if(nagesqr==0)
+		fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+	      else /* nagesqr =1 */
+		fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
 	       
-	       ij=1;
-	       for(j=3; j <=ncovmodel-nagesqr; j++){
-		 if(ij <=cptcovage) { /* Bug valgrind */
-		   if((j-2)==Tage[ij]) { /* 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++;
-		   }
-		 }
-		 else
-		   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
-	       }
-	       fprintf(ficgp,")");
-	     }
-	     fprintf(ficgp,")");
-	     if(ng ==2)
-	       fprintf(ficgp," t \"p%d%d\" ", k2,k);
-	     else /* ng= 3 */
-	       fprintf(ficgp," t \"i%d%d\" ", k2,k);
-	   }else{ /* end ng <> 1 */
-	     if( k !=k2) /* logit p11 is hard to draw */
-	       fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
-	   }
-	   if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
-	     fprintf(ficgp,",");
-	   if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
-	     fprintf(ficgp,",");
-	   i=i+ncovmodel;
-	 } /* end k */
-       } /* end k2 */
-       fprintf(ficgp,"\n set out\n");
-     } /* end jk */
-   } /* end ng */
- /* avoid: */
-   fflush(ficgp); 
+	      ij=1;
+	      for(j=3; j <=ncovmodel-nagesqr; j++){
+		if(ij <=cptcovage) { /* Bug valgrind */
+		  if((j-2)==Tage[ij]) { /* 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++;
+		  }
+		}
+		else
+		  fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+	      }
+	      fprintf(ficgp,")");
+	    }
+	    fprintf(ficgp,")");
+	    if(ng ==2)
+	      fprintf(ficgp," t \"p%d%d\" ", k2,k);
+	    else /* ng= 3 */
+	      fprintf(ficgp," t \"i%d%d\" ", k2,k);
+	  }else{ /* end ng <> 1 */
+	    if( k !=k2) /* logit p11 is hard to draw */
+	      fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
+	  }
+	  if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
+	    fprintf(ficgp,",");
+	  if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
+	    fprintf(ficgp,",");
+	  i=i+ncovmodel;
+	} /* end k */
+      } /* end k2 */
+      fprintf(ficgp,"\n set out\n");
+    } /* end jk */
+  } /* end ng */
+  /* avoid: */
+  fflush(ficgp); 
 }  /* end gnuplot */
 
 
 /*************** Moving average **************/
 /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
-int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
+ int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
    
-  int i, cpt, cptcod;
-  int modcovmax =1;
-  int mobilavrange, mob;
-  int iage=0;
-
-  double sum=0.;
-  double age;
-  double *sumnewp, *sumnewm;
-  double *agemingood, *agemaxgood; /* Currently identical for all covariates */
-  
-  
-  /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose  */
-	/* 		   a covariate has 2 modalities, should be equal to ncovcombmax  *\/ */
-
-  sumnewp = vector(1,ncovcombmax);
-  sumnewm = vector(1,ncovcombmax);
-  agemingood = vector(1,ncovcombmax);	
-  agemaxgood = vector(1,ncovcombmax);
-
-  for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
-		sumnewm[cptcod]=0.;
-		sumnewp[cptcod]=0.;
-		agemingood[cptcod]=0;
-		agemaxgood[cptcod]=0;
-	}
-  if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */
+   int i, cpt, cptcod;
+   int modcovmax =1;
+   int mobilavrange, mob;
+   int iage=0;
+
+   double sum=0.;
+   double age;
+   double *sumnewp, *sumnewm;
+   double *agemingood, *agemaxgood; /* Currently identical for all covariates */
+  
+  
+   /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose  */
+   /* 		   a covariate has 2 modalities, should be equal to ncovcombmax  *\/ */
+
+   sumnewp = vector(1,ncovcombmax);
+   sumnewm = vector(1,ncovcombmax);
+   agemingood = vector(1,ncovcombmax);	
+   agemaxgood = vector(1,ncovcombmax);
+
+   for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+     sumnewm[cptcod]=0.;
+     sumnewp[cptcod]=0.;
+     agemingood[cptcod]=0;
+     agemaxgood[cptcod]=0;
+   }
+   if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */
   
-  if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
-    if(mobilav==1) mobilavrange=5; /* default */
-    else mobilavrange=mobilav;
-    for (age=bage; age<=fage; age++)
-      for (i=1; i<=nlstate;i++)
-				for (cptcod=1;cptcod<=ncovcombmax;cptcod++)
-					mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
-    /* We keep the original values on the extreme ages bage, fage and for 
-       fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
-       we use a 5 terms etc. until the borders are no more concerned. 
-    */ 
-    for (mob=3;mob <=mobilavrange;mob=mob+2){
-      for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
-				for (i=1; i<=nlstate;i++){
-					for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
-						mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
-						for (cpt=1;cpt<=(mob-1)/2;cpt++){
-							mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
-							mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
-						}
-						mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
-					}
-				}
-      }/* end age */
-    }/* end mob */
-  }else
-    return -1;
-  for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
-    /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
-    agemingood[cptcod]=fage-(mob-1)/2;
-    for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
-      sumnewm[cptcod]=0.;
-      for (i=1; i<=nlstate;i++){
-        sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
-      }
-      if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
-				agemingood[cptcod]=age;
-      }else{ /* bad */
-				for (i=1; i<=nlstate;i++){
-					mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
-				} /* i */
-      } /* end bad */
-    }/* age */
-    sum=0.;
-    for (i=1; i<=nlstate;i++){
-      sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
-    }
-    if(fabs(sum - 1.) > 1.e-3) { /* bad */
-      printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod);
-      /* for (i=1; i<=nlstate;i++){ */
-      /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
-      /* } /\* i *\/ */
-    } /* end bad */
-    /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
-		/* From youngest, finding the oldest wrong */
-		agemaxgood[cptcod]=bage+(mob-1)/2;
-		for (age=bage+(mob-1)/2; age<=fage; age++){
-			sumnewm[cptcod]=0.;
-			for (i=1; i<=nlstate;i++){
-				sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
-			}
-			if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
-				agemaxgood[cptcod]=age;
-			}else{ /* bad */
-				for (i=1; i<=nlstate;i++){
-					mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
-				} /* i */
-			} /* end bad */
-		}/* age */
-		sum=0.;
-		for (i=1; i<=nlstate;i++){
-			sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
-		}
-		if(fabs(sum - 1.) > 1.e-3) { /* bad */
-			printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod);
-			/* for (i=1; i<=nlstate;i++){ */
-			/*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
-			/* } /\* i *\/ */
-		} /* end bad */
-		
-		for (age=bage; age<=fage; age++){
-			printf("%d %d ", cptcod, (int)age);
-			sumnewp[cptcod]=0.;
-			sumnewm[cptcod]=0.;
-			for (i=1; i<=nlstate;i++){
-				sumnewp[cptcod]+=probs[(int)age][i][cptcod];
-				sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
-				/* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */
-			}
-			/* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */
-		}
-		/* printf("\n"); */
-    /* } */
-    /* brutal averaging */
-    for (i=1; i<=nlstate;i++){
-      for (age=1; age<=bage; age++){
-				mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
-				/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
-      }	
-      for (age=fage; age<=AGESUP; age++){
-				mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
-				/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
-      }
-    } /* end i status */
-    for (i=nlstate+1; i<=nlstate+ndeath;i++){
-      for (age=1; age<=AGESUP; age++){
-				/*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
-				mobaverage[(int)age][i][cptcod]=0.;
-      }
-    }
-  }/* end cptcod */
-  free_vector(sumnewm,1, ncovcombmax);
-  free_vector(sumnewp,1, ncovcombmax);
-  free_vector(agemaxgood,1, ncovcombmax);
-  free_vector(agemingood,1, ncovcombmax);
-  return 0;
-}/* End movingaverage */
+   if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
+     if(mobilav==1) mobilavrange=5; /* default */
+     else mobilavrange=mobilav;
+     for (age=bage; age<=fage; age++)
+       for (i=1; i<=nlstate;i++)
+	 for (cptcod=1;cptcod<=ncovcombmax;cptcod++)
+	   mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+     /* We keep the original values on the extreme ages bage, fage and for 
+	fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
+	we use a 5 terms etc. until the borders are no more concerned. 
+     */ 
+     for (mob=3;mob <=mobilavrange;mob=mob+2){
+       for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
+	 for (i=1; i<=nlstate;i++){
+	   for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+	     mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
+	     for (cpt=1;cpt<=(mob-1)/2;cpt++){
+	       mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+	       mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+	     }
+	     mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
+	   }
+	 }
+       }/* end age */
+     }/* end mob */
+   }else
+     return -1;
+   for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+     /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
+     if(invalidvarcomb[cptcod]){
+       printf("\nCombination (%d) ignored because no cases \n",cptcod); 
+       continue;
+     }
+
+     agemingood[cptcod]=fage-(mob-1)/2;
+     for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+	 sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+       }
+       if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+	 agemingood[cptcod]=age;
+       }else{ /* bad */
+	 for (i=1; i<=nlstate;i++){
+	   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+	 } /* i */
+       } /* end bad */
+     }/* age */
+     sum=0.;
+     for (i=1; i<=nlstate;i++){
+       sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+     }
+     if(fabs(sum - 1.) > 1.e-3) { /* bad */
+       printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod);
+       /* for (i=1; i<=nlstate;i++){ */
+       /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+       /* } /\* i *\/ */
+     } /* end bad */
+     /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
+     /* From youngest, finding the oldest wrong */
+     agemaxgood[cptcod]=bage+(mob-1)/2;
+     for (age=bage+(mob-1)/2; age<=fage; age++){
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+	 sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+       }
+       if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+	 agemaxgood[cptcod]=age;
+       }else{ /* bad */
+	 for (i=1; i<=nlstate;i++){
+	   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+	 } /* i */
+       } /* end bad */
+     }/* age */
+     sum=0.;
+     for (i=1; i<=nlstate;i++){
+       sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+     }
+     if(fabs(sum - 1.) > 1.e-3) { /* bad */
+       printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod);
+       /* for (i=1; i<=nlstate;i++){ */
+       /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+       /* } /\* i *\/ */
+     } /* end bad */
+		
+     for (age=bage; age<=fage; age++){
+       printf("%d %d ", cptcod, (int)age);
+       sumnewp[cptcod]=0.;
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+	 sumnewp[cptcod]+=probs[(int)age][i][cptcod];
+	 sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+	 /* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */
+       }
+       /* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */
+     }
+     /* printf("\n"); */
+     /* } */
+     /* brutal averaging */
+     for (i=1; i<=nlstate;i++){
+       for (age=1; age<=bage; age++){
+	 mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+	 /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
+       }	
+       for (age=fage; age<=AGESUP; age++){
+	 mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+	 /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
+       }
+     } /* end i status */
+     for (i=nlstate+1; i<=nlstate+ndeath;i++){
+       for (age=1; age<=AGESUP; age++){
+	 /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
+	 mobaverage[(int)age][i][cptcod]=0.;
+       }
+     }
+   }/* end cptcod */
+   free_vector(sumnewm,1, ncovcombmax);
+   free_vector(sumnewp,1, ncovcombmax);
+   free_vector(agemaxgood,1, ncovcombmax);
+   free_vector(agemingood,1, ncovcombmax);
+   return 0;
+ }/* End movingaverage */
  
 
 /************** Forecasting ******************/
@@ -7148,12 +7190,13 @@ int readdata(char datafile[], int firsto
   /*-------- data file ----------*/
   FILE *fic;
   char dummy[]="                         ";
-  int i=0, j=0, n=0;
+  int i=0, j=0, n=0, iv=0;
+  int lstra;
   int linei, month, year,iout;
   char line[MAXLINE], linetmp[MAXLINE];
   char stra[MAXLINE], strb[MAXLINE];
   char *stratrunc;
-  int lstra;
+
 
 
   if((fic=fopen(datafile,"r"))==NULL)    {
@@ -7180,41 +7223,106 @@ int readdata(char datafile[], int firsto
     }
     trimbb(linetmp,line); /* Trims multiple blanks in line */
     strcpy(line, linetmp);
-  
-
+    
+    /* 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;
+				}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; 
+				}
+				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 \
+ 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																						\
+ 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																\
+ Exiting.\n",lval,linei, i,line,j);
+					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																						\
+ 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				\
+ Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
+					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 status */
-	lval=-1;
+      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 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;
-      
+
+      /* Date of Interview */
       strcpy(line,stra);
       cutv(stra, strb,line,' ');
       if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
       }
       else  if( (iout=sscanf(strb,"%s.",dummy)) != 0){
-	month=99;
-	year=9999;
+				month=99;
+				year=9999;
       }else{
-	printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
-	fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
-	return 1;
+				printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
+				fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+				return 1;
       }
       anint[j][i]= (double) year; 
       mint[j][i]= (double)month; 
       strcpy(line,stra);
-    } /* ENd Waves */
-    
+    } /* End loop on waves */
+
+    /* Date of death */
     cutv(stra, strb,line,' '); 
     if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }
@@ -7223,13 +7331,14 @@ int readdata(char datafile[], int firsto
       year=9999;
     }else{
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
-	fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
-	return 1;
+			fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
+			return 1;
     }
     andc[i]=(double) year; 
     moisdc[i]=(double) month; 
     strcpy(line,stra);
     
+    /* Date of birth */
     cutv(stra, strb,line,' '); 
     if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }
@@ -7239,18 +7348,19 @@ int readdata(char datafile[], int firsto
     }else{
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
-	return 1;
+			return 1;
     }
     if (year==9999) {
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
-	return 1;
+			return 1;
 
     }
     annais[i]=(double)(year);
     moisnais[i]=(double)(month); 
     strcpy(line,stra);
-    
+
+    /* Sample weight */
     cutv(stra, strb,line,' '); 
     errno=0;
     dval=strtod(strb,&endptr); 
@@ -7262,22 +7372,44 @@ int readdata(char datafile[], int firsto
     }
     weight[i]=dval; 
     strcpy(line,stra);
+
+    for (iv=nqv;iv>=1;iv--){  /* Loop  on fixed quantitative variables */
+      cutv(stra, strb, line, ' '); 
+      if(strb[0]=='.') { /* Missing value */
+				lval=-1;
+      }else{
+				errno=0;
+				/* what_kind_of_number(strb); */
+				dval=strtod(strb,&endptr);
+				/* 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) constant for all waves. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqv, 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) constant for all waves. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqv, maxwav);fflush(ficlog);
+					return 1;
+				}
+				coqvar[iv][i]=dval; 
+      }
+      strcpy(line,stra);
+    }/* end loop nqv */
     
+    /* Covariate values */
     for (j=ncovcol;j>=1;j--){
       cutv(stra, strb,line,' '); 
-      if(strb[0]=='.') { /* Missing status */
-	lval=-1;
+      if(strb[0]=='.') { /* Missing covariate value */
+				lval=-1;
       }else{
-	errno=0;
-	lval=strtol(strb,&endptr,10); 
-	if( strb[0]=='\0' || (*endptr != '\0')){
-	  printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
-	  fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
-	  return 1;
-	}
+				errno=0;
+				lval=strtol(strb,&endptr,10); 
+				if( strb[0]=='\0' || (*endptr != '\0')){
+					printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
+					fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
+					return 1;
+				}
       }
       if(lval <-1 || lval >1){
-	printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+				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 \
@@ -7286,7 +7418,7 @@ int readdata(char datafile[], int firsto
  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 \
  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 \
@@ -7295,7 +7427,7 @@ int readdata(char datafile[], int firsto
  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 \
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
-	return 1;
+				return 1;
       }
       covar[j][i]=(double)(lval);
       strcpy(line,stra);
@@ -7319,13 +7451,11 @@ int readdata(char datafile[], int firsto
  
   return (0);
   /* endread: */
-    printf("Exiting readdata: ");
-    fclose(fic);
-    return (1);
-
-
-
+	printf("Exiting readdata: ");
+	fclose(fic);
+	return (1);
 }
+
 void removespace(char *str) {
   char *p1 = str, *p2 = str;
   do
@@ -7456,62 +7586,62 @@ int decodemodel ( char model[], int last
         Tvar[k]=0;
       cptcovage=0;
       for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
-	cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
-					 modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ 
-	if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
-	/*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
-	/*scanf("%d",i);*/
-	if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
-	  cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
-	  if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
-	    /* covar is not filled and then is empty */
-	    cptcovprod--;
-	    cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
-	    Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
-	    cptcovage++; /* Sums the number of covariates which include age as a product */
-	    Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
-	    /*printf("stre=%s ", stre);*/
-	  } else if (strcmp(strd,"age")==0) { /* or age*Vn */
-	    cptcovprod--;
-	    cutl(stre,strb,strc,'V');
-	    Tvar[k]=atoi(stre);
-	    cptcovage++;
-	    Tage[cptcovage]=k;
-	  } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
-	    /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
-	    cptcovn++;
-	    cptcovprodnoage++;k1++;
-	    cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-	    Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
-				   because this model-covariate is a construction we invent a new column
-				   ncovcol + k1
-				   If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
-				   Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
-	    cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
-	    Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
-	    Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
-	    Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
-	    k2=k2+2;
-	    Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
-	    Tvar[cptcovt+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
-	    for (i=1; i<=lastobs;i++){
-	      /* Computes the new covariate which is a product of
-		 covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
-	      covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
-	    }
-	  } /* End age is not in the model */
-	} /* End if model includes a product */
-	else { /* no more sum */
-	  /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
-	  /*  scanf("%d",i);*/
-	  cutl(strd,strc,strb,'V');
-	  ks++; /**< Number of simple covariates */
-	  cptcovn++;
-	  Tvar[k]=atoi(strd);
-	}
-	strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
-	/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
-	  scanf("%d",i);*/
+				cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
+																				 modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ 
+				if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+				/*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+				/*scanf("%d",i);*/
+				if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
+					cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+					if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+						/* covar is not filled and then is empty */
+						cptcovprod--;
+						cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+						Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+						cptcovage++; /* Sums the number of covariates which include age as a product */
+						Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+						/*printf("stre=%s ", stre);*/
+					} else if (strcmp(strd,"age")==0) { /* or age*Vn */
+						cptcovprod--;
+						cutl(stre,strb,strc,'V');
+						Tvar[k]=atoi(stre);
+						cptcovage++;
+						Tage[cptcovage]=k;
+					} else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
+						/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+						cptcovn++;
+						cptcovprodnoage++;k1++;
+						cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+						Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
+																	 because this model-covariate is a construction we invent a new column
+																	 ncovcol + k1
+																	 If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
+																	 Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
+						cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
+						Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
+						Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+						Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+						k2=k2+2;
+						Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
+						Tvar[cptcovt+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
+						for (i=1; i<=lastobs;i++){
+							/* Computes the new covariate which is a product of
+								 covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+							covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+						}
+					} /* End age is not in the model */
+				} /* End if model includes a product */
+				else { /* no more sum */
+					/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+					/*  scanf("%d",i);*/
+					cutl(strd,strc,strb,'V');
+					ks++; /**< Number of simple covariates */
+					cptcovn++;
+					Tvar[k]=atoi(strd);
+				}
+				strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
+				/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+					scanf("%d",i);*/
       } /* end of loop + on total covariates */
     } /* end if strlen(modelsave == 0) age*age might exist */
   } /* end if strlen(model == 0) */
@@ -8131,6 +8261,10 @@ int hPijx(double *p, int bage, int fage)
     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) */
@@ -8445,13 +8579,13 @@ int main(int argc, char *argv[])
     }else
       break;
   }
-  if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
-			&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
-    if (num_filled != 8) {
-      printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
+  if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
+			&ftol, &stepm, &ncovcol, &nqv, &ntv, &nqtv, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
+    if (num_filled != 11) {
+      printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1  nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
       printf("but line=%s\n",line);
     }
-    printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
+    printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
   }
   /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
   /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
@@ -8489,8 +8623,8 @@ int main(int argc, char *argv[])
   /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
   /* numlinepar=numlinepar+3; /\* In general *\/ */
   /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
-  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
   fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){
@@ -8519,6 +8653,9 @@ int main(int argc, char *argv[])
 
    
   covar=matrix(0,NCOVMAX,1,n);  /**< used in readdata */
+  coqvar=matrix(1,nqv,1,n);  /**< used in readdata */
+  cotvar=ma3x(1,maxwav,1,ntv,1,n);  /**< used in readdata */
+  cotqvar=ma3x(1,maxwav,1,ntqv,1,n);  /**< used in readdata */
   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
   /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
      v1+v2*age+v2*v3 makes cptcovn = 3
@@ -8786,7 +8923,7 @@ Please run with mle=-1 to get a correct
 /* Main decodemodel */
 
 
-  if(decodemodel(model, lastobs) == 1)
+  if(decodemodel(model, lastobs) == 1) /* In order to get Tvar[k] V4+V3+V5 p Tvar[1]@3  = {4, 3, 5}*/
     goto end;
 
   if((double)(lastobs-imx)/(double)imx > 1.10){
@@ -9014,7 +9151,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 		 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,	\
-	      firstpass, lastpass,  stepm,  weightopt, model);
+							firstpass, lastpass,  stepm,  weightopt, model);
 
   fprintf(fichtm,"\n");
   fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
@@ -9043,7 +9180,7 @@ Interval (in months) between two waves:
     ageexmed=vector(1,n);
     agecens=vector(1,n);
     dcwave=ivector(1,n);
- 
+		
     for (i=1; i<=imx; i++){
       dcwave[i]=-1;
       for (m=firstpass; m<=lastpass; m++)
@@ -9543,9 +9680,9 @@ Please run with mle=-1 to get a correct
     ungetc(c,ficpar);
     
     fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
-    fprintf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
-    fprintf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
-    fprintf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+    fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+    fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+    fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
     /* day and month of proj2 are not used but only year anproj2.*/
     
     
@@ -9648,8 +9785,8 @@ Please run with mle=-1 to get a correct
       back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
       fclose(ficresplb);
 
-      /* hBijx(p, bage, fage, mobaverage); */
-      /* fclose(ficrespijb); */
+      hBijx(p, bage, fage, mobaverage);
+      fclose(ficrespijb);
       free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
 
       /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
@@ -9893,6 +10030,9 @@ Please run with mle=-1 to get a correct
     free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
+    free_ma3x(cotqvar,1,maxwav,1,ntqv,1,n);
+    free_ma3x(cotvar,1,maxwav,1,ntv,1,n);
+    free_matrix(coqvar,1,maxwav,1,n);
     free_matrix(covar,0,NCOVMAX,1,n);
     free_matrix(matcov,1,npar,1,npar);
     free_matrix(hess,1,npar,1,npar);