--- imach/src/imach.c	2004/05/16 15:05:56	1.98
+++ imach/src/imach.c	2004/09/15 10:38:38	1.101
@@ -1,6 +1,15 @@
-/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $
+/* $Id: imach.c,v 1.101 2004/09/15 10:38:38 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.101  2004/09/15 10:38:38  brouard
+  Fix on curr_time
+
+  Revision 1.100  2004/07/12 18:29:06  brouard
+  Add version for Mac OS X. Just define UNIX in Makefile
+
+  Revision 1.99  2004/06/05 08:57:40  brouard
+  *** empty log message ***
+
   Revision 1.98  2004/05/16 15:05:56  brouard
   New version 0.97 . First attempt to estimate force of mortality
   directly from the data i.e. without the need of knowing the health
@@ -196,7 +205,7 @@
 #include <stdlib.h>
 #include <unistd.h>
 
-#include <sys/time.h>
+/* #include <sys/time.h> */
 #include <time.h>
 #include "timeval.h"
 
@@ -224,7 +233,7 @@
 #define AGESUP 130
 #define AGEBASE 40
 #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
-#ifdef unix
+#ifdef UNIX
 #define DIRSEPARATOR '/'
 #define ODIRSEPARATOR '\\'
 #else
@@ -232,11 +241,11 @@
 #define ODIRSEPARATOR '/'
 #endif
 
-/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ */
+/* $Id: imach.c,v 1.101 2004/09/15 10:38:38 brouard Exp $ */
 /* $State: Exp $ */
 
