--- imach/src/imach.c	2002/05/07 15:53:01	1.41
+++ imach/src/imach.c	2003/06/13 07:45:28	1.41.2.2
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.41 2002/05/07 15:53:01 lievre Exp $
+/* $Id: imach.c,v 1.41.2.2 2003/06/13 07:45:28 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -60,7 +60,8 @@
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80
 /*#define DEBUG*/
-#define windows
+
+/*#define windows*/
 #define	GLOCK_ERROR_NOPATH		-1	/* empty path */
 #define	GLOCK_ERROR_GETCWD		-2	/* cannot get cwd */
 
@@ -869,6 +870,7 @@ double func( double *x)
   double **out;
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
+  int s1, s2;
   long ipmx;
   /*extern weight */
   /* We are differentiating ll according to initial status */
@@ -883,7 +885,10 @@ double func( double *x)
     for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
     for(mi=1; mi<= wav[i]-1; mi++){
       for (ii=1;ii<=nlstate+ndeath;ii++)
-	for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	for (j=1;j<=nlstate+ndeath;j++){
+	  oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+	  savm[ii][j]=(ii==j ? 1.0 : 0.0);
+	}
       for(d=0; d<dh[mi][i]; d++){
 	newm=savm;
 	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
@@ -899,17 +904,47 @@ double func( double *x)
 	
       } /* end mult */
       
-      lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
-      /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
+      s1=s[mw[mi][i]][i];
+      s2=s[mw[mi+1][i]][i];
+      if( s2 > nlstate){ 
+	/* i.e. if s2 is a death state and if the date of death is known then the contribution
+	   to the likelihood is the probability to die between last step unit time and current 
+	   step unit time, which is also the differences between probability to die before dh 
+	   and probability to die before dh-stepm . 
+	   In version up to 0.92 likelihood was computed
+	   as if date of death was unknown. Death was treated as any other
+	   health state: the date of the interview describes the actual state
+	   and not the date of a change in health state. The former idea was
+	   to consider that at each interview the state was recorded
+	   (healthy, disable or death) and IMaCh was corrected; but when we
+	   introduced the exact date of death then we should have modified
+	   the contribution of an exact death to the likelihood. This new
+	   contribution is smaller and very dependent of the step unit
+	   stepm. It is no more the probability to die between last interview
+	   and month of death but the probability to survive from last
+	   interview up to one month before death multiplied by the
+	   probability to die within a month. Thanks to Chris
+	   Jackson for correcting this bug.  Former versions increased
+	   mortality artificially. The bad side is that we add another loop
+	   which slows down the processing. The difference can be up to 10%
+	   lower mortality.
+	*/
+	lli=log(out[s1][s2] - savm[s1][s2]);
+      }else{
+	lli=log(out[s1][s2]); /* or	lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); */
+	/* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
+      }
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d lli=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],lli,weight[i],out[s1][s2],savm[s1][s2]);*/
     } /* end of wave */
   } /* end of individual */
 
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+  /*exit(0);*/
   return -l;
 }
 
@@ -2158,12 +2193,7 @@ m=pow(2,cptcoveff);
   for (cpt=1; cpt<= nlstate ; cpt ++) {
    for (k1=1; k1<= m ; k1 ++) {
 
-#ifdef windows
-    fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
-#endif
-#ifdef unix
-fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
-#endif
+     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
 
 for (i=1; i<= nlstate ; i ++) {
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
@@ -2180,9 +2210,7 @@ for (i=1; i<= nlstate ; i ++) {
   else fprintf(ficgp," \%%*lf (\%%*lf)");
 }  
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
-#ifdef unix
-fprintf(ficgp,"\nset ter gif small size 400,300");
-#endif
+
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
    }
   }
@@ -2651,7 +2679,7 @@ int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   
 
-  char version[80]="Imach version 0.8a, May 2002, INED-EUROREVES ";
+  char version[80]="Imach version 0.8a1, June 2003, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
 
@@ -3526,9 +3554,7 @@ free_matrix(mint,1,maxwav,1,n);
 
 
  end:
-#ifdef windows
   /* chdir(pathcd);*/
-#endif 
  /*system("wgnuplot graph.plt");*/
  /*system("../gp37mgw/wgnuplot graph.plt");*/
  /*system("cd ../gp37mgw");*/
@@ -3538,7 +3564,7 @@ free_matrix(mint,1,maxwav,1,n);
  strcat(plotcmd,optionfilegnuplot);
  system(plotcmd);
 
-#ifdef windows
+ /*#ifdef windows*/
   while (z[0] != 'q') {
     /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
@@ -3548,7 +3574,7 @@ free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }
-#endif 
+  /*#endif */
 }