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

Diff for /imach/src/imach.c between versions 1.25 and 1.26

version 1.25, 2002/02/26 17:11:54 version 1.26, 2002/02/27 15:42:00
Line 1176  void lubksb(double **a, int n, int *indx Line 1176  void lubksb(double **a, int n, int *indx
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2)  void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;    int i, m, jk, k1,i1, j1, bool, z1,z2,j;
Line 1236  void  freqsummary(char fileres[], int ag Line 1236  void  freqsummary(char fileres[], int ag
            }             }
          }           }
        }         }
          
          fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
         if  (cptcovn>0) {          if  (cptcovn>0) {
          fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable ");
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);           for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 1620  void varevsij(char fileres[], double *** Line 1623  void varevsij(char fileres[], double ***
         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];
       }        }
         
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
Line 1632  void varevsij(char fileres[], double *** Line 1635  void varevsij(char fileres[], double ***
         xp[i] = x[i] - (i==theta ?delti[theta]:0);          xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
    
       if (popbased==1) {        if (popbased==1) {
         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];
Line 1684  void varevsij(char fileres[], double *** Line 1687  void varevsij(char fileres[], double ***
     free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);      free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   } /* End age */    } /* End age */
     
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,npar);    free_matrix(doldm,1,nlstate,1,npar);
   free_matrix(dnewm,1,nlstate,1,nlstate);    free_matrix(dnewm,1,nlstate,1,nlstate);
Line 1876  if (i== 4) fprintf(ficresprob,"%.3e %.3e Line 1879  if (i== 4) fprintf(ficresprob,"%.3e %.3e
 }  }
  free_vector(xp,1,npar);   free_vector(xp,1,npar);
 fclose(ficresprob);  fclose(ficresprob);
  exit(0);  
 }  }
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
Line 2763  printf("Total number of individuals= %d, Line 2766  printf("Total number of individuals= %d,
     fputs(line,ficparo);      fputs(line,ficparo);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mob_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);    fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
 fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);  fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
 fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);  fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
   
  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
 /*------------ gnuplot -------------*/  /*------------ gnuplot -------------*/
  printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);   printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);