-char version[]="Imach version 0.70, May 2004, INED-EUROREVES ";
-char fullversion[]="$Revision: 1.98 $ $Date: 2004/05/16 15:05:56 $"; 
+char version[]="Imach version 0.97b, May 2004, INED-EUROREVES ";
+char fullversion[]="$Revision: 1.101 $ $Date: 2004/09/15 10:38:38 $"; 
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -356,6 +365,9 @@ double ftolhess; /* Tolerance for comput
 /**************** split *************************/
 static	int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {
+  /* From a file name with full path (either Unix or Windows) we extract the directory (dirc)
+     the name of the file (name), its extension only (ext) and its first part of the name (finame)
+  */ 
   char	*ss;				/* pointer */
   int	l1, l2;				/* length counters */
 
@@ -387,12 +399,14 @@ static	int split( char *path, char *dirc
 #endif
   */
   ss = strrchr( name, '.' );		/* find last / */
-  ss++;
-  strcpy(ext,ss);			/* save extension */
-  l1= strlen( name);
-  l2= strlen(ss)+1;
-  strncpy( finame, name, l1-l2);
-  finame[l1-l2]= 0;
+  if (ss >0){
+    ss++;
+    strcpy(ext,ss);			/* save extension */
+    l1= strlen( name);
+    l2= strlen(ss)+1;
+    strncpy( finame, name, l1-l2);
+    finame[l1-l2]= 0;
+  }
   return( 0 );				/* we're done */
 }
 
@@ -868,11 +882,11 @@ void powell(double p[], double **xi, int
     fprintf(ficrespow,"\n");fflush(ficrespow);
     if(*iter <=3){
       tm = *localtime(&curr_time.tv_sec);
-      strcpy(strcurr,asctime(&tmf));
+      strcpy(strcurr,asctime(&tm));
 /*       asctime_r(&tm,strcurr); */
-      forecast_time=curr_time;
+      forecast_time=curr_time; 
       itmp = strlen(strcurr);
-      if(strcurr[itmp-1]=='\n')
+      if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
 	strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
@@ -884,8 +898,8 @@ void powell(double p[], double **xi, int
 	itmp = strlen(strfor);
 	if(strfor[itmp-1]=='\n')
 	strfor[itmp-1]='\0';
-	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
-	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
       }
     }
     for (i=1;i<=n;i++) { 
@@ -1050,57 +1064,57 @@ double **pmij(double **ps, double *cov,
   int i,j,j1, nc, ii, jj;
 
     for(i=1; i<= nlstate; i++){
-    for(j=1; j<i;j++){
-      for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	/*s2 += param[i][j][nc]*cov[nc];*/
-	s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-	/*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/
-      }
-      ps[i][j]=s2;
-      /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/
-    }
-    for(j=i+1; j<=nlstate+ndeath;j++){
-      for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-	/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
+      for(j=1; j<i;j++){
+	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+	  /*s2 += param[i][j][nc]*cov[nc];*/
+	  s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+/* 	 printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
+	}
+	ps[i][j]=s2;
+/* 	printf("s1=%.17e, s2=%.17e\n",s1,s2); */
+      }
+      for(j=i+1; j<=nlstate+ndeath;j++){
+	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+	  s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+/* 	  printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
+	}
+	ps[i][j]=s2;
       }
-      ps[i][j]=s2;
     }
-  }
     /*ps[3][2]=1;*/
-
-  for(i=1; i<= nlstate; i++){
-     s1=0;
-    for(j=1; j<i; j++)
-      s1+=exp(ps[i][j]);
-    for(j=i+1; j<=nlstate+ndeath; j++)
-      s1+=exp(ps[i][j]);
-    ps[i][i]=1./(s1+1.);
-    for(j=1; j<i; j++)
-      ps[i][j]= exp(ps[i][j])*ps[i][i];
-    for(j=i+1; j<=nlstate+ndeath; j++)
-      ps[i][j]= exp(ps[i][j])*ps[i][i];
-    /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
-  } /* end i */
-
-  for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
-    for(jj=1; jj<= nlstate+ndeath; jj++){
-      ps[ii][jj]=0;
-      ps[ii][ii]=1;
+    
+    for(i=1; i<= nlstate; i++){
+      s1=0;
+      for(j=1; j<i; j++)
+	s1+=exp(ps[i][j]);
+      for(j=i+1; j<=nlstate+ndeath; j++)
+	s1+=exp(ps[i][j]);
+      ps[i][i]=1./(s1+1.);
+      for(j=1; j<i; j++)
+	ps[i][j]= exp(ps[i][j])*ps[i][i];
+      for(j=i+1; j<=nlstate+ndeath; j++)
+	ps[i][j]= exp(ps[i][j])*ps[i][i];
+      /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+    } /* end i */
+    
+    for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+      for(jj=1; jj<= nlstate+ndeath; jj++){
+	ps[ii][jj]=0;
+	ps[ii][ii]=1;
+      }
     }
-  }
-
+    
 
-  /*   for(ii=1; ii<= nlstate+ndeath; ii++){
-    for(jj=1; jj<= nlstate+ndeath; jj++){
-     printf("%lf ",ps[ii][jj]);
-   }
-    printf("\n ");
-    }
-    printf("\n ");printf("%lf ",cov[2]);*/
-/*
-  for(i=1; i<= npar; i++) printf("%f ",x[i]);
-  goto end;*/
+/*        for(ii=1; ii<= nlstate+ndeath; ii++){ */
+/* 	 for(jj=1; jj<= nlstate+ndeath; jj++){ */
+/* 	   printf("ddd %lf ",ps[ii][jj]); */
+/* 	 } */
+/* 	 printf("\n "); */
+/*        } */
+/*        printf("\n ");printf("%lf ",cov[2]); */
+       /*
+      for(i=1; i<= npar; i++) printf("%f ",x[i]);
+      goto end;*/
     return ps;
 }
 
@@ -1224,14 +1238,14 @@ double func( double *x)
 	} /* end mult */
       
 	/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
-	/* But now since version 0.9 we anticipate for bias and large stepm.
+	/* But now since version 0.9 we anticipate for bias at large stepm.
 	 * If stepm is larger than one month (smallest stepm) and if the exact delay 
 	 * (in months) between two waves is not a multiple of stepm, we rounded to 
 	 * the nearest (and in case of equal distance, to the lowest) interval but now
 	 * we keep into memory the bias bh[mi][i] and also the previous matrix product
-	 * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+	 * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
 	 * probability in order to take into account the bias as a fraction of the way
-	 * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+	 * from savm to out if bh is negative or even beyond if bh is positive. bh varies
 	 * -stepm/2 to stepm/2 .
 	 * For stepm=1 the results are the same as for previous versions of Imach.
 	 * For stepm > 1 the results are less biased than in previous versions. 
@@ -1239,15 +1253,15 @@ double func( double *x)
 	s1=s[mw[mi][i]][i];
 	s2=s[mw[mi+1][i]][i];
 	bbh=(double)bh[mi][i]/(double)stepm; 
-	/* bias is positive if real duration
+	/* bias bh is positive if real duration
 	 * is higher than the multiple of stepm and negative otherwise.
 	 */
 	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
 	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 . 
+	     step unit time, which is also equal to probability to die before dh 
+	     minus 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
@@ -4003,7 +4017,7 @@ int main(int argc, char *argv[])
   int *indx;
   char line[MAXLINE], linepar[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
-  char pathr[MAXLINE]; 
+  char pathr[MAXLINE], pathimach[MAXLINE]; 
   int firstobs=1, lastobs=10;
   int sdeb, sfin; /* Status at beginning and end */
   int c,  h , cpt,l;
@@ -4089,8 +4103,10 @@ int main(int argc, char *argv[])
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
   /* cutv(path,optionfile,pathtot,'\\');*/
 
+  split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
+ /*   strcpy(pathimach,argv[0]); */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);
-  printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+  printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);
   strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);
@@ -4115,12 +4131,12 @@ int main(int argc, char *argv[])
   }
   fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s\n%s",version,fullversion);
-  fprintf(ficlog,"\nEnter the parameter file name: ");
-  fprintf(ficlog,"pathtot=%s\n\
+  fprintf(ficlog,"\nEnter the parameter file name: \n");
+  fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\
  optionfile=%s\n\
  optionfilext=%s\n\
- optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
+ optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -5411,7 +5427,10 @@ Interval (in months) between two waves:
   /*------ End -----------*/
 
   chdir(path);
-  strcpy(plotcmd,GNUPLOTPROGRAM);
+  strcpy(plotcmd,"\"");
+  strcat(plotcmd,pathimach);
+  strcat(plotcmd,GNUPLOTPROGRAM);
+  strcat(plotcmd,"\"");
   strcat(plotcmd," ");
   strcat(plotcmd,optionfilegnuplot);
   printf("Starting graphs with: %s",plotcmd);fflush(stdout);
@@ -5424,7 +5443,10 @@ Interval (in months) between two waves:
     printf("\nType e to edit output files, g to graph again and q for exiting: ");
     scanf("%s",z);
 /*     if (z[0] == 'c') system("./imach"); */
-    if (z[0] == 'e') system(optionfilehtm);
+    if (z[0] == 'e') {
+      printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
+      system(optionfilehtm);
+    }
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }