]> henry.ined.fr Git - .git/commitdiff
*** empty log message ***
authorN. Brouard <brouard@ined.fr>
Sat, 5 Jun 2004 08:57:40 +0000 (08:57 +0000)
committerN. Brouard <brouard@ined.fr>
Sat, 5 Jun 2004 08:57:40 +0000 (08:57 +0000)
src/ChangeLog
src/imach.c

index 26c36d7ae678869f0267294d69be5e086f9fec6c..72f8071749de45446fe6159f771686c323eae172 100644 (file)
@@ -1,3 +1,18 @@
+2004-05-20  Brouard Nicolas  <brouard@localhost>\r
+\r
+       * imach.c (Repository): \r
+       Agnes added a direct estimation of mortality (without the need of\r
+       computing period prevalence and differential mortality). Thus here\r
+       is version 0.97a which has been distributed to some people at\r
+       REVES 16 in Brugge using an Inno setup.exe for PCs. Estimates of\r
+       mortality using covariates is not done today. Estimating direct\r
+       mortality is a very different process because it doesn't need\r
+       interpolation because it is easy to get the lx from the force of\r
+       the mortality mux in the simplest case as for a Gompertz (log mux\r
+       = a + b*x . But we have been able to incorporate the new code\r
+       within former imach program (0.96d) without deteriorating too much\r
+       the understanding of the program. \r
+\r
 2003-06-25    <brouard@BROUARD>\r
 \r
        * imach.c (Module): On windows (cygwin) function asctime_r doesn't\r
index 303c80b4cfd3481c14864e636d95e33f86390c2f..cb7b1dd7a427fc28126cd8b9b774f398ecaea7ae 100644 (file)
@@ -1,6 +1,28 @@
 /* $Id$
   $State$
   $Log$
+  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
+  state at each age, but using a Gompertz model: log u =a + b*age .
+  This is the basic analysis of mortality and should be done before any
+  other analysis, in order to test if the mortality estimated from the
+  cross-longitudinal survey is different from the mortality estimated
+  from other sources like vital statistic data.
+
+  The same imach parameter file can be used but the option for mle should be -3.
+
+  Agnès, who wrote this part of the code, tried to keep most of the
+  former routines in order to include the new code within the former code.
+
+  The output is very simple: only an estimate of the intercept and of
+  the slope with 95% confident intervals.
+
+  Current limitations:
+  A) Even if you enter covariates, i.e. with the
+  model= V1+V2 equation for example, the programm does only estimate a unique global model without covariates.
+  B) There is no computation of Life Expectancy nor Life Table.
+
   Revision 1.97  2004/02/20 13:25:42  lievre
   Version 0.96d. Population forecasting command line is (temporarily)
   suppressed.
 #include <stdlib.h>
 #include <unistd.h>
 
-#include <sys/time.h>
+/* #include <sys/time.h> */
 #include <time.h>
 #include "timeval.h"
 
 /* $Id$ */
 /* $State$ */
 
-char version[]="Imach version 0.70, May 2004, INED-EUROREVES ";
+char version[]="Imach version 0.97b, May 2004, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$"; 
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;
@@ -334,6 +356,9 @@ double ftolhess; /* Tolerance for computing hessian */
 /**************** 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 */
 
@@ -365,12 +390,14 @@ static    int split( char *path, char *dirc, char *name, char *ext, char *finame )
 #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 */
 }
 
@@ -1028,57 +1055,57 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
   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);*/
+      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); */
       }
-      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=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;
 }
 
@@ -3981,7 +4008,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;
@@ -4067,8 +4094,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);
@@ -4093,12 +4122,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);
@@ -5389,7 +5418,10 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
   /*------ 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);
@@ -5402,7 +5434,10 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
     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);
   }