Diff for /imach/src/imach.c between versions 1.2 and 1.22

version 1.2, 2001/03/13 18:10:26 version 1.22, 2002/02/22 17:54:20
Line 1 Line 1
      /* $Id$
 /*********************** Imach **************************************             Interpolate Markov Chain
   This program computes Healthy Life Expectancies from cross-longitudinal  
   data. Cross-longitudinal consist in a first survey ("cross") where    Short summary of the programme:
   individuals from different ages are interviewed on their health status   
   or degree of  disability. At least a second wave of interviews    This program computes Healthy Life Expectancies from
   ("longitudinal") should  measure each new individual health status.    cross-longitudinal data. Cross-longitudinal data consist in: -1- a
   Health expectancies are computed from the transistions observed between    first survey ("cross") where individuals from different ages are
   waves and are computed for each degree of severity of disability (number    interviewed on their health status or degree of disability (in the
   of life states). More degrees you consider, more time is necessary to    case of a health survey which is our main interest) -2- at least a
   reach the Maximum Likelihood of the parameters involved in the model.    second wave of interviews ("longitudinal") which measure each change
   The simplest model is the multinomial logistic model where pij is    (if any) in individual health status.  Health expectancies are
   the probabibility to be observed in state j at the second wave conditional    computed from the time spent in each health state according to a
   to be observed in state i at the first wave. Therefore the model is:    model. More health states you consider, more time is necessary to reach the
   log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex'    Maximum Likelihood of the parameters involved in the model.  The
   is a covariate. If you want to have a more complex model than "constant and    simplest model is the multinomial logistic model where pij is the
   age", you should modify the program where the markup    probabibility to be observed in state j at the second wave
     *Covariates have to be included here again* invites you to do it.    conditional to be observed in state i at the first wave. Therefore
   More covariates you add, less is the speed of the convergence.    the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
     'age' is age and 'sex' is a covariate. If you want to have a more
   The advantage that this computer programme claims, comes from that if the    complex model than "constant and age", you should modify the program
   delay between waves is not identical for each individual, or if some    where the markup *Covariates have to be included here again* invites
   individual missed an interview, the information is not rounded or lost, but    you to do it.  More covariates you add, slower the
   taken into account using an interpolation or extrapolation.    convergence.
   hPijx is the probability to be  
   observed in state i at age x+h conditional to the observed state i at age    The advantage of this computer programme, compared to a simple
   x. The delay 'h' can be split into an exact number (nh*stepm) of    multinomial logistic model, is clear when the delay between waves is not
   unobserved intermediate  states. This elementary transition (by month or    identical for each individual. Also, if a individual missed an
   quarter trimester, semester or year) is model as a multinomial logistic.    intermediate interview, the information is lost, but taken into
   The hPx matrix is simply the matrix product of nh*stepm elementary matrices    account using an interpolation or extrapolation.  
   and the contribution of each individual to the likelihood is simply hPijx.  
     hPijx is the probability to be observed in state i at age x+h
     conditional to the observed state i at age x. The delay 'h' can be
     split into an exact number (nh*stepm) of unobserved intermediate
     states. This elementary transition (by month or quarter trimester,
     semester or year) is model as a multinomial logistic.  The hPx
     matrix is simply the matrix product of nh*stepm elementary matrices
     and the contribution of each individual to the likelihood is simply
     hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the prevalence limits.    of the life expectancies. It also computes the prevalence limits.
Line 48 Line 56
 #include <unistd.h>  #include <unistd.h>
   
 #define MAXLINE 256  #define MAXLINE 256
   #define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
 #define windows  #define windows
   #define GLOCK_ERROR_NOPATH              -1      /* empty path */
   #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
 #define MAXPARM 30 /* Maximum number of parameters for the optimization */  #define MAXPARM 30 /* Maximum number of parameters for the optimization */
 #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */  #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
Line 59 Line 70
 #define NLSTATEMAX 8 /* Maximum number of live states (for func) */  #define NLSTATEMAX 8 /* Maximum number of live states (for func) */
 #define NDEATHMAX 8 /* Maximum number of dead states (for func) */  #define NDEATHMAX 8 /* Maximum number of dead states (for func) */
 #define NCOVMAX 8 /* Maximum number of covariates */  #define NCOVMAX 8 /* Maximum number of covariates */
 #define MAXN 80000  #define MAXN 20000
 #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
   
   
   int erreur; /* Error number */
 int nvar;  int nvar;
 static int cptcov;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 int cptcovn;  
 int npar=NPARMAX;  int npar=NPARMAX;
 int nlstate=2; /* Number of live states */  int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */  int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */  int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
   int popbased=0;
   
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */  int maxwav; /* Maxim number of waves */
   int jmin, jmax; /* min, max spacing between 2 waves */
 int mle, weightopt;  int mle, weightopt;
 int **mw; /* mw[mi][i] is number of the mi wave for this individual */  int **mw; /* mw[mi][i] is number of the mi wave for this individual */
 int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */  int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
   double jmean; /* Mean space between 2 waves */
 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;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;
 FILE *ficgp, *fichtm;  FILE *ficgp, *fichtm,*ficresprob,*ficpop;
 FILE *ficreseij;  FILE *ficreseij;
   char filerese[FILENAMELENGTH];    char filerese[FILENAMELENGTH];
  FILE  *ficresvij;   FILE  *ficresvij;
Line 89  FILE *ficreseij; Line 103  FILE *ficreseij;
  FILE  *ficresvpl;   FILE  *ficresvpl;
   char fileresvpl[FILENAMELENGTH];    char fileresvpl[FILENAMELENGTH];
   
   
   
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
 #define FTOL 1.0e-10  #define FTOL 1.0e-10
Line 125  int stepm; Line 136  int stepm;
 /* Stepm, step in month: minimum step interpolation*/  /* Stepm, step in month: minimum step interpolation*/
   
 int m,nb;  int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax;  int *num, 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;  double **pmmij, ***probs, ***mobaverage;
   double dateintmean=0;
   
 double *weight;  double *weight;
 int **s; /* Status */  int **s; /* Status */
 double *agedc, **covar, idx;  double *agedc, **covar, idx;
 int **nbcode, *Tcode, *Tvar, **codtab;  int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
   
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */  double ftol=FTOL; /* Tolerance for computing Max Likelihood */
 double ftolhess; /* Tolerance for computing hessian */  double ftolhess; /* Tolerance for computing hessian */
   
   /**************** split *************************/
   static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
   {
      char *s;                             /* pointer */
      int  l1, l2;                         /* length counters */
   
      l1 = strlen( path );                 /* length of path */
      if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   #ifdef windows
      s = strrchr( path, '\\' );           /* find last / */
   #else
      s = strrchr( path, '/' );            /* find last / */
   #endif
      if ( s == NULL ) {                   /* no directory, so use current */
   #if     defined(__bsd__)                /* get current working directory */
         extern char       *getwd( );
   
         if ( getwd( dirc ) == NULL ) {
   #else
         extern char       *getcwd( );
   
         if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
   #endif
            return( GLOCK_ERROR_GETCWD );
         }
         strcpy( name, path );             /* we've got it */
      } else {                             /* strip direcotry from path */
         s++;                              /* after this, the filename */
         l2 = strlen( s );                 /* length of filename */
         if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
         strcpy( name, s );                /* save file name */
         strncpy( dirc, path, l1 - l2 );   /* now the directory */
         dirc[l1-l2] = 0;                  /* add zero */
      }
      l1 = strlen( dirc );                 /* length of directory */
   #ifdef windows
      if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
   #else
      if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
   #endif
      s = strrchr( name, '.' );            /* find last / */
      s++;
      strcpy(ext,s);                       /* save extension */
      l1= strlen( name);
      l2= strlen( s)+1;
      strncpy( finame, name, l1-l2);
      finame[l1-l2]= 0;
      return( 0 );                         /* we're done */
   }
   
   
 /******************************************/  /******************************************/
   
