[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.98 and 1.110

version 1.98, 2004/05/16 15:05:56 version 1.110, 2006/01/25 00:51:50
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.110  2006/01/25 00:51:50  brouard
     (Module): Lots of cleaning and bugs added (Gompertz)
   
     Revision 1.109  2006/01/24 19:37:15  brouard
     (Module): Comments (lines starting with a #) are allowed in data.
   
     Revision 1.108  2006/01/19 18:05:42  lievre
     Gnuplot problem appeared...
     To be fixed
   
     Revision 1.107  2006/01/19 16:20:37  brouard
     Test existence of gnuplot in imach path
   
     Revision 1.106  2006/01/19 13:24:36  brouard
     Some cleaning and links added in html output
   
     Revision 1.105  2006/01/05 20:23:19  lievre
     *** empty log message ***
   
     Revision 1.104  2005/09/30 16:11:43  lievre
     (Module): sump fixed, loop imx fixed, and simplifications.
     (Module): If the status is missing at the last wave but we know
     that the person is alive, then we can code his/her status as -2
     (instead of missing=-1 in earlier versions) and his/her
     contributions to the likelihood is 1 - Prob of dying from last
     health status (= 1-p13= p11+p12 in the easiest case of somebody in
     the healthy state at last known wave). Version is 0.98
   
     Revision 1.103  2005/09/30 15:54:49  lievre
     (Module): sump fixed, loop imx fixed, and simplifications.
   
     Revision 1.102  2004/09/15 17:31:30  brouard
     Add the possibility to read data file including tab characters.
   
     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    Revision 1.98  2004/05/16 15:05:56  brouard
   New version 0.97 . First attempt to estimate force of mortality    New version 0.97 . First attempt to estimate force of mortality
   directly from the data i.e. without the need of knowing the health    directly from the data i.e. without the need of knowing the health
Line 194 Line 237
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
 #include <stdlib.h>  #include <stdlib.h>
   #include <string.h>
 #include <unistd.h>  #include <unistd.h>
   
 #include <sys/time.h>  #include <limits.h>
   #include <sys/types.h>
   #include <sys/stat.h>
   #include <errno.h>
   extern int errno;
   
   /* #include <sys/time.h> */
 #include <time.h>  #include <time.h>
 #include "timeval.h"  #include "timeval.h"
   
Line 204 Line 254
 /* #define _(String) gettext (String) */  /* #define _(String) gettext (String) */
   
 #define MAXLINE 256  #define MAXLINE 256
   
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 132  #define FILENAMELENGTH 132
 /*#define DEBUG*/  
 /*#define windows*/  
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
Line 224 Line 274
 #define AGESUP 130  #define AGESUP 130
 #define AGEBASE 40  #define AGEBASE 40
 #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */  #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
 #ifdef unix  #ifdef UNIX
 #define DIRSEPARATOR '/'  #define DIRSEPARATOR '/'
   #define CHARSEPARATOR "/"
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
 #else  #else
 #define DIRSEPARATOR '\\'  #define DIRSEPARATOR '\\'
   #define CHARSEPARATOR "\\"
 #define ODIRSEPARATOR '/'  #define ODIRSEPARATOR '/'
 #endif  #endif
   
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.70, May 2004, INED-EUROREVES ";  char version[]="Imach version 0.98a, January 2006, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */  int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;  int nvar;
Line 249  int popbased=0; Line 301  int popbased=0;
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */  int maxwav; /* Maxim number of waves */
 int jmin, jmax; /* min, max spacing between 2 waves */  int jmin, jmax; /* min, max spacing between 2 waves */
   int ijmin, ijmax; /* Individuals having jmin and jmax */ 
 int gipmx, gsw; /* Global variables on the number of contributions   int gipmx, gsw; /* Global variables on the number of contributions 
                    to the likelihood and the sum of weights (done by funcone)*/                     to the likelihood and the sum of weights (done by funcone)*/
 int mle, weightopt;  int mle, weightopt;
Line 301  long time_value; Line 354  long time_value;
 extern long time();  extern long time();
 char strcurr[80], strfor[80];  char strcurr[80], strfor[80];
   
   char *endptr;
   long lval;
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
 #define FTOL 1.0e-10  #define FTOL 1.0e-10
Line 349  double *weight; Line 405  double *weight;
 int **s; /* Status */  int **s; /* Status */
 double *agedc, **covar, idx;  double *agedc, **covar, idx;
 int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;  int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
   double *lsurv, *lpop, *tpop;
   
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */  double ftol=FTOL; /* Tolerance for computing Max Likelihood */
 double ftolhess; /* Tolerance for computing hessian */  double ftolhess; /* Tolerance for computing hessian */
Line 356  double ftolhess; /* Tolerance for comput Line 413  double ftolhess; /* Tolerance for comput
 /**************** split *************************/  /**************** split *************************/
 static  int split( char *path, char *dirc, char *name, char *ext, char *finame )  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 */    char  *ss;                            /* pointer */
   int   l1, l2;                         /* length counters */    int   l1, l2;                         /* length counters */
   
   l1 = strlen(path );                   /* length of path */    l1 = strlen(path );                   /* length of path */
   if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );    if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   ss= strrchr( path, DIRSEPARATOR );            /* find last / */    ss= strrchr( path, DIRSEPARATOR );            /* find last / */
   if ( ss == NULL ) {                   /* no directory, so use current */    if ( ss == NULL ) {                   /* no directory, so determine current directory */
       strcpy( name, path );               /* we got the fullname name because no directory */
     /*if(strrchr(path, ODIRSEPARATOR )==NULL)      /*if(strrchr(path, ODIRSEPARATOR )==NULL)
       printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/        printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
     /* get current working directory */      /* get current working directory */
Line 370  static int split( char *path, char *dirc Line 431  static int split( char *path, char *dirc
     if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {      if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
       return( GLOCK_ERROR_GETCWD );        return( GLOCK_ERROR_GETCWD );
     }      }
     strcpy( name, path );               /* we've got it */      /* got dirc from getcwd*/
       printf(" DIRC = %s \n",dirc);
   } else {                              /* strip direcotry from path */    } else {                              /* strip direcotry from path */
     ss++;                               /* after this, the filename */      ss++;                               /* after this, the filename */
     l2 = strlen( ss );                  /* length of filename */      l2 = strlen( ss );                  /* length of filename */
Line 378  static int split( char *path, char *dirc Line 440  static int split( char *path, char *dirc
     strcpy( name, ss );         /* save file name */      strcpy( name, ss );         /* save file name */
     strncpy( dirc, path, l1 - l2 );     /* now the directory */      strncpy( dirc, path, l1 - l2 );     /* now the directory */
     dirc[l1-l2] = 0;                    /* add zero */      dirc[l1-l2] = 0;                    /* add zero */
       printf(" DIRC2 = %s \n",dirc);
   }    }
     /* We add a separator at the end of dirc if not exists */
   l1 = strlen( dirc );                  /* length of directory */    l1 = strlen( dirc );                  /* length of directory */
   /*#ifdef windows    if( dirc[l1-1] != DIRSEPARATOR ){
   if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }      dirc[l1] =  DIRSEPARATOR;
 #else      dirc[l1+1] = 0; 
   if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }      printf(" DIRC3 = %s \n",dirc);
 #endif    }
   */  
   ss = strrchr( name, '.' );            /* find last / */    ss = strrchr( name, '.' );            /* find last / */
   ss++;    if (ss >0){
   strcpy(ext,ss);                       /* save extension */      ss++;
   l1= strlen( name);      strcpy(ext,ss);                     /* save extension */
   l2= strlen(ss)+1;      l1= strlen( name);
   strncpy( finame, name, l1-l2);      l2= strlen(ss)+1;
   finame[l1-l2]= 0;      strncpy( finame, name, l1-l2);
       finame[l1-l2]= 0;
     }
   
   return( 0 );                          /* we're done */    return( 0 );                          /* we're done */
 }  }
   
Line 425  int nbocc(char *s, char occ) Line 491  int nbocc(char *s, char occ)
   
 void cutv(char *u,char *v, char*t, char occ)  void cutv(char *u,char *v, char*t, char occ)
 {  {
   /* cuts string t into u and v where u is ended by char occ excluding it    /* cuts string t into u and v where u ends before first occurence of char 'occ' 
      and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)       and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')
      gives u="abcedf" and v="ghi2j" */       gives u="abcedf" and v="ghi2j" */
   int i,lg,j,p=0;    int i,lg,j,p=0;
   i=0;    i=0;
Line 868  void powell(double p[], double **xi, int Line 934  void powell(double p[], double **xi, int
     fprintf(ficrespow,"\n");fflush(ficrespow);      fprintf(ficrespow,"\n");fflush(ficrespow);
     if(*iter <=3){      if(*iter <=3){
       tm = *localtime(&curr_time.tv_sec);        tm = *localtime(&curr_time.tv_sec);
       strcpy(strcurr,asctime(&tmf));        strcpy(strcurr,asctime(&tm));
 /*       asctime_r(&tm,strcurr); */  /*       asctime_r(&tm,strcurr); */
       forecast_time=curr_time;        forecast_time=curr_time; 
       itmp = strlen(strcurr);        itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')        if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
         strcurr[itmp-1]='\0';          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);        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);        fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
Line 884  void powell(double p[], double **xi, int Line 950  void powell(double p[], double **xi, int
         itmp = strlen(strfor);          itmp = strlen(strfor);
         if(strfor[itmp-1]=='\n')          if(strfor[itmp-1]=='\n')
         strfor[itmp-1]='\0';          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);          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 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 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++) {       for (i=1;i<=n;i++) { 
Line 1050  double **pmij(double **ps, double *cov, Line 1116  double **pmij(double **ps, double *cov,
   int i,j,j1, nc, ii, jj;    int i,j,j1, nc, ii, jj;
   
     for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
     for(j=1; j<i;j++){        for(j=1; j<i;j++){
       for (nc=1, s2=0.;nc <=ncovmodel; nc++){          for (nc=1, s2=0.;nc <=ncovmodel; nc++){
         /*s2 += param[i][j][nc]*cov[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];            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);*/  /*       printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
       }          }
       ps[i][j]=s2;          ps[i][j]=s2;
       /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/  /*      printf("s1=%.17e, s2=%.17e\n",s1,s2); */
     }        }
     for(j=i+1; j<=nlstate+ndeath;j++){        for(j=i+1; j<=nlstate+ndeath;j++){
       for (nc=1, s2=0.;nc <=ncovmodel; nc++){          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];            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);*/  /*        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;*/      /*ps[3][2]=1;*/
       
   for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
      s1=0;        s1=0;
     for(j=1; j<i; j++)        for(j=1; j<i; j++)
       s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]);
     for(j=i+1; j<=nlstate+ndeath; j++)        for(j=i+1; j<=nlstate+ndeath; j++)
       s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]);
     ps[i][i]=1./(s1+1.);        ps[i][i]=1./(s1+1.);
     for(j=1; j<i; j++)        for(j=1; j<i; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];          ps[i][j]= exp(ps[i][j])*ps[i][i];
     for(j=i+1; j<=nlstate+ndeath; j++)        for(j=i+1; j<=nlstate+ndeath; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];          ps[i][j]= exp(ps[i][j])*ps[i][i];
     /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */        /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
   } /* end i */      } /* end i */
       
   for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){      for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){        for(jj=1; jj<= nlstate+ndeath; jj++){
       ps[ii][jj]=0;          ps[ii][jj]=0;
       ps[ii][ii]=1;          ps[ii][ii]=1;
         }
     }      }
   }      
   
   
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){  /*        for(ii=1; ii<= nlstate+ndeath; ii++){ */
     for(jj=1; jj<= nlstate+ndeath; jj++){  /*       for(jj=1; jj<= nlstate+ndeath; jj++){ */
      printf("%lf ",ps[ii][jj]);  /*         printf("ddd %lf ",ps[ii][jj]); */
    }  /*       } */
     printf("\n ");  /*       printf("\n "); */
     }  /*        } */
     printf("\n ");printf("%lf ",cov[2]);*/  /*        printf("\n ");printf("%lf ",cov[2]); */
 /*         /*
   for(i=1; i<= npar; i++) printf("%f ",x[i]);        for(i=1; i<= npar; i++) printf("%f ",x[i]);
   goto end;*/        goto end;*/
     return ps;      return ps;
 }  }
   
Line 1224  double func( double *x) Line 1290  double func( double *x)
         } /* end mult */          } /* end mult */
               
         /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */          /*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            * 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            * (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           * 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           * 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           * 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 .           * -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 the same as for previous versions of Imach.
          * For stepm > 1 the results are less biased than in previous versions.            * For stepm > 1 the results are less biased than in previous versions. 
Line 1239  double func( double *x) Line 1305  double func( double *x)
         s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;           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.           * 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]));*/          /* 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){           if( s2 > nlstate){ 
           /* i.e. if s2 is a death state and if the date of death is known then the contribution            /* i.e. if s2 is a death state and if the date of death is known 
              to the likelihood is the probability to die between last step unit time and current                then the contribution to the likelihood is the probability to 
              step unit time, which is also the differences between probability to die before dh                die between last step unit time and current  step unit time, 
              and probability to die before dh-stepm .                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               In version up to 0.92 likelihood was computed
         as if date of death was unknown. Death was treated as any other          as if date of death was unknown. Death was treated as any other
         health state: the date of the interview describes the actual state          health state: the date of the interview describes the actual state
Line 1267  double func( double *x) Line 1334  double func( double *x)
         lower mortality.          lower mortality.
           */            */
           lli=log(out[s1][s2] - savm[s1][s2]);            lli=log(out[s1][s2] - savm[s1][s2]);
         }else{  
   
           } else if  (s2==-2) {
             for (j=1,survp=0. ; j<=nlstate; j++) 
               survp += out[s1][j];
             lli= survp;
           }
           
           else if  (s2==-4) {
             for (j=3,survp=0. ; j<=nlstate; j++) 
               survp += out[s1][j];
             lli= survp;
           }
           
           else if  (s2==-5) {
             for (j=1,survp=0. ; j<=2; j++) 
               survp += out[s1][j];
             lli= survp;
           }
   
   
           else{
           lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */            lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
           /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */            /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
         }           } 
