--- imach/src/imach.c	2015/08/18 16:28:39	1.195
+++ imach/src/imach.c	2015/09/22 19:45:16	1.202
@@ -1,6 +1,32 @@
-/* $Id: imach.c,v 1.195 2015/08/18 16:28:39 brouard Exp $
+/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.202  2015/09/22 19:45:16  brouard
+  Summary: Adding some overall graph on contribution to likelihood. Might change
+
+  Revision 1.201  2015/09/15 17:34:58  brouard
+  Summary: 0.98r0
+
+  - Some new graphs like suvival functions
+  - Some bugs fixed like model=1+age+V2.
+
+  Revision 1.200  2015/09/09 16:53:55  brouard
+  Summary: Big bug thanks to Flavia
+
+  Even model=1+age+V2. did not work anymore
+
+  Revision 1.199  2015/09/07 14:09:23  brouard
+  Summary: 0.98q6 changing default small png format for graph to vectorized svg.
+
+  Revision 1.198  2015/09/03 07:14:39  brouard
+  Summary: 0.98q5 Flavia
+
+  Revision 1.197  2015/09/01 18:24:39  brouard
+  *** empty log message ***
+
+  Revision 1.196  2015/08/18 23:17:52  brouard
+  Summary: 0.98q5
+
   Revision 1.195  2015/08/18 16:28:39  brouard
   Summary: Adding a hack for testing purpose
 
@@ -621,6 +647,7 @@
 
 /* #define DEBUG */
 /* #define DEBUGBRENT */
+#define DEBUGLINMIN
 #define POWELL /* Instead of NLOPT */
 #define POWELLF1F3 /* Skip test */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
