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

Diff for /imach/src/imach.c between versions 1.5 and 1.17

version 1.5, 2001/05/02 17:42:45 version 1.17, 2002/02/20 17:15:02
Line 8 Line 8
   Health expectancies are computed from the transistions observed between    Health expectancies are computed from the transistions observed between
   waves and are computed for each degree of severity of disability (number    waves and are computed for each degree of severity of disability (number
   of life states). More degrees you consider, more time is necessary to    of life states). More degrees you consider, more time is necessary to
   reach the Maximum Likelihood of the parameters involved in the model.    reach the Maximum Likelihood of the parameters involved in the model.
   The simplest model is the multinomial logistic model where pij is    The simplest model is the multinomial logistic model where pij is
   the probabibility to be observed in state j at the second wave conditional    the probabibility to be observed in state j at the second wave conditional
   to be observed in state i at the first wave. Therefore the model is:    to be observed in state i at the first wave. Therefore the model is:
Line 22 Line 22
   delay between waves is not identical for each individual, or if some    delay between waves is not identical for each individual, or if some
   individual missed an interview, the information is not rounded or lost, but    individual missed an interview, the information is not rounded or lost, but
   taken into account using an interpolation or extrapolation.    taken into account using an interpolation or extrapolation.
   hPijx is the probability to be    hPijx is the probability to be
   observed in state i at age x+h conditional to the observed state i at age    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    x. The delay 'h' can be split into an exact number (nh*stepm) of
   unobserved intermediate  states. This elementary transition (by month or    unobserved intermediate  states. This elementary transition (by month or
Line 54 Line 54
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
   
   
 #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 70 Line 68
   
   
 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 93  FILE *ficreseij; Line 93  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 129  int stepm; Line 126  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 *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 )  static  int split( char *path, char *dirc, char *name )
 {  {
    char *s;                             /* pointer */     char *s;                             /* pointer */
Line 205  int nbocc(char *s, char occ) Line 202  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;
   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;
Line 214  void cutv(char *u,char *v, char*t, char Line 211  void cutv(char *u,char *v, char*t, char
   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]);
Line 652  double **prevalim(double **prlim, int nl Line 649  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 685  double **prevalim(double **prlim, int nl Line 694  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 708  double **pmij(double **ps, double *cov, Line 717  double **pmij(double **ps, double *cov,
         s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];          s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/          /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }        }
       ps[i][j]=s2;        ps[i][j]=(s2);
     }      }
   }    }
       /*ps[3][2]=1;*/
   
   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 732  double **pmij(double **ps, double *cov, Line 743  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 749  double **pmij(double **ps, double *cov, Line 761  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 795  double ***hpxij(double ***po, int nhstep Line 807  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 819  double ***hpxij(double ***po, int nhstep Line 835  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 829  double func( double *x) Line 845  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[1+k-1][i];          }
             }         
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;          savm=oldm;
           oldm=newm;          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 899  void hesscov(double **matcov, double p[] Line 917  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 907  void hesscov(double **matcov, double p[] Line 924  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 1018  double hessii( double x[], double delta, Line 1037  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;    return res;
     
 }  }
   
Line 1131  void lubksb(double **a, int n, int *indx Line 1150  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, int fprev1,int lprev1,double **mint,double **anint)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;    int i, m, jk, k1, k2,i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos;    double pos;
Line 1142  void  freqsummary(char fileres[], int ag Line 1161  void  freqsummary(char fileres[], int ag
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
     probs= ma3x(1,130 ,1,8, 1,8);
   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 1152  void  freqsummary(char fileres[], int ag Line 1171  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++)
Line 1167  void  freqsummary(char fileres[], int ag Line 1187  void  freqsummary(char fileres[], int ag
        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=fprev1; m<=lprev1; m++){
                k2=anint[m][i]+(mint[m][i]/12.);
                if ((k2>=1984) && (k2<=1988.5)) {
              if(agev[m][i]==0) agev[m][i]=agemax+1;               if(agev[m][i]==0) agev[m][i]=agemax+1;
              if(agev[m][i]==1) agev[m][i]=agemax+2;               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)agev[m][i]] += weight[i];
              freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];               freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
            }               }
          }              }
             }
        }         }
         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 1195  void  freqsummary(char fileres[], int ag Line 1219  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 1205  void  freqsummary(char fileres[], int ag Line 1229  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 1217  void  freqsummary(char fileres[], int ag Line 1243  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 1239  void  freqsummary(char fileres[], int ag Line 1268  void  freqsummary(char fileres[], int ag
   
 }  /* End of Freq */  }  /* End of Freq */
   
   /************ Prevalence ********************/
   void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)
   {  /* Some frequencies */
    
     int i, m, jk, k1, i1, j1, bool, z1,z2,j;
     double ***freq; /* Frequencies */
     double *pp;
     double pos;
   
     pp=vector(1,nlstate);
     probs= ma3x(1,130 ,1,8, 1,8);
    
     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=fprev1; m<=lprev1; m++){
               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]] += weight[i];
               freq[s[m][i]][s[m+1][i]][(int) agemax+3] += 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 */
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
   
 void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)  void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
