--- imach/src/imach.c	2002/07/26 12:29:55	1.58
+++ imach/src/imach.c	2002/11/18 23:01:13	1.59
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.58 2002/07/26 12:29:55 lievre Exp $
+/* $Id: imach.c,v 1.59 2002/11/18 23:01:13 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -83,7 +83,7 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES ";
+char version[80]="Imach version 0.9, November 2002, INED-EUROREVES ";
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -99,6 +99,8 @@ int jmin, jmax; /* min, max spacing betw
 int mle, weightopt;
 int **mw; /* mw[mi][i] is number of the mi wave for this individual */
 int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
+int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
+	   * wave mi and wave mi+1 is not an exact multiple of stepm. */
 double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
@@ -177,49 +179,49 @@ double ftolhess; /* Tolerance for comput
 /**************** split *************************/
 static	int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {
-   char	*s;				/* pointer */
-   int	l1, l2;				/* length counters */
+  char	*ss;				/* pointer */
+  int	l1, l2;				/* length counters */
 
-   l1 = strlen( path );			/* length of path */
-   if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
-   s= strrchr( path, DIRSEPARATOR );		/* find last / */
-   if ( s == NULL ) {			/* no directory, so use current */
-     /*if(strrchr(path, ODIRSEPARATOR )==NULL)
-       printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
+  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(strrchr(path, ODIRSEPARATOR )==NULL)
+      printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
 #if	defined(__bsd__)		/* get current working directory */
-      extern char	*getwd( );
+    extern char	*getwd( );
 
-      if ( getwd( dirc ) == NULL ) {
+    if ( getwd( dirc ) == NULL ) {
 #else
-      extern char	*getcwd( );
+    extern char	*getcwd( );
 
-      if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
+    if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
 #endif
-         return( GLOCK_ERROR_GETCWD );
-      }
-      strcpy( name, path );		/* we've got it */
-   } else {				/* strip direcotry from path */
-      s++;				/* after this, the filename */
-      l2 = strlen( s );			/* length of filename */
-      if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
-      strcpy( name, s );		/* save file name */
-      strncpy( dirc, path, l1 - l2 );	/* now the directory */
-      dirc[l1-l2] = 0;			/* add zero */
-   }
-   l1 = strlen( dirc );			/* length of directory */
+      return( GLOCK_ERROR_GETCWD );
+    }
+    strcpy( name, path );		/* we've got it */
+  } else {				/* strip direcotry from path */
+    ss++;				/* after this, the filename */
+    l2 = strlen( ss );			/* length of filename */
+    if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
+    strcpy( name, ss );		/* save file name */
+    strncpy( dirc, path, l1 - l2 );	/* now the directory */
+    dirc[l1-l2] = 0;			/* add zero */
+  }
+  l1 = strlen( dirc );			/* length of directory */
 #ifdef windows
-   if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
+  if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
 #else
-   if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
+  if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
 #endif
-   s = strrchr( name, '.' );		/* find last / */
-   s++;
-   strcpy(ext,s);			/* save extension */
-   l1= strlen( name);
-   l2= strlen( s)+1;
-   strncpy( finame, name, l1-l2);
-   finame[l1-l2]= 0;
-   return( 0 );				/* we're done */
+  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;
+  return( 0 );				/* we're done */
 }
 
 
@@ -277,7 +279,7 @@ void nrerror(char error_text[])
 {
   fprintf(stderr,"ERREUR ...\n");
   fprintf(stderr,"%s\n",error_text);
-  exit(1);
+  exit(EXIT_FAILURE);
 }
 /*********************** vector *******************/
 double *vector(int nl, int nh)
@@ -914,6 +916,8 @@ double func( double *x)
   double **out;
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
+  int s1, s2;
+  double bbh;
   long ipmx;
   /*extern weight */
   /* We are differentiating ll according to initial status */
