--- imach/src/imach.c	2004/05/16 15:05:56	1.98
+++ imach/src/imach.c	2006/03/16 17:45:01	1.121
@@ -1,6 +1,103 @@
-/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $
+/* $Id: imach.c,v 1.121 2006/03/16 17:45:01 lievre Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.121  2006/03/16 17:45:01  lievre
+  * imach.c (Module): Comments concerning covariates added
+
+  * imach.c (Module): refinements in the computation of lli if
+  status=-2 in order to have more reliable computation if stepm is
+  not 1 month. Version 0.98f
+
+  Revision 1.120  2006/03/16 15:10:38  lievre
+  (Module): refinements in the computation of lli if
+  status=-2 in order to have more reliable computation if stepm is
+  not 1 month. Version 0.98f
+
+  Revision 1.119  2006/03/15 17:42:26  brouard
+  (Module): Bug if status = -2, the loglikelihood was
+  computed as likelihood omitting the logarithm. Version O.98e
+
+  Revision 1.118  2006/03/14 18:20:07  brouard
+  (Module): varevsij Comments added explaining the second
+  table of variances if popbased=1 .
+  (Module): Covariances of eij, ekl added, graphs fixed, new html link.
+  (Module): Function pstamp added
+  (Module): Version 0.98d
+
+  Revision 1.117  2006/03/14 17:16:22  brouard
+  (Module): varevsij Comments added explaining the second
+  table of variances if popbased=1 .
+  (Module): Covariances of eij, ekl added, graphs fixed, new html link.
+  (Module): Function pstamp added
+  (Module): Version 0.98d
+
+  Revision 1.116  2006/03/06 10:29:27  brouard
+  (Module): Variance-covariance wrong links and
+  varian-covariance of ej. is needed (Saito).
+
+  Revision 1.115  2006/02/27 12:17:45  brouard
+  (Module): One freematrix added in mlikeli! 0.98c
+
+  Revision 1.114  2006/02/26 12:57:58  brouard
+  (Module): Some improvements in processing parameter
+  filename with strsep.
+
+  Revision 1.113  2006/02/24 14:20:24  brouard
+  (Module): Memory leaks checks with valgrind and:
+  datafile was not closed, some imatrix were not freed and on matrix
+  allocation too.
+
+  Revision 1.112  2006/01/30 09:55:26  brouard
+  (Module): Back to gnuplot.exe instead of wgnuplot.exe
+
+  Revision 1.111  2006/01/25 20:38:18  brouard
+  (Module): Lots of cleaning and bugs added (Gompertz)
+  (Module): Comments can be added in data file. Missing date values
+  can be a simple dot '.'.
+
+  Revision 1.110  2006/01/25 00:51:50  brouard
+  (Module): Lots of cleaning and bugs added (Gompertz)
+
+  Revision 1.109  2006/01/24 19:37:15  brouard
+  (Module): Comments (lines starting with a #) are allowed in data.
+
+  Revision 1.108  2006/01/19 18:05:42  lievre
+  Gnuplot problem appeared...
+  To be fixed
+
+  Revision 1.107  2006/01/19 16:20:37  brouard
+  Test existence of gnuplot in imach path
+
+  Revision 1.106  2006/01/19 13:24:36  brouard
+  Some cleaning and links added in html output
+
+  Revision 1.105  2006/01/05 20:23:19  lievre
+  *** empty log message ***
+
+  Revision 1.104  2005/09/30 16:11:43  lievre
+  (Module): sump fixed, loop imx fixed, and simplifications.
+  (Module): If the status is missing at the last wave but we know
+  that the person is alive, then we can code his/her status as -2
+  (instead of missing=-1 in earlier versions) and his/her
+  contributions to the likelihood is 1 - Prob of dying from last
+  health status (= 1-p13= p11+p12 in the easiest case of somebody in
+  the healthy state at last known wave). Version is 0.98
+
+  Revision 1.103  2005/09/30 15:54:49  lievre
+  (Module): sump fixed, loop imx fixed, and simplifications.
+
+  Revision 1.102  2004/09/15 17:31:30  brouard
+  Add the possibility to read data file including tab characters.
+
+  Revision 1.101  2004/09/15 10:38:38  brouard
+  Fix on curr_time
+
+  Revision 1.100  2004/07/12 18:29:06  brouard
+  Add version for Mac OS X. Just define UNIX in Makefile
+
+  Revision 1.99  2004/06/05 08:57:40  brouard
+  *** empty log message ***
+
   Revision 1.98  2004/05/16 15:05:56  brouard
   New version 0.97 . First attempt to estimate force of mortality
   directly from the data i.e. without the need of knowing the health
@@ -143,7 +240,7 @@
   hPijx.
 
   Also this programme outputs the covariance matrix of the parameters but also
-  of the life expectancies. It also computes the stable prevalence. 
+  of the life expectancies. It also computes the period (stable) prevalence. 
   
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.
@@ -172,7 +269,7 @@
       begin-prev-date,...
   open gnuplot file
   open html file
-  stable prevalence
+  period (stable) prevalence
    for age prevalim()
   h Pij x
   variance of p varprob
@@ -184,7 +281,7 @@
   varevsij() 
   if popbased==1 varevsij(,popbased)
   total life expectancies
-  Variance of stable prevalence
+  Variance of period (stable) prevalence
  end
 */
 
@@ -194,9 +291,16 @@
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <unistd.h>
 
-#include <sys/time.h>
+#include <limits.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <errno.h>
+extern int errno;
+
+/* #include <sys/time.h> */
 #include <time.h>
 #include "timeval.h"
 
@@ -204,11 +308,11 @@
 /* #define _(String) gettext (String) */
 
 #define MAXLINE 256
+
 #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 132
-/*#define DEBUG*/
-/*#define windows*/
+
 #define	GLOCK_ERROR_NOPATH		-1	/* empty path */
 #define	GLOCK_ERROR_GETCWD		-2	/* cannot get cwd */
 
@@ -224,19 +328,23 @@
 #define AGESUP 130
 #define AGEBASE 40
 #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
-#ifdef unix
+#ifdef UNIX
 #define DIRSEPARATOR '/'
+#define CHARSEPARATOR "/"
 #define ODIRSEPARATOR '\\'
 #else
 #define DIRSEPARATOR '\\'
+#define CHARSEPARATOR "\\"
 #define ODIRSEPARATOR '/'
 #endif
 
-/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ */
+/* $Id: imach.c,v 1.121 2006/03/16 17:45:01 lievre Exp $ */
 /* $State: Exp $ */
 
-char version[]="Imach version 0.70, May 2004, INED-EUROREVES ";
-char fullversion[]="$Revision: 1.98 $ $Date: 2004/05/16 15:05:56 $"; 
+char version[]="Imach version 0.98f, March 2006, INED-EUROREVES-Institut de longevite ";
+char fullversion[]="$Revision: 1.121 $ $Date: 2006/03/16 17:45:01 $"; 
+char strstart[80];
+char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -249,6 +357,7 @@ int popbased=0;
 int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */
 int jmin, jmax; /* min, max spacing between 2 waves */
+int ijmin, ijmax; /* Individuals having jmin and jmax */ 
 int gipmx, gsw; /* Global variables on the number of contributions 
 		   to the likelihood and the sum of weights (done by funcone)*/
 int mle, weightopt;
@@ -259,7 +368,7 @@ int **bh; /* bh[mi][i] is the bias (+ or
 double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
-FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
+FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficlog, *ficrespow;
 int globpr; /* Global variable for printing or not */
 double fretone; /* Only one call to likelihood */
@@ -273,13 +382,17 @@ FILE *ficresprobmorprev;
 FILE *fichtm, *fichtmcov; /* Html File */
 FILE *ficreseij;
 char filerese[FILENAMELENGTH];