Line 166  int nbocc(char *s, char occ) Line 228  int nbocc(char *s, char occ)
   
 void cutv(char *u,char *v, char*t, char occ)  void cutv(char *u,char *v, char*t, char occ)
 {  {
   int i,lg,j,p;    int i,lg,j,p=0;
   i=0;    i=0;
   if (t[0]== occ) p=0;  
   for(j=0; j<=strlen(t)-1; j++) {    for(j=0; j<=strlen(t)-1; j++) {
     if((t[j]!= occ) && (t[j+1]==occ)) p=j+1;      if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
   }    }
   
   lg=strlen(t);    lg=strlen(t);
   for(j=0; j<p; j++) {    for(j=0; j<p; j++) {
     (u[j] = t[j]);      (u[j] = t[j]);
     u[p]='\0';  
   }    }
        u[p]='\0';
   
    for(j=0; j<= lg; j++) {     for(j=0; j<= lg; j++) {
     if (j>=(p+1))(v[j-p-1] = t[j]);      if (j>=(p+1))(v[j-p-1] = t[j]);
   }    }
 }  }
   
   
 /********************** nrerror ********************/  /********************** nrerror ********************/
   
 void nrerror(char error_text[])  void nrerror(char error_text[])
Line 615  double **prevalim(double **prlim, int nl Line 675  double **prevalim(double **prlim, int nl
     for (j=1;j<=nlstate+ndeath;j++){      for (j=1;j<=nlstate+ndeath;j++){
       oldm[ii][j]=(ii==j ? 1.0 : 0.0);        oldm[ii][j]=(ii==j ? 1.0 : 0.0);
     }      }
   /* Even if hstepm = 1, at least one multiplication by the unit matrix */  
      cov[1]=1.;
    
    /* Even if hstepm = 1, at least one multiplication by the unit matrix */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
     cov[1]=1.;       cov[2]=agefin;
     cov[2]=agefin;   
     if (cptcovn>0){        for (k=1; k<=cptcovn;k++) {
       for (k=1; k<=cptcovn;k++) {cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];/*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/}          cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
     }          /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/
         }
         for (k=1; k<=cptcovage;k++)
           cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
         for (k=1; k<=cptcovprod;k++)
           cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
         /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
         /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
   
     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);      out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
   
     savm=oldm;      savm=oldm;
Line 648  double **prevalim(double **prlim, int nl Line 720  double **prevalim(double **prlim, int nl
   }    }
 }  }
   
 /*************** transition probabilities **********/  /*************** transition probabilities ***************/
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {  {
Line 674  double **pmij(double **ps, double *cov, Line 746  double **pmij(double **ps, double *cov,
       ps[i][j]=s2;        ps[i][j]=s2;
     }      }
   }    }
       /*ps[3][2]=1;*/
   
   for(i=1; i<= nlstate; i++){    for(i=1; i<= nlstate; i++){
      s1=0;       s1=0;
     for(j=1; j<i; j++)      for(j=1; j<i; j++)
Line 695  double **pmij(double **ps, double *cov, Line 769  double **pmij(double **ps, double *cov,
     }      }
   }    }
   
   
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){    /*   for(ii=1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){      for(jj=1; jj<= nlstate+ndeath; jj++){
      printf("%lf ",ps[ii][jj]);       printf("%lf ",ps[ii][jj]);
Line 712  double **pmij(double **ps, double *cov, Line 787  double **pmij(double **ps, double *cov,
   
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)  double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
 {  {
   /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times    /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */       b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
   /* in, b, out are matrice of pointers which should have been initialized    /* in, b, out are matrice of pointers which should have been initialized
      before: only the contents of out is modified. The function returns       before: only the contents of out is modified. The function returns
Line 758  double ***hpxij(double ***po, int nhstep Line 833  double ***hpxij(double ***po, int nhstep
       /* Covariates have to be included here again */        /* Covariates have to be included here again */
       cov[1]=1.;        cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
       if (cptcovn>0){        for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
       for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];        for (k=1; k<=cptcovage;k++)
     }          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
         for (k=1; k<=cptcovprod;k++)
           cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
   
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/        /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
       out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,        out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
Line 782  double ***hpxij(double ***po, int nhstep Line 861  double ***hpxij(double ***po, int nhstep
 /*************** log-likelihood *************/  /*************** log-likelihood *************/
 double func( double *x)  double func( double *x)
 {  {
   int i, ii, j, k, mi, d;    int i, ii, j, k, mi, d, kk;
   double l, ll[NLSTATEMAX], cov[NCOVMAX];    double l, ll[NLSTATEMAX], cov[NCOVMAX];
   double **out;    double **out;
   double sw; /* Sum of weights */    double sw; /* Sum of weights */
Line 792  double func( double *x) Line 871  double func( double *x)
   /* We are differentiating ll according to initial status */    /* We are differentiating ll according to initial status */
   /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/    /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
   /*for(i=1;i<imx;i++)    /*for(i=1;i<imx;i++)
 printf(" %d\n",s[4][i]);      printf(" %d\n",s[4][i]);
   */    */
     cov[1]=1.;
   
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
        for(mi=1; mi<= wav[i]-1; mi++){      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 (ii=1;ii<=nlstate+ndeath;ii++)
         for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);          for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
             for(d=0; d<dh[mi][i]; d++){        for(d=0; d<dh[mi][i]; d++){
         newm=savm;          newm=savm;
           cov[1]=1.;          cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
           cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;          for (kk=1; kk<=cptcovage;kk++) {
           if (cptcovn>0){            cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
             for (k=1; k<=cptcovn;k++) {          }
               cov[2+k]=covar[Tvar[k]][i];         
               /* printf("k=%d cptcovn=%d %lf\n",k,cptcovn,covar[Tvar[k]][i]);*/          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
             }                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             }          savm=oldm;
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,          oldm=newm;
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));         
           savm=oldm;         
           oldm=newm;  
       } /* end mult */        } /* end mult */
           
       lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);        lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
       /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/        /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
       ipmx +=1;        ipmx +=1;
Line 827  printf(" %d\n",s[4][i]); Line 907  printf(" %d\n",s[4][i]);
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */    /* 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 */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
   
   return -l;    return -l;
 }  }
   
Line 847  void mlikeli(FILE *ficres,double p[], in Line 926  void mlikeli(FILE *ficres,double p[], in
   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));     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   
 }  }
   
Line 864  void hesscov(double **matcov, double p[] Line 943  void hesscov(double **matcov, double p[]
   void lubksb(double **a, int npar, int *indx, double b[]) ;    void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;    void ludcmp(double **a, int npar, int *indx, double *d) ;
   
   
   hess=matrix(1,npar,1,npar);    hess=matrix(1,npar,1,npar);
   
   printf("\nCalculation of the hessian matrix. Wait...\n");    printf("\nCalculation of the hessian matrix. Wait...\n");
Line 872  void hesscov(double **matcov, double p[] Line 950  void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);      printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);      hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/      /*printf(" %f ",p[i]);*/
       /*printf(" %lf ",hess[i][i]);*/
   }    }
    
   for (i=1;i<=npar;i++) {    for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {      for (j=1;j<=npar;j++)  {
       if (j>i) {        if (j>i) {
         printf(".%d%d",i,j);fflush(stdout);          printf(".%d%d",i,j);fflush(stdout);
         hess[i][j]=hessij(p,delti,i,j);          hess[i][j]=hessij(p,delti,i,j);
         hess[j][i]=hess[i][j];          hess[j][i]=hess[i][j];    
           /*printf(" %lf ",hess[i][j]);*/
       }        }
     }      }
   }    }
