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

Diff for /imach/src/imach.c between versions 1.67 and 1.68

version 1.67, 2003/01/28 17:41:19 version 1.68, 2003/02/04 12:40:59
Line 919  double func( double *x) Line 919  double func( double *x)
   double sw; /* Sum of weights */    double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */    double lli; /* Individual log likelihood */
   int s1, s2;    int s1, s2;
   double bbh;    double bbh, survp;
   long ipmx;    long ipmx;
   /*extern weight */    /*extern weight */
   /* We are differentiating ll according to initial status */    /* We are differentiating ll according to initial status */
Line 972  double func( double *x) Line 972  double func( double *x)
          * is higher than the multiple of stepm and negative otherwise.           * is higher than the multiple of stepm and negative otherwise.
          */           */
         /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/          /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
           /* if s2=-2 lli=out[1][1]+out[1][2];*/
           if (s2==-2) {
             for (j=1,survp=0. ; j<=nlstate; j++) 
               survp += out[s1][j];
             lli= survp;
           }
   
           else
         lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));  /* linear interpolation */          lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));  /* linear interpolation */
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/          /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
         /*if(lli ==000.0)*/          /*if(lli ==000.0)*/
Line 1416  void  freqsummary(char fileres[], int ag Line 1424  void  freqsummary(char fileres[], int ag
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);      exit(0);
   }    }
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
       
   j=cptcoveff;    j=cptcoveff;
Line 1531  void  freqsummary(char fileres[], int ag Line 1539  void  freqsummary(char fileres[], int ag
           }            }
         }          }
                   
         for(jk=-1; jk <=nlstate+ndeath; jk++)          for(jk=-2; jk <=nlstate+ndeath; jk++)
           for(m=-1; m <=nlstate+ndeath; m++)            for(m=-2; m <=nlstate+ndeath; m++)
             if(freq[jk][m][i] !=0 ) {              if(freq[jk][m][i] !=0 ) {
             if(first==1)              if(first==1)
               printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);                printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
Line 1549  void  freqsummary(char fileres[], int ag Line 1557  void  freqsummary(char fileres[], int ag
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
       
   /* End of Freq */    /* End of Freq */
Line 1566  void prevalence(int agemin, float agemax Line 1574  void prevalence(int agemin, float agemax
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
       
   freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq=ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
       
   j=cptcoveff;    j=cptcoveff;
Line 1576  void prevalence(int agemin, float agemax Line 1584  void prevalence(int agemin, float agemax
     for(i1=1; i1<=ncodemax[k1];i1++){      for(i1=1; i1<=ncodemax[k1];i1++){
       j1++;        j1++;
               
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-2; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-2; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
             freq[i][jk][m]=0;              freq[i][jk][m]=0;
             
Line 1634  void prevalence(int agemin, float agemax Line 1642  void prevalence(int agemin, float agemax
   } /* end k1 */    } /* end k1 */
   
       
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
       
 }  /* End of Freq */  }  /* End of Freq */
Line 1664  void  concatwav(int wav[], int **dh, int Line 1672  void  concatwav(int wav[], int **dh, int
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
     while(s[m][i] <= nlstate){      while(s[m][i] <= nlstate){
       if(s[m][i]>=1)        if(s[m][i]>=1 || s[m][i]==-2)
         mw[++mi][i]=m;          mw[++mi][i]=m;
       if(m >=lastpass)        if(m >=lastpass)
         break;          break;
Line 1704  void  concatwav(int wav[], int **dh, int Line 1712  void  concatwav(int wav[], int **dh, int
           if (j <= jmin) jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           /*if (j<0) printf("j=%d num=%d \n",j,i); */            /*if (j<0) printf("j=%d num=%d \n",j,i); */
             /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
           }            }
         }          }
         else{          else{
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));            j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
             /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
           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;
Line 1742  void  concatwav(int wav[], int **dh, int Line 1752  void  concatwav(int wav[], int **dh, int
             bh[mi][i]=ju; /* At least one step */              bh[mi][i]=ju; /* At least one step */
             printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);              printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
           }            }
           if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);  
         }          }
       } /* end if mle */        } /* end if mle */
     } /* end wave */      } /* end wave */
Line 2136  void varevsij(char optionfilefiname[], d Line 2145  void varevsij(char optionfilefiname[], d
          as a weighted average of prlim.           as a weighted average of prlim.
       */        */
       for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){        for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){
         for(i=1; i<= nlstate; i++)          for(i=1,gpp[j]=0.; i<= nlstate; i++)
           gpp[j] += prlim[i][i]*p3mat[i][j][1];            gpp[j] += prlim[i][i]*p3mat[i][j][1];
       }            }    
       /* end probability of death */        /* end probability of death */
Line 2167  void varevsij(char optionfilefiname[], d Line 2176  void varevsij(char optionfilefiname[], d
          as a weighted average of prlim.           as a weighted average of prlim.
       */        */
       for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){        for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
         for(i=1; i<= nlstate; i++)          for(i=1,gmp[j]=0.; i<= nlstate; i++)
           gmp[j] += prlim[i][i]*p3mat[i][j][1];           gmp[j] += prlim[i][i]*p3mat[i][j][1];
       }            }    
       /* end probability of death */        /* end probability of death */
   
Line 2176  void varevsij(char optionfilefiname[], d Line 2185  void varevsij(char optionfilefiname[], d
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
         }          }
   
       for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */        for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
         gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];          gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
       }        }