+FILE *ficresstdeij;
+char fileresstde[FILENAMELENGTH];
+FILE *ficrescveij;
+char filerescve[FILENAMELENGTH];
 FILE  *ficresvij;
 char fileresv[FILENAMELENGTH];
 FILE  *ficresvpl;
 char fileresvpl[FILENAMELENGTH];
 char title[MAXLINE];
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
-char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
+char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
 char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH]; 
 char command[FILENAMELENGTH];
 int  outcmd=0;
@@ -301,6 +414,9 @@ long time_value;
 extern long time();
 char strcurr[80], strfor[80];
 
+char *endptr;
+long lval;
+
 #define NR_END 1
 #define FREE_ARG char*
 #define FTOL 1.0e-10
@@ -349,6 +465,7 @@ double *weight;
 int **s; /* Status */
 double *agedc, **covar, idx;
 int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
+double *lsurv, *lpop, *tpop;
 
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */
 double ftolhess; /* Tolerance for computing hessian */
@@ -356,13 +473,17 @@ double ftolhess; /* Tolerance for comput
 /**************** split *************************/
 static	int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {
+  /* From a file name with (full) path (either Unix or Windows) we extract the directory (dirc)
+     the name of the file (name), its extension only (ext) and its first part of the name (finame)
+  */ 
   char	*ss;				/* pointer */
   int	l1, l2;				/* length counters */
 
   l1 = strlen(path );			/* length of path */
   if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   ss= strrchr( path, DIRSEPARATOR );		/* find last / */
-  if ( ss == NULL ) {			/* no directory, so use current */
+  if ( ss == NULL ) {			/* no directory, so determine current directory */
+    strcpy( name, path );		/* we got the fullname name because no directory */
     /*if(strrchr(path, ODIRSEPARATOR )==NULL)
       printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
     /* get current working directory */
@@ -370,7 +491,8 @@ static	int split( char *path, char *dirc
     if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
       return( GLOCK_ERROR_GETCWD );
     }
-    strcpy( name, path );		/* we've got it */
+    /* got dirc from getcwd*/
+    printf(" DIRC = %s \n",dirc);
   } else {				/* strip direcotry from path */
     ss++;				/* after this, the filename */
     l2 = strlen( ss );			/* length of filename */
@@ -378,21 +500,25 @@ static	int split( char *path, char *dirc
     strcpy( name, ss );		/* save file name */
     strncpy( dirc, path, l1 - l2 );	/* now the directory */
     dirc[l1-l2] = 0;			/* add zero */
+    printf(" DIRC2 = %s \n",dirc);
   }
+  /* We add a separator at the end of dirc if not exists */
   l1 = strlen( dirc );			/* length of directory */
-  /*#ifdef windows
-  if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
-#else
-  if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
-#endif
-  */
+  if( dirc[l1-1] != DIRSEPARATOR ){
+    dirc[l1] =  DIRSEPARATOR;
+    dirc[l1+1] = 0; 
+    printf(" DIRC3 = %s \n",dirc);
+  }
   ss = strrchr( name, '.' );		/* find last / */
-  ss++;
-  strcpy(ext,ss);			/* save extension */
-  l1= strlen( name);
-  l2= strlen(ss)+1;
-  strncpy( finame, name, l1-l2);
-  finame[l1-l2]= 0;
+  if (ss >0){
+    ss++;
+    strcpy(ext,ss);			/* save extension */
+    l1= strlen( name);
+    l2= strlen(ss)+1;
+    strncpy( finame, name, l1-l2);
+    finame[l1-l2]= 0;
+  }
+
   return( 0 );				/* we're done */
 }
 
@@ -425,8 +551,8 @@ int nbocc(char *s, char occ)
 
 void cutv(char *u,char *v, char*t, char occ)
 {
-  /* cuts string t into u and v where u is ended by char occ excluding it
-     and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)
+  /* cuts string t into u and v where u ends before first occurence of char 'occ' 
+     and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')
      gives u="abcedf" and v="ghi2j" */
   int i,lg,j,p=0;
   i=0;
@@ -868,11 +994,11 @@ void powell(double p[], double **xi, int
     fprintf(ficrespow,"\n");fflush(ficrespow);
     if(*iter <=3){
       tm = *localtime(&curr_time.tv_sec);
-      strcpy(strcurr,asctime(&tmf));
+      strcpy(strcurr,asctime(&tm));
 /*       asctime_r(&tm,strcurr); */
-      forecast_time=curr_time;
+      forecast_time=curr_time; 
       itmp = strlen(strcurr);
-      if(strcurr[itmp-1]=='\n')
+      if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
 	strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
@@ -884,8 +1010,8 @@ void powell(double p[], double **xi, int
 	itmp = strlen(strfor);
 	if(strfor[itmp-1]=='\n')
 	strfor[itmp-1]='\0';
-	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
-	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
       }
     }
     for (i=1;i<=n;i++) { 
@@ -979,7 +1105,7 @@ void powell(double p[], double **xi, int
   } 
 } 
 
-/**** Prevalence limit (stable prevalence)  ****************/
+/**** Prevalence limit (stable or period prevalence)  ****************/
 
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
 {
@@ -1050,57 +1176,57 @@ double **pmij(double **ps, double *cov,
   int i,j,j1, nc, ii, jj;
 
     for(i=1; i<= nlstate; i++){
-    for(j=1; j<i;j++){
-      for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	/*s2 += param[i][j][nc]*cov[nc];*/
-	s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-	/*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/
-      }
-      ps[i][j]=s2;
-      /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/
-    }
-    for(j=i+1; j<=nlstate+ndeath;j++){
-      for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-	/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
+      for(j=1; j<i;j++){
+	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+	  /*s2 += param[i][j][nc]*cov[nc];*/
+	  s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+/* 	 printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
+	}
+	ps[i][j]=s2;
+/* 	printf("s1=%.17e, s2=%.17e\n",s1,s2); */
+      }
+      for(j=i+1; j<=nlstate+ndeath;j++){
+	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+	  s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+/* 	  printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
+	}
+	ps[i][j]=s2;
       }
-      ps[i][j]=s2;
     }
-  }
     /*ps[3][2]=1;*/
-
-  for(i=1; i<= nlstate; i++){
-     s1=0;
-    for(j=1; j<i; j++)
-      s1+=exp(ps[i][j]);
-    for(j=i+1; j<=nlstate+ndeath; j++)
-      s1+=exp(ps[i][j]);
-    ps[i][i]=1./(s1+1.);
-    for(j=1; j<i; j++)
-      ps[i][j]= exp(ps[i][j])*ps[i][i];
-    for(j=i+1; j<=nlstate+ndeath; j++)
-      ps[i][j]= exp(ps[i][j])*ps[i][i];
-    /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
-  } /* end i */
-
-  for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
-    for(jj=1; jj<= nlstate+ndeath; jj++){
-      ps[ii][jj]=0;
-      ps[ii][ii]=1;
+    
+    for(i=1; i<= nlstate; i++){
+      s1=0;
+      for(j=1; j<i; j++)
+	s1+=exp(ps[i][j]);
+      for(j=i+1; j<=nlstate+ndeath; j++)
+	s1+=exp(ps[i][j]);
+      ps[i][i]=1./(s1+1.);
+      for(j=1; j<i; j++)
+	ps[i][j]= exp(ps[i][j])*ps[i][i];
+      for(j=i+1; j<=nlstate+ndeath; j++)
+	ps[i][j]= exp(ps[i][j])*ps[i][i];
+      /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+    } /* end i */
+    
+    for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+      for(jj=1; jj<= nlstate+ndeath; jj++){
+	ps[ii][jj]=0;
+	ps[ii][ii]=1;
+      }
     }
-  }
-
+    
 
-  /*   for(ii=1; ii<= nlstate+ndeath; ii++){
-    for(jj=1; jj<= nlstate+ndeath; jj++){
-     printf("%lf ",ps[ii][jj]);
-   }
-    printf("\n ");
-    }
-    printf("\n ");printf("%lf ",cov[2]);*/
-/*
-  for(i=1; i<= npar; i++) printf("%f ",x[i]);
-  goto end;*/
+/*        for(ii=1; ii<= nlstate+ndeath; ii++){ */
+/* 	 for(jj=1; jj<= nlstate+ndeath; jj++){ */
+/* 	   printf("ddd %lf ",ps[ii][jj]); */
+/* 	 } */
+/* 	 printf("\n "); */
+/*        } */
+/*        printf("\n ");printf("%lf ",cov[2]); */
+       /*
+      for(i=1; i<= npar; i++) printf("%f ",x[i]);
+      goto end;*/
     return ps;
 }
 
@@ -1224,14 +1350,14 @@ double func( double *x)
 	} /* end mult */
       
 	/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
-	/* But now since version 0.9 we anticipate for bias and large stepm.
+	/* But now since version 0.9 we anticipate for bias at large stepm.
 	 * If stepm is larger than one month (smallest stepm) and if the exact delay 
 	 * (in months) between two waves is not a multiple of stepm, we rounded to 
 	 * the nearest (and in case of equal distance, to the lowest) interval but now
 	 * we keep into memory the bias bh[mi][i] and also the previous matrix product
-	 * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+	 * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
 	 * probability in order to take into account the bias as a fraction of the way
-	 * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+	 * from savm to out if bh is negative or even beyond if bh is positive. bh varies
 	 * -stepm/2 to stepm/2 .
 	 * For stepm=1 the results are the same as for previous versions of Imach.
 	 * For stepm > 1 the results are less biased than in previous versions. 
@@ -1239,15 +1365,16 @@ double func( double *x)
 	s1=s[mw[mi][i]][i];
 	s2=s[mw[mi+1][i]][i];
 	bbh=(double)bh[mi][i]/(double)stepm; 
-	/* bias is positive if real duration
+	/* bias bh is positive if real duration
 	 * is higher than the multiple of stepm and negative otherwise.
 	 */
 	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
 	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 . 
+	  /* 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 equal to probability to die before dh 
+	     minus 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
@@ -1267,7 +1394,28 @@ double func( double *x)
 	lower mortality.
 	  */
 	  lli=log(out[s1][s2] - savm[s1][s2]);
-	}else{
+
+
+	} else if  (s2==-2) {
+	  for (j=1,survp=0. ; j<=nlstate; j++) 
+	    survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+	  /*survp += out[s1][j]; */
+	  lli= log(survp);
+	}
+	
+ 	else if  (s2==-4) { 
+	  for (j=3,survp=0. ; j<=nlstate; j++)  
+	    survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ 	  lli= log(survp); 
+ 	} 
+
+ 	else if  (s2==-5) { 
+ 	  for (j=1,survp=0. ; j<=2; j++)  
+	    survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+ 	  lli= log(survp); 
+ 	} 
+	
+	else{
 	  lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
 	  /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
 	} 
@@ -1461,7 +1609,11 @@ double funcone( double *x)
        */
       if( s2 > nlstate && (mle <5) ){  /* Jackson */
 	lli=log(out[s1][s2] - savm[s1][s2]);
-      } else if (mle==1){
+      } else if  (s2==-2) {
+	for (j=1,survp=0. ; j<=nlstate; j++) 
+	  survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
+	lli= log(survp);
+      }else if (mle==1){
 	lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
       } else if(mle==2){
 	lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
@@ -1477,8 +1629,8 @@ double funcone( double *x)
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
 /*       printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){
-	fprintf(ficresilk,"%9d %6d %1d %1d %1d %1d %3d %10.6f %6.4f\
- %10.6f %10.6f %10.6f ", \
+	fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
+ %11.6f %11.6f %11.6f ", \
 		num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
 		2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
 	for(k=1,llt=0.,l=0.; k<=nlstate; k++){
@@ -1563,6 +1715,7 @@ void mlikeli(FILE *ficres,double p[], in
 
   powell(p,xi,npar,ftol,&iter,&fret,func);
 
+  free_matrix(xi,1,npar,1,npar);
   fclose(ficrespow);
   printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
@@ -1830,8 +1983,13 @@ void lubksb(double **a, int n, int *indx
   } 
 } 
 
+void pstamp(FILE *fichier)
+{
+  fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);
+}
+
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint)
+void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
 {  /* Some frequencies */
   
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;
@@ -1839,7 +1997,6 @@ void  freqsummary(char fileres[], int ia
   double ***freq; /* Frequencies */
   double *pp, **prop;
   double pos,posprop, k2, dateintsum=0,k2cpt=0;
-  FILE *ficresp;
   char fileresp[FILENAMELENGTH];
   
   pp=vector(1,nlstate);
@@ -1851,7 +2008,7 @@ void  freqsummary(char fileres[], int ia
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);
   }
-  freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);
+  freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
   j1=0;
   
   j=cptcoveff;
@@ -1864,8 +2021,8 @@ void  freqsummary(char fileres[], int ia
       j1++;
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
 	scanf("%d", i);*/
-      for (i=-1; i<=nlstate+ndeath; i++)  
-	for (jk=-1; jk<=nlstate+ndeath; jk++)  
+      for (i=-5; i<=nlstate+ndeath; i++)  
+	for (jk=-5; jk<=nlstate+ndeath; jk++)  
 	  for(m=iagemin; m <= iagemax+3; m++)
 	    freq[i][jk][m]=0;
 
@@ -1904,7 +2061,7 @@ void  freqsummary(char fileres[], int ia
       }
        
       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
-
+      pstamp(ficresp);
       if  (cptcovn>0) {
 	fprintf(ficresp, "\n#********** Variable "); 
 	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
@@ -1990,7 +2147,7 @@ void  freqsummary(char fileres[], int ia
   dateintmean=dateintsum/k2cpt; 
  
   fclose(ficresp);
-  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);
+  free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
   free_vector(pp,1,nlstate);
   free_matrix(prop,1,nlstate,iagemin, iagemax+3);
   /* End of Freq */
@@ -2099,7 +2256,7 @@ void  concatwav(int wav[], int **dh, int
     mi=0;
     m=firstpass;
     while(s[m][i] <= nlstate){
-      if(s[m][i]>=1)
+      if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)
 	mw[++mi][i]=m;
       if(m >=lastpass)
 	break;
@@ -2117,11 +2274,11 @@ void  concatwav(int wav[], int **dh, int
     if(mi==0){
       nbwarn++;
       if(first==0){
-	printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);
+	printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
 	first=1;
       }
       if(first==1){
-	fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i);
+	fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
       }
     } /* end mi==0 */
   } /* End individuals */
@@ -2139,13 +2296,19 @@ void  concatwav(int wav[], int **dh, int
 	      nberr++;
 	      printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
 	      j=1; /* Temporary Dangerous patch */
-	      printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);
+	      printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
 	      fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-	      fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);
+	      fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
 	    }
 	    k=k+1;
-	    if (j >= jmax) jmax=j;
-	    if (j <= jmin) jmin=j;
+	    if (j >= jmax){
+	      jmax=j;
+	      ijmax=i;
+	    }
+	    if (j <= jmin){
+	      jmin=j;
+	      ijmin=i;
+	    }
 	    sum=sum+j;
 	    /*if (j<0) printf("j=%d num=%d \n",j,i);*/
 	    /*	  printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
@@ -2153,10 +2316,17 @@ void  concatwav(int wav[], int **dh, int
 	}
 	else{
 	  j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
-	  /*	  printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
+/* 	  if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
+
 	  k=k+1;
-	  if (j >= jmax) jmax=j;
-	  else if (j <= jmin)jmin=j;
+	  if (j >= jmax) {
+	    jmax=j;
+	    ijmax=i;
+	  }
+	  else if (j <= jmin){
+	    jmin=j;
+	    ijmin=i;
+	  }
 	  /*	    if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
 	  /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
 	  if(j<0){
@@ -2199,8 +2369,8 @@ void  concatwav(int wav[], int **dh, int
     } /* end wave */
   }
   jmean=sum/k;
-  printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
-  fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
+  printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
+  fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
  }
 
 /*********** Tricode ****************************/
@@ -2246,7 +2416,7 @@ void tricode(int *Tvar, int **nbcode, in
  for (k=0; k< maxncov; k++) Ndum[k]=0;
 
  for (i=1; i<=ncovmodel-2; i++) { 
-   /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
+   /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
    ij=Tvar[i];
    Ndum[ij]++;
  }
@@ -2264,31 +2434,151 @@ 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,double delti[],double **matcov )
+void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
 
 {
-  /* Health expectancies */
-  int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
+  /* Health expectancies, no variances */
+  int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
   double age, agelim, hf;
-  double ***p3mat,***varhe;
+  double ***p3mat;
+  double eip;
+
+  pstamp(ficreseij);
+  fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
+  fprintf(ficreseij,"# Age");
+  for(i=1; i<=nlstate;i++){
+    for(j=1; j<=nlstate;j++){
+      fprintf(ficreseij," e%1d%1d ",i,j);
+    }
+    fprintf(ficreseij," e%1d. ",i);
+  }
+  fprintf(ficreseij,"\n");
+
+  
+  if(estepm < stepm){
+    printf ("Problem %d lower than %d\n",estepm, stepm);
+  }
+  else  hstepm=estepm;   
+  /* We compute the life expectancy from trapezoids spaced every estepm months
+   * This is mainly to measure the difference between two models: for example
+   * if stepm=24 months pijx are given only every 2 years and by summing them
+   * we are calculating an estimate of the Life Expectancy assuming a linear 
+   * progression in between and thus overestimating or underestimating according
+   * to the curvature of the survival function. If, for the same date, we 
+   * estimate the model with stepm=1 month, we can keep estepm to 24 months
+   * to compare the new estimate of Life expectancy with the same linear 
+   * hypothesis. A more precise result, taking into account a more precise
+   * curvature will be obtained if estepm is as small as stepm. */
+
+  /* For example we decided to compute the life expectancy with the smallest unit */
+  /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. 
+     nhstepm is the number of hstepm from age to agelim 
+     nstepm is the number of stepm from age to agelin. 
+     Look at hpijx to understand the reason of that which relies in memory size
+     and note for a fixed period like estepm months */
+  /* We decided (b) to get a life expectancy respecting the most precise curvature of the
+     survival function given by stepm (the optimization length). Unfortunately it
+     means that if the survival funtion is printed only each two years of age and if
+     you sum them up and add 1 year (area under the trapezoids) you won't get the same 
+     results. So we changed our mind and took the option of the best precision.
+  */
+  hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ 
+
+  agelim=AGESUP;
+  /* nhstepm age range expressed in number of stepm */
+  nstepm=(int) rint((agelim-age)*YEARM/stepm); 
+  /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
+  /* 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);
+
+  for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
+    /* 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, cij);  
+    
+    hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
+    
+    printf("%d|",(int)age);fflush(stdout);
+    fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
+    
+    /* 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]);*/
+
+	}
+
+    fprintf(ficreseij,"%3.0f",age );
+    for(i=1; i<=nlstate;i++){
+      eip=0;
+      for(j=1; j<=nlstate;j++){
+	eip +=eij[i][j][(int)age];
+	fprintf(ficreseij,"%9.4f", eij[i][j][(int)age] );
+      }
+      fprintf(ficreseij,"%9.4f", eip );
+    }
+    fprintf(ficreseij,"\n");
+    
+  }
+  free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+  printf("\n");
+  fprintf(ficlog,"\n");
+  
+}
+
+void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
+
+{
+  /* Covariances of health expectancies eij and of total life expectancies according
+   to initial status i, ei. .
+  */
+  int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
+  double age, agelim, hf;
+  double ***p3matp, ***p3matm, ***varhe;
   double **dnewm,**doldm;
-  double *xp;
+  double *xp, *xm;
   double **gp, **gm;
   double ***gradg, ***trgradg;
   int theta;
 
+  double eip, vip;
+
   varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage);
   xp=vector(1,npar);