@@ -691,7 +718,7 @@ typedef struct {
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
-#define codtabm(h,k)  1 & (h-1) >> (k-1) ;
+#define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
@@ -708,11 +735,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.195 2015/08/18 16:28:39 brouard Exp $ */
+/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ */
 /* $State: Exp $ */
-
-char version[]="Imach version 0.98q5, August 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.195 $ $Date: 2015/08/18 16:28:39 $"; 
+#include "version.h"
+char version[]=__IMACH_VERSION__;
+char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char fullversion[]="$Revision: 1.202 $ $Date: 2015/09/22 19:45:16 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -778,7 +806,7 @@ char command[FILENAMELENGTH];
 int  outcmd=0;
 
 char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
-
+char fileresu[FILENAMELENGTH]; /* fileres without r in front */
 char filelog[FILENAMELENGTH]; /* Log file */
 char filerest[FILENAMELENGTH];
 char fileregp[FILENAMELENGTH];
@@ -846,7 +874,7 @@ int estepm;
 
 int m,nb;
 long *num;
-int firstpass=0, lastpass=4,*cod, *Tage,*cens;
+int firstpass=0, lastpass=4,*cod, *cens;
 int *ncodemax;  /* ncodemax[j]= Number of modalities of the j th
 		   covariate for which somebody answered excluding 
 		   undefined. Usually 2: 0 and 1. */
@@ -866,8 +894,9 @@ double  **covar; /**< covar[j,i], value
 		  * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
 double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+int *Tage;
 int *Ndum; /** Freq of modality (tricode */
-int **codtab; /**< codtab=imatrix(1,100,1,10); */
+/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard, *Tprod, cptcovprod, *Tvaraff;
 double *lsurv, *lpop, *tpop;
 
@@ -1581,7 +1610,7 @@ void linmin(double p[], double xi[], int
   nrfunc=func; 
   for (j=1;j<=n;j++) { 
     pcom[j]=p[j]; 
-    xicom[j]=xi[j]; 
+    xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */
   } 
 
   /* axs=0.0; */
@@ -1605,6 +1634,7 @@ void linmin(double p[], double xi[], int
 
 #ifdef DEBUGLINMIN
   printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
+  fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
 #endif
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
   /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
@@ -1617,6 +1647,7 @@ void linmin(double p[], double xi[], int
 #endif
 #ifdef DEBUGLINMIN
   printf("linmin end ");
+  fprintf(ficlog,"linmin end ");
 #endif
   for (j=1;j<=n;j++) { 
     /* printf(" before xi[%d]=%12.8f", j,xi[j]); */
@@ -1628,10 +1659,14 @@ void linmin(double p[], double xi[], int
   /* printf("\n"); */
 #ifdef DEBUGLINMIN
   printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
+  fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
   for (j=1;j<=n;j++) { 
-    printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]);
-    if(j % ncovmodel == 0)
+    printf(" xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+    fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+    if(j % ncovmodel == 0){
       printf("\n");
+      fprintf(ficlog,"\n");
+    }
   }
 #endif
   free_vector(xicom,1,n); 
@@ -1666,7 +1701,7 @@ void powell(double p[], double **xi, int
   xits=vector(1,n); 
   *fret=(*func)(p); 
   for (j=1;j<=n;j++) pt[j]=p[j]; 
-    rcurr_time = time(NULL);  
+  rcurr_time = time(NULL);  
   for (*iter=1;;++(*iter)) { 
     fp=(*fret); /* From former iteration or initial value */
     ibig=0; 
@@ -1811,7 +1846,7 @@ void powell(double p[], double **xi, int
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
       t= t- del*SQR(fp-fptt);
 #endif
-      directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */
+      directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
 #ifdef DEBUG
       printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
       fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
@@ -1826,9 +1861,9 @@ void powell(double p[], double **xi, int
       if (t < 0.0) { /* Then we use it for new direction */
 #else
       if (directest*t < 0.0) { /* Contradiction between both tests */
-	printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
+	printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
         printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
-        fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
+        fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
         fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
       } 
       if (directest < 0.0) { /* Then we use it for new direction */
@@ -1836,17 +1871,23 @@ void powell(double p[], double **xi, int
 #ifdef DEBUGLINMIN
 	printf("Before linmin in direction P%d-P0\n",n);
 	for (j=1;j<=n;j++) { 
-	  printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-	  if(j % ncovmodel == 0)
+	  printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  if(j % ncovmodel == 0){
 	    printf("\n");
+	    fprintf(ficlog,"\n");
+	  }
 	}
 #endif
 	linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
 #ifdef DEBUGLINMIN
 	for (j=1;j<=n;j++) { 
 	  printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
-	  if(j % ncovmodel == 0)
+	  fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+	  if(j % ncovmodel == 0){
 	    printf("\n");
+	    fprintf(ficlog,"\n");
+	  }
 	}
 #endif
 	for (j=1;j<=n;j++) { 
@@ -1886,7 +1927,8 @@ double **prevalim(double **prlim, int nl
   /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij();
   double **newm;
-  double agefin, delaymax=50 ; /* Max number of years to converge */
+  double agefin, delaymax=100 ; /* Max number of years to converge */
+  long int ncvyear=0, ncvloop=0;
   
   for (ii=1;ii<=nlstate+ndeath;ii++)
     for (j=1;j<=nlstate+ndeath;j++){
@@ -1896,20 +1938,25 @@ double **prevalim(double **prlim, int nl
   cov[1]=1.;
   
   /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+  /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
+    ncvloop++;
     newm=savm;
     /* Covariates have to be included here again */
     cov[2]=agefin;
     if(nagesqr==1)
       cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {
-      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
-      /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
+      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+      /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
     }
     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
+    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
     for (k=1; k<=cptcovprod;k++) /* Useless */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
     
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -1928,17 +1975,21 @@ double **prevalim(double **prlim, int nl
 	sumnew=0;
 	for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
 	prlim[i][j]= newm[i][j]/(1-sumnew);
-        /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/
 	max=FMAX(max,prlim[i][j]);
 	min=FMIN(min,prlim[i][j]);
+        /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */
       }
       maxmin=max-min;
       maxmax=FMAX(maxmax,maxmin);
     } /* j loop */
     if(maxmax < ftolpl){
+      /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age-(int)agefin); */
       return prlim;
     }
   } /* age loop */
+  printf("Warning: the stable prevalence did not converge with the required precision ftolpl=6*10^5*ftol=%g. \n\
+Earliest age to start was %d-%d=%d, ncvloop=%ld, ncvyear=%d\n\
+Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin);
   return prlim; /* should not reach here */
 }
 
@@ -2082,12 +2133,15 @@ double ***hpxij(double ***po, int nhstep
       if(nagesqr==1)
 	cov[3]= agexact*agexact;
       for (k=1; k<=cptcovn;k++) 
-	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+	/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
 	/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+	/* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+	/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
 
 
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -2511,9 +2565,9 @@ 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,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
+	fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %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],
+		num[i], agexact, 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++){
 	  llt +=ll[k]*gipmx/gsw;
@@ -2545,14 +2599,14 @@ void likelione(FILE *ficres,double p[],
   int k;
 
   if(*globpri !=0){ /* Just counts and sums, no printings */
-    strcpy(fileresilk,"ilk"); 
-    strcat(fileresilk,fileres);
+    strcpy(fileresilk,"ILK_"); 
+    strcat(fileresilk,fileresu);
     if((ficresilk=fopen(fileresilk,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresilk);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
     }
     fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
-    fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight 2wlli out sav ");
+    fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav ");
     /* 	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; k<=nlstate; k++) 
       fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
@@ -2562,8 +2616,10 @@ void likelione(FILE *ficres,double p[],
   *fretone=(*funcone)(p);
   if(*globpri !=0){
     fclose(ficresilk);
-    fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
-    fflush(fichtm); 
+    fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
+    fprintf(fichtm,"<br>- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: <a href=\"%s.png\">%s.png</a><br> \
+<img src=\"%s.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+    fflush(fichtm);
   } 
   return;
 }
@@ -2595,7 +2651,7 @@ void mlikeli(FILE *ficres,double p[], in
     for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");
-  strcpy(filerespow,"pow"); 
+  strcpy(filerespow,"POW_"); 
   strcat(filerespow,fileres);
   if((ficrespow=fopen(filerespow,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", filerespow);
@@ -2906,7 +2962,7 @@ 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);
+  fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
 /************ Frequencies ********************/
@@ -2922,8 +2978,8 @@ void  freqsummary(char fileres[], int ia
   
   pp=vector(1,nlstate);
   prop=matrix(1,nlstate,iagemin,iagemax+3);
-  strcpy(fileresp,"p");
-  strcat(fileresp,fileres);
+  strcpy(fileresp,"P_");
+  strcat(fileresp,fileresu);
   if((ficresp=fopen(fileresp,"w"))==NULL) {
     printf("Problem with prevalence resultfile: %s\n", fileresp);
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
@@ -2958,13 +3014,13 @@ void  freqsummary(char fileres[], int ia
 	bool=1;
 	if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
 	  for (z1=1; z1<=cptcoveff; z1++)       
-            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){
+            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                 /* Tests if the value of each of the covariates of i is equal to filter j1 */
               bool=0;
-              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", 
-                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
-                j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/
-              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/
+              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
+                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+                j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
             } 
 	}
  
@@ -2993,10 +3049,10 @@ void  freqsummary(char fileres[], int ia
       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]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresp, "**********\n#");
 	fprintf(ficlog, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficlog, "**********\n#");
       }
       for(i=1; i<=nlstate;i++) 
@@ -3124,7 +3180,7 @@ void prevalence(double ***probs, double
 	bool=1;
 	if  (cptcovn>0) {
 	  for (z1=1; z1<=cptcoveff; z1++) 
-	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
+	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
 	      bool=0;
 	} 
 	if (bool==1) { 
@@ -3390,14 +3446,15 @@ void tricode(int *Tvar, int **nbcode, in
        nbcode[Tvar[j]][1]=0;
        nbcode[Tvar[j]][2]=1;
        nbcode[Tvar[j]][3]=2;
+       To be continued (not working yet).
     */
-    ij=0; /* ij is similar to i but can jumps over null modalities */
-    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/
-    	if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */
+    ij=0; /* ij is similar to i but can jump over null modalities */
+    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+    	if (Ndum[i] == 0) { /* If nobody responded to this modality k */
 	  break;
 	}
 	ij++;
-	nbcode[Tvar[j]][ij]=i;  /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
+	nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
 	cptcode = ij; /* New max modality for covar j */
     } /* end of loop on modality i=-1 to 1 or more */
       
@@ -3790,11 +3847,11 @@ void varevsij(char optionfilefiname[], d
 
   if(popbased==1){
     if(mobilav!=0)
-      strcpy(digitp,"-populbased-mobilav-");
-    else strcpy(digitp,"-populbased-nomobil-");
+      strcpy(digitp,"-POPULBASED-MOBILAV_");
+    else strcpy(digitp,"-POPULBASED-NOMOBIL_");
   }
   else 
-    strcpy(digitp,"-stablbased-");
+    strcpy(digitp,"-STABLBASED_");
 
   if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
@@ -3804,12 +3861,12 @@ void varevsij(char optionfilefiname[], d
     }
   }
 
-  strcpy(fileresprobmorprev,"prmorprev"); 
+  strcpy(fileresprobmorprev,"PRMORPREV-"); 
   sprintf(digit,"%-d",ij);
   /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
   strcat(fileresprobmorprev,digit); /* Tvar to be done */
   strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
-  strcat(fileresprobmorprev,fileres);
+  strcat(fileresprobmorprev,fileresu);
   if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobmorprev);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
@@ -3827,7 +3884,8 @@ void varevsij(char optionfilefiname[], d
   }  
   fprintf(ficresprobmorprev,"\n");
   fprintf(ficgp,"\n# Routine varevsij");
-  /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+  fprintf(ficgp,"\nunset title \n");
+/* 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);
 /*   } */
@@ -4039,9 +4097,11 @@ void varevsij(char optionfilefiname[], d
   free_vector(gmp,nlstate+1,nlstate+ndeath);
   free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
   free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-  fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240");
+  /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
+  fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
   /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
   fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+  fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
 /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
@@ -4049,11 +4109,11 @@ void varevsij(char optionfilefiname[], d
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
-  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
-  /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
+  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+  /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
 */
-/*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */
-  fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+/*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
+  fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
 
   free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,nlstate);
@@ -4171,20 +4231,20 @@ void varprob(char optionfilefiname[], do
   char fileresprobcor[FILENAMELENGTH];
   double ***varpij;
 
-  strcpy(fileresprob,"prob"); 
+  strcpy(fileresprob,"PROB_"); 
   strcat(fileresprob,fileres);
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
   }
-  strcpy(fileresprobcov,"probcov"); 
-  strcat(fileresprobcov,fileres);
+  strcpy(fileresprobcov,"PROBCOV_"); 
+  strcat(fileresprobcov,fileresu);
   if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcov);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
   }
-  strcpy(fileresprobcor,"probcor"); 
-  strcat(fileresprobcor,fileres);
+  strcpy(fileresprobcor,"PROBCOR_"); 
+  strcat(fileresprobcor,fileresu);
   if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcor);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
@@ -4226,10 +4286,9 @@ void varprob(char optionfilefiname[], do
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");
 
-  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\
-  file %s<br>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\
+  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
+  fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
  They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
   fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
@@ -4250,23 +4309,23 @@ To be simple, these graphs help to under
     /*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]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprob, "**********\n#\n");
 	fprintf(ficresprobcov, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcov, "**********\n#\n");
 	
 	fprintf(ficgp, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficgp, "**********\n#\n");
 	
 	
 	fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
 	
 	fprintf(ficresprobcor, "\n#********** Variable ");    
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcor, "**********\n#");    
       }
       
@@ -4279,7 +4338,8 @@ To be simple, these graphs help to under
 	if(nagesqr==1)
 	  cov[3]= age*age;
 	for (k=1; k<=cptcovn;k++) {
-	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
+	  /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
 							 * 1  1 1 1 1
 							 * 2  2 1 1 1
 							 * 3  1 2 1 1
@@ -4287,9 +4347,9 @@ To be simple, these graphs help to under
 	  /* nbcode[1][1]=0 nbcode[1][2]=1;*/
 	}
 	/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
 	for (k=1; k<=cptcovprod;k++)
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
 	
     
 	for(theta=1; theta <=npar; theta++){
@@ -4437,17 +4497,18 @@ To be simple, these graphs help to under
 		  /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
 		  if(first==1){
 		    first=0;
+ 		    fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
  		    fprintf(ficgp,"\nset parametric;unset label");
 		    fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
-		    fprintf(ficgp,"\nset ter png small size 320, 240");
+		    fprintf(ficgp,"\nset ter svg size 640, 480");
 		    fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
- :<a href=\"%s%d%1d%1d-%1d%1d.png\">\
-%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\
-			    subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\
-			    subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
-		    fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.png\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ :<a href=\"%s_%d%1d%1d-%1d%1d.svg\">\
+%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
+			    subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,\
+			    subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+		    fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 		    fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
-		    fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+		    fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 		    fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
 		    fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
 		    fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
@@ -4464,7 +4525,7 @@ To be simple, these graphs help to under
 		  }/* if first */
 		} /* age mod 5 */
 	      } /* end loop age */
-	      fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+	      fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 	      first=1;
 	    } /*l12 */
 	  } /* k12 */
@@ -4486,7 +4547,7 @@ To be simple, these graphs help to under
 
 
 /******************* Printing html file ***********/
-void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
 		  int popforecast, int estepm ,\
@@ -4499,20 +4560,20 @@ void printinghtml(char fileres[], char t
 </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"));
+	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
    fprintf(fichtm,"\
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
-	   stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
+	   stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
    fprintf(fichtm,"\
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
-	   subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
+	   subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
    fprintf(fichtm,"\
  - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates 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",
-	   estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+	   estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
    fprintf(fichtm,"\
  - Population projections by age and states: \
-   <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f"));
+   <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
 
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
@@ -4526,26 +4587,42 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
-	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
        }
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
+     /* aij, bij */
+     fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
      /* Pij */
-     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.png\">%s%d_1.png</a><br> \
-<img src=\"%s%d_1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);     
+     fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+<img src=\"%s_%d-2.svg\">",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: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \
-<img src=\"%s%d_2.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>- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
-<img src=\"%s%d_%d.png\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
-       }
+     fprintf(fichtm,"<br>\n- Iij 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,\
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
+divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
+     /* Survival functions (period) in state j */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);
+     }
+     /* State specific survival functions (period) */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Survival functions from state %d in any different live states and total.\
+ Or probability to survive in various states (1 to %d) being in state %d at different ages.\
+ <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
+     }
+     /* Period (stable) prevalence in each health state */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, 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>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
-<img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,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 in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s%d%d.svg\">%s%d%d.svg</a> <br> \
+<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
      }
    /* } /\* end i1 *\/ */
  }/* End k1 */
@@ -4554,34 +4631,41 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
  fprintf(fichtm,"\
 \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> \
- - 95%% confidence intervals and T statistics are in the log file.<br>\n", rfileres,rfileres);
+ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
+But because parameters are usually highly correlated (a higher incidence of disability \
+and a higher incidence of recovery can give very close observed transition) it might \
+be very useful to look not only at linear confidence intervals estimated from the \
+variances but at the covariance matrix. And instead of looking at the estimated coefficients \
+(parameters) of the logistic regression, it might be more meaningful to visualize the \
+covariance matrix of the one-step probabilities. \
+See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
 
  fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
+	 subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
  fprintf(fichtm,"\
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"));
+	 subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
 
  fprintf(fichtm,"\
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
+	 subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_"));
  fprintf(fichtm,"\
  - 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"));
+	   estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"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"));
+	   estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_"));
  fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> 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"));
+	 estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
  fprintf(fichtm,"\
  - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (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,"t"),subdirf2(fileres,"t"));
+	 estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
  fprintf(fichtm,"\
  - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
-	 subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
+	 subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
 
 /*  if(popforecast==1) fprintf(fichtm,"\n */
 /*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
@@ -4602,20 +4686,20 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++) 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png <br>\
-<img src=\"%s%d_%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);  
+prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\
+<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);  
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences: %s%d.png<br>\
-<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
+ observed and cahotic prevalences: %s_%d.svg<br>\
+<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
    /* } /\* end i1 *\/ */
  }/* End k1 */
  fprintf(fichtm,"</ul>");
@@ -4623,11 +4707,12 @@ true period expectancies (those weighted
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
   int ng=0;
+  int vpopbased;
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
 /*     fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
@@ -4638,76 +4723,97 @@ void printinggnuplot(char fileres[], cha
     /*#endif */
   m=pow(2,cptcoveff);
 
+  /* Contribution to likelihood */
+  /* Plot the probability implied in the likelihood */
+    fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
+    fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
+    /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
+    fprintf(ficgp,"\nset ter png size 640, 480");
+/* good for mle=4 plot by number of matrix products.
+   replot  "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
+/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */
+    /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
+    fprintf(ficgp,"\nset out \"%s.png\";",subdirf2(optionfilefiname,"ILK_"));
+    fprintf(ficgp,"\nplot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk));
+    fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk));
+    fprintf(ficgp,"\nset out\n");
+    /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
  /* 1eme*/
-  fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n");
+  fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
   for (cpt=1; cpt<= nlstate ; cpt ++) {
     for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-     fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
-     fprintf(ficgp,"\n#set out \"v%s%d_%d.png\" \n",optionfilefiname,cpt,k1);
+     fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+     fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
      fprintf(ficgp,"set xlabel \"Age\" \n\
 set ylabel \"Probability\" \n\
-set ter png small size 320, 240\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
+set ter svg size 640, 480\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
 
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else        fprintf(ficgp," %%*lf (%%*lf)");
      }
-     fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 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 lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," %%*lf (%%*lf)");
      } 
-     fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); 
+     fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," %%*lf (%%*lf)");
      }  
-     fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
-   }
-  }
+     fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+     fprintf(ficgp,"\nset out \n");
+    } /* k1 */
+  } /* cpt */
   /*2 eme*/
   fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
   for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1);
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage);
-    
-    for (i=1; i<= nlstate+1 ; i ++) {
-      k=2*i;
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
-      else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      fprintf(ficgp,"\" t\"\" w l lt 0,");
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-      else fprintf(ficgp,"\" t\"\" w l lt 0,");
-    }
-  }
-  
+    fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+    for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+      if(vpopbased==0)
+	fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+      else
+	fprintf(ficgp,"\nreplot ");
+      for (i=1; i<= nlstate+1 ; i ++) {
+	k=2*i;
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+	else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	fprintf(ficgp,"\" t\"\" w l lt 0,");
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+	else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+      } /* state */
+    } /* vpopbased */
+    fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+  } /* k1 */
   /*3eme*/
   
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       /*       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 size 320, 240\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt);
+      fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+      fprintf(ficgp,"set ter svg size 640, 480\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
 	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
 	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
@@ -4717,39 +4823,98 @@ 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+i,cpt,i+1);
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"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);
+      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
     }
   }
   
-  /* CV preval stable (period) */
+  /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       k=3;
+      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f]  ", ageminpar, agemaxpar);
+      for (i=1; i<= nlstate ; i ++){
+	if(i==1)
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	else
+	  fprintf(ficgp,", '' ");
+	l=(nlstate+ndeath)*(i-1)+1;
+	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+	for (j=2; j<= nlstate+ndeath ; j ++)
+	  fprintf(ficgp,"+$%d",k+l+j-1);
+	fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+      } /* nlstate */
+      fprintf(ficgp,"\nset out\n");
+    } /* end cpt state*/ 
+  } /* end covariate */  
+
+  /* Survival functions (period) from state i in state j by final state j */
+  for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+      k=3;
+      fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f]  ", ageminpar, agemaxpar);
+      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	if(j==1)
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	else
+	  fprintf(ficgp,", '' ");
+	l=(nlstate+ndeath)*(cpt-1) +j;
+	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+	/* for (i=2; i<= nlstate+ndeath ; i ++) */
+	/*   fprintf(ficgp,"+$%d",k+l+i-1); */
+	fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+      } /* nlstate */
+      fprintf(ficgp,", '' ");
+      fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	l=(nlstate+ndeath)*(cpt-1) +j;
+	if(j < nlstate)
+	  fprintf(ficgp,"$%d +",k+l);
+	else
+	  fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+      }
+      fprintf(ficgp,"\nset out\n");
+    } /* end cpt state*/ 
+  } /* end covariate */  
+
+  /* CV preval stable (period) for each covariate */
+  for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+      k=3;
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
-      fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter png small size 320, 240\n\
+set ter svg size 640, 480\n\
 unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);
       for (i=1; i<= nlstate ; i ++){
 	if(i==1)
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij"));
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
 	else
 	  fprintf(ficgp,", '' ");
 	l=(nlstate+ndeath)*(i-1)+1;
 	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
-	for (j=1; j<= (nlstate-1) ; j ++)
-	  fprintf(ficgp,"+$%d",k+l+j);
+	for (j=2; j<= nlstate ; j ++)
+	  fprintf(ficgp,"+$%d",k+l+j-1);
 	fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
       } /* nlstate */
