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

Diff for /imach/src/imach.c between versions 1.13 and 1.14

version 1.13, 2002/02/20 17:02:08 version 1.14, 2002/02/20 17:05:44
Line 73  int npar=NPARMAX; Line 73  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, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
   int popbased=0, fprev,lprev;
   
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */  int maxwav; /* Maxim number of waves */
Line 1191  void  freqsummary(char fileres[], int ag Line 1192  void  freqsummary(char fileres[], int ag
                bool=0;                 bool=0;
          }           }
           if (bool==1) {            if (bool==1) {
            for(m=firstpass; m<=lastpass-1; m++){             for(m=fprev; m<=lprev; m++){
              if(agev[m][i]==0) agev[m][i]=agemax+1;               if(agev[m][i]==0) agev[m][i]=agemax+1;
              if(agev[m][i]==1) agev[m][i]=agemax+2;               if(agev[m][i]==1) agev[m][i]=agemax+2;
              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];               freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
Line 1215  void  freqsummary(char fileres[], int ag Line 1216  void  freqsummary(char fileres[], int ag
       printf("Age %d", i);        printf("Age %d", i);
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)        for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }      }
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pos=0; m <=0 ; m++)        for(m=-1, pos=0; m <=0 ; m++)
Line 1225  void  freqsummary(char fileres[], int ag Line 1226  void  freqsummary(char fileres[], int ag
       else        else
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);          printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     }      }
     for(jk=1; jk <=nlstate ; jk++){  
       for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)       for(jk=1; jk <=nlstate ; jk++){
         for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }       }
   
     for(jk=1,pos=0; jk <=nlstate ; jk++)      for(jk=1,pos=0; jk <=nlstate ; jk++)
       pos += pp[jk];        pos += pp[jk];
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
Line 1492  void varevsij(char fileres[], double *** Line 1495  void varevsij(char fileres[], double ***
       }        }
       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) {
           for(i=1; i<=nlstate;i++)
             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 1503  void varevsij(char fileres[], double *** Line 1512  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) {
           for(i=1; i<=nlstate;i++)
             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, gm[h][j]=0.;i<=nlstate;i++)            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
             gm[h][j] += prlim[i][i]*p3mat[i][j][h];              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
         }          }
       }        }
   
       for(j=1; j<= nlstate; j++)        for(j=1; j<= nlstate; j++)
         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];
Line 1773  int main() Line 1789  int main()
   int c,  h , cpt,l;    int c,  h , cpt,l;
   int ju,jl, mi;    int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
      int mobilav=0, fprevfore=1, lprevfore=1;
   int hstepm, nhstepm;    int hstepm, nhstepm;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
   double **prlim;    double **prlim;
Line 1852  split(pathtot, path,optionfile); Line 1869  split(pathtot, path,optionfile);
   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 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 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 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 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 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);
   while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
    while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);
    
   covar=matrix(0,NCOVMAX,1,n);    covar=matrix(0,NCOVMAX,1,n);
   cptcovn=0;    cptcovn=0;
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;    if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
Line 2014  split(pathtot, path,optionfile); Line 2049  split(pathtot, path,optionfile);
     if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;      if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;      if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;      if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
   }      }
   for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/      for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
Line 2525  Interval (in months) between two waves: Line 2560  Interval (in months) between two waves:
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>          - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>          - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>          - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);          - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
           - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
   <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
   
  fprintf(fichtm," <li>Graphs</li><p>");   fprintf(fichtm," <li>Graphs</li><p>");
   
Line 2663  fclose(fichtm); Line 2700  fclose(fichtm);
   
   fclose(ficrespij);    fclose(ficrespij);
   
   exit(0);  
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   
   strcpy(fileresf,"f");    strcpy(fileresf,"f");
Line 2682  fclose(fichtm); Line 2718  fclose(fichtm);
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
   mobaverage= ma3x(1,130 ,1,8, 1,8);    if (mobilav==1) {
   for (agedeb=bage+3; agedeb<=fage-2; agedeb++)      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     for (i=1; i<=nlstate;i++)      for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
       for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)        for (i=1; i<=nlstate;i++)
         mobaverage[(int)agedeb][i][cptcod]=0.;          for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
              mobaverage[(int)agedeb][i][cptcod]=0.;
   for (agedeb=bage+4; agedeb<=fage; agedeb++){     
     for (i=1; i<=nlstate;i++){      for (agedeb=bage+4; agedeb<=fage; agedeb++){
       for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){        for (i=1; i<=nlstate;i++){
         for (cpt=0;cpt<=4;cpt++){          for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
           mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];            for (cpt=0;cpt<=4;cpt++){
               mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
           }            }
           mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;            mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
           }
       }        }
     }      }  
   }    }
   
 /* if (cptcod==2) printf("m=%f p=%f %d age=%d ",mobaverage[(int)agedeb-2][i][cptcod],probs[(int)agedeb-cpt][i][cptcod],cpt,(int)agedeb-2);*/  
   
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
   if (stepm<=24) stepsize=2;    if (stepm<=24) stepsize=2;
   
Line 2725  fclose(fichtm); Line 2760  fclose(fichtm);
   
       for (agedeb=fage; agedeb>=bage; agedeb--){        for (agedeb=fage; agedeb>=bage; agedeb--){
         fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);          fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
          if (mobilav==1) {
         for(j=1; j<=nlstate;j++)          for(j=1; j<=nlstate;j++)
           fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]);            fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);
       }          }
           else {
             for(j=1; j<=nlstate;j++)
             fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]);
           }
          
       for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");        for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");
         }
       for (cpt=1; cpt<=8;cpt++)          for (cpt=1; cpt<=NCOVMAX;cpt++)  
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */        for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
                 
         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */          nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
Line 2749  fclose(fichtm); Line 2790  fclose(fichtm);
           for(j=1; j<=nlstate+ndeath;j++) {            for(j=1; j<=nlstate+ndeath;j++) {
             kk1=0.;              kk1=0.;
             for(i=1; i<=nlstate;i++) {                      for(i=1; i<=nlstate;i++) {        
               /*   kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];*/                if (mobilav==1)
                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];                kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
             }                else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];
                          }    
             if (h*hstepm/YEARM*stepm==cpt)              if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);
               fprintf(ficresf," %.5f ", kk1);  
           }  
           }            }
         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);
   fclose(ficresf);    fclose(ficresf);
   
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");    strcpy(filerest,"t");
Line 2822  fclose(fichtm); Line 2862  fclose(fichtm);
       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);
           if (popbased==1) {
             for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];
           }
          
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %.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++) {
Removed from v.1.13  
changed lines
  Added in v.1.14

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