Line 983  double hessii( double x[], double delta, Line 1063  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;      return res;
    
 }  }
   
 double hessij( double x[], double delti[], int thetai,int thetaj)  double hessij( double x[], double delti[], int thetai,int thetaj)
Line 1095  void lubksb(double **a, int n, int *indx Line 1176  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 *Tvar, int **nbcode, int *ncodemax)  void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2)
 {  /* 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;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos;    double pos, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;    FILE *ficresp;
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
     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 1116  void  freqsummary(char fileres[], int ag Line 1197  void  freqsummary(char fileres[], int ag
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
   
   j=cptcovn;    j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
   for(k1=1; k1<=j;k1++){    for(k1=1; k1<=j;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){     for(i1=1; i1<=ncodemax[k1];i1++){
        j1++;         j1++;
          /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
            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=agemin; m <= agemax+3; m++)
              freq[i][jk][m]=0;               freq[i][jk][m]=0;
          
           dateintsum=0;
           k2cpt=0;
        for (i=1; i<=imx; i++) {         for (i=1; i<=imx; i++) {
          bool=1;           bool=1;
          if  (cptcovn>0) {           if  (cptcovn>0) {
            for (z1=1; z1<=cptcovn; z1++)             for (z1=1; z1<=cptcoveff; z1++)
              if (covar[Tvar[z1]][i]!= nbcode[Tvar[z1]][codtab[j1][z1]]) bool=0;               if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
                  bool=0;
          }           }
           if (bool==1) {           if (bool==1) {
            for(m=firstpass; m<=lastpass-1; m++){             for(m=firstpass; m<=lastpass; m++){
              if(agev[m][i]==0) agev[m][i]=agemax+1;               k2=anint[m][i]+(mint[m][i]/12.);
              if(agev[m][i]==1) agev[m][i]=agemax+2;               if ((k2>=dateprev1) && (k2<=dateprev2)) {
              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                 if(agev[m][i]==0) agev[m][i]=agemax+1;
              freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                 if(agev[m][i]==1) agev[m][i]=agemax+2;
                  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                  freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
                  if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
                    dateintsum=dateintsum+k2;
                    k2cpt++;
                  }
   
                }
            }             }
          }           }
        }         }
         if  (cptcovn>0) {          if  (cptcovn>0) {
          fprintf(ficresp, "\n#Variable");           fprintf(ficresp, "\n#********** Variable ");
          for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvar[z1],nbcode[Tvar[z1]][codtab[j1][z1]]);           for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
        }         fprintf(ficresp, "**********\n#");
        fprintf(ficresp, "\n#");          }
        for(i=1; i<=nlstate;i++)         for(i=1; i<=nlstate;i++)
          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");
Line 1159  void  freqsummary(char fileres[], int ag Line 1252  void  freqsummary(char fileres[], int ag
       printf("Age %d", i);        printf("Age %d", i);
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)        for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }      }
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pos=0; m <=0 ; m++)        for(m=-1, pos=0; m <=0 ; m++)
Line 1169  void  freqsummary(char fileres[], int ag Line 1262  void  freqsummary(char fileres[], int ag
       else        else
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);          printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     }      }
     for(jk=1; jk <=nlstate ; jk++){  
       for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)       for(jk=1; jk <=nlstate ; jk++){
         for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }       }
   
     for(jk=1,pos=0; jk <=nlstate ; jk++)      for(jk=1,pos=0; jk <=nlstate ; jk++)
       pos += pp[jk];        pos += pp[jk];
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
Line 1181  void  freqsummary(char fileres[], int ag Line 1276  void  freqsummary(char fileres[], int ag
       else        else
         printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);          printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
       if( i <= (int) agemax){        if( i <= (int) agemax){
         if(pos>=1.e-5)          if(pos>=1.e-5){
           fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);            fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, 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]);*/
           }
       else        else
           fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);            fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
       }        }
