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

Diff for /imach/src/imach.c between versions 1.53 and 1.54

version 1.53, 2002/07/23 23:59:37 version 1.54, 2002/07/24 09:07:45
Line 39 Line 39
   hPijx.    hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the prevalence limits.     of the life expectancies. It also computes the stable prevalence. 
       
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).    Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.             Institut national d'études démographiques, Paris.
Line 83 Line 83
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
 #endif  #endif
   
 char version[80]="Imach version 0.8j, July 2002, INED-EUROREVES ";  char version[80]="Imach version 0.8k, July 2002, 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 701  void powell(double p[], double **xi, int Line 701  void powell(double p[], double **xi, int
         printf("\n");          printf("\n");
         fprintf(ficlog,"\n");          fprintf(ficlog,"\n");
 #endif  #endif
       }         }
     }       } 
   }     } 
 }   } 
   
 /**** Prevalence limit ****************/  /**** Prevalence limit (stable prevalence)  ****************/
   
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)  double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
 {  {
Line 1827  void varevsij(char optionfilefiname[], d Line 1827  void varevsij(char optionfilefiname[], d
     strcpy(digitp,"-populbased-");      strcpy(digitp,"-populbased-");
   else    else
     strcpy(digitp,"-stablbased-");      strcpy(digitp,"-stablbased-");
   if(mobilav==1)    if(mobilav!=0)
     strcat(digitp,"mobilav-");      strcat(digitp,"mobilav-");
   else    else
     strcat(digitp,"nomobil-");      strcat(digitp,"nomobil-");
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, bage, fage, mobaverage);      if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   strcpy(fileresprobmorprev,"prmorprev");     strcpy(fileresprobmorprev,"prmorprev"); 
Line 1928  void varevsij(char optionfilefiname[], d Line 1931  void varevsij(char optionfilefiname[], d
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
   
       if (popbased==1) {        if (popbased==1) {
         if(mobilav !=1){          if(mobilav ==0){
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];              prlim[i][i]=probs[(int)age][i][ij];
         }else{ /* mobilav=1 */           }else{ /* mobilav */ 
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=mobaverage[(int)age][i][ij];              prlim[i][i]=mobaverage[(int)age][i][ij];
         }          }
Line 1956  void varevsij(char optionfilefiname[], d Line 1959  void varevsij(char optionfilefiname[], d
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
     
       if (popbased==1) {        if (popbased==1) {
         if(mobilav !=1){          if(mobilav ==0){
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];              prlim[i][i]=probs[(int)age][i][ij];
         }else{ /* mobilav=1 */           }else{ /* mobilav */ 
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=mobaverage[(int)age][i][ij];              prlim[i][i]=mobaverage[(int)age][i][ij];
         }          }
Line 2025  void varevsij(char optionfilefiname[], d Line 2028  void varevsij(char optionfilefiname[], d
     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);      prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
     
     if (popbased==1) {      if (popbased==1) {
       if(mobilav !=1){        if(mobilav ==0){
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=probs[(int)age][i][ij];            prlim[i][i]=probs[(int)age][i][ij];
       }else{ /* mobilav=1 */         }else{ /* mobilav */ 
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=mobaverage[(int)age][i][ij];            prlim[i][i]=mobaverage[(int)age][i][ij];
       }        }
Line 2106  void varprevlim(char fileres[], double * Line 2109  void varprevlim(char fileres[], double *
   double age,agelim;    double age,agelim;
   int theta;    int theta;
         
   fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");    fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
   fprintf(ficresvpl,"# Age");    fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);        fprintf(ficresvpl," %1d-%1d",i,i);
Line 2599  void printinggnuplot(char fileres[], dou Line 2602  void printinggnuplot(char fileres[], dou
     fprintf(ficlog,"Problem with file %s",optionfilegnuplot);      fprintf(ficlog,"Problem with file %s",optionfilegnuplot);
   }    }
   
 #ifdef windows    /*#ifdef windows */
     fprintf(ficgp,"cd \"%s\" \n",pathc);      fprintf(ficgp,"cd \"%s\" \n",pathc);
 #endif      /*#endif */
 m=pow(2,cptcoveff);  m=pow(2,cptcoveff);
       
  /* 1eme*/   /* 1eme*/
Line 2614  m=pow(2,cptcoveff); Line 2617  m=pow(2,cptcoveff);
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," \%%*lf (\%%*lf)");
      }       }
      fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);       fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," \%%*lf (\%%*lf)");
Line 2768  m=pow(2,cptcoveff); Line 2771  m=pow(2,cptcoveff);
   
   
 /*************** Moving average **************/  /*************** Moving average **************/
 void movingaverage(double ***probs, double bage,double fage, double ***mobaverage){  int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
   
   int i, cpt, cptcod;    int i, cpt, cptcod;
     int mobilavrange, mob;
   double age;    double age;
   for (age=bage; age<=fage; age++)    if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
     for (i=1; i<=nlstate;i++)      if(mobilav==1) mobilavrange=5; /* default */
       for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)      else mobilavrange=mobilav;
         mobaverage[(int)age][i][cptcod]=0.;      for (age=bage; age<=fage; age++)
           for (i=1; i<=nlstate;i++)
   for (age=bage+4; age<=fage; age++){          for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
     for (i=1; i<=nlstate;i++){            mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
       for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      /* We keep the original values on the extreme ages bage, fage and for 
         for (cpt=0;cpt<=4;cpt++){         fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
           mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]+probs[(int)age-cpt][i][cptcod];         we use a 5 terms etc. until the borders are no more concerned. 
       */ 
       for (mob=3;mob <=mobilavrange;mob=mob+2){
         for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
           for (i=1; i<=nlstate;i++){
             for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
               mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
                 for (cpt=1;cpt<=(mob-1)/2;cpt++){
                   mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
                   mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
                 }
               mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
             }
         }          }
         mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]/5;        }/* end age */
       }      }/* end mob */
     }    }else return -1;
   }    return 0;
     }/* End movingaverage */
 }  
   
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
Line 2818  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2833  calagedate=(anproj1+mproj1/12.+jproj1/36
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar,fage, mobaverage);      if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
Line 2891  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2909  calagedate=(anproj1+mproj1/12.+jproj1/36
     }      }
   }    }
                 
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   fclose(ficresf);    fclose(ficresf);
 }  }
