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

Diff for /imach/src/imach.c between versions 1.31 and 1.34

version 1.31, 2002/03/10 13:43:02 version 1.34, 2002/03/13 17:19:16
Line 82  int cptcovn, cptcovage=0, cptcoveff=0,cp Line 82  int cptcovn, cptcovage=0, cptcoveff=0,cp
 int npar=NPARMAX;  int npar=NPARMAX;
 int nlstate=2; /* Number of live states */  int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */  int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */  int ncovmodel, ncovcol;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
 int popbased=0;  int popbased=0;
   
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
Line 1513  void tricode(int *Tvar, int **nbcode, in Line 1513  void tricode(int *Tvar, int **nbcode, in
   
  ij=1;   ij=1;
  for (i=1; i<=10; i++) {   for (i=1; i<=10; i++) {
    if((Ndum[i]!=0) && (i<=ncov)){     if((Ndum[i]!=0) && (i<=ncovcol)){
      Tvaraff[ij]=i;       Tvaraff[ij]=i;
      ij++;       ij++;
    }     }
Line 1540  void evsij(char fileres[], double ***eij Line 1540  void evsij(char fileres[], double ***eij
   
   k=1;             /* For example stepm=6 months */    k=1;             /* For example stepm=6 months */
   hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */    hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */
   hstepm=1;   /* or (b) We decided to compute the life expectancy with the smallest unit */    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.    /* 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       nhstepm is the number of hstepm from age to agelim
      nstepm is the number of stepm from age to agelin.       nstepm is the number of stepm from age to agelin.
      Look at hpijx to understand the reason of that which relies in memory size       Look at hpijx to understand the reason of that which relies in memory size
      and note for a fixed period like k years */       and note for a fixed period like k years */
   /* We decided (b) to get a life expectancy respecting the most precise curvature of the    /* We decided (b) to get a life expectancy respecting the most precise curvature of the
      survival function given par stepm (the optimization length). Unfortunately it       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       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       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.       results. So we changed our mind and took the option of the best precision.
Line 1559  void evsij(char fileres[], double ***eij Line 1559  void evsij(char fileres[], double ***eij
     /* nhstepm age range expressed in number of stepm */      /* nhstepm age range expressed in number of stepm */
     nstepm=(int) rint((agelim-age)*YEARM/stepm);      nstepm=(int) rint((agelim-age)*YEARM/stepm);
     /* Typically if 20 years nstepm = 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 = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=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/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++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
Line 1903  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.71a </font> <hr size=\"2\" color=\"#EC5E5E\">   fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.8 </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 2466  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.71a, March 2002, INED-EUROREVES ";    char version[80]="Imach version 0.8, 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 2527  int main(int argc, char *argv[]) Line 2527  int main(int argc, char *argv[])
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf 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);    fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
   printf("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);    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"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(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
 while((c=getc(ficpar))=='#' && c!= EOF){  while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
Line 2686  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2686  while((c=getc(ficpar))=='#' && c!= EOF){
         cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);          cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
   
         cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);          cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
         for (j=ncov;j>=1;j--){          for (j=ncovcol;j>=1;j--){
           cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);            cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
         }          }
         num[i]=atol(stra);          num[i]=atol(stra);
Line 2755  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2755  while((c=getc(ficpar))=='#' && c!= EOF){
         }          }
         else {          else {
           cutv(strb,stre,strc,'V');            cutv(strb,stre,strc,'V');
           Tvar[i]=ncov+k1;            Tvar[i]=ncovcol+k1;
           cutv(strb,strc,strd,'V');            cutv(strb,strc,strd,'V');
           Tprod[k1]=i;            Tprod[k1]=i;
           Tvard[k1][1]=atoi(strc);            Tvard[k1][1]=atoi(strc);
Line 2763  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2763  while((c=getc(ficpar))=='#' && c!= EOF){
           Tvar[cptcovn+k2]=Tvard[k1][1];            Tvar[cptcovn+k2]=Tvard[k1][1];
           Tvar[cptcovn+k2+1]=Tvard[k1][2];            Tvar[cptcovn+k2+1]=Tvard[k1][2];
           for (k=1; k<=lastobs;k++)            for (k=1; k<=lastobs;k++)
             covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];              covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
           k1++;            k1++;
           k2=k2+2;            k2=k2+2;
         }          }
Line 2917  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= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, weightopt,model);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
     
   
    jk=1;     jk=1;
    fprintf(ficres,"# Parameters\n");     fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
    printf("# Parameters\n");     printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
    for(i=1,jk=1; i <=nlstate; i++){     for(i=1,jk=1; i <=nlstate; i++){
      for(k=1; k <=(nlstate+ndeath); k++){       for(k=1; k <=(nlstate+ndeath); k++){
        if (k != i)         if (k != i)
Line 2944  printf("Total number of individuals= %d, Line 2944  printf("Total number of individuals= %d,
     ftolhess=ftol; /* Usually correct */      ftolhess=ftol; /* Usually correct */
     hesscov(matcov, p, npar, delti, ftolhess, func);      hesscov(matcov, p, npar, delti, ftolhess, func);
  }   }
     fprintf(ficres,"# Scales\n");      fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
     printf("# Scales\n");      printf("# Scales (for hessian or gradient estimation)\n");
      for(i=1,jk=1; i <=nlstate; i++){       for(i=1,jk=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath; j++){        for(j=1; j <=nlstate+ndeath; j++){
         if (j!=i) {          if (j!=i) {
Line 2963  printf("Total number of individuals= %d, Line 2963  printf("Total number of individuals= %d,
      }       }
         
     k=1;      k=1;
     fprintf(ficres,"# Covariance\n");      fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
     printf("# Covariance\n");      printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
     for(i=1;i<=npar;i++){      for(i=1;i<=npar;i++){
       /*  if (k>nlstate) k=1;        /*  if (k>nlstate) k=1;
       i1=(i-1)/(ncovmodel*nlstate)+1;        i1=(i-1)/(ncovmodel*nlstate)+1;
Line 3175  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3175  while((c=getc(ficpar))=='#' && c!= EOF){
   
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   if((stepm == 1) && (model==".")){    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 3243  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 3253  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");
       }        }
Line 3323  if (popforecast==1) populforecast(filere Line 3321  if (popforecast==1) populforecast(filere
   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)
     printf("End of Imach with error %d\n",erreur);      printf("End of Imach with error or warning %d\n",erreur);
   else   printf("End of Imach\n");    else   printf("End of Imach\n");
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
     
Removed from v.1.31  
changed lines
  Added in v.1.34

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