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

Diff for /imach/src/imach.c between versions 1.39 and 1.42

version 1.39, 2002/04/05 15:45:00 version 1.42, 2002/05/21 18:44:41
Line 56 Line 56
 #include <unistd.h>  #include <unistd.h>
   
 #define MAXLINE 256  #define MAXLINE 256
 #define GNUPLOTPROGRAM "wgnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
Line 1330  void prevalence(int agemin, float agemax Line 1330  void prevalence(int agemin, float agemax
   j=cptcoveff;    j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
     
  for(k1=1; k1<=j;k1++){    for(k1=1; k1<=j;k1++){
     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=-1; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
Line 1352  void prevalence(int agemin, float agemax Line 1352  void prevalence(int agemin, float agemax
             if ((k2>=dateprev1) && (k2<=dateprev2)) {              if ((k2>=dateprev1) && (k2<=dateprev2)) {
               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;
               if (m<lastpass) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];                if (m<lastpass) {
               /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  */                  if (calagedate>0)
                     freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
                   else
                     freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                   freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];
                 }
             }              }
           }            }
         }          }
       }        }
         for(i=(int)agemin; i <= (int)agemax+3; i++){        for(i=(int)agemin; i <= (int)agemax+3; 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++)
             pos += freq[jk][m][i];              pos += freq[jk][m][i];
         }          }
                 
          for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){
            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)            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++) pos += pp[jk];  
   
          for(jk=1; jk <=nlstate ; jk++){            
            if( i <= (int) agemax){  
              if(pos>=1.e-5){  
                probs[i][jk][j1]= pp[jk]/pos;  
              }  
            }  
          }  
            
         }          }
          
           for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
          
           for(jk=1; jk <=nlstate ; jk++){    
             if( i <= (int) agemax){
               if(pos>=1.e-5){
                 probs[i][jk][j1]= pp[jk]/pos;
               }
             }
           }
          
         }
     }      }
   }    }
    
     
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
Line 1531  void tricode(int *Tvar, int **nbcode, in Line 1536  void tricode(int *Tvar, int **nbcode, in
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm)  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )
   
 {  {
   /* Health expectancies */    /* Health expectancies */
   int i, j, nhstepm, hstepm, h, nstepm;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
   double age, agelim, hf;    double age, agelim, hf;
   double ***p3mat;    double ***p3mat,***varhe;
     double **dnewm,**doldm;
     double *xp;
     double **gp, **gm;
     double ***gradg, ***trgradg;
     int theta;
   
     varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);
     xp=vector(1,npar);
     dnewm=matrix(1,nlstate*2,1,npar);
     doldm=matrix(1,nlstate*2,1,nlstate*2);
     
   fprintf(ficreseij,"# Health expectancies\n");    fprintf(ficreseij,"# Health expectancies\n");
   fprintf(ficreseij,"# Age");    fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)      for(j=1; j<=nlstate;j++)
       fprintf(ficreseij," %1d-%1d",i,j);        fprintf(ficreseij," %1d-%1d (SE)",i,j);
   fprintf(ficreseij,"\n");    fprintf(ficreseij,"\n");
   
   if(estepm < stepm){    if(estepm < stepm){
Line 1582  void evsij(char fileres[], double ***eij Line 1598  void evsij(char fileres[], double ***eij
     /* 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);
       gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);
       gp=matrix(0,nhstepm,1,nlstate*2);
       gm=matrix(0,nhstepm,1,nlstate*2);
   
     /* 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);  
     
     /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/  
   
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
       /* Computing Variances of health expectancies */
   
        for(theta=1; theta <=npar; theta++){
         for(i=1; i<=npar; i++){
           xp[i] = x[i] + (i==theta ?delti[theta]:0);
         }
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
    
         cptj=0;
         for(j=1; j<= nlstate; j++){
           for(i=1; i<=nlstate; i++){
             cptj=cptj+1;
             for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){
               gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
             }
           }
         }
        
        
         for(i=1; i<=npar; i++)
           xp[i] = x[i] - (i==theta ?delti[theta]:0);
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
        
         cptj=0;
         for(j=1; j<= nlstate; j++){
           for(i=1;i<=nlstate;i++){
             cptj=cptj+1;
             for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
               gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
             }
           }
         }
        
      
   
         for(j=1; j<= nlstate*2; j++)
           for(h=0; h<=nhstepm-1; h++){
             gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
           }
   
        }
      
   /* End theta */
   
        trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);
   
        for(h=0; h<=nhstepm-1; h++)
         for(j=1; j<=nlstate*2;j++)
           for(theta=1; theta <=npar; theta++)
           trgradg[h][j][theta]=gradg[h][theta][j];
   
   
        for(i=1;i<=nlstate*2;i++)
         for(j=1;j<=nlstate*2;j++)
           varhe[i][j][(int)age] =0.;
   
        printf("%d||",(int)age);fflush(stdout);
       for(h=0;h<=nhstepm-1;h++){
         for(k=0;k<=nhstepm-1;k++){
           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]);
           for(i=1;i<=nlstate*2;i++)
             for(j=1;j<=nlstate*2;j++)
               varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
         }
       }
   
        
       /* Computing expectancies */
     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++){
           eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;            eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
           /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/           
   /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
   
         }          }
   
     fprintf(ficreseij,"%3.0f",age );      fprintf(ficreseij,"%3.0f",age );
       cptj=0;
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);          cptj++;
           fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficreseij,"\n");
      
       free_matrix(gm,0,nhstepm,1,nlstate*2);
       free_matrix(gp,0,nhstepm,1,nlstate*2);
       free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);
       free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }    }
     free_vector(xp,1,npar);
     free_matrix(dnewm,1,nlstate*2,1,npar);
     free_matrix(doldm,1,nlstate*2,1,nlstate*2);
     free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);
 }  }
   
 /************ Variance ******************/  /************ Variance ******************/