Line 1831  void lubksb(double **a, int n, int *indx Line 1919  void lubksb(double **a, int n, int *indx
 }   } 
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint)  void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
 {  /* Some frequencies */  {  /* Some frequencies */
       
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;    int i, m, jk, k1,i1, j1, bool, z1,z2,j;
Line 1851  void  freqsummary(char fileres[], int ia Line 1939  void  freqsummary(char fileres[], int ia
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);      exit(0);
   }    }
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);    freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
   j1=0;    j1=0;
       
   j=cptcoveff;    j=cptcoveff;
Line 1864  void  freqsummary(char fileres[], int ia Line 1952  void  freqsummary(char fileres[], int ia
       j1++;        j1++;
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
         scanf("%d", i);*/          scanf("%d", i);*/
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-5; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-5; jk<=nlstate+ndeath; jk++)  
           for(m=iagemin; m <= iagemax+3; m++)            for(m=iagemin; m <= iagemax+3; m++)
             freq[i][jk][m]=0;              freq[i][jk][m]=0;
   
Line 1904  void  freqsummary(char fileres[], int ia Line 1992  void  freqsummary(char fileres[], int ia
       }        }
                 
       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
   fprintf(ficresp, "#Local time at start: %s", strstart);
       if  (cptcovn>0) {        if  (cptcovn>0) {
         fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable "); 
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 1990  void  freqsummary(char fileres[], int ia Line 2078  void  freqsummary(char fileres[], int ia
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);    free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   free_matrix(prop,1,nlstate,iagemin, iagemax+3);    free_matrix(prop,1,nlstate,iagemin, iagemax+3);
   /* End of Freq */    /* End of Freq */
Line 2099  void  concatwav(int wav[], int **dh, int Line 2187  void  concatwav(int wav[], int **dh, int
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
     while(s[m][i] <= nlstate){      while(s[m][i] <= nlstate){
       if(s[m][i]>=1)        if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)
         mw[++mi][i]=m;          mw[++mi][i]=m;
       if(m >=lastpass)        if(m >=lastpass)
         break;          break;
Line 2139  void  concatwav(int wav[], int **dh, int Line 2227  void  concatwav(int wav[], int **dh, int
               nberr++;                nberr++;
               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               j=1; /* Temporary Dangerous patch */                j=1; /* Temporary Dangerous patch */
               printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);                printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
               fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);                fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
             }              }
             k=k+1;              k=k+1;
             if (j >= jmax) jmax=j;              if (j >= jmax){
             if (j <= jmin) jmin=j;                jmax=j;
                 ijmax=i;
               }
               if (j <= jmin){
                 jmin=j;
                 ijmin=i;
               }
             sum=sum+j;              sum=sum+j;
             /*if (j<0) printf("j=%d num=%d \n",j,i);*/              /*if (j<0) printf("j=%d num=%d \n",j,i);*/
             /*    printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/              /*    printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
Line 2153  void  concatwav(int wav[], int **dh, int Line 2247  void  concatwav(int wav[], int **dh, int
         }          }
         else{          else{
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));            j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
           /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/  /*        if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
   
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) {
           else if (j <= jmin)jmin=j;              jmax=j;
               ijmax=i;
             }
             else if (j <= jmin){
               jmin=j;
               ijmin=i;
             }
           /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */            /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/            /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
           if(j<0){            if(j<0){
Line 2199  void  concatwav(int wav[], int **dh, int Line 2300  void  concatwav(int wav[], int **dh, int
     } /* end wave */      } /* end wave */
   }    }
   jmean=sum/k;    jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);    printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
   fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);    fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
  }   }
   
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
Line 2246  void tricode(int *Tvar, int **nbcode, in Line 2347  void tricode(int *Tvar, int **nbcode, in
  for (k=0; k< maxncov; k++) Ndum[k]=0;   for (k=0; k< maxncov; k++) Ndum[k]=0;
   
  for (i=1; i<=ncovmodel-2; i++) {    for (i=1; i<=ncovmodel-2; i++) { 
    /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/     /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
    ij=Tvar[i];     ij=Tvar[i];
    Ndum[ij]++;     Ndum[ij]++;
  }   }
Line 2264  void tricode(int *Tvar, int **nbcode, in Line 2365  void tricode(int *Tvar, int **nbcode, in
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov,char strstart[] )
   
 {  {
   /* Health expectancies */    /* Health expectancies */
Line 2282  void evsij(char fileres[], double ***eij Line 2383  void evsij(char fileres[], double ***eij
   dnewm=matrix(1,nlstate*nlstate,1,npar);    dnewm=matrix(1,nlstate*nlstate,1,npar);
   doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);    doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);
       
     fprintf(ficreseij,"# Local time at start: %s", strstart);
   fprintf(ficreseij,"# Health expectancies\n");    fprintf(ficreseij,"# Health expectancies\n");
   fprintf(ficreseij,"# Age");    fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
Line 2436  void evsij(char fileres[], double ***eij Line 2538  void evsij(char fileres[], double ***eij
 }  }
   
 /************ Variance ******************/  /************ Variance ******************/
 void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav)  void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
 {  {
   /* Variance of health expectancies */    /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
Line 2487  void varevsij(char optionfilefiname[], d Line 2589  void varevsij(char optionfilefiname[], d
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
   }    }
   printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);    printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
    
   fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
     fprintf(ficresprobmorprev, "#Local time at start: %s", strstart);
   fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
   fprintf(ficresprobmorprev,"# Age cov=%-d",ij);    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
   for(j=nlstate+1; j<=(nlstate+ndeath);j++){    for(j=nlstate+1; j<=(nlstate+ndeath);j++){
Line 2497  void varevsij(char optionfilefiname[], d Line 2601  void varevsij(char optionfilefiname[], d
   }      }  
   fprintf(ficresprobmorprev,"\n");    fprintf(ficresprobmorprev,"\n");
   fprintf(ficgp,"\n# Routine varevsij");    fprintf(ficgp,"\n# Routine varevsij");
     /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
   fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
 /*   } */  /*   } */
   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
    fprintf(ficresvij, "#Local time at start: %s", strstart);
   fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");
   fprintf(ficresvij,"# Age");    fprintf(ficresvij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
Line 2734  void varevsij(char optionfilefiname[], d Line 2839  void varevsij(char optionfilefiname[], d
 }  /* end varevsij */  }  /* end varevsij */
   
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
 void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)  void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[])
 {  {
   /* Variance of prevalence limit */    /* Variance of prevalence limit */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
Line 2747  void varprevlim(char fileres[], double * Line 2852  void varprevlim(char fileres[], double *
   double **gradg, **trgradg;    double **gradg, **trgradg;
   double age,agelim;    double age,agelim;
   int theta;    int theta;
        fprintf(ficresvpl, "#Local time at start: %s", strstart); 
   fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");    fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
   fprintf(ficresvpl,"# Age");    fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
Line 2817  void varprevlim(char fileres[], double * Line 2922  void varprevlim(char fileres[], double *
 }  }
   
 /************ Variance of one-step probabilities  ******************/  /************ Variance of one-step probabilities  ******************/
 void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)  void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[])
 {  {
   int i, j=0,  i1, k1, l1, t, tj;    int i, j=0,  i1, k1, l1, t, tj;
   int k2, l2, j1,  z1;    int k2, l2, j1,  z1;
Line 2862  void varprob(char optionfilefiname[], do Line 2967  void varprob(char optionfilefiname[], do
   fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);    fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
   printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);    printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
   fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);    fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
       fprintf(ficresprob, "#Local time at start: %s", strstart);
   fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");    fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");    fprintf(ficresprob,"# Age");
     fprintf(ficresprobcov, "#Local time at start: %s", strstart);
   fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");    fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
   fprintf(ficresprobcov,"# Age");    fprintf(ficresprobcov,"# Age");
     fprintf(ficresprobcor, "#Local time at start: %s", strstart);
   fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");    fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
   fprintf(ficresprobcov,"# Age");    fprintf(ficresprobcov,"# Age");
   
Line 3137  void printinghtml(char fileres[], char t Line 3244  void printinghtml(char fileres[], char t
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   
    fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
      <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
   </ul>");
      fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",   - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
Line 3188  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3298  fprintf(fichtm," \n<ul><li><b>Graphs</b>
   
   
  fprintf(fichtm,"\   fprintf(fichtm,"\
 \n<br><li><h4> Result files (second order: variances)</h4>\n\  \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);
   
  fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",   fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
Line 3902  double gompertz(double x[]) Line 4012  double gompertz(double x[])
 {   { 
   double A,B,L=0.0,sump=0.,num=0.;    double A,B,L=0.0,sump=0.,num=0.;
   int i,n=0; /* n is the size of the sample */    int i,n=0; /* n is the size of the sample */
   
   for (i=0;i<=imx-1 ; i++) {    for (i=0;i<=imx-1 ; i++) {
     sump=sump+weight[i];      sump=sump+weight[i];
     sump=sump+1;      /*    sump=sump+1;*/
     num=num+1;      num=num+1;
   }    }
     
     
   /* for (i=1; i<=imx; i++)     /* for (i=0; i<=imx; i++) 
      if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/       if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
   
   for (i=0;i<=imx-1 ; i++)    for (i=1;i<=imx ; i++)
     {      {
       if (cens[i]==1 & wav[i]>1)        if (cens[i] == 1 && wav[i]>1)
         A=-x[1]/(x[2])*          A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
           (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)));  
               
       if (cens[i]==0 & wav[i]>1)        if (cens[i] == 0 && wav[i]>1)
         A=-x[1]/(x[2])*          A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
              (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)))               +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);  
           +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM);        
               
       if (wav[i]>1 & agecens[i]>15) {        /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
         if (wav[i] > 1 ) { /* ??? */
         L=L+A*weight[i];          L=L+A*weight[i];
         /*      printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/          /*      printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
       }        }
Line 3937  double gompertz(double x[]) Line 4047  double gompertz(double x[])
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \  void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
                   int imx,  double p[],double **matcov){                    int imx,  double p[],double **matcov,double agemortsup){
   int i;    int i,k;
   
   fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");    fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");
   fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);    fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);
Line 3946  void printinghtmlmort(char fileres[], ch Line 4056  void printinghtmlmort(char fileres[], ch
     fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));      fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
   fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");    fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");
   fprintf(fichtm,"</ul>");    fprintf(fichtm,"</ul>");
   
   fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>");
   
    fprintf(fichtm,"\nAge   l<inf>x</inf>     q<inf>x</inf> d(x,x+1)    L<inf>x</inf>     T<inf>x</inf>     e<infx</inf><br>");
   
    for (k=agegomp;k<(agemortsup-2);k++) 
      fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf<br>\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
   
    
   fflush(fichtm);    fflush(fichtm);
 }  }
   
Line 3983  int main(int argc, char *argv[]) Line 4102  int main(int argc, char *argv[])
 {  {
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);    int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
     int linei, month, year,iout;
   int jj, ll, li, lj, lk, imk;    int jj, ll, li, lj, lk, imk;
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
   int itimes;    int itimes;
Line 3991  int main(int argc, char *argv[]) Line 4111  int main(int argc, char *argv[])
   char ca[32], cb[32], cc[32];    char ca[32], cb[32], cc[32];
   /*  FILE *fichtm; *//* Html File */    /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */    /* FILE *ficgp;*/ /*Gnuplot File */
     struct stat info;
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;    double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
   
Line 4003  int main(int argc, char *argv[]) Line 4124  int main(int argc, char *argv[])
   int *indx;    int *indx;
   char line[MAXLINE], linepar[MAXLINE];    char line[MAXLINE], linepar[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[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 firstobs=1, lastobs=10;
   int sdeb, sfin; /* Status at beginning and end */    int sdeb, sfin; /* Status at beginning and end */
   int c,  h , cpt,l;    int c,  h , cpt,l;
Line 4013  int main(int argc, char *argv[]) Line 4134  int main(int argc, char *argv[])
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */    int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;    int mobilav=0,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
     int agemortsup;
     float  sumlpop=0.;
   double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;    double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
   double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;    double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
   
Line 4089  int main(int argc, char *argv[]) Line 4212  int main(int argc, char *argv[])
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/      printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
   /* cutv(path,optionfile,pathtot,'\\');*/    /* cutv(path,optionfile,pathtot,'\\');*/
   
     /* Split argv[0], imach program to get pathimach */
     printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]);
     split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
     printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
    /*   strcpy(pathimach,argv[0]); */
     /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);    chdir(path);
   strcpy(command,"mkdir ");    strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);    strcat(command,optionfilefiname);
Line 4115  int main(int argc, char *argv[]) Line 4244  int main(int argc, char *argv[])
   }    }
   fprintf(ficlog,"Log filename:%s\n",filelog);    fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s\n%s",version,fullversion);    fprintf(ficlog,"\n%s\n%s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: ");    fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathtot=%s\n\    fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\   path=%s \n\
  optionfile=%s\n\   optionfile=%s\n\
  optionfilext=%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);    printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);    fprintf(ficlog,"Local time (at start): %s",strstart);
Line 4257  int main(int argc, char *argv[]) Line 4386  int main(int argc, char *argv[])
     }        }  
     fflush(ficlog);      fflush(ficlog);
   
   
     p=param[1][1];      p=param[1][1];
           
     /* Reads comments: lines beginning with '#' */      /* Reads comments: lines beginning with '#' */
Line 4381  int main(int argc, char *argv[]) Line 4509  int main(int argc, char *argv[])
   ncodemax=ivector(1,8);    ncodemax=ivector(1,8);
   
   i=1;    i=1;
   while (fgets(line, MAXLINE, fic) != NULL)    {    linei=0;
     if ((i >= firstobs) && (i <=lastobs)) {    while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs)))    {
               linei=linei+1;
       for(j=strlen(line); j>=0;j--){  /* Untabifies line */
           if(line[j] == '\t')
             line[j] = ' ';
         }
         for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
           ;
         };
         line[j+1]=0;  /* Trims blanks at end of line */
         if(line[0]=='#'){
           fprintf(ficlog,"Comment line\n%s\n",line);
           printf("Comment line\n%s\n",line);
           continue;
         }
       for (j=maxwav;j>=1;j--){        for (j=maxwav;j>=1;j--){
         cutv(stra, strb,line,' '); s[j][i]=atoi(strb);           cutv(stra, strb,line,' '); 
           errno=0;
           lval=strtol(strb,&endptr,10); 
           /*      if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
           if( strb[0]=='\0' || (*endptr != '\0')){
             printf("Error reading data around '%d' at line number %d %s for individual %d\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n",lval, i,line,linei,j,maxwav);
             exit(1);
           }
           s[j][i]=lval;
   
         strcpy(line,stra);          strcpy(line,stra);
         cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);          cutv(stra, strb,line,' ');
         cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);          if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
       }          }
           else  if(iout=sscanf(strb,".") != 0){
             month=99;
             year=9999;
           }else{
             printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
             exit(1);
           }
           anint[j][i]= (double) year; 
           mint[j][i]= (double)month; 
           strcpy(line,stra);
         } /* ENd Waves */
                   
       cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);        cutv(stra, strb,line,' '); 
       cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);        if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
         }
       cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);        else  if(iout=sscanf(strb,".") != 0){
       cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);          month=99;
           year=9999;
         }else{
           printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
           exit(1);
         }
         andc[i]=(double) year; 
         moisdc[i]=(double) month; 
         strcpy(line,stra);
   
         cutv(stra, strb,line,' '); 
         if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
         }
         else  if(iout=sscanf(strb,".") != 0){
           month=99;
           year=9999;
         }else{
           printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
           exit(1);
         }
         annais[i]=(double)(year);
         moisnais[i]=(double)(month); 
         strcpy(line,stra);
   
         cutv(stra, strb,line,' '); 
         errno=0;
         lval=strtol(strb,&endptr,10); 
         if( strb[0]=='\0' || (*endptr != '\0')){
           printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight.  Exiting.\n",lval, i,line,linei);
           exit(1);
         }
         weight[i]=(double)(lval); 
         strcpy(line,stra);
   
       cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);  
       for (j=ncovcol;j>=1;j--){        for (j=ncovcol;j>=1;j--){
         cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);          cutv(stra, strb,line,' '); 
           errno=0;
           lval=strtol(strb,&endptr,10); 
           if( strb[0]=='\0' || (*endptr != '\0')){
             printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a covar (meaning 0 for the reference or 1).  Exiting.\n",lval, i,line,linei);
             exit(1);
           }
           if(lval <-1 || lval >1){
             printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self).  Exiting.\n",lval,i,line,linei,j);
             exit(1);
           }
           covar[j][i]=(double)(lval);
           strcpy(line,stra);
       }         } 
       lstra=strlen(stra);        lstra=strlen(stra);
   
       if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */        if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
         stratrunc = &(stra[lstra-9]);          stratrunc = &(stra[lstra-9]);
         num[i]=atol(stratrunc);          num[i]=atol(stratrunc);
       }        }
       else        else
         num[i]=atol(stra);          num[i]=atol(stra);
           
       /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){        /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
         printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/          printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
   
       i=i+1;        i=i+1;
     }    } /* End loop reading  data */
   }  
   /* printf("ii=%d", ij);    /* printf("ii=%d", ij);
      scanf("%d",i);*/       scanf("%d",i);*/
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
Line 4428  int main(int argc, char *argv[]) Line 4631  int main(int argc, char *argv[])
      if (s[4][i]==9)  s[4][i]=-1;        if (s[4][i]==9)  s[4][i]=-1; 
      printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}*/       printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}*/
       
  for (i=1; i<=imx; i++)    /* for (i=1; i<=imx; i++) */
     
    /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;     /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;
      else weight[i]=1;*/       else weight[i]=1;*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameters from char model */
   Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */    Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
   Tprod=ivector(1,15);     Tprod=ivector(1,15); 
   Tvaraff=ivector(1,15);     Tvaraff=ivector(1,15); 