Line 2190  void varevsij(char optionfilefiname[], d Line 2200  void varevsij(char optionfilefiname[], d
           trgradg[h][j][theta]=gradg[h][theta][j];            trgradg[h][j][theta]=gradg[h][theta][j];
   
     for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
       for(theta=1; theta <=npar; theta++)        for(theta=1; theta <=npar; theta++) {
         trgradgp[j][theta]=gradgp[theta][j];          trgradgp[j][theta]=gradgp[theta][j];
         }
   
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
     for(i=1;i<=nlstate;i++)      for(i=1;i<=nlstate;i++)
Line 2211  void varevsij(char optionfilefiname[], d Line 2222  void varevsij(char optionfilefiname[], d
     /* pptj */      /* pptj */
     matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);      matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
     matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);      matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
     for(j=nlstate+1;j<=nlstate+ndeath;j++)   
       for(i=nlstate+1;i<=nlstate+ndeath;i++)     for(j=nlstate+1;j<=nlstate+ndeath;j++)
        for(i=nlstate+1;i<=nlstate+ndeath;i++){
         varppt[j][i]=doldmp[j][i];          varppt[j][i]=doldmp[j][i];
        }
   
     /* end ppptj */      /* end ppptj */
     /*  x centered again */      /*  x centered again */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);        hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);  
Line 2233  void varevsij(char optionfilefiname[], d Line 2247  void varevsij(char optionfilefiname[], d
        computed over hstepm (estepm) matrices product = hstepm*stepm months)          computed over hstepm (estepm) matrices product = hstepm*stepm months) 
        as a weighted average of prlim.         as a weighted average of prlim.
     */      */
     for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){      for(j=nlstate+1;j<=nlstate+ndeath;j++){
       for(i=1; i<= nlstate; i++)        for(i=1,gmp[j]=0.;i<= nlstate; i++) 
         gmp[j] += prlim[i][i]*p3mat[i][j][1];           gmp[j] += prlim[i][i]*p3mat[i][j][1]; 
     }          }    
     /* end probability of death */      /* end probability of death */
Line 3704  int main(int argc, char *argv[]) Line 3718  int main(int argc, char *argv[])
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);      agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
     for(m=1; (m<= maxwav); m++){      for(m=1; (m<= maxwav); m++){
       if(s[m][i] >0){        if(s[m][i] >0 || s[m][i]==-2){
         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(moisdc[i]!=99 && 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 (andc[i]!=9999){
Line 3770  int main(int argc, char *argv[]) Line 3783  int main(int argc, char *argv[])
   dh=imatrix(1,lastpass-firstpass+1,1,imx);    dh=imatrix(1,lastpass-firstpass+1,1,imx);
   bh=imatrix(1,lastpass-firstpass+1,1,imx);    bh=imatrix(1,lastpass-firstpass+1,1,imx);
   mw=imatrix(1,lastpass-firstpass+1,1,imx);    mw=imatrix(1,lastpass-firstpass+1,1,imx);
        
   
   /* Concatenates waves */    /* Concatenates waves */
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   
Line 3979  int main(int argc, char *argv[]) Line 3993  int main(int argc, char *argv[])
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);    fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);    fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   
   
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
   /*------------ gnuplot -------------*/    /*------------ gnuplot -------------*/
Line 4141  Interval (in months) between two waves: Line 4156  Interval (in months) between two waves:
   
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   if((stepm == 1) && (strcmp(model,".")==0)){     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);
   }     } 
Line 4149  Interval (in months) between two waves: Line 4164  Interval (in months) between two waves:
     erreur=108;      erreur=108;
     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);      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);
     fprintf(ficlog,"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);      fprintf(ficlog,"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);
   }      }
       
   
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
Line 4173  Interval (in months) between two waves: Line 4188  Interval (in months) between two waves:
   printf("Computing Health Expectancies: result on file '%s' \n", filerese);    printf("Computing Health Expectancies: result on file '%s' \n", filerese);
   fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);    fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
   
   
   strcpy(fileresv,"v");    strcpy(fileresv,"v");
   strcat(fileresv,fileres);    strcat(fileresv,fileres);
   if((ficresvij=fopen(fileresv,"w"))==NULL) {    if((ficresvij=fopen(fileresv,"w"))==NULL) {
Removed from v.1.67  
changed lines
  Added in v.1.68

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