+  xm=vector(1,npar);
   dnewm=matrix(1,nlstate*nlstate,1,npar);
   doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);
   
-  fprintf(ficreseij,"# Health expectancies\n");
-  fprintf(ficreseij,"# Age");
-  for(i=1; i<=nlstate;i++)
+  pstamp(ficresstdeij);
+  fprintf(ficresstdeij,"# Health expectancies with standard errors\n");
+  fprintf(ficresstdeij,"# Age");
+  for(i=1; i<=nlstate;i++){
     for(j=1; j<=nlstate;j++)
-      fprintf(ficreseij," %1d-%1d (SE)",i,j);
-  fprintf(ficreseij,"\n");
+      fprintf(ficresstdeij," e%1d%1d (SE)",i,j);
+    fprintf(ficresstdeij," e%1d. ",i);
+  }
+  fprintf(ficresstdeij,"\n");
 
+  pstamp(ficrescveij);
+  fprintf(ficrescveij,"# Subdiagonal matrix of covariances of health expectancies by age: cov(eij,ekl)\n");
+  fprintf(ficrescveij,"# Age");
+  for(i=1; i<=nlstate;i++)
+    for(j=1; j<=nlstate;j++){
+      cptj= (j-1)*nlstate+i;
+      for(i2=1; i2<=nlstate;i2++)
+	for(j2=1; j2<=nlstate;j2++){
+	  cptj2= (j2-1)*nlstate+i2;
+	  if(cptj2 <= cptj)
+	    fprintf(ficrescveij,"  %1d%1d,%1d%1d",i,j,i2,j2);
+	}
+    }
+  fprintf(ficrescveij,"\n");
+  
   if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);
   }