Line 1196  void  freqsummary(char fileres[], int ag Line 1294  void  freqsummary(char fileres[], int ag
     }      }
     }      }
  }   }
     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,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   
     /* End of Freq */
   }
   
   /************ Prevalence ********************/
   void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate)
   {  /* Some frequencies */
    
     int i, m, jk, k1, i1, j1, bool, z1,z2,j;
     double ***freq; /* Frequencies */
     double *pp;
     double pos, k2;
   
     pp=vector(1,nlstate);
     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
    
     freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
     j1=0;
    
     j=cptcoveff;
     if (cptcovn<1) {j=1;ncodemax[1]=1;}
    
    for(k1=1; k1<=j;k1++){
       for(i1=1; i1<=ncodemax[k1];i1++){
         j1++;
    
         for (i=-1; i<=nlstate+ndeath; i++)  
           for (jk=-1; jk<=nlstate+ndeath; jk++)  
             for(m=agemin; m <= agemax+3; m++)
               freq[i][jk][m]=0;
        
         for (i=1; i<=imx; i++) {
           bool=1;
           if  (cptcovn>0) {
             for (z1=1; z1<=cptcoveff; z1++)
               if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
                 bool=0;
           }
           if (bool==1) {
             for(m=firstpass; m<=lastpass; m++){
               k2=anint[m][i]+(mint[m][i]/12.);
               if ((k2>=dateprev1) && (k2<=dateprev2)) {
                 if(agev[m][i]==0) agev[m][i]=agemax+1;
                 if(agev[m][i]==1) agev[m][i]=agemax+2;
                 freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
                 freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  
               }
             }
           }
         }
        
           for(i=(int)agemin; i <= (int)agemax+3; i++){
             for(jk=1; jk <=nlstate ; jk++){
               for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
                 pp[jk] += freq[jk][m][i];
             }
             for(jk=1; jk <=nlstate ; jk++){
               for(m=-1, pos=0; m <=0 ; m++)
               pos += freq[jk][m][i];
           }
          
            for(jk=1; jk <=nlstate ; jk++){
              for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
                pp[jk] += freq[jk][m][i];
            }
            
            for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
   
            for(jk=1; jk <=nlstate ; jk++){          
              if( i <= (int) agemax){
                if(pos>=1.e-5){
                  probs[i][jk][j1]= pp[jk]/pos;
                }
              }
            }
            
           }
       }
     }
    
    
     free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
     free_vector(pp,1,nlstate);
    
 }  /* End of Freq */  }  /* End of Freq */
   
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
Line 1215  void  concatwav(int wav[], int **dh, int Line 1397  void  concatwav(int wav[], int **dh, int
      */       */
   
   int i, mi, m;    int i, mi, m;
   int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
 float sum=0.;       double sum=0., jmean=0.;*/
   
     int j, k=0,jk, ju, jl;
     double sum=0.;
     jmin=1e+5;
     jmax=-1;
     jmean=0.;
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
Line 1247  float sum=0.; Line 1434  float sum=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) {
             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;
             if (j >= jmax) jmax=j;
             if (j <= jmin) jmin=j;
             sum=sum+j;
             /* if (j<10) printf("j=%d num=%d ",j,i); */
             }
         }          }
         else{          else{
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));            j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
           /*printf("i=%d agevi+1=%lf agevi=%lf j=%d\n", i,agev[mw[mi+1][i]][i],agev[mw[mi][i]][i],j);*/  
   
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           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); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1271  float sum=0.; Line 1464  float sum=0.;
       }        }
     }      }
   }    }
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k);    jmean=sum/k;
 }    printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
    }
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
   int Ndum[80],ij, k, j, i;    int Ndum[20],ij=1, k, j, i;
   int cptcode=0;    int cptcode=0;
   for (k=0; k<79; k++) Ndum[k]=0;    cptcoveff=0;
    
     for (k=0; k<19; k++) Ndum[k]=0;
   for (k=1; k<=7; k++) ncodemax[k]=0;    for (k=1; k<=7; k++) ncodemax[k]=0;
    
   for (j=1; j<=cptcovn; j++) {    for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
     for (i=1; i<=imx; i++) {      for (i=1; i<=imx; i++) {
       ij=(int)(covar[Tvar[j]][i]);        ij=(int)(covar[Tvar[j]][i]);
       Ndum[ij]++;        Ndum[ij]++;
         /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
       if (ij > cptcode) cptcode=ij;        if (ij > cptcode) cptcode=ij;
     }      }
     /*printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);*/  
     for (i=0; i<=cptcode; i++) {      for (i=0; i<=cptcode; i++) {
       if(Ndum[i]!=0) ncodemax[j]++;        if(Ndum[i]!=0) ncodemax[j]++;
     }      }
    
     ij=1;      ij=1;
   
   
     for (i=1; i<=ncodemax[j]; i++) {      for (i=1; i<=ncodemax[j]; i++) {
       for (k=0; k<=79; k++) {        for (k=0; k<=19; k++) {
         if (Ndum[k] != 0) {          if (Ndum[k] != 0) {
           nbcode[Tvar[j]][ij]=k;            nbcode[Tvar[j]][ij]=k;
           ij++;            ij++;
Line 1302  void tricode(int *Tvar, int **nbcode, in Line 1500  void tricode(int *Tvar, int **nbcode, in
         if (ij > ncodemax[j]) break;          if (ij > ncodemax[j]) break;
       }          }  
     }      }
   }      }  
   
   }   for (k=0; k<19; k++) Ndum[k]=0;
   
    for (i=1; i<=ncovmodel-2; i++) {
         ij=Tvar[i];
         Ndum[ij]++;
       }
   
    ij=1;
    for (i=1; i<=10; i++) {
      if((Ndum[i]!=0) && (i<=ncov)){
        Tvaraff[ij]=i;
        ij++;
      }
    }
    
       cptcoveff=ij-1;
   }
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
Line 1365  void varevsij(char fileres[], double *** Line 1579  void varevsij(char fileres[], double ***
   double **dnewm,**doldm;    double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h;
   int k, cptcode;    int k, cptcode;
    double *xp;    double *xp;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   double ***p3mat;    double ***p3mat;
Line 1401  void varevsij(char fileres[], double *** Line 1615  void varevsij(char fileres[], double ***
       }        }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
   
         if (popbased==1) {
           for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];
         }
        
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
Line 1412  void varevsij(char fileres[], double *** Line 1632  void varevsij(char fileres[], double ***
         xp[i] = x[i] - (i==theta ?delti[theta]:0);          xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
   
         if (popbased==1) {
           for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];
         }
   
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gm[h][j]=0.;i<=nlstate;i++)            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
             gm[h][j] += prlim[i][i]*p3mat[i][j][h];              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
         }          }
       }        }
   
       for(j=1; j<= nlstate; j++)        for(j=1; j<= nlstate; j++)
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; 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 1547  void varprevlim(char fileres[], double * Line 1774  void varprevlim(char fileres[], double *
   
 }  }
   
   /************ Variance of one-step probabilities  ******************/
   void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
   {
     int i, j;
     int k=0, cptcode;
     double **dnewm,**doldm;
     double *xp;
     double *gp, *gm;
     double **gradg, **trgradg;
     double age,agelim, cov[NCOVMAX];
     int theta;
     char fileresprob[FILENAMELENGTH];
   
     strcpy(fileresprob,"prob");
     strcat(fileresprob,fileres);
     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprob);
     }
     printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);
    
   
     xp=vector(1,npar);
     dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
    
     cov[1]=1;
     for (age=bage; age<=fage; age ++){
       cov[2]=age;
       gradg=matrix(1,npar,1,9);
       trgradg=matrix(1,9,1,npar);
       gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
       gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
      
       for(theta=1; theta <=npar; theta++){
         for(i=1; i<=npar; i++)
           xp[i] = x[i] + (i==theta ?delti[theta]:0);
        
         pmij(pmmij,cov,ncovmodel,xp,nlstate);
      
         k=0;
         for(i=1; i<= (nlstate+ndeath); i++){
           for(j=1; j<=(nlstate+ndeath);j++){
              k=k+1;
             gp[k]=pmmij[i][j];
           }
         }
   
         for(i=1; i<=npar; i++)
           xp[i] = x[i] - (i==theta ?delti[theta]:0);
      
   
         pmij(pmmij,cov,ncovmodel,xp,nlstate);
         k=0;
         for(i=1; i<=(nlstate+ndeath); i++){
           for(j=1; j<=(nlstate+ndeath);j++){
             k=k+1;
             gm[k]=pmmij[i][j];
           }
         }
        
          for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
              gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
       }
   
        for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
         for(theta=1; theta <=npar; theta++)
         trgradg[j][theta]=gradg[theta][j];
    
        matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
        matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
   
        pmij(pmmij,cov,ncovmodel,x,nlstate);
   
        k=0;
        for(i=1; i<=(nlstate+ndeath); i++){
          for(j=1; j<=(nlstate+ndeath);j++){
            k=k+1;
            gm[k]=pmmij[i][j];
           }
        }
        
        /*printf("\n%d ",(int)age);
        for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
          
   
          printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
        }*/
   
     fprintf(ficresprob,"\n%d ",(int)age);
   
     for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
       if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
   if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
     }
   
       free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
       free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
       free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
       free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   }
    free_vector(xp,1,npar);
   fclose(ficresprob);
    exit(0);
   }
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
 /***********************************************/  /***********************************************/
   
 /*int main(int argc, char *argv[])*/  int main(int argc, char *argv[])
 int main()  
 {  {
   
   int i,j, k, n=MAXN,iter,m,size,cptcode, aaa, cptcod;    int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
   double agemin=1.e20, agemax=-1.e20;    double agemin=1.e20, agemax=-1.e20;
   
Line 1569  int main() Line 1898  int main()
   int *indx;    int *indx;
   char line[MAXLINE], linepar[MAXLINE];    char line[MAXLINE], linepar[MAXLINE];
   char title[MAXLINE];    char title[MAXLINE];
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];    char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];    char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
    
     char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];;
   
   char filerest[FILENAMELENGTH];    char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];    char fileregp[FILENAMELENGTH];
     char popfile[FILENAMELENGTH];
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];    char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
   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;
   int ju,jl, mi;    int ju,jl, mi;
   int i1,j1, k1,jk,aa,bb, stepsize;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
      int mobilav=0,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
     int *popage;/*boolprev=0 if date and zero if wave*/
     double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
   double **prlim;    double **prlim;
