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

Diff for /imach/src/imach.c between versions 1.59 and 1.63

version 1.59, 2002/11/18 23:01:13 version 1.63, 2002/11/20 17:35:59
Line 928  double func( double *x) Line 928  double func( double *x)
   cov[1]=1.;    cov[1]=1.;
   
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){  
     for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];    if(mle==1){
     for(mi=1; mi<= wav[i]-1; mi++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (ii=1;ii<=nlstate+ndeath;ii++)        for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for (j=1;j<=nlstate+ndeath;j++){        for(mi=1; mi<= wav[i]-1; mi++){
           oldm[ii][j]=(ii==j ? 1.0 : 0.0);          for (ii=1;ii<=nlstate+ndeath;ii++)
           savm[ii][j]=(ii==j ? 1.0 : 0.0);            for (j=1;j<=nlstate+ndeath;j++){
         }              oldm[ii][j]=(ii==j ? 1.0 : 0.0);
       for(d=0; d<dh[mi][i]; d++){              savm[ii][j]=(ii==j ? 1.0 : 0.0);
         newm=savm;            }
         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;          for(d=0; d<dh[mi][i]; d++){
         for (kk=1; kk<=cptcovage;kk++) {            newm=savm;
           cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];            cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
         }            for (kk=1; kk<=cptcovage;kk++) {
                       cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            }
                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));          
         savm=oldm;            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
         oldm=newm;                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             savm=oldm;
             oldm=newm;
                   
                   
       } /* end mult */          } /* end mult */
               
       /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */          /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
       /* But now since version 0.9 we anticipate for bias and large stepm.          /* But now since version 0.9 we anticipate for bias and large stepm.
        * If stepm is larger than one month (smallest stepm) and if the exact delay            * If stepm is larger than one month (smallest stepm) and if the exact delay 
        * (in months) between two waves is not a multiple of stepm, we rounded to            * (in months) between two waves is not a multiple of stepm, we rounded to 
        * the nearest (and in case of equal distance, to the lowest) interval but now           * the nearest (and in case of equal distance, to the lowest) interval but now
        * we keep into memory the bias bh[mi][i] and also the previous matrix product           * we keep into memory the bias bh[mi][i] and also the previous matrix product
        * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the           * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
        * probability in order to take into account the bias as a fraction of the way           * probability in order to take into account the bias as a fraction of the way
        * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies           * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
        * -stepm/2 to stepm/2 .           * -stepm/2 to stepm/2 .
        * For stepm=1 the results are the same as for previous versions of Imach.           * For stepm=1 the results are the same as for previous versions of Imach.
        * For stepm > 1 the results are less biased than in previous versions.            * For stepm > 1 the results are less biased than in previous versions. 
        */           */
       s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
       s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
       bbh=(double)bh[mi][i]/(double)stepm;          bbh=(double)bh[mi][i]/(double)stepm;
       lli= (savm[s1][s2]>(double)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]));*/          /*      lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
       /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/  
       /*if(lli ==000.0)*/          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 */
       /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */  
       ipmx +=1;          /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/
       sw += weight[i];          /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          /*if(lli ==000.0)*/
     } /* end of wave */          /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
   } /* end of individual */          ipmx +=1;
           sw += weight[i];
           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
         } /* end of wave */
       } /* end of individual */
     }  else{ 
       for (i=1,ipmx=0, sw=0.; i<=imx; i++){
         for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         for(mi=1; mi<= wav[i]-1; mi++){
           for (ii=1;ii<=nlstate+ndeath;ii++)
             for (j=1;j<=nlstate+ndeath;j++){
               oldm[ii][j]=(ii==j ? 1.0 : 0.0);
               savm[ii][j]=(ii==j ? 1.0 : 0.0);
             }
           for(d=0; d<dh[mi][i]; d++){
             newm=savm;
             cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
             for (kk=1; kk<=cptcovage;kk++) {
               cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
             }
           
             out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                          1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
             savm=oldm;
             oldm=newm;
           } /* end mult */
         
           lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
           ipmx +=1;
           sw += weight[i];
           ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
         } /* end of wave */
       } /* end of individual */
     } /* End of if */
   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 */
Line 3122  populforecast(char fileres[], double anp Line 3156  populforecast(char fileres[], double anp
   
 int main(int argc, char *argv[])  int main(int argc, char *argv[])
 {  {
     int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   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 3160  int main(int argc, char *argv[]) Line 3194  int main(int argc, char *argv[])
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
   /*int *movingaverage; */  
   
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
Line 3652  int main(int argc, char *argv[]) Line 3685  int main(int argc, char *argv[])
   /* 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'. */
   
       pmmij= 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 */
       savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
       oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
           
         
   /* For Powell, parameters are in a vector p[] starting at p[1]    /* For Powell, parameters are in a vector p[] starting at p[1]
      so we point p on param[1][1] so that p[1] maps on param[1][1][1] */       so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
   
   if(mle==1){    if(mle>=1){ /* Could be 1 or 2 */
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);      mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   }    }
           
Line 3884  Interval (in months) between two waves: Line 3922  Interval (in months) between two waves:
   fprintf(ficrespl,"\n");    fprintf(ficrespl,"\n");
       
   prlim=matrix(1,nlstate,1,nlstate);    prlim=matrix(1,nlstate,1,nlstate);
   pmmij= 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 */  
   savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
   oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */  
   
   agebase=ageminpar;    agebase=ageminpar;
   agelim=agemaxpar;    agelim=agemaxpar;
Removed from v.1.59  
changed lines
  Added in v.1.63

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