Line 4546  int main(int argc, char *argv[]) Line 4749  int main(int argc, char *argv[])
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);      agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
     for(m=firstpass; (m<= lastpass); m++){      for(m=firstpass; (m<= lastpass); m++){
       if(s[m][i] >0){        if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0)            if(agedc[i]>0)
             if((int)moisdc[i]!=99 && (int)andc[i]!=9999)              if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
Line 4562  int main(int argc, char *argv[]) Line 4765  int main(int argc, char *argv[])
             }              }
         }          }
         else if(s[m][i] !=9){ /* Standard case, age in fractional          else if(s[m][i] !=9){ /* Standard case, age in fractional
                                  years but with the precision of a                                   years but with the precision of a month */
                                  month */  
           agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);            agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
           if((int)mint[m][i]==99 || (int)anint[m][i]==9999)            if((int)mint[m][i]==99 || (int)anint[m][i]==9999)
             agev[m][i]=1;              agev[m][i]=1;
Line 4720  Title=%s <br>Datafile=%s Firstpass=%d La Line 4922  Title=%s <br>Datafile=%s Firstpass=%d La
       
   /* Calculates basic frequencies. Computes observed prevalence at single age    /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
   fprintf(fichtm,"<br>Total number of observations=%d <br>\n\    fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
Line 4739  Interval (in months) between two waves: Line 4941  Interval (in months) between two waves:
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
   
   globpr=0; /* To get the number ipmx of contributions and the sum of weights*/    globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
   
   if (mle==-3){    if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);      ximort=matrix(1,NDIM,1,NDIM);
     cens=ivector(1,n);      cens=ivector(1,n);