Line 1251  void  concatwav(int wav[], int **dh, int Line 1357  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 1283  float sum=0.; Line 1394  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));
           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 1305  float sum=0.; Line 1424  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 1336  void tricode(int *Tvar, int **nbcode, in Line 1460  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 1399  void varevsij(char fileres[], double *** Line 1539  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 1435  void varevsij(char fileres[], double *** Line 1575  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 1446  void varevsij(char fileres[], double *** Line 1592  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 1581  void varprevlim(char fileres[], double * Line 1734  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 *****************/
Line 1591  void varprevlim(char fileres[], double * Line 1847  void varprevlim(char fileres[], double *
 int main()  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 1603  int main() Line 1859  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 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,k2,k3,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, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
     int *popage;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
   double **prlim;    double **prlim;
Line 1628  int main() Line 1887  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;
   
     char version[80]="Imach version 64b, May 2001, 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>
Line 1641  int main() Line 1904  int main()
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   
   
   printf("\nIMACH, Version 0.64a");    printf("\nIMACH, Version 0.7");
   printf("\nEnter the parameter file name: ");    printf("\nEnter the parameter file name: ");
   
 #ifdef windows  #ifdef windows
Line 1689  split(pathtot, path,optionfile); Line 1952  split(pathtot, path,optionfile);
   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);
    
     fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
     fprintf(ficparo,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased);
    
    while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
   }    }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);
     fprintf(ficparo,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);
        
    
   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\n",&popforecast,popfile);
    
     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 1727  split(pathtot, path,optionfile); Line 2018  split(pathtot, path,optionfile);
       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 1788  split(pathtot, path,optionfile); Line 2080  split(pathtot, path,optionfile);
   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 1818  split(pathtot, path,optionfile); Line 2109  split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);      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 1840  split(pathtot, path,optionfile); Line 2131  split(pathtot, path,optionfile);
           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;}*/
   
         i=i+1;          i=i+1;
       }        }
     }      }
       /* printf("ii=%d", ij);
     /*scanf("%d",i);*/         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++) 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);*/
         if (strchr(strb,'*')) {
           cutv(strd,strc,strb,'*');
           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 {
             cutv(strb,stre,strc,'V');
             Tvar[i]=ncov+k1;
           cutv(strb,strc,strd,'V');            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++)            for (k=1; k<=lastobs;k++)
             covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];              covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
         }            k1++;
         else {cutv(strd,strc,strb,'V');            k2=k2+2;
         Tvar[i+1]=atoi(strc);  
         }          }
         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]);
Line 1899  split(pathtot, path,optionfile); Line 2246  split(pathtot, path,optionfile);
             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 */
Line 1944  printf("Total number of individuals= %d, Line 2293  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 1958  printf("Total number of individuals= %d, Line 2307  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;
          }           }
        }         }
      }       }
    }     }
   
    /*for(i=1; i <=m ;i++){  
      for(k=1; k <=cptcovn; k++){  
        printf("i=%d k=%d %d ",i,k,codtab[i][k]);  
      }  
      printf("\n");  
    }*/  
    /*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);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev,mint,anint);
    
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n);
    
     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) */
      
     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);
         fprintf(ficres,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased);
      fprintf(ficres,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);
   
    jk=1;     jk=1;
    fprintf(ficres,"# Parameters\n");     fprintf(ficres,"# Parameters\n");
    printf("# Parameters\n");     printf("# Parameters\n");