@@ -2318,77 +2608,64 @@ void evsij(char fileres[], double ***eij
   */
   hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ 
 
+  /* If stepm=6 months */
+  /* nhstepm age range expressed in number of stepm */
   agelim=AGESUP;
-  for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
-    /* nhstepm age range expressed in number of stepm */
-    nstepm=(int) rint((agelim-age)*YEARM/stepm); 
-    /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
-    /* 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*nlstate);
-    gp=matrix(0,nhstepm,1,nlstate*nlstate);
-    gm=matrix(0,nhstepm,1,nlstate*nlstate);
+  nstepm=(int) rint((agelim-age)*YEARM/stepm); 
+  /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
+  /* if (stepm >= YEARM) hstepm=1;*/
+  nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
+  
+  p3matp=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+  p3matm=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+  gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);
+  trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
+  gp=matrix(0,nhstepm,1,nlstate*nlstate);
+  gm=matrix(0,nhstepm,1,nlstate*nlstate);
+
+  for (age=bage; age<=fage; age ++){ 
 
     /* 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++){
+    /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
+       decrease memory allocation */
+    for(theta=1; theta <=npar; theta++){
       for(i=1; i<=npar; i++){ 
 	xp[i] = x[i] + (i==theta ?delti[theta]:0);
+	xm[i] = x[i] - (i==theta ?delti[theta]:0);
       }
-      hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
+      hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);  
+      hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);  
   
-      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(h=0; h<=nhstepm-1; h++){
+	    gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
+	    gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[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*nlstate; j++)
+      for(ij=1; ij<= nlstate*nlstate; ij++)
 	for(h=0; h<=nhstepm-1; h++){
-	  gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
+	  gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
 	}
-     } 
-   
-/* End theta */
-
-     trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
-
-     for(h=0; h<=nhstepm-1; h++)
+    }/* End theta */
+    
+    
+    for(h=0; h<=nhstepm-1; h++)
       for(j=1; j<=nlstate*nlstate;j++)
 	for(theta=1; theta <=npar; theta++)
 	  trgradg[h][j][theta]=gradg[h][theta][j];
-     
+    
 