@@ -928,7 +932,10 @@ double func( double *x)
     for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
     for(mi=1; mi<= wav[i]-1; mi++){
       for (ii=1;ii<=nlstate+ndeath;ii++)
-	for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	for (j=1;j<=nlstate+ndeath;j++){
+	  oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	  savm[ii][j]=(ii==j ? 1.0 : 0.0);
+	}
       for(d=0; d<dh[mi][i]; d++){
 	newm=savm;
 	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
@@ -944,8 +951,27 @@ double func( double *x)
 	
       } /* end mult */
       
-      lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
-      /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
+      /*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.
+       * 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
+       * 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
+       * -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. 
+       */
+      s1=s[mw[mi][i]][i];
+      s2=s[mw[mi+1][i]][i];
+      bbh=(double)bh[mi][i]/(double)stepm;
+      lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));
+      /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/
+      /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+      /*if(lli ==000.0)*/
+      /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
@@ -1237,7 +1263,7 @@ void lubksb(double **a, int n, int *indx
 } 
 
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
+void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
 {  /* Some frequencies */
   
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;
@@ -1482,12 +1508,12 @@ void prevalence(int agemin, float agemax
 
 /************* Waves Concatenation ***************/
 
-void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
+void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
 {
   /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
      Death is a valid wave (if date is known).
      mw[mi][i] is the mi (mi=1 to wav[i])  effective wave of individual i
-     dh[m][i] of dh[mw[mi][i][i] is the delay between two effective waves m=mw[mi][i]
+     dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
      and mw[mi+1][i]. dh depends on stepm.
      */
 
@@ -1558,12 +1584,20 @@ void  concatwav(int wav[], int **dh, int
 	jk= j/stepm;
 	jl= j -jk*stepm;
 	ju= j -(jk+1)*stepm;
-	if(jl <= -ju)
+	if(jl <= -ju){
 	  dh[mi][i]=jk;
-	else
+	  bh[mi][i]=jl;
+	}
+	else{
 	  dh[mi][i]=jk+1;
-	if(dh[mi][i]==0)
+	  bh[mi][i]=ju;
+	}
+	if(dh[mi][i]==0){
 	  dh[mi][i]=1; /* At least one step */
+	  bh[mi][i]=ju; /* At least one step */
+	  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
+	}
+	if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
       }
     }
   }
@@ -2102,7 +2136,7 @@ void varevsij(char optionfilefiname[], d
 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)
 {
   /* Variance of prevalence limit */
-  /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
+  /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
   double **newm;
   double **dnewm,**doldm;
   int i, j, nhstepm, hstepm;
@@ -2274,7 +2308,6 @@ void varprob(char optionfilefiname[], do
 
   }
 
- 
   cov[1]=1;
   tj=cptcoveff;
   if (cptcovn<1) {tj=1;ncodemax[1]=1;}
@@ -2282,7 +2315,6 @@ void varprob(char optionfilefiname[], do
   for(t=1; t<=tj;t++){
     for(i1=1; i1<=ncodemax[t];i1++){ 
       j1++;
-      
       if  (cptcovn>0) {
 	fprintf(ficresprob, "\n#********** Variable "); 
 	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
@@ -2355,7 +2387,11 @@ void varprob(char optionfilefiname[], do
 	
 	matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); 
 	matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
-	
+	free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+	free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+	free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+	free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+
 	pmij(pmmij,cov,ncovmodel,x,nlstate);
 	
 	k=0;
@@ -2370,10 +2406,10 @@ void varprob(char optionfilefiname[], do
 	    varpij[i][j][(int)age] = doldm[i][j];
 
 	/*printf("\n%d ",(int)age);
-     for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
-       printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
-       fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
-     }*/
+	  for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+	  printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+	  fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
+	  }*/
 
 	fprintf(ficresprob,"\n%d ",(int)age);
 	fprintf(ficresprobcov,"\n%d ",(int)age);
@@ -2401,13 +2437,13 @@ void varprob(char optionfilefiname[], do
 
       /* Confidence intervalle of pij  */
       /*
-      fprintf(ficgp,"\nset noparametric;unset label");
-      fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
-      fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
-      fprintf(fichtm,"\n<br>Probability with  confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
-      fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
-      fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
-      fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
+	fprintf(ficgp,"\nset noparametric;unset label");
+	fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
+	fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+	fprintf(fichtm,"\n<br>Probability with  confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
+	fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
+	fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
+	fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
       */
 
       /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
@@ -2479,13 +2515,9 @@ void varprob(char optionfilefiname[], do
 	} /*l1 */
       }/* k1 */
     } /* loop covariates */
-    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
-    free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
-    free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
-    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
-    free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
-    free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   }
+  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_vector(xp,1,npar);
   fclose(ficresprob);
   fclose(ficresprobcov);
@@ -2929,7 +2961,7 @@ populforecast(char fileres[], double anp
   
   int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
   int *popage;
-  double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+  double calagedate, agelim, kk1, kk2;
   double *popeffectif,*popcount;
   double ***p3mat,***tabpop,***tabpopprev;
   double ***mobaverage;
@@ -3109,7 +3141,7 @@ int main(int argc, char *argv[])
   int c,  h , cpt,l;
   int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
-  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; 
+  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab; 
   int mobilav=0,popforecast=0;
   int hstepm, nhstepm;
   double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
@@ -3128,7 +3160,7 @@ int main(int argc, char *argv[])
   double *epj, vepp;
   double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
-  
+  /*int *movingaverage; */
 
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
@@ -3139,9 +3171,9 @@ int main(int argc, char *argv[])
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
  
   /* long total_usecs;
-  struct timeval start_time, end_time;
+     struct timeval start_time, end_time;
   
-  gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
+     gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   getcwd(pathcd, size);
 
   printf("\n%s",version);
@@ -3158,11 +3190,11 @@ int main(int argc, char *argv[])
   /* cutv(path,optionfile,pathtot,'\\');*/
 
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);
-   printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+  printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);
   replace(pathc,path);
 
-/*-------- arguments in the command line --------*/
+  /*-------- arguments in the command line --------*/
 
   /* Log file */
   strcat(filelog, optionfilefiname);
@@ -3210,7 +3242,7 @@ int main(int argc, char *argv[])
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-while((c=getc(ficpar))=='#' && c!= EOF){
+  while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);
     puts(line);
@@ -3237,7 +3269,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   ungetc(c,ficpar);
   
   param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
-    for(i=1; i <=nlstate; i++)
+  for(i=1; i <=nlstate; i++)
     for(j=1; j <=nlstate+ndeath-1; j++){
       fscanf(ficpar,"%1d%1d",&i1,&j1);
       fprintf(ficparo,"%1d%1d",i1,j1);
@@ -3261,7 +3293,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
       fprintf(ficparo,"\n");
     }
   
-    npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
+  npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
 
   p=param[1][1];
   
@@ -3334,73 +3366,72 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   fprintf(ficlog,"\n");
 
 
-    /*-------- Rewriting paramater file ----------*/
-     strcpy(rfileres,"r");    /* "Rparameterfile */
-     strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
-     strcat(rfileres,".");    /* */
-     strcat(rfileres,optionfilext);    /* Other files have txt extension */
-    if((ficres =fopen(rfileres,"w"))==NULL) {
-      printf("Problem writing new parameter file: %s\n", fileres);goto end;
-      fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
-    }
-    fprintf(ficres,"#%s\n",version);
+  /*-------- Rewriting paramater file ----------*/
+  strcpy(rfileres,"r");    /* "Rparameterfile */
+  strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
+  strcat(rfileres,".");    /* */
+  strcat(rfileres,optionfilext);    /* Other files have txt extension */
+  if((ficres =fopen(rfileres,"w"))==NULL) {
+    printf("Problem writing new parameter file: %s\n", fileres);goto end;
+    fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
+  }
+  fprintf(ficres,"#%s\n",version);
     
-    /*-------- 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;
-    }
-
-    n= lastobs;
-    severity = vector(1,maxwav);
-    outcome=imatrix(1,maxwav+1,1,n);
-    num=ivector(1,n);
-    moisnais=vector(1,n);
-    annais=vector(1,n);
-    moisdc=vector(1,n);
-    andc=vector(1,n);
-    agedc=vector(1,n);
-    cod=ivector(1,n);
-    weight=vector(1,n);
-    for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
-    mint=matrix(1,maxwav,1,n);
-    anint=matrix(1,maxwav,1,n);
-    s=imatrix(1,maxwav+1,1,n);
-    adl=imatrix(1,maxwav+1,1,n);    
-    tab=ivector(1,NCOVMAX);
-    ncodemax=ivector(1,8);
-
-    i=1;
-    while (fgets(line, MAXLINE, fic) != NULL)    {
-      if ((i >= firstobs) && (i <=lastobs)) {
+  /*-------- 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;
+  }
+
+  n= lastobs;
+  severity = vector(1,maxwav);
+  outcome=imatrix(1,maxwav+1,1,n);
+  num=ivector(1,n);
+  moisnais=vector(1,n);
+  annais=vector(1,n);
+  moisdc=vector(1,n);
+  andc=vector(1,n);
+  agedc=vector(1,n);
+  cod=ivector(1,n);
+  weight=vector(1,n);
+  for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
+  mint=matrix(1,maxwav,1,n);
+  anint=matrix(1,maxwav,1,n);
+  s=imatrix(1,maxwav+1,1,n);
+  tab=ivector(1,NCOVMAX);
+  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);
-	}
+      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);
+      }
 	
-	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,'/'); 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,'/'); 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);
-	} 
-	num[i]=atol(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);
+      } 
+      num[i]=atol(stra);
 	
-	/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
-	  printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
+      /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+	printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
 
-	i=i+1;
-      }
-    } 
-    /* printf("ii=%d", ij);
-       scanf("%d",i);*/
+      i=i+1;
+    }
+  }
+  /* printf("ii=%d", ij);
+     scanf("%d",i);*/
   imx=i-1; /* Number of individuals */
 
   /* for (i=1; i<=imx; i++){
@@ -3434,11 +3465,11 @@ while((c=getc(ficpar))=='#' && c!= EOF){
       goto end;
     }
     
-    /* This loop fill the array Tvar from the string 'model'.*/
+    /* This loop fills the array Tvar from the string 'model'.*/
 
     for(i=(j+1); i>=1;i--){
       cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ 
-      if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */
+      if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
       /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
       /*scanf("%d",i);*/
       if (strchr(strb,'*')) {  /* Model includes a product */
@@ -3496,269 +3527,265 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   fclose(fic);*/
 
     /*  if(mle==1){*/
-    if (weightopt != 1) { /* Maximisation without weights*/
-      for(i=1;i<=n;i++) weight[i]=1.0;
-    }
+  if (weightopt != 1) { /* Maximisation without weights*/
+    for(i=1;i<=n;i++) weight[i]=1.0;
+  }
     /*-calculation of age at interview from date of interview and age at death -*/
-    agev=matrix(1,maxwav,1,imx);
+  agev=matrix(1,maxwav,1,imx);
 
-    for (i=1; i<=imx; i++) {
-      for(m=2; (m<= maxwav); m++) {
-       if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
-	 anint[m][i]=9999;
-	 s[m][i]=-1;
-       }
-     if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
-      }
+  for (i=1; i<=imx; i++) {
+    for(m=2; (m<= maxwav); m++) {
+      if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
+	anint[m][i]=9999;
+	s[m][i]=-1;
+      }
+      if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
     }
+  }
 
-    for (i=1; i<=imx; i++)  {
-      agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
-      for(m=1; (m<= maxwav); m++){
-	if(s[m][i] >0){
-	  if (s[m][i] >= nlstate+1) {
-	    if(agedc[i]>0)
-	      if(moisdc[i]!=99 && andc[i]!=9999)
-		agev[m][i]=agedc[i];
-	    /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
-	   else {
+  for (i=1; i<=imx; i++)  {
+    agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
+    for(m=1; (m<= maxwav); m++){
+      if(s[m][i] >0){
+	if (s[m][i] >= nlstate+1) {
+	  if(agedc[i]>0)
+	    if(moisdc[i]!=99 && andc[i]!=9999)
+	      agev[m][i]=agedc[i];
+	  /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
+	    else {
 	      if (andc[i]!=9999){
-	      printf("Warning negative age at death: %d line:%d\n",num[i],i);
-	      fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
-	      agev[m][i]=-1;
+		printf("Warning negative age at death: %d line:%d\n",num[i],i);
+		fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
+		agev[m][i]=-1;
 	      }
 	    }
-	  }
-	  else if(s[m][i] !=9){ /* Should no more exist */
-	    agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
-	    if(mint[m][i]==99 || anint[m][i]==9999)
-	      agev[m][i]=1;
-	    else if(agev[m][i] <agemin){ 
-	      agemin=agev[m][i];
-	      /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
-	    }
-	    else if(agev[m][i] >agemax){
-	      agemax=agev[m][i];
-	     /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
-	    }
-	    /*agev[m][i]=anint[m][i]-annais[i];*/
-	    /*	 agev[m][i] = age[i]+2*m;*/
-	  }
-	  else { /* =9 */
+	}
+	else if(s[m][i] !=9){ /* Should no more exist */
+	  agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
+	  if(mint[m][i]==99 || anint[m][i]==9999)
 	    agev[m][i]=1;
-	    s[m][i]=-1;
+	  else if(agev[m][i] <agemin){ 
+	    agemin=agev[m][i];
+	    /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
+	  }
+	  else if(agev[m][i] >agemax){
+	    agemax=agev[m][i];
+	    /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
 	  }
+	  /*agev[m][i]=anint[m][i]-annais[i];*/
+	  /*	 agev[m][i] = age[i]+2*m;*/
 	}
-	else /*= 0 Unknown */
+	else { /* =9 */
 	  agev[m][i]=1;
+	  s[m][i]=-1;
+	}
       }
-    
+      else /*= 0 Unknown */
+	agev[m][i]=1;
     }
-    for (i=1; i<=imx; i++)  {
-      for(m=1; (m<= maxwav); m++){
-	if (s[m][i] > (nlstate+ndeath)) {
-	  printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	
-	  fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	
-	  goto end;
-	}
+    
+  }
+  for (i=1; i<=imx; i++)  {
+    for(m=1; (m<= maxwav); m++){
+      if (s[m][i] > (nlstate+ndeath)) {
+	printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	
+	fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	
+	goto end;
       }
     }
+  }
 
-    printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
-    fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); 
-
-    free_vector(severity,1,maxwav);
-    free_imatrix(outcome,1,maxwav+1,1,n);
-    free_vector(moisnais,1,n);
-    free_vector(annais,1,n);
-    /* free_matrix(mint,1,maxwav,1,n);
-       free_matrix(anint,1,maxwav,1,n);*/
-    free_vector(moisdc,1,n);
-    free_vector(andc,1,n);
+  printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
+  fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); 
+
+  free_vector(severity,1,maxwav);
+  free_imatrix(outcome,1,maxwav+1,1,n);
+  free_vector(moisnais,1,n);
+  free_vector(annais,1,n);
+  /* free_matrix(mint,1,maxwav,1,n);
+     free_matrix(anint,1,maxwav,1,n);*/
+  free_vector(moisdc,1,n);
+  free_vector(andc,1,n);
 
    
-    wav=ivector(1,imx);
-    dh=imatrix(1,lastpass-firstpass+1,1,imx);
-    mw=imatrix(1,lastpass-firstpass+1,1,imx);
+  wav=ivector(1,imx);
+  dh=imatrix(1,lastpass-firstpass+1,1,imx);
+  bh=imatrix(1,lastpass-firstpass+1,1,imx);
+  mw=imatrix(1,lastpass-firstpass+1,1,imx);
    
-    /* Concatenates waves */
-      concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
+  /* Concatenates waves */
+  concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
 
-      /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
+  /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
 
-      Tcode=ivector(1,100);
-      nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
-      ncodemax[1]=1;
-      if (cptcovn > 0) tricode(Tvar,nbcode,imx);
+  Tcode=ivector(1,100);
+  nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
+  ncodemax[1]=1;
+  if (cptcovn > 0) tricode(Tvar,nbcode,imx);
       
-      codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of 
-				     the estimations*/
-   h=0;
-   m=pow(2,cptcoveff);
+  codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of 
+				 the estimations*/
+  h=0;
+  m=pow(2,cptcoveff);
  
-   for(k=1;k<=cptcoveff; k++){
-     for(i=1; i <=(m/pow(2,k));i++){
-       for(j=1; j <= ncodemax[k]; j++){
-	 for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
-	   h++;
-	   if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
-	   /*  printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
-	 } 
-       }
-     }
-   } 
-   /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); 
-      codtab[1][2]=1;codtab[2][2]=2; */
-   /* for(i=1; i <=m ;i++){ 
-      for(k=1; k <=cptcovn; k++){
-      printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
-      }
-      printf("\n");
+  for(k=1;k<=cptcoveff; k++){
+    for(i=1; i <=(m/pow(2,k));i++){
+      for(j=1; j <= ncodemax[k]; j++){
+	for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
+	  h++;
+	  if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
+	  /*  printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
+	} 
       }
-      scanf("%d",i);*/
+    }
+  } 
+  /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); 
+     codtab[1][2]=1;codtab[2][2]=2; */
+  /* for(i=1; i <=m ;i++){ 
+     for(k=1; k <=cptcovn; k++){
+     printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
+     }
+     printf("\n");
+     }
+     scanf("%d",i);*/
     
-   /* Calculates basic frequencies. Computes observed prevalence at single age
-       and prints on file fileres'p'. */
+  /* Calculates basic frequencies. Computes observed prevalence at single age
+     and prints on file fileres'p'. */
 
     
    
-    pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
-    oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
-    newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
-    savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
-    oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
-     
-    /* For Powell, parameters are in a vector p[] starting at p[1]
-       so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
-    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
+  /* For Powell, parameters are in a vector p[] starting at p[1]
+     so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
+  p=param[1][1]; /* *(*(*(param +1)+1)+0) */
 
-    if(mle==1){
+  if(mle==1){
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
-    }
+  }
     
-    /*--------- results files --------------*/
-    fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
+  /*--------- results files --------------*/
+  fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
   
 
-   jk=1;
-   fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
-   printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
-   fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
-   for(i=1,jk=1; i <=nlstate; i++){
-     for(k=1; k <=(nlstate+ndeath); k++){
-       if (k != i) 
-	 {
-	   printf("%d%d ",i,k);
-	   fprintf(ficlog,"%d%d ",i,k);
-	   fprintf(ficres,"%1d%1d ",i,k);
-	   for(j=1; j <=ncovmodel; j++){
-	     printf("%f ",p[jk]);
-	     fprintf(ficlog,"%f ",p[jk]);
-	     fprintf(ficres,"%f ",p[jk]);
-	     jk++; 
-	   }
-	   printf("\n");
-	   fprintf(ficlog,"\n");
-	   fprintf(ficres,"\n");
-	 }
-     }
-   }
-   if(mle==1){
-     /* Computing hessian and covariance matrix */
-     ftolhess=ftol; /* Usually correct */
-     hesscov(matcov, p, npar, delti, ftolhess, func);
-   }
-   fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
-   printf("# Scales (for hessian or gradient estimation)\n");
-   fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
-   for(i=1,jk=1; i <=nlstate; i++){
-     for(j=1; j <=nlstate+ndeath; j++){
-       if (j!=i) {
-	 fprintf(ficres,"%1d%1d",i,j);
-	 printf("%1d%1d",i,j);
-	 fprintf(ficlog,"%1d%1d",i,j);
-	 for(k=1; k<=ncovmodel;k++){
-	   printf(" %.5e",delti[jk]);
-	   fprintf(ficlog," %.5e",delti[jk]);
-	   fprintf(ficres," %.5e",delti[jk]);
-	   jk++;
-	 }
-	 printf("\n");
-	 fprintf(ficlog,"\n");
-	 fprintf(ficres,"\n");
-       }
-     }
-   }
+  jk=1;
+  fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+  printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+  fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+  for(i=1,jk=1; i <=nlstate; i++){
+    for(k=1; k <=(nlstate+ndeath); k++){
+      if (k != i) 
+	{
+	  printf("%d%d ",i,k);
+	  fprintf(ficlog,"%d%d ",i,k);
+	  fprintf(ficres,"%1d%1d ",i,k);
+	  for(j=1; j <=ncovmodel; j++){
+	    printf("%f ",p[jk]);
+	    fprintf(ficlog,"%f ",p[jk]);
+	    fprintf(ficres,"%f ",p[jk]);
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	  fprintf(ficres,"\n");
+	}
+    }
+  }
+  if(mle==1){
+    /* Computing hessian and covariance matrix */
+    ftolhess=ftol; /* Usually correct */
+    hesscov(matcov, p, npar, delti, ftolhess, func);
+  }
+  fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
+  printf("# Scales (for hessian or gradient estimation)\n");
+  fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
+  for(i=1,jk=1; i <=nlstate; i++){
+    for(j=1; j <=nlstate+ndeath; j++){
+      if (j!=i) {
+	fprintf(ficres,"%1d%1d",i,j);
+	printf("%1d%1d",i,j);
+	fprintf(ficlog,"%1d%1d",i,j);
+	for(k=1; k<=ncovmodel;k++){
+	  printf(" %.5e",delti[jk]);
+	  fprintf(ficlog," %.5e",delti[jk]);
+	  fprintf(ficres," %.5e",delti[jk]);
+	  jk++;
+	}
+	printf("\n");
+	fprintf(ficlog,"\n");
+	fprintf(ficres,"\n");
+      }
+    }
+  }
    
-   k=1;
-   fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
-   if(mle==1)
-     printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
-   fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
-   for(i=1;i<=npar;i++){
-     /*  if (k>nlstate) k=1;
-	 i1=(i-1)/(ncovmodel*nlstate)+1; 
-	 fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
-	 printf("%s%d%d",alph[k],i1,tab[i]);*/
-     fprintf(ficres,"%3d",i);
-     if(mle==1)
-       printf("%3d",i);
-     fprintf(ficlog,"%3d",i);
-     for(j=1; j<=i;j++){
-       fprintf(ficres," %.5e",matcov[i][j]);
-       if(mle==1)
-	 printf(" %.5e",matcov[i][j]);
-       fprintf(ficlog," %.5e",matcov[i][j]);
-     }
-     fprintf(ficres,"\n");
-     if(mle==1)
-       printf("\n");
-     fprintf(ficlog,"\n");
-     k++;
-   }
+  fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
+  if(mle==1)
+    printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
+  fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
+  for(i=1,k=1;i<=npar;i++){
+    /*  if (k>nlstate) k=1;
+	i1=(i-1)/(ncovmodel*nlstate)+1; 
+	fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
+	printf("%s%d%d",alph[k],i1,tab[i]);
+    */
+    fprintf(ficres,"%3d",i);
+    if(mle==1)
+      printf("%3d",i);
+    fprintf(ficlog,"%3d",i);
+    for(j=1; j<=i;j++){
+      fprintf(ficres," %.5e",matcov[i][j]);
+      if(mle==1)
+	printf(" %.5e",matcov[i][j]);
+      fprintf(ficlog," %.5e",matcov[i][j]);
+    }
+    fprintf(ficres,"\n");
+    if(mle==1)
+      printf("\n");
+    fprintf(ficlog,"\n");
+    k++;
+  }
    
-   while((c=getc(ficpar))=='#' && c!= EOF){
-     ungetc(c,ficpar);
-     fgets(line, MAXLINE, ficpar);
-     puts(line);
-     fputs(line,ficparo);
-   }
-   ungetc(c,ficpar);
-   estepm=0;
-   fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
-   if (estepm==0 || estepm < stepm) estepm=stepm;
-   if (fage <= 2) {
-     bage = ageminpar;
-     fage = agemaxpar;
-   }
+  while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
+
+  estepm=0;
+  fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
+  if (estepm==0 || estepm < stepm) estepm=stepm;
+  if (fage <= 2) {
+    bage = ageminpar;
+    fage = agemaxpar;
+  }
    
-   fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
-   fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
-   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+  fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
+  fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+  fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
    
-   while((c=getc(ficpar))=='#' && c!= EOF){
-     ungetc(c,ficpar);
-     fgets(line, MAXLINE, ficpar);
-     puts(line);
-     fputs(line,ficparo);
-   }
-   ungetc(c,ficpar);
+  while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
   
-   fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
-   fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-   fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+  fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
+  fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+  fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
    
-   while((c=getc(ficpar))=='#' && c!= EOF){
-     ungetc(c,ficpar);
-     fgets(line, MAXLINE, ficpar);
-     puts(line);
-     fputs(line,ficparo);
-   }
-   ungetc(c,ficpar);
+  while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
  
 
-   dateprev1=anprev1+mprev1/12.+jprev1/365.;
-   dateprev2=anprev2+mprev2/12.+jprev2/365.;
+  dateprev1=anprev1+mprev1/12.+jprev1/365.;
+  dateprev2=anprev2+mprev2/12.+jprev2/365.;
 
   fscanf(ficpar,"pop_based=%d\n",&popbased);
   fprintf(ficparo,"pop_based=%d\n",popbased);   
@@ -3773,11 +3800,11 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   ungetc(c,ficpar);
 
   fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
-fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
-fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
+  fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
+  fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
 
 
-while((c=getc(ficpar))=='#' && c!= EOF){
+  while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);
     puts(line);
@@ -3789,22 +3816,22 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
 
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
 
-/*------------ gnuplot -------------*/
- strcpy(optionfilegnuplot,optionfilefiname);
- strcat(optionfilegnuplot,".gp");
- if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
-   printf("Problem with file %s",optionfilegnuplot);
- }
- else{
-   fprintf(ficgp,"\n# %s\n", version); 
-   fprintf(ficgp,"# %s\n", optionfilegnuplot); 
-   fprintf(ficgp,"set missing 'NaNq'\n");
-}
- fclose(ficgp);
- printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
-/*--------- index.htm --------*/
+  /*------------ gnuplot -------------*/
+  strcpy(optionfilegnuplot,optionfilefiname);
+  strcat(optionfilegnuplot,".gp");
+  if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+    printf("Problem with file %s",optionfilegnuplot);
+  }
+  else{
+    fprintf(ficgp,"\n# %s\n", version); 
+    fprintf(ficgp,"# %s\n", optionfilegnuplot); 
+    fprintf(ficgp,"set missing 'NaNq'\n");
+  }
+  fclose(ficgp);
+  printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
+  /*--------- index.htm --------*/
 
   strcpy(optionfilehtm,optionfile);
   strcat(optionfilehtm,".htm");
@@ -3824,19 +3851,21 @@ Interval (in months) between two waves:
  - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
   fclose(fichtm);
 
- printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
+  printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
  
-/*------------ free_vector  -------------*/
- chdir(path);
+  /*------------ 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);
+  free_ivector(wav,1,imx);
+  free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
+  free_imatrix(bh,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,0,NCOVMAX,1,n);
+  /*free_matrix(covar,1,NCOVMAX,1,n);*/
+  fclose(ficparo);
+  fclose(ficres);
 
 
   /*--------------- Prevalence limit  (stable prevalence) --------------*/
@@ -3860,38 +3889,38 @@ Interval (in months) between two waves:
   newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
   savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
   oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
-  k=0;
+
   agebase=ageminpar;
   agelim=agemaxpar;
   ftolpl=1.e-10;
   i1=cptcoveff;
   if (cptcovn < 1){i1=1;}
 
-  for(cptcov=1;cptcov<=i1;cptcov++){
+  for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-	k=k+1;
-	/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
-	fprintf(ficrespl,"\n#******");
-	printf("\n#******");
-	fprintf(ficlog,"\n#******");
-	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	}
-	fprintf(ficrespl,"******\n");
-	printf("******\n");
-	fprintf(ficlog,"******\n");
+      k=k+1;
+      /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
+      fprintf(ficrespl,"\n#******");
+      printf("\n#******");
+      fprintf(ficlog,"\n#******");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+      }
+      fprintf(ficrespl,"******\n");
+      printf("******\n");
+      fprintf(ficlog,"******\n");
 	
-	for (age=agebase; age<=agelim; age++){
-	  prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
-	  fprintf(ficrespl,"%.0f",age );
-	  for(i=1; i<=nlstate;i++)
+      for (age=agebase; age<=agelim; age++){
+	prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+	fprintf(ficrespl,"%.0f",age );
+	for(i=1; i<=nlstate;i++)
 	  fprintf(ficrespl," %.5f", prlim[i][i]);
-	  fprintf(ficrespl,"\n");
-	}
+	fprintf(ficrespl,"\n");
       }
     }
+  }
   fclose(ficrespl);
 
   /*------------- h Pij x at various ages ------------*/
@@ -3913,39 +3942,38 @@ Interval (in months) between two waves:
 
   /* hstepm=1;   aff par mois*/
 
-  k=0;
-  for(cptcov=1;cptcov<=i1;cptcov++){
+  for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;
-	fprintf(ficrespij,"\n#****** ");
-	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	fprintf(ficrespij,"******\n");
+      fprintf(ficrespij,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++) 
+	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+      fprintf(ficrespij,"******\n");
 	
-	for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
-	  nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
-	  nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
-
-	  /*	  nhstepm=nhstepm*YEARM; aff par mois*/
-
-	  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);  
-	  fprintf(ficrespij,"# Age");
+      for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
+	nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+
+	/*	  nhstepm=nhstepm*YEARM; aff par mois*/
+
+	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);  
+	fprintf(ficrespij,"# Age");
+	for(i=1; i<=nlstate;i++)
+	  for(j=1; j<=nlstate+ndeath;j++)
+	    fprintf(ficrespij," %1d-%1d",i,j);
+	fprintf(ficrespij,"\n");
+	for (h=0; h<=nhstepm; h++){
+	  fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
 	  for(i=1; i<=nlstate;i++)
 	    for(j=1; j<=nlstate+ndeath;j++)
-	      fprintf(ficrespij," %1d-%1d",i,j);
-	  fprintf(ficrespij,"\n");
-	   for (h=0; h<=nhstepm; h++){
-	    fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
-	    for(i=1; i<=nlstate;i++)
-	      for(j=1; j<=nlstate+ndeath;j++)
-		fprintf(ficrespij," %.5f", p3mat[i][j][h]);
-	    fprintf(ficrespij,"\n");
-	     }
-	  free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	      fprintf(ficrespij," %.5f", p3mat[i][j][h]);
 	  fprintf(ficrespij,"\n");
 	}
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	fprintf(ficrespij,"\n");
+      }
     }
   }
 
