--- imach/src/imach.c	2003/06/13 07:45:28	1.41.2.2
+++ imach/src/imach.c	2002/05/21 18:44:41	1.42
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.41.2.2 2003/06/13 07:45:28 brouard Exp $
+/* $Id: imach.c,v 1.42 2002/05/21 18:44:41 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -56,12 +56,11 @@
 #include <unistd.h>
 
 #define MAXLINE 256
-#define GNUPLOTPROGRAM "wgnuplot"
+#define GNUPLOTPROGRAM "gnuplot"
 /*#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 */
 
@@ -870,7 +869,6 @@ 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 */
@@ -885,10 +883,7 @@ 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);
-	  savm[ii][j]=(ii==j ? 1.0 : 0.0);
-	}
+	for (j=1;j<=nlstate+ndeath;j++) oldm[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;
@@ -904,47 +899,17 @@ double func( double *x)
 	
       } /* end mult */
       
-      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]]);*/
-      }
+      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;
 }
 
@@ -1365,10 +1330,10 @@ void prevalence(int agemin, float agemax
   j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
- for(k1=1; k1<=j;k1++){
+  for(k1=1; k1<=j;k1++){
     for(i1=1; i1<=ncodemax[k1];i1++){
       j1++;
- 
+      
       for (i=-1; i<=nlstate+ndeath; i++)  
 	for (jk=-1; jk<=nlstate+ndeath; jk++)  
 	  for(m=agemin; m <= agemax+3; m++)
@@ -1387,41 +1352,43 @@ void prevalence(int agemin, float agemax
 	    if ((k2>=dateprev1) && (k2<=dateprev2)) {
 	      if(agev[m][i]==0) agev[m][i]=agemax+1;
 	      if(agev[m][i]==1) agev[m][i]=agemax+2;
-	      if (m<lastpass) 
-	      	if (calagedate>0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
-	      else
-	       freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
-	       freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; 
+	      if (m<lastpass) {
+	      	if (calagedate>0) 
+		  freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
+		else
+		  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
+		freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; 
+	      }
 	    }
 	  }
 	}
       }
-	for(i=(int)agemin; i <= (int)agemax+3; i++){ 
-	  for(jk=1; jk <=nlstate ; jk++){
-	    for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-	      pp[jk] += freq[jk][m][i]; 
-	  }
-	  for(jk=1; jk <=nlstate ; jk++){
-	    for(m=-1, pos=0; m <=0 ; m++)
+      for(i=(int)agemin; i <= (int)agemax+3; i++){ 
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+	    pp[jk] += freq[jk][m][i]; 
+	}
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pos=0; m <=0 ; m++)
 	    pos += freq[jk][m][i];
 	}
 	
-	 for(jk=1; jk <=nlstate ; jk++){
-	   for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
-	     pp[jk] += freq[jk][m][i];
-	 }
-	 
-	 for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
-
-	 for(jk=1; jk <=nlstate ; jk++){	   
-	   if( i <= (int) agemax){
-	     if(pos>=1.e-5){
-	       probs[i][jk][j1]= pp[jk]/pos;
-	     }
-	   }
-	 }
-	 
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
+	    pp[jk] += freq[jk][m][i];
+	}
+	
+	for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
+	
+	for(jk=1; jk <=nlstate ; jk++){	   
+	  if( i <= (int) agemax){
+	    if(pos>=1.e-5){
+	      probs[i][jk][j1]= pp[jk]/pos;
+	    }
+	  }
 	}
+	
+      }
     }
   }
 
@@ -1698,6 +1665,7 @@ void evsij(char fileres[], double ***eij
       for(j=1;j<=nlstate*2;j++)
 	varhe[i][j][(int)age] =0.;
 
+     printf("%d||",(int)age);fflush(stdout);
     for(h=0;h<=nhstepm-1;h++){
       for(k=0;k<=nhstepm-1;k++){
 	matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
@@ -2149,24 +2117,24 @@ fprintf(fichtm," <li>Graphs</li><p>");
 	   fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
 	 fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        }
-       fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
-<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
+       fprintf(fichtm,"<br>- Probabilities: pe%s%d.png<br>
+<img src=\"pe%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
        for(cpt=1; cpt<nlstate;cpt++){
-	 fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
-<img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
+	 fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.png<br>
+<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }
     for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
-interval) in state (%d): v%s%d%d.gif <br>
-<img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
+interval) in state (%d): v%s%d%d.png <br>
+<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
-        fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
-<img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
+        fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
+<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and
-health expectancies in states (1) and (2): e%s%d.gif<br>
-<img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
+health expectancies in states (1) and (2): e%s%d.png<br>
+<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
 fprintf(fichtm,"\n</body>");
    }
    }
