--- imach/src/imach.c	2002/02/22 18:10:15	1.24
+++ imach/src/imach.c	2002/02/26 17:11:54	1.25
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.24 2002/02/22 18:10:15 lievre Exp $
+/* $Id: imach.c,v 1.25 2002/02/26 17:11:54 lievre Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -95,7 +95,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;
-FILE *ficgp, *fichtm,*ficresprob,*ficpop;
+FILE *ficgp,*ficresprob,*ficpop;
 FILE *ficreseij;
   char filerese[FILENAMELENGTH];
  FILE  *ficresvij;
@@ -123,7 +123,7 @@ FILE *ficreseij;
 static double maxarg1,maxarg2;
 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))
 #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))
- 
+  
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
 #define rint(a) floor(a+0.5)
 
@@ -1879,6 +1879,270 @@ fclose(ficresprob);
  exit(0);
 }
 
+/******************* 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[] ){
+  int jj1, k1, i1, cpt;
+  FILE *fichtm;
+  /*char optionfilehtm[FILENAMELENGTH];*/
+
+  strcpy(optionfilehtm,optionfile);
+  strcat(optionfilehtm,".htm");
+  if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
+    printf("Problem with %s \n",optionfilehtm), exit(0);
+  }
+
+ fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\"> 
+Title=%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>
+        - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
+        - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
+        - 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>
+        - Prevalences and population 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><p>");
+
+ m=cptcoveff;
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}
+
+ 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.gif<br>
+<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
+       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);
+       }
+    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);  
+     }
+     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>- 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);
+fprintf(fichtm,"\n</body>");
+   }
+   }
+fclose(fichtm);
+}
+
+/******************* Gnuplot file **************/
+void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){
+
+  int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
+
+  strcpy(optionfilegnuplot,optionfilefiname);
+  strcat(optionfilegnuplot,".plt");
+  if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+    printf("Problem with file %s",optionfilegnuplot);
+  }
+
+#ifdef windows
+    fprintf(ficgp,"cd \"%s\" \n",pathc);
+#endif
+m=pow(2,cptcoveff);
+  
+ /* 1eme*/
+  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",agemin,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",agemin,fage,fileres);
+#endif
+
+for (i=1; i<= nlstate ; i ++) {
+  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+  else fprintf(ficgp," \%%*lf (\%%*lf)");
+}
+    fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
+    for (i=1; i<= nlstate ; i ++) {
+  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+  else fprintf(ficgp," \%%*lf (\%%*lf)");
+} 
+  fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1); 
+     for (i=1; i<= nlstate ; i ++) {
+  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
+  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);
+   }
+  }
+  /*2 eme*/
+
+  for (k1=1; k1<= m ; k1 ++) { 
+    fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
+    
+    for (i=1; i<= nlstate+1 ; i ++) {
+      k=2*i;
+      fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
+      for (j=1; j<= nlstate+1 ; j ++) {
+  if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+  else fprintf(ficgp," \%%*lf (\%%*lf)");
+}   
+      if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
+      else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
+    fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
+      for (j=1; j<= nlstate+1 ; j ++) {
+	if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+	else fprintf(ficgp," \%%*lf (\%%*lf)");
+}   
+      fprintf(ficgp,"\" t\"\" w l 0,");
+     fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
+      for (j=1; j<= nlstate+1 ; j ++) {
+  if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
+  else fprintf(ficgp," \%%*lf (\%%*lf)");
+}   
+      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",agemin,fage,fileres,k1-1,k1-1,k,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,"\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",agemin,agemax,fileres,k1,k+cpt+1,k+1);
+
+      for (i=1; i< nlstate ; i ++)
+	fprintf(ficgp,"+$%d",k+i+1);
+      fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
+      
+      l=3+(nlstate+ndeath)*cpt;
+      fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);
+      for (i=1; i< nlstate ; i ++) {
+	l=3+(nlstate+ndeath)*cpt;
+	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);
+    } 
+  }  
+  
+  /* proba elementaires */
+   for(i=1,jk=1; i <=nlstate; i++){
+    for(k=1; k <=(nlstate+ndeath); k++){
+      if (k != i) {
+	for(j=1; j <=ncovmodel; j++){
+	
+	  fprintf(ficgp,"p%d=%f ",jk,p[jk]);
+	  jk++; 
+	  fprintf(ficgp,"\n");
+	}
+      }
+    }
+    }
+
+    for(jk=1; jk <=m; jk++) {
+  fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);
+   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;
+       }
+     }
+   }
+   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
+   }
+   
+  fclose(ficgp);
+}  /* end gnuplot */
+
+
+/*************** Moving average **************/
+void movingaverage(double agedeb, double fage,double agemin, double ***mobaverage){
+
+  int i, cpt, cptcod;
+    for (agedeb=agemin; agedeb<=fage; agedeb++)
+      for (i=1; i<=nlstate;i++)
+	for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
+	  mobaverage[(int)agedeb][i][cptcod]=0.;
+    
+    for (agedeb=agemin+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;
+	}
+      }
+    }
+    
+}
+
 /***********************************************/
 /**************** Main Program *****************/
 /***********************************************/