Line 4748  Interval (in months) between two waves: Line 4951  Interval (in months) between two waves:
     
     for (i=1; i<=imx; i++){      for (i=1; i<=imx; i++){
       dcwave[i]=-1;        dcwave[i]=-1;
       for (j=1; j<=lastpass; j++)        for (m=firstpass; m<=lastpass; m++)
         if (s[j][i]>nlstate) {          if (s[m][i]>nlstate) {
           dcwave[i]=j;            dcwave[i]=m;
           /*    printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/            /*    printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
           break;            break;
         }          }
Line 4759  Interval (in months) between two waves: Line 4962  Interval (in months) between two waves:
     for (i=1; i<=imx; i++) {      for (i=1; i<=imx; i++) {
       if (wav[i]>0){        if (wav[i]>0){
         ageexmed[i]=agev[mw[1][i]][i];          ageexmed[i]=agev[mw[1][i]][i];
         j=wav[i];agecens[i]=1.;           j=wav[i];
         if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];          agecens[i]=1.; 
         cens[i]=1;  
                   if (ageexmed[i]> 1 && wav[i] > 0){
         if (ageexmed[i]<1) cens[i]=-1;            agecens[i]=agev[mw[j][i]][i];
         if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;            cens[i]= 1;
           }else if (ageexmed[i]< 1) 
             cens[i]= -1;
           if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
             cens[i]=0 ;
       }        }
       else cens[i]=-1;        else cens[i]=-1;
     }      }
Line 4773  Interval (in months) between two waves: Line 4980  Interval (in months) between two waves:
       for (j=1;j<=NDIM;j++)        for (j=1;j<=NDIM;j++)
         ximort[i][j]=(i == j ? 1.0 : 0.0);          ximort[i][j]=(i == j ? 1.0 : 0.0);
     }      }
       
     p[1]=0.1; p[2]=0.1;      p[1]=0.1; p[NDIM]=0.1;
     /*printf("%lf %lf", p[1], p[2]);*/      /*printf("%lf %lf", p[1], p[2]);*/
           
           
   printf("Powell\n");  fprintf(ficlog,"Powell\n");      printf("Powell\n");  fprintf(ficlog,"Powell\n");
   strcpy(filerespow,"pow-mort");       strcpy(filerespow,"pow-mort"); 
   strcat(filerespow,fileres);      strcat(filerespow,fileres);
   if((ficrespow=fopen(filerespow,"w"))==NULL) {      if((ficrespow=fopen(filerespow,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", filerespow);        printf("Problem with resultfile: %s\n", filerespow);
     fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);        fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
   }      }
   fprintf(ficrespow,"# Powell\n# iter -2*LL");      fprintf(ficrespow,"# Powell\n# iter -2*LL");
   /*  for (i=1;i<=nlstate;i++)      /*  for (i=1;i<=nlstate;i++)
     for(j=1;j<=nlstate+ndeath;j++)          for(j=1;j<=nlstate+ndeath;j++)
       if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);          if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
   */      */
   fprintf(ficrespow,"\n");      fprintf(ficrespow,"\n");
       
     powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);      powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
     fclose(ficrespow);      fclose(ficrespow);
           
     hesscov(matcov, p, NDIM,delti, 1e-4, gompertz);       hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); 
   
     for(i=1; i <=NDIM; i++)      for(i=1; i <=NDIM; i++)
       for(j=i+1;j<=NDIM;j++)        for(j=i+1;j<=NDIM;j++)
Line 4812  Interval (in months) between two waves: Line 5019  Interval (in months) between two waves:
     printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);      printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
     for (i=1;i<=NDIM;i++)       for (i=1;i<=NDIM;i++) 
       printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));        printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
   
       lsurv=vector(1,AGESUP);
       lpop=vector(1,AGESUP);
       tpop=vector(1,AGESUP);
       lsurv[agegomp]=100000;
       
       for (k=agegomp;k<=AGESUP;k++) {
         agemortsup=k;
         if (p[1]*exp(p[2]*(k-agegomp))>1) break;
       }
       
       for (k=agegomp;k<agemortsup;k++)
         lsurv[k+1]=lsurv[k]-lsurv[k]*(p[1]*exp(p[2]*(k-agegomp)));
       
       for (k=agegomp;k<agemortsup;k++){
         lpop[k]=(lsurv[k]+lsurv[k+1])/2.;
         sumlpop=sumlpop+lpop[k];
       }
       
       tpop[agegomp]=sumlpop;
       for (k=agegomp;k<(agemortsup-3);k++){
         /*  tpop[k+1]=2;*/
         tpop[k+1]=tpop[k]-lpop[k];
       }
       
       
       printf("\nAge   lx     qx    dx    Lx     Tx     e(x)\n");
       for (k=agegomp;k<(agemortsup-2);k++) 
         printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
       
       
     replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
     printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
           
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \      printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
                      stepm, weightopt,\                       stepm, weightopt,\
                      model,imx,p,matcov);                       model,imx,p,matcov,agemortsup);
       
       free_vector(lsurv,1,AGESUP);
       free_vector(lpop,1,AGESUP);
       free_vector(tpop,1,AGESUP);
   } /* Endof if mle==-3 */    } /* Endof if mle==-3 */
     
   else{ /* For mle >=1 */    else{ /* For mle >=1 */
       
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */      likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
Line 5077  Interval (in months) between two waves: Line 5319  Interval (in months) between two waves:
     }      }
     printf("Computing stable prevalence: result on file '%s' \n", filerespl);      printf("Computing stable prevalence: result on file '%s' \n", filerespl);
     fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);      fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
       fprintf(ficrespl, "#Local time at start: %s", strstart);
     fprintf(ficrespl,"#Stable prevalence \n");      fprintf(ficrespl,"#Stable prevalence \n");
     fprintf(ficrespl,"#Age ");      fprintf(ficrespl,"#Age ");
     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);      for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
