[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.6

version 1.5, 2001/05/02 17:42:45 version 1.6, 2001/05/02 17:47:10
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 71 Line 69
   
 int nvar;  int nvar;
 static int cptcov;  static int cptcov;
 int cptcovn;  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 93  FILE *ficreseij; Line 91  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 124  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;
   
Line 205  int nbocc(char *s, char occ) Line 200  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 209  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 647  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][k]];
     }  
   /*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/
         }
         for (k=1; k<=cptcovage;k++)
           cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
      
     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 819  double ***hpxij(double ***po, int nhstep Line 822  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 831  double func( double *x) Line 834  double func( double *x)
   /*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 (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[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];                   /*printf("%d %d",kk,Tage[kk]);*/
             }                }
                 /*cov[4]=covar[1][i]*cov[2];scanf("%d", i);*/
                 /*cov[3]=pow(cov[2],2)/1000.;*/
   
           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;
Line 1310  float sum=0.; Line 1318  float sum=0.;
 /*********** 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[80],ij=1, k, j, i, Ntvar[20];
   int cptcode=0;    int cptcode=0;
   for (k=0; k<79; k++) Ndum[k]=0;    for (k=0; k<79; 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; 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]);
Line 1337  void tricode(int *Tvar, int **nbcode, in Line 1345  void tricode(int *Tvar, int **nbcode, in
       }          }  
     }      }
   }      }  
    
   }  }
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
Line 1690  split(pathtot, path,optionfile); Line 1698  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(1,NCOVMAX,1,n);        covar=matrix(0,NCOVMAX,1,n);    
   if (strlen(model)<=1) cptcovn=0;    if (strlen(model)<=1) cptcovn=0;
   else {    else {
     j=0;      j=0;
Line 1788  split(pathtot, path,optionfile); Line 1796  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 1851  split(pathtot, path,optionfile); Line 1858  split(pathtot, path,optionfile);
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
   
   /* 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);    
     Tage=ivector(1,15);      
         
   if (strlen(model) >1){    if (strlen(model) >1){
     j=0;      j=0, j1=0;
     j=nbocc(model,'+');      j=nbocc(model,'+');
       j1=nbocc(model,'*');
     cptcovn=j+1;      cptcovn=j+1;
         
     strcpy(modelsav,model);      strcpy(modelsav,model);
     if (j==0) {      if (j==0) {
       cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb);        if (j1==0){
          cutv(stra,strb,modelsav,'V');
          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 {      else {
       for(i=j; i>=1;i--){        for(i=j; i>=1;i--){
         cutv(stra,strb,modelsav,'+');          cutv(stra,strb,modelsav,'+');
           /*printf("%s %s %s\n", stra,strb,modelsav);*/
         if (strchr(strb,'*')) {          if (strchr(strb,'*')) {
           cutv(strd,strc,strb,'*');            cutv(strd,strc,strb,'*');
           cutv(strb,stre,strc,'V');Tvar[i+1]=ncov+1;            if (strcmp(strc,"age")==0) {
           cutv(strb,strc,strd,'V');              cutv(strb,stre,strd,'V');
           for (k=1; k<=lastobs;k++)              Tvar[i+1]=atoi(stre);
             covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];              cptcovage++;
               Tage[cptcovage]=i+1;
               printf("stre=%s ", stre);
             }
             else if (strcmp(strd,"age")==0) {
               cutv(strb,stre,strc,'V');
               Tvar[i+1]=atoi(stre);
               cptcovage++;
               Tage[cptcovage]=i+1;
             }
             else {
               cutv(strb,stre,strc,'V');
               Tvar[i+1]=ncov+1;
               cutv(strb,strc,strd,'V');
               for (k=1; k<=lastobs;k++)
                 covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
             }
         }          }
         else {cutv(strd,strc,strb,'V');          else {
             cutv(strd,strc,strb,'V');
             /* printf("%s %s %s", strd,strc,strb);*/
   
         Tvar[i+1]=atoi(strc);          Tvar[i+1]=atoi(strc);
         }          }
         strcpy(modelsav,stra);            strcpy(modelsav,stra);  
Line 1881  split(pathtot, path,optionfile); Line 1927  split(pathtot, path,optionfile);
       Tvar[1]=atoi(strc);        Tvar[1]=atoi(strc);
     }      }
   }    }
   /*printf("tvar=%d ",Tvar[1]);  
   scanf("%d ",i);*/    /* printf("tvar=%d %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);
        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;
     }      }
Line 1958  printf("Total number of individuals= %d, Line 2006  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(1,nvar,1,8);  
    ncodemax[1]=1;     ncodemax[1]=1;
    if (cptcovn > 0) tricode(Tvar,nbcode,imx);     if (cptcovn > 0) tricode(Tvar,nbcode,imx);
Line 1978  Tcode=ivector(1,100); Line 2026  Tcode=ivector(1,100);
      }       }
    }     }
   
    /*for(i=1; i <=m ;i++){     /* for(i=1; i <=m ;i++){
      for(k=1; k <=cptcovn; k++){       for(k=1; k <=cptcovn; k++){
        printf("i=%d k=%d %d ",i,k,codtab[i][k]);         printf("i=%d k=%d %d ",i,k,codtab[i][k]);
      }       }
      printf("\n");       printf("\n");
    }*/     }
    /*scanf("%d",i);*/     scanf("%d",i);*/
         
    /* Calculates basic frequencies. Computes observed prevalence at single age     /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */         and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);     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 */
Line 2214  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2262  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);          fprintf(ficgp," exp(p%d+p%d*x",i,i+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]]);            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+1; k1=k1+2){            for(k1=1; k1 <=nlstate; k1++){  
             fprintf(ficgp,"+exp(p%d+p%d*x",k1+k3-1,k1+k3);            fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
             for(j=3; j <=ncovmodel; j++)
             for(j=3; j <=ncovmodel; j++)              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",k2,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);
         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);
         
Line 2354  fclose(fichtm); Line 2401  fclose(fichtm);
     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<=cptcovn;j++)
           fprintf(ficrespl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);            fprintf(ficrespl," V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
         fprintf(ficrespl,"******\n");          fprintf(ficrespl,"******\n");
                 
         for (age=agebase; age<=agelim; age++){          for (age=agebase; age<=agelim; age++){
Removed from v.1.5  
changed lines
  Added in v.1.6

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