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

Diff for /imach064/src/imach.c between versions 1.2 and 1.5

version 1.2, 2001/03/13 18:10:26 version 1.5, 2001/05/02 17:42:45
Line 51 Line 51
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
 #define windows  #define windows
   #define GLOCK_ERROR_NOPATH              -1      /* empty path */
   #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
   
   
 #define MAXPARM 30 /* Maximum number of parameters for the optimization */  #define MAXPARM 30 /* Maximum number of parameters for the optimization */
 #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */  #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
Line 59 Line 63
 #define NLSTATEMAX 8 /* Maximum number of live states (for func) */  #define NLSTATEMAX 8 /* Maximum number of live states (for func) */
 #define NDEATHMAX 8 /* Maximum number of dead states (for func) */  #define NDEATHMAX 8 /* Maximum number of dead states (for func) */
 #define NCOVMAX 8 /* Maximum number of covariates */  #define NCOVMAX 8 /* Maximum number of covariates */
 #define MAXN 80000  #define MAXN 20000
 #define YEARM 12. /* Number of months per year */  #define YEARM 12. /* Number of months per year */
 #define AGESUP 130  #define AGESUP 130
 #define AGEBASE 40  #define AGEBASE 40
Line 138  double ftol=FTOL; /* Tolerance for compu Line 142  double ftol=FTOL; /* Tolerance for compu
 double ftolhess; /* Tolerance for computing hessian */  double ftolhess; /* Tolerance for computing hessian */
   
   
   static  int split( char *path, char *dirc, char *name )
   {
      char *s;                             /* pointer */
      int  l1, l2;                         /* length counters */
   
      l1 = strlen( path );                 /* length of path */
      if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
      s = strrchr( path, '\\' );           /* find last / */
      if ( s == NULL ) {                   /* no directory, so use current */
   #if     defined(__bsd__)                /* get current working directory */
         extern char       *getwd( );
   
         if ( getwd( dirc ) == NULL ) {
   #else
         extern char       *getcwd( );
   
         if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
   #endif
            return( GLOCK_ERROR_GETCWD );
         }
         strcpy( name, path );             /* we've got it */
      } else {                             /* strip direcotry from path */
         s++;                              /* after this, the filename */
         l2 = strlen( s );                 /* length of filename */
         if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
         strcpy( name, s );                /* save file name */
         strncpy( dirc, path, l1 - l2 );   /* now the directory */
         dirc[l1-l2] = 0;                  /* add zero */
      }
      l1 = strlen( dirc );                 /* length of directory */
      if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
      return( 0 );                         /* we're done */
   }
   
   
 /******************************************/  /******************************************/
   
 void replace(char *s, char*t)  void replace(char *s, char*t)
Line 168  void cutv(char *u,char *v, char*t, char Line 207  void cutv(char *u,char *v, char*t, char
 {  {
   int i,lg,j,p;    int i,lg,j,p;
   i=0;    i=0;
   if (t[0]== occ) p=0;  
   for(j=0; j<=strlen(t)-1; j++) {    for(j=0; j<=strlen(t)-1; j++) {
     if((t[j]!= occ) && (t[j+1]==occ)) p=j+1;      if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
   }    }
   
   lg=strlen(t);    lg=strlen(t);
Line 184  void cutv(char *u,char *v, char*t, char Line 222  void cutv(char *u,char *v, char*t, char
   }    }
 }  }
   
   
 /********************** nrerror ********************/  /********************** nrerror ********************/
   
 void nrerror(char error_text[])  void nrerror(char error_text[])
Line 805  printf(" %d\n",s[4][i]); Line 842  printf(" %d\n",s[4][i]);
           cov[1]=1.;            cov[1]=1.;
           cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;            cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
           if (cptcovn>0){            if (cptcovn>0){
             for (k=1; k<=cptcovn;k++) {              for (k=1; k<=cptcovn;k++) cov[2+k]=covar[1+k-1][i];
               cov[2+k]=covar[Tvar[k]][i];  
               /* printf("k=%d cptcovn=%d %lf\n",k,cptcovn,covar[Tvar[k]][i]);*/  
             }  
             }              }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            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]]);
Line 827  printf(" %d\n",s[4][i]); Line 863  printf(" %d\n",s[4][i]);
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */    /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
   
   return -l;    return -l;
 }  }
   
