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

Diff for /imach/src/imach.c between versions 1.73 and 1.85

version 1.73, 2003/04/08 14:06:50 version 1.85, 2003/06/17 13:12:43
Line 1 Line 1
 /* $Id$  /* $Id$
     $State$
     $Log$
     Revision 1.85  2003/06/17 13:12:43  brouard
     * imach.c (Repository): Check when date of death was earlier that
     current date of interview. It may happen when the death was just
     prior to the death. In this case, dh was negative and likelihood
     was wrong (infinity). We still send an "Error" but patch by
     assuming that the date of death was just one stepm after the
     interview.
     (Repository): Because some people have very long ID (first column)
     we changed int to long in num[] and we added a new lvector for
     memory allocation. But we also truncated to 8 characters (left
     truncation)
     (Repository): No more line truncation errors.
   
     Revision 1.84  2003/06/13 21:44:43  brouard
     * imach.c (Repository): Replace "freqsummary" at a correct
     place. It differs from routine "prevalence" which may be called
     many times. Probs is memory consuming and must be used with
     parcimony.
     Version 0.95a2 (should output exactly the same maximization than 0.8a2)
   
     Revision 1.83  2003/06/10 13:39:11  lievre
     *** empty log message ***
   
     Revision 1.82  2003/06/05 15:57:20  brouard
     Add log in  imach.c and  fullversion number is now printed.
   
   */
   /*
    Interpolated Markov Chain     Interpolated Markov Chain
   
   Short summary of the programme:    Short summary of the programme:
Line 48 Line 78
   It is copyrighted identically to a GNU software product, ie programme and    It is copyrighted identically to a GNU software product, ie programme and
   software can be distributed freely for non commercial use. Latest version    software can be distributed freely for non commercial use. Latest version
   can be accessed at http://euroreves.ined.fr/imach .    can be accessed at http://euroreves.ined.fr/imach .
   
     Help to debug: LD_PRELOAD=/usr/local/lib/libnjamd.so ./imach foo.imach
     or better on gdb : set env LD_PRELOAD=/usr/local/lib/libnjamd.so
     
   **********************************************************************/    **********************************************************************/
   /*
     main
     read parameterfile
     read datafile
     concatwav
     freqsummary
     if (mle >= 1)
       mlikeli
     print results files
     if mle==1 
        computes hessian
     read end of parameter file: agemin, agemax, bage, fage, estepm
         begin-prev-date,...
     open gnuplot file
     open html file
     stable prevalence
      for age prevalim()
     h Pij x
     variance of p varprob
     forecasting if prevfcast==1 prevforecast call prevalence()
     health expectancies
     Variance-covariance of DFLE
     prevalence()
      movingaverage()
     varevsij() 
     if popbased==1 varevsij(,popbased)
     total life expectancies
     Variance of stable prevalence
    end
   */
   
   
   
     
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 58 Line 125
 #define MAXLINE 256  #define MAXLINE 256
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 132
 /*#define DEBUG*/  /*#define DEBUG*/
 #define windows  /*#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 75 Line 142
 #define YEARM 12. /* Number of months per year */  #define YEARM 12. /* Number of months per year */
 #define AGESUP 130  #define AGESUP 130
 #define AGEBASE 40  #define AGEBASE 40
 #ifdef windows  #ifdef unix
 #define DIRSEPARATOR '\\'  
 #define ODIRSEPARATOR '/'  
 #else  
 #define DIRSEPARATOR '/'  #define DIRSEPARATOR '/'
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
   #else
   #define DIRSEPARATOR '\\'
   #define ODIRSEPARATOR '/'
 #endif  #endif
   
 char version[80]="Imach version 0.94, February 2003, INED-EUROREVES ";  /* $Id$ */
   /* $State$ */
   
   char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES ";
   char fullversion[]="$Revision$ $Date$"; 
 int erreur; /* Error number */  int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
Line 105  double jmean; /* Mean space between 2 wa Line 176  double jmean; /* Mean space between 2 wa
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficlog;  FILE *ficlog, *ficrespow;
   int globpr; /* Global variable for printing or not */
   double fretone; /* Only one call to likelihood */
   long ipmx; /* Number of contributions */
   double sw; /* Sum of weights */
   char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */
   FILE *ficresilk;
 FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficresprobmorprev;  FILE *ficresprobmorprev;
 FILE *fichtm; /* Html File */  FILE *fichtm; /* Html File */
Line 163  int estepm; Line 240  int estepm;
 /* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/  /* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/
   
 int m,nb;  int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;  long *num;
   int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs;  double **pmmij, ***probs;
 double dateintmean=0;  double dateintmean=0;
Line 188  static int split( char *path, char *dirc Line 266  static int split( char *path, char *dirc
   if ( ss == NULL ) {                   /* no directory, so use current */    if ( ss == NULL ) {                   /* no directory, so use current */
     /*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);*/
 #if     defined(__bsd__)                /* get current working directory */      /* get current working directory */
     extern char *getwd( );      /*    extern  char* getcwd ( char *buf , int len);*/
   
     if ( getwd( dirc ) == NULL ) {  
 #else  
     extern char *getcwd( );  
   
     if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {      if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
 #endif  
       return( GLOCK_ERROR_GETCWD );        return( GLOCK_ERROR_GETCWD );
     }      }
     strcpy( name, path );               /* we've got it */      strcpy( name, path );               /* we've got it */
Line 209  static int split( char *path, char *dirc Line 281  static int split( char *path, char *dirc
     dirc[l1-l2] = 0;                    /* add zero */      dirc[l1-l2] = 0;                    /* add zero */
   }    }
   l1 = strlen( dirc );                  /* length of directory */    l1 = strlen( dirc );                  /* length of directory */
 #ifdef windows    /*#ifdef windows
   if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }    if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
 #else  #else
   if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }    if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
 #endif  #endif
     */
   ss = strrchr( name, '.' );            /* find last / */    ss = strrchr( name, '.' );            /* find last / */
   ss++;    ss++;
   strcpy(ext,ss);                       /* save extension */    strcpy(ext,ss);                       /* save extension */
Line 311  void free_ivector(int *v, long nl, long Line 384  void free_ivector(int *v, long nl, long
   free((FREE_ARG)(v+nl-NR_END));    free((FREE_ARG)(v+nl-NR_END));
 }  }
   
   /************************lvector *******************************/
   long *lvector(long nl,long nh)
   {
     long *v;
     v=(long *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(long)));
     if (!v) nrerror("allocation failure in ivector");
     return v-nl+NR_END;
   }
   
   /******************free lvector **************************/
   void free_lvector(long *v, long nl, long nh)
   {
     free((FREE_ARG)(v+nl-NR_END));
   }
   
 /******************* imatrix *******************************/  /******************* imatrix *******************************/
 int **imatrix(long nrl, long nrh, long ncl, long nch)   int **imatrix(long nrl, long nrh, long ncl, long nch) 
      /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */        /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ 
Line 365  double **matrix(long nrl, long nrh, long Line 453  double **matrix(long nrl, long nrh, long
   
   for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;    for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
   return m;    return m;
     /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1]) 
      */
 }  }
   
 /*************************free matrix ************************/  /*************************free matrix ************************/
