--- imach/src/imach.c	2015/09/01 18:24:39	1.197
+++ imach/src/imach.c	2015/09/22 19:45:16	1.202
@@ -1,6 +1,26 @@
-/* $Id: imach.c,v 1.197 2015/09/01 18:24: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 ***
 
@@ -627,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 *\/ */
@@ -714,12 +735,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.197 2015/09/01 18:24:39 brouard Exp $ */
+/* $Id: imach.c,v 1.202 2015/09/22 19:45:16 brouard Exp $ */
 /* $State: Exp $ */
 #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.197 $ $Date: 2015/09/01 18:24:39 $"; 
+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  */
@@ -785,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];
@@ -875,7 +896,7 @@ 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;
 
@@ -1589,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; */
@@ -1613,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]) */
@@ -1625,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]); */
@@ -1636,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); 
@@ -1674,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; 
@@ -1819,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);
@@ -1834,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 */
@@ -1844,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++) { 
@@ -1894,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++){
@@ -1904,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]);*/
@@ -1936,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 */
 }
 
@@ -2090,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);*/
@@ -2519,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;
@@ -2553,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);
@@ -2570,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;
 }
@@ -2603,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);
@@ -2930,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);
@@ -2966,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*/
             } 
 	}
  
@@ -3001,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++) 
@@ -3132,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) { 
@@ -3799,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);
@@ -3813,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);
@@ -3836,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);
 /*   } */
@@ -4048,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); */
@@ -4058,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);
@@ -4180,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);
@@ -4235,7 +4286,7 @@ 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 one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
+  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.\
@@ -4258,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#");    
       }
       
@@ -4287,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
@@ -4295,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++){
@@ -4445,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",\
@@ -4472,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 */
@@ -4494,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 ,\
@@ -4507,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>");
 
@@ -4534,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 */
@@ -4572,31 +4641,31 @@ covariance matrix of the one-step probab
 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 */
@@ -4617,20 +4686,20 @@ See page 'Matrix of variance-covariance
      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>");
@@ -4638,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); */
@@ -4653,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);
@@ -4732,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++){
@@ -4783,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");
@@ -4797,71 +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++) {
 	       /* 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]][codtab[jk][Tvar[j-2]]]);
+		   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(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]][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: */
@@ -4929,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);
@@ -4979,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");
@@ -5003,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++) {
@@ -5053,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);
@@ -5101,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");
@@ -5407,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;
@@ -5416,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>");
@@ -5431,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];
 
@@ -5445,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);
 
@@ -5758,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 */
@@ -6138,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 */
 
@@ -6171,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");
@@ -6186,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");
@@ -6196,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");
@@ -6218,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;
@@ -6243,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 */
@@ -6479,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 --------*/
 
@@ -6494,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);
@@ -6544,7 +6725,7 @@ int main(int argc, char *argv[])
     }
     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 */
@@ -6557,8 +6738,10 @@ int main(int argc, char *argv[])
     }else
       break;
   }
-  if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){
-    if (num_filled != 1) {
+  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';
@@ -6568,10 +6751,10 @@ int main(int argc, char *argv[])
       if (model[0]=='+'){
 	for(i=1; i<=strlen(model);i++)
 	  modeltemp[i-1]=model[i];
+	strcpy(model,modeltemp); 
       }
-      strcpy(model,modeltemp); 
     }
-    printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout);
+    /* 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 *\/ */
@@ -6803,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 */
@@ -6914,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;
 
@@ -6947,14 +7130,14 @@ Please run with mle=-1 to get a correct
 	   *    15 i=8 1     2     2     2
 	   *    16     2     2     2     2
 	   */
-  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(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*\/ */
@@ -6987,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) {
@@ -7006,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);
@@ -7122,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);
@@ -7277,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);
     
@@ -7548,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);
       
@@ -7575,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
@@ -7602,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{ */
@@ -7632,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);
@@ -7646,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");
 
@@ -7663,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;
@@ -7673,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);
@@ -7682,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);
@@ -7691,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);
@@ -7706,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);
@@ -7735,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 ");
@@ -7743,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++)
@@ -7760,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++)
@@ -7799,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);
@@ -7813,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);
@@ -7850,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);