--- imach096d/src/imach.c	2001/05/02 17:42:45	1.5
+++ imach096d/src/imach.c	2002/02/20 17:05:44	1.14
@@ -8,7 +8,7 @@
   Health expectancies are computed from the transistions observed between 
   waves and are computed for each degree of severity of disability (number
   of life states). More degrees you consider, more time is necessary to
-  reach the Maximum Likelihood of the parameters involved in the model.
+  reach the Maximum Likelihood of the parameters involved in the model. 
   The simplest model is the multinomial logistic model where pij is
   the probabibility to be observed in state j at the second wave conditional
   to be observed in state i at the first wave. Therefore the model is:
@@ -22,7 +22,7 @@
   delay between waves is not identical for each individual, or if some 
   individual missed an interview, the information is not rounded or lost, but
   taken into account using an interpolation or extrapolation.
-  hPijx is the probability to be
+  hPijx is the probability to be 
   observed in state i at age x+h conditional to the observed state i at age 
   x. The delay 'h' can be split into an exact number (nh*stepm) of 
   unobserved intermediate  states. This elementary transition (by month or 
@@ -54,8 +54,6 @@
 #define	GLOCK_ERROR_NOPATH		-1	/* empty path */
 #define	GLOCK_ERROR_GETCWD		-2	/* cannot get cwd */
 
-
-
 #define MAXPARM 30 /* Maximum number of parameters for the optimization */
 #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
 
@@ -70,22 +68,24 @@
 
 
 int nvar;
-static int cptcov;
-int cptcovn;
+int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 int npar=NPARMAX;
 int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int popbased=0, fprev,lprev;
 
 int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */
+int jmin, jmax; /* min, max spacing between 2 waves */
 int mle, weightopt;
 int **mw; /* mw[mi][i] is number of the mi wave for this individual */
 int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
+double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
-FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest;
-FILE *ficgp, *fichtm;
+FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;
+FILE *ficgp, *fichtm,*ficresprob;
 FILE *ficreseij;
   char filerese[FILENAMELENGTH];
  FILE  *ficresvij;
@@ -93,9 +93,6 @@ FILE *ficreseij;
  FILE  *ficresvpl;
   char fileresvpl[FILENAMELENGTH];
 
-
-
-
 #define NR_END 1
 #define FREE_ARG char*
 #define FTOL 1.0e-10
@@ -129,19 +126,19 @@ int stepm;
 /* Stepm, step in month: minimum step interpolation*/
 
 int m,nb;
-int *num, firstpass=0, lastpass=4,*cod, *ncodemax;
+int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
-double **pmmij;
+double **pmmij, ***probs, ***mobaverage;
 
 double *weight;
 int **s; /* Status */
 double *agedc, **covar, idx;
-int **nbcode, *Tcode, *Tvar, **codtab;
+int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
 
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */
 double ftolhess; /* Tolerance for computing hessian */
 
-
+/**************** split *************************/
 static	int split( char *path, char *dirc, char *name )
 {
    char	*s;				/* pointer */
@@ -205,7 +202,7 @@ int nbocc(char *s, char occ)
 
 void cutv(char *u,char *v, char*t, char occ)
 {
-  int i,lg,j,p;
+  int i,lg,j,p=0;
   i=0;
   for(j=0; j<=strlen(t)-1; j++) {
     if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
@@ -214,8 +211,8 @@ void cutv(char *u,char *v, char*t, char
   lg=strlen(t);
   for(j=0; j<p; j++) {
     (u[j] = t[j]);
-    u[p]='\0';
   }
+     u[p]='\0';
 
    for(j=0; j<= lg; j++) {
     if (j>=(p+1))(v[j-p-1] = t[j]);
@@ -652,15 +649,27 @@ double **prevalim(double **prlim, int nl
     for (j=1;j<=nlstate+ndeath;j++){
       oldm[ii][j]=(ii==j ? 1.0 : 0.0);
     }
-  /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+
+   cov[1]=1.;
+ 
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;
     /* Covariates have to be included here again */
-    cov[1]=1.;
-    cov[2]=agefin;
-    if (cptcovn>0){
-      for (k=1; k<=cptcovn;k++) {cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];/*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/}
-    }
+     cov[2]=agefin;
+  
+      for (k=1; k<=cptcovn;k++) {
+	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+	/*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/
+      }
+      for (k=1; k<=cptcovage;k++)
+	cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+      for (k=1; k<=cptcovprod;k++)
+	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+
+      /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
+      /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
+
     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
 
     savm=oldm;
@@ -685,7 +694,7 @@ double **prevalim(double **prlim, int nl
   }
 }
 
-/*************** transition probabilities **********/ 
+/*************** transition probabilities ***************/ 
 
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {
@@ -708,9 +717,11 @@ double **pmij(double **ps, double *cov,
 	s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
 	/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }
-      ps[i][j]=s2;
+      ps[i][j]=(s2);
     }
   }
+    /*ps[3][2]=1;*/
+
   for(i=1; i<= nlstate; i++){
      s1=0;
     for(j=1; j<i; j++)
@@ -732,6 +743,7 @@ double **pmij(double **ps, double *cov,
     }
   }
 
+
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){
      printf("%lf ",ps[ii][jj]);
@@ -749,7 +761,7 @@ double **pmij(double **ps, double *cov,
 
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
 {
-  /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
+  /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
   /* in, b, out are matrice of pointers which should have been initialized 
      before: only the contents of out is modified. The function returns
@@ -795,9 +807,13 @@ double ***hpxij(double ***po, int nhstep
       /* Covariates have to be included here again */
       cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
-      if (cptcovn>0){
-      for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];
-    }
+      for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+      for (k=1; k<=cptcovage;k++)
+	cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+      for (k=1; k<=cptcovprod;k++)
+	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[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]);*/
       out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, 
@@ -819,7 +835,7 @@ double ***hpxij(double ***po, int nhstep
 /*************** log-likelihood *************/
 double func( double *x)
 {
-  int i, ii, j, k, mi, d;
+  int i, ii, j, k, mi, d, kk;
   double l, ll[NLSTATEMAX], cov[NCOVMAX];
   double **out;
   double sw; /* Sum of weights */
@@ -829,29 +845,31 @@ double func( double *x)
   /* We are differentiating ll according to initial status */
   /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
   /*for(i=1;i<imx;i++) 
-printf(" %d\n",s[4][i]);
+    printf(" %d\n",s[4][i]);
   */
+  cov[1]=1.;
 
   for(k=1; k<=nlstate; k++) ll[k]=0.;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){
-       for(mi=1; mi<= wav[i]-1; mi++){
+    for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+    for(mi=1; mi<= wav[i]-1; mi++){
       for (ii=1;ii<=nlstate+ndeath;ii++)
 	for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
-            for(d=0; d<dh[mi][i]; d++){
+      for(d=0; d<dh[mi][i]; d++){
 	newm=savm;
-	  cov[1]=1.;
-	  cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
-	  if (cptcovn>0){
-	    for (k=1; k<=cptcovn;k++) cov[2+k]=covar[1+k-1][i];
-	    }
-	  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
-		       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
-	  savm=oldm;
-	  oldm=newm;
-
-
+	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+	for (kk=1; kk<=cptcovage;kk++) {
+	  cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+	}
+	
+	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+		     1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+	savm=oldm;
+	oldm=newm;
+	
+	
       } /* end mult */
-   
+      
       lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
       /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
       ipmx +=1;
@@ -899,7 +917,6 @@ void hesscov(double **matcov, double p[]
   void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;
 
-
   hess=matrix(1,npar,1,npar);
 
   printf("\nCalculation of the hessian matrix. Wait...\n");
@@ -907,14 +924,16 @@ void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/
+    /*printf(" %lf ",hess[i][i]);*/
   }
-
+  
   for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {
       if (j>i) { 
 	printf(".%d%d",i,j);fflush(stdout);
 	hess[i][j]=hessij(p,delti,i,j);
-	hess[j][i]=hess[i][j];
+	hess[j][i]=hess[i][j];    
+	/*printf(" %lf ",hess[i][j]);*/
       }
     }
   }
@@ -1018,7 +1037,7 @@ double hessii( double x[], double delta,
     }
   }
   delti[theta]=delts;
-  return res;
+  return res; 
   
 }
 
@@ -1142,7 +1161,7 @@ void  freqsummary(char fileres[], int ag
   char fileresp[FILENAMELENGTH];
 
   pp=vector(1,nlstate);
-
+ probs= ma3x(1,130 ,1,8, 1,8);
   strcpy(fileresp,"p");
   strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {
@@ -1152,13 +1171,14 @@ void  freqsummary(char fileres[], int ag
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;
 
-  j=cptcovn;
+  j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
 
   for(k1=1; k1<=j;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){
        j1++;
-
+       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
+	 scanf("%d", i);*/
         for (i=-1; i<=nlstate+ndeath; i++)  
 	 for (jk=-1; jk<=nlstate+ndeath; jk++)  
 	   for(m=agemin; m <= agemax+3; m++)
@@ -1167,11 +1187,12 @@ void  freqsummary(char fileres[], int ag
        for (i=1; i<=imx; i++) {
 	 bool=1;
 	 if  (cptcovn>0) {
-	   for (z1=1; z1<=cptcovn; z1++) 
-	     if (covar[Tvar[z1]][i]!= nbcode[Tvar[z1]][codtab[j1][z1]]) bool=0;
+	   for (z1=1; z1<=cptcoveff; z1++) 
+	     if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
+	       bool=0;
 	 }
 	  if (bool==1) {
-	   for(m=firstpass; m<=lastpass-1; m++){
+	   for(m=fprev; m<=lprev; m++){
 	     if(agev[m][i]==0) agev[m][i]=agemax+1;
 	     if(agev[m][i]==1) agev[m][i]=agemax+2;
 	     freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
@@ -1180,10 +1201,10 @@ void  freqsummary(char fileres[], int ag
 	 }
        }
         if  (cptcovn>0) {
-	 fprintf(ficresp, "\n#Variable"); 
-	 for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvar[z1],nbcode[Tvar[z1]][codtab[j1][z1]]);
-       }
-       fprintf(ficresp, "\n#");
+	 fprintf(ficresp, "\n#********** Variable "); 
+	 for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+       fprintf(ficresp, "**********\n#");
+	}
        for(i=1; i<=nlstate;i++) 
 	 fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
        fprintf(ficresp, "\n");
@@ -1195,7 +1216,7 @@ void  freqsummary(char fileres[], int ag
       printf("Age %d", i);
     for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-	pp[jk] += freq[jk][m][i];
+	pp[jk] += freq[jk][m][i]; 
     }
     for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pos=0; m <=0 ; m++)
@@ -1205,10 +1226,12 @@ void  freqsummary(char fileres[], int ag
       else
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     }
-    for(jk=1; jk <=nlstate ; jk++){
-      for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)
+
+     for(jk=1; jk <=nlstate ; jk++){
+      for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
 	pp[jk] += freq[jk][m][i];
-    }
+     }
+
     for(jk=1,pos=0; jk <=nlstate ; jk++)
       pos += pp[jk];
     for(jk=1; jk <=nlstate ; jk++){
@@ -1217,8 +1240,11 @@ void  freqsummary(char fileres[], int ag
       else
 	printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
       if( i <= (int) agemax){
-	if(pos>=1.e-5)
+	if(pos>=1.e-5){
 	  fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
+	  probs[i][jk][j1]= pp[jk]/pos;
+	  /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
+	}
       else
 	  fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
       }
@@ -1251,9 +1277,14 @@ void  concatwav(int wav[], int **dh, int
      */
 
   int i, mi, m;
-  int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
-float sum=0.;
+  /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
+     double sum=0., jmean=0.;*/
 
+  int j, k=0,jk, ju, jl;
+  double sum=0.;
+  jmin=1e+5;
+  jmax=-1;
+  jmean=0.;
   for(i=1; i<=imx; i++){
     mi=0;
     m=firstpass;
@@ -1283,14 +1314,22 @@ float sum=0.;
 	dh[mi][i]=1;
       else{
 	if (s[mw[mi+1][i]][i] > nlstate) {
+	  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 */
+	  if(j==0) j=1;  /* Survives at least one month after exam */
+	  k=k+1;
+	  if (j >= jmax) jmax=j;
+	  if (j <= jmin) jmin=j;
+	  sum=sum+j;
+	  /* if (j<10) printf("j=%d num=%d ",j,i); */
+	  }
 	}
 	else{
 	  j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
 	  k=k+1;
 	  if (j >= jmax) jmax=j;
 	  else if (j <= jmin)jmin=j;
+	  /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
 	  sum=sum+j;
 	}
 	jk= j/stepm;
@@ -1305,30 +1344,35 @@ float sum=0.;
       }
     }
   }
-  printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k);
-}
+  jmean=sum/k;
+  printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
+ }
 /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)
 {
-  int Ndum[80],ij, k, j, i;
+  int Ndum[20],ij=1, k, j, i;
   int cptcode=0;
-  for (k=0; k<79; k++) Ndum[k]=0;
+  cptcoveff=0; 
+ 
+  for (k=0; k<19; k++) Ndum[k]=0;
   for (k=1; k<=7; k++) ncodemax[k]=0;
-  
-  for (j=1; j<=cptcovn; j++) {
+
+  for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
     for (i=1; i<=imx; i++) {
       ij=(int)(covar[Tvar[j]][i]);
       Ndum[ij]++; 
+      /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
       if (ij > cptcode) cptcode=ij; 
     }
-    /*printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);*/
+
     for (i=0; i<=cptcode; i++) {
       if(Ndum[i]!=0) ncodemax[j]++;
     }
-  
     ij=1; 
+
+
     for (i=1; i<=ncodemax[j]; i++) {
-      for (k=0; k<=79; k++) {
+      for (k=0; k<=19; k++) {
 	if (Ndum[k] != 0) {
 	  nbcode[Tvar[j]][ij]=k; 
 	  ij++;
@@ -1336,9 +1380,25 @@ void tricode(int *Tvar, int **nbcode, in
 	if (ij > ncodemax[j]) break; 
       }  
     } 
-  }   
+  }  
 
-  }
+ for (k=0; k<19; k++) Ndum[k]=0;
+
+ for (i=1; i<=ncovmodel-2; i++) {
+      ij=Tvar[i];
+      Ndum[ij]++; 
+    }
+
+ ij=1;
+ for (i=1; i<=10; i++) {
+   if((Ndum[i]!=0) && (i<=ncov)){
+     Tvaraff[ij]=i; 
+     ij++;
+   }
+ }
+ 
+    cptcoveff=ij-1;
+}
 
 /*********** Health Expectancies ****************/
 
@@ -1399,7 +1459,7 @@ void varevsij(char fileres[], double ***
   double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;
   int k, cptcode;
-   double *xp;
+  double *xp;
   double **gp, **gm;
   double ***gradg, ***trgradg;
   double ***p3mat;
@@ -1435,6 +1495,12 @@ void varevsij(char fileres[], double ***
       }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+      if (popbased==1) {
+	for(i=1; i<=nlstate;i++)
+	  prlim[i][i]=probs[(int)age][i][ij];
+      }
+      
       for(j=1; j<= nlstate; j++){
 	for(h=0; h<=nhstepm; h++){
 	  for(i=1, gp[h][j]=0.;i<=nlstate;i++)
@@ -1446,12 +1512,19 @@ void varevsij(char fileres[], double ***
 	xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+      if (popbased==1) {
+	for(i=1; i<=nlstate;i++)
+	  prlim[i][i]=probs[(int)age][i][ij];
+      }
+
       for(j=1; j<= nlstate; j++){
 	for(h=0; h<=nhstepm; h++){
 	  for(i=1, gm[h][j]=0.;i<=nlstate;i++)
 	    gm[h][j] += prlim[i][i]*p3mat[i][j][h];
 	}
       }
+
       for(j=1; j<= nlstate; j++)
 	for(h=0; h<=nhstepm; h++){
 	  gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
@@ -1581,7 +1654,110 @@ void varprevlim(char fileres[], double *
 
 }
 
+/************ Variance of one-step probabilities  ******************/
+void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
+{
+  int i, j;
+  int k=0, cptcode;
+  double **dnewm,**doldm;
+  double *xp;
+  double *gp, *gm;
+  double **gradg, **trgradg;
+  double age,agelim, cov[NCOVMAX];
+  int theta;
+  char fileresprob[FILENAMELENGTH];
+
+  strcpy(fileresprob,"prob"); 
+  strcat(fileresprob,fileres);
+  if((ficresprob=fopen(fileresprob,"w"))==NULL) {
+    printf("Problem with resultfile: %s\n", fileresprob);
+  }
+  printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);
+  
+
+  xp=vector(1,npar);
+  dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+  doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
+  
+  cov[1]=1;
+  for (age=bage; age<=fage; age ++){ 
+    cov[2]=age;
+    gradg=matrix(1,npar,1,9);
+    trgradg=matrix(1,9,1,npar);
+    gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
+    gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
+    
+    for(theta=1; theta <=npar; theta++){
+      for(i=1; i<=npar; i++)
+	xp[i] = x[i] + (i==theta ?delti[theta]:0);
+     
+      pmij(pmmij,cov,ncovmodel,xp,nlstate);
+   
+      k=0;
+      for(i=1; i<= (nlstate+ndeath); 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]:0);
+    
+
+      pmij(pmmij,cov,ncovmodel,xp,nlstate);
+      k=0;
+      for(i=1; i<=(nlstate+ndeath); i++){
+	for(j=1; j<=(nlstate+ndeath);j++){
+	  k=k+1;
+	  gm[k]=pmmij[i][j];
+	}
+      }
+     
+       for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) 
+	   gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
+    }
+
+     for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
+      for(theta=1; theta <=npar; theta++)
+      trgradg[j][theta]=gradg[theta][j];
+ 
+     matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
+     matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
+
+     pmij(pmmij,cov,ncovmodel,x,nlstate);
+
+     k=0;
+     for(i=1; i<=(nlstate+ndeath); i++){
+       for(j=1; j<=(nlstate+ndeath);j++){
+	 k=k+1;
+	 gm[k]=pmmij[i][j];
+	}
+     }
+     
+     /*printf("\n%d ",(int)age);
+     for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
+       
+
+       printf("%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);
+
+  for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
+    if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
+if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
+  }
+
+    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);
+}
+ free_vector(xp,1,npar);
+fclose(ficresprob);
+ exit(0);
+}
 
 /***********************************************/
 /**************** Main Program *****************/
@@ -1591,7 +1767,7 @@ void varprevlim(char fileres[], double *
 int main()
 {
 
-  int i,j, k, n=MAXN,iter,m,size,cptcode, aaa, cptcod;
+  int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   double agedeb, agefin,hf;
   double agemin=1.e20, agemax=-1.e20;
 
@@ -1603,8 +1779,8 @@ int main()
   int *indx;
   char line[MAXLINE], linepar[MAXLINE];
   char title[MAXLINE];
-  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
-  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
+  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
+  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
   char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
@@ -1612,10 +1788,11 @@ int main()
   int sdeb, sfin; /* Status at beginning and end */
   int c,  h , cpt,l;
   int ju,jl, mi;
-  int i1,j1, k1,k2,k3,jk,aa,bb, stepsize;
-  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
-  
+  int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
+  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; 
+  int mobilav=0, fprevfore=1, lprevfore=1;
   int hstepm, nhstepm;
+
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
   double **prlim;
@@ -1628,9 +1805,12 @@ int main()
   double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;
-  char version[80]="Imach version 62c, May 1999, INED-EUROREVES ";
+  double kk1;
+
+  char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
+
   char z[1]="c", occ;
 #include <sys/time.h>
 #include <time.h>
@@ -1641,7 +1821,7 @@ int main()
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
 
 
-  printf("\nIMACH, Version 0.64a");
+  printf("\nIMACH, Version 0.64b");
   printf("\nEnter the parameter file name: ");
 
 #ifdef windows
@@ -1689,14 +1869,28 @@ split(pathtot, path,optionfile);
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
-
-  covar=matrix(1,NCOVMAX,1,n);    
-  if (strlen(model)<=1) cptcovn=0;
-  else {
-    j=0;
-    j=nbocc(model,'+');
-    cptcovn=j+1;
+while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
+  
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
   }
+  ungetc(c,ficpar);
+  
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);
+ 
+  covar=matrix(0,NCOVMAX,1,n); 
+  cptcovn=0; 
+  if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
 
   ncovmodel=2+cptcovn;
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
@@ -1727,7 +1921,8 @@ split(pathtot, path,optionfile);
       fprintf(ficparo,"\n");
     }
   
-  npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
+    npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
+
   p=param[1][1];
   
   /* Reads comments: lines beginning with '#' */
@@ -1788,7 +1983,6 @@ split(pathtot, path,optionfile);
   printf("\n");
 
 
-   if(mle==1){
     /*-------- data file ----------*/
     if((ficres =fopen(fileres,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileres);goto end;
@@ -1818,7 +2012,7 @@ split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);
 
-    i=1; 
+    i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {
       if ((i >= firstobs) && (i <=lastobs)) {
 	
@@ -1840,56 +2034,112 @@ split(pathtot, path,optionfile);
 	  cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
 	} 
 	num[i]=atol(stra);
-
-	/*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
+	
+	/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+	  printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
 
 	i=i+1;
       }
     } 
-
-    /*scanf("%d",i);*/
+    /* printf("ii=%d", ij);
+       scanf("%d",i);*/
   imx=i-1; /* Number of individuals */
 
+  /* for (i=1; i<=imx; i++){
+    if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
+    if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
+    if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
+    }
+    for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
+
   /* Calculation of the number of parameter from char model*/
-  Tvar=ivector(1,8);    
+  Tvar=ivector(1,15); 
+  Tprod=ivector(1,15); 
+  Tvaraff=ivector(1,15); 
+  Tvard=imatrix(1,15,1,2);
+  Tage=ivector(1,15);      
    
   if (strlen(model) >1){
-    j=0;
+    j=0, j1=0, k1=1, k2=1;
     j=nbocc(model,'+');
+    j1=nbocc(model,'*');
     cptcovn=j+1;
+    cptcovprod=j1;
+    
     
     strcpy(modelsav,model); 
-    if (j==0) {
-      cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb);
+    if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
+      printf("Error. Non available option model=%s ",model);
+      goto end;
     }
-    else {
-      for(i=j; i>=1;i--){
-	cutv(stra,strb,modelsav,'+');
-	if (strchr(strb,'*')) {
-	  cutv(strd,strc,strb,'*');
-	  cutv(strb,stre,strc,'V');Tvar[i+1]=ncov+1;
+    
+    for(i=(j+1); i>=1;i--){
+      cutv(stra,strb,modelsav,'+');
+      if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); 
+      /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+      /*scanf("%d",i);*/
+      if (strchr(strb,'*')) {
+	cutv(strd,strc,strb,'*');
+	if (strcmp(strc,"age")==0) {
+	  cptcovprod--;
+	  cutv(strb,stre,strd,'V');
+	  Tvar[i]=atoi(stre);
+	  cptcovage++;
+	    Tage[cptcovage]=i;
+	    /*printf("stre=%s ", stre);*/
+	}
+	else if (strcmp(strd,"age")==0) {
+	  cptcovprod--;
+	  cutv(strb,stre,strc,'V');
+	  Tvar[i]=atoi(stre);
+	  cptcovage++;
+	  Tage[cptcovage]=i;
+	}
+	else {
+	  cutv(strb,stre,strc,'V');
+	  Tvar[i]=ncov+k1;
 	  cutv(strb,strc,strd,'V'); 
+	  Tprod[k1]=i;
+	  Tvard[k1][1]=atoi(strc);
+	  Tvard[k1][2]=atoi(stre);
+	  Tvar[cptcovn+k2]=Tvard[k1][1];
+	  Tvar[cptcovn+k2+1]=Tvard[k1][2]; 
 	  for (k=1; k<=lastobs;k++) 
-	    covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
-	}
-	else {cutv(strd,strc,strb,'V');
-	Tvar[i+1]=atoi(strc);
+	    covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
+	  k1++;
+	  k2=k2+2;
 	}
-	strcpy(modelsav,stra);   
       }
-      cutv(strd,strc,stra,'V');
-      Tvar[1]=atoi(strc);
+      else {
+	/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+       /*  scanf("%d",i);*/
+      cutv(strd,strc,strb,'V');
+      Tvar[i]=atoi(strc);
+      }
+      strcpy(modelsav,stra);  
+      /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+	scanf("%d",i);*/
     }
-  }
-  /*printf("tvar=%d ",Tvar[1]);
-  scanf("%d ",i);*/ 
+}
+  
+  /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
+  printf("cptcovprod=%d ", cptcovprod);
+  scanf("%d ",i);*/
     fclose(fic);
 
+    /*  if(mle==1){*/
     if (weightopt != 1) { /* Maximisation without weights*/
       for(i=1;i<=n;i++) weight[i]=1.0;
     }
     /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);
+
+   for (i=1; i<=imx; i++) 
+     for(m=2; (m<= maxwav); m++)
+       if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
+	 anint[m][i]=9999;
+	 s[m][i]=-1;
+       }
     
     for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
@@ -1899,9 +2149,11 @@ split(pathtot, path,optionfile);
 	    if(agedc[i]>0)
 	      if(moisdc[i]!=99 && andc[i]!=9999)
 	      agev[m][i]=agedc[i];
-	    else{
+	    else {
+	      if (andc[i]!=9999){
 	      printf("Warning negative age at death: %d line:%d\n",num[i],i);
 	      agev[m][i]=-1;
+	      }
 	    }
 	  }
 	  else if(s[m][i] !=9){ /* Should no more exist */
@@ -1958,19 +2210,19 @@ printf("Total number of individuals= %d,
       concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
 
 
-Tcode=ivector(1,100);
-   nbcode=imatrix(1,nvar,1,8);  
-   ncodemax[1]=1;
-   if (cptcovn > 0) tricode(Tvar,nbcode,imx);
- 
+      Tcode=ivector(1,100);
+      nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
+      ncodemax[1]=1;
+      if (cptcovn > 0) tricode(Tvar,nbcode,imx);
+      
    codtab=imatrix(1,100,1,10);
    h=0;
-   m=pow(2,cptcovn);
+   m=pow(2,cptcoveff);
  
-   for(k=1;k<=cptcovn; k++){
+   for(k=1;k<=cptcoveff; k++){
      for(i=1; i <=(m/pow(2,k));i++){
        for(j=1; j <= ncodemax[k]; j++){
-	 for(cpt=1; cpt <=(m/pow(2,cptcovn+1-k)); cpt++){
+	 for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
 	   h++;
 	   if (h>m) h=1;codtab[h][k]=j;
 	 } 
@@ -1978,13 +2230,14 @@ Tcode=ivector(1,100);
      }
    } 
 
+
    /*for(i=1; i <=m ;i++){ 
      for(k=1; k <=cptcovn; k++){
-       printf("i=%d k=%d %d ",i,k,codtab[i][k]);
+       printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
      }
      printf("\n");
-   }*/
-   /*scanf("%d",i);*/
+   }
+   scanf("%d",i);*/
     
    /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */
@@ -1995,13 +2248,14 @@ Tcode=ivector(1,100);
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
-    
+     
     /* For Powell, parameters are in a vector p[] starting at p[1]
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */
-    
-    mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
 
+    if(mle==1){
+    mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
+    }
     
     /*--------- results files --------------*/
     fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
@@ -2025,10 +2279,11 @@ Tcode=ivector(1,100);
 	 }
      }
    }
-
+ if(mle==1){
     /* Computing hessian and covariance matrix */
     ftolhess=ftol; /* Usually correct */
     hesscov(matcov, p, npar, delti, ftolhess, func);
+ }
     fprintf(ficres,"# Scales\n");
     printf("# Scales\n");
      for(i=1,jk=1; i <=nlstate; i++){
@@ -2083,6 +2338,8 @@ Tcode=ivector(1,100);
 
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
+
+   
 /*------------ gnuplot -------------*/
 chdir(pathcd);
   if((ficgp=fopen("graph.plt","w"))==NULL) {
@@ -2091,7 +2348,7 @@ chdir(pathcd);
 #ifdef windows
   fprintf(ficgp,"cd \"%s\" \n",pathc);
 #endif
-m=pow(2,cptcovn);
+m=pow(2,cptcoveff);
   
  /* 1eme*/
   for (cpt=1; cpt<= nlstate ; cpt ++) {
@@ -2187,7 +2444,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);   
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
     } 
-  }
+  }  
 
   /* proba elementaires */
    for(i=1,jk=1; i <=nlstate; i++){
@@ -2212,26 +2469,38 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
      for(k=1; k<=(nlstate+ndeath); k++) {
        if (k != k2){
     	fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-
-	for(j=3; j <=ncovmodel; j++) 
-	  fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
-	fprintf(ficgp,")/(1");
-
-	for(k1=1; k1 <=nlstate+1; k1=k1+2){   
-	    fprintf(ficgp,"+exp(p%d+p%d*x",k1+k3-1,k1+k3);
-
-	    for(j=3; j <=ncovmodel; j++)
-	      fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
-	    fprintf(ficgp,")");
+ij=1;
+	for(j=3; j <=ncovmodel; j++) {
+	  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+	    fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+	    ij++;
+	  }
+	  else
+	  fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+	}
+	  fprintf(ficgp,")/(1");
+	
+	for(k1=1; k1 <=nlstate; k1++){   
+	  fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ij=1;
+	  for(j=3; j <=ncovmodel; j++){
+	  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+	    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+	    ij++;
+	  }
+	  else
+	    fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+	  }
+	  fprintf(ficgp,")");
 	}
 	fprintf(ficgp,") t \"p%d%d\" ", k2,k);
 	if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
-    i=i+ncovmodel;
+	i=i+ncovmodel;
        }
      }
    }
-  fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
-   }
+   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
+  }
    
   fclose(ficgp);
    
@@ -2250,7 +2519,7 @@ chdir(path);
     /*free_matrix(covar,1,NCOVMAX,1,n);*/
     fclose(ficparo);
     fclose(ficres);
-   }
+    /*  }*/
    
    /*________fin mle=1_________*/
    
@@ -2271,11 +2540,18 @@ chdir(path);
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 /*--------- index.htm --------*/
 
-  if((fichtm=fopen("index.htm","w"))==NULL)    {
-    printf("Problem with index.htm \n");goto end;
+  strcpy(optionfilehtm,optionfile);
+  strcat(optionfilehtm,".htm");
+  if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
+    printf("Problem with %s \n",optionfilehtm);goto end;
   }
 
- fprintf(fichtm,"<body><ul> Imach, Version 0.64a<hr> <li>Outputs files<br><br>\n
+ fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\"> 
+Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
+Total number of observations=%d <br>
+Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
+<hr  size=\"2\" color=\"#EC5E5E\"> 
+<li>Outputs files<br><br>\n
         - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
 - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
         - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
@@ -2284,11 +2560,13 @@ chdir(path);
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
-        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
+        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
+        - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
+<br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
 
- fprintf(fichtm," <li>Graphs</li>\n<p>");
+ fprintf(fichtm," <li>Graphs</li><p>");
 
- m=cptcovn;
+ m=cptcoveff;
  if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
  j1=0;
@@ -2296,10 +2574,10 @@ chdir(path);
    for(i1=1; i1<=ncodemax[k1];i1++){
        j1++;
        if (cptcovn > 0) {
-	 fprintf(fichtm,"<hr>************ Results for covariates");
-	 for (cpt=1; cpt<=cptcovn;cpt++) 
-	   fprintf(fichtm," V%d=%d ",Tvar[cpt],nbcode[Tvar[cpt]][codtab[j1][cpt]]);
-	 fprintf(fichtm," ************\n<hr>");
+	 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]][codtab[j1][cpt]]);
+	 fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        }
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);     
@@ -2347,16 +2625,16 @@ fclose(fichtm);
   agebase=agemin;
   agelim=agemax;
   ftolpl=1.e-10;
-  i1=cptcovn;
+  i1=cptcoveff;
   if (cptcovn < 1){i1=1;}
 
   for(cptcov=1;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
 	k=k+1;
 	/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
-	fprintf(ficrespl,"\n#****** ");
-	for(j=1;j<=cptcovn;j++) 
-	  fprintf(ficrespl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+	fprintf(ficrespl,"\n#******");
+	for(j=1;j<=cptcoveff;j++) 
+	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
 	fprintf(ficrespl,"******\n");
 	
 	for (age=agebase; age<=agelim; age++){
@@ -2369,6 +2647,7 @@ fclose(fichtm);
       }
     }
   fclose(ficrespl);
+
   /*------------- h Pij x at various ages ------------*/
   
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);
@@ -2378,7 +2657,7 @@ fclose(fichtm);
   printf("Computing pij: result on file '%s' \n", filerespij);
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;
-  if (stepm<=24) stepsize=2;
+  /*if (stepm<=24) stepsize=2;*/
 
   agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */
@@ -2389,8 +2668,8 @@ fclose(fichtm);
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;
 	fprintf(ficrespij,"\n#****** ");
-	for(j=1;j<=cptcovn;j++) 
-	  fprintf(ficrespij,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+	for(j=1;j<=cptcoveff;j++) 
+	  fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
 	fprintf(ficrespij,"******\n");
 	
 	for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -2417,8 +2696,113 @@ fclose(fichtm);
     }
   }
 
+  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
+
   fclose(ficrespij);
 
+  /*---------- Forecasting ------------------*/
+
+  strcpy(fileresf,"f"); 
+  strcat(fileresf,fileres);
+  if((ficresf=fopen(fileresf,"w"))==NULL) {
+    printf("Problem with forecast resultfile: %s\n", fileresf);goto end;
+  }
+  printf("Computing forecasting: result on file '%s' \n", fileresf);
+
+  /* Mobile average */
+
+  /* for (agedeb=bage; agedeb<=fage; agedeb++)
+    for (i=1; i<=nlstate;i++)
+      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)
+      printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/
+
+  if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+  if (mobilav==1) {
+    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+    for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
+      for (i=1; i<=nlstate;i++)
+	for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
+	  mobaverage[(int)agedeb][i][cptcod]=0.;
+    
+    for (agedeb=bage+4; agedeb<=fage; agedeb++){
+      for (i=1; i<=nlstate;i++){
+	for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+	  for (cpt=0;cpt<=4;cpt++){
+	    mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
+	  }
+	  mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
+	}
+      }
+    }   
+  }
+
+  stepsize=(int) (stepm+YEARM-1)/YEARM;
+  if (stepm<=24) stepsize=2;
+
+  agelim=AGESUP;
+  hstepm=stepsize*YEARM; /* Every year of age */
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
+  hstepm=12;
+   k=0;
+  for(cptcov=1;cptcov<=i1;cptcov++){
+    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+      k=k+1;
+      fprintf(ficresf,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+      }
+      
+      fprintf(ficresf,"******\n");
+
+      fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");
+      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
+
+      for (agedeb=fage; agedeb>=bage; agedeb--){ 
+	fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
+       if (mobilav==1) {
+	for(j=1; j<=nlstate;j++) 
+	  fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);
+	}
+	else {
+	  for(j=1; j<=nlstate;j++) 
+	  fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]);
+	}
+       
+      for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");
+      }
+      for (cpt=1; cpt<=NCOVMAX;cpt++)   
+      for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+	
+	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
+	nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+	/*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/
+
+	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	oldm=oldms;savm=savms;
+	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
+		
+	for (h=0; h<=nhstepm; h++){
+	
+	  if (h*hstepm/YEARM*stepm==cpt)
+ fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);
+	  
+	  for(j=1; j<=nlstate+ndeath;j++) {
+	    kk1=0.;
+	    for(i=1; i<=nlstate;i++) {	      
+	      if (mobilav==1)
+	      kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
+	      else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];
+	    }    
+	    if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);
+	  }
+	}
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+      }
+    }
+  }
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  fclose(ficresf);
   /*---------- Health expectancies and variances ------------*/
 
   strcpy(filerest,"t");
@@ -2448,17 +2832,17 @@ fclose(fichtm);
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;
       fprintf(ficrest,"\n#****** ");
-      for(j=1;j<=cptcovn;j++) 
-	fprintf(ficrest,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+      for(j=1;j<=cptcoveff;j++) 
+	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficrest,"******\n");
 
       fprintf(ficreseij,"\n#****** ");
-      for(j=1;j<=cptcovn;j++) 
+      for(j=1;j<=cptcoveff;j++) 
 	fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
       fprintf(ficreseij,"******\n");
 
       fprintf(ficresvij,"\n#****** ");
-      for(j=1;j<=cptcovn;j++) 
+      for(j=1;j<=cptcoveff;j++) 
 	fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
       fprintf(ficresvij,"******\n");
 
@@ -2478,6 +2862,11 @@ fclose(fichtm);
       epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){
 	prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+	if (popbased==1) {
+	  for(i=1; i<=nlstate;i++)
+	    prlim[i][i]=probs[(int)age][i][k];
+	}
+	
 	fprintf(ficrest," %.0f",age);
 	for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
 	  for(i=1, epj[j]=0.;i <=nlstate;i++) {
@@ -2497,6 +2886,9 @@ fclose(fichtm);
     }
   }
 	
+       
+
+
  fclose(ficreseij);
  fclose(ficresvij);
   fclose(ficrest);
@@ -2519,8 +2911,8 @@ strcpy(fileresvpl,"vpl");
    for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
      k=k+1;
      fprintf(ficresvpl,"\n#****** ");
-     for(j=1;j<=cptcovn;j++) 
-       fprintf(ficresvpl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+     for(j=1;j<=cptcoveff;j++) 
+       fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
      fprintf(ficresvpl,"******\n");
      
      varpl=matrix(1,nlstate,(int) bage, (int) fage);
@@ -2542,7 +2934,7 @@ strcpy(fileresvpl,"vpl");
   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_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);
   
@@ -2555,11 +2947,13 @@ strcpy(fileresvpl,"vpl");
   /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/
 
+
  end:
 #ifdef windows
  chdir(pathcd);
 #endif 
- system("wgnuplot graph.plt");
+ 
+ system("..\\gp37mgw\\wgnuplot graph.plt");
 
 #ifdef windows
   while (z[0] != 'q') {
@@ -2569,7 +2963,7 @@ strcpy(fileresvpl,"vpl");
     if (z[0] == 'c') system("./imach");
     else if (z[0] == 'e') {
       chdir(path);
-      system("index.htm");
+      system(optionfilehtm);
     }
     else if (z[0] == 'q') exit(0);
   }