Line 404  double ***ma3x(long nrl, long nrh, long Line 494  double ***ma3x(long nrl, long nrh, long
     for (j=ncl+1; j<=nch; j++)       for (j=ncl+1; j<=nch; j++) 
       m[i][j]=m[i][j-1]+nlay;        m[i][j]=m[i][j-1]+nlay;
   }    }
   return m;    return m; 
     /*  gdb: p *(m+1) <=> p m[1] and p (m+1) <=> p (m+1) <=> p &(m[1])
              &(m[i][j][k]) <=> *((*(m+i) + j)+k)
     */
 }  }
   
 /*************************free ma3x ************************/  /*************************free ma3x ************************/
Line 612  void powell(double p[], double **xi, int Line 705  void powell(double p[], double **xi, int
     del=0.0;       del=0.0; 
     printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
     for (i=1;i<=n;i++)       fprintf(ficrespow,"%d %.12f",*iter,*fret);
       for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);        printf(" %d %.12f",i, p[i]);
     fprintf(ficlog," %d %.12f",i, p[i]);        fprintf(ficlog," %d %.12lf",i, p[i]);
         fprintf(ficrespow," %.12lf", p[i]);
       }
     printf("\n");      printf("\n");
     fprintf(ficlog,"\n");      fprintf(ficlog,"\n");
       fprintf(ficrespow,"\n");
     for (i=1;i<=n;i++) {       for (i=1;i<=n;i++) { 
       for (j=1;j<=n;j++) xit[j]=xi[j][i];         for (j=1;j<=n;j++) xit[j]=xi[j][i]; 
       fptt=(*fret);         fptt=(*fret); 
Line 1109  double func( double *x) Line 1206  double func( double *x)
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */        } /* end of wave */
     } /* end of individual */      } /* end of individual */
   }else{  /* ml=4 no inter-extrapolation */    }else if (mle==4){  /* ml=4 no inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];        for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
Line 1131  double func( double *x) Line 1228  double func( double *x)
           oldm=newm;            oldm=newm;
         } /* end mult */          } /* end mult */
               
           s1=s[mw[mi][i]][i];
           s2=s[mw[mi+1][i]][i];
           if( s2 > nlstate){ 
             lli=log(out[s1][s2] - savm[s1][s2]);
           }else{
             lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
           }
           ipmx +=1;
           sw += weight[i];
           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
   /*      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
         } /* end of wave */
       } /* end of individual */
     }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
       for (i=1,ipmx=0, sw=0.; i<=imx; i++){
         for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for(mi=1; mi<= wav[i]-1; mi++){
           for (ii=1;ii<=nlstate+ndeath;ii++)
             for (j=1;j<=nlstate+ndeath;j++){
               oldm[ii][j]=(ii==j ? 1.0 : 0.0);
               savm[ii][j]=(ii==j ? 1.0 : 0.0);
             }
           for(d=0; d<dh[mi][i]; d++){
             newm=savm;
             cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
             for (kk=1; kk<=cptcovage;kk++) {
               cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
             }
           
             out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                          1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             savm=oldm;
             oldm=newm;
           } /* end mult */
         
           s1=s[mw[mi][i]][i];
           s2=s[mw[mi+1][i]][i];
         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 */
         ipmx +=1;          ipmx +=1;
         sw += weight[i];          sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
           /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
       } /* end of wave */        } /* end of wave */
     } /* end of individual */      } /* end of individual */
   } /* End of if */    } /* End of if */
Line 1144  double func( double *x) Line 1279  double func( double *x)
   return -l;    return -l;
 }  }
   
   /*************** log-likelihood *************/
   double funcone( double *x)
   {
     int i, ii, j, k, mi, d, kk;
     double l, ll[NLSTATEMAX], cov[NCOVMAX];
     double **out;
     double lli; /* Individual log likelihood */
     int s1, s2;
     double bbh, survp;
     /*extern weight */
     /* We are differentiating ll according to initial status */
     /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
     /*for(i=1;i<imx;i++) 
       printf(" %d\n",s[4][i]);
     */
     cov[1]=1.;
   
     for(k=1; k<=nlstate; k++) ll[k]=0.;
   
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){
         for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++){
             oldm[ii][j]=(ii==j ? 1.0 : 0.0);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }
         for(d=0; d<dh[mi][i]; d++){
           newm=savm;
           cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
           for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
           }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;
           oldm=newm;
         } /* end mult */
         
         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
          * is higher than the multiple of stepm and negative otherwise.
          */
         if( s2 > nlstate && (mle <5) ){  /* Jackson */
           lli=log(out[s1][s2] - savm[s1][s2]);
         } else if (mle==1){
           lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
         } else if(mle==2){
           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 */
         } else if(mle==3){  /* exponential inter-extrapolation */
           lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
         } else if (mle==4){  /* mle=4 no inter-extrapolation */
           lli=log(out[s1][s2]); /* Original formula */
         } else{  /* ml>=5 no inter-extrapolation no jackson =0.8a */
           lli=log(out[s1][s2]); /* Original formula */
         } /* End of if */
         ipmx +=1;
         sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
   /*       printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
         if(globpr){
           fprintf(ficresilk,"%6d %1d %1d %1d %1d %3d %10.6f %6.4f %10.6f %10.6f %10.6f ", \
                   i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
           for(k=1,l=0.; k<=nlstate; k++) 
             fprintf(ficresilk," %10.6f",ll[k]);
           fprintf(ficresilk,"\n");
         }
       } /* end of wave */
     } /* end of individual */
     for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
     /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
     l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
     return -l;
   }
   
   
   void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpr, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
   {
     /* This routine should help understanding what is done with the selection of individuals/waves and
        to check the exact contribution to the likelihood.
        Plotting could be done.
      */
     int k;
     if(globpr !=0){ /* Just counts and sums no printings */
       strcpy(fileresilk,"ilk"); 
       strcat(fileresilk,fileres);
       if((ficresilk=fopen(fileresilk,"w"))==NULL) {
         printf("Problem with resultfile: %s\n", fileresilk);
         fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
       }
       fprintf(ficresilk, "# individual(line's record) s1 s2 wave# effective_wave# number_of_product_matrix pij weight 2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state");
       fprintf(ficresilk, "# i s1 s2 mi mw dh likeli weight out sav ");
       /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
       for(k=1; k<=nlstate; k++) 
         fprintf(ficresilk," ll[%d]",k);
       fprintf(ficresilk,"\n");
     }
   
     *fretone=(*funcone)(p);
     if(globpr !=0)
       fclose(ficresilk);
     return;
   }
   
 /*********** Maximum Likelihood Estimation ***************/  /*********** Maximum Likelihood Estimation ***************/
   
 void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))  void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
 {  {
   int i,j, iter;    int i,j, iter;
   double **xi,*delti;    double **xi;
   double fret;    double fret;
     double fretone; /* Only one call to likelihood */
     char filerespow[FILENAMELENGTH];
   xi=matrix(1,npar,1,npar);    xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)    for (i=1;i<=npar;i++)
     for (j=1;j<=npar;j++)      for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);        xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");    printf("Powell\n");  fprintf(ficlog,"Powell\n");
     strcpy(filerespow,"pow"); 
     strcat(filerespow,fileres);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
     }
     fprintf(ficrespow,"# Powell\n# iter -2*LL");
     for (i=1;i<=nlstate;i++)
       for(j=1;j<=nlstate+ndeath;j++)
         if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
     fprintf(ficrespow,"\n");
   
   powell(p,xi,npar,ftol,&iter,&fret,func);    powell(p,xi,npar,ftol,&iter,&fret,func);
   
    printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));    fclose(ficrespow);
     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));    fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   