@@ -2193,7 +2161,14 @@ m=pow(2,cptcoveff);
   for (cpt=1; cpt<= nlstate ; cpt ++) {
    for (k1=1; k1<= m ; k1 ++) {
 
-     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);
+#ifdef windows
+     fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
+#endif
+#ifdef unix
+fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
+#endif
 
 for (i=1; i<= nlstate ; i ++) {
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
@@ -2210,14 +2185,16 @@ 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));
-
-fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
+#ifdef unix
+fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
+#endif
    }
   }
   /*2 eme*/
 
   for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);
+    fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1);
+    fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
     
     for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;
@@ -2242,7 +2219,6 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
       else fprintf(ficgp,"\" t\"\" w l 0,");
     }
-    fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
   }
  
   /*3eme*/
@@ -2250,7 +2226,8 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       k=2+nlstate*(2*cpt-2);
-      fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
+      fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+      fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-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);
@@ -2263,15 +2240,15 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
 	fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);
 
       } 
-      fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
-    }
     }
+  }
  
   /* CV preval stat */
     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<nlstate ; cpt ++) {
       k=3;
-      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
+      fprintf(ficgp,"set out \"p%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
 
       for (i=1; i< nlstate ; i ++)
 	fprintf(ficgp,"+$%d",k+i+1);
@@ -2284,7 +2261,6 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
 	fprintf(ficgp,"+$%d",l+i+1);
       }
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);   
-      fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
     } 
   }  
   
@@ -2300,50 +2276,50 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
 	}
       }
     }
-    }
+   }
 
-    for(jk=1; jk <=m; jk++) {
-  fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
-   i=1;
-   for(k2=1; k2<=nlstate; k2++) {
-     k3=i;
-     for(k=1; k<=(nlstate+ndeath); k++) {
-       if (k != k2){
-    	fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-ij=1;
-	for(j=3; j <=ncovmodel; j++) {
-	  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
-	    fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-	    ij++;
-	  }
-	  else
-	  fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
-	}
-	  fprintf(ficgp,")/(1");
-	
-	for(k1=1; k1 <=nlstate; k1++){   
-	  fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
-ij=1;
-	  for(j=3; j <=ncovmodel; j++){
-	  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
-	    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-	    ij++;
-	  }
-	  else
-	    fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
-	  }
-	  fprintf(ficgp,")");
-	}
-	fprintf(ficgp,") t \"p%d%d\" ", k2,k);
-	if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
-	i=i+ncovmodel;
+   for(jk=1; jk <=m; jk++) {
+     fprintf(ficgp,"\nset out \"pe%s%d.png\" \n\n",strtok(optionfile, "."),jk); 
+     fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+     i=1;
+     for(k2=1; k2<=nlstate; k2++) {
+       k3=i;
+       for(k=1; k<=(nlstate+ndeath); k++) {
+	 if (k != k2){
+	   fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+	   ij=1;
+	   for(j=3; j <=ncovmodel; j++) {
+	     if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+	       fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+	       ij++;
+	     }
+	     else
+	       fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+	   }
+	   fprintf(ficgp,")/(1");
+	   
+	   for(k1=1; k1 <=nlstate; k1++){   
+	     fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+	     ij=1;
+	     for(j=3; j <=ncovmodel; j++){
+	       if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+		 fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+		 ij++;
+	       }
+	       else
+		 fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+	     }
+	     fprintf(ficgp,")");
+	   }
+	   fprintf(ficgp,") t \"p%d%d\" ", k2,k);
+	   if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
+	   i=i+ncovmodel;
+	 }
        }
      }
    }
-   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); 
-   }
    
-  fclose(ficgp);
+   fclose(ficgp);
 }  /* end gnuplot */
 
 
@@ -2679,7 +2655,7 @@ int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   
 
-  char version[80]="Imach version 0.8a1, June 2003, INED-EUROREVES ";
+  char version[80]="Imach version 0.8d, May 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
 
@@ -3554,7 +3530,9 @@ 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");*/
@@ -3564,7 +3542,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: ");
@@ -3574,7 +3552,7 @@ free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }
-  /*#endif */
+#endif 
 }