--- imach/src/imach.c	2002/05/21 18:44:41	1.42
+++ imach/src/imach.c	2002/05/30 17:44:35	1.46
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.42 2002/05/21 18:44:41 brouard Exp $
+/* $Id: imach.c,v 1.46 2002/05/30 17:44:35 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -96,7 +96,7 @@ double jmean; /* Mean space between 2 wa
 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,*ficresf,*ficrespop;
-FILE *ficgp,*ficresprob,*ficpop;
+FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficreseij;
   char filerese[FILENAMELENGTH];
  FILE  *ficresvij;
@@ -1641,14 +1641,10 @@ void evsij(char fileres[], double ***eij
 	  }
 	}
       }
-     
-   
-
       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 */
@@ -1658,15 +1654,15 @@ void evsij(char fileres[], double ***eij
      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];
-
+	  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++){
+     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]);
@@ -1724,7 +1720,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++)
@@ -1859,7 +1855,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);
@@ -1931,7 +1927,7 @@ void varprevlim(char fileres[], double *
 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, i1, k1, j1, z1;
-  int k=0, cptcode;
+  int k=0,l, cptcode;
   double **dnewm,**doldm;
   double *xp;
   double *gp, *gm;
@@ -1939,24 +1935,43 @@ void varprob(char fileres[], double **ma
   double age,agelim, cov[NCOVMAX];
   int theta;
   char fileresprob[FILENAMELENGTH];
+  char fileresprobcov[FILENAMELENGTH];
+  char fileresprobcor[FILENAMELENGTH];
 
   strcpy(fileresprob,"prob"); 
   strcat(fileresprob,fileres);
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);
   }
+  strcpy(fileresprobcov,"probcov"); 
+  strcat(fileresprobcov,fileres);
+  if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
+    printf("Problem with resultfile: %s\n", fileresprobcov);
+  }
+  strcpy(fileresprobcor,"probcor"); 
+  strcat(fileresprobcor,fileres);
+  if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
+    printf("Problem with resultfile: %s\n", fileresprobcor);
+  }
   printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+  printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+  printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
   
-fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
+  fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");
+  fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
+  fprintf(ficresprobcov,"# Age");
+  fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
+  fprintf(ficresprobcov,"# Age");
   for(i=1; i<=nlstate;i++)
-    for(j=1; j<=(nlstate+ndeath);j++)
+    for(j=1; j<=(nlstate+ndeath);j++){
       fprintf(ficresprob," p%1d-%1d (SE)",i,j);
-
-
+      fprintf(ficresprobcov," p%1d-%1d ",i,j);
+      fprintf(ficresprobcor," p%1d-%1d ",i,j);
+    }  
   fprintf(ficresprob,"\n");
-
-
+  fprintf(ficresprobcov,"\n");
+  fprintf(ficresprobcor,"\n");
   xp=vector(1,npar);
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
@@ -1971,15 +1986,20 @@ fprintf(ficresprob,"#One-step probabilit
 
     if  (cptcovn>0) {
       fprintf(ficresprob, "\n#********** Variable "); 
+      fprintf(ficresprobcov, "\n#********** Variable "); 
+      fprintf(ficresprobcor, "\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 (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+      fprintf(ficresprobcov, "**********\n#");
+      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+      fprintf(ficresprobcor, "**********\n#");
     }
     
       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 (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
 	for (k=1; k<=cptcovprod;k++)
@@ -2037,16 +2057,33 @@ fprintf(ficresprob,"#One-step probabilit
 	  }
 	}
      
-     /*printf("\n%d ",(int)age);
+	/*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);
+	fprintf(ficresprobcov,"\n%d ",(int)age);
+	fprintf(ficresprobcor,"\n%d ",(int)age);
 
 	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
-	  fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
-  
+	  fprintf(ficresprob,"%12.3e (%12.3e) ",gm[i],sqrt(doldm[i][j]));
+	for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+	  fprintf(ficresprobcov,"%12.3e ",gm[i]);
+	  fprintf(ficresprobcor,"%12.3e ",gm[i]);
+	}
+	i=0;
+	for (k=1; k<=(nlstate);k++){
+ 	  for (l=1; l<=(nlstate+ndeath);l++){ 
+ 	    i=i++;
+	    fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
+	    fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
+	    for (j=1; j<=i;j++){
+	      fprintf(ficresprobcov," %12.3e",doldm[i][j]);
+	      fprintf(ficresprobcor," %12.3e",doldm[i][j]/sqrt(doldm[i][i])/sqrt(doldm[j][j]));
+	    }
+	  }
+	}
       }
     }
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
@@ -2056,15 +2093,19 @@ fprintf(ficresprob,"#One-step probabilit
   }
   free_vector(xp,1,npar);
   fclose(ficresprob);
-  
+  fclose(ficresprobcov);
+  fclose(ficresprobcor);
 }
 
+
 /******************* 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];*/
@@ -2075,26 +2116,32 @@ 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
-  - 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,fileres,fileres, 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
+ - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
+ - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%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,fileres,fileres,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
@@ -2110,17 +2157,22 @@ fprintf(fichtm," <li>Graphs</li><p>");
  jj1=0;
  for(k1=1; k1<=m;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){
-       jj1++;
-       if (cptcovn > 0) {
-	 fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
-	 for (cpt=1; cpt<=cptcoveff;cpt++) 
-	   fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
-	 fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
-       }
-       fprintf(fichtm,"<br>- Probabilities: pe%s%d.png<br>
-<img src=\"pe%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
+     jj1++;
+     if (cptcovn > 0) {
+       fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
+       for (cpt=1; cpt<=cptcoveff;cpt++) 
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+       fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
+     }
+     /* 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\">",stepm,strtok(optionfile, "."),jj1,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\">",stepm,strtok(optionfile, "."),jj1,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.png<br>
+	 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++) {
@@ -2137,7 +2189,7 @@ health expectancies in states (1) and (2
 <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
 fprintf(fichtm,"\n</body>");
    }
-   }
+ }
 fclose(fichtm);
 }
 
@@ -2145,9 +2197,9 @@ 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");
+  strcat(optionfilegnuplot,".gp");
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     printf("Problem with file %s",optionfilegnuplot);
   }
@@ -2162,11 +2214,11 @@ m=pow(2,cptcoveff);
    for (k1=1; k1<= m ; k1 ++) {
 
 #ifdef windows
-     fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+     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\n",strtok(optionfile, "."),cpt,k1);
+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
 
@@ -2193,7 +2245,7 @@ fprintf(ficgp,"\nset ter png small\nset
   /*2 eme*/
 
   for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1);
+    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 ++) {
@@ -2226,7 +2278,7 @@ fprintf(ficgp,"\nset ter png small\nset
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       k=2+nlstate*(2*cpt-2);
-      fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+      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) ");
@@ -2247,7 +2299,7 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<nlstate ; cpt ++) {
       k=3;
-      fprintf(ficgp,"set out \"p%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+      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 ++)
@@ -2278,47 +2330,55 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
     }
    }
 
-   for(jk=1; jk <=m; jk++) {
-     fprintf(ficgp,"\nset out \"pe%s%d.png\" \n\n",strtok(optionfile, "."),jk); 
-     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){
-	   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++;
-	     }
+   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",YEARM/stepm,i,i+1);
 	     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);
+	       fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
 	     ij=1;
-	     for(j=3; j <=ncovmodel; j++){
+	     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]]]);
+		 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",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
 	     }
-	     fprintf(ficgp,")");
+	     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,") t \"p%d%d\" ", k2,k);
-	   if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
-	   i=i+ncovmodel;
 	 }
        }
      }
    }
-   
    fclose(ficgp);
 }  /* end gnuplot */
 
@@ -2655,7 +2715,7 @@ int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   
 
-  char version[80]="Imach version 0.8d, May 2002, INED-EUROREVES ";
+  char version[80]="Imach version 0.8g, May 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
 
@@ -3266,7 +3326,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 --------------*/