/* $Id$
$State$
$Log$
+ Revision 1.266 2017/05/13 07:26:12 brouard
+ Summary: Version 0.99r13 (improvements and bugs fixed)
+
Revision 1.265 2017/04/26 16:22:11 brouard
Summary: imach 0.99r13 Some bugs fixed
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
- - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k, nres);
Computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
}
for(i=1; i<=nlstate+ndeath; i++)
for(j=1;j<=nlstate+ndeath;j++) {
- po[i][j][h]=newm[i][j];
- /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+ po[i][j][h]=newm[i][j];
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
}
/*printf("h=%d ",h);*/
} /* end h */
- /* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d \n",h); */
return po;
}
/************* Higher Back Matrix Product ***************/
/* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
-double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
+double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij, int nres )
{
/* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
/* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
-
}
- for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
- /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
- for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
+ for (k=1; k<=nsq;k++) { /* For single varying covariates only */
+ /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
+ cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
+ /* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
+ }
+ for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */
+ if(Dummy[Tvar[Tage[k]]]){
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ } else{
+ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
+ }
+ /* printf("hBxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
+ }
+ for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
-
+ }
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+
/* Careful transposed matrix */
/* age is in cov[2], prevacurrent at beginning of transition. */
/* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
/************** 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 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
anproj2 year of en of projection (same day and month as proj1).
*/
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
}
/* /\************** Back Forecasting ******************\/ */
-/* void prevbackforecast(char fileres[], 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). */
-/* *\/ */
-/* int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */
-/* double agec; /\* generic age *\/ */
-/* double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */
-/* double *popeffectif,*popcount; */
-/* double ***p3mat; */
-/* /\* double ***mobaverage; *\/ */
-/* char fileresfb[FILENAMELENGTH]; */
-
-/* agelim=AGESUP; */
-/* /\* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people */
-/* in each health status at the date of interview (if between dateprev1 and dateprev2). */
-/* We still use firstpass and lastpass as another selection. */
-/* *\/ */
-/* /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */
-/* /\* firstpass, lastpass, stepm, weightopt, model); *\/ */
-/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-
-/* strcpy(fileresfb,"FB_"); */
-/* strcat(fileresfb,fileresu); */
-/* if((ficresfb=fopen(fileresfb,"w"))==NULL) { */
-/* printf("Problem with back forecast resultfile: %s\n", fileresfb); */
-/* fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */
-/* } */
-/* printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-/* fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-
-/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
-
-/* /\* if (mobilav!=0) { *\/ */
-/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-/* /\* 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; */
-/* if (stepm<=12) stepsize=1; */
-/* if(estepm < stepm){ */
-/* printf ("Problem %d lower than %d\n",estepm, stepm); */
-/* } */
-/* else hstepm=estepm; */
-
-/* hstepm=hstepm/stepm; */
-/* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
-/* fractional in yp1 *\/ */
-/* anprojmean=yp; */
-/* yp2=modf((yp1*12),&yp); */
-/* mprojmean=yp; */
-/* yp1=modf((yp2*30.5),&yp); */
-/* jprojmean=yp; */
-/* if(jprojmean==0) jprojmean=1; */
-/* if(mprojmean==0) jprojmean=1; */
-
-/* i1=cptcoveff; */
-/* if (cptcovn < 1){i1=1;} */
+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).
+ */
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ double agec; /* generic age */
+ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double *popeffectif,*popcount;
+ double ***p3mat;
+ /* double ***mobaverage; */
+ char fileresfb[FILENAMELENGTH];
+
+ agelim=AGESUP;
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
+ /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
+ /* firstpass, lastpass, stepm, weightopt, model); */
+
+ /*Do we need to compute prevalence again?*/
+
+ /* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-/* fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); */
+ strcpy(fileresfb,"FB_");
+ strcat(fileresfb,fileresu);
+ if((ficresfb=fopen(fileresfb,"w"))==NULL) {
+ printf("Problem with back forecast resultfile: %s\n", fileresfb);
+ fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
+ }
+ printf("\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
+ fprintf(ficlog,"\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
-/* 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; */
-/* fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */
-/* for(j=1;j<=cptcoveff;j++) { */
-/* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* } */
-/* fprintf(ficresfb," yearbproj age"); */
-/* for(j=1; j<=nlstate+ndeath;j++){ */
-/* for(i=1; i<=nlstate;i++) */
-/* fprintf(ficresfb," p%d%d",i,j); */
-/* fprintf(ficresfb," p.%d",j); */
-/* } */
-/* for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { */
-/* /\* 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); */
-/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
-/* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
-/* nhstepm = nhstepm/hstepm; */
-/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* oldm=oldms;savm=savms; */
-/* hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k); */
-/* for (h=0; h<=nhstepm; h++){ */
-/* if (h*hstepm/YEARM*stepm ==yearp) { */
-/* fprintf(ficresfb,"\n"); */
-/* for(j=1;j<=cptcoveff;j++) */
-/* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */
-/* } */
-/* for(j=1; j<=nlstate+ndeath;j++) { */
-/* ppij=0.; */
-/* for(i=1; i<=nlstate;i++) { */
-/* if (mobilav==1) */
-/* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; */
-/* else { */
-/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */
-/* } */
-/* if (h*hstepm/YEARM*stepm== yearp) { */
-/* fprintf(ficresfb," %.3f", p3mat[i][j][h]); */
-/* } */
-/* } /\* end i *\/ */
-/* if (h*hstepm/YEARM*stepm==yearp) { */
-/* fprintf(ficresfb," %.3f", ppij); */
-/* } */
-/* }/\* end j *\/ */
-/* } /\* end h *\/ */
-/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* } /\* end agec *\/ */
-/* } /\* end yearp *\/ */
-/* } /\* end cptcod *\/ */
-/* } /\* end cptcov *\/ */
-
-/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-
-/* fclose(ficresfb); */
-/* printf("End of Computing Back forecasting \n"); */
-/* fprintf(ficlog,"End of Computing Back forecasting\n"); */
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+
+ hstepm=hstepm/stepm;
+ yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
+ fractional in yp1 */
+ anprojmean=yp;
+ yp2=modf((yp1*12),&yp);
+ mprojmean=yp;
+ yp1=modf((yp2*30.5),&yp);
+ jprojmean=yp;
+ if(jprojmean==0) jprojmean=1;
+ if(mprojmean==0) jprojmean=1;
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+
+ 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)
+ continue;
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) projection ignored because no cases \n",k);
+ continue;
+ }
+ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+ fprintf(ficresfb," yearbproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresfb," p%d%d",i,j);
+ fprintf(ficresfb," p.%d",j);
+ }
+ for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
+ /* 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);
+ /* for (agec=fage; agec>=(ageminpar-1); agec--){ */
+ /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
+ for (agec=fage; agec>=fage-20; agec--){ /* testing up to 10 */
+ nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresfb,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
+ }
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ /* if (mobilav==1) */
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+ /* else { */
+ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
+ /* } */
+ if (h*hstepm/YEARM*stepm== yearp) {
+ fprintf(ficresfb," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresfb," %.3f", ppij);
+ }
+ }/* end j */
+ } /* end h */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end k */
+
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+
+ fclose(ficresfb);
+ printf("End of Computing Back forecasting \n");
+ fprintf(ficlog,"End of Computing Back forecasting\n");
-/* } */
+}
/************** 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){ */
/* oldm=oldms;savm=savms; */
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
- hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+ hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
for(i=1; i<=nlstate;i++)
fclose(ficrespijb);
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
- /* prevbackforecast(fileresu, 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);
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);
git://henry.ined.fr / .git / commitdiff