Line 1423  void lubksb(double **a, int n, int *indx Line 1678  void lubksb(double **a, int n, int *indx
 }   } 
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)  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)
 {  /* 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 1435  void  freqsummary(char fileres[], int ag Line 1690  void  freqsummary(char fileres[], int ag
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
       
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   prop=matrix(1,nlstate,agemin,agemax+3);    prop=matrix(1,nlstate,iagemin,iagemax+3);
   probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);  
   strcpy(fileresp,"p");    strcpy(fileresp,"p");
   strcat(fileresp,fileres);    strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {    if((ficresp=fopen(fileresp,"w"))==NULL) {
Line 1444  void  freqsummary(char fileres[], int ag Line 1698  void  freqsummary(char fileres[], int ag
     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,agemin,agemax+3);    freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);
   j1=0;    j1=0;
       
   j=cptcoveff;    j=cptcoveff;
Line 1459  void  freqsummary(char fileres[], int ag Line 1713  void  freqsummary(char fileres[], int ag
         scanf("%d", i);*/          scanf("%d", i);*/
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-1; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=iagemin; m <= iagemax+3; m++)
             freq[i][jk][m]=0;              freq[i][jk][m]=0;
   
     for (i=1; i<=nlstate; i++)        for (i=1; i<=nlstate; i++)  
       for(m=agemin; m <= agemax+3; m++)        for(m=iagemin; m <= iagemax+3; m++)
         prop[i][m]=0;          prop[i][m]=0;
               
       dateintsum=0;        dateintsum=0;
Line 1478  void  freqsummary(char fileres[], int ag Line 1732  void  freqsummary(char fileres[], int ag
         if (bool==1){          if (bool==1){
           for(m=firstpass; m<=lastpass; m++){            for(m=firstpass; m<=lastpass; m++){
             k2=anint[m][i]+(mint[m][i]/12.);              k2=anint[m][i]+(mint[m][i]/12.);
             if ((k2>=dateprev1) && (k2<=dateprev2)) {              /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
               if(agev[m][i]==0) agev[m][i]=agemax+1;                if(agev[m][i]==0) agev[m][i]=iagemax+1;
               if(agev[m][i]==1) agev[m][i]=agemax+2;                if(agev[m][i]==1) agev[m][i]=iagemax+2;
               if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];                if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
               if (m<lastpass) {                if (m<lastpass) {
                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                 freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                  freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];
               }                }
                               
               if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {                if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) {
                 dateintsum=dateintsum+k2;                  dateintsum=dateintsum+k2;
                 k2cpt++;                  k2cpt++;
               }                }
             }                /*}*/
           }            }
         }          }
       }        }
                 
       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);*/
   
       if  (cptcovn>0) {        if  (cptcovn>0) {
         fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable "); 
Line 1507  void  freqsummary(char fileres[], int ag Line 1761  void  freqsummary(char fileres[], int ag
         fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
       fprintf(ficresp, "\n");        fprintf(ficresp, "\n");
               
       for(i=(int)agemin; i <= (int)agemax+3; i++){        for(i=iagemin; i <= iagemax+3; i++){
         if(i==(int)agemax+3){          if(i==iagemax+3){
           fprintf(ficlog,"Total");            fprintf(ficlog,"Total");
         }else{          }else{
           if(first==1){            if(first==1){
Line 1554  void  freqsummary(char fileres[], int ag Line 1808  void  freqsummary(char fileres[], int ag
               printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);                printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
             fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);              fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
           }            }
           if( i <= (int) agemax){            if( i <= iagemax){
             if(pos>=1.e-5){              if(pos>=1.e-5){
               fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);                fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
               probs[i][jk][j1]= pp[jk]/pos;                /*probs[i][jk][j1]= pp[jk]/pos;*/
               /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/                /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
             }              }
             else              else
Line 1572  void  freqsummary(char fileres[], int ag Line 1826  void  freqsummary(char fileres[], int ag
               printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);                printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
               fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);                fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
             }              }
         if(i <= (int) agemax)          if(i <= iagemax)
           fprintf(ficresp,"\n");            fprintf(ficresp,"\n");
         if(first==1)          if(first==1)
           printf("Others in log...\n");            printf("Others in log...\n");
Line 1583  void  freqsummary(char fileres[], int ag Line 1837  void  freqsummary(char fileres[], int ag
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3);    free_matrix(prop,1,nlstate,iagemin, iagemax+3);
   /* End of Freq */    /* End of Freq */
 }  }
   
 /************ Prevalence ********************/  /************ Prevalence ********************/
 void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)  void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
 {    {  
   /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people    /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
      in each health status at the date of interview (if between dateprev1 and dateprev2).       in each health status at the date of interview (if between dateprev1 and dateprev2).
Line 1602  void prevalence(int agemin, float agemax Line 1856  void prevalence(int agemin, float agemax
   double *pp, **prop;    double *pp, **prop;
   double pos,posprop;     double pos,posprop; 
   double  y2; /* in fractional years */    double  y2; /* in fractional years */
     int iagemin, iagemax;
   
   pp=vector(1,nlstate);    iagemin= (int) agemin;
   prop=matrix(1,nlstate,agemin,agemax+3);     iagemax= (int) agemax;
   freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    /*pp=vector(1,nlstate);*/
     prop=matrix(1,nlstate,iagemin,iagemax+3); 
     /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
   j1=0;    j1=0;
       
   j=cptcoveff;    j=cptcoveff;
Line 1616  void prevalence(int agemin, float agemax Line 1873  void prevalence(int agemin, float agemax
       j1++;        j1++;
               
       for (i=1; i<=nlstate; i++)          for (i=1; i<=nlstate; i++)  
         for(m=agemin; m <= agemax+3; m++)          for(m=iagemin; m <= iagemax+3; m++)
           prop[i][m]=0;            prop[i][m]=0.0;
             
       for (i=1; i<=imx; i++) { /* Each individual */        for (i=1; i<=imx; i++) { /* Each individual */
         bool=1;          bool=1;
Line 1630  void prevalence(int agemin, float agemax Line 1887  void prevalence(int agemin, float agemax
           for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/            for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/
             y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */              y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
             if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */              if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
               if(agev[m][i]==0) agev[m][i]=agemax+1;                if(agev[m][i]==0) agev[m][i]=iagemax+1;
               if(agev[m][i]==1) agev[m][i]=agemax+2;                if(agev[m][i]==1) agev[m][i]=iagemax+2;
               if (s[m][i]>0 && s[m][i]<=nlstate) {                if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); 
                 prop[s[m][i]][(int)agev[m][i]] += weight[i];                if (s[m][i]>0 && s[m][i]<=nlstate) { 
                 prop[s[m][i]][(int)(agemax+3)] += weight[i];                  /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
               }                  prop[s[m][i]][(int)agev[m][i]] += weight[i];
                   prop[s[m][i]][iagemax+3] += weight[i]; 
                 } 
             }              }
           } /* end selection of waves */            } /* end selection of waves */
         }          }
       }        }
       for(i=(int)agemin; i <= (int)agemax+3; i++){         for(i=iagemin; i <= iagemax+3; i++){  
                   
         for(jk=1,posprop=0; jk <=nlstate ; jk++) {          for(jk=1,posprop=0; jk <=nlstate ; jk++) { 
           posprop += prop[jk][i];            posprop += prop[jk][i]; 
         }          } 
           
         for(jk=1; jk <=nlstate ; jk++){              for(jk=1; jk <=nlstate ; jk++){     
           if( i <= (int) agemax){            if( i <=  iagemax){ 
             if(posprop>=1.e-5){              if(posprop>=1.e-5){ 
              probs[i][jk][j1]= prop[jk][i]/posprop;                probs[i][jk][j1]= prop[jk][i]/posprop;
             }              } 
           }            } 
         }/* end jk */          }/* end jk */ 
       }/* end i */        }/* end i */ 
     } /* end i1 */      } /* end i1 */
   } /* end k1 */    } /* end k1 */
   
       
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
   free_vector(pp,1,nlstate);    /*free_vector(pp,1,nlstate);*/
   free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3);    free_matrix(prop,1,nlstate, iagemin,iagemax+3);
 }  /* End of Freq */  }  /* End of prevalence */
   
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
   