-     for(i=1;i<=nlstate*nlstate;i++)
-      for(j=1;j<=nlstate*nlstate;j++)
-	varhe[i][j][(int)age] =0.;
+     for(ij=1;ij<=nlstate*nlstate;ij++)
+      for(ji=1;ji<=nlstate*nlstate;ji++)
+	varhe[ij][ji][(int)age] =0.;
 
      printf("%d|",(int)age);fflush(stdout);
      fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
@@ -2396,39 +2673,60 @@ void evsij(char fileres[], double ***eij
       for(k=0;k<=nhstepm-1;k++){
 	matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
 	matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
-	for(i=1;i<=nlstate*nlstate;i++)
-	  for(j=1;j<=nlstate*nlstate;j++)
-	    varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
+	for(ij=1;ij<=nlstate*nlstate;ij++)
+	  for(ji=1;ji<=nlstate*nlstate;ji++)
+	    varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }
     }
     /* Computing expectancies */
+    hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
     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;
+	  eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[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;
+    fprintf(ficresstdeij,"%3.0f",age );
+    for(i=1; i<=nlstate;i++){
+      eip=0.;
+      vip=0.;
+      for(j=1; j<=nlstate;j++){
+	eip += eij[i][j][(int)age];
+	for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
+	  vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
+	fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
+      }
+      fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
+    }
+    fprintf(ficresstdeij,"\n");
+
+    fprintf(ficrescveij,"%3.0f",age );
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){
-	cptj++;
-	fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
+	cptj= (j-1)*nlstate+i;
+	for(i2=1; i2<=nlstate;i2++)
+	  for(j2=1; j2<=nlstate;j2++){
+	    cptj2= (j2-1)*nlstate+i2;
+	    if(cptj2 <= cptj)
+	      fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
+	  }
       }
-    fprintf(ficreseij,"\n");
+    fprintf(ficrescveij,"\n");
    
-    free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
-    free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
-    free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
-    free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
-    free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }
+  free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
+  free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
+  free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
+  free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
+  free_ma3x(p3matm,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+  free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   printf("\n");
   fprintf(ficlog,"\n");
 
+  free_vector(xm,1,npar);
   free_vector(xp,1,npar);
   free_matrix(dnewm,1,nlstate*nlstate,1,npar);
   free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
@@ -2436,7 +2734,7 @@ void evsij(char fileres[], double ***eij
 }
 
 /************ Variance ******************/
-void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav)
+void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
 {
   /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
@@ -2487,7 +2785,9 @@ void varevsij(char optionfilefiname[], d
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
   }
   printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+ 
   fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+  pstamp(ficresprobmorprev);
   fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
   fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
   for(j=nlstate+1; j<=(nlstate+ndeath);j++){
@@ -2497,16 +2797,21 @@ void varevsij(char optionfilefiname[], d
   }  
   fprintf(ficresprobmorprev,"\n");
   fprintf(ficgp,"\n# Routine varevsij");
+  /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
   fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
 /*   } */
   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-
-  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");
+  pstamp(ficresvij);
+  fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
+  if(popbased==1)
+    fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)");
+  else
+    fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
   fprintf(ficresvij,"# Age");
   for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)
-      fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);
+      fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
   fprintf(ficresvij,"\n");
 
   xp=vector(1,npar);
@@ -2734,7 +3039,7 @@ void varevsij(char optionfilefiname[], d
 }  /* end varevsij */
 
 /************ Variance of prevlim ******************/
