--- imach096d/src/imach.c	2002/04/03 12:19:36	1.38
+++ imach096d/src/imach.c	2003/06/13 07:45:28	1.41.2.2
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.38 2002/04/03 12:19:36 lievre Exp $
+/* $Id: imach.c,v 1.41.2.2 2003/06/13 07:45:28 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -14,7 +14,7 @@
   model. More health states you consider, more time is necessary to 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
+  probability 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: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
   'age' is age and 'sex' is a covariate. If you want to have a more
@@ -60,7 +60,8 @@
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80
 /*#define DEBUG*/
-#define windows
+
+/*#define windows*/
 #define	GLOCK_ERROR_NOPATH		-1	/* empty path */
 #define	GLOCK_ERROR_GETCWD		-2	/* cannot get cwd */
 
@@ -869,6 +870,7 @@ double func( double *x)
   double **out;
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
+  int s1, s2;
   long ipmx;
   /*extern weight */
   /* We are differentiating ll according to initial status */
@@ -883,7 +885,10 @@ double func( double *x)
     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 (j=1;j<=nlstate+ndeath;j++){
+	  oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	  savm[ii][j]=(ii==j ? 1.0 : 0.0);
+	}
       for(d=0; d<dh[mi][i]; d++){
 	newm=savm;
 	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
@@ -899,17 +904,47 @@ double func( double *x)
 	
       } /* 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]]);*/
+      s1=s[mw[mi][i]][i];
+      s2=s[mw[mi+1][i]][i];
+      if( s2 > nlstate){ 
+	/* i.e. if s2 is a death state and if the date of death is known then the contribution
+	   to the likelihood is the probability to die between last step unit time and current 
+	   step unit time, which is also the differences between probability to die before dh 
+	   and probability to die before dh-stepm . 
+	   In version up to 0.92 likelihood was computed
+	   as if date of death was unknown. Death was treated as any other
+	   health state: the date of the interview describes the actual state
+	   and not the date of a change in health state. The former idea was
+	   to consider that at each interview the state was recorded
+	   (healthy, disable or death) and IMaCh was corrected; but when we
+	   introduced the exact date of death then we should have modified
+	   the contribution of an exact death to the likelihood. This new
+	   contribution is smaller and very dependent of the step unit
+	   stepm. It is no more the probability to die between last interview
+	   and month of death but the probability to survive from last
+	   interview up to one month before death multiplied by the
+	   probability to die within a month. Thanks to Chris
+	   Jackson for correcting this bug.  Former versions increased
+	   mortality artificially. The bad side is that we add another loop
+	   which slows down the processing. The difference can be up to 10%
+	   lower mortality.
+	*/
+	lli=log(out[s1][s2] - savm[s1][s2]);
+      }else{
+	lli=log(out[s1][s2]); /* or	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;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d lli=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],lli,weight[i],out[s1][s2],savm[s1][s2]);*/
     } /* end of wave */
   } /* end of individual */
 
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+  /*exit(0);*/
   return -l;
 }
 
@@ -1352,8 +1387,11 @@ void prevalence(int agemin, float agemax
 	    if ((k2>=dateprev1) && (k2<=dateprev2)) {
 	      if(agev[m][i]==0) agev[m][i]=agemax+1;
 	      if(agev[m][i]==1) agev[m][i]=agemax+2;
-	      if (m<lastpass) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
-	      /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  */ 
+	      if (m<lastpass) 
+	      	if (calagedate>0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
+	      else
+	       freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
+	       freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; 
 	    }
 	  }
 	}
@@ -1386,7 +1424,7 @@ void prevalence(int agemin, float agemax
 	}
     }
   }
-  
+
   
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
@@ -1503,7 +1541,7 @@ void tricode(int *Tvar, int **nbcode, in
       for (k=0; k<=19; k++) {
 	if (Ndum[k] != 0) {
 	  nbcode[Tvar[j]][ij]=k; 
-	  /*	 printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/
+	  
 	  ij++;
 	}
 	if (ij > ncodemax[j]) break; 
@@ -1531,18 +1569,29 @@ void tricode(int *Tvar, int **nbcode, in
 
 /*********** Health Expectancies ****************/
 
-void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm)
+void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )
+
 {
   /* Health expectancies */
-  int i, j, nhstepm, hstepm, h, nstepm;
+  int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
   double age, agelim, hf;
-  double ***p3mat;
+  double ***p3mat,***varhe;
+  double **dnewm,**doldm;
+  double *xp;
+  double **gp, **gm;
+  double ***gradg, ***trgradg;
+  int theta;
+
+  varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);
+  xp=vector(1,npar);
+  dnewm=matrix(1,nlstate*2,1,npar);
+  doldm=matrix(1,nlstate*2,1,nlstate*2);
   
   fprintf(ficreseij,"# Health expectancies\n");
   fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)