Line 5137  Interval (in months) between two waves: Line 5380  Interval (in months) between two waves:
     hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
     /* hstepm=1;   aff par mois*/      /* hstepm=1;   aff par mois*/
       fprintf(ficrespij, "#Local time at start: %s", strstart);
     fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");      fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
     for(cptcov=1,k=0;cptcov<=i1;cptcov++){      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
Line 5174  Interval (in months) between two waves: Line 5417  Interval (in months) between two waves:
       }        }
     }      }
   
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
   
     fclose(ficrespij);      fclose(ficrespij);
   
Line 5263  Interval (in months) between two waves: Line 5506  Interval (in months) between two waves:
   
         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);            evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);  
     
         vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);          varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart);
         if(popbased==1){          if(popbased==1){
           varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav);            varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
         }          }
   
            fprintf(ficrest, "#Local time at start: %s", strstart);
         fprintf(ficrest,"#Total LEs with variances: e.. (std) ");          fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
Line 5345  Interval (in months) between two waves: Line 5588  Interval (in months) between two waves:
               
         varpl=matrix(1,nlstate,(int) bage, (int) fage);          varpl=matrix(1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);          varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
         free_matrix(varpl,1,nlstate,(int) bage, (int)fage);          free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
       }        }
     }      }
Line 5394  Interval (in months) between two waves: Line 5637  Interval (in months) between two waves:
   tm = *localtime(&end_time.tv_sec);    tm = *localtime(&end_time.tv_sec);
   tmg = *gmtime(&end_time.tv_sec);    tmg = *gmtime(&end_time.tv_sec);
   strcpy(strtend,asctime(&tm));    strcpy(strtend,asctime(&tm));
   printf("Local time at start %s\nLocaltime at end   %s",strstart, strtend);     printf("Local time at start %s\nLocal time at end   %s",strstart, strtend); 
   fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend);     fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend); 
   printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));    printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
   
