--- imach/src/imach.c	2016/09/15 15:01:13	1.251
+++ imach/src/imach.c	2017/04/26 06:01:29	1.264
@@ -1,6 +1,47 @@
-/* $Id: imach.c,v 1.251 2016/09/15 15:01:13 brouard Exp $
+/* $Id: imach.c,v 1.264 2017/04/26 06:01:29 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.264  2017/04/26 06:01:29  brouard
+  Summary: Labels in graphs
+
+  Revision 1.263  2017/04/24 15:23:15  brouard
+  Summary: to save
+
+  Revision 1.262  2017/04/18 16:48:12  brouard
+  *** empty log message ***
+
+  Revision 1.261  2017/04/05 10:14:09  brouard
+  Summary: Bug in E_ as well as in T_ fixed nres-1 vs k1-1
+
+  Revision 1.260  2017/04/04 17:46:59  brouard
+  Summary: Gnuplot indexations fixed (humm)
+
+  Revision 1.259  2017/04/04 13:01:16  brouard
+  Summary: Some errors to warnings only if date of death is unknown but status is death we could set to pi3
+
+  Revision 1.258  2017/04/03 10:17:47  brouard
+  Summary: Version 0.99r12
+
+  Some cleanings, conformed with updated documentation.
+
+  Revision 1.257  2017/03/29 16:53:30  brouard
+  Summary: Temp
+
+  Revision 1.256  2017/03/27 05:50:23  brouard
+  Summary: Temporary
+
+  Revision 1.255  2017/03/08 16:02:28  brouard
+  Summary: IMaCh version 0.99r10 bugs in gnuplot fixed
+
+  Revision 1.254  2017/03/08 07:13:00  brouard
+  Summary: Fixing data parameter line
+
+  Revision 1.253  2016/12/15 11:59:41  brouard
+  Summary: 0.99 in progress
+
+  Revision 1.252  2016/09/15 21:15:37  brouard
+  *** empty log message ***
+
   Revision 1.251  2016/09/15 15:01:13  brouard
   Summary: not working
 
@@ -128,9 +169,7 @@
   Author: Nicolas Brouard
 
   Revision 1.210  2015/11/18 17:41:20  brouard
-  Summary: Start working on projected prevalences
-
-  Revision 1.209  2015/11/17 22:12:03  brouard
+  Summary: Start working on projected prevalences  Revision 1.209  2015/11/17 22:12:03  brouard
   Summary: Adding ftolpl parameter
   Author: N Brouard
 
@@ -955,12 +994,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.251 2016/09/15 15:01:13 brouard Exp $ */
+/* $Id: imach.c,v 1.264 2017/04/26 06:01:29 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.251 $ $Date: 2016/09/15 15:01:13 $"; 
+char fullversion[]="$Revision: 1.264 $ $Date: 2017/04/26 06:01:29 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1161,6 +1200,7 @@ int *TvarsQind;
 
 #define MAXRESULTLINES 10
 int nresult=0;
+int parameterline=0; /* # of the parameter (type) line */
 int TKresult[MAXRESULTLINES];
 int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
 int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
@@ -2393,8 +2433,8 @@ void powell(double p[], double **xi, int
 	    flatd++;
 	}
 	if(flatd >0){
-	  printf("%d flat directions\n",flatd);
-	  fprintf(ficlog,"%d flat directions\n",flatd);
+	  printf("%d flat directions: ",flatd);
+	  fprintf(ficlog,"%d flat directions :",flatd);
 	  for (j=1;j<=n;j++) { 
 	    if(flatdir[j]>0){
 	      printf("%d ",j);
@@ -2578,7 +2618,7 @@ Earliest age to start was %d-%d=%d, ncvl
  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
   double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres)
 {
-  /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
+  /* Computes the prevalence limit in each live state at age x and for covariate combination ij (<=2**cptcoveff) by left multiplying the unit
      matrix by transitions matrix until convergence is reached with precision ftolpl */
   /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
   /* Wx is row vector: population in state 1, population in state 2, population dead */
@@ -2639,7 +2679,7 @@ Earliest age to start was %d-%d=%d, ncvl
     for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
  			/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
       cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
-      /* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
+      /* printf("bprevalim Dummy agefin=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agefin,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
     }
     /* for (k=1; k<=cptcovn;k++) { */
     /*   /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
@@ -4122,7 +4162,16 @@ void ludcmp(double **a, int n, int *indx
     big=0.0; 
     for (j=1;j<=n;j++) 
       if ((temp=fabs(a[i][j])) > big) big=temp; 
-    if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); 
+    if (big == 0.0){
+      printf(" Singular Hessian matrix at row %d:\n",i);
+      for (j=1;j<=n;j++) {
+	printf(" a[%d][%d]=%f,",i,j,a[i][j]);
+	fprintf(ficlog," a[%d][%d]=%f,",i,j,a[i][j]);
+      }
+      fflush(ficlog);
+      fclose(ficlog);
+      nrerror("Singular matrix in routine ludcmp"); 
+    }
     vv[i]=1.0/big; 
   } 
   for (j=1;j<=n;j++) { 
@@ -4188,17 +4237,84 @@ void pstamp(FILE *fichier)
   fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
+int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
+
+  /* y=a+bx regression */
+  double   sumx = 0.0;                        /* sum of x                      */
+  double   sumx2 = 0.0;                       /* sum of x**2                   */
+  double   sumxy = 0.0;                       /* sum of x * y                  */
+  double   sumy = 0.0;                        /* sum of y                      */
+  double   sumy2 = 0.0;                       /* sum of y**2                   */
+  double   sume2; /* sum of square or residuals */
+  double yhat;
+  
+  double denom=0;
+  int i;
+  int ne=*no;
+  
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    sumx  += x[i];       
+    sumx2 += x[i]*x[i];  
+    sumxy += x[i] * y[i];
+    sumy  += y[i];      
+    sumy2 += y[i]*y[i]; 
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  
+  denom = (ne * sumx2 - sumx*sumx);
+  if (denom == 0) {
+    // vertical, slope m is infinity
+    *b = INFINITY;
+    *a = 0;
+    if (r) *r = 0;
+    return 1;
+  }
+  
+  *b = (ne * sumxy  -  sumx * sumy) / denom;
+  *a = (sumy * sumx2  -  sumx * sumxy) / denom;
+  if (r!=NULL) {
+    *r = (sumxy - sumx * sumy / ne) /          /* compute correlation coeff     */
+      sqrt((sumx2 - sumx*sumx/ne) *
+	   (sumy2 - sumy*sumy/ne));
+  }
+  *no=ne;
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    yhat = y[i] - *a -*b* x[i];
+    sume2  += yhat * yhat ;       
+    
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne));
+  *sa= *sb * sqrt(sumx2/ne);
+  
+  return 0; 
+}
+
 /************ Frequencies ********************/
 void  freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
 		  int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],	\
 		  int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies as well as proposing some starting values */
   