Line 1705  void  concatwav(int wav[], int **dh, int Line 1963  void  concatwav(int wav[], int **dh, int
     wav[i]=mi;      wav[i]=mi;
     if(mi==0){      if(mi==0){
       if(first==0){        if(first==0){
         printf("Warning, no any valid information for:%d line=%d and may be others, see log file\n",num[i],i);          printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);
         first=1;          first=1;
       }        }
       if(first==1){        if(first==1){
         fprintf(ficlog,"Warning, no any valid information for:%d line=%d\n",num[i],i);          fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i);
       }        }
     } /* end mi==0 */      } /* end mi==0 */
   }    } /* End individuals */
   
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     for(mi=1; mi<wav[i];mi++){      for(mi=1; mi<wav[i];mi++){
       if (stepm <=0)        if (stepm <=0)
         dh[mi][i]=1;          dh[mi][i]=1;
       else{        else{
         if (s[mw[mi+1][i]][i] > nlstate) {          if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
           if (agedc[i] < 2*AGESUP) {            if (agedc[i] < 2*AGESUP) {
           j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);               j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
           if(j==0) j=1;  /* Survives at least one month after exam */              if(j==0) j=1;  /* Survives at least one month after exam */
           k=k+1;              else if(j<0){
           if (j >= jmax) jmax=j;                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]);
           if (j <= jmin) jmin=j;                j=1; /* Careful Patch */
           sum=sum+j;                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 fixe the contradiction between dates.\n",stepm);
           /*if (j<0) printf("j=%d num=%d \n",j,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]);
           /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,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);
           /*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]);*/              }
               k=k+1;
               if (j >= jmax) jmax=j;
               if (j <= jmin) jmin=j;
               sum=sum+j;
               /*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);*/
           }            }
         }          }
         else{          else{
Line 1740  void  concatwav(int wav[], int **dh, int Line 2004  void  concatwav(int wav[], int **dh, int
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
           /*        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){
               printf("Error! Negative delay (%d) 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) 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]);
             }
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
         jl= j -jk*stepm;          jl= j -jk*stepm;
         ju= j -(jk+1)*stepm;          ju= j -(jk+1)*stepm;
         if(mle <=1){           if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
           if(jl==0){            if(jl==0){
             dh[mi][i]=jk;              dh[mi][i]=jk;
             bh[mi][i]=0;              bh[mi][i]=0;
Line 1770  void  concatwav(int wav[], int **dh, int Line 2038  void  concatwav(int wav[], int **dh, int
             bh[mi][i]=ju; /* At least one step */              bh[mi][i]=ju; /* At least one step */
             /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/              /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
           }            }
         }          } /* end if mle */
       } /* end if mle */        }
     } /* end wave */      } /* end wave */
   }    }
   jmean=sum/k;    jmean=sum/k;
Line 1853  void evsij(char fileres[], double ***eij Line 2121  void evsij(char fileres[], double ***eij
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   int theta;    int theta;
   
   varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);    varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage);
   xp=vector(1,npar);    xp=vector(1,npar);
   dnewm=matrix(1,nlstate*2,1,npar);    dnewm=matrix(1,nlstate*nlstate,1,npar);
   doldm=matrix(1,nlstate*2,1,nlstate*2);    doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);
       
   fprintf(ficreseij,"# Health expectancies\n");    fprintf(ficreseij,"# Health expectancies\n");
   fprintf(ficreseij,"# Age");    fprintf(ficreseij,"# Age");
Line 1902  void evsij(char fileres[], double ***eij Line 2170  void evsij(char fileres[], double ***eij
     /* if (stepm >= YEARM) hstepm=1;*/      /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);      gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);
     gp=matrix(0,nhstepm,1,nlstate*2);      gp=matrix(0,nhstepm,1,nlstate*nlstate);
     gm=matrix(0,nhstepm,1,nlstate*2);      gm=matrix(0,nhstepm,1,nlstate*nlstate);
   
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
Line 1941  void evsij(char fileres[], double ***eij Line 2209  void evsij(char fileres[], double ***eij
         for(i=1;i<=nlstate;i++){          for(i=1;i<=nlstate;i++){
           cptj=cptj+1;            cptj=cptj+1;
           for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){            for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
   
             gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;              gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
           }            }
         }          }
       }        }
       for(j=1; j<= nlstate*2; j++)        for(j=1; j<= nlstate*nlstate; j++)
         for(h=0; h<=nhstepm-1; h++){          for(h=0; h<=nhstepm-1; h++){
           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
         }          }
Line 1953  void evsij(char fileres[], double ***eij Line 2222  void evsij(char fileres[], double ***eij
         
 /* End theta */  /* End theta */
   
      trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);       trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
   
      for(h=0; h<=nhstepm-1; h++)       for(h=0; h<=nhstepm-1; h++)
       for(j=1; j<=nlstate*2;j++)        for(j=1; j<=nlstate*nlstate;j++)
         for(theta=1; theta <=npar; theta++)          for(theta=1; theta <=npar; theta++)
           trgradg[h][j][theta]=gradg[h][theta][j];            trgradg[h][j][theta]=gradg[h][theta][j];
             
   
      for(i=1;i<=nlstate*2;i++)       for(i=1;i<=nlstate*nlstate;i++)
       for(j=1;j<=nlstate*2;j++)        for(j=1;j<=nlstate*nlstate;j++)
         varhe[i][j][(int)age] =0.;          varhe[i][j][(int)age] =0.;
   
      printf("%d|",(int)age);fflush(stdout);       printf("%d|",(int)age);fflush(stdout);
      fprintf(ficlog,"%d|",(int)age);fflush(ficlog);       fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
      for(h=0;h<=nhstepm-1;h++){       for(h=0;h<=nhstepm-1;h++){
       for(k=0;k<=nhstepm-1;k++){        for(k=0;k<=nhstepm-1;k++){
         matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);          matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
         for(i=1;i<=nlstate*2;i++)          for(i=1;i<=nlstate*nlstate;i++)
           for(j=1;j<=nlstate*2;j++)            for(j=1;j<=nlstate*nlstate;j++)
             varhe[i][j][(int)age] += doldm[i][j]*hf*hf;              varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
       }        }
     }      }
