From 555dc6e0bf3faa3625a561e835ff9c98ae838d1a Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Tue, 17 Nov 2015 14:31:57 +0000 Subject: [PATCH] Summary: temporary --- src/imach.c | 305 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 191 insertions(+), 114 deletions(-) diff --git a/src/imach.c b/src/imach.c index 4577bfe..99ac78b 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,9 @@ /* $Id$ $State$ $Log$ + Revision 1.207 2015/10/27 17:36:57 brouard + *** empty log message *** + Revision 1.206 2015/10/24 07:14:11 brouard *** empty log message *** @@ -1986,7 +1989,7 @@ double **prevalim(double **prlim, int nlstate, double x[], double age, double ** /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **pmij(); double **newm; - double agefin, delaymax=100 ; /* Max number of years to converge */ + double agefin, delaymax=100. ; /* 100 Max number of years to converge */ int ncvloop=0; for (ii=1;ii<=nlstate+ndeath;ii++) @@ -2036,21 +2039,47 @@ double **prevalim(double **prlim, int nlstate, double x[], double age, double ** prlim[i][j]= newm[i][j]/(1-sumnew); max=FMAX(max,prlim[i][j]); min=FMIN(min,prlim[i][j]); - printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); + /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */ } maxmin=(max-min)/(max+min)*2; maxmax=FMAX(maxmax,maxmin); + /* for(i=1; i<=nlstate; i++) { */ + /* sumnew=0.; */ + /* sumnew+=prlim[i][j]; */ + /* } */ + /* prmimj = sumnew/(float)nlstate; /\* Means of various prevalence limits. */ } /* j loop */ *ncvyear= (int)age- (int)agefin; - printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); + /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ if(maxmax < ftolpl){ /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ return prlim; } } /* age loop */ - printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\ -Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); -/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ + /* After some age loop it doesn't converge */ + printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g within %.0f years. \n\ +Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + /* printf(" age= %d newm\n",(int)age); */ + /* for(i=1; i<=nlstate+ndeath; i++) { */ + /* for(j=1;j<=nlstate+ndeath;j++){ */ + /* printf(" %lf", newm[i][j]); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* printf("\n"); */ + /* printf("prlim\n"); */ + /* for(i=1; i<=nlstate; i++) { */ + /* sumnew=0; */ + /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */ + /* for(j=1;j<=nlstate;j++){ */ + /* prlim[i][j]= newm[i][j]/(1-sumnew); */ + /* printf(" %lf", prlim[i][j]); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* printf("\n"); */ + + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ return prlim; /* should not reach here */ } @@ -2683,16 +2712,16 @@ void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, lon fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ -",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); - fprintf(fichtm,"
- and by state of destination %s-dest.png
\ -",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); - fflush(fichtm); for (k=1; k<= nlstate ; k++) { fprintf(fichtm,"
- Probability p%dj by origin %d and destination j %s-p%dj.png
\ ",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); } + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fflush(fichtm); } return; } @@ -2957,6 +2986,8 @@ double hessij( double x[], double **hess, double delti[], int thetai,int thetaj, double p2[MAXPARM+1]; int k, kmax=1; double v1, v2, cv12, lc1, lc2; + + int firstime=0; fx=func(x); for (k=1; k<=kmax; k=k+10) { @@ -2978,13 +3009,14 @@ double hessij( double x[], double **hess, double delti[], int thetai,int thetaj, k4=func(p2)-fx; res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ if(k1*k2*k3*k4 <0.){ + firstime=1; kmax=kmax+10; - if(kmax >=10){ + } + if(kmax >=10 || firstime ==1){ printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); - } } #ifdef DEBUGHESSIJ v1=hess[thetai][thetai]; @@ -4044,7 +4076,6 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); } printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); - fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); pstamp(ficresprobmorprev); fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); @@ -4055,6 +4086,7 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j); } fprintf(ficresprobmorprev,"\n"); + fprintf(ficgp,"\n# Routine varevsij"); fprintf(ficgp,"\nunset title \n"); /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ @@ -4092,9 +4124,9 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa /* For example we decided to compute the life expectancy with the smallest unit */ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. nhstepm is the number of hstepm from age to agelim - nstepm is the number of stepm from age to agelin. - Look at function hpijx to understand why (it is linked to memory size questions) */ - /* We decided (b) to get a life expectancy respecting the most precise curvature of the + nstepm is the number of stepm from age to agelim. + Look at function hpijx to understand why (it is linked to memory size questions) + we decided (b) to get a life expectancy respecting the most precise curvature of the survival function given by stepm (the optimization length). Unfortunately it means that if the survival funtion is printed every two years of age and if you sum them up and add 1 year (area under the trapezoids) you won't get the same @@ -4310,6 +4342,7 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mgm, **mgp; double age,agelim; int theta; @@ -4331,7 +4364,10 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ if (stepm >= YEARM) hstepm=1; nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); gradg=matrix(1,npar,1,nlstate); + mgp=matrix(1,npar,1,nlstate); + mgm=matrix(1,npar,1,nlstate); gp=vector(1,nlstate); gm=vector(1,nlstate); @@ -4339,16 +4375,19 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Missing or not useful because 1 year */ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); - for(i=1;i<=nlstate;i++) + for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; - + mgp[theta][i] = prlim[i][i]; + } for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); - for(i=1;i<=nlstate;i++) + for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; - + mgm[theta][i] = prlim[i][i]; + } for(i=1;i<=nlstate;i++) gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; } /* End theta */ @@ -4358,10 +4397,28 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[j][theta]=gradg[theta][j]; + if((int)age==68 ||(int)age== 69 ){ + printf("\nmgm mgp %d ",(int)age); + for(j=1; j<=nlstate;j++){ + printf("%d ",j); + for(theta=1; theta <=npar; theta++) + printf("%d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); + printf("\n "); + } + } + if((int)age==68 ||(int)age== 69 ){ + printf("\n gradg %d ",(int)age); + for(j=1; j<=nlstate;j++){ + printf("%d ",j); + for(theta=1; theta <=npar; theta++) + printf("%d %lf ",theta,gradg[theta][j]); + printf("\n "); + } + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] =0.; - if((int)age==67 ||(int)age== 66 ){ + if((int)age==68 ||(int)age== 69 ){ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); }else{ @@ -4377,6 +4434,8 @@ void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fa fprintf(ficresvpl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); + free_matrix(mgm,1,npar,1,nlstate); + free_matrix(mgp,1,npar,1,nlstate); free_matrix(gradg,1,npar,1,nlstate); free_matrix(trgradg,1,nlstate,1,npar); } /* End age */ @@ -4788,7 +4847,7 @@ divided by h: hPij/h : %s_%d-3.svg
\ } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Survival functions from state %d in any different live states and total.\ + fprintf(fichtm,"
\n- Survival functions from state %d in each live state and total.\ Or probability to survive in various states (1 to %d) being in state %d at different ages.\ %s%d_%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); } @@ -8039,8 +8098,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); } - printf("Computing Health Expectancies: result on file '%s' \n", filerese); - fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout); + fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog); /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ @@ -8059,7 +8118,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa /*}*/ } fclose(ficreseij); - + printf("done evsij\n");fflush(stdout); + fprintf(ficlog,"done evsij\n");fflush(ficlog); /*---------- Health expectancies and variances ------------*/ @@ -8070,8 +8130,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa 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); + 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_"); @@ -8080,8 +8140,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); } - printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); - fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); strcpy(filerescve,"CVE_"); strcat(filerescve,fileresu); @@ -8089,8 +8149,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); } - printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); - fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); strcpy(fileresv,"V_"); strcat(fileresv,fileresu); @@ -8098,97 +8158,109 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa printf("Problem with variance resultfile: %s\n", fileresv);exit(0); fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); } - printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout); + fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog); /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ for (k=1; k <= (int) pow(2,cptcoveff); k++){ - fprintf(ficrest,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficrest,"******\n"); - - fprintf(ficresstdeij,"\n#****** "); - fprintf(ficrescveij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - } - fprintf(ficresstdeij,"******\n"); - fprintf(ficrescveij,"******\n"); - - fprintf(ficresvij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficresvij,"******\n"); - - eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); - /* - */ - /* goto endfree; */ - - vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); - pstamp(ficrest); - - - 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 */ - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ - fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); - if(vpopbased==1) - fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); - else - fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); - fprintf(ficrest,"# Age popbased mobilav e.. (std) "); - 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); - for(age=bage; age <=fage ;age++){ - prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */ - if (vpopbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][k]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][k]; - } - } - - fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); - /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ - /* printf(" age %4.0f ",age); */ - for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ - for(i=1, epj[j]=0.;i <=nlstate;i++) { - epj[j] += prlim[i][i]*eij[i][j][(int)age]; - /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ - /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ - } - epj[nlstate+1] +=epj[j]; + fprintf(ficrest,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrest,"******\n"); + + fprintf(ficresstdeij,"\n#****** "); + fprintf(ficrescveij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + } + fprintf(ficresstdeij,"******\n"); + fprintf(ficrescveij,"******\n"); + + fprintf(ficresvij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresvij,"******\n"); + + eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + printf(" cvevsij %d, ",k); + fprintf(ficlog, " cvevsij %d, ",k); + cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); + printf(" end cvevsij \n "); + fprintf(ficlog, " end cvevsij \n "); + + /* + */ + /* goto endfree; */ + + vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); + pstamp(ficrest); + + + 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 */ + printf("varevsij %d \n",vpopbased); + fprintf(ficlog, "varevsij %d \n",vpopbased); + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); + if(vpopbased==1) + fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); + else + fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); + fprintf(ficrest,"# Age popbased mobilav e.. (std) "); + 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++){ + prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */ + if (vpopbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][k]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][k]; } - /* printf(" age %4.0f \n",age); */ - - for(i=1, vepp=0.;i <=nlstate;i++) - for(j=1;j <=nlstate;j++) - vepp += vareij[i][j][(int)age]; - fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); - for(j=1;j <=nlstate;j++){ - fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); + } + + fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); + /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ + /* printf(" age %4.0f ",age); */ + for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ + for(i=1, epj[j]=0.;i <=nlstate;i++) { + epj[j] += prlim[i][i]*eij[i][j][(int)age]; + /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ + /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ } - fprintf(ficrest,"\n"); + epj[nlstate+1] +=epj[j]; + } + /* printf(" age %4.0f \n",age); */ + + for(i=1, vepp=0.;i <=nlstate;i++) + for(j=1;j <=nlstate;j++) + vepp += vareij[i][j][(int)age]; + fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); + for(j=1;j <=nlstate;j++){ + fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); } + fprintf(ficrest,"\n"); } - 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); + } /* End vpopbased */ + 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 \n");fflush(stdout); + fprintf(ficlog,"done\n");fflush(ficlog); + /*}*/ - } + } /* End k */ free_vector(weight,1,n); free_imatrix(Tvard,1,NCOVMAX,1,2); free_imatrix(s,1,maxwav+1,1,n); @@ -8200,6 +8272,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa fclose(ficrescveij); fclose(ficresvij); fclose(ficrest); + printf("done Health expectancies\n");fflush(stdout); + fprintf(ficlog,"done Health expectancies\n");fflush(ficlog); fclose(ficpar); /*------- Variance of period (stable) prevalence------*/ @@ -8210,7 +8284,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl); + 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++){*/ @@ -8229,6 +8304,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa } fclose(ficresvpl); + printf("done variance-covariance of period prevalence\n");fflush(stdout); + fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); /*---------- End : free ----------------*/ if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); -- 2.43.0