-void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)
+void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[])
 {
   /* Variance of prevalence limit */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
@@ -2747,8 +3052,9 @@ void varprevlim(char fileres[], double *
   double **gradg, **trgradg;
   double age,agelim;
   int theta;
-   
-  fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
+  
+  pstamp(ficresvpl);
+  fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
   fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);
@@ -2817,7 +3123,7 @@ void varprevlim(char fileres[], double *
 }
 
 /************ Variance of one-step probabilities  ******************/
-void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
+void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[])
 {
   int i, j=0,  i1, k1, l1, t, tj;
   int k2, l2, j1,  z1;
@@ -2862,13 +3168,15 @@ void varprob(char optionfilefiname[], do
   fprintf(ficlog,"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(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
-  
+  pstamp(ficresprob);
   fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");
+  pstamp(ficresprobcov);
   fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
   fprintf(ficresprobcov,"# Age");
+  pstamp(ficresprobcor);
   fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
-  fprintf(ficresprobcov,"# Age");
+  fprintf(ficresprobcor,"# Age");
 
 
   for(i=1; i<=nlstate;i++)
@@ -3119,6 +3427,8 @@ To be simple, these graphs help to under
   }
   free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
   free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
+  free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+  free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
   free_vector(xp,1,npar);
   fclose(ficresprob);
   fclose(ficresprobcov);
@@ -3137,20 +3447,24 @@ void printinghtml(char fileres[], char t
 		  double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;
 
-   fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \
+   fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
+   <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
+</ul>");
+   fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
 	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
    fprintf(fichtm,"\
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
 	   stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
    fprintf(fichtm,"\
- - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
+ - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
 	   subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
    fprintf(fichtm,"\
- - Life expectancies by age and initial health status (estepm=%2d months): \
+ - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age:  ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",
 	   estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
 
+
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
  m=cptcoveff;
@@ -3167,20 +3481,20 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
        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: %s%d1.png<br> \
-<img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);     
+     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d1.png\">%s%d1.png</a><br> \
+<img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),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: %s%d2.png<br> \
-<img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); 
-       /* Stable prevalence in each health state */
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d2.png\">%s%d2.png</a><br> \
+<img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); 
+       /* Period (stable) prevalence in each health state */
        for(cpt=1; cpt<nlstate;cpt++){
-	 fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br> \
-<img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+	 fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \
+<img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
        }
      for(cpt=1; cpt<=nlstate;cpt++) {
-        fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \
-<img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+        fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
+<img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
      }
    } /* end i1 */
  }/* End k1 */
@@ -3188,7 +3502,7 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
 
 
  fprintf(fichtm,"\
-\n<br><li><h4> Result files (second order: variances)</h4>\n\
+\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);
 
  fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
@@ -3201,13 +3515,21 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
 	 subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
  fprintf(fichtm,"\
- - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
+ - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
+   <a href=\"%s\">%s</a> <br>\n</li>",
+	   estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve"));
+ fprintf(fichtm,"\
+ - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
+   <a href=\"%s\">%s</a> <br>\n</li>",
+	   estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
+ fprintf(fichtm,"\
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
 	 estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
  fprintf(fichtm,"\
- - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n",
+ - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n",
 	 subdirf2(fileres,"t"),subdirf2(fileres,"t"));
  fprintf(fichtm,"\
- - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\
+ - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
 	 subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
 
 /*  if(popforecast==1) fprintf(fichtm,"\n */
@@ -3279,7 +3601,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");
      }
-     fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+     fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");
@@ -3327,7 +3649,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
-      k=2+nlstate*(2*cpt-2);
+      /*       k=2+nlstate*(2*cpt-2); */
+      k=2+(nlstate+1)*(cpt-1);
       fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
       fprintf(ficgp,"set ter png small\n\
 set size 0.65,0.65\n\
@@ -3341,9 +3664,11 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	
       */
       for (i=1; i< nlstate ; i ++) {
-	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1);
+	/*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
 	
       } 
+      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt);
     }
   }
   
@@ -3902,28 +4227,28 @@ double gompertz(double x[])
 { 
   double A,B,L=0.0,sump=0.,num=0.;
   int i,n=0; /* n is the size of the sample */
+
   for (i=0;i<=imx-1 ; i++) {
     sump=sump+weight[i];
-    sump=sump+1;
+    /*    sump=sump+1;*/
     num=num+1;
   }
  
  
-  /* for (i=1; i<=imx; i++) 
+  /* for (i=0; i<=imx; i++) 
      if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
 
-  for (i=0;i<=imx-1 ; i++)
+  for (i=1;i<=imx ; i++)
     {
-      if (cens[i]==1 & wav[i]>1)
-	A=-x[1]/(x[2])*
-	  (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)));
+      if (cens[i] == 1 && wav[i]>1)
+	A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
       
-      if (cens[i]==0 & wav[i]>1)
-	A=-x[1]/(x[2])*
-             (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)))
-	  +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM);      
+      if (cens[i] == 0 && wav[i]>1)
+	A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
+	     +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);  
       
-      if (wav[i]>1 & agecens[i]>15) {
+      /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
+      if (wav[i] > 1 ) { /* ??? */
 	L=L+A*weight[i];
 	/* 	printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
       }
@@ -3937,8 +4262,8 @@ double gompertz(double x[])
 /******************* Printing html file ***********/
 void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
-		  int imx,  double p[],double **matcov){
-  int i;
+		  int imx,  double p[],double **matcov,double agemortsup){
+  int i,k;
 
   fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");
   fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);
@@ -3946,6 +4271,15 @@ void printinghtmlmort(char fileres[], ch
     fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
   fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");
   fprintf(fichtm,"</ul>");
+
+fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>");
+
+ fprintf(fichtm,"\nAge   l<inf>x</inf>     q<inf>x</inf> d(x,x+1)    L<inf>x</inf>     T<inf>x</inf>     e<infx</inf><br>");
+
+ for (k=agegomp;k<(agemortsup-2);k++) 
+   fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf<br>\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
+
+ 
   fflush(fichtm);
 }
 
@@ -3975,6 +4309,7 @@ void printinggnuplotmort(char fileres[],
 
 
 
+
 /***********************************************/
 /**************** Main Program *****************/
 /***********************************************/
@@ -3983,14 +4318,17 @@ int main(int argc, char *argv[])
 {
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
+  int linei, month, year,iout;
   int jj, ll, li, lj, lk, imk;
   int numlinepar=0; /* Current linenumber of parameter file */
   int itimes;
   int NDIM=2;
 
   char ca[32], cb[32], cc[32];
+  char dummy[]="                         ";
   /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */
+  struct stat info;
   double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
 
@@ -4003,7 +4341,8 @@ int main(int argc, char *argv[])
   int *indx;
   char line[MAXLINE], linepar[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
-  char pathr[MAXLINE]; 
+  char pathr[MAXLINE], pathimach[MAXLINE]; 
+  char **bp, *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
   int sdeb, sfin; /* Status at beginning and end */
   int c,  h , cpt,l;
@@ -4013,6 +4352,8 @@ int main(int argc, char *argv[])
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
   int hstepm, nhstepm;
+  int agemortsup;
+  float  sumlpop=0.;
   double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
   double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
 
@@ -4037,7 +4378,7 @@ int main(int argc, char *argv[])
   char z[1]="c", occ;
 
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
-  char strstart[80], *strt, strtend[80];
+  char  *strt, strtend[80];
   char *stratrunc;
   int lstra;
 
@@ -4079,7 +4420,17 @@ int main(int argc, char *argv[])
   printf("\n%s\n%s",version,fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
-    scanf("%s",pathtot);
+    fgets(pathr,FILENAMELENGTH,stdin);
+    i=strlen(pathr);
+    if(pathr[i-1]=='\n')
+      pathr[i-1]='\0';
+   for (tok = pathr; tok != NULL; ){
+      printf("Pathr |%s|\n",pathr);
+      while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
+      printf("val= |%s| pathr=%s\n",val,pathr);
+      strcpy (pathtot, val);
+      if(pathr[0] == '\0') break; /* Dirty */
+    }
   }
   else{
     strcpy(pathtot,argv[1]);
@@ -4089,9 +4440,17 @@ int main(int argc, char *argv[])
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
   /* cutv(path,optionfile,pathtot,'\\');*/
 
+  /* Split argv[0], imach program to get pathimach */
+  printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]);
+  split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
+  printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
+ /*   strcpy(pathimach,argv[0]); */
+  /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);
-  printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
-  chdir(path);
+  printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+  chdir(path); /* Can be a relative path */
+  if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */
+    printf("Current directory %s!\n",pathcd);
   strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);
   if((outcmd=system(command)) != 0){
@@ -4115,12 +4474,12 @@ int main(int argc, char *argv[])
   }
   fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s\n%s",version,fullversion);
-  fprintf(ficlog,"\nEnter the parameter file name: ");
-  fprintf(ficlog,"pathtot=%s\n\
+  fprintf(ficlog,"\nEnter the parameter file name: \n");
+  fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\
  optionfile=%s\n\
  optionfilext=%s\n\
- optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+ optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -4257,7 +4616,6 @@ int main(int argc, char *argv[])
     }  
     fflush(ficlog);
 
-
     p=param[1][1];
     
     /* Reads comments: lines beginning with '#' */
@@ -4358,8 +4716,8 @@ int main(int argc, char *argv[])
 
   /*-------- data file ----------*/
   if((fic=fopen(datafile,"r"))==NULL)    {
-    printf("Problem with datafile: %s\n", datafile);goto end;
-    fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;
+    printf("Problem while opening datafile: %s\n", datafile);goto end;
+    fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end;
   }
 
   n= lastobs;
@@ -4381,40 +4739,117 @@ int main(int argc, char *argv[])
   ncodemax=ivector(1,8);
 
   i=1;
-  while (fgets(line, MAXLINE, fic) != NULL)    {
-    if ((i >= firstobs) && (i <=lastobs)) {
-	
-      for (j=maxwav;j>=1;j--){
-	cutv(stra, strb,line,' '); s[j][i]=atoi(strb); 
-	strcpy(line,stra);
-	cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
-	cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
+  linei=0;
+  while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
+    linei=linei+1;
+    for(j=strlen(line); j>=0;j--){  /* Untabifies line */
+      if(line[j] == '\t')
+	line[j] = ' ';
+    }
+    for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
+      ;
+    };
+    line[j+1]=0;  /* Trims blanks at end of line */
+    if(line[0]=='#'){
+      fprintf(ficlog,"Comment line\n%s\n",line);
+      printf("Comment line\n%s\n",line);
+      continue;
+    }
+
+    for (j=maxwav;j>=1;j--){
+      cutv(stra, strb,line,' '); 
+      errno=0;
+      lval=strtol(strb,&endptr,10); 
+      /*	if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+      if( strb[0]=='\0' || (*endptr != '\0')){
+	printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+	exit(1);
       }
-	
-      cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);
-      cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);
-
-      cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);
-      cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
-
-      cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
-      for (j=ncovcol;j>=1;j--){
-	cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
-      } 
-      lstra=strlen(stra);
-      if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
-	stratrunc = &(stra[lstra-9]);
-	num[i]=atol(stratrunc);
+      s[j][i]=lval;
+      
+      strcpy(line,stra);
+      cutv(stra, strb,line,' ');
+      if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+      }
+      else  if(iout=sscanf(strb,"%s.") != 0){
+	month=99;
+	year=9999;
+      }else{
+	printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
+	exit(1);
+      }
+      anint[j][i]= (double) year; 
+      mint[j][i]= (double)month; 
+      strcpy(line,stra);
+    } /* ENd Waves */
+    
+    cutv(stra, strb,line,' '); 
+    if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+    }
+    else  if(iout=sscanf(strb,"%s.",dummy) != 0){
+      month=99;
+      year=9999;
+    }else{
+      printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
+      exit(1);
+    }
+    andc[i]=(double) year; 
+    moisdc[i]=(double) month; 
+    strcpy(line,stra);
+    
+    cutv(stra, strb,line,' '); 
+    if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+    }
+    else  if(iout=sscanf(strb,"%s.") != 0){
+      month=99;
+      year=9999;
+    }else{
+      printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line,j);
+      exit(1);
+    }
+    annais[i]=(double)(year);
+    moisnais[i]=(double)(month); 
+    strcpy(line,stra);
+    
+    cutv(stra, strb,line,' '); 
+    errno=0;
+    lval=strtol(strb,&endptr,10); 
+    if( strb[0]=='\0' || (*endptr != '\0')){
+      printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight.  Exiting.\n",lval, i,line,linei);
+      exit(1);
+    }
+    weight[i]=(double)(lval); 
+    strcpy(line,stra);
+    
+    for (j=ncovcol;j>=1;j--){
+      cutv(stra, strb,line,' '); 
+      errno=0;
+      lval=strtol(strb,&endptr,10); 
+      if( strb[0]=='\0' || (*endptr != '\0')){
+	printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1).  Exiting.\n",lval, linei,i, line);
+	exit(1);
+      }
+      if(lval <-1 || lval >1){
+	printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self).  Exiting.\n",lval,linei, i,line,j);
+	exit(1);
       }
-      else
-	num[i]=atol(stra);
-	
-      /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
-	printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
-
-      i=i+1;
-    }
-  }
+      covar[j][i]=(double)(lval);
+      strcpy(line,stra);
+    } 
+    lstra=strlen(stra);
+    
+    if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
+      stratrunc = &(stra[lstra-9]);
+      num[i]=atol(stratrunc);
+    }
+    else
+      num[i]=atol(stra);
+    /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+      printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
+    
+    i=i+1;
+  } /* End loop reading  data */
+  fclose(fic);
   /* printf("ii=%d", ij);
      scanf("%d",i);*/
   imx=i-1; /* Number of individuals */
@@ -4428,12 +4863,12 @@ int main(int argc, char *argv[])
      if (s[4][i]==9)  s[4][i]=-1; 
      printf("%ld %.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]));}*/
   
- for (i=1; i<=imx; i++)
+  /* for (i=1; i<=imx; i++) */
  
    /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;
      else weight[i]=1;*/
 
-  /* Calculation of the number of parameter from char model*/
+  /* Calculation of the number of parameters from char model */
   Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
   Tprod=ivector(1,15); 
   Tvaraff=ivector(1,15); 
@@ -4512,8 +4947,7 @@ int main(int argc, char *argv[])
   printf("cptcovprod=%d ", cptcovprod);
   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
 
-  scanf("%d ",i);
-  fclose(fic);*/
+  scanf("%d ",i);*/
 
     /*  if(mle==1){*/
   if (weightopt != 1) { /* Maximisation without weights*/
@@ -4546,7 +4980,7 @@ int main(int argc, char *argv[])
   for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
     for(m=firstpass; (m<= lastpass); m++){
-      if(s[m][i] >0){
+      if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
 	if (s[m][i] >= nlstate+1) {
 	  if(agedc[i]>0)
 	    if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
@@ -4562,8 +4996,7 @@ int main(int argc, char *argv[])
 	    }
 	}
 	else if(s[m][i] !=9){ /* Standard case, age in fractional
-				 years but with the precision of a
-				 month */
+				 years but with the precision of a month */
 	  agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
 	  if((int)mint[m][i]==99 || (int)anint[m][i]==9999)
 	    agev[m][i]=1;
@@ -4705,11 +5138,13 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 \n\
 <hr  size=\"2\" color=\"#EC5E5E\">\
  <ul><li><h4>Parameter files</h4>\n\
+ - Parameter file: <a href=\"%s.%s\">%s.%s</a><br>\n\
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\
  - Log file of the run: <a href=\"%s\">%s</a><br>\n\
  - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
  - Date and time at start: %s</ul>\n",\
 	  fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+	  optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
 	  fileres,fileres,\
 	  filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
   fflush(fichtm);
@@ -4720,7 +5155,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   
   /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);
 
   fprintf(fichtm,"\n");
   fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
@@ -4739,6 +5174,7 @@ Interval (in months) between two waves:
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */
 
   globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
+
   if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);
     cens=ivector(1,n);
@@ -4748,9 +5184,9 @@ Interval (in months) between two waves:
  
     for (i=1; i<=imx; i++){
       dcwave[i]=-1;
-      for (j=1; j<=lastpass; j++)
-	if (s[j][i]>nlstate) {
-	  dcwave[i]=j;
+      for (m=firstpass; m<=lastpass; m++)
+	if (s[m][i]>nlstate) {
+	  dcwave[i]=m;
 	  /*	printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
 	  break;
 	}
@@ -4759,12 +5195,16 @@ Interval (in months) between two waves:
     for (i=1; i<=imx; i++) {
       if (wav[i]>0){
 	ageexmed[i]=agev[mw[1][i]][i];
-	j=wav[i];agecens[i]=1.; 
-	if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];
-	cens[i]=1;
-	
-	if (ageexmed[i]<1) cens[i]=-1;
-	if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;
+	j=wav[i];
+	agecens[i]=1.; 
+
+	if (ageexmed[i]> 1 && wav[i] > 0){
+	  agecens[i]=agev[mw[j][i]][i];
+	  cens[i]= 1;
+	}else if (ageexmed[i]< 1) 
+	  cens[i]= -1;
+	if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+	  cens[i]=0 ;
       }
       else cens[i]=-1;
     }
@@ -4773,29 +5213,29 @@ Interval (in months) between two waves:
       for (j=1;j<=NDIM;j++)
 	ximort[i][j]=(i == j ? 1.0 : 0.0);
     }
-
-    p[1]=0.1; p[2]=0.1;
+    
+    p[1]=0.0268; p[NDIM]=0.083;
     /*printf("%lf %lf", p[1], p[2]);*/
     
     
-  printf("Powell\n");  fprintf(ficlog,"Powell\n");
-  strcpy(filerespow,"pow-mort"); 
-  strcat(filerespow,fileres);
-  if((ficrespow=fopen(filerespow,"w"))==NULL) {
-    printf("Problem with resultfile: %s\n", filerespow);
-    fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
-  }
-  fprintf(ficrespow,"# Powell\n# iter -2*LL");
-  /*  for (i=1;i<=nlstate;i++)
-    for(j=1;j<=nlstate+ndeath;j++)
-      if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
-  */
-  fprintf(ficrespow,"\n");
-
+    printf("Powell\n");  fprintf(ficlog,"Powell\n");
+    strcpy(filerespow,"pow-mort"); 
+    strcat(filerespow,fileres);
+    if((ficrespow=fopen(filerespow,"w"))==NULL) {
+      printf("Problem with resultfile: %s\n", filerespow);
+      fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
+    }
+    fprintf(ficrespow,"# Powell\n# iter -2*LL");
+    /*  for (i=1;i<=nlstate;i++)
+	for(j=1;j<=nlstate+ndeath;j++)
+	if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
+    */
+    fprintf(ficrespow,"\n");
+    
     powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
     fclose(ficrespow);
     
-    hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); 
+    hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); 
 
     for(i=1; i <=NDIM; i++)
       for(j=i+1;j<=NDIM;j++)
@@ -4812,14 +5252,49 @@ Interval (in months) between two waves:
     printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
     for (i=1;i<=NDIM;i++) 
       printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
-    replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
+
+    lsurv=vector(1,AGESUP);
+    lpop=vector(1,AGESUP);
+    tpop=vector(1,AGESUP);
+    lsurv[agegomp]=100000;
+    
+    for (k=agegomp;k<=AGESUP;k++) {
+      agemortsup=k;
+      if (p[1]*exp(p[2]*(k-agegomp))>1) break;
+    }
+    
+    for (k=agegomp;k<agemortsup;k++)
+      lsurv[k+1]=lsurv[k]-lsurv[k]*(p[1]*exp(p[2]*(k-agegomp)));
+    
+    for (k=agegomp;k<agemortsup;k++){
+      lpop[k]=(lsurv[k]+lsurv[k+1])/2.;
+      sumlpop=sumlpop+lpop[k];
+    }
+    
+    tpop[agegomp]=sumlpop;
+    for (k=agegomp;k<(agemortsup-3);k++){
+      /*  tpop[k+1]=2;*/
+      tpop[k+1]=tpop[k]-lpop[k];
+    }
+    
+    
+    printf("\nAge   lx     qx    dx    Lx     Tx     e(x)\n");
+    for (k=agegomp;k<(agemortsup-2);k++) 
+      printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
+    
+    
+    replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
-		     model,imx,p,matcov);
+		     model,imx,p,matcov,agemortsup);
+    
+    free_vector(lsurv,1,AGESUP);
+    free_vector(lpop,1,AGESUP);
+    free_vector(tpop,1,AGESUP);
   } /* Endof if mle==-3 */
-
+  
   else{ /* For mle >=1 */
   
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
@@ -5045,7 +5520,7 @@ Interval (in months) between two waves:
     /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/
     /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
     
-    replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
+    replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
     printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
@@ -5067,17 +5542,18 @@ Interval (in months) between two waves:
     fclose(ficres);
 
 
-    /*--------------- Prevalence limit  (stable prevalence) --------------*/
+    /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
   
     strcpy(filerespl,"pl");
     strcat(filerespl,fileres);
     if((ficrespl=fopen(filerespl,"w"))==NULL) {
-      printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
-      fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
+      printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
+      fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
     }
-    printf("Computing stable prevalence: result on file '%s' \n", filerespl);
-    fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
-    fprintf(ficrespl,"#Stable prevalence \n");
+    printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+    fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+    pstamp(ficrespl);
+    fprintf(ficrespl,"# Period (stable) prevalence \n");
     fprintf(ficrespl,"#Age ");
     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
     fprintf(ficrespl,"\n");
@@ -5137,7 +5613,7 @@ Interval (in months) between two waves:
     hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
 
     /* hstepm=1;   aff par mois*/
-
+    pstamp(ficrespij);
     fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
@@ -5174,7 +5650,7 @@ Interval (in months) between two waves:
       }
     }
 