Line 1594  int main() Line 1930  int main()
   double ***eij, ***vareij;    double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */    double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;    double *epj, vepp;
   char version[80]="Imach version 62c, May 1999, INED-EUROREVES ";    double kk1, kk2;
     double *popeffectif,*popcount;
     double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
     double yp,yp1,yp2;
   
     char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
   char z[1]="c", occ;    char z[1]="c", occ;
 #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];
    
   /* long total_usecs;    /* long total_usecs;
   struct timeval start_time, end_time;    struct timeval start_time, end_time;
     
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   
   
   printf("\nIMACH, Version 0.63");    printf("\n%s",version);
   printf("\nEnter the parameter file name: ");    if(argc <=1){
       printf("\nEnter the parameter file name: ");
       scanf("%s",pathtot);
     }
     else{
       strcpy(pathtot,argv[1]);
     }
     /*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/
     /*cygwin_split_path(pathtot,path,optionfile);
       printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
     /* cutv(path,optionfile,pathtot,'\\');*/
   
 #ifdef windows    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   scanf("%s",pathtot);     printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   getcwd(pathcd, size);  
   cutv(path,optionfile,pathtot,'\\');  
   chdir(path);    chdir(path);
   replace(pathc,path);    replace(pathc,path);
 #endif  
 #ifdef unix  
   scanf("%s",optionfile);  
 #endif  
   
 /*-------- arguments in the command line --------*/  /*-------- arguments in the command line --------*/
   
   strcpy(fileres,"r");    strcpy(fileres,"r");
   strcat(fileres, optionfile);    strcat(fileres, optionfilefiname);
     strcat(fileres,".txt");    /* Other files have txt extension */
   
   /*---------arguments file --------*/    /*---------arguments file --------*/
   