@@ -4008,8 +4036,7 @@ Interval (in months) between two waves:
     }
   }
 
-  k=0;
-  for(cptcov=1;cptcov<=i1;cptcov++){
+  for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1; 
       fprintf(ficrest,"\n#****** ");
@@ -4036,7 +4063,7 @@ Interval (in months) between two waves:
       varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);
       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);
-       }
+      }
 
  
       fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
@@ -4074,16 +4101,22 @@ Interval (in months) between two waves:
 	}
 	fprintf(ficrest,"\n");
       }
+      free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+      free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+      free_vector(epj,1,nlstate+1);
     }
   }
-free_matrix(mint,1,maxwav,1,n);
-    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
-    free_vector(weight,1,n);
+  free_vector(weight,1,n);
+  free_imatrix(Tvard,1,15,1,2);
+  free_imatrix(s,1,maxwav+1,1,n);
+  free_matrix(anint,1,maxwav,1,n); 
+  free_matrix(mint,1,maxwav,1,n);
+  free_ivector(cod,1,n);
+  free_ivector(tab,1,NCOVMAX);
   fclose(ficreseij);
   fclose(ficresvij);
   fclose(ficrest);
   fclose(ficpar);
-  free_vector(epj,1,nlstate+1);
   
   /*------- Variance of stable prevalence------*/   
 
