--- imach/src/imach.c 2002/05/07 15:53:01 1.41 +++ imach/src/imach.c 2002/05/24 13:01:48 1.44 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.41 2002/05/07 15:53:01 lievre Exp $ +/* $Id: imach.c,v 1.44 2002/05/24 13:01:48 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -56,7 +56,7 @@ #include #define MAXLINE 256 -#define GNUPLOTPROGRAM "wgnuplot" +#define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ @@ -1330,10 +1330,10 @@ void prevalence(int agemin, float agemax j=cptcoveff; 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++){ j1++; - + for (i=-1; i<=nlstate+ndeath; i++) for (jk=-1; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) @@ -1352,41 +1352,43 @@ void prevalence(int agemin, float agemax if ((k2>=dateprev1) && (k2<=dateprev2)) { if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m0) 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]; + if (m0) + 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(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) + for(i=(int)agemin; i <= (int)agemax+3; i++){ + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) + pp[jk] += freq[jk][m][i]; + } + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pos=0; m <=0 ; m++) pos += freq[jk][m][i]; } - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - 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; jk <=nlstate ; jk++){ + for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) + 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; + } + } } + + } } } @@ -1663,6 +1665,7 @@ void evsij(char fileres[], double ***eij 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); @@ -1721,7 +1724,7 @@ void varevsij(char fileres[], double *** double age,agelim, hf; int theta; - fprintf(ficresvij,"# Covariances of life expectancies\n"); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); fprintf(ficresvij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1856,7 +1859,7 @@ void varprevlim(char fileres[], double * double age,agelim; int theta; - fprintf(ficresvpl,"# Standard deviation of prevalences limit\n"); + fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -2058,10 +2061,12 @@ fprintf(ficresprob,"#One-step probabilit /******************* Printing html file ***********/ void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ - int lastpass, int stepm, int weightopt, char model[],\ - int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ - char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ - char version[], int popforecast, int estepm ){ + int lastpass, int stepm, int weightopt, char model[],\ + int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ + char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ + char version[], int popforecast, int estepm ,/* \ */ + double jprev1, double mprev1,double anprev1, \ + double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ @@ -2072,26 +2077,30 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm," %s
\n + fprintf(fichtm," %s
\n Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n \n Total number of observations=%d
\n Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n -
-
  • Outputs files
    \n +
    +
    • Parameter files
      \n - Copy of the parameter file: o%s
      \n - - Gnuplot file name: %s
      \n - - Observed prevalence in each state: p%s
      \n - - Stationary prevalence in each state: pl%s
      \n - - Transition probabilities: pij%s
      \n - - Life expectancies by age and initial health status (estepm=%2d months): e%s
      \n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres); - - fprintf(fichtm,"\n - - Parameter file with estimated parameters and the covariance matrix: %s
      \n - - Variance of one-step probabilities: prob%s
      \n - - Variances of life expectancies by age and initial health status (estepm=%d months): v%s
      \n - - Health expectancies with their variances: t%s
      \n - - Standard deviation of stationary prevalences: vpl%s
      \n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); + - Gnuplot file name: %s
    \n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot); + + fprintf(fichtm,"
    • Result files (first order: no variance)
      \n + - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): p%s
      \n + - Estimated transition probabilities over %d (stepm) months: pij%s
      \n + - Stable prevalence in each health state: pl%s
      \n + - Life expectancies by age and initial health status (estepm=%2d months): + e%s
      \n
    • ", \ + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres); + + fprintf(fichtm,"\n
    • Result files (second order: variances)
      \n + - Parameter file with estimated parameters and covariance matrix: %s
      \n + - Variance of one-step probabilities: prob%s
      \n + - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): v%s
      \n + - Health expectancies with their variances (no covariance): t%s
      \n + - Standard deviation of stable prevalences: vpl%s
      \n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); if(popforecast==1) fprintf(fichtm,"\n - Prevalences forecasting: f%s
      \n @@ -2114,24 +2123,29 @@ fprintf(fichtm,"
    • Graphs
    • "); fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); fprintf(fichtm," ************\n


      "); } - fprintf(fichtm,"
      - Probabilities: pe%s%d.gif
      -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); + /* Pij */ + fprintf(fichtm,"
      - Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png
      +",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1); + /* Quasi-incidences */ + fprintf(fichtm,"
      - 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
      +",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1); + /* Stable prevalence in each health state */ for(cpt=1; cpt- Prevalence of disability : p%s%d%d.gif
      -",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"
      - Stable prevalence in each health state : p%s%d%d.png
      +",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
      - Observed and stationary prevalence (with confident -interval) in state (%d): v%s%d%d.gif
      -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); +interval) in state (%d): v%s%d%d.png
      +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
      - Health life expectancies by age and initial health state (%d): exp%s%d%d.gif
      -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"\n
      - Health life expectancies by age and initial health state (%d): exp%s%d%d.png
      +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } fprintf(fichtm,"\n
      - Total life expectancy by age and -health expectancies in states (1) and (2): e%s%d.gif
      -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); +health expectancies in states (1) and (2): e%s%d.png
      +",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); fprintf(fichtm,"\n"); } } @@ -2142,7 +2156,7 @@ fclose(fichtm); void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; - + int ng; strcpy(optionfilegnuplot,optionfilefiname); strcat(optionfilegnuplot,".gp.txt"); if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { @@ -2159,9 +2173,11 @@ m=pow(2,cptcoveff); for (k1=1; k1<= m ; k1 ++) { #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",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 #ifdef unix +fprintf(ficgp,"\nset out \"v%s%d%d.png\" \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); #endif @@ -2181,15 +2197,15 @@ 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)); #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 -fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); } } /*2 eme*/ 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",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 ++) { k=2*i; @@ -2214,7 +2230,6 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" if (i== (nlstate+1)) 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*/ @@ -2222,7 +2237,8 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { 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",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); @@ -2235,15 +2251,15 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt); 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 */ for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt