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

Diff for /imach/src/imach.c between versions 1.43 and 1.46

version 1.43, 2002/05/24 13:00:54 version 1.46, 2002/05/30 17:44:35
Line 96  double jmean; /* Mean space between 2 wa Line 96  double jmean; /* Mean space between 2 wa
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficgp,*ficresprob,*ficpop;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficreseij;  FILE *ficreseij;
   char filerese[FILENAMELENGTH];    char filerese[FILENAMELENGTH];
  FILE  *ficresvij;   FILE  *ficresvij;
Line 1641  void evsij(char fileres[], double ***eij Line 1641  void evsij(char fileres[], double ***eij
           }            }
         }          }
       }        }
        
      
   
       for(j=1; j<= nlstate*2; j++)        for(j=1; j<= nlstate*2; j++)
         for(h=0; h<=nhstepm-1; h++){          for(h=0; h<=nhstepm-1; 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];
         }          }
   
      }       }
         
 /* End theta */  /* End theta */
Line 1658  void evsij(char fileres[], double ***eij Line 1654  void evsij(char fileres[], double ***eij
      for(h=0; h<=nhstepm-1; h++)       for(h=0; h<=nhstepm-1; h++)
       for(j=1; j<=nlstate*2;j++)        for(j=1; j<=nlstate*2;j++)
         for(theta=1; theta <=npar; theta++)          for(theta=1; theta <=npar; theta++)
         trgradg[h][j][theta]=gradg[h][theta][j];            trgradg[h][j][theta]=gradg[h][theta][j];
        
   
      for(i=1;i<=nlstate*2;i++)       for(i=1;i<=nlstate*2;i++)
       for(j=1;j<=nlstate*2;j++)        for(j=1;j<=nlstate*2;j++)
         varhe[i][j][(int)age] =0.;          varhe[i][j][(int)age] =0.;
   
      printf("%d|",(int)age);fflush(stdout);       printf("%d|",(int)age);fflush(stdout);
     for(h=0;h<=nhstepm-1;h++){       for(h=0;h<=nhstepm-1;h++){
       for(k=0;k<=nhstepm-1;k++){        for(k=0;k<=nhstepm-1;k++){
         matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);          matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
Line 1931  void varprevlim(char fileres[], double * Line 1927  void varprevlim(char fileres[], double *
 void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)  void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
 {  {
   int i, j, i1, k1, j1, z1;    int i, j, i1, k1, j1, z1;
   int k=0, cptcode;    int k=0,l, cptcode;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double *xp;    double *xp;
   double *gp, *gm;    double *gp, *gm;
Line 1939  void varprob(char fileres[], double **ma Line 1935  void varprob(char fileres[], double **ma
   double age,agelim, cov[NCOVMAX];    double age,agelim, cov[NCOVMAX];
   int theta;    int theta;
   char fileresprob[FILENAMELENGTH];    char fileresprob[FILENAMELENGTH];
     char fileresprobcov[FILENAMELENGTH];
     char fileresprobcor[FILENAMELENGTH];
   
   strcpy(fileresprob,"prob");    strcpy(fileresprob,"prob");
   strcat(fileresprob,fileres);    strcat(fileresprob,fileres);
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {    if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);      printf("Problem with resultfile: %s\n", fileresprob);
   }    }
     strcpy(fileresprobcov,"probcov");
     strcat(fileresprobcov,fileres);
     if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprobcov);
     }
     strcpy(fileresprobcor,"probcor");
     strcat(fileresprobcor,fileres);
     if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprobcor);
     }
   printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);    printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
     printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
     printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
     
 fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");    fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");    fprintf(ficresprob,"# Age");
     fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
     fprintf(ficresprobcov,"# Age");
     fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
     fprintf(ficresprobcov,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
     for(j=1; j<=(nlstate+ndeath);j++)      for(j=1; j<=(nlstate+ndeath);j++){
       fprintf(ficresprob," p%1d-%1d (SE)",i,j);        fprintf(ficresprob," p%1d-%1d (SE)",i,j);
         fprintf(ficresprobcov," p%1d-%1d ",i,j);
         fprintf(ficresprobcor," p%1d-%1d ",i,j);
       }  
   fprintf(ficresprob,"\n");    fprintf(ficresprob,"\n");
     fprintf(ficresprobcov,"\n");
     fprintf(ficresprobcor,"\n");
   xp=vector(1,npar);    xp=vector(1,npar);
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);    dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));    doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
Line 1971  fprintf(ficresprob,"#One-step probabilit Line 1986  fprintf(ficresprob,"#One-step probabilit
   
     if  (cptcovn>0) {      if  (cptcovn>0) {
       fprintf(ficresprob, "\n#********** Variable ");        fprintf(ficresprob, "\n#********** Variable ");
         fprintf(ficresprobcov, "\n#********** Variable ");
         fprintf(ficresprobcor, "\n#********** Variable ");
       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
       fprintf(ficresprob, "**********\n#");        fprintf(ficresprob, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficresprobcov, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficresprobcor, "**********\n#");
     }      }
         
       for (age=bage; age<=fage; age ++){        for (age=bage; age<=fage; age ++){
         cov[2]=age;          cov[2]=age;
         for (k=1; k<=cptcovn;k++) {          for (k=1; k<=cptcovn;k++) {
           cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];            cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
            
         }          }
         for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
         for (k=1; k<=cptcovprod;k++)          for (k=1; k<=cptcovprod;k++)
Line 2037  fprintf(ficresprob,"#One-step probabilit Line 2057  fprintf(ficresprob,"#One-step probabilit
           }            }
         }          }
             
      /*printf("\n%d ",(int)age);          /*printf("\n%d ",(int)age);
      for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){       for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
        printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));         printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
      }*/       }*/
   
         fprintf(ficresprob,"\n%d ",(int)age);          fprintf(ficresprob,"\n%d ",(int)age);
           fprintf(ficresprobcov,"\n%d ",(int)age);
           fprintf(ficresprobcor,"\n%d ",(int)age);
   
         for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)          for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
           fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));            fprintf(ficresprob,"%12.3e (%12.3e) ",gm[i],sqrt(doldm[i][j]));
            for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
             fprintf(ficresprobcov,"%12.3e ",gm[i]);
             fprintf(ficresprobcor,"%12.3e ",gm[i]);
           }
           i=0;
           for (k=1; k<=(nlstate);k++){
             for (l=1; l<=(nlstate+ndeath);l++){
               i=i++;
               fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
               fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
               for (j=1; j<=i;j++){
                 fprintf(ficresprobcov," %12.3e",doldm[i][j]);
                 fprintf(ficresprobcor," %12.3e",doldm[i][j]/sqrt(doldm[i][i])/sqrt(doldm[j][j]));
               }
             }
           }
       }        }
     }      }
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));      free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
Line 2056  fprintf(ficresprob,"#One-step probabilit Line 2093  fprintf(ficresprob,"#One-step probabilit
   }    }
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   fclose(ficresprob);    fclose(ficresprob);
      fclose(ficresprobcov);
     fclose(ficresprobcor);
 }  }
   
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
                   int imx,int jmin, int jmax, double jmeanint,char optionfile[], \                    int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
                   char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\                    char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
                   char version[], int popforecast, int estepm ,/* \ */                    char version[], int popforecast, int estepm ,\
                   double jprev1, double mprev1,double anprev1, \                    double jprev1, double mprev1,double anprev1, \
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
Line 2098  Interval (in months) between two waves: Line 2137  Interval (in months) between two waves:
  fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n   fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
  - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n   - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
    - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
    - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n
  - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n   - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
  - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n   - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
  - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);   - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
   
  if(popforecast==1) fprintf(fichtm,"\n   if(popforecast==1) fprintf(fichtm,"\n
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n   - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
Line 2116  fprintf(fichtm," <li>Graphs</li><p>"); Line 2157  fprintf(fichtm," <li>Graphs</li><p>");
  jj1=0;   jj1=0;
  for(k1=1; k1<=m;k1++){   for(k1=1; k1<=m;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){     for(i1=1; i1<=ncodemax[k1];i1++){
        jj1++;       jj1++;
        if (cptcovn > 0) {       if (cptcovn > 0) {
          fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
          for (cpt=1; cpt<=cptcoveff;cpt++)         for (cpt=1; cpt<=cptcoveff;cpt++)
            fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);           fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
          fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        }       }
        /* Pij */       /* Pij */
        fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>
 <img src=\"pe%s%d1.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);      <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
        /* Quasi-incidences */       /* Quasi-incidences */
        fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>
 <img src=\"pe%s%d2.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);      <img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
        /* Stable prevalence in each health state */         /* Stable prevalence in each health state */
        for(cpt=1; cpt<nlstate;cpt++){         for(cpt=1; cpt<nlstate;cpt++){
          fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>           fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
Line 2148  health expectancies in states (1) and (2 Line 2189  health expectancies in states (1) and (2
 <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);  <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
 fprintf(fichtm,"\n</body>");  fprintf(fichtm,"\n</body>");
    }     }
    }   }
 fclose(fichtm);  fclose(fichtm);
 }  }
   