Line 1995  void evsij(char fileres[], double ***eij Line 2264  void evsij(char fileres[], double ***eij
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficreseij,"\n");
         
     free_matrix(gm,0,nhstepm,1,nlstate*2);      free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
     free_matrix(gp,0,nhstepm,1,nlstate*2);      free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
     free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);      free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
     free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);      free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }    }
   printf("\n");    printf("\n");
   fprintf(ficlog,"\n");    fprintf(ficlog,"\n");
   
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(dnewm,1,nlstate*2,1,npar);    free_matrix(dnewm,1,nlstate*nlstate,1,npar);
   free_matrix(doldm,1,nlstate*2,1,nlstate*2);    free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
   free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);    free_ma3x(varhe,1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int)fage);
 }  }
   
 /************ Variance ******************/  /************ Variance ******************/
Line 2318  void varevsij(char optionfilefiname[], d Line 2587  void varevsij(char optionfilefiname[], d
   fclose(ficresprobmorprev);    fclose(ficresprobmorprev);
   fclose(ficgp);    fclose(ficgp);
   fclose(fichtm);    fclose(fichtm);
 }    }  /* 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)
Line 2542  void varprob(char optionfilefiname[], do Line 2811  void varprob(char optionfilefiname[], do
           
         for(theta=1; theta <=npar; theta++){          for(theta=1; theta <=npar; theta++){
           for(i=1; i<=npar; i++)            for(i=1; i<=npar; i++)
             xp[i] = x[i] + (i==theta ?delti[theta]:0);              xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);
                       
           pmij(pmmij,cov,ncovmodel,xp,nlstate);            pmij(pmmij,cov,ncovmodel,xp,nlstate);
                       
Line 2555  void varprob(char optionfilefiname[], do Line 2824  void varprob(char optionfilefiname[], do
           }            }
                       
           for(i=1; i<=npar; i++)            for(i=1; i<=npar; i++)
             xp[i] = x[i] - (i==theta ?delti[theta]:0);              xp[i] = x[i] - (i==theta ?delti[theta]:(double)0);
           
           pmij(pmmij,cov,ncovmodel,xp,nlstate);            pmij(pmmij,cov,ncovmodel,xp,nlstate);
           k=0;            k=0;
Line 2567  void varprob(char optionfilefiname[], do Line 2836  void varprob(char optionfilefiname[], do
           }            }
             
           for(i=1; i<= (nlstate)*(nlstate+ndeath); i++)             for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) 
             gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];                gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta];  
         }          }
   
         for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)          for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
Line 2730  void printinghtml(char fileres[], char t Line 2999  void printinghtml(char fileres[], char t
     fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);      fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
    fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n     fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n   - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \
  - Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n   - Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n \
  - Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n   - Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n \
  - Life expectancies by age and initial health status (estepm=%2d months):    - Life expectancies by age and initial health status (estepm=%2d months): \
    <a href=\"e%s\">e%s</a> <br>\n</li>", \     <a href=\"e%s\">e%s</a> <br>\n</li>", \
   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);    jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
   
Line 2754  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3023  fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: pe%s%d1.png<br> \
 <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);       <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
    before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br> \
 <img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);   <img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); 
        /* Stable prevalence in each health state */         /* Stable prevalence in each health state */
        for(cpt=1; cpt<nlstate;cpt++){         for(cpt=1; cpt<nlstate;cpt++){
          fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>           fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br> \
 <img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }         }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>          fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br> \
 <img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2): e%s%d.png<br>  health expectancies in states (1) and (2): e%s%d.png<br>\
 <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);  <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
   
   
  fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n   fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n\
  - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n   - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n\
  - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n   - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n\
  - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n   - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n\
  - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n    - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n\
  - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n   - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n\
  - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);   - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
   
  if(popforecast==1) fprintf(fichtm,"\n  /*  if(popforecast==1) fprintf(fichtm,"\n */
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n  /*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n  /*  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n */
         <br>",fileres,fileres,fileres,fileres);  /*      <br>",fileres,fileres,fileres,fileres); */
  else   /*  else  */
    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);  /*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
 fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
   
  m=cptcoveff;   m=cptcoveff;
Line 2807  fprintf(fichtm," <ul><li><b>Graphs</b></ Line 3077  fprintf(fichtm," <ul><li><b>Graphs</b></
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident         fprintf(fichtm,"<br>- Observed and period prevalence (with confident\
 interval) in state (%d): v%s%d%d.png <br>  interval) in state (%d): v%s%d%d.png <br>\
 <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);    <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }       }
    } /* end i1 */     } /* end i1 */