-  int i, m, jk, j1, bool, z1,j, k, iv, jj=0;
+  int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0;
   int iind=0, iage=0;
   int mi; /* Effective wave */
   int first;
   double ***freq; /* Frequencies */
+  double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */
+  int no;
   double *meanq;
   double **meanqt;
   double *pp, **prop, *posprop, *pospropt;
@@ -4251,6 +4367,8 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   }
   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
   
+  y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
+  x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   j1=0;
   
@@ -4268,7 +4386,15 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   dateintsum=0;
   k2cpt=0;
 
-  for (j = 0; j <= cptcoveff; j+=cptcoveff){   /* j= 0 constant model */
+  if(cptcoveff == 0 )
+    nl=1;  /* Constant model only */
+  else
+    nl=2;
+  for (nj = 1; nj <= nl; nj++){   /* nj= 1 constant model, nl number of loops. */
+    if(nj==1)
+      j=0;  /* First pass for the constant */
+    else
+      j=cptcoveff; /* Other passes for the covariate values */
     first=1;
     for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
       posproptt=0.;
@@ -4496,6 +4622,8 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	  if( iage <= iagemax){
 	    if(pos>=1.e-5){
 	      fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+	    /* fprintf(ficresp, "%d %d %d %.5f %.0f %.0f",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)],iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); */
+
 	      fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
 	      /*probs[iage][jk][j1]= pp[jk]/pos;*/
 	      /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
@@ -4514,6 +4642,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	      if(first==1){
 		printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
 	      }
+	      /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); */
 	      fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
 	    }
 	    if(jk!=0 && m!=0)
@@ -4547,7 +4676,8 @@ Title=%s <br>Datafile=%s Firstpass=%d La
       if(posproptt < 1.e-5){
 	fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
 	fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
-	fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);
+	fprintf(ficlog,"#  This combination (%d) is not valid and no result will be produced\n",j1);
+	printf("#  This combination (%d) is not valid and no result will be produced\n",j1);
 	invalidvarcomb[j1]=1;
       }else{
 	fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
@@ -4560,22 +4690,33 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	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(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */
-		if(jj==1){  /* Constant case */
+		if(jj==1){  /* Constant case (in fact cste + age) */
 		  if(j1==1){ /* All dummy covariates to zero */
 		    freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */
 		    freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */
+		    printf("%d%d ",i,k);
+		    fprintf(ficlog,"%d%d ",i,k);
+		    printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
+		    fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+		    pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
+		  }
+		}else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */
+		  for(iage=iagemin; iage <= iagemax+3; iage++){
+		    x[iage]= (double)iage;
+		    y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
+		    /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */
 		  }
-		  printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
-		  fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
-		  pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
-		}else if( (log(j1-1)/log(2)+1 == jj -2 -nagesqr)  && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ 
+		  linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
+		  pstart[jk]=b;
+		  pstart[jk-1]=a;
+		}else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj)  && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ 
+		  printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
+		  printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
 		  pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]));
+		  printf("%d%d ",i,k);
+		  fprintf(ficlog,"%d%d ",i,k);
 		  printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4]));
-		}else if(jj==2 || nagesqr==1){ /* age or age*age parameter */
-		  ;
 		}else{ /* Other cases, like quantitative fixed or varying covariates */
 		  ;
 		}
@@ -4583,8 +4724,6 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 		/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
 		jk++; 
 	      } /* end jj */
-	      printf("\n");
-	      fprintf(ficlog,"\n");
 	    } /* end k!= i */
 	  } /* end k */
 	} /* end i, jk */
@@ -4599,7 +4738,9 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	    printf("%d%d ",i,k);
 	    fprintf(ficlog,"%d%d ",i,k);
 	    for(jj=1; jj <=ncovmodel; jj++){
-	      if(jj==1){
+	      pstart[jk]=p[jk]; /* Setting pstart to p values by default */
+	      if(jj==1){ /* Age has to be done */
+		pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
 		printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
 		fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
 	      }
@@ -4630,6 +4771,28 @@ Title=%s <br>Datafile=%s Firstpass=%d La
       fprintf(ficlog,"\n");
     } /* end j=0 */
   } /* end j */
+
+  if(mle == -2){  /* We want to use these values as starting values */
+    for(i=1, jk=1; i <=nlstate; i++){
+      for(j=1; j <=nlstate+ndeath; j++){
+	if(j!=i){
+	  /*ca[0]= k+'a'-1;ca[1]='\0';*/
+	  printf("%1d%1d",i,j);
+	  fprintf(ficparo,"%1d%1d",i,j);
+	  for(k=1; k<=ncovmodel;k++){
+	    /* 	  printf(" %lf",param[i][j][k]); */
+	    /* 	  fprintf(ficparo," %lf",param[i][j][k]); */
+	    p[jk]=pstart[jk];
+	    printf(" %f ",pstart[jk]);
+	    fprintf(ficparo," %f ",pstart[jk]);
+	    jk++;
+	  }
+	  printf("\n");
+	  fprintf(ficparo,"\n");
+	}
+      }
+    }
+  } /* end mle=-2 */
   dateintmean=dateintsum/k2cpt; 
   
   fclose(ficresp);
@@ -4637,6 +4800,8 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   fclose(ficresphtmfr);
   free_vector(meanq,1,nqfveff);
   free_matrix(meanqt,1,lastpass,1,nqtveff);
+  free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
+  free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pospropt,1,nlstate);
   free_vector(posprop,1,nlstate);
@@ -4794,10 +4959,10 @@ void  concatwav(int wav[], int **dh, int
 #else
 	if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
 	  if(firsthree == 0){
-	    printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+	    printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
 	    firsthree=1;
 	  }
-	  fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+	  fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p%d%d .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath);
 	  mw[++mi][i]=m;
 	  mli=m;
 	}
@@ -4824,13 +4989,9 @@ void  concatwav(int wav[], int **dh, int
       /* if(mi==0)  never been interviewed correctly before death */
       /* Only death is a correct wave */
       mw[mi][i]=m;
-    }
+    } /* else not in a death state */
 #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE
