--- imach/src/imach.c 2017/05/18 20:09:32 1.268 +++ imach/src/imach.c 2017/05/23 08:39:25 1.269 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $ +/* $Id: imach.c,v 1.269 2017/05/23 08:39:25 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.269 2017/05/23 08:39:25 brouard + Summary: Code into subroutine, cleanings + Revision 1.268 2017/05/18 20:09:32 brouard Summary: backprojection and confidence intervals of backprevalence @@ -1007,12 +1010,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $ */ +/* $Id: imach.c,v 1.269 2017/05/23 08:39:25 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.268 $ $Date: 2017/05/18 20:09:32 $"; +char fullversion[]="$Revision: 1.269 $ $Date: 2017/05/23 08:39:25 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1084,10 +1087,7 @@ FILE *ficrescveij; char filerescve[FILENAMELENGTH]; FILE *ficresvij; char fileresv[FILENAMELENGTH]; -FILE *ficresvpl; -char fileresvpl[FILENAMELENGTH]; -FILE *ficresvbl; -char fileresvbl[FILENAMELENGTH]; + char title[MAXLINE]; char model[MAXLINE]; /**< The model line */ char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH]; @@ -2932,14 +2932,15 @@ double **pmij(double **ps, double *cov, /* Diag(w_x) */ /* Problem with prevacurrent which can be zero */ sumnew=0.; - /* for (ii=1;ii<=nlstate+ndeath;ii++){ */ + /*for (ii=1;ii<=nlstate+ndeath;ii++){*/ for (ii=1;ii<=nlstate;ii++){ /* Only on live states */ + /* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]); */ sumnew+=prevacurrent[(int)agefin][ii][ij]; } if(sumnew >0.01){ /* At least some value in the prevalence */ for (ii=1;ii<=nlstate+ndeath;ii++){ for (j=1;j<=nlstate+ndeath;j++) - doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0); + doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0); } }else{ for (ii=1;ii<=nlstate+ndeath;ii++){ @@ -5694,7 +5695,8 @@ void concatwav(int wav[], int **dh, int /* 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]);*/ } - + + /* Standard deviation of expectancies ij */ fprintf(ficresstdeij,"%3.0f",age ); for(i=1; i<=nlstate;i++){ eip=0.; @@ -5709,6 +5711,7 @@ void concatwav(int wav[], int **dh, int } fprintf(ficresstdeij,"\n"); + /* Variance of expectancies ij */ fprintf(ficrescveij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ @@ -6062,7 +6065,7 @@ void concatwav(int wav[], int **dh, int } /* end varevsij */ /************ Variance of prevlim ******************/ - void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres) + void varprevlim(char fileresvpl[], FILE *ficresvpl, double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres) { /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ @@ -6187,7 +6190,7 @@ void concatwav(int wav[], int **dh, int /************ Variance of backprevalence limit ******************/ - void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres) + void varbrevlim(char fileresvbl[], FILE *ficresvbl, double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres) { /* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ @@ -8112,7 +8115,7 @@ set ter svg size 640, 480\nunset log y\n /************** Forecasting ******************/ -void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ + void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ /* proj1, year, month, day of starting projection agemin, agemax range of age dateprev1 dateprev2 range of dates during which prevalence is computed @@ -8218,7 +8221,7 @@ void prevforecast(char fileres[], double ppij=0.; for(i=1; i<=nlstate;i++) { /* if (mobilav>=1) */ - ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k]; + ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k]; /* else { */ /* even if mobilav==-1 we use mobaverage */ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ /* } */ @@ -8239,12 +8242,13 @@ void prevforecast(char fileres[], double } -/* /\************** Back Forecasting ******************\/ */ -void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ +/************** Back Forecasting ******************/ + void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ /* back1, year, month, day of starting backection agemin, agemax range of age dateprev1 dateprev2 range of dates during which prevalence is computed - anback2 year of en of backection (same day and month as back1). + anback2 year of end of backprojection (same day and month as back1). + prevacurrent and prev are prevalences. */ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; double agec; /* generic age */ @@ -8304,10 +8308,6 @@ void prevbackforecast(char fileres[], do fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); - /* if (h==(int)(YEARM*yearp)){ */ - /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */ - /* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */ - /* k=k+1; */ for(nres=1; nres <= nresult; nres++) /* For each resultline */ for(k=1; k<=i1;k++){ if(i1 != 1 && TKresult[nres]!= k) @@ -8333,9 +8333,7 @@ void prevbackforecast(char fileres[], do /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */ fprintf(ficresfb,"\n"); fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); - printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); - /* for (agec=fage; agec>=(ageminpar-1); agec--){ */ - /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */ + /* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */ for (agec=bage; agec<=agemax-1; agec++){ /* testing */ /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/ nhstepm=(int) rint((agec-agelim)*YEARM/stepm); @@ -8360,8 +8358,10 @@ void prevbackforecast(char fileres[], do ppij=0.;ppi=0.; for(j=1; j<=nlstate;j++) { /* if (mobilav==1) */ - ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; - ppi=ppi+mobaverage[(int)agec][j][k]; + ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k]; + ppi=ppi+prevacurrent[(int)agec][j][k]; + /* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; */ + /* ppi=ppi+mobaverage[(int)agec][j][k]; */ /* else { */ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ /* } */ @@ -8386,6 +8386,123 @@ void prevbackforecast(char fileres[], do } +/* Variance of prevalence limit: varprlim */ + void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ + /*------- Variance of period (stable) prevalence------*/ + + char fileresvpl[FILENAMELENGTH]; + FILE *ficresvpl; + double **oldm, **savm; + double **varpl; /* Variances of prevalence limits by age */ + int i1, k, nres, j ; + + strcpy(fileresvpl,"VPL_"); + strcat(fileresvpl,fileresu); + if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { + printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); + exit(0); + } + printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); + fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); + + /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ + + i1=pow(2,cptcoveff); + if (cptcovn < 1){i1=1;} + + for(nres=1; nres <= nresult; nres++) /* For each resultline */ + for(k=1; k<=i1;k++){ + if(i1 != 1 && TKresult[nres]!= k) + continue; + fprintf(ficresvpl,"\n#****** "); + printf("\n#****** "); + fprintf(ficlog,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } + for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ + printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + } + fprintf(ficresvpl,"******\n"); + printf("******\n"); + fprintf(ficlog,"******\n"); + + varpl=matrix(1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + varprevlim(fileresvpl, ficresvpl, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyearp, k, strstart, nres); + free_matrix(varpl,1,nlstate,(int) bage, (int)fage); + /*}*/ + } + + fclose(ficresvpl); + printf("done variance-covariance of period prevalence\n");fflush(stdout); + fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); + + } +/* Variance of back prevalence: varbprlim */ + void varbprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **bprlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ + /*------- Variance of back (stable) prevalence------*/ + + char fileresvbl[FILENAMELENGTH]; + FILE *ficresvbl; + + double **oldm, **savm; + double **varbpl; /* Variances of back prevalence limits by age */ + int i1, k, nres, j ; + + strcpy(fileresvbl,"VBL_"); + strcat(fileresvbl,fileresu); + if((ficresvbl=fopen(fileresvbl,"w"))==NULL) { + printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl); + exit(0); + } + printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout); + fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog); + + + i1=pow(2,cptcoveff); + if (cptcovn < 1){i1=1;} + + for(nres=1; nres <= nresult; nres++) /* For each resultline */ + for(k=1; k<=i1;k++){ + if(i1 != 1 && TKresult[nres]!= k) + continue; + fprintf(ficresvbl,"\n#****** "); + printf("\n#****** "); + fprintf(ficlog,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } + for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ + printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + } + fprintf(ficresvbl,"******\n"); + printf("******\n"); + fprintf(ficlog,"******\n"); + + varbpl=matrix(1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + + varbrevlim(fileresvbl, ficresvbl, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, ncvyearp, k, strstart, nres); + free_matrix(varbpl,1,nlstate,(int) bage, (int)fage); + /*}*/ + } + + fclose(ficresvbl); + printf("done variance-covariance of back prevalence\n");fflush(stdout); + fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog); + + } /* End of varbprlim */ + /************** Forecasting *****not tested NB*************/ /* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */ @@ -10494,7 +10611,7 @@ int main(int argc, char *argv[]) double *delti; /* Scale */ double ***eij, ***vareij; double **varpl; /* Variances of prevalence limits by age */ - double **varbpl; /* Variances of back prevalence limits by age */ + double *epj, vepp; double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; @@ -11994,16 +12111,17 @@ Please run with mle=-1 to get a correct k=1; varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); - /* Prevalence for each covariates in probs[age][status][cov] */ - probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); - for(i=1;i<=AGESUP;i++) + /* Prevalence for each covariate combination in probs[age][status][cov] */ + probs= ma3x(AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); + for(i=AGEINF;i<=AGESUP;i++) for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */ for(k=1;k<=ncovcombmax;k++) probs[i][j][k]=0.; - prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, + ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); if (mobilav!=0 ||mobilavproj !=0 ) { - mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); - for(i=1;i<=AGESUP;i++) + mobaverages= ma3x(AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); + for(i=AGEINF;i<=AGESUP;i++) for(j=1;j<=nlstate+ndeath;j++) for(k=1;k<=ncovcombmax;k++) mobaverages[i][j][k]=0.; @@ -12015,31 +12133,27 @@ Please run with mle=-1 to get a correct fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); printf(" Error in movingaverage mobilav=%d\n",mobilav); } - } - /* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */ - /* for(i=1;i<=AGESUP;i++) */ - /* for(j=1;j<=nlstate;j++) */ - /* for(k=1;k<=ncovcombmax;k++) */ - /* mobaverages[i][j][k]=probs[i][j][k]; */ - /* /\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */ - /* /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */ - /* } */ - else if (mobilavproj !=0) { + } else if (mobilavproj !=0) { printf("Movingaveraging projected observed prevalence\n"); fprintf(ficlog,"Movingaveraging projected observed prevalence\n"); if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); } + }else{ + printf("Internal error moving average\n"); + fflush(stdout); + exit(1); } }/* end if moving average */ /*---------- Forecasting ------------------*/ - /*if((stepm == 1) && (strcmp(model,".")==0)){*/ if(prevfcast==1){ /* if(stepm ==1){*/ - prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); + prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); } + + /* Backcasting */ if(backcast==1){ ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); @@ -12048,71 +12162,23 @@ Please run with mle=-1 to get a correct /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ bprlim=matrix(1,nlstate,1,nlstate); + back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); fclose(ficresplb); hBijx(p, bage, fage, mobaverage); fclose(ficrespijb); - /* free_matrix(bprlim,1,nlstate,1,nlstate); /\*here or after loop ? *\/ */ - prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, - bage, fage, firstpass, lastpass, anback2, p, cptcoveff); + prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, + mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); + varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); - /*------- Variance of back (stable) prevalence------*/ - strcpy(fileresvbl,"VBL_"); - strcat(fileresvbl,fileresu); - if((ficresvbl=fopen(fileresvbl,"w"))==NULL) { - printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl); - exit(0); - } - printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout); - fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog); - - /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ - - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} - - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ - if(i1 != 1 && TKresult[nres]!= k) - continue; - fprintf(ficresvbl,"\n#****** "); - printf("\n#****** "); - fprintf(ficlog,"\n#****** "); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - } - for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - } - fprintf(ficresvbl,"******\n"); - printf("******\n"); - fprintf(ficlog,"******\n"); - - varbpl=matrix(1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - - varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres); - free_matrix(varbpl,1,nlstate,(int) bage, (int)fage); - /*}*/ - } - - fclose(ficresvbl); - printf("done variance-covariance of back prevalence\n");fflush(stdout); - fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog); - + free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath); free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath); free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath); - } - free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ + } /* end Backcasting */ /* ------ Other prevalence ratios------------ */ @@ -12165,10 +12231,10 @@ Please run with mle=-1 to get a correct fclose(ficreseij); printf("done evsij\n");fflush(stdout); fprintf(ficlog,"done evsij\n");fflush(ficlog); + /*---------- State-specific expectancies and variances ------------*/ - strcpy(filerest,"T_"); strcat(filerest,fileresu); if((ficrest=fopen(filerest,"w"))==NULL) { @@ -12177,8 +12243,6 @@ Please run with mle=-1 to get a correct } printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout); fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog); - - strcpy(fileresstde,"STDE_"); strcat(fileresstde,fileresu); if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { @@ -12206,9 +12270,6 @@ Please run with mle=-1 to get a correct printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout); fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog); - /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ - i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ if (cptcovn < 1){i1=1;} @@ -12270,7 +12331,7 @@ Please run with mle=-1 to get a correct vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); pstamp(ficrest); - + epj=vector(1,nlstate+1); for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ oldm=oldms;savm=savms; /* ZZ Segmentation fault */ cptcod= 0; /* To be deleted */ @@ -12286,7 +12347,6 @@ Please run with mle=-1 to get a correct for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ - epj=vector(1,nlstate+1); printf("Computing age specific period (stable) prevalences in each health state \n"); fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); for(age=bage; age <=fage ;age++){ @@ -12324,66 +12384,19 @@ Please run with mle=-1 to get a correct fprintf(ficrest,"\n"); } } /* End vpopbased */ + free_vector(epj,1,nlstate+1); free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); - free_vector(epj,1,nlstate+1); printf("done selection\n");fflush(stdout); fprintf(ficlog,"done selection\n");fflush(ficlog); - /*}*/ } /* End k selection */ printf("done State-specific expectancies\n");fflush(stdout); fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); - /*------- Variance of period (stable) prevalence------*/ - - strcpy(fileresvpl,"VPL_"); - strcat(fileresvpl,fileresu); - if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); - exit(0); - } - printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); - fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); - - /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ - - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} - - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ - if(i1 != 1 && TKresult[nres]!= k) - continue; - fprintf(ficresvpl,"\n#****** "); - printf("\n#****** "); - fprintf(ficlog,"\n#****** "); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - } - for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); - } - fprintf(ficresvpl,"******\n"); - printf("******\n"); - fprintf(ficlog,"******\n"); - - varpl=matrix(1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres); - free_matrix(varpl,1,nlstate,(int) bage, (int)fage); - /*}*/ - } - - fclose(ficresvpl); - printf("done variance-covariance of period prevalence\n");fflush(stdout); - fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); + /* variance-covariance of period prevalence*/ + varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); free_vector(weight,1,n); @@ -12402,8 +12415,8 @@ Please run with mle=-1 to get a correct /*---------- End : free ----------------*/ if (mobilav!=0 ||mobilavproj !=0) - free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */ - free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); + free_ma3x(mobaverages,AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */ + free_ma3x(probs,AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); } /* mle==-3 arrives here for freeing */ @@ -12420,6 +12433,7 @@ Please run with mle=-1 to get a correct /*free_vector(delti,1,npar);*/ free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); free_matrix(agev,1,maxwav,1,imx); + free_ma3x(paramstart,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); free_ivector(ncodemax,1,NCOVMAX);