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

Diff for /imach/src/imach.c between versions 1.6 and 1.12

version 1.6, 2001/05/02 17:47:10 version 1.12, 2002/02/20 16:57:00
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 68 Line 68
   
   
 int nvar;  int nvar;
 static int cptcov;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 int cptcovn, cptcovage=0;  
 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 */
Line 77  int ncovmodel, ncov;     /* Total number Line 76  int ncovmodel, ncov;     /* Total number
   
 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;
Line 131  double **pmmij; Line 132  double **pmmij;
 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 657  double **prevalim(double **prlim, int nl Line 658  double **prevalim(double **prlim, int nl
      cov[2]=agefin;       cov[2]=agefin;
     
       for (k=1; k<=cptcovn;k++) {        for (k=1; k<=cptcovn;k++) {
         cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];          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]);*/
 /*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/  
       }        }
       for (k=1; k<=cptcovage;k++)        for (k=1; k<=cptcovage;k++)
         cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          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 688  double **prevalim(double **prlim, int nl Line 693  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 711  double **pmij(double **ps, double *cov, Line 716  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 735  double **pmij(double **ps, double *cov, Line 742  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 798  double ***hpxij(double ***po, int nhstep Line 806  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 832  double func( double *x) Line 844  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.;    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 (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
        for(mi=1; mi<= wav[i]-1; mi++){      for(mi=1; mi<= wav[i]-1; mi++){
       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[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++) {          for (kk=1; kk<=cptcovage;kk++) {
                  cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];            cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
                  /*printf("%d %d",kk,Tage[kk]);*/          }
               }         
               /*cov[4]=covar[1][i]*cov[2];scanf("%d", i);*/          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
               /*cov[3]=pow(cov[2],2)/1000.;*/                       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 907  void hesscov(double **matcov, double p[] Line 916  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 915  void hesscov(double **matcov, double p[] Line 923  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 1026  double hessii( double x[], double delta, Line 1036  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;    return res;
     
 }  }
   
Line 1160  void  freqsummary(char fileres[], int ag Line 1170  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 1175  void  freqsummary(char fileres[], int ag Line 1186  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=firstpass; m<=lastpass-1; m++){
Line 1188  void  freqsummary(char fileres[], int ag Line 1200  void  freqsummary(char fileres[], int ag
          }           }
        }         }
         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 1259  void  concatwav(int wav[], int **dh, int Line 1271  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 1291  float sum=0.; Line 1308  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 1313  float sum=0.; Line 1338  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=1, k, j, i, Ntvar[20];    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 1344  void tricode(int *Tvar, int **nbcode, in Line 1374  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 1407  void varevsij(char fileres[], double *** Line 1453  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 1599  void varprevlim(char fileres[], double * Line 1645  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 1611  int main() Line 1657  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];
   char filerest[FILENAMELENGTH];    char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];    char fileregp[FILENAMELENGTH];
Line 1620  int main() Line 1666  int main()
   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 hstepm, nhstepm;    int hstepm, nhstepm;
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 1636  int main() Line 1682  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 ";    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;
Line 1649  int main() Line 1695  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.64b");
   printf("\nEnter the parameter file name: ");    printf("\nEnter the parameter file name: ");
   
 #ifdef windows  #ifdef windows
Line 1698  split(pathtot, path,optionfile); Line 1744  split(pathtot, path,optionfile);
   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);
   
   covar=matrix(0,NCOVMAX,1,n);        covar=matrix(0,NCOVMAX,1,n);
   if (strlen(model)<=1) cptcovn=0;    cptcovn=0;
   else {    if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
     j=0;  
     j=nbocc(model,'+');  
     cptcovn=j+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 1735  split(pathtot, path,optionfile); Line 1777  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 1825  split(pathtot, path,optionfile); Line 1868  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 1847  split(pathtot, path,optionfile); Line 1890  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,15);        Tvar=ivector(1,15);
     Tprod=ivector(1,15);
     Tvaraff=ivector(1,15);
     Tvard=imatrix(1,15,1,2);
   Tage=ivector(1,15);          Tage=ivector(1,15);      
         
   if (strlen(model) >1){    if (strlen(model) >1){
     j=0, j1=0;      j=0, j1=0, k1=1, k2=1;
     j=nbocc(model,'+');      j=nbocc(model,'+');
     j1=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)){
       if (j1==0){        printf("Error. Non available option model=%s ",model);
        cutv(stra,strb,modelsav,'V');        goto end;
        Tvar[1]=atoi(strb);  
       }  
       else if (j1==1) {  
        cutv(stra,strb,modelsav,'*');  
        /*      printf("stra=%s strb=%s modelsav=%s ",stra,strb,modelsav);*/  
        Tage[1]=1; cptcovage++;  
        if (strcmp(stra,"age")==0) {  
          cutv(strd,strc,strb,'V');  
          Tvar[1]=atoi(strc);  
        }  
        else if (strcmp(strb,"age")==0) {  
          cutv(strd,strc,stra,'V');  
          Tvar[1]=atoi(strc);  
        }  
        else {printf("Error"); exit(0);}  
       }  
     }      }
     else {     
       for(i=j; i>=1;i--){      for(i=(j+1); i>=1;i--){
         cutv(stra,strb,modelsav,'+');        cutv(stra,strb,modelsav,'+');
         /*printf("%s %s %s\n", stra,strb,modelsav);*/        if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);
         if (strchr(strb,'*')) {        /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
           cutv(strd,strc,strb,'*');        /*scanf("%d",i);*/
           if (strcmp(strc,"age")==0) {        if (strchr(strb,'*')) {
             cutv(strb,stre,strd,'V');          cutv(strd,strc,strb,'*');
             Tvar[i+1]=atoi(stre);          if (strcmp(strc,"age")==0) {
             cptcovage++;            cptcovprod--;
             Tage[cptcovage]=i+1;            cutv(strb,stre,strd,'V');
             printf("stre=%s ", stre);            Tvar[i]=atoi(stre);
           }            cptcovage++;
           else if (strcmp(strd,"age")==0) {              Tage[cptcovage]=i;
             cutv(strb,stre,strc,'V');              /*printf("stre=%s ", stre);*/
             Tvar[i+1]=atoi(stre);          }
             cptcovage++;          else if (strcmp(strd,"age")==0) {
             Tage[cptcovage]=i+1;            cptcovprod--;
           }            cutv(strb,stre,strc,'V');
           else {            Tvar[i]=atoi(stre);
             cutv(strb,stre,strc,'V');            cptcovage++;
             Tvar[i+1]=ncov+1;            Tage[cptcovage]=i;
             cutv(strb,strc,strd,'V');  
             for (k=1; k<=lastobs;k++)  
               covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];  
           }  
         }          }
         else {          else {
           cutv(strd,strc,strb,'V');            cutv(strb,stre,strc,'V');
           /* printf("%s %s %s", strd,strc,strb);*/            Tvar[i]=ncov+k1;
             cutv(strb,strc,strd,'V');
         Tvar[i+1]=atoi(strc);            Tprod[k1]=i;
         }            Tvard[k1][1]=atoi(strc);
         strcpy(modelsav,stra);              Tvard[k1][2]=atoi(stre);
       }            Tvar[cptcovn+k2]=Tvard[k1][1];
       cutv(strd,strc,stra,'V');            Tvar[cptcovn+k2+1]=Tvard[k1][2];
       Tvar[1]=atoi(strc);            for (k=1; k<=lastobs;k++)
               covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
             k1++;
             k2=k2+2;
           }
         }
         else {
           /*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 %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);    /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
      scanf("%d ",i);*/    printf("cptcovprod=%d ", cptcovprod);
     scanf("%d ",i);*/
     fclose(fic);      fclose(fic);
   
    if(mle==1){      /*  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 1947  split(pathtot, path,optionfile); Line 2005  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 2007  printf("Total number of individuals= %d, Line 2067  printf("Total number of individuals= %d,
   
   
       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 2026  printf("Total number of individuals= %d, Line 2086  printf("Total number of individuals= %d,
      }       }
    }     }
   
    /* 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");
    }     }
Line 2036  printf("Total number of individuals= %d, Line 2097  printf("Total number of individuals= %d,
         
    /* 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);
   
     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) */
      
     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,"\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);
Line 2073  printf("Total number of individuals= %d, Line 2135  printf("Total number of individuals= %d,
          }           }
      }       }
    }     }
    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 2131  printf("Total number of individuals= %d, Line 2194  printf("Total number of individuals= %d,
   
     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 2139  chdir(pathcd); Line 2204  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 2260  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2325  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++) {
             if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
               fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
               ij++;
             }
             else
           fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);            fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
         fprintf(ficgp,")/(1");          }
             fprintf(ficgp,")/(1");
                 
         for(k1=1; k1 <=nlstate; k1++){            for(k1=1; k1 <=nlstate; k1++){  
           fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);            fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
           for(j=3; j <=ncovmodel; j++)  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,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
             }
           fprintf(ficgp,")");            fprintf(ficgp,")");
         }          }
         fprintf(ficgp,") t \"p%d%d\" ", k2,k);          fprintf(ficgp,") t \"p%d%d\" ", k2,k);
Line 2297  chdir(path); Line 2375  chdir(path);
     /*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 2318  chdir(path); Line 2396  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.64b </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 2331  chdir(path); Line 2416  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><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);
   
  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 2343  chdir(path); Line 2428  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 2394  fclose(fichtm); Line 2479  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++){
Line 2402  fclose(fichtm); Line 2487  fclose(fichtm);
         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 2436  fclose(fichtm); Line 2521  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 2495  fclose(fichtm); Line 2580  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 2566  strcpy(fileresvpl,"vpl"); Line 2651  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 2602  strcpy(fileresvpl,"vpl"); Line 2687  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("wgnuplot graph.plt");*/
    /*system("../gp37mgw/wgnuplot graph.plt");*/
    /*system("cd ../gp37mgw");*/
    system("..\\gp37mgw\\wgnuplot graph.plt");
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
Line 2616  strcpy(fileresvpl,"vpl"); Line 2705  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.6  
changed lines
  Added in v.1.12

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