@@ -1898,8 +2162,8 @@ int main(int argc, char *argv[])
   int *indx;
   char line[MAXLINE], linepar[MAXLINE];
   char title[MAXLINE];
-  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
-  char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
+  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
+  char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
   
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
 
@@ -2472,9 +2736,9 @@ printf("Total number of individuals= %d,
   }
   ungetc(c,ficpar);
   
-  fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
-  fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
- fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+  fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
+  fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
+ fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
      
   while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);
@@ -2499,197 +2763,27 @@ printf("Total number of individuals= %d,
     fputs(line,ficparo);
   }
   ungetc(c,ficpar);
-  fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
-fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
-fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
+  fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mob_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
+fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
+fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
 
  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
 
-    
-    /*------------ gnuplot -------------*/
-    /*chdir(pathcd);*/
-    strcpy(optionfilegnuplot,optionfilefiname);
-    strcat(optionfilegnuplot,".plt");
-    if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
-      printf("Problem with file %s",optionfilegnuplot);goto end;
-    }
-#ifdef windows
-    fprintf(ficgp,"cd \"%s\" \n",pathc);
-#endif
-m=pow(2,cptcoveff);
-  
- /* 1eme*/
-  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",agemin,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",agemin,fage,fileres);
-#endif
-
-for (i=1; i<= nlstate ; i ++) {
-  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
-  else fprintf(ficgp," \%%*lf (\%%*lf)");
-}
-    fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
-    for (i=1; i<= nlstate ; i ++) {
-  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
-  else fprintf(ficgp," \%%*lf (\%%*lf)");
-} 
-  fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1); 
-     for (i=1; i<= nlstate ; i ++) {
-  if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
-  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);
-   }
-  }
-  /*2 eme*/
-
-  for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
-    
-    for (i=1; i<= nlstate+1 ; i ++) {
-      k=2*i;
-      fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-  if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
-  else fprintf(ficgp," \%%*lf (\%%*lf)");
-}   
-      if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
-      else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
-    fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
-	else fprintf(ficgp," \%%*lf (\%%*lf)");
-}   
-      fprintf(ficgp,"\" t\"\" w l 0,");
-     fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-  if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
-  else fprintf(ficgp," \%%*lf (\%%*lf)");
-}   
-      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);
-  }
+/*------------ gnuplot -------------*/
+ printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);
  
-  /*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",agemin,fage,fileres,k1-1,k1-1,k,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,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
-    }
-  }
+/*------------ free_vector  -------------*/
+ chdir(path);
+ 
+ free_ivector(wav,1,imx);
+ free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
+ free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   
+ free_ivector(num,1,n);
+ free_vector(agedc,1,n);
+ /*free_matrix(covar,1,NCOVMAX,1,n);*/
+ fclose(ficparo);
+ fclose(ficres);
  