-    else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */
-      /* m++; */
-      /* mi++; */
-      /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */
-      /* mw[mi][i]=m; */
+    else if ((int) andc[i] != 9999) {  /* Date of death is known */
       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
 	if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
 	  nbwarn++;
@@ -4843,12 +5004,12 @@ void  concatwav(int wav[], int **dh, int
 	}else{ /* Death occured afer last wave potential bias */
 	  nberr++;
 	  if(firstwo==0){
-	    printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+	    printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
 	    firstwo=1;
 	  }
-	  fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+	  fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
 	}
-      }else{ /* end date of interview is known */
+      }else{ /* if date of interview is unknown */
 	/* death is known but not confirmed by death status at any wave */
 	if(firstfour==0){
 	  printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
@@ -4879,7 +5040,7 @@ void  concatwav(int wav[], int **dh, int
       if (stepm <=0)
 	dh[mi][i]=1;
       else{
-	if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
+	if (s[mw[mi+1][i]][i] > nlstate) { /* A death, but what if date is unknown? */
 	  if (agedc[i] < 2*AGESUP) {
 	    j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
 	    if(j==0) j=1;  /* Survives at least one month after exam */
@@ -6241,7 +6402,7 @@ To be simple, these graphs help to under
 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 prevfcast, int backcast, int estepm , \
+		  int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
 		  double jprev1, double mprev1,double anprev1, double dateprev1, \
 		  double jprev2, double mprev2,double anprev2, double dateprev2){
   int jj1, k1, i1, cpt, k4, nres;
@@ -6279,21 +6440,66 @@ void printinghtml(char fileresu[], char
    <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
    }
 
-   fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
    m=pow(2,cptcoveff);
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}
 
+   fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
+
+   jj1=0;
+
+   fprintf(fichtm," \n<ul>");
+   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
+     if(m != 1 && TKresult[nres]!= k1)
+       continue;
+     jj1++;
+     if (cptcovn > 0) {
+       fprintf(fichtm,"\n<li><a  size=\"1\" color=\"#EC5E5E\" href=\"#rescov");
+       for (cpt=1; cpt<=cptcoveff;cpt++){ 
+	 fprintf(fichtm,"_V%d=%d_",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+       }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	 fprintf(fichtm,"_V%d=%f_",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
+       fprintf(fichtm,"\">");
+       
+       /* if(nqfveff+nqtveff 0) */ /* Test to be done */
+       fprintf(fichtm,"************ Results for covariates");
+       for (cpt=1; cpt<=cptcoveff;cpt++){ 
+	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+       }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	 fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
+       if(invalidvarcomb[k1]){
+	 fprintf(fichtm," Warning Combination (%d) ignored because no cases ",k1); 
+	 continue;
+       }
+       fprintf(fichtm,"</a></li>");
+     } /* cptcovn >0 */
+   }
+     fprintf(fichtm," \n</ul>");
+
    jj1=0;
 
    for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
 
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
      if (cptcovn > 0) {
+       fprintf(fichtm,"\n<p><a name=\"rescov");
+       for (cpt=1; cpt<=cptcoveff;cpt++){ 
+	 fprintf(fichtm,"_V%d=%d_",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
+       }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+	 fprintf(fichtm,"_V%d=%f_",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
+       fprintf(fichtm,"\"</a>");
+ 
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
 	 fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
@@ -6315,7 +6521,7 @@ void printinghtml(char fileresu[], char
        }
      }
      /* aij, bij */
-     fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
+     fprintf(fichtm,"<br>- Logit model (yours is: logit(pij)=log(pij/pii)= aij+ bij age+%s) as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
 <img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
      /* Pij */
      fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \
@@ -6339,21 +6545,21 @@ divided by h: <sub>h</sub>P<sub>ij</sub>
      }
      /* 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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }
      if(backcast==1){
        /* Period (stable) back prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Convergence to period (stable) back 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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+	 fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
        }
      }
      if(prevfcast==1){
        /* Projection of prevalence up to period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up 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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up 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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }
      }
 	 
@@ -6420,7 +6626,7 @@ See page 'Matrix of variance-covariance
 
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
    for(k1=1; k1<=m;k1++){
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
@@ -6441,9 +6647,9 @@ See page 'Matrix of variance-covariance
        }
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
-       fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
+       fprintf(fichtm,"\n<br>- Observed (cross-sectional with mov_average=%d) and period (incidence based) \
 prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\
-<img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
+<img src=\"%s_%d-%d-%d.svg\">",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
      }
      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) \
@@ -6462,7 +6668,7 @@ true period expectancies (those weighted
 void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
-  char gplotcondition[132];
+  char gplotcondition[132], gplotlabel[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;
   int ng=0;
@@ -6516,11 +6722,12 @@ void printinggnuplot(char fileresu[], ch
     for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
 	/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-	if(TKresult[nres]!= k1)
+	if(m != 1 && TKresult[nres]!= k1)
 	  continue;
 	/* We are interested in selected combination by the resultline */
 	/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
 	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
+	strcpy(gplotlabel,"(");
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -6530,32 +6737,38 @@ void printinggnuplot(char fileresu[], ch
 	  /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
 	  /* printf(" V%d=%d ",Tvaraff[k],vlv); */
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-	}	
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	}
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	/* printf("\n#\n"); */
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
+          /*k1=k1-1;*/ /* To be checked */
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
 	  continue;
 	}
       
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
 	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
-	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
-      
+	fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
+	/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
+      /* k1-1 error should be nres-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==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
+	fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
 	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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); 
+	fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); 
 	for (i=1; i<= nlstate ; i ++) {
 	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	  else fprintf(ficgp," %%*lf (%%*lf)");
@@ -6589,7 +6802,7 @@ void printinggnuplot(char fileresu[], ch
 	    } /* end covariate */
 	  } /* end if no covariate */
 	} /* end if backcast */
-	fprintf(ficgp,"\nset out \n");
+	fprintf(ficgp,"\nset out ;unset label;\n");
       } /* nres */
     } /* k1 */
   } /* cpt */
@@ -6598,9 +6811,10 @@ void printinggnuplot(char fileresu[], ch
   /*2 eme*/
   for (k1=1; k1<= m ; k1 ++){  
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+      strcpy(gplotlabel,"(");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -6608,12 +6822,15 @@ void printinggnuplot(char fileresu[], ch
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
       }
       /* for(k=1; k <= ncovds; k++){ */
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }
+      strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){
 	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6622,26 +6839,27 @@ void printinggnuplot(char fileresu[], ch
 			
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-	if(vpopbased==0)
+	fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);
+	if(vpopbased==0){
 	  fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-	else
+	}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);
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-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);
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-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);
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
 	  for (j=1; j<= nlstate+1 ; j ++) {
 	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
 	    else fprintf(ficgp," %%*lf (%%*lf)");
@@ -6650,7 +6868,7 @@ void printinggnuplot(char fileresu[], ch
 	  else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
 	} /* state */
       } /* vpopbased */
-      fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
+      fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
     } /* end nres */
   } /* k1 end 2 eme*/
 	
@@ -6658,11 +6876,12 @@ void printinggnuplot(char fileresu[], ch
   /*3eme*/
   for (k1=1; k1<= m ; k1 ++){
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
 
       for (cpt=1; cpt<= nlstate ; cpt ++) {
-	fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
+	fprintf(ficgp,"\n\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
+	strcpy(gplotlabel,"(");
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -6670,10 +6889,13 @@ void printinggnuplot(char fileresu[], ch
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6683,8 +6905,9 @@ void printinggnuplot(char fileresu[], ch
 	/*       k=2+nlstate*(2*cpt-2); */
 	k=2+(nlstate+1)*(cpt-1);
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
+	fprintf(ficgp,"set label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel);
 	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);
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),nres-1,nres-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);
@@ -6694,12 +6917,13 @@ 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(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(fileresu,"E_"),nres-1,nres-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(fileresu,"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_"),nres-1,nres-1,k+nlstate,cpt);
       }
+      fprintf(ficgp,"\nunset label;\n");
     } /* end nres */
   } /* end kl 3eme */
   
@@ -6707,9 +6931,10 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
+	strcpy(gplotlabel,"(");
 	fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
@@ -6718,10 +6943,13 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6729,6 +6957,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	}
       
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+	fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3;
@@ -6744,7 +6973,7 @@ set ter svg size 640, 480\nunset log y\n
 	    fprintf(ficgp,"+$%d",k+l+j-1);
 	  fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
 	} /* nlstate */