Line 2924  populforecast(char fileres[], double anp Line 2942  populforecast(char fileres[], double anp
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar, fage, mobaverage);      if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
Line 3039  populforecast(char fileres[], double anp Line 3060  populforecast(char fileres[], double anp
    }      } 
   }    }
     
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   if (popforecast==1) {    if (popforecast==1) {
     free_ivector(popage,0,AGESUP);      free_ivector(popage,0,AGESUP);
Line 3750  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3771  while((c=getc(ficpar))=='#' && c!= EOF){
   
  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
 /*------------ gnuplot -------------*/  /*------------ gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);   strcpy(optionfilegnuplot,optionfilefiname);
   strcat(optionfilegnuplot,".gp");   strcat(optionfilegnuplot,".gp");
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     printf("Problem with file %s",optionfilegnuplot);     printf("Problem with file %s",optionfilegnuplot);
   }   }
   fclose(ficgp);   else{
      fprintf(ficgp,"\n# %s\n", version); 
      fprintf(ficgp,"# %s\n", optionfilegnuplot); 
      fprintf(ficgp,"set missing 'NaNq'\n");
   }
    fclose(ficgp);
  printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);   printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
Line 3792  Interval (in months) between two waves: Line 3818  Interval (in months) between two waves:
  fclose(ficres);   fclose(ficres);
   
   
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit  (stable prevalence) --------------*/
       
   strcpy(filerespl,"pl");    strcpy(filerespl,"pl");
   strcat(filerespl,fileres);    strcat(filerespl,fileres);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {    if((ficrespl=fopen(filerespl,"w"))==NULL) {
     printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;      printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
     fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;      fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
   }    }
   printf("Computing prevalence limit: result on file '%s' \n", filerespl);    printf("Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);    fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficrespl,"#Prevalence limit\n");    fprintf(ficrespl,"#Stable prevalence \n");
   fprintf(ficrespl,"#Age ");    fprintf(ficrespl,"#Age ");
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");    fprintf(ficrespl,"\n");
Line 3950  Interval (in months) between two waves: Line 3976  Interval (in months) between two waves:
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   calagedate=-1;    calagedate=-1;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);    prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar, fage, mobaverage);      if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   k=0;    k=0;
Line 3994  Interval (in months) between two waves: Line 4023  Interval (in months) between two waves:
       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);
         if (popbased==1) {          if (popbased==1) {
           if(mobilav !=1){            if(mobilav ==0){
             for(i=1; i<=nlstate;i++)              for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];                prlim[i][i]=probs[(int)age][i][k];
           }else{ /* mobilav=1 */             }else{ /* mobilav */ 
             for(i=1; i<=nlstate;i++)              for(i=1; i<=nlstate;i++)
               prlim[i][i]=mobaverage[(int)age][i][k];                prlim[i][i]=mobaverage[(int)age][i][k];
           }            }
Line 4032  free_matrix(mint,1,maxwav,1,n); Line 4061  free_matrix(mint,1,maxwav,1,n);
   fclose(ficpar);    fclose(ficpar);
   free_vector(epj,1,nlstate+1);    free_vector(epj,1,nlstate+1);
       
   /*------- Variance limit prevalence------*/       /*------- Variance of stable prevalence------*/   
   
   strcpy(fileresvpl,"vpl");    strcpy(fileresvpl,"vpl");
   strcat(fileresvpl,fileres);    strcat(fileresvpl,fileres);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);      printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);
     exit(0);      exit(0);
   }    }
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);    printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
   
   k=0;    k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
Line 4075  free_matrix(mint,1,maxwav,1,n); Line 4104  free_matrix(mint,1,maxwav,1,n);
   free_vector(delti,1,npar);    free_vector(delti,1,npar);
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   fprintf(fichtm,"\n</body>");    fprintf(fichtm,"\n</body>");
   fclose(fichtm);    fclose(fichtm);
Line 4109  free_matrix(mint,1,maxwav,1,n); Line 4138  free_matrix(mint,1,maxwav,1,n);
  strcpy(plotcmd,GNUPLOTPROGRAM);   strcpy(plotcmd,GNUPLOTPROGRAM);
  strcat(plotcmd," ");   strcat(plotcmd," ");
  strcat(plotcmd,optionfilegnuplot);   strcat(plotcmd,optionfilegnuplot);
    printf("Starting: %s\n",plotcmd);fflush(stdout);
  system(plotcmd);   system(plotcmd);
   
 #ifdef windows   /*#ifdef windows*/
   while (z[0] != 'q') {    while (z[0] != 'q') {
     /* chdir(path); */      /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");      printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
Line 4121  free_matrix(mint,1,maxwav,1,n); Line 4151  free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);      else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }
 #endif     /*#endif */
 }  }
   
   
Removed from v.1.53  
changed lines
  Added in v.1.54

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