--- imach/src/imach.c	2017/04/24 15:23:15	1.263
+++ imach/src/imach.c	2017/04/26 06:01:29	1.264
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.263 2017/04/24 15:23:15 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
 
@@ -991,12 +994,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.263 2017/04/24 15:23:15 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.263 $ $Date: 2017/04/24 15:23:15 $"; 
+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  */
@@ -2615,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 */
@@ -2676,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])]; *\/ */
@@ -4619,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]);*/
@@ -6435,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 */
@@ -6450,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]);
@@ -6495,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);
        }
      }
@@ -6618,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;
@@ -6677,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 */
@@ -6686,11 +6737,14 @@ 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]){
@@ -6701,6 +6755,7 @@ void printinggnuplot(char fileresu[], ch
       
 	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 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*/
@@ -6747,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 */
@@ -6759,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 */
@@ -6766,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); 
@@ -6780,9 +6839,10 @@ 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;
@@ -6808,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*/
 	
@@ -6821,6 +6881,7 @@ void printinggnuplot(char fileresu[], ch
 
       for (cpt=1; cpt<= nlstate ; 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 */
@@ -6828,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); 
@@ -6841,6 +6905,7 @@ 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_"),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);
@@ -6858,6 +6923,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	} 
 	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 */
   
@@ -6868,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 */
@@ -6876,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); 
@@ -6887,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;
@@ -6902,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 */  
@@ -6914,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 */
@@ -6922,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); 
@@ -6933,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;
@@ -6956,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 */
@@ -6968,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 */
@@ -6977,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); 
@@ -6988,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 */
@@ -7002,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 */  
   
@@ -7015,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 */
@@ -7023,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); 
@@ -7034,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 */
@@ -7052,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 */
@@ -7066,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 */
@@ -7074,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); 
@@ -7086,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 */
@@ -7148,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 */
@@ -7191,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 */
@@ -7251,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++;
 		}
@@ -7260,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 */
@@ -7276,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 */
@@ -7294,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,")");
 	    }
@@ -7330,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); 
@@ -9737,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;
       }