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

Diff for /imach/src/imach.c between versions 1.75 and 1.76

version 1.75, 2003/05/03 01:18:24 version 1.76, 2003/05/16 10:44:42
Line 119 Line 119
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
 #endif  #endif
   
 char version[80]="Imach version 0.95, May 2003, INED-EUROREVES ";  char version[80]="Imach version 0.95a, May 2003, INED-EUROREVES ";
 int erreur; /* Error number */  int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
Line 141  double jmean; /* Mean space between 2 wa Line 141  double jmean; /* Mean space between 2 wa
 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,*ficresf,*ficrespop;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficlog;  FILE *ficlog, *ficrespow;
 FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficresprobmorprev;  FILE *ficresprobmorprev;
 FILE *fichtm; /* Html File */  FILE *fichtm; /* Html File */
Line 327  void free_vector(double*v, int nl, int n Line 327  void free_vector(double*v, int nl, int n
 }  }
   
 /************************ivector *******************************/  /************************ivector *******************************/
   char *cvector(long nl,long nh)
   {
     char *v;
     v=(char *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(char)));
     if (!v) nrerror("allocation failure in cvector");
     return v-nl+NR_END;
   }
   
   /******************free ivector **************************/
   void free_cvector(char *v, long nl, long nh)
   {
     free((FREE_ARG)(v+nl-NR_END));
   }
   
   /************************ivector *******************************/
 int *ivector(long nl,long nh)  int *ivector(long nl,long nh)
 {  {
   int *v;    int *v;
Line 647  void powell(double p[], double **xi, int Line 662  void powell(double p[], double **xi, int
     del=0.0;       del=0.0; 
     printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
     for (i=1;i<=n;i++)       fprintf(ficrespow,"%d %.12f",*iter,*fret);
       for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);        printf(" %d %.12f",i, p[i]);
     fprintf(ficlog," %d %.12f",i, p[i]);        fprintf(ficlog," %d %.12lf",i, p[i]);
         fprintf(ficrespow," %.12lf", p[i]);
       }
     printf("\n");      printf("\n");
     fprintf(ficlog,"\n");      fprintf(ficlog,"\n");
       fprintf(ficrespow,"\n");
     for (i=1;i<=n;i++) {       for (i=1;i<=n;i++) { 
       for (j=1;j<=n;j++) xit[j]=xi[j][i];         for (j=1;j<=n;j++) xit[j]=xi[j][i]; 
       fptt=(*fret);         fptt=(*fret); 
Line 1187  void mlikeli(FILE *ficres,double p[], in Line 1206  void mlikeli(FILE *ficres,double p[], in
   int i,j, iter;    int i,j, iter;
   double **xi;    double **xi;
   double fret;    double fret;
     char filerespow[FILENAMELENGTH];
   xi=matrix(1,npar,1,npar);    xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)    for (i=1;i<=npar;i++)
     for (j=1;j<=npar;j++)      for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);        xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");    printf("Powell\n");  fprintf(ficlog,"Powell\n");
     strcpy(filerespow,"pow"); 
     strcat(filerespow,fileres);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
     }
     fprintf(ficrespow,"# Powell\n# iter -2*LL");
     for (i=1;i<=nlstate;i++)
       for(j=1;j<=nlstate+ndeath;j++)
         if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
     fprintf(ficrespow,"\n");
   powell(p,xi,npar,ftol,&iter,&fret,func);    powell(p,xi,npar,ftol,&iter,&fret,func);
   
    printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));    fclose(ficrespow);
     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));    fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   
Line 2792  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 2824  fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: pe%s%d1.png<br>
 <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);       <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>