Line 1650  int main() Line 1999  int main()
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);    fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);    fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
   while((c=getc(ficpar))=='#' && c!= EOF){
   covar=matrix(1,NCOVMAX,1,n);          ungetc(c,ficpar);
   if (strlen(model)<=1) cptcovn=0;      fgets(line, MAXLINE, ficpar);
   else {      puts(line);
     j=0;      fputs(line,ficparo);
     j=nbocc(model,'+');  
     cptcovn=j+1;  
   }    }
     ungetc(c,ficpar);
    
      
     covar=matrix(0,NCOVMAX,1,n);
     cptcovn=0;
     if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
   
   ncovmodel=2+cptcovn;    ncovmodel=2+cptcovn;
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */    nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
Line 1688  int main() Line 2041  int main()
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
     
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel;      npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
   
   p=param[1][1];    p=param[1][1];
     
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
Line 1749  int main() Line 2103  int main()
   printf("\n");    printf("\n");
   
   
    if(mle==1){  
     /*-------- data file ----------*/      /*-------- data file ----------*/
     if((ficres =fopen(fileres,"w"))==NULL) {      if((ficres =fopen(fileres,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileres);goto end;        printf("Problem with resultfile: %s\n", fileres);goto end;
Line 1777  int main() Line 2130  int main()
     s=imatrix(1,maxwav+1,1,n);      s=imatrix(1,maxwav+1,1,n);
     adl=imatrix(1,maxwav+1,1,n);          adl=imatrix(1,maxwav+1,1,n);    
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,NCOVMAX);      ncodemax=ivector(1,8);
   
     i=1;      i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {      while (fgets(line, MAXLINE, fic) != NULL)    {
       if ((i >= firstobs) && (i <=lastobs)) {        if ((i >= firstobs) && (i <=lastobs)) {
                 
Line 1801  int main() Line 2154  int main()
           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);          num[i]=atol(stra);
          
         /* 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]));*/          /*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("%d %.lf %.lf %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),(covar[3][i]), (covar[4][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]));*/  
   
         i=i+1;          i=i+1;
       }        }
     }      }
     /*scanf("%d",i);*/      /* printf("ii=%d", ij);
          scanf("%d",i);*/
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
    
     /* for (i=1; i<=imx; i++){
       if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
       if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
       if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
       }
   
       for (i=1; i<=imx; i++)
       if (covar[1][i]==0) 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]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,8);        Tvar=ivector(1,15);
     Tprod=ivector(1,15);
     Tvaraff=ivector(1,15);
     Tvard=imatrix(1,15,1,2);
     Tage=ivector(1,15);      
         
   if (strlen(model) >1){    if (strlen(model) >1){
     j=0;      j=0, j1=0, k1=1, k2=1;
     j=nbocc(model,'+');      j=nbocc(model,'+');
       j1=nbocc(model,'*');
     cptcovn=j+1;      cptcovn=j+1;
        cptcovprod=j1;
      
      
     strcpy(modelsav,model);      strcpy(modelsav,model);
     if (j==0) {      if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
       cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb);        printf("Error. Non available option model=%s ",model);
         goto end;
     }      }
     else {     
       for(i=j; i>=1;i--){      for(i=(j+1); i>=1;i--){
         cutv(stra,strb,modelsav,'+');        cutv(stra,strb,modelsav,'+');
         if (strchr(strb,'*')) {        if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);
           cutv(strd,strc,strb,'*');        /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
           cutv(strb,stre,strc,'V');Tvar[i+1]=ncov+1;        /*scanf("%d",i);*/
           cutv(strb,strc,strd,'V');        if (strchr(strb,'*')) {
           for (k=1; k<=lastobs;k++)          cutv(strd,strc,strb,'*');
             covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];          if (strcmp(strc,"age")==0) {
             cptcovprod--;
             cutv(strb,stre,strd,'V');
             Tvar[i]=atoi(stre);
             cptcovage++;
               Tage[cptcovage]=i;
               /*printf("stre=%s ", stre);*/
           }
           else if (strcmp(strd,"age")==0) {
             cptcovprod--;
             cutv(strb,stre,strc,'V');
             Tvar[i]=atoi(stre);
             cptcovage++;
             Tage[cptcovage]=i;
         }          }
         else {          else {
           cutv(strd,strc,strb,'V');            cutv(strb,stre,strc,'V');
           Tvar[i+1]=atoi(strc);            Tvar[i]=ncov+k1;
             cutv(strb,strc,strd,'V');
             Tprod[k1]=i;
             Tvard[k1][1]=atoi(strc);
             Tvard[k1][2]=atoi(stre);
             Tvar[cptcovn+k2]=Tvard[k1][1];
             Tvar[cptcovn+k2+1]=Tvard[k1][2];
             for (k=1; k<=lastobs;k++)
               covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
             k1++;
             k2=k2+2;
         }          }
         strcpy(modelsav,stra);    
       }        }
       /*cutv(strd,strc,stra,'V');*/        else {
       Tvar[1]=atoi(strc);          /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
          /*  scanf("%d",i);*/
         cutv(strd,strc,strb,'V');
         Tvar[i]=atoi(strc);
         }
         strcpy(modelsav,stra);  
         /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
           scanf("%d",i);*/
     }      }
   }  }
   /*printf("tvar=%d ",Tvar[1]);*/   
   /*scanf("%d ",i);*/    /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
     printf("cptcovprod=%d ", cptcovprod);
     scanf("%d ",i);*/
     fclose(fic);      fclose(fic);
   
       /*  if(mle==1){*/
     if (weightopt != 1) { /* Maximisation without weights*/      if (weightopt != 1) { /* Maximisation without weights*/
       for(i=1;i<=n;i++) weight[i]=1.0;        for(i=1;i<=n;i++) weight[i]=1.0;
     }      }
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);      agev=matrix(1,maxwav,1,imx);
   
      for (i=1; i<=imx; i++)
        for(m=2; (m<= maxwav); m++)
          if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
            anint[m][i]=9999;
            s[m][i]=-1;
          }
         
     for (i=1; i<=imx; i++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
       for(m=1; (m<= maxwav); m++){        for(m=1; (m<= maxwav); m++){
         if (mint[m][i]==99 || anint[m][i]==9999) s[m][i]=-1;    
         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(moisdc[i]!=99 && andc[i]!=9999)
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
             else{              else {
                 if (andc[i]!=9999){
               printf("Warning negative age at death: %d line:%d\n",num[i],i);                printf("Warning negative age at death: %d line:%d\n",num[i],i);
               agev[m][i]=-1;                agev[m][i]=-1;
                 }
             }              }
           }            }
           else if(s[m][i] !=9){ /* Should no more exist */            else if(s[m][i] !=9){ /* Should no more exist */
             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(mint[m][i]==99 || anint[m][i]==9999)
               agev[m][i]=1;                agev[m][i]=1;
               /* printf("i=%d m=%d agev=%lf \n",i,m, agev[m][i]);    */  
             }  
             else if(agev[m][i] <agemin){              else if(agev[m][i] <agemin){
               agemin=agev[m][i];                agemin=agev[m][i];
               /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/                /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
Line 1911  printf("Total number of individuals= %d, Line 2318  printf("Total number of individuals= %d,
     free_imatrix(outcome,1,maxwav+1,1,n);      free_imatrix(outcome,1,maxwav+1,1,n);
     free_vector(moisnais,1,n);      free_vector(moisnais,1,n);
     free_vector(annais,1,n);      free_vector(annais,1,n);
     free_matrix(mint,1,maxwav,1,n);      /* free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n);         free_matrix(anint,1,maxwav,1,n);*/
     free_vector(moisdc,1,n);      free_vector(moisdc,1,n);
     free_vector(andc,1,n);      free_vector(andc,1,n);
   
Line 1925  printf("Total number of individuals= %d, Line 2332  printf("Total number of individuals= %d,
       concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);        concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   
   
 Tcode=ivector(1,100);        Tcode=ivector(1,100);
    nbcode=imatrix(1,nvar,1,8);          nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
    ncodemax[1]=1;        ncodemax[1]=1;
    if (cptcovn > 0) tricode(Tvar,nbcode,imx);        if (cptcovn > 0) tricode(Tvar,nbcode,imx);
         
    codtab=imatrix(1,100,1,10);     codtab=imatrix(1,100,1,10);
    h=0;     h=0;
    m=pow(2,cptcovn);     m=pow(2,cptcoveff);
     
    for(k=1;k<=cptcovn; k++){     for(k=1;k<=cptcoveff; k++){
      for(i=1; i <=(m/pow(2,k));i++){       for(i=1; i <=(m/pow(2,k));i++){
        for(j=1; j <= ncodemax[k]; j++){         for(j=1; j <= ncodemax[k]; j++){
          for(cpt=1; cpt <=(m/pow(2,cptcovn+1-k)); cpt++){           for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
            h++;             h++;
            if (h>m) h=1;codtab[h][k]=j;             if (h>m) h=1;codtab[h][k]=j;
          }           }
Line 1945  Tcode=ivector(1,100); Line 2352  Tcode=ivector(1,100);
      }       }
    }     }
   
   
    /*for(i=1; i <=m ;i++){     /*for(i=1; i <=m ;i++){
      for(k=1; k <=cptcovn; k++){       for(k=1; k <=cptcovn; k++){
        printf("i=%d k=%d %d ",i,k,codtab[i][k]);         printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
      }       }
      printf("\n");       printf("\n");
    }     }
   scanf("%d",i);*/     scanf("%d",i);*/
         
    /* 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,Tvar,nbcode, ncodemax);  
   
      
      
     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 */
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */      oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
           
     /* For Powell, parameters are in a vector p[] starting at p[1]      /* For Powell, parameters are in a vector p[] starting at p[1]
        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) */
     /*scanf("%d",i);*/  
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);  
   
       if(mle==1){
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
       }
         
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
       
   
    jk=1;     jk=1;
    fprintf(ficres,"# Parameters\n");     fprintf(ficres,"# Parameters\n");
    printf("# Parameters\n");     printf("# Parameters\n");
Line 1992  Tcode=ivector(1,100); Line 2403  Tcode=ivector(1,100);
          }           }
      }       }
    }     }
    if(mle==1){
     /* 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);
    }
     fprintf(ficres,"# Scales\n");      fprintf(ficres,"# Scales\n");
     printf("# Scales\n");      printf("# Scales\n");
      for(i=1,jk=1; i <=nlstate; i++){       for(i=1,jk=1; i <=nlstate; i++){
Line 2012  Tcode=ivector(1,100); Line 2424  Tcode=ivector(1,100);
           fprintf(ficres,"\n");            fprintf(ficres,"\n");
         }          }
       }        }
       }       }
         
     k=1;      k=1;
     fprintf(ficres,"# Covariance\n");      fprintf(ficres,"# Covariance\n");
Line 2047  Tcode=ivector(1,100); Line 2459  Tcode=ivector(1,100);
       bage = agemin;        bage = agemin;
       fage = agemax;        fage = agemax;
     }      }
      
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);      fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 /*------------ gnuplot -------------*/      fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 chdir(pathcd);   
   if((ficgp=fopen("graph.plt","w"))==NULL) {      while((c=getc(ficpar))=='#' && c!= EOF){
     printf("Problem with file graph.gp");goto end;      ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
   }    }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
     fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
        
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
   
      dateprev1=anprev1+mprev1/12.+jprev1/365.;
      dateprev2=anprev2+mprev2/12.+jprev2/365.;
   
     fscanf(ficpar,"pop_based=%d\n",&popbased);
      fprintf(ficparo,"pop_based=%d\n",popbased);  
      fprintf(ficres,"pop_based=%d\n",popbased);  
   
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
     fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
   fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
   fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
   
    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
   
      
       /*------------ gnuplot -------------*/
       /*chdir(pathcd);*/
       strcpy(optionfilegnuplot,optionfilefiname);
       strcat(optionfilegnuplot,".plt");
       if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
         printf("Problem with file %s",optionfilegnuplot);goto end;
       }
 #ifdef windows  #ifdef windows
   fprintf(ficgp,"cd \"%s\" \n",pathc);      fprintf(ficgp,"cd \"%s\" \n",pathc);
 #endif  #endif
 m=pow(2,cptcovn);  m=pow(2,cptcoveff);
     
  /* 1eme*/   /* 1eme*/
   for (cpt=1; cpt<= nlstate ; cpt ++) {    for (cpt=1; cpt<= nlstate ; cpt ++) {
Line 2154  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2612  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);          fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);  
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);        fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
     }      }
   }    }  
   
   /* proba elementaires */    /* proba elementaires */
   for(i=1,jk=1; i <=nlstate; i++){     for(i=1,jk=1; i <=nlstate; i++){
     for(k=1; k <=(nlstate+ndeath); k++){      for(k=1; k <=(nlstate+ndeath); k++){
       if (k != i) {        if (k != i) {
         /*  fprintf(ficgp,"%1d%1d ",i,k);*/  
         for(j=1; j <=ncovmodel; j++){          for(j=1; j <=ncovmodel; j++){
           fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);            /*fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/
             /*fprintf(ficgp,"%s",alph[1]);*/
             fprintf(ficgp,"p%d=%f ",jk,p[jk]);
           jk++;            jk++;
           fprintf(ficgp,"\n");            fprintf(ficgp,"\n");
         }          }
       }        }
     }      }
   }      }
   
   for(jk=1; jk <=m; jk++) {    for(jk=1; jk <=m; jk++) {
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);    fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);
   for(i=1; i <=nlstate; i++) {     i=1;
     for(k=1; k <=(nlstate+ndeath); k++){     for(k2=1; k2<=nlstate; k2++) {
       if (k != i) {       k3=i;
         fprintf(ficgp," exp(a%d%d+b%d%d*x",i,k,i,k);       for(k=1; k<=(nlstate+ndeath); k++) {
         for(j=3; j <=ncovmodel; j++)         if (k != k2){
           fprintf(ficgp,"+%s%d%d*%d",alph[j],i,k,nbcode[Tvar[j-2]][codtab[jk][j-2]]);          fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
         fprintf(ficgp,")/(1");  ij=1;
         for(k1=1; k1 <=(nlstate+ndeath); k1++)          for(j=3; j <=ncovmodel; j++) {
           if (k1 != i) {            if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
             fprintf(ficgp,"+exp(a%d%d+b%d%d*x",i,k1,i,k1);              fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
             for(j=3; j <=ncovmodel; j++)              ij++;
               fprintf(ficgp,"+%s%d%d*%d",alph[j],i,k,nbcode[Tvar[j-2]][codtab[jk][j-2]]);  
             fprintf(ficgp,")");  
           }            }
         fprintf(ficgp,") t \"p%d%d\" ", i,k);            else
       if ((i+k)!= (nlstate*2+ndeath)) fprintf(ficgp,",");            fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
       }          }
     }            fprintf(ficgp,")/(1");
   }         
 fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);            for(k1=1; k1 <=nlstate; k1++){  
             fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
   ij=1;
             for(j=3; j <=ncovmodel; j++){
             if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
               fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
               ij++;
             }
             else
               fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
             }
             fprintf(ficgp,")");
           }
           fprintf(ficgp,") t \"p%d%d\" ", k2,k);
           if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
           i=i+ncovmodel;
          }
        }
      }
      fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
   }    }
      
   fclose(ficgp);    fclose(ficgp);
     /* end gnuplot */
         
 chdir(path);  chdir(path);
     free_matrix(agev,1,maxwav,1,imx);     
     free_ivector(wav,1,imx);      free_ivector(wav,1,imx);
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);      free_imatrix(dh,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_imatrix(s,1,maxwav+1,1,n);  
      
      
     free_ivector(num,1,n);      free_ivector(num,1,n);
     free_vector(agedc,1,n);      free_vector(agedc,1,n);
     free_vector(weight,1,n);  
     /*free_matrix(covar,1,NCOVMAX,1,n);*/      /*free_matrix(covar,1,NCOVMAX,1,n);*/
     fclose(ficparo);      fclose(ficparo);
     fclose(ficres);      fclose(ficres);
    }      /*  }*/
         
    /*________fin mle=1_________*/     /*________fin mle=1_________*/
         
