--- imach/src/imach.c	2017/04/03 10:17:47	1.258
+++ imach/src/imach.c	2017/04/26 06:01:29	1.264
@@ -1,6 +1,24 @@
-/* $Id: imach.c,v 1.258 2017/04/03 10:17:47 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
 
@@ -976,12 +994,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.258 2017/04/03 10:17:47 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.258 $ $Date: 2017/04/03 10:17:47 $"; 
+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  */
@@ -2600,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 */
@@ -2661,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])]; *\/ */
@@ -4604,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]);*/
@@ -4656,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);
@@ -4938,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;
 	}
@@ -5019,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 */
@@ -6419,11 +6440,47 @@ 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 */
@@ -6434,6 +6491,15 @@ void printinghtml(char fileresu[], char
      /* 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]);
@@ -6455,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> \
@@ -6479,13 +6545,13 @@ 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 some years earlier, knowing that we will be 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 mixed (stable) back prevalence in state %d. Or probability to be in state %d at a younger age, knowing that we will be in state (1 to %d) at different older 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);
        }
      }
@@ -6602,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;
@@ -6661,6 +6727,7 @@ void printinggnuplot(char fileresu[], ch
 	/* 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 */
@@ -6670,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)");
@@ -6729,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 */
@@ -6741,6 +6814,7 @@ void printinggnuplot(char fileresu[], ch
       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 */
@@ -6748,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); 
@@ -6762,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)");
@@ -6790,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*/
 	
@@ -6802,7 +6880,8 @@ void printinggnuplot(char fileresu[], ch
 	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 */
@@ -6810,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); 
@@ -6823,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);
@@ -6834,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 */
   
@@ -6850,6 +6934,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
       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 */
@@ -6858,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); 
@@ -6869,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;
@@ -6884,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 */  
@@ -6896,6 +6985,7 @@ set ter svg size 640, 480\nunset log y\n
       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 */
@@ -6904,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); 
@@ -6915,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;
@@ -6938,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 */
@@ -6950,7 +7044,7 @@ set ter svg size 640, 480\nunset log y\n
     if(m != 1 && TKresult[nres]!= k1)
       continue;
     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 */
@@ -6959,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); 
@@ -6970,6 +7067,7 @@ 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 */
@@ -6984,7 +7082,7 @@ set ter svg size 640, 480\nunset log y\n
 	  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 */  
   
@@ -6997,7 +7095,8 @@ set ter svg size 640, 480\nunset log y\n
       if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */
-	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
+	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 */
@@ -7005,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); 
@@ -7016,6 +7118,7 @@ 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 */
@@ -7034,7 +7137,7 @@ set ter svg size 640, 480\nunset log y\n
 	  /* 	/\* 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 */
@@ -7048,6 +7151,7 @@ set ter svg size 640, 480\nunset log y\n
       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 */
@@ -7056,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); 
@@ -7068,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 */
@@ -7130,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 */
@@ -7173,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(m != 1 && 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 */
@@ -7233,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++;
 		}
@@ -7242,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 */
@@ -7258,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 */
@@ -7276,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,")");
 	    }
@@ -7312,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); 
@@ -9108,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 */
       }
     }
   }
@@ -9719,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;
       }
@@ -10065,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';
@@ -11233,8 +11358,8 @@ Please run with mle=-1 to get a correct
 	/*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 finloal-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 finloal-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    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);
@@ -11268,10 +11393,11 @@ Please run with mle=-1 to get a correct
 	  fprintf(ficlog,"result: %s\n",resultline);
 	  break;
 	case 14: 
-	  if(ncovmodel >2){
-	    printf("ERROR: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
+	  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 */