@@ -4095,8 +4128,7 @@ free_matrix(mint,1,maxwav,1,n);
   }
   printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
 
-  k=0;
-  for(cptcov=1;cptcov<=i1;cptcov++){
+  for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;
       fprintf(ficresvpl,"\n#****** ");
@@ -4106,19 +4138,14 @@ free_matrix(mint,1,maxwav,1,n);
       
       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);
+      free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
     }
- }
+  }
 
   fclose(ficresvpl);
 
   /*---------- End : free ----------------*/
-  free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
-  
-  free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-  free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-  
-  
   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);
@@ -4129,6 +4156,12 @@ free_matrix(mint,1,maxwav,1,n);
   free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  free_ivector(ncodemax,1,8);
+  free_ivector(Tvar,1,15);
+  free_ivector(Tprod,1,15);
+  free_ivector(Tvaraff,1,15);
+  free_ivector(Tage,1,15);
+  free_ivector(Tcode,1,100);
 
   fprintf(fichtm,"\n</body>");
   fclose(fichtm);
@@ -4150,8 +4183,7 @@ free_matrix(mint,1,maxwav,1,n);
   /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/
 
-
- end:
+  end:
 #ifdef windows
   /* chdir(pathcd);*/
 #endif 
@@ -4159,11 +4191,11 @@ free_matrix(mint,1,maxwav,1,n);
  /*system("../gp37mgw/wgnuplot graph.plt");*/
  /*system("cd ../gp37mgw");*/
  /* system("..\\gp37mgw\\wgnuplot graph.plt");*/
- strcpy(plotcmd,GNUPLOTPROGRAM);
- strcat(plotcmd," ");
- strcat(plotcmd,optionfilegnuplot);
- printf("Starting: %s\n",plotcmd);fflush(stdout);
- system(plotcmd);
+  strcpy(plotcmd,GNUPLOTPROGRAM);
+  strcat(plotcmd," ");
+  strcat(plotcmd,optionfilegnuplot);
+  printf("Starting: %s\n",plotcmd);fflush(stdout);
+  system(plotcmd);
 
  /*#ifdef windows*/
   while (z[0] != 'q') {