Line 1842  void varprob(char fileres[], double **ma Line 1945  void varprob(char fileres[], double **ma
   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);
   }    }
   printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);    printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
     
   fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
     fprintf(ficresprob,"# Age");
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=(nlstate+ndeath);j++)
         fprintf(ficresprob," p%1d-%1d (SE)",i,j);
   
   
     fprintf(ficresprob,"\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);
Line 1933  void varprob(char fileres[], double **ma Line 2045  void varprob(char fileres[], double **ma
         fprintf(ficresprob,"\n%d ",(int)age);          fprintf(ficresprob,"\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],doldm[i][i]);            fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
     
       }        }
     }      }
Line 1979  Interval (in months) between two waves: Line 2091  Interval (in months) between two waves:
   
  fprintf(fichtm,"\n   fprintf(fichtm,"\n
  - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n   - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n
  - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n    - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
    - Variances 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: <a href=\"t%s\">t%s</a> <br>\n   - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n
  - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);   - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,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 2004  fprintf(fichtm," <li>Graphs</li><p>"); Line 2117  fprintf(fichtm," <li>Graphs</li><p>");
            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\">");
        }         }
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>         fprintf(fichtm,"<br>- Probabilities: pe%s%d.png<br>
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);      <img src=\"pe%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
        for(cpt=1; cpt<nlstate;cpt++){         for(cpt=1; cpt<nlstate;cpt++){
          fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>           fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.png<br>
 <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }         }
     for(cpt=1; cpt<=nlstate;cpt++) {      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident         fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
 interval) in state (%d): v%s%d%d.gif <br>  interval) in state (%d): v%s%d%d.png <br>
 <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);    <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>          fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
 <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and       fprintf(fichtm,"\n<br>- Total life expectancy by age and
 health expectancies in states (1) and (2): e%s%d.gif<br>  health expectancies in states (1) and (2): e%s%d.png<br>
 <img src=\"e%s%d.gif\">",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>");
    }     }
    }     }
Line 2049  m=pow(2,cptcoveff); Line 2162  m=pow(2,cptcoveff);
    for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) {
   
 #ifdef windows  #ifdef windows
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);       fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
 #endif  #endif
 #ifdef unix  #ifdef unix
   fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
 fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);  fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
 #endif  #endif
   
Line 2071  for (i=1; i<= nlstate ; i ++) { Line 2186  for (i=1; i<= nlstate ; i ++) {
 }    }  
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));       fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
 #ifdef unix  #ifdef unix
 fprintf(ficgp,"\nset ter gif small size 400,300");  fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
 #endif  #endif
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
    }     }
   }    }
   /*2 eme*/    /*2 eme*/
   
   for (k1=1; k1<= m ; k1 ++) {    for (k1=1; k1<= m ; k1 ++) {
     fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);      fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1);
       fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
         
     for (i=1; i<= nlstate+1 ; i ++) {      for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;        k=2*i;
Line 2104  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2219  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");        if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
       else fprintf(ficgp,"\" t\"\" w l 0,");        else fprintf(ficgp,"\" t\"\" w l 0,");
     }      }
     fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);  
   }    }
     
   /*3eme*/    /*3eme*/
   
   for (k1=1; k1<= m ; k1 ++) {    for (k1=1; k1<= m ; k1 ++) {
     for (cpt=1; cpt<= nlstate ; cpt ++) {      for (cpt=1; cpt<= nlstate ; cpt ++) {
       k=2+nlstate*(cpt-1);        k=2+nlstate*(2*cpt-2);
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);        fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
         /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
    for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
   fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
   fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
    for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
   fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
   
   */
       for (i=1; i< nlstate ; i ++) {        for (i=1; i< nlstate ; i ++) {
         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);          fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);
   
       }        }
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
     }      }
     }
     
   /* CV preval stat */    /* CV preval stat */
     for (k1=1; k1<= m ; k1 ++) {      for (k1=1; k1<= m ; k1 ++) {
     for (cpt=1; cpt<nlstate ; cpt ++) {      for (cpt=1; cpt<nlstate ; cpt ++) {
       k=3;        k=3;
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);        fprintf(ficgp,"set out \"p%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
   
       for (i=1; i< nlstate ; i ++)        for (i=1; i< nlstate ; i ++)
         fprintf(ficgp,"+$%d",k+i+1);          fprintf(ficgp,"+$%d",k+i+1);
Line 2137  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2261  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
         fprintf(ficgp,"+$%d",l+i+1);          fprintf(ficgp,"+$%d",l+i+1);
       }        }
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);          fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);  
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }      }
   }      }  
     