Line 2903  fprintf(ficres,"popforecast=%d popfile=% Line 2906  fprintf(ficres,"popforecast=%d popfile=%
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;    calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
   
   /*printf("calage= %f", calagedate);*/   
     
   prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);    prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
    
   strcpy(fileresf,"f");    strcpy(fileresf,"f");
   strcat(fileresf,fileres);    strcat(fileresf,fileres);
   if((ficresf=fopen(fileresf,"w"))==NULL) {    if((ficresf=fopen(fileresf,"w"))==NULL) {
     printf("Problem with forecast resultfile: %s\n", fileresf);goto end;      printf("Problem with forecast resultfile: %s\n", fileresf);goto end;
   }    }
   printf("Computing forecasting: result on file '%s' \n", fileresf);    printf("Computing forecasting: result on file '%s' \n", fileresf);
    
   free_matrix(mint,1,maxwav,1,n);    free_matrix(mint,1,maxwav,1,n);
   free_matrix(anint,1,maxwav,1,n);    free_matrix(anint,1,maxwav,1,n);
   free_matrix(agev,1,maxwav,1,imx);  
   
   /* Mobile average */  
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
Line 2930  fprintf(ficres,"popforecast=%d popfile=% Line 2930  fprintf(ficres,"popforecast=%d popfile=%
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
   if (stepm<=12) stepsize=1;    if (stepm<=12) stepsize=1;
    
   agelim=AGESUP;    agelim=AGESUP;
   /*hstepm=stepsize*YEARM; *//* Every year of age */   
   hstepm=1;    hstepm=1;
   hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */    hstepm=hstepm/stepm;
   yp1=modf(dateintmean,&yp);    yp1=modf(dateintmean,&yp);
   anprojmean=yp;    anprojmean=yp;
   yp2=modf((yp1*12),&yp);    yp2=modf((yp1*12),&yp);
Line 2943  fprintf(ficres,"popforecast=%d popfile=% Line 2943  fprintf(ficres,"popforecast=%d popfile=%
   jprojmean=yp;    jprojmean=yp;
   if(jprojmean==0) jprojmean=1;    if(jprojmean==0) jprojmean=1;
   if(mprojmean==0) jprojmean=1;    if(mprojmean==0) jprojmean=1;
    
   fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);    fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);
    
   if (popforecast==1) {    if (popforecast==1) {
     if((ficpop=fopen(popfile,"r"))==NULL)    {      if((ficpop=fopen(popfile,"r"))==NULL)    {
       printf("Problem with population file : %s\n",popfile);goto end;        printf("Problem with population file : %s\n",popfile);goto end;
Line 2953  fprintf(ficres,"popforecast=%d popfile=% Line 2953  fprintf(ficres,"popforecast=%d popfile=%
     popage=ivector(0,AGESUP);      popage=ivector(0,AGESUP);
     popeffectif=vector(0,AGESUP);      popeffectif=vector(0,AGESUP);
     popcount=vector(0,AGESUP);      popcount=vector(0,AGESUP);
      
     i=1;        i=1;  
     while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)      while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)
       {        {
Line 2963  fprintf(ficres,"popforecast=%d popfile=% Line 2963  fprintf(ficres,"popforecast=%d popfile=%
         
     for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];      for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
   }    }
    
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;        k=k+1;
Line 2976  fprintf(ficres,"popforecast=%d popfile=% Line 2976  fprintf(ficres,"popforecast=%d popfile=%
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);        for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
       if (popforecast==1)  fprintf(ficresf," [Population]");        if (popforecast==1)  fprintf(ficresf," [Population]");
             
      for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {        for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {
         fprintf(ficresf,"\n");          fprintf(ficresf,"\n");
         fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);            fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
   
         for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){  
         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);  
         nhstepm = nhstepm/hstepm;  
                 
         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
         oldm=oldms;savm=savms;            nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);              nhstepm = nhstepm/hstepm;
            
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
                 
         for (h=0; h<=nhstepm; h++){            for (h=0; h<=nhstepm; h++){
           if (h==(int) (calagedate+YEARM*cpt)) {              if (h==(int) (calagedate+YEARM*cpt)) {
             fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);                fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
           }              }
           for(j=1; j<=nlstate+ndeath;j++) {              for(j=1; j<=nlstate+ndeath;j++) {
             kk1=0.;kk2=0;                kk1=0.;kk2=0;
             for(i=1; i<=nlstate;i++) {                        for(i=1; i<=nlstate;i++) {              
               if (mobilav==1)                  if (mobilav==1)
                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];                    kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
               else {                  else {
                 kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];                    kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
                                  }
                   if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
               }                }
                 if (h==(int)(calagedate+12*cpt)){
               if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];                  fprintf(ficresf," %.3f", kk1);
             }  
            
             if (h==(int)(calagedate+12*cpt)){  
               fprintf(ficresf," %.3f", kk1);  
                             
               if (popforecast==1) fprintf(ficresf," [%.f]", kk2);                  if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
                 }
             }              }
           }            }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         }          }
         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);        }
         }  
      }  
     }      }
   }    }
          
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   if (popforecast==1) {    if (popforecast==1) {
Line 3024  fprintf(ficres,"popforecast=%d popfile=% Line 3022  fprintf(ficres,"popforecast=%d popfile=%
     free_vector(popeffectif,0,AGESUP);      free_vector(popeffectif,0,AGESUP);
     free_vector(popcount,0,AGESUP);      free_vector(popcount,0,AGESUP);
   }    }
    
   free_imatrix(s,1,maxwav+1,1,n);    free_imatrix(s,1,maxwav+1,1,n);
   free_vector(weight,1,n);    free_vector(weight,1,n);
   fclose(ficresf);    fclose(ficresf);
   }/* End forecasting */    }
   
     /* End forecasting */
   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("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);
   
   }    }
    
     
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
Line 3083  fprintf(ficres,"popforecast=%d popfile=% Line 3083  fprintf(ficres,"popforecast=%d popfile=%
       evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);          evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);  
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
       varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);         varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
           
   
    
       fprintf(ficrest,"#Total LEs with variances: e.. (std) ");        fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
       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;        hf=1;
       if (stepm >= YEARM) hf=stepm/YEARM;        if (stepm >= YEARM) hf=stepm/YEARM;
       epj=vector(1,nlstate+1);        epj=vector(1,nlstate+1);
Line 3127  fprintf(ficres,"popforecast=%d popfile=% Line 3129  fprintf(ficres,"popforecast=%d popfile=%
   fclose(ficpar);    fclose(ficpar);
   free_vector(epj,1,nlstate+1);    free_vector(epj,1,nlstate+1);
   /*  scanf("%d ",i); */    /*  scanf("%d ",i); */
    
   /*------- Variance limit prevalence------*/      /*------- Variance limit prevalence------*/  
   
 strcpy(fileresvpl,"vpl");  strcpy(fileresvpl,"vpl");
Line 3169  strcpy(fileresvpl,"vpl"); Line 3171  strcpy(fileresvpl,"vpl");
     
   free_matrix(matcov,1,npar,1,npar);    free_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);    free_vector(delti,1,npar);
      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(erreur >0)    if(erreur >0)
Removed from v.1.25  
changed lines
  Added in v.1.26

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