Line 2025  Tcode=ivector(1,100); Line 2372  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 2083  Tcode=ivector(1,100); Line 2431  Tcode=ivector(1,100);
   
     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 -------------*/  /*------------ gnuplot -------------*/
 chdir(pathcd);  chdir(pathcd);
   if((ficgp=fopen("graph.plt","w"))==NULL) {    if((ficgp=fopen("graph.plt","w"))==NULL) {
Line 2091  chdir(pathcd); Line 2441  chdir(pathcd);
 #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 2187  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2537  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++){
Line 2212  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2562  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
      for(k=1; k<=(nlstate+ndeath); k++) {       for(k=1; k<=(nlstate+ndeath); k++) {
        if (k != k2){         if (k != k2){
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);          fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
   ij=1;
         for(j=3; j <=ncovmodel; j++)          for(j=3; j <=ncovmodel; j++) {
           fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);            if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
         fprintf(ficgp,")/(1");              fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
               ij++;
         for(k1=1; k1 <=nlstate+1; k1=k1+2){              }
             fprintf(ficgp,"+exp(p%d+p%d*x",k1+k3-1,k1+k3);            else
             fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
             for(j=3; j <=ncovmodel; j++)          }
               fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);            fprintf(ficgp,")/(1");
             fprintf(ficgp,")");         
           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);          fprintf(ficgp,") t \"p%d%d\" ", k2,k);
         if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");          if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
     i=i+ncovmodel;          i=i+ncovmodel;
        }         }
      }       }
    }     }
   fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);     fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
    }    }
         
   fclose(ficgp);    fclose(ficgp);
         
 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 2271  chdir(path); Line 2628  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.64a<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 2284  chdir(path); Line 2648  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 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 2296  chdir(path); Line 2662  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 2347  fclose(fichtm); Line 2713  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 2369  fclose(fichtm); Line 2735  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 2378  fclose(fichtm); Line 2745  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 2389  fclose(fichtm); Line 2756  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 2417  fclose(fichtm); Line 2784  fclose(fichtm);
     }      }
   }    }
   
     /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
   
   fclose(ficrespij);    fclose(ficrespij);
   
     /*---------- Forecasting ------------------*/
   
     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);
   
     prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore);
   
    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=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
    
     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 Horizon(in years)");
         for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
         if (popforecast==1)  fprintf(ficresf," [Population]");
   
         for (agedeb=fage; agedeb>=bage; agedeb--){
           fprintf(ficresf,"\n%.f %.f 0",agedeb, agedeb);
          if (mobilav==1) {
           for(j=1; j<=nlstate;j++)
             fprintf(ficresf," %.3f",mobaverage[(int)agedeb][j][cptcod]);
           }
           else {
             for(j=1; j<=nlstate;j++)
             fprintf(ficresf," %.3f",probs[(int)agedeb][j][cptcod]);
           }  
   
          for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");
          if (popforecast==1) fprintf(ficresf," [%.f] ",popeffectif[(int)agedeb]);
         }
        
         for (cpt=1; cpt<=nforecast;cpt++) {
           fprintf(ficresf,"\n");
         for (agedeb=fage; agedeb>=bage; 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*hstepm/YEARM*stepm==cpt)
               fprintf(ficresf,"\n%.f %.f %.f",agedeb, agedeb+ h*hstepm/YEARM*stepm, 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][i][cptcod];
                else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];
                if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
               }
              if (h*hstepm/YEARM*stepm==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);
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");    strcpy(filerest,"t");
Line 2448  fclose(fichtm); Line 2950  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 2478  fclose(fichtm); Line 2980  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 2497  fclose(fichtm); Line 3004  fclose(fichtm);
     }      }
   }    }
                 
          
   
   
  fclose(ficreseij);   fclose(ficreseij);
  fclose(ficresvij);   fclose(ficresvij);
   fclose(ficrest);    fclose(ficrest);
Line 2519  strcpy(fileresvpl,"vpl"); Line 3029  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 2542  strcpy(fileresvpl,"vpl"); Line 3052  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);
     
Line 2555  strcpy(fileresvpl,"vpl"); Line 3065  strcpy(fileresvpl,"vpl");
   /*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 graph.plt");   
    system("..\\gp37mgw\\wgnuplot graph.plt");
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
Line 2569  strcpy(fileresvpl,"vpl"); Line 3081  strcpy(fileresvpl,"vpl");
     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.5  
changed lines
  Added in v.1.17

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