-	fprintf(ficgp,"\nset out\n");
+	fprintf(ficgp,"\nset out; unset label;\n");
       } /* end cpt state*/ 
     } /* end nres */
   } /* end covariate k1 */  
@@ -6753,9 +6982,10 @@ set ter svg size 640, 480\nunset log y\n
   /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+	strcpy(gplotlabel,"(");
 	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);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
@@ -6764,10 +6994,13 @@ set ter svg size 640, 480\nunset log y\n
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6775,6 +7008,7 @@ set ter svg size 640, 480\nunset log y\n
 	}
       
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+	fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3;
@@ -6798,7 +7032,7 @@ set ter svg size 640, 480\nunset log y\n
 	  else
 	    fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
 	}
-	fprintf(ficgp,"\nset out\n");
+	fprintf(ficgp,"\nset out; unset label;\n");
       } /* end cpt state*/ 
     } /* end covariate */  
   } /* end nres */
@@ -6807,10 +7041,10 @@ set ter svg size 640, 480\nunset log y\n
   /* CV preval stable (period) for each covariate */
   for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-    if(TKresult[nres]!= k1)
+    if(m != 1 && TKresult[nres]!= k1)
       continue;
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-      
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */
+      strcpy(gplotlabel,"(");      
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
@@ -6819,10 +7053,13 @@ set ter svg size 640, 480\nunset log y\n
 	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	vlv= nbcode[Tvaraff[k]][lv];
 	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
       }
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }	
+      strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){
 	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6830,21 +7067,22 @@ set ter svg size 640, 480\nunset log y\n
       }
       
       fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+      fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */
-      for (i=1; i<= nlstate ; i ++){
+      for (i=1; i<= nlstate ; i ++){ /* State of origin */
 	if(i==1)
 	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
 	else
 	  fprintf(ficgp,", '' ");
-	l=(nlstate+ndeath)*(i-1)+1;
+	l=(nlstate+ndeath)*(i-1)+1; /* 1, 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
 	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
 	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,"\nset out\n");
+      fprintf(ficgp,"\nset out; unset label;\n");
     } /* end cpt state*/ 
   } /* end covariate */  
   
@@ -6854,10 +7092,11 @@ set ter svg size 640, 480\nunset log y\n
     /* CV back preval stable (period) for each covariate */
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
-      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */
+	strcpy(gplotlabel,"(");      
+ 	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -6865,10 +7104,13 @@ set ter svg size 640, 480\nunset log y\n
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6876,25 +7118,26 @@ set ter svg size 640, 480\nunset log y\n
 	}
 	
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+	fprintf(ficgp,"set label \"Ending alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3; /* Offset */
-	for (i=1; i<= nlstate ; i ++){
+	for (i=1; i<= nlstate ; i ++){ /* State of origin */
 	  if(i==1)
 	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
 	  else
 	    fprintf(ficgp,", '' ");
 	  /* l=(nlstate+ndeath)*(i-1)+1; */
-	  l=(nlstate+ndeath)*(cpt-1)+1;
+	  l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
 	  /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
 	  /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
-	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
+	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */
 	  /* for (j=2; j<= nlstate ; j ++) */
 	  /* 	fprintf(ficgp,"+$%d",k+l+j-1); */
 	  /* 	/\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
 	  fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
 	} /* nlstate */
-	fprintf(ficgp,"\nset out\n");
+	fprintf(ficgp,"\nset out; unset label;\n");
       } /* end cpt state*/ 
     } /* end covariate */  
   } /* End if backcast */
@@ -6905,9 +7148,10 @@ set ter svg size 640, 480\nunset log y\n
     
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+	strcpy(gplotlabel,"(");      
 	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
@@ -6916,10 +7160,13 @@ set ter svg size 640, 480\nunset log y\n
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv];
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
+	strcpy(gplotlabel+strlen(gplotlabel),")");
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -6928,6 +7175,7 @@ set ter svg size 640, 480\nunset log y\n
 	
 	fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+	fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
@@ -6990,7 +7238,7 @@ set ter svg size 640, 480\nunset log y\n
 	    }
 	  } /* end if covariate */
 	} /* nlstate */
-	fprintf(ficgp,"\nset out\n");
+	fprintf(ficgp,"\nset out; unset label;\n");
       } /* end cpt state*/
     } /* end covariate */
   } /* End if prevfcast */
@@ -7033,17 +7281,31 @@ set ter svg size 640, 480\nunset log y\n
     fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
     fprintf(ficgp,"#model=%s \n",model);
     fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
-    fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
-    for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
+    fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
+    for(k1=1; k1 <=m; k1++)  /* For each combination of covariate */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= jk)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
-      fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
+      fprintf(ficgp,"\n\n# Combination of dummy  k1=%d which is ",k1);
+      strcpy(gplotlabel,"(");
+      sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);
+      for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
+	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+	/* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+	/* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+	vlv= nbcode[Tvaraff[k]][lv];
+	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
+      }
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }	
+      strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
-      fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres);
+      fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres);
+      fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel);
       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 */
