--- imach/src/imach.c 2002/07/23 23:59:37 1.53 +++ imach/src/imach.c 2002/07/24 09:07:45 1.54 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.53 2002/07/23 23:59:37 brouard Exp $ +/* $Id: imach.c,v 1.54 2002/07/24 09:07:45 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -39,7 +39,7 @@ hPijx. Also this programme outputs the covariance matrix of the parameters but also - of the life expectancies. It also computes the prevalence limits. + of the life expectancies. It also computes the stable prevalence. Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). Institut national d'études démographiques, Paris. @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.8j, July 2002, INED-EUROREVES "; +char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES "; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -701,12 +701,12 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } + } } } } -/**** Prevalence limit ****************/ +/**** Prevalence limit (stable prevalence) ****************/ double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) { @@ -1827,13 +1827,16 @@ void varevsij(char optionfilefiname[], d strcpy(digitp,"-populbased-"); else strcpy(digitp,"-stablbased-"); - if(mobilav==1) + if(mobilav!=0) strcat(digitp,"mobilav-"); else strcat(digitp,"nomobil-"); - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, bage, fage, mobaverage); + if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } strcpy(fileresprobmorprev,"prmorprev"); @@ -1928,10 +1931,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -1956,10 +1959,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -2025,10 +2028,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -2106,7 +2109,7 @@ void varprevlim(char fileres[], double * double age,agelim; int theta; - fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n"); + fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -2599,9 +2602,9 @@ void printinggnuplot(char fileres[], dou fprintf(ficlog,"Problem with file %s",optionfilegnuplot); } -#ifdef windows + /*#ifdef windows */ fprintf(ficgp,"cd \"%s\" \n",pathc); -#endif + /*#endif */ m=pow(2,cptcoveff); /* 1eme*/ @@ -2614,7 +2617,7 @@ m=pow(2,cptcoveff); 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); + fprintf(ficgp,"\" t\"Stable 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)"); @@ -2768,27 +2771,39 @@ m=pow(2,cptcoveff); /*************** Moving average **************/ -void movingaverage(double ***probs, double bage,double fage, double ***mobaverage){ +int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){ int i, cpt, cptcod; + int mobilavrange, mob; double age; - for (age=bage; age<=fage; age++) - for (i=1; i<=nlstate;i++) - for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) - mobaverage[(int)age][i][cptcod]=0.; - - for (age=bage+4; age<=fage; age++){ - for (i=1; i<=nlstate;i++){ - for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ - for (cpt=0;cpt<=4;cpt++){ - mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]+probs[(int)age-cpt][i][cptcod]; + if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ + if(mobilav==1) mobilavrange=5; /* default */ + else mobilavrange=mobilav; + for (age=bage; age<=fage; age++) + for (i=1; i<=nlstate;i++) + for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) + mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; + /* We keep the original values on the extreme ages bage, fage and for + fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 + we use a 5 terms etc. until the borders are no more concerned. + */ + for (mob=3;mob <=mobilavrange;mob=mob+2){ + for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ + for (i=1; i<=nlstate;i++){ + for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ + mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; + for (cpt=1;cpt<=(mob-1)/2;cpt++){ + mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; + mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; + } + mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; + } } - mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]/5; - } - } - } - -} + }/* end age */ + }/* end mob */ + }else return -1; + return 0; +}/* End movingaverage */ /************** Forecasting ******************/ @@ -2818,9 +2833,12 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 if (cptcoveff==0) ncodemax[cptcoveff]=1; - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar,fage, mobaverage); + if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -2891,7 +2909,7 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 } } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficresf); } @@ -2924,9 +2942,12 @@ populforecast(char fileres[], double anp if (cptcoveff==0) ncodemax[cptcoveff]=1; - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar, fage, mobaverage); + if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -3039,7 +3060,7 @@ populforecast(char fileres[], double anp } } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); if (popforecast==1) { free_ivector(popage,0,AGESUP); @@ -3750,12 +3771,17 @@ while((c=getc(ficpar))=='#' && c!= EOF){ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ - strcpy(optionfilegnuplot,optionfilefiname); - strcat(optionfilegnuplot,".gp"); - if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { - printf("Problem with file %s",optionfilegnuplot); - } - fclose(ficgp); + strcpy(optionfilegnuplot,optionfilefiname); + strcat(optionfilegnuplot,".gp"); + if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { + printf("Problem with file %s",optionfilegnuplot); + } + else{ + fprintf(ficgp,"\n# %s\n", version); + fprintf(ficgp,"# %s\n", optionfilegnuplot); + fprintf(ficgp,"set missing 'NaNq'\n"); +} + fclose(ficgp); printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p); /*--------- index.htm --------*/ @@ -3792,17 +3818,17 @@ Interval (in months) between two waves: fclose(ficres); - /*--------------- Prevalence limit --------------*/ + /*--------------- Prevalence limit (stable prevalence) --------------*/ strcpy(filerespl,"pl"); strcat(filerespl,fileres); if((ficrespl=fopen(filerespl,"w"))==NULL) { - printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end; - fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end; + printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; + fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; } - printf("Computing prevalence limit: result on file '%s' \n", filerespl); - fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl); - fprintf(ficrespl,"#Prevalence limit\n"); + printf("Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficrespl,"#Stable prevalence \n"); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -3950,9 +3976,12 @@ Interval (in months) between two waves: fprintf(ficlog,"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); - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar, fage, mobaverage); + if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } k=0; @@ -3994,10 +4023,10 @@ Interval (in months) between two waves: for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][k]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][k]; } @@ -4032,15 +4061,15 @@ free_matrix(mint,1,maxwav,1,n); fclose(ficpar); free_vector(epj,1,nlstate+1); - /*------- Variance limit prevalence------*/ + /*------- Variance of stable prevalence------*/ strcpy(fileresvpl,"vpl"); strcat(fileresvpl,fileres); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance prev lim resultfile: %s\n", fileresvpl); + printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl); + printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); k=0; for(cptcov=1;cptcov<=i1;cptcov++){ @@ -4075,7 +4104,7 @@ free_matrix(mint,1,maxwav,1,n); free_vector(delti,1,npar); free_matrix(agev,1,maxwav,1,imx); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fprintf(fichtm,"\n"); fclose(fichtm); @@ -4109,9 +4138,10 @@ free_matrix(mint,1,maxwav,1,n); strcpy(plotcmd,GNUPLOTPROGRAM); strcat(plotcmd," "); strcat(plotcmd,optionfilegnuplot); + printf("Starting: %s\n",plotcmd);fflush(stdout); system(plotcmd); -#ifdef windows + /*#ifdef windows*/ while (z[0] != 'q') { /* chdir(path); */ printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); @@ -4121,7 +4151,7 @@ free_matrix(mint,1,maxwav,1,n); else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } -#endif + /*#endif */ }