Line 983  double hessii( double x[], double delta, Line 1018  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;      return res;
    
 }  }
   
 double hessij( double x[], double delti[], int thetai,int thetaj)  double hessij( double x[], double delti[], int thetai,int thetaj)
Line 1252  float sum=0.; Line 1288  float sum=0.;
         }          }
         else{          else{
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));            j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
           /*printf("i=%d agevi+1=%lf agevi=%lf j=%d\n", i,agev[mw[mi+1][i]][i],agev[mw[mi][i]][i],j);*/  
   
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
Line 1578  int main() Line 1612  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,jk,aa,bb, stepsize;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize;
   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;
Line 1596  int main() Line 1630  int main()
   double *epj, vepp;    double *epj, vepp;
   char version[80]="Imach version 62c, May 1999, INED-EUROREVES ";    char version[80]="Imach version 62c, May 1999, 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 1606  int main() Line 1641  int main()
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   
   
   printf("\nIMACH, Version 0.63");    printf("\nIMACH, Version 0.64a");
   printf("\nEnter the parameter file name: ");    printf("\nEnter the parameter file name: ");
   
 #ifdef windows  #ifdef windows
   scanf("%s",pathtot);    scanf("%s",pathtot);
   getcwd(pathcd, size);    getcwd(pathcd, size);
   cutv(path,optionfile,pathtot,'\\');    /*cygwin_split_path(pathtot,path,optionfile);
       printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
     /* cutv(path,optionfile,pathtot,'\\');*/
   
   split(pathtot, path,optionfile);
   chdir(path);    chdir(path);
   replace(pathc,path);    replace(pathc,path);
 #endif  #endif
Line 1777  int main() Line 1816  int main()
     s=imatrix(1,maxwav+1,1,n);      s=imatrix(1,maxwav+1,1,n);
     adl=imatrix(1,maxwav+1,1,n);          adl=imatrix(1,maxwav+1,1,n);    
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,NCOVMAX);      ncodemax=ivector(1,8);
   
     i=1;      i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {      while (fgets(line, MAXLINE, fic) != NULL)    {
Line 1802  int main() Line 1841  int main()
         }          }
         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]));*/          /*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
         /*printf("%d %.lf %.lf %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),(covar[3][i]), (covar[4][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]));*/  
   
         i=i+1;          i=i+1;
       }        }
     }      }
     /*scanf("%d",i);*/  
   
       /*scanf("%d",i);*/
   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,8);    
         
Line 1820  int main() Line 1857  int main()
     j=0;      j=0;
     j=nbocc(model,'+');      j=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);        cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb);
Line 1835  int main() Line 1872  int main()
           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+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
         }          }
         else {          else {cutv(strd,strc,strb,'V');
           cutv(strd,strc,strb,'V');          Tvar[i+1]=atoi(strc);
           Tvar[i+1]=atoi(strc);  
         }          }
         strcpy(modelsav,stra);            strcpy(modelsav,stra);  
       }        }
       /*cutv(strd,strc,stra,'V');*/        cutv(strd,strc,stra,'V');
       Tvar[1]=atoi(strc);        Tvar[1]=atoi(strc);
     }      }
   }    }
   /*printf("tvar=%d ",Tvar[1]);*/    /*printf("tvar=%d ",Tvar[1]);
   /*scanf("%d ",i);*/    scanf("%d ",i);*/
     fclose(fic);      fclose(fic);
   
     if (weightopt != 1) { /* Maximisation without weights*/      if (weightopt != 1) { /* Maximisation without weights*/
Line 1858  int main() Line 1894  int main()
     for (i=1; i<=imx; i++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
       for(m=1; (m<= maxwav); m++){        for(m=1; (m<= maxwav); m++){
         if (mint[m][i]==99 || anint[m][i]==9999) s[m][i]=-1;    
         if(s[m][i] >0){          if(s[m][i] >0){
           if (s[m][i] == nlstate+1) {            if (s[m][i] == nlstate+1) {
             if(agedc[i]>0)              if(agedc[i]>0)
Line 1871  int main() Line 1906  int main()
           }            }
           else if(s[m][i] !=9){ /* Should no more exist */            else if(s[m][i] !=9){ /* Should no more exist */
             agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);              agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
             if(mint[m][i]==99 || anint[m][i]==9999){              if(mint[m][i]==99 || anint[m][i]==9999)
               agev[m][i]=1;                agev[m][i]=1;
               /* printf("i=%d m=%d agev=%lf \n",i,m, agev[m][i]);    */  
             }  
             else if(agev[m][i] <agemin){              else if(agev[m][i] <agemin){
               agemin=agev[m][i];                agemin=agev[m][i];
               /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/                /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
Line 1950  Tcode=ivector(1,100); Line 1983  Tcode=ivector(1,100);
        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'. */
Line 1966  Tcode=ivector(1,100); Line 1999  Tcode=ivector(1,100);
     /* For Powell, parameters are in a vector p[] starting at p[1]      /* For Powell, parameters are in a vector p[] starting at p[1]
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */         so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */      p=param[1][1]; /* *(*(*(param +1)+1)+0) */
     /*scanf("%d",i);*/     
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);      mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   
         
Line 2157  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2190  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
   }    }
   
   /* proba elementaires */    /* proba elementaires */
   for(i=1,jk=1; i <=nlstate; i++){     for(i=1,jk=1; i <=nlstate; i++){
     for(k=1; k <=(nlstate+ndeath); k++){      for(k=1; k <=(nlstate+ndeath); k++){
       if (k != i) {        if (k != i) {
         /*  fprintf(ficgp,"%1d%1d ",i,k);*/  
         for(j=1; j <=ncovmodel; j++){          for(j=1; j <=ncovmodel; j++){
           fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);            /*fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/
             /*fprintf(ficgp,"%s",alph[1]);*/
             fprintf(ficgp,"p%d=%f ",jk,p[jk]);
           jk++;            jk++;
           fprintf(ficgp,"\n");            fprintf(ficgp,"\n");
         }          }
       }        }
     }      }
   }      }
   
   for(jk=1; jk <=m; jk++) {    for(jk=1; jk <=m; jk++) {
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);    fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);
   for(i=1; i <=nlstate; i++) {     i=1;
     for(k=1; k <=(nlstate+ndeath); k++){     for(k2=1; k2<=nlstate; k2++) {
       if (k != i) {       k3=i;
         fprintf(ficgp," exp(a%d%d+b%d%d*x",i,k,i,k);       for(k=1; k<=(nlstate+ndeath); k++) {
          if (k != k2){
           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,"+%s%d%d*%d",alph[j],i,k,nbcode[Tvar[j-2]][codtab[jk][j-2]]);            fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
         fprintf(ficgp,")/(1");          fprintf(ficgp,")/(1");
         for(k1=1; k1 <=(nlstate+ndeath); k1++)  
           if (k1 != i) {          for(k1=1; k1 <=nlstate+1; k1=k1+2){  
             fprintf(ficgp,"+exp(a%d%d+b%d%d*x",i,k1,i,k1);              fprintf(ficgp,"+exp(p%d+p%d*x",k1+k3-1,k1+k3);
   
             for(j=3; j <=ncovmodel; j++)              for(j=3; j <=ncovmodel; j++)
               fprintf(ficgp,"+%s%d%d*%d",alph[j],i,k,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\" ", i,k);          fprintf(ficgp,") t \"p%d%d\" ", k2,k);
       if ((i+k)!= (nlstate*2+ndeath)) fprintf(ficgp,",");          if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
       }      i=i+ncovmodel;
     }         }
   }       }
 fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);       }
   }    fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
      }
      
   fclose(ficgp);    fclose(ficgp);
         
 chdir(path);  chdir(path);
Line 2234  chdir(path); Line 2275  chdir(path);
     printf("Problem with index.htm \n");goto end;      printf("Problem with index.htm \n");goto end;
   }    }
   
  fprintf(fichtm,"<body><ul> Imach, Version 0.63<hr> <li>Outputs files<br><br>\n   fprintf(fichtm,"<body><ul> Imach, Version 0.64a<hr> <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 2518  strcpy(fileresvpl,"vpl"); Line 2559  strcpy(fileresvpl,"vpl");
 #ifdef windows  #ifdef windows
  chdir(pathcd);   chdir(pathcd);
 #endif  #endif
  system("wgnuplot ../gp37mgw/graph.plt");   system("wgnuplot graph.plt");
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
Removed from v.1.2  
changed lines
  Added in v.1.5

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