Line 2908  m=pow(2,cptcoveff); Line 3178  m=pow(2,cptcoveff);
     }      }
   }    }
       
   /* CV preval stat */    /* CV preval stable (period) */
   for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<nlstate ; cpt ++) {      for (cpt=1; cpt<=nlstate ; cpt ++) {
       k=3;        k=3;
       fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);        fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
Line 3054  prevforecast(char fileres[], double anpr Line 3324  prevforecast(char fileres[], double anpr
   char fileresf[FILENAMELENGTH];    char fileresf[FILENAMELENGTH];
   
   agelim=AGESUP;    agelim=AGESUP;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);    prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     
   strcpy(fileresf,"f");     strcpy(fileresf,"f"); 
   strcat(fileresf,fileres);    strcat(fileresf,fileres);
Line 3077  prevforecast(char fileres[], double anpr Line 3347  prevforecast(char fileres[], double anpr
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
   if (stepm<=12) stepsize=1;    if (stepm<=12) stepsize=1;
       if(estepm < stepm){
   hstepm=1;      printf ("Problem %d lower than %d\n",estepm, stepm);
     }
     else  hstepm=estepm;   
   
   hstepm=hstepm/stepm;     hstepm=hstepm/stepm; 
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and    yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
                                fractional in yp1 */                                 fractional in yp1 */
Line 3097  prevforecast(char fileres[], double anpr Line 3370  prevforecast(char fileres[], double anpr
       
   fprintf(ficresf,"#****** Routine prevforecast **\n");    fprintf(ficresf,"#****** Routine prevforecast **\n");
   
   /*            if (h==(int)(YEARM*yearp)){ */
   for(cptcov=1, k=0;cptcov<=i1;cptcov++){    for(cptcov=1, k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;        k=k+1;
Line 3111  prevforecast(char fileres[], double anpr Line 3385  prevforecast(char fileres[], double anpr
           fprintf(ficresf," p%d%d",i,j);            fprintf(ficresf," p%d%d",i,j);
         fprintf(ficresf," p.%d",j);          fprintf(ficresf," p.%d",j);
       }        }
       for (yearp=0; yearp<=(anproj2-anproj1);yearp++) {         for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { 
         fprintf(ficresf,"\n");          fprintf(ficresf,"\n");
         fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);             fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
   
Line 3123  prevforecast(char fileres[], double anpr Line 3397  prevforecast(char fileres[], double anpr
           hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);              hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);  
                   
           for (h=0; h<=nhstepm; h++){            for (h=0; h<=nhstepm; h++){
             if (h==(int) (YEARM*yearp)) {              if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");                fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++)                 for(j=1;j<=cptcoveff;j++) 
                 fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);                  fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
Line 3137  prevforecast(char fileres[], double anpr Line 3411  prevforecast(char fileres[], double anpr
                 else {                  else {
                   ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];                    ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
                 }                  }
                 if (h==(int)(YEARM*yearp))                  if (h*hstepm/YEARM*stepm== yearp) {
                   fprintf(ficresf," %.3f", p3mat[i][j][h]);                    fprintf(ficresf," %.3f", p3mat[i][j][h]);
               }                  }
               if (h==(int)(YEARM*yearp)){                } /* end i */
                 if (h*hstepm/YEARM*stepm==yearp) {
                 fprintf(ficresf," %.3f", ppij);                  fprintf(ficresf," %.3f", ppij);
               }                }
             }              }/* end j */
           }            } /* end h */
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         }          } /* end agec */
       }        } /* end yearp */
     }      } /* end cptcod */
   }    } /* end  cptcov */
                 
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
Line 3172  populforecast(char fileres[], double anp Line 3447  populforecast(char fileres[], double anp
   agelim=AGESUP;    agelim=AGESUP;
   calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;    calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
       
   prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);    prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
       
       
   strcpy(filerespop,"pop");     strcpy(filerespop,"pop"); 
Line 3314  populforecast(char fileres[], double anp Line 3589  populforecast(char fileres[], double anp
   free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   fclose(ficrespop);    fclose(ficrespop);
 }  } /* End of popforecast */
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
Line 3323  populforecast(char fileres[], double anp Line 3598  populforecast(char fileres[], double anp
 int main(int argc, char *argv[])  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,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
     int jj;
     int numlinepar=0; /* Current linenumber of parameter file */
   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 3335  int main(int argc, char *argv[]) Line 3612  int main(int argc, char *argv[])
   double ***mobaverage;    double ***mobaverage;
   int *indx;    int *indx;
   char line[MAXLINE], linepar[MAXLINE];    char line[MAXLINE], linepar[MAXLINE];
   char path[80],pathc[80],pathcd[80],pathtot[80],model[80];    char path[132],pathc[132],pathcd[132],pathtot[132],model[132];
   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 3345  int main(int argc, char *argv[]) Line 3622  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;
   double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;    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 bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 3360  int main(int argc, char *argv[]) Line 3638  int main(int argc, char *argv[])
   double **varpl; /* Variances of prevalence limits by age */    double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
   
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
Line 3369  int main(int argc, char *argv[]) Line 3647  int main(int argc, char *argv[])
 #include <sys/time.h>  #include <sys/time.h>
 #include <time.h>  #include <time.h>
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];    char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
     char *strt, *strtend;
     char *stratrunc;
     int lstra;
   
     long total_usecs;
     struct timeval start_time, end_time, curr_time;
     struct timezone tzp;
     extern int gettimeofday();
     struct tm tmg, tm, *gmtime(), *localtime();
     long time_value;
     extern long time();
     
   /* long total_usecs;    /*   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
      struct timeval start_time, end_time;    (void) gettimeofday(&start_time,&tzp);
       tm = *localtime(&start_time.tv_sec);
      gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    tmg = *gmtime(&start_time.tv_sec);
     strt=asctime(&tm);
   /*  printf("Localtime (at start)=%s",strt); */
   /*  tp.tv_sec = tp.tv_sec +86400; */
   /*  tm = *localtime(&start_time.tv_sec); */
   /*   tmg.tm_year=tmg.tm_year +dsign*dyear; */
   /*   tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */
   /*   tmg.tm_hour=tmg.tm_hour + 1; */
   /*   tp.tv_sec = mktime(&tmg); */
   /*   strt=asctime(&tmg); */
   /*   printf("Time(after) =%s",strt);  */
   /*  (void) time (&time_value);
   *  printf("time=%d,t-=%d\n",time_value,time_value-86400);
   *  tm = *localtime(&time_value);
   *  strt=asctime(&tm);
   *  printf("tim_value=%d,asctime=%s\n",time_value,strt); 
   */
   
   getcwd(pathcd, size);    getcwd(pathcd, size);
   
   printf("\n%s",version);    printf("\n%s\n%s",version,fullversion);
   if(argc <=1){    if(argc <=1){
     printf("\nEnter the parameter file name: ");      printf("\nEnter the parameter file name: ");
     scanf("%s",pathtot);      scanf("%s",pathtot);
Line 3384  int main(int argc, char *argv[]) Line 3690  int main(int argc, char *argv[])
   else{    else{
     strcpy(pathtot,argv[1]);      strcpy(pathtot,argv[1]);
   }    }
   /*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/    /*if(getcwd(pathcd, 132)!= NULL)printf ("Error pathcd\n");*/
   /*cygwin_split_path(pathtot,path,optionfile);    /*cygwin_split_path(pathtot,path,optionfile);
     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(pathtot,path,optionfile,optionfilext,optionfilefiname);    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);    chdir(path);
   replace(pathc,path);    replace(pathc,path);
   
Line 3404  int main(int argc, char *argv[]) Line 3710  int main(int argc, char *argv[])
     goto end;      goto end;
   }    }
   fprintf(ficlog,"Log filename:%s\n",filelog);    fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s",version);    fprintf(ficlog,"\n%s\n%s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: ");    fprintf(ficlog,"\nEnter the parameter file name: ");
   fprintf(ficlog,"pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    fprintf(ficlog,"pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
     printf("Localtime (at start)=%s",strt); 
     fprintf(ficlog,"Localtime (at start)=%s",strt); 
   fflush(ficlog);    fflush(ficlog);
   
   /* */    /* */
Line 3419  int main(int argc, char *argv[]) Line 3727  int main(int argc, char *argv[])
   if((ficpar=fopen(optionfile,"r"))==NULL)    {    if((ficpar=fopen(optionfile,"r"))==NULL)    {
     printf("Problem with optionfile %s\n",optionfile);      printf("Problem with optionfile %s\n",optionfile);
     fprintf(ficlog,"Problem with optionfile %s\n",optionfile);      fprintf(ficlog,"Problem with optionfile %s\n",optionfile);
       fflush(ficlog);
     goto end;      goto end;
   }    }
   
Line 3427  int main(int argc, char *argv[]) Line 3736  int main(int argc, char *argv[])
   if((ficparo=fopen(filereso,"w"))==NULL) {    if((ficparo=fopen(filereso,"w"))==NULL) {
     printf("Problem with Output resultfile: %s\n", filereso);      printf("Problem with Output resultfile: %s\n", filereso);
     fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);      fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
       fflush(ficlog);
     goto end;      goto end;
   }    }
   
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
     numlinepar=0;
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
       numlinepar++;
     puts(line);      puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
       fputs(line,ficlog);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);    fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
     numlinepar++;
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);    fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
     fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
     fflush(ficlog);
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
       numlinepar++;
     puts(line);      puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
       fputs(line,ficlog);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
     
         
   covar=matrix(0,NCOVMAX,1,n);     covar=matrix(0,NCOVMAX,1,n); 
   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/    cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/
Line 3463  int main(int argc, char *argv[]) Line 3781  int main(int argc, char *argv[])
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
       numlinepar++;
     puts(line);      puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
       fputs(line,ficlog);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
     
   param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);    param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
   for(i=1; i <=nlstate; i++)    for(i=1; i <=nlstate; i++){
     for(j=1; j <=nlstate+ndeath-1; j++){      j=0;
       for(jj=1; jj <=nlstate+ndeath; jj++){
         if(jj==i) continue;
         j++;
       fscanf(ficpar,"%1d%1d",&i1,&j1);        fscanf(ficpar,"%1d%1d",&i1,&j1);
         if ((i1 != i) && (j1 != j)){
           printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
           exit(1);
         }
       fprintf(ficparo,"%1d%1d",i1,j1);        fprintf(ficparo,"%1d%1d",i1,j1);
       if(mle==1)        if(mle==1)
         printf("%1d%1d",i,j);          printf("%1d%1d",i,j);
Line 3487  int main(int argc, char *argv[]) Line 3814  int main(int argc, char *argv[])
         fprintf(ficparo," %lf",param[i][j][k]);          fprintf(ficparo," %lf",param[i][j][k]);
       }        }
       fscanf(ficpar,"\n");        fscanf(ficpar,"\n");
         numlinepar++;
       if(mle==1)        if(mle==1)
         printf("\n");          printf("\n");
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
       }  
     fflush(ficlog);
   
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/    npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
   
   p=param[1][1];    p=param[1][1];
