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

Diff for /imach/src/imach.c between versions 1.86 and 1.87

version 1.86, 2003/06/17 20:04:08 version 1.87, 2003/06/18 12:26:01
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.87  2003/06/18 12:26:01  brouard
     Version 0.96
   
   Revision 1.86  2003/06/17 20:04:08  brouard    Revision 1.86  2003/06/17 20:04:08  brouard
   (Module): Change position of html and gnuplot routines and added    (Module): Change position of html and gnuplot routines and added
   routine fileappend.    routine fileappend.
Line 161 Line 164
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES ";  char version[]="Imach version 0.96, June 2003, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 int erreur; /* Error number */  int erreur; /* Error number */
 int nvar;  int nvar;
Line 175  int popbased=0; Line 178  int popbased=0;
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */  int maxwav; /* Maxim number of waves */
 int jmin, jmax; /* min, max spacing between 2 waves */  int jmin, jmax; /* min, max spacing between 2 waves */
   int gipmx, gsw; /* Global variables on the number of contributions 
                      to the likelihood and the sum of weights (done by funcone)*/
 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 */
Line 1290  double func( double *x) Line 1295  double func( double *x)
 /*************** log-likelihood *************/  /*************** log-likelihood *************/
 double funcone( double *x)  double funcone( double *x)
 {  {
     /* Same as likeli but slower because of a lot of printf and if */
   int i, ii, j, k, mi, d, kk;    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 lli; /* Individual log likelihood */    double lli; /* Individual log likelihood */
     double llt;
   int s1, s2;    int s1, s2;
   double bbh, survp;    double bbh, survp;
   /*extern weight */    /*extern weight */
Line 1354  double funcone( double *x) Line 1361  double funcone( double *x)
  %10.6f %10.6f %10.6f ", \   %10.6f %10.6f %10.6f ", \
                 num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],                  num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
         for(k=1,l=0.; k<=nlstate; k++)           for(k=1,llt=0.,l=0.; k<=nlstate; k++){
           fprintf(ficresilk," %10.6f",ll[k]);            llt +=ll[k]*gipmx/gsw;
         fprintf(ficresilk,"\n");            fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
           }
           fprintf(ficresilk," %10.6f\n", -llt);
       }        }
     } /* end of wave */      } /* end of wave */
   } /* end of individual */    } /* end of individual */
   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 */
     if(globpr==0){ /* First time we count the contributions and weights */
       gipmx=ipmx;
       gsw=sw;
     }
   return -l;    return -l;
 }  }
   
   
 void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpr, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))  void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
 {  {
   /* This routine should help understanding what is done with the selection of individuals/waves and    /* This routine should help understanding what is done with 
        the selection of individuals/waves and
      to check the exact contribution to the likelihood.       to check the exact contribution to the likelihood.
      Plotting could be done.       Plotting could be done.
    */     */
   int k;    int k;
   if(globpr !=0){ /* Just counts and sums no printings */  
     if(*globpri !=0){ /* Just counts and sums no printings */
     strcpy(fileresilk,"ilk");       strcpy(fileresilk,"ilk"); 
     strcat(fileresilk,fileres);      strcat(fileresilk,fileres);
     if((ficresilk=fopen(fileresilk,"w"))==NULL) {      if((ficresilk=fopen(fileresilk,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresilk);        printf("Problem with resultfile: %s\n", fileresilk);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
     }      }
     fprintf(ficresilk, "#individual(line's record) s1 s2 wave# effective_wave# number_of_product_matrix pij weight 2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state\n");      fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
     fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight out sav ");      fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight out sav ");
     /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
     for(k=1; k<=nlstate; k++)       for(k=1; k<=nlstate; k++) 
       fprintf(ficresilk," ll[%d]",k);        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
     fprintf(ficresilk,"\n");      fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n");
   }    }
   
   *fretone=(*funcone)(p);    *fretone=(*funcone)(p);
   if(globpr !=0){    if(*globpri !=0){
     fclose(ficresilk);      fclose(ficresilk);
     if((fichtm=fopen(optionfilehtm,"a"))==NULL) {      fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk,fileresilk);
       printf("Problem with html file: %s\n", optionfilehtm);      fflush(fichtm); 
       fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);    } 
       exit(0);  
     }  
     else{  
       fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk);  
       fclose(fichtm);  
     }  
   }  
   return;    return;
 }  }
   
Line 2368  void varevsij(char optionfilefiname[], d Line 2376  void varevsij(char optionfilefiname[], d
   else{    else{
     fprintf(ficgp,"\n# Routine varevsij");      fprintf(ficgp,"\n# Routine varevsij");
   }    }
   if((fichtm=fopen(optionfilehtm,"a"))==NULL) {  /*   if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
     printf("Problem with html file: %s\n", optionfilehtm);  /*     printf("Problem with html file: %s\n", optionfilehtm); */
     fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);  /*     fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */
     exit(0);  /*     exit(0); */
   }  /*   } */
   else{  /*   else{ */
     fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
     fprintf(fichtm,"\n<br>%s  <br>\n",digitp);    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
   }  /*   } */
   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
   
   fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");
Line 2606  void varevsij(char optionfilefiname[], d Line 2614  void varevsij(char optionfilefiname[], d
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   fclose(ficresprobmorprev);    fclose(ficresprobmorprev);
   fclose(ficgp);    fclose(ficgp);
   fclose(fichtm);  /*   fclose(fichtm); */
 }  /* end varevsij */  }  /* end varevsij */
   
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
Line 2771  void varprob(char optionfilefiname[], do Line 2779  void varprob(char optionfilefiname[], do
   else{    else{
     fprintf(ficgp,"\n# Routine varprob");      fprintf(ficgp,"\n# Routine varprob");
   }    }
   if((fichtm=fopen(optionfilehtm,"a"))==NULL) {  /*   if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
     printf("Problem with html file: %s\n", optionfilehtm);  /*     printf("Problem with html file: %s\n", optionfilehtm); */
     fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);  /*     fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */
     exit(0);  /*     exit(0); */
   }  /*   } */
   else{  /*   else{ */
     fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");      fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
     fprintf(fichtm,"\n");      fprintf(fichtm,"\n");
   
Line 2784  void varprob(char optionfilefiname[], do Line 2792  void varprob(char optionfilefiname[], do
     fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");      fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
     fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n");      fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n");
   
   }  /*   } */
   
   cov[1]=1;    cov[1]=1;
   tj=cptcoveff;    tj=cptcoveff;
Line 3001  void varprob(char optionfilefiname[], do Line 3009  void varprob(char optionfilefiname[], do
   fclose(ficresprobcov);    fclose(ficresprobcov);
   fclose(ficresprobcor);    fclose(ficresprobcor);
   fclose(ficgp);    fclose(ficgp);
   fclose(fichtm);  
 }  }
   
   
Line 3014  void printinghtml(char fileres[], char t Line 3021  void printinghtml(char fileres[], char t
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   /*char optionfilehtm[FILENAMELENGTH];*/    /*char optionfilehtm[FILENAMELENGTH];*/
   if((fichtm=fopen(optionfilehtm,"a"))==NULL)    {  /*   if((fichtm=fopen(optionfilehtm,"a"))==NULL)    { */
     printf("Problem with %s \n",optionfilehtm), exit(0);  /*     printf("Problem with %s \n",optionfilehtm), exit(0); */
     fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);  /*     fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */
   }  /*   } */
   
    fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \     fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \   - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \
Line 3104  interval) in state (%d): v%s%d%d.png <br Line 3111  interval) in state (%d): v%s%d%d.png <br
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
 fclose(fichtm);   fflush(fichtm);
 }  }
   
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
Line 3611  populforecast(char fileres[], double anp Line 3618  populforecast(char fileres[], double anp
   fclose(ficrespop);    fclose(ficrespop);
 } /* End of popforecast */  } /* End of popforecast */
   
 int fileappend(FILE *fichier, char *optionfile)  int fileappend(FILE *fichier, char *optionfich)
 {  {
   if((fichier=fopen(optionfile,"a"))==NULL) {    if((fichier=fopen(optionfich,"a"))==NULL) {
     printf("Problem with file: %s\n", optionfile);      printf("Problem with file: %s\n", optionfich);
     fprintf(ficlog,"Problem with file: %s\n", optionfile);      fprintf(ficlog,"Problem with file: %s\n", optionfich);
     return (1);      return (0);
   }    }
     fflush(fichier);
     return (1);
 }  }
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
Line 3630  int main(int argc, char *argv[]) Line 3638  int main(int argc, char *argv[])
   int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
   int jj;    int jj;
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
     /*  FILE *fichtm; *//* Html File */
     /* FILE *ficgp;*/ /*Gnuplot File */
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;    double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
   
Line 4268  int main(int argc, char *argv[]) Line 4278  int main(int argc, char *argv[])
   fclose(ficgp);    fclose(ficgp);
   /*--------- index.htm --------*/    /*--------- index.htm --------*/
   
   strcpy(optionfilehtm,optionfile);    strcpy(optionfilehtm,optionfilefiname);
   strcat(optionfilehtm,".htm");    strcat(optionfilehtm,".htm");
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {    if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
   fprintf(fichtm,"<body> <font size=\"2\">%s <br> %s</font> \    fprintf(fichtm,"<body>\n<title>IMaCh %s</title>\n <font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 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\
Line 4282  Title=%s <br>Datafile=%s Firstpass=%d La Line 4292  Title=%s <br>Datafile=%s Firstpass=%d La
  <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>\n\   - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
  - Date and time at start: %s</ul>\n",\   - Date and time at start: %s</ul>\n",\
           version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\            fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\
           model,fileres,fileres,\            model,fileres,fileres,\
           filelog,filelog,optionfilegnuplot,optionfilegnuplot,strt);            filelog,filelog,optionfilegnuplot,optionfilegnuplot,strt);
   fclose(fichtm);    /*fclose(fichtm);*/
     fflush(fichtm);
   
   /* 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,Tvaraff,nbcode, ncodemax,mint,anint);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
   
   if(fileappend(fichtm, optionfilehtm)){  /*   if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */
   /*     printf("Problem with file: %s\n", optionfilehtm); */
   /*     fprintf(ficlog,"Problem with file: %s\n", optionfilehtm); */
   /*   } */
   
   
   /*   if(fileappend(fichtm, optionfilehtm)){ */
       fprintf(fichtm,"\n");
     fprintf(fichtm,"<br>Total number of observations=%d <br>\n\      fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
 Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\  Youngest age at first (selected) 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",\
         imx,agemin,agemax,jmin,jmax,jmean);          imx,agemin,agemax,jmin,jmax,jmean);
     fclose(fichtm);  /*     fclose(fichtm); */
   }  /*   } */
   
     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 */
Line 4837  ageminpar, agemax, s[lastpass][imx], age Line 4856  ageminpar, agemax, s[lastpass][imx], age
   printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);    printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
   fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);    fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);
   /*  printf("Total time was %d uSec.\n", total_usecs);*/    /*  printf("Total time was %d uSec.\n", total_usecs);*/
   if(fileappend(fichtm,optionfilehtm)){  /*   if(fileappend(fichtm,optionfilehtm)){ */
     fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend);    fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend);
     fclose(fichtm);    fclose(fichtm);
   }  
   /*------ End -----------*/    /*------ End -----------*/
   
   end:    end:
Removed from v.1.86  
changed lines
  Added in v.1.87

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