From d609b475ce9bccc5332da1a170ae92490311ee5f Mon Sep 17 00:00:00 2001 From: =?utf8?q?Agn=C3=A8s=20Li=C3=A8vre?= Date: Tue, 26 Feb 2002 17:11:54 +0000 Subject: [PATCH] Cleaning of the main function: gnuplot and html file are now subroutines; the bug in the forecasting is corrected; the moving average statement moved to the last command line. --- src/imach.c | 584 +++++++++++++++++++++++++++------------------------- 1 file changed, 304 insertions(+), 280 deletions(-) diff --git a/src/imach.c b/src/imach.c index 7ce5471..2049be2 100644 --- a/src/imach.c +++ b/src/imach.c @@ -95,7 +95,7 @@ double jmean; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf; -FILE *ficgp, *fichtm,*ficresprob,*ficpop; +FILE *ficgp,*ficresprob,*ficpop; FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -123,7 +123,7 @@ FILE *ficreseij; static double maxarg1,maxarg2; #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2)) #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2)) - + #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) #define rint(a) floor(a+0.5) @@ -1879,6 +1879,270 @@ fclose(ficresprob); exit(0); } +/******************* 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[] ){ + int jj1, k1, i1, cpt; + FILE *fichtm; + /*char optionfilehtm[FILENAMELENGTH];*/ + + strcpy(optionfilehtm,optionfile); + strcat(optionfilehtm,".htm"); + if((fichtm=fopen(optionfilehtm,"w"))==NULL) { + printf("Problem with %s \n",optionfilehtm), exit(0); + } + + fprintf(fichtm,"

Imach, Version 0.7

Outputs files

\n + - Observed prevalence in each state: p%s
\n +- Estimated parameters and the covariance matrix: %s
+ - Stationary prevalence in each state: pl%s
+ - Transition probabilities: pij%s
+ - Copy of the parameter file: o%s
+ - Life expectancies by age and initial health status: e%s
+ - Variances of life expectancies by age and initial health status: v%s
+ - Health expectancies with their variances: t%s
+ - Standard deviation of stationary prevalences: vpl%s
+ - Prevalences and population forecasting: f%s
+
",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,"
Graphs

"); + + m=cptcoveff; + if (cptcovn < 1) {m=1;ncodemax[1]=1;} + + jj1=0; + for(k1=1; k1<=m;k1++){ + for(i1=1; i1<=ncodemax[k1];i1++){ + jj1++; + if (cptcovn > 0) { + fprintf(fichtm,"

- Prevalence of disability : p%s%d%d.gif
+

",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); + } + 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
- Total life expectancy by age and +health expectancies in states (1) and (2): e%s%d.gif
+

",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); +fprintf(fichtm,"\n"); + } + } +fclose(fichtm); +} + +/******************* Gnuplot file **************/ +void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){ + + int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; + + strcpy(optionfilegnuplot,optionfilefiname); + strcat(optionfilegnuplot,".plt"); + if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { + printf("Problem with file %s",optionfilegnuplot); + } + +#ifdef windows + fprintf(ficgp,"cd \"%s\" \n",pathc); +#endif +m=pow(2,cptcoveff); + + /* 1eme*/ + for (cpt=1; cpt<= nlstate ; cpt ++) { + 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",agemin,fage,fileres,k1-1,k1-1); +#endif +#ifdef unix +fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres); +#endif + +for (i=1; i<= nlstate ; i ++) { + if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1); + for (i=1; i<= nlstate ; i ++) { + if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1); + for (i=1; i<= nlstate ; i ++) { + if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + 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"); +#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] ",agemin,fage); + + for (i=1; i<= nlstate+1 ; i ++) { + k=2*i; + fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1); + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,"); + else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1); + fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1); + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + fprintf(ficgp,"\" t\"\" w l 0,"); + fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1); + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); +} + 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*/ + + for (k1=1; k1<= m ; k1 ++) { + for (cpt=1; cpt<= nlstate ; cpt ++) { + k=2+nlstate*(cpt-1); + 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",agemin,fage,fileres,k1-1,k1-1,k,cpt); + 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,"\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

Imach, Version 0.7

Outputs files

\n - - Observed prevalence in each state: p%s
\n -- Estimated parameters and the covariance matrix: %s
- - Stationary prevalence in each state: pl%s
- - Transition probabilities: pij%s
- - Copy of the parameter file: o%s
- - Life expectancies by age and initial health status: e%s
- - Variances of life expectancies by age and initial health status: v%s
- - Health expectancies with their variances: t%s
- - Standard deviation of stationary prevalences: vpl%s
- - Prevalences and population forecasting: f%s
-
",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,"
Graphs

"); - - m=cptcoveff; - if (cptcovn < 1) {m=1;ncodemax[1]=1;} - - j1=0; - for(k1=1; k1<=m;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - if (cptcovn > 0) { - fprintf(fichtm,"

- Prevalence of disability : p%s%d%d.gif
-

",strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1); - } - 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,j1,strtok(optionfile, "."),cpt,j1); - } - 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,j1,strtok(optionfile, "."),cpt,j1); - } - fprintf(fichtm,"\n
- Total life expectancy by age and -health expectancies in states (1) and (2): e%s%d.gif
-

",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1); -fprintf(fichtm,"\n"); - } - } -fclose(fichtm); + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm); + /*--------------- Prevalence limit --------------*/ strcpy(filerespl,"pl"); @@ -2883,27 +2918,14 @@ fclose(fichtm); free_matrix(mint,1,maxwav,1,n); free_matrix(anint,1,maxwav,1,n); free_matrix(agev,1,maxwav,1,imx); + /* Mobile average */ if (cptcoveff==0) ncodemax[cptcoveff]=1; if (mobilav==1) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - for (agedeb=bage+3; agedeb<=fage-2; agedeb++) - for (i=1; i<=nlstate;i++) - 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 (cptcod=1;cptcod<=ncodemax[cptcoveff];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; - } - } - } + movingaverage(agedeb, fage, agemin, mobaverage); } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -2953,16 +2975,15 @@ fclose(fichtm); fprintf(ficresf,"# StartingAge FinalAge"); for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j); if (popforecast==1) fprintf(ficresf," [Population]"); - - for (cpt=0; cpt<4;cpt++) { + + for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { fprintf(ficresf,"\n"); fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */ + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; - /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/ - + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); oldm=oldms;savm=savms; hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); @@ -2978,7 +2999,7 @@ fclose(fichtm); kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; else { kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; - /* fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h], probs[(int)(agedeb)+1][i][cptcod]);*/ + } if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb]; @@ -2991,26 +3012,29 @@ fclose(fichtm); } } } - /* 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); + if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (popforecast==1) { free_ivector(popage,0,AGESUP); free_vector(popeffectif,0,AGESUP); free_vector(popcount,0,AGESUP); } free_imatrix(s,1,maxwav+1,1,n); - free_vector(weight,1,n);*/ + free_vector(weight,1,n); fclose(ficresf); }/* End forecasting */ else{ 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); + } + /*---------- Health expectancies and variances ------------*/ strcpy(filerest,"t"); -- 2.43.0