--- imach096d/src/imach.c	2002/04/05 15:45:00	1.39
+++ imach096d/src/imach.c	2002/05/24 13:00:54	1.43
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.39 2002/04/05 15:45:00 lievre Exp $
+/* $Id: imach.c,v 1.43 2002/05/24 13:00:54 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -56,7 +56,7 @@
 #include <unistd.h>
 
 #define MAXLINE 256
-#define GNUPLOTPROGRAM "wgnuplot"
+#define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80
 /*#define DEBUG*/
@@ -1330,10 +1330,10 @@ void prevalence(int agemin, float agemax
   j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
- for(k1=1; k1<=j;k1++){
+  for(k1=1; k1<=j;k1++){
     for(i1=1; i1<=ncodemax[k1];i1++){
       j1++;
- 
+      
       for (i=-1; i<=nlstate+ndeath; i++)  
 	for (jk=-1; jk<=nlstate+ndeath; jk++)  
 	  for(m=agemin; m <= agemax+3; m++)
@@ -1352,41 +1352,46 @@ 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]; 
+	      }
 	    }
 	  }
 	}
       }
-	for(i=(int)agemin; i <= (int)agemax+3; i++){ 
-	  for(jk=1; jk <=nlstate ; jk++){
-	    for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-	      pp[jk] += freq[jk][m][i]; 
-	  }
-	  for(jk=1; jk <=nlstate ; jk++){
-	    for(m=-1, pos=0; m <=0 ; m++)
+      for(i=(int)agemin; i <= (int)agemax+3; i++){ 
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+	    pp[jk] += freq[jk][m][i]; 
+	}
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pos=0; m <=0 ; m++)
 	    pos += freq[jk][m][i];
 	}
 	
-	 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++){	   
-	   if( i <= (int) agemax){
-	     if(pos>=1.e-5){
-	       probs[i][jk][j1]= pp[jk]/pos;
-	     }
-	   }
-	 }
-	 
+	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++){	   
+	  if( i <= (int) agemax){
+	    if(pos>=1.e-5){
+	      probs[i][jk][j1]= pp[jk]/pos;
+	    }
+	  }
+	}
+	
+      }
     }
   }
-  
+
   
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
@@ -1531,18 +1536,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,27 +1598,114 @@ 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);  
-  
-    /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/
+ 
 
     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.;
+
+     printf("%d|",(int)age);fflush(stdout);
+    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 ******************/
@@ -1621,7 +1724,7 @@ void varevsij(char fileres[], double ***
   double age,agelim, hf;
   int theta;
 
-   fprintf(ficresvij,"# Covariances of life expectancies\n");
+  fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");
   fprintf(ficresvij,"# Age");
   for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)
@@ -1756,7 +1859,7 @@ void varprevlim(char fileres[], double *
   double age,agelim;
   int theta;
    
-  fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");
+  fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");
   fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);
@@ -1842,8 +1945,17 @@ 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);
@@ -1933,7 +2045,7 @@ void varprob(char fileres[], double **ma
 	fprintf(ficresprob,"\n%d ",(int)age);
 
 	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
-	  fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]);
+	  fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
   
       }
     }
@@ -1949,10 +2061,12 @@ void varprob(char fileres[], double **ma
 
 /******************* Printing html file ***********/
 void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
- int lastpass, int stepm, int weightopt, char model[],\
- int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
- char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
- char version[], int popforecast, int estepm ){
+		  int lastpass, int stepm, int weightopt, char model[],\
+		  int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
+		  char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
+		  char version[], int popforecast, int estepm ,/* \ */
+		  double jprev1, double mprev1,double anprev1, \
+		  double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;
   FILE *fichtm;
   /*char optionfilehtm[FILENAMELENGTH];*/
@@ -1963,25 +2077,30 @@ void printinghtml(char fileres[], char t
     printf("Problem with %s \n",optionfilehtm), exit(0);
   }
 
- fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
+  fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
 \n
 Total number of observations=%d <br>\n
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
-<hr  size=\"2\" color=\"#EC5E5E\"> 
- <ul><li>Outputs files<br>\n
+<hr  size=\"2\" color=\"#EC5E5E\">
+ <ul><li>Parameter files<br>\n
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
- - Gnuplot file name: <a href=\"%s\">%s</a><br>\n
- - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
- - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>\n
- - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>\n
- - Life expectancies by age and initial health status (estepm=%2d months): <a href=\"e%s\">e%s</a> <br>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
-
- 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
- - 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);
+ - Gnuplot file name: <a href=\"%s\">%s</a><br></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);
+
+   fprintf(fichtm,"<ul><li>Result files (first order: no variance)<br>\n
+ - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n
+ - Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
+ - Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
+ - Life expectancies by age and initial health status (estepm=%2d months): 
+   <a href=\"e%s\">e%s</a> <br>\n</li>", \
+  jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
+
+ fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n
+ - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
+ - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
+ - Variances and covariances 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 (no covariance): <a href=\"t%s\">t%s</a> <br>\n
+ - Standard deviation of stable 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
@@ -2004,24 +2123,29 @@ fprintf(fichtm," <li>Graphs</li><p>");
 	   fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][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, "."),jj1,strtok(optionfile, "."),jj1);     
+       /* Pij */
+       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>
+<img src=\"pe%s%d1.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);     
+       /* Quasi-incidences */
+       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>
+<img src=\"pe%s%d2.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);     
+       /* Stable prevalence in each health state */
        for(cpt=1; cpt<nlstate;cpt++){
-	 fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
-<img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
+	 fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
+<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }
     for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
-interval) in state (%d): v%s%d%d.gif <br>
-<img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
+interval) in state (%d): v%s%d%d.png <br>
+<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
-        fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
-<img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
+        fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
+<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and
-health expectancies in states (1) and (2): e%s%d.gif<br>
-<img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
+health expectancies in states (1) and (2): e%s%d.png<br>
+<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
 fprintf(fichtm,"\n</body>");
    }
    }
@@ -2032,7 +2156,7 @@ fclose(fichtm);
 void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
-
+  int ng;
   strcpy(optionfilegnuplot,optionfilefiname);
   strcat(optionfilegnuplot,".gp.txt");
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
@@ -2049,9 +2173,11 @@ m=pow(2,cptcoveff);
    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);
+     fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
+     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
 #endif
 #ifdef unix
+fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
 fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
 #endif
 
@@ -2071,15 +2197,15 @@ for (i=1; i<= nlstate ; i ++) {
 }  
      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");
+fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
 #endif
-fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
    }
   }
   /*2 eme*/
 
   for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);
+    fprintf(ficgp,"\nset out \"e%s%d.png\" \n",strtok(optionfile, "."),k1);
+    fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
     
     for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;
@@ -2104,27 +2230,36 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
       else fprintf(ficgp,"\" t\"\" w l 0,");
     }
-    fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
   }
  
   /*3eme*/
 
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
-      k=2+nlstate*(cpt-1);
-      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);
+      k=2+nlstate*(2*cpt-2);
+      fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
+      fprintf(ficgp,"set ter png small\nset size 0.65,0.65\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);
-    }
     }
+  }
  
   /* CV preval stat */
     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<nlstate ; cpt ++) {
       k=3;
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
+      fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
 
       for (i=1; i< nlstate ; i ++)
 	fprintf(ficgp,"+$%d",k+i+1);
@@ -2137,7 +2272,6 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
 	fprintf(ficgp,"+$%d",l+i+1);
       }
       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);
     } 
   }  
   
@@ -2153,50 +2287,58 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
 	}
       }
     }
-    }
+   }
 
-    for(jk=1; jk <=m; jk++) {
-  fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
-   i=1;
-   for(k2=1; k2<=nlstate; k2++) {
-     k3=i;
-     for(k=1; k<=(nlstate+ndeath); k++) {
-       if (k != k2){
-    	fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-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;
+   for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+     for(jk=1; jk <=m; jk++) {
+       fprintf(ficgp,"\nset out \"pe%s%d%d.png\" \n",strtok(optionfile, "."),jk,ng); 
+       if (ng==2)
+	 fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
+       else
+	 fprintf(ficgp,"\nset title \"Probability\"\n");
+       fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+       i=1;
+       for(k2=1; k2<=nlstate; k2++) {
+	 k3=i;
+	 for(k=1; k<=(nlstate+ndeath); k++) {
+	   if (k != k2){
+	     if(ng==2)
+	       fprintf(ficgp," %f*exp(p%d+p%d*x",stepm/YEARM,i,i+1);
+	     else
+	       fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+	     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;
+	   }
+	 }
        }
      }
    }
-   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
-   }
-   
-  fclose(ficgp);
+   fclose(ficgp);
 }  /* end gnuplot */
 
 
@@ -2514,7 +2656,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;
@@ -2532,7 +2674,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.8d, May 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
 
@@ -3143,7 +3285,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 
 /*--------- index.htm --------*/
 
-  printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm);
+  printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
 
   
   /*--------------- Prevalence limit --------------*/
@@ -3227,13 +3369,13 @@ 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");
 	}
@@ -3249,9 +3391,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);
@@ -3281,6 +3421,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++){
@@ -3303,7 +3445,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);
@@ -3326,9 +3469,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];
@@ -3340,7 +3485,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);