Line 5411  Interval (in months) between two waves: Line 5654  Interval (in months) between two waves:
   /*------ End -----------*/    /*------ End -----------*/
   
   chdir(path);    chdir(path);
   strcpy(plotcmd,GNUPLOTPROGRAM);  #ifndef UNIX
     /*  strcpy(plotcmd,"\""); */
   #endif
     strcpy(plotcmd,pathimach);
     /*strcat(plotcmd,CHARSEPARATOR);*/
     strcat(plotcmd,GNUPLOTPROGRAM);
   #ifndef UNIX
     strcat(plotcmd,".exe");
     /*  strcat(plotcmd,"\"");*/
   #endif
     if(stat(plotcmd,&info)){
       printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
     }
   
   #ifndef UNIX
     strcpy(plotcmd,"\"");
   #endif
     strcat(plotcmd,pathimach);
     strcat(plotcmd,GNUPLOTPROGRAM);
   #ifndef UNIX
     strcat(plotcmd,".exe");
     strcat(plotcmd,"\"");
   #endif
   strcat(plotcmd," ");    strcat(plotcmd," ");
   strcat(plotcmd,optionfilegnuplot);    strcat(plotcmd,optionfilegnuplot);
   printf("Starting graphs with: %s",plotcmd);fflush(stdout);    printf("Starting graphs with: %s",plotcmd);fflush(stdout);
   
   if((outcmd=system(plotcmd)) != 0){    if((outcmd=system(plotcmd)) != 0){
     printf(" Problem with gnuplot\n");      printf("\n Problem with gnuplot\n");
   }    }
   printf(" Wait...");    printf(" Wait...");
   while (z[0] != 'q') {    while (z[0] != 'q') {
Line 5424  Interval (in months) between two waves: Line 5690  Interval (in months) between two waves:
     printf("\nType e to edit output files, g to graph again and q for exiting: ");      printf("\nType e to edit output files, g to graph again and q for exiting: ");
     scanf("%s",z);      scanf("%s",z);
 /*     if (z[0] == 'c') system("./imach"); */  /*     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] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }
Removed from v.1.98  
changed lines
  Added in v.1.110

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>