-    varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);
+    varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
 
     fclose(ficrespij);
 
@@ -5207,8 +5683,8 @@ Interval (in months) between two waves:
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;
       fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
     }
-    printf("Computing Total LEs with variances: file '%s' \n", filerest); 
-    fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); 
+    printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
+    fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
 
 
     strcpy(filerese,"e");
@@ -5220,6 +5696,24 @@ Interval (in months) between two waves:
     printf("Computing Health Expectancies: result on file '%s' \n", filerese);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
 
+    strcpy(fileresstde,"stde");
+    strcat(fileresstde,fileres);
+    if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
+      printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+      fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+    }
+    printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+    fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+
+    strcpy(filerescve,"cve");
+    strcat(filerescve,fileres);
+    if((ficrescveij=fopen(filerescve,"w"))==NULL) {
+      printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
+      fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
+    }
+    printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+    fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+
     strcpy(fileresv,"v");
     strcat(fileresv,fileres);
     if((ficresvij=fopen(fileresv,"w"))==NULL) {
@@ -5252,9 +5746,16 @@ Interval (in months) between two waves:
 	fprintf(ficrest,"******\n");
 
 	fprintf(ficreseij,"\n#****** ");
-	for(j=1;j<=cptcoveff;j++) 
+	fprintf(ficresstdeij,"\n#****** ");
+	fprintf(ficrescveij,"\n#****** ");
+	for(j=1;j<=cptcoveff;j++) {
 	  fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	}
 	fprintf(ficreseij,"******\n");
+	fprintf(ficresstdeij,"******\n");
+	fprintf(ficrescveij,"******\n");
 
 	fprintf(ficresvij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
@@ -5263,17 +5764,18 @@ Interval (in months) between two waves:
 
 	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, delti, matcov);  
+	evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
+	cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);  
  
 	vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
 	oldm=oldms;savm=savms;
