--- imach/src/imach.c 2006/03/20 10:52:43 1.123 +++ imach/src/imach.c 2006/04/28 18:11:50 1.127 @@ -1,6 +1,32 @@ -/* $Id: imach.c,v 1.123 2006/03/20 10:52:43 brouard Exp $ +/* $Id: imach.c,v 1.127 2006/04/28 18:11:50 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.127 2006/04/28 18:11:50 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + (Module): In order to speed up (in case of numerous covariates) we + compute health expectancies (without variances) in a first step + and then all the health expectancies with variances or standard + deviation (needs data from the Hessian matrices) which slows the + computation. + In the future we should be able to stop the program is only health + expectancies and graph are needed without standard deviations. + + Revision 1.126 2006/04/28 17:23:28 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + Version 0.98h + + Revision 1.125 2006/04/04 15:20:31 lievre + Errors in calculation of health expectancies. Age was not initialized. + Forecasting file added. + + Revision 1.124 2006/03/22 17:13:53 lievre + Parameters are printed with %lf instead of %f (more numbers after the comma). + The log-likelihood is printed in the log file + Revision 1.123 2006/03/20 10:52:43 brouard * imach.c (Module): changed, corresponds to .htm file name. <head> headers where missing. @@ -357,11 +383,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.123 2006/03/20 10:52:43 brouard Exp $ */ +/* $Id: imach.c,v 1.127 2006/04/28 18:11:50 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.123 $ $Date: 2006/03/20 10:52:43 $"; +char version[]="Imach version 0.98h, April 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.127 $ $Date: 2006/04/28 18:11:50 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1001,9 +1027,8 @@ void powell(double p[], double **xi, int last_time=curr_time; (void) gettimeofday(&curr_time,&tzp); printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); - /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); - fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); - */ + fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog); +/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */ for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); @@ -2454,11 +2479,12 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) +void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) { /* Health expectancies, no variances */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; + int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; double ***p3mat; double eip; @@ -2505,18 +2531,28 @@ void evsij(char fileres[], double ***eij hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ agelim=AGESUP; - /* nhstepm age range expressed in number of stepm */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepm matrices, stored + in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + +/* nhstepm age range expressed in number of stepm */ + nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ /* Typically if 20 years nstepm = 20*12/6=40 stepm */ /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - /* Computed by stepm unit matrices, product of hstepm matrices, stored - in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + for (age=bage; age<=fage; age ++){ + nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */ + + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepma matrices, stored + in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij); hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ @@ -2551,13 +2587,14 @@ void evsij(char fileres[], double ***eij } -void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) +void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) { /* Covariances of health expectancies eij and of total life expectancies according to initial status i, ei. . */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; + int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; double ***p3matp, ***p3matm, ***varhe; double **dnewm,**doldm; @@ -2631,7 +2668,7 @@ void cvevsij(char fileres[], double ***e /* If stepm=6 months */ /* nhstepm age range expressed in number of stepm */ agelim=AGESUP; - nstepm=(int) rint((agelim-age)*YEARM/stepm); + nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Typically if 20 years nstepm = 20*12/6=40 stepm */ /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ @@ -2644,10 +2681,15 @@ void cvevsij(char fileres[], double ***e gm=matrix(0,nhstepm,1,nlstate*nlstate); for (age=bage; age<=fage; age ++){ - - /* Computed by stepm unit matrices, product of hstepm matrices, stored - in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - + nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */ + + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepma matrices, stored + in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ /* Computing Variances of health expectancies */ @@ -2698,6 +2740,7 @@ void cvevsij(char fileres[], double ***e varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; } } + /* Computing expectancies */ hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); for(i=1; i<=nlstate;i++) @@ -3481,9 +3524,11 @@ void printinghtml(char fileres[], char t subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); fprintf(fichtm,"\ - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ - <a href=\"%s\">%s</a> <br>\n</li>", + <a href=\"%s\">%s</a> <br>\n", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); - + fprintf(fichtm,"\ + - Population projections by age and states: \ + <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f")); fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); @@ -5357,9 +5402,9 @@ Interval (in months) between two waves: fprintf(ficlog,"%d%d ",i,k); fprintf(ficres,"%1d%1d ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%f ",p[jk]); - fprintf(ficlog,"%f ",p[jk]); - fprintf(ficres,"%f ",p[jk]); + printf("%lf ",p[jk]); + fprintf(ficlog,"%lf ",p[jk]); + fprintf(ficres,"%lf ",p[jk]); jk++; } printf("\n"); @@ -5706,18 +5751,24 @@ Interval (in months) between two waves: } - /*---------- Health expectancies and variances ------------*/ + /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ - strcpy(filerest,"t"); - strcat(filerest,fileres); - if((ficrest=fopen(filerest,"w"))==NULL) { - printf("Problem with total LE resultfile: %s\n", filerest);goto end; - fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; + prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ + ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); + */ + + if (mobilav!=0) { + mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + 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); + } } - printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); - fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); + /*---------- Health expectancies, no variances ------------*/ + strcpy(filerese,"e"); strcat(filerese,fileres); if((ficreseij=fopen(filerese,"w"))==NULL) { @@ -5726,6 +5777,37 @@ Interval (in months) between two waves: } printf("Computing Health Expectancies: result on file '%s' \n", filerese); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + fprintf(ficreseij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + fprintf(ficreseij,"******\n"); + + eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart); + + free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); + } + } + fclose(ficreseij); + + + /*---------- Health expectancies and variances ------------*/ + + + strcpy(filerest,"t"); + strcat(filerest,fileres); + if((ficrest=fopen(filerest,"w"))==NULL) { + printf("Problem with total LE resultfile: %s\n", filerest);goto end; + fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; + } + printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); + fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); + strcpy(fileresstde,"stde"); strcat(fileresstde,fileres); @@ -5754,20 +5836,6 @@ Interval (in months) between two waves: printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ - prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); - /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ - ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); - */ - - if (mobilav!=0) { - mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - 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); - } - } - for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; @@ -5776,15 +5844,12 @@ Interval (in months) between two waves: fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); fprintf(ficrest,"******\n"); - fprintf(ficreseij,"\n#****** "); fprintf(ficresstdeij,"\n#****** "); fprintf(ficrescveij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); } - fprintf(ficreseij,"******\n"); fprintf(ficresstdeij,"******\n"); fprintf(ficrescveij,"******\n"); @@ -5795,8 +5860,7 @@ Interval (in months) between two waves: eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart); - cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); + cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; @@ -5853,7 +5917,6 @@ Interval (in months) between two waves: free_matrix(mint,1,maxwav,1,n); free_ivector(cod,1,n); free_ivector(tab,1,NCOVMAX); - fclose(ficreseij); fclose(ficresstdeij); fclose(ficrescveij); fclose(ficresvij);