--- imach/src/imach.c	2022/09/02 14:26:02	1.337
+++ imach/src/imach.c	2022/09/11 07:53:11	1.340
@@ -1,6 +1,28 @@
-/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $
+/* $Id: imach.c,v 1.340 2022/09/11 07:53:11 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.340  2022/09/11 07:53:11  brouard
+  Summary: Version imach 0.99r37
+
+  * imach.c (Module): Adding timevarying products of any kinds,
+  should work before shifting cotvar from ncovcol+nqv columns in
+  order to have a correspondance between the column of cotvar and
+  the id of column.
+
+  Revision 1.339  2022/09/09 17:55:22  brouard
+  Summary: version 0.99r37
+
+  * imach.c (Module): Many improvements for fixing products of fixed
+  timevarying as well as fixed * fixed, and test with quantitative
+  covariate.
+
+  Revision 1.338  2022/09/04 17:40:33  brouard
+  Summary: 0.99r36
+
+  * imach.c (Module): Now the easy runs i.e. without result or
+  model=1+age only did not work. The defautl combination should be 1
+  and not 0 because everything hasn't been tranformed yet.
+
   Revision 1.337  2022/09/02 14:26:02  brouard
   Summary: version 0.99r35
 
@@ -1294,12 +1316,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $ */
+/* $Id: imach.c,v 1.340 2022/09/11 07:53:11 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="September 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.337 $ $Date: 2022/09/02 14:26:02 $"; 
+char fullversion[]="$Revision: 1.340 $ $Date: 2022/09/11 07:53:11 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1315,6 +1337,7 @@ int cptcovprodnoage=0; /**< Number of co
 int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
 int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
 int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
+int ncovvt=0; /* Total number of effective (wave) varying covariates (dummy or quantitative or products [without age]) in the model */
 int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
 int nsd=0; /**< Total number of single dummy variables (output) */
 int nsq=0; /**< Total number of single quantitative variables (output) */
@@ -1330,7 +1353,8 @@ int npar=NPARMAX; /* Number of parameter
 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 nqv=0, ntv=0, nqtv=0;    /* Total number of quantitative variables, time variable (dummy), quantitative and time variable*/
+int ncovcolt=0; /* ncovcolt=ncovcol+nqv+ntv+nqtv; total of covariates in the data, not in the model equation*/ 
 int popbased=0;
 
 int *wav; /* Number of waves for this individuual 0 is possible */
@@ -1517,7 +1541,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv
 /*TnsdVar[Tvar]   1   2                              3 */ 
 /*Tvaraff[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
 /*TvarsD[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
-/*TvarsDind[k]    2   3                              9 */ /* position K of single dummy cova */
+/*TvarsDind[nsd]  2   3                              9 */ /* position K of single dummy cova */
 /*    nsq      1                     2                 */ /* Counting single quantit tv */
 /* TvarsQ[k]   5                     2                 */ /* Number of single quantitative cova */
 /* TvarsQind   1                     6                 */ /* position K of single quantitative cova */
@@ -1572,7 +1596,13 @@ int *TvarVD; /* TvarVD[1]=V5 in V5+V4+V3
 int *TvarVDind; /* TvarVDind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
 int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
 int *TvarVQind; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
-
+int *TvarVV; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */
+int *TvarVVind; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */
+      /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+      /* model V1+V3+age*V1+age*V3+V1*V3 */
+      /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+      /* TvarVV={3,1,3}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */	       
+      /* TvarVVind={2,5,5}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */	       
 int *Tvarsel; /**< Selected covariates for output */
 double *Tvalsel; /**< Selected modality value of covariate for output */
 int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
@@ -2510,7 +2540,8 @@ void powell(double p[], double **xi, int
   xits=vector(1,n); 
   *fret=(*func)(p); 
   for (j=1;j<=n;j++) pt[j]=p[j]; 
-  rcurr_time = time(NULL);  
+  rcurr_time = time(NULL);
+  fp=(*fret); /* Initialisation */
   for (*iter=1;;++(*iter)) { 
     ibig=0; 
     del=0.0; 
@@ -2824,7 +2855,7 @@ void powell(double p[], double **xi, int
   
   double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres)
   {
-    /**< Computes the prevalence limit in each live state at age x and for covariate combination ij 
+    /**< Computes the prevalence limit in each live state at age x and for covariate combination ij . Nicely done
      *   (and selected quantitative values in nres)
      *  by left multiplying the unit
      *  matrix by transitions matrix until convergence is reached with precision ftolpl 
@@ -3272,12 +3303,12 @@ double **pmij(double **ps, double *cov,
   for(i=1; i<= nlstate; i++){
     s1=0;
     for(j=1; j<i; j++){
+      /* printf("debug1 %d %d ps=%lf exp(ps)=%lf \n",i,j,ps[i][j],exp(ps[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++){
+      /* printf("debug2 %d %d ps=%lf exp(ps)=%lf \n",i,j,ps[i][j],exp(ps[i][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.);
@@ -3855,6 +3886,9 @@ double func( double *x)
 {
   int i, ii, j, k, mi, d, kk, kf=0;
   int ioffset=0;
+  int ipos=0,iposold=0,ncovv=0;
+
+  double cotvarv, cotvarvold;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
   double lli; /* Individual log likelihood */
@@ -3902,16 +3936,45 @@ double func( double *x)
 	 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]
+	 And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv-ncovcol-nqv][i] because (-ncovcol-nqv) in the data
 	 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++){  /* Varying with waves */
       /* Wave varying (but not age varying) */
-	for(k=1; k <= ncovv ; k++){ /* Varying  covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3)  Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/
-	  /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? */
-	  cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
-	}
+	/* for(k=1; k <= ncovv ; k++){ /\* Varying  covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3)  Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*\/ */
+	/*   /\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? *\/ */
+	/*   cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */
+	/* } */
+	for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/
+	  itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
+	  ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+	  if(TvarFind[itv]==0){ /* Not a fixed covariate */
+	    cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]-ncovcol-nqv][i];  /* cotvar[wav][ntv+iv][i] */
+	  }else{ /* fixed covariate */
+	    cotvarv=covar[Tvar[TvarFind[itv]]][i];
+	  }
+	  if(ipos!=iposold){ /* Not a product or first of a product */
+	    cotvarvold=cotvarv;
+	  }else{ /* A second product */
+	    cotvarv=cotvarv*cotvarvold;
+	  }
+	  iposold=ipos;
+	  cov[ioffset+ipos]=cotvarv;
+	}
+	/* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
+	/*   iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
+	/*   cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
+	/*   k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
+	/*   cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
+	/*   printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
+	/* } */
+	/* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */
+	/*   iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
+	/*   /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
+	/*   cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */
+	/* } */
+	/* for products of time varying to be done */
 	for (ii=1;ii<=nlstate+ndeath;ii++)
 	  for (j=1;j<=nlstate+ndeath;j++){
 	    oldm[ii][j]=(ii==j ? 1.0 : 0.0);
@@ -3930,7 +3993,7 @@ double func( double *x)
 	    if(!FixedV[Tvar[Tage[kk]]])
 	      cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
 	    else
-	      cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
+	      cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ 
 	  }
 	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 		       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