Line 2823  health expectancies in states (1) and (2 Line 2855  health expectancies in states (1) and (2
  - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n   - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
  - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);   - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
   
  if(popforecast==1) fprintf(fichtm,"\n  /*  if(popforecast==1) fprintf(fichtm,"\n */
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n  /*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n  /*  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n */
         <br>",fileres,fileres,fileres,fileres);  /*      <br>",fileres,fileres,fileres,fileres); */
  else   /*  else  */
    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);  /*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
 fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
   
  m=cptcoveff;   m=cptcoveff;
Line 2845  fprintf(fichtm," <ul><li><b>Graphs</b></ Line 2877  fprintf(fichtm," <ul><li><b>Graphs</b></
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident         fprintf(fichtm,"<br>- Observed and period prevalence (with confident
 interval) in state (%d): v%s%d%d.png <br>  interval) in state (%d): v%s%d%d.png <br>
 <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);    <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }       }
Line 2946  m=pow(2,cptcoveff); Line 2978  m=pow(2,cptcoveff);
     }      }
   }    }
       
   /* CV preval stat */    /* CV preval stable (period) */
   for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<nlstate ; cpt ++) {      for (cpt=1; cpt<=nlstate ; cpt ++) {
       k=3;        k=3;
       fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);        fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
               
       for (i=1; i< nlstate ; i ++)        for (i=1; i<= nlstate ; i ++)
         fprintf(ficgp,"+$%d",k+i+1);          fprintf(ficgp,"+$%d",k+i+1);
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);        fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
               
Line 3786  int main(int argc, char *argv[]) Line 3818  int main(int argc, char *argv[])
   
   for (i=1; i<=imx; i++) {    for (i=1; i<=imx; i++) {
     for(m=2; (m<= maxwav); m++) {      for(m=2; (m<= maxwav); m++) {
       if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){        if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
         anint[m][i]=9999;          anint[m][i]=9999;
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
           printf("Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown %d, set an arbitrary year of death\n",(int)moisdc[i],(int)andc[i],num[i],i);
           fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown %d, set an arbitrary year of death\n",(int)moisdc[i],(int)andc[i],num[i],i);
           s[m][i]=-1;
         }
         if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
           printf("Warning! Month of death of individual %d on line %d was unknown %2d, we set it to June\n",num[i],i,(int)moisdc[i]); 
           fprintf(ficlog,"Warning! Month of death of individual %d on line %d was unknown %f, we set it to June\n",num[i],i,moisdc[i]); 
           moisdc[i]=6;
           s[m][i]=-1;
         }
     }      }
   }    }
   
Line 3800  int main(int argc, char *argv[]) Line 3842  int main(int argc, char *argv[])
       if(s[m][i] >0){        if(s[m][i] >0){
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0)            if(agedc[i]>0)
             if(moisdc[i]!=99 && andc[i]!=9999)              if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/            /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             else {              else {
               if (andc[i]!=9999){                if ((int)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);
                 fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);                  fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
                 agev[m][i]=-1;                  agev[m][i]=-1;
Line 3815  int main(int argc, char *argv[]) Line 3857  int main(int argc, char *argv[])
                                  years but with the precision of a                                   years but with the precision of a
                                  month */                                   month */
           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((int)mint[m][i]==99 || (int)anint[m][i]==9999)
             agev[m][i]=1;              agev[m][i]=1;
           else if(agev[m][i] <agemin){             else if(agev[m][i] <agemin){ 
             agemin=agev[m][i];              agemin=agev[m][i];
Line 4113  int main(int argc, char *argv[]) Line 4155  int main(int argc, char *argv[])
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
 \n  \n
 Total number of observations=%d <br>\n  Total number of observations=%d <br>\n
   Youngest age at first pass %.2f, oldest age %.2f<br>\n
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n  Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 <hr  size=\"2\" color=\"#EC5E5E\">  <hr  size=\"2\" color=\"#EC5E5E\">
  <ul><li><h4>Parameter files</h4>\n   <ul><li><h4>Parameter files</h4>\n
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n   - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
  - Log file of the run: <a href=\"%s\">%s</a><br>\n   - Log file of the run: <a href=\"%s\">%s</a><br>\n
  - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);   - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,agemin,agemax,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
    fclose(fichtm);     fclose(fichtm);
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
Removed from v.1.75  
changed lines
  Added in v.1.76

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