-      fprintf(ficreseij," %1d-%1d",i,j);
+      fprintf(ficreseij," %1d-%1d (SE)",i,j);
   fprintf(ficreseij,"\n");
 
   if(estepm < stepm){
@@ -1582,24 +1631,113 @@ void evsij(char fileres[], double ***eij
     /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+    gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);
+    gp=matrix(0,nhstepm,1,nlstate*2);
+    gm=matrix(0,nhstepm,1,nlstate*2);
+
     /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
+ 
+
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
+
+    /* Computing Variances of health expectancies */
+
+     for(theta=1; theta <=npar; theta++){
+      for(i=1; i<=npar; i++){ 
+	xp[i] = x[i] + (i==theta ?delti[theta]:0);
+      }
+      hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
+  
+      cptj=0;
+      for(j=1; j<= nlstate; j++){
+	for(i=1; i<=nlstate; i++){
+	  cptj=cptj+1;
+	  for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){
+	    gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
+	  }
+	}
+      }
+     
+     
+      for(i=1; i<=npar; i++) 
+	xp[i] = x[i] - (i==theta ?delti[theta]:0);
+      hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
+      
+      cptj=0;
+      for(j=1; j<= nlstate; j++){
+	for(i=1;i<=nlstate;i++){
+	  cptj=cptj+1;
+	  for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
+	    gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
+	  }
+	}
+      }
+     
+   
+
+      for(j=1; j<= nlstate*2; j++)
+	for(h=0; h<=nhstepm-1; h++){
+	  gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
+	}
+
+     } 
+   
+/* End theta */
+
+     trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);
+
+     for(h=0; h<=nhstepm-1; h++)
+      for(j=1; j<=nlstate*2;j++)
+	for(theta=1; theta <=npar; theta++)
+	trgradg[h][j][theta]=gradg[h][theta][j];
+
+
+     for(i=1;i<=nlstate*2;i++)
+      for(j=1;j<=nlstate*2;j++)
+	varhe[i][j][(int)age] =0.;
+
+    for(h=0;h<=nhstepm-1;h++){
+      for(k=0;k<=nhstepm-1;k++){
+	matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
+	matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
+	for(i=1;i<=nlstate*2;i++)
+	  for(j=1;j<=nlstate*2;j++)
+	    varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
+      }
+    }
+
+     
+    /* Computing expectancies */
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)
 	for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
 	  eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
-	  /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
+	  
+/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
+
 	}
+
     fprintf(ficreseij,"%3.0f",age );
+    cptj=0;
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){
-	fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);
+	cptj++;
+	fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
       }
     fprintf(ficreseij,"\n");
+   
+    free_matrix(gm,0,nhstepm,1,nlstate*2);
+    free_matrix(gp,0,nhstepm,1,nlstate*2);
+    free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);
+    free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }
+  free_vector(xp,1,npar);
+  free_matrix(dnewm,1,nlstate*2,1,npar);
+  free_matrix(doldm,1,nlstate*2,1,nlstate*2);
+  free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);
 }
 
 /************ Variance ******************/
@@ -1822,9 +1960,9 @@ 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)
+void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
 {
-  int i, j;
+  int i, j, i1, k1, j1, z1;
   int k=0, cptcode;
   double **dnewm,**doldm;
   double *xp;
@@ -1839,91 +1977,118 @@ void varprob(char fileres[], double **ma
   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);
+  printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
   
+fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
+  fprintf(ficresprob,"# Age");
+  for(i=1; i<=nlstate;i++)
+    for(j=1; j<=(nlstate+ndeath);j++)
+      fprintf(ficresprob," p%1d-%1d (SE)",i,j);
+
+
+  fprintf(ficresprob,"\n");
+
 
   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));
+  j=cptcoveff;
+  if (cptcovn<1) {j=1;ncodemax[1]=1;}
+  j1=0;
+  for(k1=1; k1<=1;k1++){
+    for(i1=1; i1<=ncodemax[k1];i1++){ 
+    j1++;
+
+    if  (cptcovn>0) {
+      fprintf(ficresprob, "\n#********** Variable "); 
+      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+      fprintf(ficresprob, "**********\n#");
+    }
     
-    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 (age=bage; age<=fage; age ++){ 
+	cov[2]=age;
+	for (k=1; k<=cptcovn;k++) {
+	  cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
+	  
 	}
-      }
-
-      for(i=1; i<=npar; i++)
-	xp[i] = x[i] - (i==theta ?delti[theta]:0);
+	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]]];
+	
+	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));
     
-
-      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(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);
+	  for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) 
+	    gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
+	}
 
-     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(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]);
-  }
+	fprintf(ficresprob,"\n%d ",(int)age);
 
+	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
+	  fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(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);
-
+  }
+  free_vector(xp,1,npar);
+  fclose(ficresprob);
+  
 }
 
 /******************* Printing html file ***********/
@@ -1958,9 +2123,10 @@ Interval (in months) between two waves:
 
  fprintf(fichtm,"\n
  - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n
- - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
+  - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
+ - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n 
  - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n
- - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
+ - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
 
  if(popforecast==1) fprintf(fichtm,"\n
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
@@ -2027,12 +2193,7 @@ m=pow(2,cptcoveff);
   for (cpt=1; cpt<= nlstate ; cpt ++) {
    for (k1=1; k1<= m ; k1 ++) {
 
-#ifdef windows
-    fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
-#endif
-#ifdef unix
-fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
-#endif
+     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
 
 for (i=1; i<= nlstate ; i ++) {
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
@@ -2049,9 +2210,7 @@ for (i=1; i<= nlstate ; i ++) {
   else fprintf(ficgp," \%%*lf (\%%*lf)");
 }  
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
-#ifdef unix
-fprintf(ficgp,"\nset ter gif small size 400,300");
-#endif
+
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
    }
   }
@@ -2090,10 +2249,19 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
 
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
-      k=2+nlstate*(cpt-1);
+      k=2+nlstate*(2*cpt-2);
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
+      /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+
+*/
       for (i=1; i< nlstate ; i ++) {
-	fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);
+	fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);
+
       } 
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
     }
@@ -2493,7 +2661,7 @@ int main(int argc, char *argv[])
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; 
   int mobilav=0,popforecast=0;
   int hstepm, nhstepm;
-  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;
+  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
 
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
@@ -2511,7 +2679,7 @@ int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   
 
-  char version[80]="Imach version 0.8b, March 2002, INED-EUROREVES ";
+  char version[80]="Imach version 0.8a1, June 2003, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
 
@@ -2751,11 +2919,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     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++){
+   /*  for (i=1; i<=imx; i++){
      if (s[4][i]==9)  s[4][i]=-1; 
-     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]));}
-  */
+     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,15); 
@@ -3207,20 +3374,20 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 	    for(j=1; j<=nlstate+ndeath;j++)
 	      fprintf(ficrespij," %1d-%1d",i,j);
 	  fprintf(ficrespij,"\n");
-	  for (h=0; h<=nhstepm; h++){
+	   for (h=0; h<=nhstepm; h++){
 	    fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
 	    for(i=1; i<=nlstate;i++)
 	      for(j=1; j<=nlstate+ndeath;j++)
 		fprintf(ficrespij," %.5f", p3mat[i][j][h]);
 	    fprintf(ficrespij,"\n");
-	  }
+	     }
 	  free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
 	  fprintf(ficrespij,"\n");
 	}
     }
   }
 
-  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
+  varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
 
   fclose(ficrespij);
 
@@ -3229,9 +3396,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   if((stepm == 1) && (strcmp(model,".")==0)){
     prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
     if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
-    free_matrix(mint,1,maxwav,1,n);
-    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
-    free_vector(weight,1,n);}
+  } 
   else{
     erreur=108;
     printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
@@ -3261,6 +3426,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
   }
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
+  calagedate=-1;
+prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
 
   k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){
@@ -3283,7 +3450,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;
-      evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm);  
+      evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  
+ 
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;
        varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
@@ -3306,9 +3474,11 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 	for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
 	  for(i=1, epj[j]=0.;i <=nlstate;i++) {
 	    epj[j] += prlim[i][i]*eij[i][j][(int)age];
+	    /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
 	  }
 	  epj[nlstate+1] +=epj[j];
 	}
+
 	for(i=1, vepp=0.;i <=nlstate;i++)
 	  for(j=1;j <=nlstate;j++)
 	    vepp += vareij[i][j][(int)age];
@@ -3320,7 +3490,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){
       }
     }
   }
-
+free_matrix(mint,1,maxwav,1,n);
+    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
+    free_vector(weight,1,n);
   fclose(ficreseij);
   fclose(ficresvij);
   fclose(ficrest);
@@ -3382,9 +3554,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 
 
  end:
-#ifdef windows
   /* chdir(pathcd);*/
-#endif 
  /*system("wgnuplot graph.plt");*/
  /*system("../gp37mgw/wgnuplot graph.plt");*/
  /*system("cd ../gp37mgw");*/
@@ -3394,7 +3564,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
  strcat(plotcmd,optionfilegnuplot);
  system(plotcmd);
 
-#ifdef windows
+ /*#ifdef windows*/
   while (z[0] != 'q') {
     /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
@@ -3404,7 +3574,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }
-#endif 
+  /*#endif */
 }