Line 2158  void printinggnuplot(char fileres[],char Line 2199  void printinggnuplot(char fileres[],char
   int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;    int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
   int ng;    int ng;
   strcpy(optionfilegnuplot,optionfilefiname);    strcpy(optionfilegnuplot,optionfilefiname);
   strcat(optionfilegnuplot,".gp.txt");    strcat(optionfilegnuplot,".gp");
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {    if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     printf("Problem with file %s",optionfilegnuplot);      printf("Problem with file %s",optionfilegnuplot);
   }    }
Line 2303  fprintf(ficgp,"\" t \"e%d1\" w l",cpt); Line 2344  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
          for(k=1; k<=(nlstate+ndeath); k++) {           for(k=1; k<=(nlstate+ndeath); k++) {
            if (k != k2){             if (k != k2){
              if(ng==2)               if(ng==2)
                fprintf(ficgp," %f*exp(p%d+p%d*x",stepm/YEARM,i,i+1);                 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
              else               else
                fprintf(ficgp," exp(p%d+p%d*x",i,i+1);                 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
              ij=1;               ij=1;
Line 2674  int main(int argc, char *argv[]) Line 2715  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.8d, May 2002, INED-EUROREVES ";    char version[80]="Imach version 0.8g, May 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
Line 3285  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3326  while((c=getc(ficpar))=='#' && c!= EOF){
   
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
   
     
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
Removed from v.1.43  
changed lines
  Added in v.1.46

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