@@ -7093,7 +7355,7 @@ set ter svg size 640, 480\nunset log y\n
 		    fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
 		  }else{ /* quantitative */
 		    fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
-		    /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		    /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
 		  }
 		  ij++;
 		}
@@ -7102,10 +7364,10 @@ set ter svg size 640, 480\nunset log y\n
 		if(ijp <=cptcovprod) { /* Product */
 		  if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
 		    if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
-		      /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],nbcode[Tvard[ijp][2]][codtabm(jk,j)]); */
+		      /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
 		      fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
 		    }else{ /* Vn is dummy and Vm is quanti */
-		      /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(jk,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+		      /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
 		      fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
 		    }
 		  }else{ /* Vn*Vm Vn is quanti */
@@ -7118,12 +7380,12 @@ set ter svg size 640, 480\nunset log y\n
 		  ijp++;
   		}
 	      } else{  /* simple covariate */
-		/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
+		/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
 		if(Dummy[j]==0){
 		  fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
 		}else{ /* quantitative */
 		  fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* */
-		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
 		}
 	      } /* end simple */
 	    } /* end j */
@@ -7136,23 +7398,23 @@ set ter svg size 640, 480\nunset log y\n
 	  if(ng != 1){
 	    fprintf(ficgp,")/(1");
 	    
-	    for(k1=1; k1 <=nlstate; k1++){ 
+	    for(cpt=1; cpt <=nlstate; cpt++){ 
 	      if(nagesqr==0)
-		fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+		fprintf(ficgp,"+exp(p%d+p%d*x",k3+(cpt-1)*ncovmodel,k3+(cpt-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);
+		fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(cpt-1)*ncovmodel,k3+(cpt-1)*ncovmodel+1,k3+(cpt-1)*ncovmodel+1+nagesqr);
 	       
 	      ij=1;
 	      for(j=3; j <=ncovmodel-nagesqr; j++){
 		if((j-2)==Tage[ij]) { /* Bug valgrind */
 		  if(ij <=cptcovage) { /* 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])]); */
+		    fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
+		    /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,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)]);/* Valgrind bug nbcode */
+		  fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
 	      }
 	      fprintf(ficgp,")");
 	    }
@@ -7172,8 +7434,8 @@ set ter svg size 640, 480\nunset log y\n
 	  i=i+ncovmodel;
 	} /* end k */
       } /* end k2 */
-      fprintf(ficgp,"\n set out\n");
-    } /* end jk */
+      fprintf(ficgp,"\n set out; unset label;\n");
+    } /* end k1 */
   } /* end ng */
   /* avoid: */
   fflush(ficgp); 
@@ -7398,7 +7660,7 @@ set ter svg size 640, 480\nunset log y\n
 /* 	      if (h==(int)(YEARM*yearp)){ */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */
   for(k=1; k<=i1;k++){
-    if(TKresult[nres]!= k)
+    if(i1 != 1 && TKresult[nres]!= k)
       continue;
     if(invalidvarcomb[k]){
       printf("\nCombination (%d) projection ignored because no cases \n",k); 
@@ -8349,6 +8611,11 @@ int decoderesult ( char resultline[], in
   if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
   }
+  if(j == 0){ /* Resultline but no = */
+    TKresult[nres]=0; /* Combination for the nresult and the model */
+    return (0);
+  }
+    
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
     printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
     fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
@@ -8963,16 +9230,16 @@ int calandcheckages(int imx, int maxwav,
 	*nberr = *nberr + 1;
 	if(firstone == 0){
 	  firstone=1;
-	printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\nOther similar cases in log file\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
+	printf("Warning (#%d)! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown but status is a death state %d at wave %d. If you don't know the vital status, please enter -2. If he/she is still alive but don't know the state, please code with '-1 or '.'. Here, we do not believe in a death, skipped.\nOther similar cases in log file\n", *nberr,(int)moisdc[i],(int)andc[i],num[i],i,s[m][i],m);
 	}
-	fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
-	s[m][i]=-1;
+	fprintf(ficlog,"Warning (#%d)! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown but status is a death state %d at wave %d. If you don't know the vital status, please enter -2. If he/she is still alive but don't know the state, please code with '-1 or '.'. Here, we do not believe in a death, skipped.\n", *nberr,(int)moisdc[i],(int)andc[i],num[i],i,s[m][i],m);
+	s[m][i]=-1;  /* Droping the death status */
       }
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
 	(*nberr)++;
-	printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
-	fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
-	s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
+	printf("Error (#%d)! Month of death of individual %ld on line %d was unknown (%2d) (year of death is %4d) and status is a death state %d at wave %d. Please impute an arbitrary (or not) month and rerun. Currently this transition to death will be skipped (status is set to -2).\nOther similar cases in log file\n", *nberr, num[i],i,(int)moisdc[i],(int)andc[i],s[m][i],m);
+	fprintf(ficlog,"Error (#%d)! Month of death of individual %ld on line %d was unknown (%2d) (year of death is %4d) and status is a death state %d at wave %d. Please impute an arbitrary (or not) month and rerun. Currently this transition to death will be skipped (status is set to -2).\n", *nberr, num[i],i,(int)moisdc[i],(int)andc[i],s[m][i],m);
+	s[m][i]=-2; /* We prefer to skip it (and to skip it in version 0.8a1 too */
       }
     }
   }