-	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);
+	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart);
 	if(popbased==1){
-	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav);
+	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
 	}
 
- 
-	fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
+	pstamp(ficrest);
+	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) ");
 	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	fprintf(ficrest,"\n");
 
@@ -5321,19 +5823,21 @@ Interval (in months) between two waves:
     free_ivector(cod,1,n);
     free_ivector(tab,1,NCOVMAX);
     fclose(ficreseij);
+    fclose(ficresstdeij);
+    fclose(ficrescveij);
     fclose(ficresvij);
     fclose(ficrest);
     fclose(ficpar);
   
-    /*------- Variance of stable prevalence------*/   
+    /*------- Variance of period (stable) prevalence------*/   
 
     strcpy(fileresvpl,"vpl");
     strcat(fileresvpl,fileres);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
-      printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);
+      printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
       exit(0);
     }
-    printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
+    printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
 
     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
@@ -5345,7 +5849,7 @@ Interval (in months) between two waves:
       
 	varpl=matrix(1,nlstate,(int) bage, (int) fage);
 	oldm=oldms;savm=savms;
-	varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
+	varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
 	free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
       }
     }
@@ -5357,11 +5861,11 @@ Interval (in months) between two waves:
     free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   }  /* mle==-3 arrives here for freeing */
+  free_matrix(prlim,1,nlstate,1,nlstate);
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
-  
     free_matrix(covar,0,NCOVMAX,1,n);
     free_matrix(matcov,1,npar,1,npar);
     /*free_vector(delti,1,npar);*/
@@ -5376,7 +5880,8 @@ Interval (in months) between two waves:
     free_ivector(Tage,1,15);
     free_ivector(Tcode,1,100);
 
-
+    free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
+    free_imatrix(codtab,1,100,1,10);
   fflush(fichtm);
   fflush(ficgp);
   
@@ -5394,7 +5899,7 @@ Interval (in months) between two waves:
   tm = *localtime(&end_time.tv_sec);
   tmg = *gmtime(&end_time.tv_sec);
   strcpy(strtend,asctime(&tm));
-  printf("Local time at start %s\nLocaltime at end   %s",strstart, strtend); 
+  printf("Local time at start %s\nLocal time at end   %s",strstart, strtend); 
   fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend); 
   printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
 
@@ -5410,13 +5915,38 @@ Interval (in months) between two waves:
   fclose(ficlog);
   /*------ End -----------*/
 
-  chdir(path);
-  strcpy(plotcmd,GNUPLOTPROGRAM);
-  strcat(plotcmd," ");
-  strcat(plotcmd,optionfilegnuplot);
-  printf("Starting graphs with: %s",plotcmd);fflush(stdout);
+
+   printf("Before Current directory %s!\n",pathcd);
+   if(chdir(pathcd) != 0)
+    printf("Can't move to directory %s!\n",path);
+  if(getcwd(pathcd,MAXLINE) > 0)
+    printf("Current directory %s!\n",pathcd);
+  /*strcat(plotcmd,CHARSEPARATOR);*/
+  sprintf(plotcmd,"gnuplot");
+#ifndef UNIX
+  sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
+#endif
+  if(!stat(plotcmd,&info)){
+    printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+    if(!stat(getenv("GNUPLOTBIN"),&info)){
+      printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
+    }else
+      strcpy(pplotcmd,plotcmd);
+#ifdef UNIX
+    strcpy(plotcmd,GNUPLOTPROGRAM);
+    if(!stat(plotcmd,&info)){
+      printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+    }else
+      strcpy(pplotcmd,plotcmd);
+#endif
+  }else
+    strcpy(pplotcmd,plotcmd);
+  
+  sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
+  printf("Starting graphs with: %s\n",plotcmd);fflush(stdout);
+
   if((outcmd=system(plotcmd)) != 0){
-    printf(" Problem with gnuplot\n");
+    printf("\n Problem with gnuplot\n");
   }
   printf(" Wait...");
   while (z[0] != 'q') {
@@ -5424,7 +5954,10 @@ Interval (in months) between two waves:
     printf("\nType e to edit output files, g to graph again and q for exiting: ");
     scanf("%s",z);
 /*     if (z[0] == 'c') system("./imach"); */
-    if (z[0] == 'e') system(optionfilehtm);
+    if (z[0] == 'e') {
+      printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
+      system(optionfilehtm);
+    }
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }