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

Diff for /imach/src/imach.c between versions 1.28 and 1.33

version 1.28, 2002/03/01 17:56:23 version 1.33, 2002/03/12 22:13:38
Line 1443  void  concatwav(int wav[], int **dh, int Line 1443  void  concatwav(int wav[], int **dh, int
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           if (j <= jmin) jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           /* if (j<10) printf("j=%d num=%d ",j,i); */            /*if (j<0) printf("j=%d num=%d \n",j,i); */
           }            }
         }          }
         else{          else{
Line 1451  void  concatwav(int wav[], int **dh, int Line 1451  void  concatwav(int wav[], int **dh, int
           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); */            /*        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 1527  void tricode(int *Tvar, int **nbcode, in Line 1527  void tricode(int *Tvar, int **nbcode, in
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)
 {  {
   /* Health expectancies */    /* Health expectancies */
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h, nstepm, k;
   double age, agelim,hf;    double age, agelim,hf;
   double ***p3mat;    double ***p3mat;
     
Line 1538  void evsij(char fileres[], double ***eij Line 1538  void evsij(char fileres[], double ***eij
       fprintf(ficreseij," %1d-%1d",i,j);        fprintf(ficreseij," %1d-%1d",i,j);
   fprintf(ficreseij,"\n");    fprintf(ficreseij,"\n");
   
   hstepm=1*YEARM; /*  Every j years of age (in month) */    k=1;             /* For example stepm=6 months */
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */    hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */
     hstepm=stepm;   /* or (b) We decided to compute the life expectancy with the smallest unit */
     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
        nhstepm is the number of hstepm from age to agelim
        nstepm is the number of stepm from age to agelin.
        Look at hpijx to understand the reason of that which relies in memory size
        and note for a fixed period like k years */
     /* We decided (b) to get a life expectancy respecting the most precise curvature of the
        survival function given by stepm (the optimization length). Unfortunately it
        means that if the survival funtion is printed only each two years of age and if
        you sum them up and add 1 year (area under the trapezoids) you won't get the same
        results. So we changed our mind and took the option of the best precision.
     */
     hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */
   
   agelim=AGESUP;    agelim=AGESUP;
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
     /* nhstepm age range expressed in number of stepm */      /* nhstepm age range expressed in number of stepm */
     nhstepm=(int) rint((agelim-age)*YEARM/stepm);      nstepm=(int) rint((agelim-age)*YEARM/stepm);
     /* Typically if 20 years = 20*12/6=40 stepm */      /* Typically if 20 years nstepm = 20*12/6=40 stepm */
     if (stepm >= YEARM) hstepm=1;      /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);        hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
       hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm; h++){          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
           eij[i][j][(int)age] +=p3mat[i][j][h];            eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
             /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
         }          }
          fprintf(ficreseij,"%3.0f",age );
     hf=1;  
     if (stepm >= YEARM) hf=stepm/YEARM;  
     fprintf(ficreseij,"%.0f",age );  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         fprintf(ficreseij," %.4f", hf*eij[i][j][(int)age]);          fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficreseij,"\n");
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 1882  fclose(ficresprob); Line 1892  fclose(ficresprob);
 }  }
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[] ){  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[],char rfileres[] ){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   FILE *fichtm;    FILE *fichtm;
   /*char optionfilehtm[FILENAMELENGTH];*/    /*char optionfilehtm[FILENAMELENGTH];*/
Line 1893  void printinghtml(char fileres[], char t Line 1903  void printinghtml(char fileres[], char t
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
  fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">   fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.71e </font> <hr size=\"2\" color=\"#EC5E5E\">
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
   
 Total number of observations=%d <br>  Total number of observations=%d <br>
Line 1911  Interval (in months) between two waves: Line 1921  Interval (in months) between two waves:
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>          - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
         - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>          - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
         - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>          - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%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,fileres,fileres,fileres,fileres);          <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
     
 fprintf(fichtm," <li>Graphs</li><p>");  fprintf(fichtm," <li>Graphs</li><p>");
   
Line 2424  int main(int argc, char *argv[]) Line 2434  int main(int argc, char *argv[])
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];    char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
   char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];    char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
     
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];    char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
   
   char filerest[FILENAMELENGTH];    char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];    char fileregp[FILENAMELENGTH];
Line 2456  int main(int argc, char *argv[]) Line 2466  int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
     
   
   char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";    char version[80]="Imach version 0.71e, March 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
Line 2624  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2634  while((c=getc(ficpar))=='#' && c!= EOF){
   printf("\n");    printf("\n");
   
   
     /*-------- data file ----------*/      /*-------- Rewriting paramater file ----------*/
     if((ficres =fopen(fileres,"w"))==NULL) {       strcpy(rfileres,"r");    /* "Rparameterfile */
       printf("Problem with resultfile: %s\n", fileres);goto end;       strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
        strcat(rfileres,".");    /* */
        strcat(rfileres,optionfilext);    /* Other files have txt extension */
       if((ficres =fopen(rfileres,"w"))==NULL) {
         printf("Problem writing new parameter file: %s\n", fileres);goto end;
     }      }
     fprintf(ficres,"#%s\n",version);      fprintf(ficres,"#%s\n",version);
         
       /*-------- data file ----------*/
     if((fic=fopen(datafile,"r"))==NULL)    {      if((fic=fopen(datafile,"r"))==NULL)    {
       printf("Problem with datafile: %s\n", datafile);goto end;        printf("Problem with datafile: %s\n", datafile);goto end;
     }      }
Line 2902  printf("Total number of individuals= %d, Line 2917  printf("Total number of individuals= %d,
     }      }
         
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"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(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, weightopt,model);
     
   
    jk=1;     jk=1;
Line 3057  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3072  while((c=getc(ficpar))=='#' && c!= EOF){
   
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres);
   
     
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
Line 3160  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3175  while((c=getc(ficpar))=='#' && c!= EOF){
   
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   if(stepm == 1) {    if((stepm == 1) && (strcmp(model,".")==0)){
     prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);      prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
 if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);      if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
     free_matrix(mint,1,maxwav,1,n);      free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);      free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
     free_vector(weight,1,n);}      free_vector(weight,1,n);}
   else{    else{
     erreur=108;      erreur=108;
     printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);      printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
   }    }
     
   
Line 3228  if (popforecast==1) populforecast(filere Line 3243  if (popforecast==1) populforecast(filere
       for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);        for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
       fprintf(ficrest,"\n");        fprintf(ficrest,"\n");
   
       hf=1;  
       if (stepm >= YEARM) hf=stepm/YEARM;  
       epj=vector(1,nlstate+1);        epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){        for(age=bage; age <=fage ;age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);          prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
Line 3238  if (popforecast==1) populforecast(filere Line 3251  if (popforecast==1) populforecast(filere
             prlim[i][i]=probs[(int)age][i][k];              prlim[i][i]=probs[(int)age][i][k];
         }          }
                 
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %4.0f",age);
         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){          for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
           for(i=1, epj[j]=0.;i <=nlstate;i++) {            for(i=1, epj[j]=0.;i <=nlstate;i++) {
             epj[j] += prlim[i][i]*hf*eij[i][j][(int)age];              epj[j] += prlim[i][i]*eij[i][j][(int)age];
           }            }
           epj[nlstate+1] +=epj[j];            epj[nlstate+1] +=epj[j];
         }          }
         for(i=1, vepp=0.;i <=nlstate;i++)          for(i=1, vepp=0.;i <=nlstate;i++)
           for(j=1;j <=nlstate;j++)            for(j=1;j <=nlstate;j++)
             vepp += vareij[i][j][(int)age];              vepp += vareij[i][j][(int)age];
         fprintf(ficrest," %.2f (%.2f)", epj[nlstate+1],hf*sqrt(vepp));          fprintf(ficrest," %7.2f (%7.2f)", epj[nlstate+1],sqrt(vepp));
         for(j=1;j <=nlstate;j++){          for(j=1;j <=nlstate;j++){
           fprintf(ficrest," %.2f (%.2f)", epj[j],hf*sqrt(vareij[j][j][(int)age]));            fprintf(ficrest," %7.2f (%7.2f)", epj[j],sqrt(vareij[j][j][(int)age]));
         }          }
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
       }        }
Removed from v.1.28  
changed lines
  Added in v.1.33

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