Line 3501  int main(int argc, char *argv[]) Line 3831  int main(int argc, char *argv[])
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
       numlinepar++;
     puts(line);      puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
       fputs(line,ficlog);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   
   delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);    delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
   delti=vector(1,npar); /* Scale of each paramater (output from hesscov) */    /* delti=vector(1,npar); *//* Scale of each paramater (output from hesscov) */
   for(i=1; i <=nlstate; i++){    for(i=1; i <=nlstate; i++){
     for(j=1; j <=nlstate+ndeath-1; j++){      for(j=1; j <=nlstate+ndeath-1; j++){
       fscanf(ficpar,"%1d%1d",&i1,&j1);        fscanf(ficpar,"%1d%1d",&i1,&j1);
         if ((i1-i)*(j1-j)!=0){
           printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
           exit(1);
         }
       printf("%1d%1d",i,j);        printf("%1d%1d",i,j);
       fprintf(ficparo,"%1d%1d",i1,j1);        fprintf(ficparo,"%1d%1d",i1,j1);
         fprintf(ficlog,"%1d%1d",i1,j1);
       for(k=1; k<=ncovmodel;k++){        for(k=1; k<=ncovmodel;k++){
         fscanf(ficpar,"%le",&delti3[i][j][k]);          fscanf(ficpar,"%le",&delti3[i][j][k]);
         printf(" %le",delti3[i][j][k]);          printf(" %le",delti3[i][j][k]);
         fprintf(ficparo," %le",delti3[i][j][k]);          fprintf(ficparo," %le",delti3[i][j][k]);
           fprintf(ficlog," %le",delti3[i][j][k]);
       }        }
       fscanf(ficpar,"\n");        fscanf(ficpar,"\n");
         numlinepar++;
       printf("\n");        printf("\n");
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
         fprintf(ficlog,"\n");
     }      }
   }    }
     fflush(ficlog);
   
   delti=delti3[1][1];    delti=delti3[1][1];
   
   
     /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */
       
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
       numlinepar++;
     puts(line);      puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
       fputs(line,ficlog);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
       
Line 3545  int main(int argc, char *argv[]) Line 3892  int main(int argc, char *argv[])
       fscanf(ficpar," %le",&matcov[i][j]);        fscanf(ficpar," %le",&matcov[i][j]);
       if(mle==1){        if(mle==1){
         printf(" %.5le",matcov[i][j]);          printf(" %.5le",matcov[i][j]);
         fprintf(ficlog," %.5le",matcov[i][j]);  
       }        }
       else        fprintf(ficlog," %.5le",matcov[i][j]);
         fprintf(ficlog," %.5le",matcov[i][j]);  
       fprintf(ficparo," %.5le",matcov[i][j]);        fprintf(ficparo," %.5le",matcov[i][j]);
     }      }
     fscanf(ficpar,"\n");      fscanf(ficpar,"\n");
       numlinepar++;
     if(mle==1)      if(mle==1)
       printf("\n");        printf("\n");
     fprintf(ficlog,"\n");      fprintf(ficlog,"\n");
Line 3565  int main(int argc, char *argv[]) Line 3911  int main(int argc, char *argv[])
     printf("\n");      printf("\n");
   fprintf(ficlog,"\n");    fprintf(ficlog,"\n");
   
     fflush(ficlog);
   
   /*-------- Rewriting paramater file ----------*/    /*-------- Rewriting paramater file ----------*/
   strcpy(rfileres,"r");    /* "Rparameterfile */    strcpy(rfileres,"r");    /* "Rparameterfile */
Line 3586  int main(int argc, char *argv[]) Line 3933  int main(int argc, char *argv[])
   n= lastobs;    n= lastobs;
   severity = vector(1,maxwav);    severity = vector(1,maxwav);
   outcome=imatrix(1,maxwav+1,1,n);    outcome=imatrix(1,maxwav+1,1,n);
   num=ivector(1,n);    num=lvector(1,n);
   moisnais=vector(1,n);    moisnais=vector(1,n);
   annais=vector(1,n);    annais=vector(1,n);
   moisdc=vector(1,n);    moisdc=vector(1,n);
Line 3622  int main(int argc, char *argv[]) Line 3969  int main(int argc, char *argv[])
       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,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
       }         } 
       num[i]=atol(stra);        lstra=strlen(stra);
         if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
           stratrunc = &(stra[lstra-9]);
           num[i]=atol(stratrunc);
         }
         else
           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("%d %.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;
     }      }
Line 3641  int main(int argc, char *argv[]) Line 3994  int main(int argc, char *argv[])
     }*/      }*/
    /*  for (i=1; i<=imx; i++){     /*  for (i=1; i<=imx; i++){
      if (s[4][i]==9)  s[4][i]=-1;        if (s[4][i]==9)  s[4][i]=-1; 
      printf("%d %.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++)
     
Line 3739  int main(int argc, char *argv[]) Line 4092  int main(int argc, char *argv[])
   
   for (i=1; i<=imx; i++) {    for (i=1; i<=imx; i++) {
     for(m=2; (m<= maxwav); m++) {      for(m=2; (m<= maxwav); m++) {
       if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){        if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
         anint[m][i]=9999;          anint[m][i]=9999;
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
           printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
           fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
           s[m][i]=-1;
         }
         if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
           printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
           fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
           s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
         }
     }      }
   }    }
   
Line 3753  int main(int argc, char *argv[]) Line 4115  int main(int argc, char *argv[])
       if(s[m][i] >0){        if(s[m][i] >0){
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0)            if(agedc[i]>0)
             if(moisdc[i]!=99 && andc[i]!=9999)              if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/            /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             else {              else {
               if (andc[i]!=9999){                if ((int)andc[i]!=9999){
                 printf("Warning negative age at death: %d line:%d\n",num[i],i);                  printf("Warning negative age at death: %ld line:%d\n",num[i],i);
                 fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);                  fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);
                 agev[m][i]=-1;                  agev[m][i]=-1;
               }                }
             }              }
Line 3768  int main(int argc, char *argv[]) Line 4130  int main(int argc, char *argv[])
                                  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(mint[m][i]==99 || anint[m][i]==9999)            if((int)mint[m][i]==99 || (int)anint[m][i]==9999)
             agev[m][i]=1;              agev[m][i]=1;
           else if(agev[m][i] <agemin){             else if(agev[m][i] <agemin){ 
             agemin=agev[m][i];              agemin=agev[m][i];
Line 3803  int main(int argc, char *argv[]) Line 4165  int main(int argc, char *argv[])
   
   /*for (i=1; i<=imx; i++){    /*for (i=1; i<=imx; i++){
   for (m=firstpass; (m<lastpass); m++){    for (m=firstpass; (m<lastpass); m++){
      printf("%d %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);       printf("%ld %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
 }  }
   
 }*/  }*/
Line 3864  int main(int argc, char *argv[]) Line 4226  int main(int argc, char *argv[])
           
   /* 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);
   
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
Line 3876  int main(int argc, char *argv[]) Line 4239  int main(int argc, char *argv[])
      so we point p on param[1][1] so that p[1] maps on param[1][1][1] */       so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
   
     globpr=0; /* To get ipmx number of contributions and sum of weights*/
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
     for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);
     printf("\n");
     globpr=1; /* to print the contributions */
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
     for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);
     printf("\n");
   if(mle>=1){ /* Could be 1 or 2 */    if(mle>=1){ /* Could be 1 or 2 */
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);      mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   }    }
Line 3907  int main(int argc, char *argv[]) Line 4282  int main(int argc, char *argv[])
         }          }
     }      }
   }    }
   if(mle==1){    if(mle!=0){
     /* Computing hessian and covariance matrix */      /* Computing hessian and covariance matrix */
     ftolhess=ftol; /* Usually correct */      ftolhess=ftol; /* Usually correct */
     hesscov(matcov, p, npar, delti, ftolhess, func);      hesscov(matcov, p, npar, delti, ftolhess, func);
Line 4038  int main(int argc, char *argv[]) Line 4413  int main(int argc, char *argv[])
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);    fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);    fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
     /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
   
   /*------------ gnuplot -------------*/    /*------------ gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);    strcpy(optionfilegnuplot,optionfilefiname);
Line 4061  int main(int argc, char *argv[]) Line 4437  int main(int argc, char *argv[])
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
   fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n    fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\
 \n  \n\
 Total number of observations=%d <br>\n  Total number of observations=%d <br>\n\
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n  Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\
 <hr  size=\"2\" color=\"#EC5E5E\">  Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n\
  <ul><li><h4>Parameter files</h4>\n  <hr  size=\"2\" color=\"#EC5E5E\">\
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n   <ul><li><h4>Parameter files</h4>\n\
  - Log file of the run: <a href=\"%s\">%s</a><br>\n   - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\
  - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);   - Log file of the run: <a href=\"%s\">%s</a><br>\n\
    - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",\
             version,title,datafile,firstpass,lastpass,stepm, weightopt,\
             model,imx,agemin,agemax,jmin,jmax,jmean,fileres,fileres,\
             filelog,filelog,optionfilegnuplot,optionfilegnuplot);
   fclose(fichtm);    fclose(fichtm);
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
                  jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
     
   /*------------ free_vector  -------------*/    /*------------ free_vector  -------------*/
   chdir(path);    chdir(path);
Line 4082  Interval (in months) between two waves: Line 4464  Interval (in months) between two waves:
   free_imatrix(dh,1,lastpass-firstpass+1,1,imx);    free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
   free_imatrix(bh,1,lastpass-firstpass+1,1,imx);    free_imatrix(bh,1,lastpass-firstpass+1,1,imx);
   free_imatrix(mw,1,lastpass-firstpass+1,1,imx);       free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   
   free_ivector(num,1,n);    free_lvector(num,1,n);
   free_vector(agedc,1,n);    free_vector(agedc,1,n);
   /*free_matrix(covar,0,NCOVMAX,1,n);*/    /*free_matrix(covar,0,NCOVMAX,1,n);*/
   /*free_matrix(covar,1,NCOVMAX,1,n);*/    /*free_matrix(covar,1,NCOVMAX,1,n);*/
Line 4197  Interval (in months) between two waves: Line 4579  Interval (in months) between two waves:
     }      }
   }    }
   
   varprob(optionfilefiname, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);    varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);
   
   fclose(ficrespij);    fclose(ficrespij);
   
     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/    /*if((stepm == 1) && (strcmp(model,".")==0)){*/
   if(prevfcast==1){    if(prevfcast==1){
     if(stepm ==1){      /*    if(stepm ==1){*/
       prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);        prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);        /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
     }   /*      }  */
     else{  /*      else{ */
       erreur=108;  /*        erreur=108; */
       printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);  /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
       fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);  /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
     }  /*      } */
   }    }
       
   
Line 4247  Interval (in months) between two waves: Line 4630  Interval (in months) between two waves:
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   
   prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);    /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
     prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
   ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
     */
   
   if (mobilav!=0) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
Line 4374  Interval (in months) between two waves: Line 4761  Interval (in months) between two waves:
       
   free_matrix(covar,0,NCOVMAX,1,n);    free_matrix(covar,0,NCOVMAX,1,n);
   free_matrix(matcov,1,npar,1,npar);    free_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);    /*free_vector(delti,1,npar);*/
     free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   free_ivector(ncodemax,1,8);    free_ivector(ncodemax,1,8);
   free_ivector(Tvar,1,15);    free_ivector(Tvar,1,15);
   free_ivector(Tprod,1,15);    free_ivector(Tprod,1,15);
Line 4385  Interval (in months) between two waves: Line 4775  Interval (in months) between two waves:
   free_ivector(Tage,1,15);    free_ivector(Tage,1,15);
   free_ivector(Tcode,1,100);    free_ivector(Tcode,1,100);
   
   fprintf(fichtm,"\n</body>");    /*  fclose(fichtm);*/
   fclose(fichtm);    /*  fclose(ficgp);*/ /* ALready done */
   fclose(ficgp);  
       
   
   if(erreur >0){    if(erreur >0){
Line 4400  Interval (in months) between two waves: Line 4789  Interval (in months) between two waves:
   printf("See log file on %s\n",filelog);    printf("See log file on %s\n",filelog);
   fclose(ficlog);    fclose(ficlog);
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
       (void) gettimeofday(&end_time,&tzp);
   /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/    tm = *localtime(&end_time.tv_sec);
   /*printf("Total time was %d uSec.\n", total_usecs);*/    tmg = *gmtime(&end_time.tv_sec);
     strtend=asctime(&tm);
     printf("Localtime at start %s and at end=%s",strt, strtend); 
     fprintf(ficlog,"Localtime at start %s and at end=%s",strt, strtend); 
     /*  printf("Total time used %d Sec\n", asc_time(end_time.tv_sec -start_time.tv_sec);*/
   
     printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
     fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
     /*  printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/    /*------ End -----------*/
   
   end:    end:
Line 4416  Interval (in months) between two waves: Line 4813  Interval (in months) between two waves:
   strcpy(plotcmd,GNUPLOTPROGRAM);    strcpy(plotcmd,GNUPLOTPROGRAM);
   strcat(plotcmd," ");    strcat(plotcmd," ");
   strcat(plotcmd,optionfilegnuplot);    strcat(plotcmd,optionfilegnuplot);
   printf("Starting: %s\n",plotcmd);fflush(stdout);    printf("Starting graphs with: %s",plotcmd);fflush(stdout);
   system(plotcmd);    system(plotcmd);
     printf(" Wait...");
   
  /*#ifdef windows*/   /*#ifdef windows*/
   while (z[0] != 'q') {    while (z[0] != 'q') {
Removed from v.1.73  
changed lines
  Added in v.1.85

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