--- imach/src/imach.c	2024/04/24 21:21:17	1.359
+++ imach/src/imach.c	2024/04/30 10:59:22	1.360
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $
+/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.360  2024/04/30 10:59:22  brouard
+  Summary: Version 0.99s4 and estimation of std of e.j/e..
+
   Revision 1.359  2024/04/24 21:21:17  brouard
   Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
 
@@ -1390,12 +1393,12 @@ double gnuplotversion=GNUPLOTVERSION;
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ */
+/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
-char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.359 $ $Date: 2024/04/24 21:21:17 $"; 
+char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
+char fullversion[]="$Revision: 1.360 $ $Date: 2024/04/30 10:59:22 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -8812,7 +8815,7 @@ void  concatwav(int wav[], int **dh, int
      }
 		
      /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
-      * p.j overall mortality formula 49 but computed directly because
+      * p.j overall mortality formula 19 but computed directly because
       * we compute the grad (wix pijx) instead of grad (pijx),even if
       * wix is independent of theta.
       */
@@ -9910,7 +9913,10 @@ prevalence (with 95%% confidence interva
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \
-health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \
+health expectancies in each live state (1 to %d) with confidence intervals \
+on left y-scale as well as proportions of time spent in each live state \
+(with confidence intervals) on right y-scale 0 to 100%%.\
+ If popbased=1 the smooth (due to the model)				\
 true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\
  observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
@@ -10264,18 +10270,18 @@ void printinggnuplot(char fileresu[], ch
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	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);
+	  fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage);
 	}else
 	  fprintf(ficgp,"\nreplot ");
-	for (i=1; i<= nlstate+1 ; i ++) {
+	for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
 	  k=2*i;
-	  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)");
+	  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 fixed variables age, popbased, mobilav */
+	  for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
+	    else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
 	  }   
 	  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);
+	  else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
 	  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)");
@@ -10287,9 +10293,39 @@ void printinggnuplot(char fileresu[], ch
 	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
 	    else fprintf(ficgp," %%*lf (%%*lf)");
 	  }   
-	  if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+	  if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */
 	  else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
 	} /* state */
+	/* again for the percentag spent in state i-1=1 to i-1=nlstate */
+	for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
+	  k=2*i;
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d &&  ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
+	  for (j=1; j<= nlstate ; j ++)
+	    fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+	  for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
+	    else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
+	  }   
+	  if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */
+	  else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+	  for (j=1; j<= nlstate ; j ++)
+	    fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+	  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 axis x1y2,");
+	  fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+	  for (j=1; j<= nlstate ; j ++)
+	    fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+	  for (j=1; j<= nlstate+1 ; j ++) {
+	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	    else fprintf(ficgp," %%*lf (%%*lf)");
+	  }   
+	  if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2");
+	  else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n");
+	} /* state for percent */
       } /* vpopbased */
       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 */
@@ -16740,16 +16776,19 @@ Please run with mle=-1 to get a correct
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	oldm=oldms;savm=savms; /* ZZ Segmentation fault */
 	cptcod= 0; /* To be deleted */
-	printf("varevsij vpopbased=%d \n",vpopbased);
-	fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);
+	printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
+	fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
 	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */
-	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
+	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\
+#  (these are weighted average of eij where weights are ");
 	if(vpopbased==1)
-	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally)\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
 	else
-	  fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
+	  fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n");
+	fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n ");
 	fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
 	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+	for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (std) ",i);
 	fprintf(ficrest,"\n");
 	/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
 	printf("Computing age specific forward period (stable) prevalences in each health state \n");
@@ -16783,9 +16822,19 @@ Please run with mle=-1 to get a correct
 	    for(j=1;j <=nlstate;j++)
 	      vepp += vareij[i][j][(int)age];
 	  fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+	  /* vareij[j][i] is the variance of epj */
 	  for(j=1;j <=nlstate;j++){
 	    fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
 	  }
+	  /* And proportion of time spent in state j */
+	  /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */
+	  /* \sigma^2_x/\mu_y^2 +\sigma^2_y \mu^2x/\mu_y^4 */
+	  /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp */
+	  /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstata+1]^4 */
+	  for(j=1;j <=nlstate;j++){
+	    /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
+	    fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] ));
+	  }
 	  fprintf(ficrest,"\n");
 	}
       } /* End vpopbased */