Line 2153  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2276  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
         }          }
       }        }
     }      }
     }     }
   
     for(jk=1; jk <=m; jk++) {     for(jk=1; jk <=m; jk++) {
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);       fprintf(ficgp,"\nset out \"pe%s%d.png\" \n\n",strtok(optionfile, "."),jk);
    i=1;       fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
    for(k2=1; k2<=nlstate; k2++) {       i=1;
      k3=i;       for(k2=1; k2<=nlstate; k2++) {
      for(k=1; k<=(nlstate+ndeath); k++) {         k3=i;
        if (k != k2){         for(k=1; k<=(nlstate+ndeath); k++) {
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);           if (k != k2){
 ij=1;             fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
         for(j=3; j <=ncovmodel; j++) {             ij=1;
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {             for(j=3; j <=ncovmodel; j++) {
             fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);               if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
             ij++;                 fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
           }                 ij++;
           else               }
           fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);               else
         }                 fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
           fprintf(ficgp,")/(1");             }
                     fprintf(ficgp,")/(1");
         for(k1=1; k1 <=nlstate; k1++){               
           fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);             for(k1=1; k1 <=nlstate; k1++){  
 ij=1;               fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
           for(j=3; j <=ncovmodel; j++){               ij=1;
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {               for(j=3; j <=ncovmodel; j++){
             fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);                 if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
             ij++;                   fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
           }                   ij++;
           else                 }
             fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);                 else
           }                   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
           fprintf(ficgp,")");               }
         }               fprintf(ficgp,")");
         fprintf(ficgp,") t \"p%d%d\" ", k2,k);             }
         if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");             fprintf(ficgp,") t \"p%d%d\" ", k2,k);
         i=i+ncovmodel;             if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
              i=i+ncovmodel;
            }
        }         }
      }       }
    }     }
    fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);  
    }  
         
   fclose(ficgp);     fclose(ficgp);
 }  /* end gnuplot */  }  /* end gnuplot */
   
   
Line 2514  int main(int argc, char *argv[]) Line 2637  int main(int argc, char *argv[])
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
   int mobilav=0,popforecast=0;    int mobilav=0,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
   double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;    double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 2532  int main(int argc, char *argv[]) Line 2655  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.8b, March 2002, INED-EUROREVES ";    char version[80]="Imach version 0.8d, 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 3227  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3350  while((c=getc(ficpar))=='#' && c!= EOF){
             for(j=1; j<=nlstate+ndeath;j++)              for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespij," %1d-%1d",i,j);                fprintf(ficrespij," %1d-%1d",i,j);
           fprintf(ficrespij,"\n");            fprintf(ficrespij,"\n");
           for (h=0; h<=nhstepm; h++){             for (h=0; h<=nhstepm; h++){
             fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );              fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
             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(ficrespij," %.5f", p3mat[i][j][h]);                  fprintf(ficrespij," %.5f", p3mat[i][j][h]);
             fprintf(ficrespij,"\n");              fprintf(ficrespij,"\n");
           }               }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespij,"\n");            fprintf(ficrespij,"\n");
         }          }
Line 3249  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3372  while((c=getc(ficpar))=='#' && c!= EOF){
   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);
     free_matrix(mint,1,maxwav,1,n);    }
     free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);  
     free_vector(weight,1,n);}  
   else{    else{
     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);
Line 3281  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3402  while((c=getc(ficpar))=='#' && c!= EOF){
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);      printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
   }    }
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
     calagedate=-1;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
   k=0;    k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
Line 3303  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3426  while((c=getc(ficpar))=='#' && c!= EOF){
   
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
       evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm);          evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  
    
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
        varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);         varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
Line 3326  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3450  while((c=getc(ficpar))=='#' && c!= EOF){
         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]*eij[i][j][(int)age];              epj[j] += prlim[i][i]*eij[i][j][(int)age];
               /*  printf("%lf %lf ", 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];
Line 3340  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3466  while((c=getc(ficpar))=='#' && c!= EOF){
       }        }
     }      }
   }    }
   free_matrix(mint,1,maxwav,1,n);
       free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
       free_vector(weight,1,n);
   fclose(ficreseij);    fclose(ficreseij);
   fclose(ficresvij);    fclose(ficresvij);
   fclose(ficrest);    fclose(ficrest);
Removed from v.1.39  
changed lines
  Added in v.1.42

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