-  /* 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",agemin,agemax,fileres,k1,k+cpt+1,k+1);
-      for (i=1; i< nlstate ; i ++)
-	fprintf(ficgp,"+$%d",k+i+1);
-      fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
-      
-      l=3+(nlstate+ndeath)*cpt;
-      fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);
-      for (i=1; i< nlstate ; i ++) {
-	l=3+(nlstate+ndeath)*cpt;
-	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);
-    } 
-  }  
-
-  /* proba elementaires */
-   for(i=1,jk=1; i <=nlstate; i++){
-    for(k=1; k <=(nlstate+ndeath); k++){
-      if (k != i) {
-	for(j=1; j <=ncovmodel; j++){
-	  /*fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/
-	  /*fprintf(ficgp,"%s",alph[1]);*/
-	  fprintf(ficgp,"p%d=%f ",jk,p[jk]);
-	  jk++; 
-	  fprintf(ficgp,"\n");
-	}
-      }
-    }
-    }
-
-  for(jk=1; jk <=m; jk++) {
-  fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);
-   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;
-       }
-     }
-   }
-   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
-  }
-   
-  fclose(ficgp);
-  /* end gnuplot */
-   
-chdir(path);
-   
-    free_ivector(wav,1,imx);
-    free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
-    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   
-    free_ivector(num,1,n);
-    free_vector(agedc,1,n);
-    /*free_matrix(covar,1,NCOVMAX,1,n);*/
-    fclose(ficparo);
-    fclose(ficres);
-    /*  }*/
-   
-   /*________fin mle=1_________*/
-   
-
-  
-    /* No more information from the sample is required now */
   /* Reads comments: lines beginning with '#' */
   while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);
@@ -2704,68 +2798,9 @@ chdir(path);
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 /*--------- index.htm --------*/
 
-  strcpy(optionfilehtm,optionfile);
-  strcat(optionfilehtm,".htm");
-  if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
-    printf("Problem with %s \n",optionfilehtm);goto end;
-  }
-
- fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </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>
-        - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
-        - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
-        - 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>
-        - Prevalences and population 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><p>");
-
- m=cptcoveff;
- if (cptcovn < 1) {m=1;ncodemax[1]=1;}
-
- j1=0;
- for(k1=1; k1<=m;k1++){
-   for(i1=1; i1<=ncodemax[k1];i1++){
-       j1++;
-       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[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);     
-       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,j1,strtok(optionfile, "."),cpt,j1);
-       }
-    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,j1,strtok(optionfile, "."),cpt,j1);  
-     }
-     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,j1,strtok(optionfile, "."),cpt,j1);
-     }
-     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, "."),j1,strtok(optionfile, "."),j1);
-fprintf(fichtm,"\n</body>");
-   }
- }
-fclose(fichtm);
+  printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);
 
+  
   /*--------------- Prevalence limit --------------*/
   
   strcpy(filerespl,"pl");
@@ -2883,27 +2918,14 @@ fclose(fichtm);
   free_matrix(mint,1,maxwav,1,n);
   free_matrix(anint,1,maxwav,1,n);
   free_matrix(agev,1,maxwav,1,imx);
+
   /* Mobile average */
 
   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;
-	}
-      }
-    }   
+    movingaverage(agedeb, fage, agemin, mobaverage);
   }
 
   stepsize=(int) (stepm+YEARM-1)/YEARM;
@@ -2953,16 +2975,15 @@ fclose(fichtm);
       fprintf(ficresf,"# StartingAge FinalAge");
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
       if (popforecast==1)  fprintf(ficresf," [Population]");
-   
-      for (cpt=0; cpt<4;cpt++) { 
+      
+     for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { 
 	fprintf(ficresf,"\n");
 	fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   
 
-	for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */
+       	for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ 
 	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); 
 	nhstepm = nhstepm/hstepm; 
-	/*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,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);  
@@ -2978,7 +2999,7 @@ fclose(fichtm);
 		kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
 	      else {
 		kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
-	        /* fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h], probs[(int)(agedeb)+1][i][cptcod]);*/
+	        
 	      }
 
 	      if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
@@ -2991,26 +3012,29 @@ fclose(fichtm);
 	    }
 	  }
 	}
-	/*	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);*/
-      }
-      }
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	}
+     }
     }
   }
-  /*  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+
   if (popforecast==1) {
     free_ivector(popage,0,AGESUP);
     free_vector(popeffectif,0,AGESUP);
     free_vector(popcount,0,AGESUP);
   }
   free_imatrix(s,1,maxwav+1,1,n);
-  free_vector(weight,1,n);*/
+  free_vector(weight,1,n);
   fclose(ficresf);
   }/* End forecasting */
   else{
     erreur=108;
     printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);
+
   }
 
+ 
   /*---------- Health expectancies and variances ------------*/
 
   strcpy(filerest,"t");