@@ -9285,7 +9552,7 @@ int prevalence_limit(double *p, double *
 
   for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
 
       /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
@@ -9382,7 +9649,7 @@ int back_prevalence_limit(double *p, dou
   
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+     if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
       fprintf(ficresplb,"#******");
@@ -9441,6 +9708,7 @@ int back_prevalence_limit(double *p, dou
 	fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
       } /* Age */
       /* was end of cptcod */
+      /*fprintf(ficresplb,"\n");*/ /* Seems to be necessary for gnuplot only if two result lines and no covariate. */
     } /* end of any combination */
   } /* end of nres */  
   /* hBijx(p, bage, fage); */
@@ -9485,7 +9753,7 @@ int hPijx(double *p, int bage, int fage)
 		/*  	k=k+1;  */
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
@@ -9557,14 +9825,14 @@ int hPijx(double *p, int bage, int fage)
   
   /* hstepm=1;   aff par mois*/
   pstamp(ficrespijb);
-  fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
+  fprintf(ficrespijb,"#****** h Bij x Back probability to be in state i at age x-h being in j at x: B1j+B2j+...=1 ");
   i1= pow(2,cptcoveff);
   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
   /*  	k=k+1;  */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficrespijb,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)
@@ -9573,7 +9841,7 @@ int hPijx(double *p, int bage, int fage)
 	fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
       }
       fprintf(ficrespijb,"******\n");
-      if(invalidvarcomb[k]){
+      if(invalidvarcomb[k]){  /* Is it necessary here? */
 	fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); 
 	continue;
       }
@@ -9591,7 +9859,7 @@ int hPijx(double *p, int bage, int fage)
 	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
 	hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
 	/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
-	fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
+	fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
 	for(i=1; i<=nlstate;i++)
 	  for(j=1; j<=nlstate+ndeath;j++)
 	    fprintf(ficrespijb," %1d-%1d",i,j);
@@ -9637,6 +9905,7 @@ int main(int argc, char *argv[])
   int NDIM=2;
   int vpopbased=0;
   int nres=0;
+  int endishere=0;
 
   char ca[32], cb[32];
   /*  FILE *fichtm; *//* Html File */
@@ -9918,9 +10187,12 @@ int main(int argc, char *argv[])
       break;
   }
   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
-    if (num_filled == 0)
-            model[0]='\0';
-    else if (num_filled != 1){
+    if (num_filled == 0){
+      printf("ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);
+      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);
+      model[0]='\0';
+      goto end;
+    } 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';
@@ -10007,11 +10279,6 @@ int main(int argc, char *argv[])
     fclose (ficlog);
     goto end;
     exit(0);
-  } else if(mle==-2) { /* Guessing from means */
-    prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
-    printf(" You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
-    fprintf(ficlog," You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
-   
   }  else if(mle==-5) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
@@ -10997,70 +11264,31 @@ Please run with mle=-1 to get a correct
     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
 		
     /* Other stuffs, more or less useful */    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficparo);
-    }
-    ungetc(c,ficpar);
-    
-    fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
-    fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficparo);
-    }
-    ungetc(c,ficpar);
-    
-    
-    dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
-    dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
-    
-    fscanf(ficpar,"pop_based=%d\n",&popbased);
-    fprintf(ficlog,"pop_based=%d\n",popbased);
-    fprintf(ficparo,"pop_based=%d\n",popbased);   
-    fprintf(ficres,"pop_based=%d\n",popbased);   
-    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficres);
-      fputs(line,ficparo);
+    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;
     }