-      fprintf(ficgp,"\n");
+      fprintf(ficgp,"\nset out\n");
     } /* end cpt state*/ 
   } /* end covariate */  
-  
+
   /* proba elementaires */
   fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
   for(i=1,jk=1; i <=nlstate; i++){
@@ -4768,7 +4933,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"##############\n#\n");
 
   /*goto avoid;*/
-  fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n");
+  fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
   fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
   fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -4782,66 +4947,100 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
   fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
   fprintf(ficgp,"#\n");
-   for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+   for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
      fprintf(ficgp,"# ng=%d\n",ng);
      fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);
      for(jk=1; jk <=m; jk++) {
        fprintf(ficgp,"#    jk=%d\n",jk);
-       fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); 
-       if (ng==2)
+       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
+       fprintf(ficgp,"\nset ter svg size 640, 480 ");
+       if (ng==1){
+	 fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
+	 fprintf(ficgp,"\nunset log y");
+       }else if (ng==2){
+	 fprintf(ficgp,"\nset ylabel \"Probability\"\n");
+	 fprintf(ficgp,"\nset log y");
+       }else if (ng==3){
 	 fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
-       else
-	 fprintf(ficgp,"\nset title \"Probability\"\n");
-       fprintf(ficgp,"\nset ter png small size 320, 240\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+	 fprintf(ficgp,"\nset log y");
+       }else
+	 fprintf(ficgp,"\nunset title ");
+       fprintf(ficgp,"\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
        i=1;
        for(k2=1; k2<=nlstate; k2++) {
 	 k3=i;
 	 for(k=1; k<=(nlstate+ndeath); k++) {
 	   if (k != k2){
-	     if(ng==2)
+	     switch( ng) {
+	     case 1:
 	       if(nagesqr==0)
-		 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+		 fprintf(ficgp," p%d+p%d*x",i,i+1);
 	       else /* nagesqr =1 */
-		 fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
-	     else
+		 fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	       break;
+	     case 2: /* ng=2 */
 	       if(nagesqr==0)
 		 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
 	       else /* nagesqr =1 */
-		 fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+		   fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	       break;
+	     case 3:
+	       if(nagesqr==0)
+		 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+	       else /* nagesqr =1 */
+		 fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
+	       break;
+	     }
 	     ij=1;/* To be checked else nbcode[0][0] wrong */
 	     for(j=3; j <=ncovmodel-nagesqr; j++) {
-	       if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
-	       	 fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-	       	 ij++;
+	       /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
+	       if(ij <=cptcovage) { /* Bug valgrind */
+		 if((j-2)==Tage[ij]) { /* Bug valgrind */
+		   fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+		   /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		   ij++;
+		 }
 	       }
 	       else
-		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 	     }
-	     fprintf(ficgp,")/(1");
+	     if(ng != 1){
+	       fprintf(ficgp,")/(1");
 	     
-	     for(k1=1; k1 <=nlstate; k1++){ 
-	       if(nagesqr==0)
-		 fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
-	       else /* nagesqr =1 */
-		 fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-  
-	       ij=1;
-	       for(j=3; j <=ncovmodel-nagesqr; j++){
-		 if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
-		   fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-		   ij++;
+	       for(k1=1; k1 <=nlstate; k1++){ 
+		 if(nagesqr==0)
+		   fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+		 else /* nagesqr =1 */
+		   fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+		 
+		 ij=1;
+		 for(j=3; j <=ncovmodel-nagesqr; j++){
+		   if(ij <=cptcovage) { /* Bug valgrind */
+		     if((j-2)==Tage[ij]) { /* Bug valgrind */
+		       fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+		       /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		       ij++;
+		     }
+		   }
+		   else
+		     fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 		 }
-		 else
-		   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,")");
 	       }
 	       fprintf(ficgp,")");
+	       if(ng ==2)
+		 fprintf(ficgp," t \"p%d%d\" ", k2,k);
+	       else /* ng= 3 */
+		 fprintf(ficgp," t \"i%d%d\" ", k2,k);
+	     }else{ /* end ng <> 1 */
+	       fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
 	     }
-	     fprintf(ficgp,") t \"p%d%d\" ", k2,k);
 	     if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
 	     i=i+ncovmodel;
 	   }
 	 } /* end k */
        } /* end k2 */
+       fprintf(ficgp,"\n set out\n");
      } /* end jk */
    } /* end ng */
  /* avoid: */
@@ -4909,8 +5108,8 @@ void prevforecast(char fileres[], double
   agelim=AGESUP;
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
  
-  strcpy(fileresf,"f"); 
-  strcat(fileresf,fileres);
+  strcpy(fileresf,"F_"); 
+  strcat(fileresf,fileresu);
   if((ficresf=fopen(fileresf,"w"))==NULL) {
     printf("Problem with forecast resultfile: %s\n", fileresf);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
@@ -4959,7 +5158,7 @@ void prevforecast(char fileres[], double
       k=k+1;
       fprintf(ficresf,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficresf,"******\n");
       fprintf(ficresf,"# Covariate valuofcovar yearproj age");
@@ -4983,7 +5182,7 @@ void prevforecast(char fileres[], double
 	    if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
-                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	      fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
 	    } 
 	    for(j=1; j<=nlstate+ndeath;j++) {
@@ -5033,8 +5232,8 @@ void populforecast(char fileres[], doubl
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   
   
-  strcpy(filerespop,"pop"); 
-  strcat(filerespop,fileres);
+  strcpy(filerespop,"POP_"); 
+  strcat(filerespop,fileresu);
   if((ficrespop=fopen(filerespop,"w"))==NULL) {
     printf("Problem with forecast resultfile: %s\n", filerespop);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
@@ -5081,7 +5280,7 @@ void populforecast(char fileres[], doubl
       k=k+1;
       fprintf(ficrespop,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficrespop,"******\n");
       fprintf(ficrespop,"# Age");
@@ -5387,7 +5586,7 @@ double gompertz_f(const gsl_vector *v, v
 #endif
 
 /******************* Printing html file ***********/
-void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtmlmort(char fileresu[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,  double p[],double **matcov,double agemortsup){
   int i,k;
@@ -5396,7 +5595,7 @@ void printinghtmlmort(char fileres[], ch
   fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);
   for (i=1;i<=2;i++) 
     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,"<br><br><img src=\"graphmort.svg\">");
   fprintf(fichtm,"</ul>");
 
 fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>");
@@ -5411,7 +5610,7 @@ fprintf(fichtm,"<ul><li><h4>Life table</
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplotmort(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
 
@@ -5425,9 +5624,9 @@ void printinggnuplotmort(char fileres[],
 
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
-  fprintf(ficgp,"set out \"graphmort.png\"\n "); 
+  fprintf(ficgp,"set out \"graphmort.svg\"\n "); 
   fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); 
-  fprintf(ficgp, "set ter png small size 320, 240\n set log y\n"); 
+  fprintf(ficgp, "set ter svg size 640, 480\n set log y\n"); 
   /* fprintf(ficgp, "set size 0.65,0.65\n"); */
   fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
 
@@ -5738,9 +5937,9 @@ int decodemodel ( char model[], int last
       /* 	k=1 Tvar[1]=2 (from V2) */
       /* 	k=5 Tvar[5] */
       /* for (k=1; k<=cptcovn;k++) { */
-      /* 	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
+      /* 	cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       /* 	} */
-      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
+      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */
       /*
        * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
       for(k=cptcovt; k>=1;k--) /**< Number of covariates */
@@ -6118,25 +6317,25 @@ void syscompilerinfo(int logged)
 
  }
 
-int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl){
   /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
   int i, j, k, i1 ;
-  double ftolpl = 1.e-10;
+  /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;
 
-    strcpy(filerespl,"pl");
-    strcat(filerespl,fileres);
-    if((ficrespl=fopen(filerespl,"w"))==NULL) {
-      printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
-      fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
-    }
-    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");
+  strcpy(filerespl,"PL_");
+  strcat(filerespl,fileresu);
+  if((ficrespl=fopen(filerespl,"w"))==NULL) {
+    printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+    fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+  }
+  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");
   
     /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
 
@@ -6151,14 +6350,14 @@ int prevalence_limit(double *p, double *
       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
 	k=k+1;
 	/* to clean */
-	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
-	fprintf(ficrespl,"\n#******");
-	printf("\n#******");
-	fprintf(ficlog,"\n#******");
+	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+	fprintf(ficrespl,"#******");
+	printf("#******");
+	fprintf(ficlog,"#******");
 	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," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficrespl,"******\n");
 	printf("******\n");
@@ -6166,7 +6365,7 @@ int prevalence_limit(double *p, double *
 
 	fprintf(ficrespl,"#Age ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
 	fprintf(ficrespl,"\n");
@@ -6176,7 +6375,7 @@ int prevalence_limit(double *p, double *
 	  prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
 	  fprintf(ficrespl,"%.0f ",age );
 	  for(j=1;j<=cptcoveff;j++)
-	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	  for(i=1; i<=nlstate;i++)
 	    fprintf(ficrespl," %.5f", prlim[i][i]);
 	  fprintf(ficrespl,"\n");
@@ -6198,7 +6397,7 @@ int hPijx(double *p, int bage, int fage)
   double agedeb;
   double ***p3mat;
 
-    strcpy(filerespij,"pij");  strcat(filerespij,fileres);
+    strcpy(filerespij,"PIJ_");  strcat(filerespij,fileresu);
     if((ficrespij=fopen(filerespij,"w"))==NULL) {
       printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
       fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
@@ -6223,7 +6422,7 @@ int hPijx(double *p, int bage, int fage)
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
-	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       fprintf(ficrespij,"******\n");
       
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -6275,6 +6474,7 @@ int main(int argc, char *argv[])
 
   int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */
+  int num_filled;
   int itimes;
   int NDIM=2;
   int vpopbased=0;
@@ -6293,7 +6493,9 @@ int main(int argc, char *argv[])
   double ***mobaverage;
 
   char line[MAXLINE];
-  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
+  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
+
+  char model[MAXLINE], modeltemp[MAXLINE];
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
@@ -6373,7 +6575,7 @@ int main(int argc, char *argv[])
   getcwd(pathcd, size);
 #endif
   syscompilerinfo(0);
-  printf("\n%s\n%s",version,fullversion);
+  printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
     fgets(pathr,FILENAMELENGTH,stdin);
@@ -6437,7 +6639,7 @@ int main(int argc, char *argv[])
     goto end;
   }
   fprintf(ficlog,"Log filename:%s\n",filelog);
-  fprintf(ficlog,"\n%s\n%s",version,fullversion);
+  fprintf(ficlog,"Version %s %s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\
@@ -6445,7 +6647,7 @@ int main(int argc, char *argv[])
  optionfilext=%s\n\
  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
-  syscompilerinfo(0);
+  syscompilerinfo(1);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -6456,7 +6658,9 @@ int main(int argc, char *argv[])
   /* */
   strcpy(fileres,"r");
   strcat(fileres, optionfilefiname);
+  strcat(fileresu, optionfilefiname); /* Without r in front */
   strcat(fileres,".txt");    /* Other files have txt extension */
+  strcat(fileresu,".txt");    /* Other files have txt extension */
 
   /* Main ---------arguments file --------*/
 
@@ -6471,7 +6675,7 @@ int main(int argc, char *argv[])
 
 
   strcpy(filereso,"o");
-  strcat(filereso,fileres);
+  strcat(filereso,fileresu);
   if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
     printf("Problem with Output resultfile: %s\n", filereso);
     fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
@@ -6481,23 +6685,84 @@ int main(int argc, char *argv[])
 
   /* Reads comments: lines beginning with '#' */
   numlinepar=0;
-  while((c=getc(ficpar))=='#' && c!= EOF){
-    ungetc(c,ficpar);
-    fgets(line, MAXLINE, ficpar);
+
+    /* First parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \
+			title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
+    if (num_filled != 5) {
+      printf("Should be 5 parameters\n");
+    }
     numlinepar++;
-    fputs(line,stdout);
-    fputs(line,ficparo);
-    fputs(line,ficlog);
+    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
   }
-  ungetc(c,ficpar);
-
-  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=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
-  numlinepar=numlinepar+3; /* In general */
-  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=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+  /* Second parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
+			&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
+    if (num_filled != 8) {
+      printf("Not 8\n");
+    }
+    printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
+  }
+  ftolpl=6*ftol*1.e5; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
+  /* Third parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
+    if (num_filled == 0)
+            model[0]='\0';
+    else if (num_filled != 1){
+      printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+      model[0]='\0';
+      goto end;
+    }
+    else{
+      if (model[0]=='+'){
+	for(i=1; i<=strlen(model);i++)
+	  modeltemp[i-1]=model[i];
+	strcpy(model,modeltemp); 
+      }
+    }
+    /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
+  }
+  /* 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=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
+  /* numlinepar=numlinepar+3; /\* In general *\/ */
+  /* 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=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
   if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
     model[strlen(model)-1]='\0';
-  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=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-  fprintf(ficlog,"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=1+age+%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=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficlog,"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=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
   fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){
@@ -6721,8 +6986,8 @@ Please run with mle=-1 to get a correct
     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;
+      printf("Problem writing new parameter file: %s\n", rfileres);goto end;
+      fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
     }
     fprintf(ficres,"#%s\n",version);
   }    /* End of mle != -3 */
@@ -6832,8 +7097,8 @@ Please run with mle=-1 to get a correct
      V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
   /* 1 to ncodemax[j] is the maximum value of this jth covariate */
 
-  codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
-  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/
+  /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
+  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
   h=0;
 
@@ -6843,13 +7108,6 @@ Please run with mle=-1 to get a correct
  
   m=pow(2,cptcoveff);
  
-  for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
-    for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ 
-      for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/
-	for(cpt=1; cpt <=pow(2,k-1); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ 
-	  h++;
-	  if (h>m) 
-	    h=1;
 	  /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
 	   * For k=4 covariates, h goes from 1 to 2**k
 	   * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
@@ -6863,32 +7121,47 @@ Please run with mle=-1 to get a correct
 	   *     6     2     1     2     1
 	   *     7 i=4 1     2     2     1
 	   *     8     2     2     2     1
-	   *     9 i=5 1 i=3 1 i=2 1     1
-	   *    10     2     1     1     1
-	   *    11 i=6 1     2     1     1
-	   *    12     2     2     1     1
-	   *    13 i=7 1 i=4 1     2     1    
-	   *    14     2     1     2     1
-	   *    15 i=8 1     2     2     1
-	   *    16     2     2     2     1
+	   *     9 i=5 1 i=3 1 i=2 1     2
+	   *    10     2     1     1     2
+	   *    11 i=6 1     2     1     2
+	   *    12     2     2     1     2
+	   *    13 i=7 1 i=4 1     2     2    
+	   *    14     2     1     2     2
+	   *    15 i=8 1     2     2     2
+	   *    16     2     2     2     2
 	   */
-	  codtab[h][k]=j;
-	  /* codtab[12][3]=1; */
-	  /*codtab[h][Tvar[k]]=j;*/
-	  printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
-	} 
-      }
-    }
-  } 
+  /* /\* for(h=1; h <=100 ;h++){  *\/ */
+  /*   /\* printf("h=%2d ", h); *\/ */
+  /*    /\* for(k=1; k <=10; k++){ *\/ */
+  /*      /\* printf("k=%d %d ",k,codtabm(h,k)); *\/ */
+  /*    /\*   codtab[h][k]=codtabm(h,k); *\/ */
+  /*    /\* } *\/ */
+  /*    /\* printf("\n"); *\/ */
+  /* } */
+  /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
+  /*   for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/  */
+  /*     for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */
+  /* 	for(cpt=1; cpt <=pow(2,k-1); cpt++){  /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/  */
+  /* 	  h++; */
+  /* 	  if (h>m)  */
+  /* 	    h=1; */
+  /* 	  codtab[h][k]=j; */
+  /* 	  /\* codtab[12][3]=1; *\/ */
+  /* 	  /\*codtab[h][Tvar[k]]=j;*\/ */
+  /* 	  /\* printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][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");
-     }
-     scanf("%d",i);*/
+  /* 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);*/
 
  free_ivector(Ndum,-1,NCOVMAX);
 
@@ -6897,7 +7170,7 @@ Please run with mle=-1 to get a correct
   /* Initialisation of ----------- gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);
   if(mle==-3)
-    strcat(optionfilegnuplot,"-mort");
+    strcat(optionfilegnuplot,"-MORT_");
   strcat(optionfilegnuplot,".gp");
 
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
@@ -6916,7 +7189,7 @@ Please run with mle=-1 to get a correct
 
   strcpy(optionfilehtm,optionfilefiname); /* Main html file */
   if(mle==-3)
-    strcat(optionfilehtm,"-mort");
+    strcat(optionfilehtm,"-MORT_");
   strcat(optionfilehtm,".htm");
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
     printf("Problem with %s \n",optionfilehtm);
@@ -7032,8 +7305,8 @@ Interval (in months) between two waves:
 #else
     printf("Powell\n");  fprintf(ficlog,"Powell\n");
 #endif
-    strcpy(filerespow,"pow-mort"); 
-    strcat(filerespow,fileres);
+    strcpy(filerespow,"POW-MORT_"); 
+    strcat(filerespow,fileresu);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
@@ -7187,8 +7460,8 @@ Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else
-      printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
-    printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
+      printinggnuplotmort(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+    printinghtmlmort(fileresu,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
 		     model,imx,p,matcov,agemortsup);
     
@@ -7251,16 +7524,16 @@ Please run with mle=-1 to get a correct
       ftolhess=ftol; /* Usually correct */
       hesscov(matcov, p, npar, delti, ftolhess, func);
     }
-    printf("Parameters and 95%% confidence intervals\n");
-    fprintf(ficlog, "Parameters, T and confidence intervals\n");
+    printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+    fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\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);
 	  for(j=1; j <=ncovmodel; j++){
-	    printf("%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk]));
-	    fprintf(ficlog,"%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk]));
+	    printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+	    fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
 	    jk++; 
 	  }
 	  printf("\n");
@@ -7458,9 +7731,9 @@ Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else
-      printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
-    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
+    printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
 		 model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
 		 jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
       
@@ -7485,7 +7758,7 @@ Please run with mle=-1 to get a correct
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
     prlim=matrix(1,nlstate,1,nlstate);
-    prevalence_limit(p, prlim,  ageminpar, agemaxpar);
+    prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl);
     fclose(ficrespl);
 
 #ifdef FREEEXIT2
@@ -7512,7 +7785,7 @@ Please run with mle=-1 to get a correct
     /*if((stepm == 1) && (strcmp(model,".")==0)){*/
     if(prevfcast==1){
       /*    if(stepm ==1){*/
-      prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+      prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
       /*      }  */
       /*      else{ */
@@ -7542,8 +7815,8 @@ Please run with mle=-1 to get a correct
 
     /*---------- Health expectancies, no variances ------------*/
 
-    strcpy(filerese,"e");
-    strcat(filerese,fileres);
+    strcpy(filerese,"E_");
+    strcat(filerese,fileresu);
     if((ficreseij=fopen(filerese,"w"))==NULL) {
       printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
       fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
@@ -7556,7 +7829,7 @@ Please run with mle=-1 to get a correct
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
 	fprintf(ficreseij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficreseij,"******\n");
 
@@ -7573,8 +7846,8 @@ Please run with mle=-1 to get a correct
     /*---------- Health expectancies and variances ------------*/
 
 
-    strcpy(filerest,"t");
-    strcat(filerest,fileres);
+    strcpy(filerest,"T_");
+    strcat(filerest,fileresu);
     if((ficrest=fopen(filerest,"w"))==NULL) {
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;
       fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
@@ -7583,8 +7856,8 @@ Please run with mle=-1 to get a correct
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
 
 
-    strcpy(fileresstde,"stde");
-    strcat(fileresstde,fileres);
+    strcpy(fileresstde,"STDE_");
+    strcat(fileresstde,fileresu);
     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);
@@ -7592,8 +7865,8 @@ Please run with mle=-1 to get a correct
     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);
+    strcpy(filerescve,"CVE_");
+    strcat(filerescve,fileresu);
     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);
@@ -7601,8 +7874,8 @@ Please run with mle=-1 to get a correct
     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);
+    strcpy(fileresv,"V_");
+    strcat(fileresv,fileresu);
     if((ficresvij=fopen(fileresv,"w"))==NULL) {
       printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
       fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
@@ -7616,21 +7889,21 @@ Please run with mle=-1 to get a correct
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
     	fprintf(ficrest,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficrest,"******\n");
 
 	fprintf(ficresstdeij,"\n#****** ");
 	fprintf(ficrescveij,"\n#****** ");
 	for(j=1;j<=cptcoveff;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(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficresstdeij,"******\n");
 	fprintf(ficrescveij,"******\n");
 
 	fprintf(ficresvij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficresvij,"******\n");
 
 	eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
@@ -7645,7 +7918,7 @@ Please run with mle=-1 to get a correct
 
 
 	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-	  oldm=oldms;savm=savms; /* Segmentation fault */
+	  oldm=oldms;savm=savms; /* ZZ Segmentation fault */
 	  cptcod= 0; /* To be deleted */
 	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
 	  fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
@@ -7653,13 +7926,13 @@ Please run with mle=-1 to get a correct
 	    fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
 	  else
 	    fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
-	  fprintf(ficrest,"# Age e.. (std) ");
+	  fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
 	  for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	  fprintf(ficrest,"\n");
-
+	  /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
 	  epj=vector(1,nlstate+1);
 	  for(age=bage; age <=fage ;age++){
-	    prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+	    prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); /*ZZ Is it the correct prevalim */
 	    if (vpopbased==1) {
 	      if(mobilav ==0){
 		for(i=1; i<=nlstate;i++)
@@ -7670,14 +7943,17 @@ Please run with mle=-1 to get a correct
 	      }
 	    }
 	
-	    fprintf(ficrest," %4.0f",age);
+	    fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+	    /* printf(" age %4.0f ",age); */
 	    for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
 	      for(i=1, epj[j]=0.;i <=nlstate;i++) {
 		epj[j] += prlim[i][i]*eij[i][j][(int)age];
-		/*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+		/*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+	        /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
 	      }
 	      epj[nlstate+1] +=epj[j];
 	    }
+	    /* printf(" age %4.0f \n",age); */
 
 	    for(i=1, vepp=0.;i <=nlstate;i++)
 	      for(j=1;j <=nlstate;j++)
@@ -7709,8 +7985,8 @@ Please run with mle=-1 to get a correct
   
     /*------- Variance of period (stable) prevalence------*/   
 
-    strcpy(fileresvpl,"vpl");
-    strcat(fileresvpl,fileres);
+    strcpy(fileresvpl,"VPL_");
+    strcat(fileresvpl,fileresu);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
       printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
       exit(0);
@@ -7723,7 +7999,7 @@ Please run with mle=-1 to get a correct
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
     	fprintf(ficresvpl,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficresvpl,"******\n");
       
 	varpl=matrix(1,nlstate,(int) bage, (int) fage);
@@ -7760,7 +8036,7 @@ Please run with mle=-1 to get a correct
     free_ivector(Tage,1,NCOVMAX);
 
     free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
-    free_imatrix(codtab,1,100,1,10);
+    /* free_imatrix(codtab,1,100,1,10); */
   fflush(fichtm);
   fflush(ficgp);