@@ -4026,7 +4089,7 @@ double func( double *x)
 	} 
 	/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
 	/*if(lli ==000.0)*/
-	/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+	/* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
 	ipmx +=1;
 	sw += weight[i];
 	ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
@@ -4043,7 +4106,7 @@ double func( double *x)
 	cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];
       for(mi=1; mi<= wav[i]-1; mi++){
 	for(k=1; k <= ncovv ; k++){
-	  cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
+	  cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; /* Cotvar starts at ntv */
 	}
 	for (ii=1;ii<=nlstate+ndeath;ii++)
 	  for (j=1;j<=nlstate+ndeath;j++){
@@ -4090,7 +4153,10 @@ double func( double *x)
 	  if(nagesqr==1)
 	    cov[3]= agexact*agexact;
 	  for (kk=1; kk<=cptcovage;kk++) {
-	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+	    if(!FixedV[Tvar[Tage[kk]]])
+	      cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+	    else
+	      cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ 
 	  }
 	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 		       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
@@ -4146,7 +4212,7 @@ double func( double *x)
 	ipmx +=1;
 	sw += weight[i];
 	ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-/* 	printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+	/* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       } /* end of wave */
     } /* end of individual */
   }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
@@ -4165,7 +4231,10 @@ double func( double *x)
 	  if(nagesqr==1)
 	    cov[3]= agexact*agexact;
 	  for (kk=1; kk<=cptcovage;kk++) {
-	    cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+	    if(!FixedV[Tvar[Tage[kk]]])
+	      cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+	    else
+	      cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ 
 	  }
 	
 	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
@@ -4196,6 +4265,9 @@ double funcone( double *x)
   /* Same as func but slower because of a lot of printf and if */
   int i, ii, j, k, mi, d, kk, kf=0;
   int ioffset=0;
+  int ipos=0,iposold=0,ncovv=0;
+
+  double cotvarv, cotvarvold;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
   double lli; /* Individual log likelihood */
@@ -4228,6 +4300,9 @@ double funcone( double *x)
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
     /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
     for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+      /* printf("Debug3 TvarFind[%d]=%d",kf, TvarFind[kf]); */
+      /* printf(" Tvar[TvarFind[kf]]=%d", Tvar[TvarFind[kf]]); */
+      /* printf(" i=%d covar[Tvar[TvarFind[kf]]][i]=%f\n",i,covar[Tvar[TvarFind[kf]]][i]); */
       cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
 /*    cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i];  */
 /*    cov[2+6]=covar[Tvar[6]][i];  */
@@ -4247,7 +4322,7 @@ double funcone( double *x)
 	 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]
+	 And the iv th varying covariate in the DATA 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]
       */
@@ -4261,17 +4336,59 @@ double funcone( double *x)
     
 
     for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */
-    /* Wave varying (but not age varying) */
-      for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/
-	/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
-	cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
-      }
-      /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
-      /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
-      /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
-      /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
-      /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
-      /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
+      /* Wave varying (but not age varying) *//* V1+V3+age*V1+age*V3+V1*V3 with V4 tv and V5 tvq k= 1 to 5 and extra at V(5+1)=6 for V1*V3 */
+      /* for(k=1; k <= ncovv ; k++){ /\* Varying  covariates (single and product but no age )*\/ */
+      /* 	/\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; *\/ */
+      /* 	cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */
+      /* } */
+      
+      /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+      /* model V1+V3+age*V1+age*V3+V1*V3 */
+      /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+      /*  TvarVV[1]=V3 (first time varying in the model equation, TvarVV[2]=V1 (in V1*V3) TvarVV[3]=3(V3)  */
+      /* We need the position of the time varying or product in the model */
+      /* TvarVVind={2,5,5}, for V3 at position 2 and then the product V1*V3 is decomposed into V1 and V3 but at same position 5 */	       
+      /* TvarVV gives the variable name */
+      /* Other example V1 + V3 + V5 + age*V1  + age*V3 + age*V5 + V1*V3  + V3*V5  + V1*V5 
+      *      k=         1   2     3     4         5        6        7       8        9
+      *  varying            1     2                                 3       4        5
+      *  ncovv              1     2                                3 4     5 6      7 8
+      *  TvarVV            V3     5                                1 3     3 5      1 5
+      * TvarVVind           2     3                                7 7     8 8      9 9
+      * TvarFind[k]     1   0     0     0         0        0        0       0        0
+      * cotvar starts at ntv=2 (because of V3 V4)
+      */
+      for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age) including individual from products */
+	itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
+	ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+	if(TvarFind[itv]==0){ /* Not a fixed covariate */
+	  cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]-ncovcol-nqv][i];  /* cotvar[wav][ntv+iv][i] */
+	}else{ /* fixed covariate */
+	  cotvarv=covar[Tvar[TvarFind[itv]]][i];
+	}
+	if(ipos!=iposold){ /* Not a product or first of a product */
+	  cotvarvold=cotvarv;
+	}else{ /* A second product */
+	  cotvarv=cotvarv*cotvarvold;
+	}
+	iposold=ipos;
+	cov[ioffset+ipos]=cotvarv;
+	/* For products */
+      }
+      /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */
+      /* 	iv=TvarVDind[itv]; /\* iv, position in the model equation of time varying covariate itv *\/ */
+      /* 	/\*         "V1+V3+age*V1+age*V3+V1*V3" with V3 time varying *\/ */
+      /* 	/\*           1  2   3      4      5                         *\/ */
+      /* 	/\*itv           1                                           *\/ */
+      /* 	/\* TvarVInd[1]= 2                                           *\/ */
+      /* 	/\* iv= Tvar[Tmodelind[itv]]-ncovcol-nqv;  /\\* Counting the # varying covariate from 1 to ntveff *\\/ *\/ */
+      /* 	/\* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; *\/ */
+      /* 	/\* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; *\/ */
+      /* 	/\* k=ioffset-2-nagesqr-cptcovage+itv; /\\* position in simple model *\\/ *\/ */
+      /* 	/\* cov[ioffset+iv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; *\/ */
+      /* 	cov[ioffset+iv]=cotvar[mw[mi][i]][itv][i]; */
+      /* 	/\* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][itv][i]=%f\n", i, mi, itv, TvarVDind[itv],cotvar[mw[mi][i]][itv][i]); *\/ */
+      /* } */
       /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */
       /* 	iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
       /* 	/\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
@@ -4353,7 +4470,26 @@ double funcone( double *x)
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-      /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
+      printf("Funcone num[i]=%ld i=%6d V= ", num[i], i);
+      for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+	printf("%g",covar[Tvar[TvarFind[kf]]][i]);
+      }
+      for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age) including individual from products */
+	ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+	if(ipos!=iposold){ /* Not a product or first of a product */
+	  printf(" %g",cov[ioffset+ipos]);
+	}else{
+	  printf("*");
+	}
+	iposold=ipos;
+      }
+      for (kk=1; kk<=cptcovage;kk++) {
+	if(!FixedV[Tvar[Tage[kk]]])
+	  printf(" %g*age",covar[Tvar[Tage[kk]]][i]);
+	else
+	  printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]);
+      }
+      printf(" s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
       if(globpr){
 	fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \
@@ -5129,7 +5265,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 		/* }else  */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/
 		/* if( iind >=imx-3) printf("Searching error iind=%d Tvaraff[z1]=%d covar[Tvaraff[z1]][iind]=%.f TnsdVar[Tvaraff[z1]]=%d, cptcoveff=%d, cptcovs=%d \n",iind, Tvaraff[z1], covar[Tvaraff[z1]][iind],TnsdVar[Tvaraff[z1]],cptcoveff, cptcovs); */
 		if(Tvaraff[z1]<1 || Tvaraff[z1]>=NCOVMAX)
-		  printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
+		  printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=1+age+%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
 		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */
 		  /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */
 		  bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
@@ -5150,7 +5286,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	      if(anyvaryingduminmodel==1){ /* Some are varying covariates */
 		for (z1=1; z1<=cptcoveff; z1++) {
 		  if( Fixed[Tmodelind[z1]]==1){
-		    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+		    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; /* Good */
 		    if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's 
 										      value is -1, we don't select. It differs from the 
 										      constant and age model which counts them. */
@@ -5231,7 +5367,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
 	fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
 	fprintf(ficlog, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcovs; z1++){
+	for (z1=1; z1<=cptcoveff; z1++){
 	  if(!FixedV[Tvaraff[z1]]){
 	    printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
 	    fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
@@ -5439,7 +5575,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	printf("#  This combination (%d) is not valid and no result will be produced\n",j1);
 	invalidvarcomb[j1]=1;
       }else{
-	fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
+	fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced (or no resultline).</p>",j1);
 	invalidvarcomb[j1]=0;
       }
       fprintf(ficresphtmfr,"</table>\n");
@@ -5996,12 +6132,14 @@ void  concatwav(int wav[], int **dh, int
    /* Loop on covariates without age and products and no quantitative variable */
    for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
      for (j=-1; (j < maxncov); j++) Ndum[j]=0;
-     if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
+     printf("Testing k=%d, cptcovt=%d\n",k, cptcovt);
+     if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */ 
        switch(Fixed[k]) {
        case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
 	 modmaxcovj=0;
 	 modmincovj=0;
 	 for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
+	   /* printf("Waiting for error tricode Tvar[%d]=%d i=%d (int)(covar[Tvar[k]][i]=%d\n",k,Tvar[k], i, (int)(covar[Tvar[k]][i])); */
 	   ij=(int)(covar[Tvar[k]][i]);
 	   /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
 	    * If product of Vn*Vm, still boolean *:
@@ -7208,9 +7346,9 @@ To be simple, these graphs help to under
 
        /* Including quantitative variables of the resultline to be done */
        for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */
-	 printf("Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
-	 fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
-	 /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */
+	 printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);
+	 fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);
+	 /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */
 	 if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline  */
 	   if(Fixed[modelresult[nres][z1]]==0){ /* Fixed referenced to model equation */
 	     fprintf(ficresprob,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
@@ -7581,7 +7719,8 @@ void printinghtml(char fileresu[], char
    fprintf(fichtm," \n<ul>");
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */
-     k1= TKresult[nres];
+     k1=TKresult[nres];
+     if(TKresult[nres]==0)k1=1; /* To be checked for no result */
    /* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
    /*   if(m != 1 && TKresult[nres]!= k1) */
    /*     continue; */
@@ -7624,7 +7763,8 @@ void printinghtml(char fileresu[], char
 
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */
-     k1= TKresult[nres];
+     k1=TKresult[nres];
+     if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
    /* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
    /*   if(m != 1 && TKresult[nres]!= k1) */
    /*     continue; */
@@ -7649,7 +7789,7 @@ void printinghtml(char fileresu[], char
 	 /* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
        }
        /* if(nqfveff+nqtveff 0) */ /* Test to be done */
-       fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
+       fprintf(fichtm," (model=1+age+%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
        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); 
@@ -7685,14 +7825,15 @@ divided by h: <sub>h</sub>P<sub>ij</sub>
      /* Period (forward 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 for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
-       fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"P_"),subdirf2(optionfilefiname,"P_"));
+       fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
       fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }
      if(prevbcast==1){
        /* Backward prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+	 fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+	 fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJB_"),subdirf2(optionfilefiname,"PIJB_"));
+	 fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
        }
      }
      if(prevfcast==1){
@@ -7769,7 +7910,7 @@ See page 'Matrix of variance-covariance
 /*  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n */
 /* 	<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); */
+/*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=1+age+%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
    fflush(fichtm);
 
    m=pow(2,cptcoveff);
@@ -7782,7 +7923,7 @@ See page 'Matrix of variance-covariance
    fprintf(fichtm," \n<ul>");
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */
-     k1= TKresult[nres];
+     k1=TKresult[nres];
      /* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
      /* if(m != 1 && TKresult[nres]!= k1) */
      /*   continue; */
@@ -7812,7 +7953,8 @@ See page 'Matrix of variance-covariance
 
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */
-     k1= TKresult[nres];
+     k1=TKresult[nres];
+     if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
      /* for(k1=1; k1<=m;k1++){ */
      /* if(m != 1 && TKresult[nres]!= k1) */
      /*   continue; */
@@ -7832,7 +7974,7 @@ See page 'Matrix of variance-covariance
 	 /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
        }
 
-       fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
+       fprintf(fichtm," (model=1+age+%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
 
        if(invalidvarcomb[k1]){
 	 fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); 
@@ -7933,6 +8075,7 @@ void printinggnuplot(char fileresu[], ch
     /* for (k1=1; k1<= m ; k1 ++){ /\* For each valid combination of covariate *\/ */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
 	k1=TKresult[nres];
+	if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
 	/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
 	/* if(m != 1 && TKresult[nres]!= k1) */
 	/*   continue; */
@@ -7973,7 +8116,7 @@ void printinggnuplot(char fileresu[], ch
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
 	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
 	/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
-	fprintf(ficgp,"set title \"Alive state %d %s model=%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model);
+	fprintf(ficgp,"set title \"Alive state %d %s model=1+age+%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
 	/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
       /* k1-1 error should be nres-1*/
@@ -8080,6 +8223,7 @@ void printinggnuplot(char fileresu[], ch
   /* for (k1=1; k1<= m ; k1 ++){   */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
@@ -8151,6 +8295,7 @@ void printinggnuplot(char fileresu[], ch
   /* for (k1=1; k1<= m ; k1 ++){ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
 
@@ -8212,6 +8357,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   /* for (k1=1; k1<=m; k1++){    /\* For each covariate and each value *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
@@ -8269,6 +8415,7 @@ set ter svg size 640, 480\nunset log y\n
   /* for (k1=1; k1<= m ; k1++){ /\* For each covariate combination if any *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
@@ -8334,6 +8481,7 @@ set ter svg size 640, 480\nunset log y\n
   /* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      k1=TKresult[nres];
+     if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
      /* if(m != 1 && TKresult[nres]!= k1) */
      /*  continue; */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */
@@ -8391,6 +8539,7 @@ set ter svg size 640, 480\nunset log y\n
     /* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */
@@ -8453,6 +8602,7 @@ set ter svg size 640, 480\nunset log y\n
     /* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       /* if(m != 1 && TKresult[nres]!= k1) */
       /* 	continue; */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
@@ -8574,6 +8724,7 @@ set ter svg size 640, 480\nunset log y\n
     /* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      k1=TKresult[nres];
+     if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
        /* if(m != 1 && TKresult[nres]!= k1) */
        /* 	continue; */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
@@ -8651,7 +8802,7 @@ set ter svg size 640, 480\nunset log y\n
 	    kl=0;
 	    strcpy(gplotcondition,"(");
 	    for (k=1; k<=cptcovs; k++){    /* For each covariate k of the resultline, get corresponding value lv for combination k1 */
-	      if(Dummy[Tvresult[nres][k]]==0){  /* To be verified */
+	      if(Dummy[modelresult[nres][k]]==0){  /* To be verified */
 		/* for (k=1; k<=cptcoveff; k++){    /\* For each covariate writing the chain of conditions *\/ */
 		/* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
 		/* lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
@@ -8666,7 +8817,7 @@ set ter svg size 640, 480\nunset log y\n
 		/* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
 		sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,Tvresult[nres][k], kl+1,TinvDoQresult[nres][Tvresult[nres][k]]);
 		kl++;
-		if(k <cptcoveff && cptcoveff>1)
+		if(k <cptcovs && cptcovs>1)
 		  sprintf(gplotcondition+strlen(gplotcondition)," && ");
 	      }
 	    }
@@ -8733,13 +8884,14 @@ set ter svg size 640, 480\nunset log y\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,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
-    fprintf(ficgp,"#model=%s \n",model);
+    fprintf(ficgp,"#model=1+age+%s \n",model);
     fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
     fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
     /* for(k1=1; k1 <=m; k1++)  /\* For each combination of covariate *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */
-      k1= TKresult[nres];
+      k1=TKresult[nres];
+      if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
       fprintf(ficgp,"\n\n# Resultline k1=%d ",k1);
       strcpy(gplotlabel,"(");
       /*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/
@@ -9555,6 +9707,7 @@ void prevforecast(char fileres[], double
 
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
        k=TKresult[nres];
+       if(TKresult[nres]==0) k=1; /* To be checked for noresult */
      /* for(k=1; k<=i1;k++){ /\* We find the combination equivalent to result line values of dummies *\/ */
       if(i1 != 1 && TKresult[nres]!= k)
 	continue;
@@ -9615,6 +9768,7 @@ void prevforecast(char fileres[], double
    
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      k=TKresult[nres];
+     if(TKresult[nres]==0) k=1; /* To be checked for noresult */
     /* for(k=1; k<=i1;k++){ */
     /*    if(i1 != 1 && TKresult[nres]!= k) */
     /* 	 continue; */
@@ -9622,9 +9776,9 @@ void prevforecast(char fileres[], double
        printf("\n#****** ");
        fprintf(ficlog,"\n#****** ");
        for (j=1; j<= cptcovs; j++){ /* For each selected (single) quantitative value */
-	 printf(" V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
-	 fprintf(ficresvbl," V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
-	 fprintf(ficlog," V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+	 printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+	 fprintf(ficresvbl," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+	 fprintf(ficlog," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
        /* for(j=1;j<=cptcoveff;j++) { */
        /* 	 fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
        /* 	 fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
@@ -10120,7 +10274,9 @@ int readdata(char datafile[], int firsto
     printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
     fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
   }
-
+  
+  ncovcolt=ncovcol+nqv+ntv+nqtv; /* total of covariates in the data, not in the model equation */
+  
   if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
     fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
@@ -10520,8 +10676,8 @@ int decoderesult( char resultline[], int
 	}
       }
       if(match == 0){
-	printf("Error in result line (Dummy single): V%d is missing in result: %s according to model=%s. Tvar[k1=%d]=%d is different from Tvarsel[k2=%d]=%d.\n",Tvar[k1], resultline, model,k1, Tvar[k1], k2, Tvarsel[k2]);
-	fprintf(ficlog,"Error in result line (Dummy single): V%d is missing in result: %s according to model=%s\n",Tvar[k1], resultline, model);
+	printf("Error in result line (Dummy single): V%d is missing in result: %s according to model=1+age+%s. Tvar[k1=%d]=%d is different from Tvarsel[k2=%d]=%d.\n",Tvar[k1], resultline, model,k1, Tvar[k1], k2, Tvarsel[k2]);
+	fprintf(ficlog,"Error in result line (Dummy single): V%d is missing in result: %s according to model=1+age+%s\n",Tvar[k1], resultline, model);
 	return 1;
       }
     }else if(Typevar[k1]==1){ /* Product with age We want to get the position k2 in the resultline of the product k1 in the model line*/
@@ -10537,8 +10693,8 @@ int decoderesult( char resultline[], int
 	}
       }
       if(match == 0){
-	printf("Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
-	fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
+	printf("Error in result line (Product with age): V%d is missing in result: %s according to model=1+age+%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
+	fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=1+age+%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
       return 1;
       }
     }else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/
@@ -10553,8 +10709,8 @@ int decoderesult( char resultline[], int
 	}
       }
       if(match == 0){
-	printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
-	fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
+	printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
+	fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
 	return 1;
       }
       match=0;
@@ -10567,8 +10723,8 @@ int decoderesult( char resultline[], int
 	}
       }
       if(match == 0){
-	printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][2], resultline, model);
-	fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=%s\n",Tvardk[k1][2], resultline, model);
+	printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
+	fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
 	return 1;
       }
     }/* End of testing */
@@ -10588,12 +10744,12 @@ int decoderesult( char resultline[], int
       }
     }
     if(match == 0){
-      printf("Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
-      fprintf(ficlog,"Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
+      printf("Error in result line: variable V%d is missing in model; result: %s, model=1+age+%s\n",Tvarsel[k2], resultline, model);
+      fprintf(ficlog,"Error in result line: variable V%d is missing in model; result: %s, model=1+age+%s\n",Tvarsel[k2], resultline, model);
       return 1;
     }else if(match > 1){
-      printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
-      fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
+      printf("Error in result line: %d doubled; result: %s, model=1+age+%s\n",k2, resultline, model);
+      fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=1+age+%s\n",k2, resultline, model);
       return 1;
     }
   }
@@ -10701,8 +10857,9 @@ int decodemodel( char model[], int lasto
 	* - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
 	* - cptcovage number of covariates with age*products =2
 	* - cptcovs number of simple covariates
+	* ncovcolt=ncovcol+nqv+ntv+nqtv total of covariates in the data, not in the model equation
 	* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
-	*     which is a new column after the 9 (ncovcol) variables. 
+	*     which is a new column after the 9 (ncovcol+nqv+ntv+nqtv) variables. 
 	* - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual
 	* - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
 	*    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
@@ -10725,18 +10882,18 @@ int decodemodel( char model[], int lasto
       return 1;
     }
     if (strstr(model,"v") !=0){
-      printf("Error. 'v' must be in upper case 'V' model=%s ",model);
-      fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog);
+      printf("Error. 'v' must be in upper case 'V' model=1+age+%s ",model);
+      fprintf(ficlog,"Error. 'v' must be in upper case model=1+age+%s ",model);fflush(ficlog);
       return 1;
     }
     strcpy(modelsav,model); 
     if ((strpt=strstr(model,"age*age")) !=0){
-      printf(" strpt=%s, model=%s\n",strpt, model);
+      printf(" strpt=%s, model=1+age+%s\n",strpt, model);
       if(strpt != model){
-	printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+	printf("Error in model: 'model=1+age+%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model);
-	fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+	fprintf(ficlog,"Error in model: 'model=1+age+%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model); fflush(ficlog);
 	return 1;
@@ -10846,12 +11003,12 @@ int decodemodel( char model[], int lasto
 	    cptcovn++;
 	    cptcovprodnoage++;k1++;
 	    cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-	    Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
+	    Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* ncovcolt+k1; For model-covariate k tells which data-covariate to use but
 						because this model-covariate is a construction we invent a new column
 						which is after existing variables ncovcol+nqv+ntv+nqtv + k1
 						If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
 						thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
-						Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
+						Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=3 etc */
 	    /* Please remark that the new variables are model dependent */
 	    /* If we have 4 variable but the model uses only 3, like in
 	     * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
@@ -10860,7 +11017,7 @@ int decodemodel( char model[], int lasto
 	     * Tage[kk]    [1]= 2           [2]=5      [3]=6                  kk=1 to cptcovage=3
 	     * Tvar[Tage[kk]][1]=2          [2]=2      [3]=3
 	     */
-	    Typevar[k]=2;  /* 2 for double fixed dummy covariates */
+	    Typevar[k]=2;  /* 2 for product */
 	    cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
 	    Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
 	    Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
@@ -10873,11 +11030,13 @@ int decodemodel( char model[], int lasto
 	    /* Tvar[cptcovt+k2+1]=Tvard[k1][2];  /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
             /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
 	    /*                     1  2   3      4     5 | Tvar[5+1)=1, Tvar[7]=2   */
-	    for (i=1; i<=lastobs;i++){
+	    if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+	      for (i=1; i<=lastobs;i++){/* For fixed product */
 	      /* 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];
-	    }
+	      covar[ncovcolt+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+	      }
+	    } /*End of FixedV */
 	  } /* End age is not in the model */
 	} /* End if model includes a product */
 	else { /* not a product */
@@ -10929,7 +11088,7 @@ Typevar: 0 for simple covariate (dummy,
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
   for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
-  for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
+  for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt */
     if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;
       Dummy[k]= 0;
@@ -10944,7 +11103,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarFind[ncovf]=k;
       TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-    }else if( Tvar[k] <=ncovcol &&  Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */
+    /* }else if( Tvar[k] <=ncovcol &&  Typevar[k]==2){ /\* Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol *\/ */
+    }else if( Tposprod[k]>0  &&  Typevar[k]==2 && FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* Needs a fixed product Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol */
       Fixed[k]= 0;
       Dummy[k]= 0;
       ncoveff++;
@@ -10970,6 +11130,13 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
       TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
     }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */
+      /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+      /* model V1+V3+age*V1+age*V3+V1*V3 */
+      /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+      ncovvt++;
+      TvarVV[ncovvt]=Tvar[k];  /*  TvarVV[1]=V3 (first time varying in the model equation  */
+      TvarVVind[ncovvt]=k;    /*  TvarVVind[1]=2 (second position in the model equation  */
+
       Fixed[k]= 1;
       Dummy[k]= 0;
       ntveff++; /* Only simple time varying dummy variable */
@@ -10987,6 +11154,13 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
       printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv);
     }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv  && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/
+      /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+      /* model V1+V3+age*V1+age*V3+V1*V3 */
+      /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+      ncovvt++;
+      TvarVV[ncovvt]=Tvar[k];  /*  TvarVV[1]=V3 (first time varying in the model equation  */
+      TvarVVind[ncovvt]=k;  /*  TvarVV[1]=V3 (first time varying in the model equation  */
+      
       Fixed[k]= 1;
       Dummy[k]= 1;
       nqtveff++;
@@ -11033,10 +11207,21 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       	modell[k].subtype= APVQ;		/*	Product age * varying quantitative */
       	/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
       }
-    }else if (Typevar[k] == 2) {  /* product without age */
-      k1=Tposprod[k];
-      if(Tvard[k1][1] <=ncovcol){
-    	if(Tvard[k1][2] <=ncovcol){
+    }else if (Typevar[k] == 2) {  /* product Vn * Vm without age, V1+V3+age*V1+age*V3+V1*V3 looking at V1*V3, Typevar={0, 0, 1, 1, 2}, k=5, V1 is fixed, V3 is timevary, V5 is a product  */
+      /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+      /* model V1+V3+age*V1+age*V3+V1*V3 */
+      /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+      k1=Tposprod[k];  /* Position in the products of product k, Tposprod={0, 0, 0, 0, 1} k1=1 first product but second time varying because of V3 */
+      ncovvt++;
+      TvarVV[ncovvt]=Tvard[k1][1];  /*  TvarVV[2]=V1 (because TvarVV[1] was V3, first time varying covariates */
+      TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+      ncovvt++;
+      TvarVV[ncovvt]=Tvard[k1][2];  /*  TvarVV[3]=V3 */
+      TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+
+
+      if(Tvard[k1][1] <=ncovcol){ /* Vn is dummy fixed, (Tvard[1][1]=V1), (Tvard[1][1]=V3 time varying) */
+    	if(Tvard[k1][2] <=ncovcol){ /* Vm is dummy fixed */
     	  Fixed[k]= 1;
     	  Dummy[k]= 0;
     	  modell[k].maintype= FTYPE;
@@ -11044,23 +11229,23 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
     	  ncovf++; /* Fixed variables without age */
     	  TvarF[ncovf]=Tvar[k];
     	  TvarFind[ncovf]=k;
-    	}else if(Tvard[k1][2] <=ncovcol+nqv){
-    	  Fixed[k]= 0;  /* or 2 ?*/
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){ /* Vm is quanti fixed */
+    	  Fixed[k]= 0;  /* Fixed product */
     	  Dummy[k]= 1;
     	  modell[k].maintype= FTYPE;
     	  modell[k].subtype= FPDQ;		/*	Product fixed dummy * fixed quantitative */
     	  ncovf++; /* Varying variables without age */
     	  TvarF[ncovf]=Tvar[k];
     	  TvarFind[ncovf]=k;
-    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is a time varying dummy covariate */
     	  Fixed[k]= 1;
     	  Dummy[k]= 0;
     	  modell[k].maintype= VTYPE;
     	  modell[k].subtype= VPDD;		/*	Product fixed dummy * varying dummy */
     	  ncovv++; /* Varying variables without age */
-    	  TvarV[ncovv]=Tvar[k];
-    	  TvarVind[ncovv]=k;
-    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  TvarV[ncovv]=Tvar[k];  /* TvarV[1]=Tvar[5]=5 because there is a V4 */
+    	  TvarVind[ncovv]=k;/* TvarVind[1]=5 */ 
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is a time varying quantitative covariate */
     	  Fixed[k]= 1;
     	  Dummy[k]= 1;
     	  modell[k].maintype= VTYPE;
@@ -11069,16 +11254,16 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
     	  TvarV[ncovv]=Tvar[k];
     	  TvarVind[ncovv]=k;
     	}
-      }else if(Tvard[k1][1] <=ncovcol+nqv){
-    	if(Tvard[k1][2] <=ncovcol){
-    	  Fixed[k]= 0;  /* or 2 ?*/
+      }else if(Tvard[k1][1] <=ncovcol+nqv){ /* Vn is fixed quanti  */
+    	if(Tvard[k1][2] <=ncovcol){ /* Vm is fixed dummy */
+    	  Fixed[k]= 0;  /*  Fixed product */
     	  Dummy[k]= 1;
     	  modell[k].maintype= FTYPE;
     	  modell[k].subtype= FPDQ;		/*	Product fixed quantitative * fixed dummy */
     	  ncovf++; /* Fixed variables without age */
     	  TvarF[ncovf]=Tvar[k];
     	  TvarFind[ncovf]=k;
-    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is time varying */
     	  Fixed[k]= 1;
     	  Dummy[k]= 1;
     	  modell[k].maintype= VTYPE;
@@ -11086,7 +11271,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
     	  ncovv++; /* Varying variables without age */
     	  TvarV[ncovv]=Tvar[k];
     	  TvarVind[ncovv]=k;
-    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is time varying quanti */
     	  Fixed[k]= 1;
     	  Dummy[k]= 1;
     	  modell[k].maintype= VTYPE;
@@ -11098,7 +11283,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
     	  TvarV[ncovv]=Tvar[k];
     	  TvarVind[ncovv]=k;
     	}
-      }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
+      }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ /* Vn is time varying dummy */
     	if(Tvard[k1][2] <=ncovcol){
     	  Fixed[k]= 1;
     	  Dummy[k]= 1;
@@ -11132,7 +11317,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
     	  TvarV[ncovv]=Tvar[k];
     	  TvarVind[ncovv]=k;
     	}
-      }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
+      }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ /* Vn is time varying quanti */
     	if(Tvard[k1][2] <=ncovcol){
     	  Fixed[k]= 1;
     	  Dummy[k]= 1;
@@ -11185,16 +11370,16 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       if((Typevar[k1]==Typevar[k2]) && (Fixed[k1]==Fixed[k2]) && (Dummy[k1]==Dummy[k2] )){
 	if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
 	  if(Tvar[k1]==Tvar[k2]){
-	    printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]);
-	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog);
+	    printf("Error duplication in the model=1+age+%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]);
+	    fprintf(ficlog,"Error duplication in the model=1+age+%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog);
 	    return(1);
 	  }
 	}else if (Typevar[k1] ==2){
 	  k3=Tposprod[k1];
 	  k4=Tposprod[k2];
 	  if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
-	    printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
-	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+	    printf("Error duplication in the model=1+age+%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+	    fprintf(ficlog,"Error duplication in the model=1+age+%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
 	    return(1);
 	  }
 	}
@@ -11553,6 +11738,7 @@ int prevalence_limit(double *p, double *
   /* for(k=1; k<=i1;k++){ /\* For each combination k of dummy covariates in the model *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       k=TKresult[nres];
+      if(TKresult[nres]==0) k=1; /* To be checked for noresult */
       /* if(i1 != 1 && TKresult[nres]!= k) /\* We found the combination k corresponding to the resultline value of dummies *\/ */
       /* 	continue; */
 
@@ -11658,23 +11844,30 @@ int back_prevalence_limit(double *p, dou
   if (cptcovn < 1){i1=1;}
   
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-     if(i1 != 1 && TKresult[nres]!= k)
-	continue;
-     /*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*/
+    /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
+      k=TKresult[nres];
+      if(TKresult[nres]==0) k=1; /* To be checked for noresult */
+     /* if(i1 != 1 && TKresult[nres]!= k) */
+     /* 	continue; */
+     /* /\*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*\/ */
       fprintf(ficresplb,"#******");
       printf("#******");
       fprintf(ficlog,"#******");
-      for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-	fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
-      }
-      for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
-	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+      for(j=1;j<=cptcovs ;j++) {/**< cptcovs number of SIMPLE covariates in the model or resultline V2+V1 =2 (dummy or quantit or time varying) */
+	printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+	fprintf(ficresplb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+	fprintf(ficlog," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
       }
+      /* for(j=1;j<=cptcoveff ;j++) {/\* all covariates *\/ */
+      /* 	fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+      /* 	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+      /* 	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+      /* } */
+      /* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */
+      /* 	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
+      /* 	fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
+      /* 	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
+      /* } */
       fprintf(ficresplb,"******\n");
       printf("******\n");
       fprintf(ficlog,"******\n");
@@ -11686,8 +11879,8 @@ int back_prevalence_limit(double *p, dou
       }
     
       fprintf(ficresplb,"#Age ");
-      for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+      for(j=1;j<=cptcovs;j++) {
+	fprintf(ficresplb,"V%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
       }
       for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
       fprintf(ficresplb,"Total Years_to_converge\n");
@@ -11710,8 +11903,8 @@ int back_prevalence_limit(double *p, dou
           /* exit(1); */
 	}
 	fprintf(ficresplb,"%.0f ",age );
-	for(j=1;j<=cptcoveff;j++)
-	  fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+	for(j=1;j<=cptcovs;j++)
+	  fprintf(ficresplb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
 	tot=0.;
 	for(i=1; i<=nlstate;i++){
 	  tot +=  bprlim[i][i];
@@ -11721,7 +11914,7 @@ int back_prevalence_limit(double *p, dou
       } /* Age */
       /* was end of cptcod */
       /*fprintf(ficresplb,"\n");*/ /* Seems to be necessary for gnuplot only if two result lines and no covariate. */
-    } /* end of any combination */
+    /* } /\* end of any combination *\/ */
   } /* end of nres */  
   /* hBijx(p, bage, fage); */
   /* fclose(ficrespijb); */
@@ -11765,6 +11958,7 @@ int hPijx(double *p, int bage, int fage)
   /*  	k=k+1;  */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     k=TKresult[nres];
+    if(TKresult[nres]==0) k=1; /* To be checked for noresult */
     /* for(k=1; k<=i1;k++){ */
     /* if(i1 != 1 && TKresult[nres]!= k) */
     /* 	continue; */
@@ -11846,12 +12040,13 @@ int hPijx(double *p, int bage, int fage)
   /*  	k=k+1;  */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     k=TKresult[nres];
+    if(TKresult[nres]==0) k=1; /* To be checked for noresult */
     /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
     /*    if(i1 != 1 && TKresult[nres]!= k) */
     /* 	continue; */
     fprintf(ficrespijb,"\n#****** ");
     for(j=1;j<=cptcovs;j++){
-      fprintf(ficrespij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+      fprintf(ficrespijb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
       /* for(j=1;j<=cptcoveff;j++) */
       /* 	fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
       /* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */
@@ -12291,7 +12486,7 @@ int main(int argc, char *argv[])
 	strcpy(model,modeltemp); 
       }
     }
-    /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
+    /* printf(" model=1+age%s modeltemp= %s, model=1+age+%s\n",model, modeltemp, model);fflush(stdout); */
     printf("model=1+age+%s\n",model);fflush(stdout);
     fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout);
     fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout);
@@ -12333,6 +12528,7 @@ int main(int argc, char *argv[])
   if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs);  /**< Fixed quantitative covariate */
   if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs);  /**< Time varying quantitative covariate */
   if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /**< Time varying covariate (dummy and quantitative)*/
+  /* if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);  /\**< Might be better *\/ */
   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
@@ -12613,6 +12809,8 @@ Please run with mle=-1 to get a correct
   TvarVDind=ivector(1,NCOVMAX); /*  */
   TvarVQ=ivector(1,NCOVMAX); /*  */
   TvarVQind=ivector(1,NCOVMAX); /*  */
+  TvarVV=ivector(1,NCOVMAX); /*  */
+  TvarVVind=ivector(1,NCOVMAX); /*  */
 
   Tvalsel=vector(1,NCOVMAX); /*  */
   Tvarsel=ivector(1,NCOVMAX); /*  */
@@ -12722,8 +12920,8 @@ Please run with mle=-1 to get a correct
   }
   
   ncovcombmax=pow(2,cptcoveff);
-  invalidvarcomb=ivector(1, ncovcombmax); 
-  for(i=1;i<ncovcombmax;i++)
+  invalidvarcomb=ivector(0, ncovcombmax); 
+  for(i=0;i<ncovcombmax;i++)
     invalidvarcomb[i]=0;
   
   /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
@@ -14173,12 +14371,15 @@ Please run with mle=-1 to get a correct
   free_ivector(TvarVDind,1,NCOVMAX);
   free_ivector(TvarVQ,1,NCOVMAX);
   free_ivector(TvarVQind,1,NCOVMAX);
+  free_ivector(TvarVV,1,NCOVMAX);
+  free_ivector(TvarVVind,1,NCOVMAX);
+  
   free_ivector(Tvarsel,1,NCOVMAX);
   free_vector(Tvalsel,1,NCOVMAX);
   free_ivector(Tposprod,1,NCOVMAX);
   free_ivector(Tprod,1,NCOVMAX);
   free_ivector(Tvaraff,1,NCOVMAX);
-  free_ivector(invalidvarcomb,1,ncovcombmax);
+  free_ivector(invalidvarcomb,0,ncovcombmax);
   free_ivector(Tage,1,NCOVMAX);
   free_ivector(Tmodelind,1,NCOVMAX);
   free_ivector(TmodelInvind,1,NCOVMAX);