Line 2230  chdir(path); Line 2704  chdir(path);
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);    fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
   if((fichtm=fopen("index.htm","w"))==NULL)    {    strcpy(optionfilehtm,optionfile);
     printf("Problem with index.htm \n");goto end;    strcat(optionfilehtm,".htm");
     if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
       printf("Problem with %s \n",optionfilehtm);goto end;
   }    }
   
  fprintf(fichtm,"<body><ul> Imach, Version 0.63<hr> <li>Outputs files<br><br>\n   fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">
   Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
   Total number of observations=%d <br>
   Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
   <hr  size=\"2\" color=\"#EC5E5E\">
   <li>Outputs files<br><br>\n
         - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n          - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
 - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>  - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
         - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>          - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
Line 2243  chdir(path); Line 2724  chdir(path);
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>          - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>          - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>          - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);          - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
           - Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>
   <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
   
  fprintf(fichtm," <li>Graphs</li>\n<p>");   fprintf(fichtm," <li>Graphs</li><p>");
   
  m=cptcovn;   m=cptcoveff;
  if (cptcovn < 1) {m=1;ncodemax[1]=1;}   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   
  j1=0;   j1=0;
Line 2255  chdir(path); Line 2738  chdir(path);
    for(i1=1; i1<=ncodemax[k1];i1++){     for(i1=1; i1<=ncodemax[k1];i1++){
        j1++;         j1++;
        if (cptcovn > 0) {         if (cptcovn > 0) {
          fprintf(fichtm,"<hr>************ Results for covariates");           fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
          for (cpt=1; cpt<=cptcovn;cpt++)           for (cpt=1; cpt<=cptcoveff;cpt++)
            fprintf(fichtm," V%d=%d ",Tvar[cpt],nbcode[Tvar[cpt]][codtab[j1][cpt]]);             fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);
          fprintf(fichtm," ************\n<hr>");           fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        }         }
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>         fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);      <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);    
Line 2306  fclose(fichtm); Line 2789  fclose(fichtm);
   agebase=agemin;    agebase=agemin;
   agelim=agemax;    agelim=agemax;
   ftolpl=1.e-10;    ftolpl=1.e-10;
   i1=cptcovn;    i1=cptcoveff;
   if (cptcovn < 1){i1=1;}    if (cptcovn < 1){i1=1;}
   
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
         k=k+1;          k=k+1;
         /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/          /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
         fprintf(ficrespl,"\n#****** ");          fprintf(ficrespl,"\n#******");
         for(j=1;j<=cptcovn;j++)          for(j=1;j<=cptcoveff;j++)
           fprintf(ficrespl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);            fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrespl,"******\n");          fprintf(ficrespl,"******\n");
                 
         for (age=agebase; age<=agelim; age++){          for (age=agebase; age<=agelim; age++){
Line 2328  fclose(fichtm); Line 2811  fclose(fichtm);
       }        }
     }      }
   fclose(ficrespl);    fclose(ficrespl);
   
   /*------------- h Pij x at various ages ------------*/    /*------------- h Pij x at various ages ------------*/
     
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);    strcpy(filerespij,"pij");  strcat(filerespij,fileres);
Line 2337  fclose(fichtm); Line 2821  fclose(fichtm);
   printf("Computing pij: result on file '%s' \n", filerespij);    printf("Computing pij: result on file '%s' \n", filerespij);
     
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
   if (stepm<=24) stepsize=2;    /*if (stepm<=24) stepsize=2;*/
   
   agelim=AGESUP;    agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */    hstepm=stepsize*YEARM; /* Every year of age */
Line 2348  fclose(fichtm); Line 2832  fclose(fichtm);
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;        k=k+1;
         fprintf(ficrespij,"\n#****** ");          fprintf(ficrespij,"\n#****** ");
         for(j=1;j<=cptcovn;j++)          for(j=1;j<=cptcoveff;j++)
           fprintf(ficrespij,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);            fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrespij,"******\n");          fprintf(ficrespij,"******\n");
                 
         for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */          for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