-    ungetc(c,ficpar);
-    
-    fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);
-    fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-    printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-    fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-    fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-    /* day and month of proj2 are not used but only year anproj2.*/
-    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficparo);
-      fputs(line,ficres);
+
+    if((num_filled=sscanf(line,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav)) !=EOF){
+      
+      if (num_filled != 7) {
+	printf("Error: Not 7 (data)parameters in line but %d, for example:begin-prev-date=1/1/1990 end-prev-date=1/6/2004  mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	fprintf(ficlog,"Error: Not 7 (data)parameters in line but %d, for example:begin-prev-date=1/1/1990 end-prev-date=1/6/2004  mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	goto end;
+      }
+      printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
     }
-    ungetc(c,ficpar);
-    
-    fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
-    fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
-    fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
-    fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
-    /* day and month of proj2 are not used but only year anproj2.*/
-    
-    /* Results */
-    nresult=0;
+
     while(fgets(line, MAXLINE, ficpar)) {
       /* If line starts with a # it is a comment */
       if (line[0] == '#') {
@@ -11068,49 +11296,114 @@ Please run with mle=-1 to get a correct
 	fputs(line,stdout);
 	fputs(line,ficparo);
 	fputs(line,ficlog);
-	fputs(line,ficres);
 	continue;
       }else
 	break;
     }
-    if (!feof(ficpar))
-    while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
-      if (num_filled == 0){
-	resultline[0]='\0';
-      break;
-      } else if (num_filled != 1){
-	printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
-      }
-      nresult++; /* Sum of resultlines */
-      printf("Result %d: result=%s\n",nresult, resultline);
-      if(nresult > MAXRESULTLINES){
-	printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
-	fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+    
+    
+    dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
+    dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
+    
+    if((num_filled=sscanf(line,"pop_based=%d\n",&popbased)) !=EOF){
+      if (num_filled != 1) {
+	printf("Error: Not 1 (data)parameters in line but %d, for example:pop_based=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	fprintf(ficlog,"Error: Not 1 (data)parameters in line but %d, for example: pop_based=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
 	goto end;
       }
-      decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
-      fprintf(ficparo,"result: %s\n",resultline);
-      fprintf(ficres,"result: %s\n",resultline);
-      fprintf(ficlog,"result: %s\n",resultline);
-      while(fgets(line, MAXLINE, ficpar)) {
+      printf("pop_based=%d\n",popbased);
+      fprintf(ficlog,"pop_based=%d\n",popbased);
+      fprintf(ficparo,"pop_based=%d\n",popbased);   
+      fprintf(ficres,"pop_based=%d\n",popbased);   
+    }
+     
+    /* Results */
+    nresult=0;
+    do{
+      if(!fgets(line, MAXLINE, ficpar)){
+	endishere=1;
+	parameterline=14;
+      }else if (line[0] == '#') {
 	/* If line starts with a # it is a comment */
-	if (line[0] == '#') {
-	  numlinepar++;
-	  fputs(line,stdout);
-	  fputs(line,ficparo);
-	  fputs(line,ficres);
-	  fputs(line,ficlog);
-	  continue;
-	}else
-	  break;
+	numlinepar++;
+	fputs(line,stdout);
+	fputs(line,ficparo);
+	fputs(line,ficlog);
+	continue;
+      }else if(sscanf(line,"prevforecast=%[^\n]\n",modeltemp))
+	parameterline=11;
+      else if(sscanf(line,"backcast=%[^\n]\n",modeltemp))
+	parameterline=12;
+      else if(sscanf(line,"result:%[^\n]\n",modeltemp))
+	parameterline=13;
+      else{
+	parameterline=14;
       }
-      if (feof(ficpar))
+      switch (parameterline){ 
+      case 11:
+	if((num_filled=sscanf(line,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj)) !=EOF){
+	  if (num_filled != 8) {
+	    printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    goto end;
+	  }
+	  fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+	  printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+	  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+	  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+	  /* day and month of proj2 are not used but only year anproj2.*/
+	}
 	break;
-      else{ /* Processess output results for this combination of covariate values */
-      }				   
-    } /* end while */
-
-
+      case 12:
+	/*fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);*/
+	if((num_filled=sscanf(line,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj)) !=EOF){
+	  if (num_filled != 8) {
+	    printf("Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    goto end;
+	  }
+	  printf("backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+	  fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+	  fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+	  fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+	  /* day and month of proj2 are not used but only year anproj2.*/
+	}
+	break;
+      case 13:
+	if((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
+	  if (num_filled == 0){
+	    resultline[0]='\0';
+	    printf("Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
+	    fprintf(ficlog,"Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
+	    break;
+	  } else if (num_filled != 1){
+	    printf("ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
+	    fprintf(ficlog,"ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
+	  }
+	  nresult++; /* Sum of resultlines */
+	  printf("Result %d: result=%s\n",nresult, resultline);
+	  if(nresult > MAXRESULTLINES){
+	    printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	    fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	    goto end;
+	  }
+	  decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+	  fprintf(ficparo,"result: %s\n",resultline);
+	  fprintf(ficres,"result: %s\n",resultline);
+	  fprintf(ficlog,"result: %s\n",resultline);
+	  break;
+	case 14: 
+	  if(ncovmodel >2 && nresult==0 ){
+	    printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
+	    goto end;
+	  }
+	  break;
+	default:
+	  nresult=1;
+	  decoderesult(".",nresult ); /* No covariate */
+	}
+      } /* End switch parameterline */
+    }while(endishere==0); /* End do */
     
     /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
@@ -11127,7 +11420,7 @@ Please run with mle=-1 to get a correct
       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p);
     }
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
-		 model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
+		 model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
 		 jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
 		
     /*------------ free_vector  -------------*/
@@ -11180,6 +11473,7 @@ Please run with mle=-1 to get a correct
       mobaverage=mobaverages;
       if (mobilav!=0) {
 	printf("Movingaveraging observed prevalence\n");
+	fprintf(ficlog,"Movingaveraging observed prevalence\n");
 	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
 	  fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
 	  printf(" Error in movingaverage mobilav=%d\n",mobilav);
@@ -11189,6 +11483,7 @@ Please run with mle=-1 to get a correct
       /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
       else if (mobilavproj !=0) {
 	printf("Movingaveraging projected observed prevalence\n");
+	fprintf(ficlog,"Movingaveraging projected observed prevalence\n");
 	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
 	  fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
 	  printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
@@ -11251,7 +11546,7 @@ Please run with mle=-1 to get a correct
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficreseij,"\n#****** ");
       printf("\n#****** ");
@@ -11324,7 +11619,7 @@ Please run with mle=-1 to get a correct
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       printf("\n#****** Result for:");
       fprintf(ficrest,"\n#****** Result for:");
@@ -11465,7 +11760,7 @@ Please run with mle=-1 to get a correct
 
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficresvpl,"\n#****** ");
       printf("\n#****** ");