Line 2376  fclose(fichtm); Line 2860  fclose(fichtm);
     }      }
   }    }
   
     /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
   
   fclose(ficrespij);    fclose(ficrespij);
   
     if(stepm == 1) {
     /*---------- Forecasting ------------------*/
     calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
   
     /*printf("calage= %f", calagedate);*/
    
     prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
   
     strcpy(fileresf,"f");
     strcat(fileresf,fileres);
     if((ficresf=fopen(fileresf,"w"))==NULL) {
       printf("Problem with forecast resultfile: %s\n", fileresf);goto end;
     }
     printf("Computing forecasting: result on file '%s' \n", fileresf);
   
     free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n);
     free_matrix(agev,1,maxwav,1,imx);
     /* Mobile average */
   
     if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
     if (mobilav==1) {
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
       for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
         for (i=1; i<=nlstate;i++)
           for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
             mobaverage[(int)agedeb][i][cptcod]=0.;
      
       for (agedeb=bage+4; agedeb<=fage; agedeb++){
         for (i=1; i<=nlstate;i++){
           for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
             for (cpt=0;cpt<=4;cpt++){
               mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
             }
             mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
           }
         }
       }  
     }
   
     stepsize=(int) (stepm+YEARM-1)/YEARM;
     if (stepm<=12) stepsize=1;
   
     agelim=AGESUP;
     /*hstepm=stepsize*YEARM; *//* Every year of age */
     hstepm=1;
     hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
     yp1=modf(dateintmean,&yp);
     anprojmean=yp;
     yp2=modf((yp1*12),&yp);
     mprojmean=yp;
     yp1=modf((yp2*30.5),&yp);
     jprojmean=yp;
     if(jprojmean==0) jprojmean=1;
     if(mprojmean==0) jprojmean=1;
   
     fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);
   
     if (popforecast==1) {
       if((ficpop=fopen(popfile,"r"))==NULL)    {
         printf("Problem with population file : %s\n",popfile);goto end;
       }
       popage=ivector(0,AGESUP);
       popeffectif=vector(0,AGESUP);
       popcount=vector(0,AGESUP);
   
       i=1;  
       while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)
         {
           i=i+1;
         }
       imx=i;
      
       for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
     }
   
     for(cptcov=1;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
         k=k+1;
         fprintf(ficresf,"\n#******");
         for(j=1;j<=cptcoveff;j++) {
           fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         }
         fprintf(ficresf,"******\n");
         fprintf(ficresf,"# StartingAge FinalAge");
         for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
         if (popforecast==1)  fprintf(ficresf," [Population]");
      
         for (cpt=0; cpt<4;cpt++) {
           fprintf(ficresf,"\n");
           fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
   
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
           nhstepm = nhstepm/hstepm;
           /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/
   
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;
           hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
          
           for (h=0; h<=nhstepm; h++){
             if (h==(int) (calagedate+YEARM*cpt)) {
               fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
             }
             for(j=1; j<=nlstate+ndeath;j++) {
               kk1=0.;kk2=0;
               for(i=1; i<=nlstate;i++) {        
                 if (mobilav==1)
                   kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
                 else {
                   kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
                   /* fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h], probs[(int)(agedeb)+1][i][cptcod]);*/
                 }
   
                 if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
               }
            
               if (h==(int)(calagedate+12*cpt)){
                 fprintf(ficresf," %.3f", kk1);
                
                 if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
               }
             }
           }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         }
         }
       }
     }
     if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     if (popforecast==1) {
       free_ivector(popage,0,AGESUP);
       free_vector(popeffectif,0,AGESUP);
       free_vector(popcount,0,AGESUP);
     }
     free_imatrix(s,1,maxwav+1,1,n);
     free_vector(weight,1,n);
     fclose(ficresf);
     }/* End forecasting */
     else{
       erreur=108;
       printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);
     }
   
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");    strcpy(filerest,"t");
Line 2407  fclose(fichtm); Line 3040  fclose(fichtm);
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;        k=k+1;
       fprintf(ficrest,"\n#****** ");        fprintf(ficrest,"\n#****** ");
       for(j=1;j<=cptcovn;j++)        for(j=1;j<=cptcoveff;j++)
         fprintf(ficrest,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);          fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficrest,"******\n");        fprintf(ficrest,"******\n");
   
       fprintf(ficreseij,"\n#****** ");        fprintf(ficreseij,"\n#****** ");
       for(j=1;j<=cptcovn;j++)        for(j=1;j<=cptcoveff;j++)
         fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);          fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
       fprintf(ficreseij,"******\n");        fprintf(ficreseij,"******\n");
   
       fprintf(ficresvij,"\n#****** ");        fprintf(ficresvij,"\n#****** ");
       for(j=1;j<=cptcovn;j++)        for(j=1;j<=cptcoveff;j++)
         fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);          fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
       fprintf(ficresvij,"******\n");        fprintf(ficresvij,"******\n");
   
Line 2437  fclose(fichtm); Line 3070  fclose(fichtm);
       epj=vector(1,nlstate+1);        epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){        for(age=bage; age <=fage ;age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);          prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
           if (popbased==1) {
             for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];
           }
          
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %.0f",age);
         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){          for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
           for(i=1, epj[j]=0.;i <=nlstate;i++) {            for(i=1, epj[j]=0.;i <=nlstate;i++) {
Line 2456  fclose(fichtm); Line 3094  fclose(fichtm);
     }      }
   }    }
                 
          
   
   
  fclose(ficreseij);   fclose(ficreseij);
  fclose(ficresvij);   fclose(ficresvij);
   fclose(ficrest);    fclose(ficrest);
Line 2478  strcpy(fileresvpl,"vpl"); Line 3119  strcpy(fileresvpl,"vpl");
    for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
      k=k+1;       k=k+1;
      fprintf(ficresvpl,"\n#****** ");       fprintf(ficresvpl,"\n#****** ");
      for(j=1;j<=cptcovn;j++)       for(j=1;j<=cptcoveff;j++)
        fprintf(ficresvpl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);         fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
      fprintf(ficresvpl,"******\n");       fprintf(ficresvpl,"******\n");
             
      varpl=matrix(1,nlstate,(int) bage, (int) fage);       varpl=matrix(1,nlstate,(int) bage, (int) fage);
Line 2501  strcpy(fileresvpl,"vpl"); Line 3142  strcpy(fileresvpl,"vpl");
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
     
   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(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
   printf("End of Imach\n");    if(erreur >0)
       printf("End of Imach with error %d\n",erreur);
     else   printf("End of Imach\n");
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
     
   /* 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);*/    /* 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);*/
   /*printf("Total time was %d uSec.\n", total_usecs);*/    /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/    /*------ End -----------*/
   
   
  end:   end:
 #ifdef windows  #ifdef windows
  chdir(pathcd);    /* chdir(pathcd);*/
 #endif  #endif
  system("wgnuplot ../gp37mgw/graph.plt");   /*system("wgnuplot graph.plt");*/
    /*system("../gp37mgw/wgnuplot graph.plt");*/
    /*system("cd ../gp37mgw");*/
    /* system("..\\gp37mgw\\wgnuplot graph.plt");*/
    strcpy(plotcmd,GNUPLOTPROGRAM);
    strcat(plotcmd," ");
    strcat(plotcmd,optionfilegnuplot);
    system(plotcmd);
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
     chdir(pathcd);      chdir(path);
     printf("\nType e to edit output files, c to start again, and q for exiting: ");      printf("\nType e to edit output files, c to start again, and q for exiting: ");
     scanf("%s",z);      scanf("%s",z);
     if (z[0] == 'c') system("./imach");      if (z[0] == 'c') system("./imach");
     else if (z[0] == 'e') {      else if (z[0] == 'e') {
       chdir(path);        chdir(path);
       system("index.htm");        system(optionfilehtm);
     }      }
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }

Removed from v.1.2  
changed lines
  Added in v.1.22


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