--- imach/src/imach.c 2002/12/11 16:58:19 1.65 +++ imach/src/imach.c 2003/02/05 12:40:38 1.70 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.65 2002/12/11 16:58:19 lievre Exp $ +/* $Id: imach.c,v 1.70 2003/02/05 12:40:38 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -32,8 +32,8 @@ hPijx is the probability to be observed in state i at age x+h conditional to the observed state i at age x. The delay 'h' can be split into an exact number (nh*stepm) of unobserved intermediate - states. This elementary transition (by month or quarter trimester, - semester or year) is model as a multinomial logistic. The hPx + states. This elementary transition (by month, quarter, + semester or year) is modelled as a multinomial logistic. The hPx matrix is simply the matrix product of nh*stepm elementary matrices and the contribution of each individual to the likelihood is simply hPijx. @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.9, November 2002, INED-EUROREVES "; +char version[80]="Imach version 0.92, February 2003, INED-EUROREVES "; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -856,11 +856,13 @@ double **matprod2(double **out, double * double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) { - /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month - duration (i.e. until - age (in years) age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices. + /* 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 + nhstepm*hstepm matrices. Output is stored in matrix po[i][j][h] for h every 'hstepm' step - (typically every 2 years instead of every month which is too big). + (typically every 2 years instead of every month which is too big + for the memory). Model is determined by parameters x and covariates have to be included manually here. @@ -917,7 +919,7 @@ double func( double *x) double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ int s1, s2; - double bbh; + double bbh, survp; long ipmx; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -1554,13 +1556,18 @@ void freqsummary(char fileres[], int ag } /************ Prevalence ********************/ -void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate) -{ /* Some frequencies */ +void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) +{ + /* 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. + */ int i, m, jk, k1, i1, j1, bool, z1,z2,j; double ***freq; /* Frequencies */ double *pp; - double pos, k2; + double pos; + double y2; /* in fractional years */ pp=vector(1,nlstate); @@ -1579,7 +1586,7 @@ void prevalence(int agemin, float agemax for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; - for (i=1; i<=imx; i++) { + for (i=1; i<=imx; i++) { /* Each individual */ bool=1; if (cptcovn>0) { for (z1=1; z1<=cptcoveff; z1++) @@ -1587,20 +1594,17 @@ void prevalence(int agemin, float agemax bool=0; } if (bool==1) { - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - if ((k2>=dateprev1) && (k2<=dateprev2)) { + for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ + y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ + if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; if (m0) - freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i]; - else - freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; + freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; } } - } + } /* end selection of waves */ } } for(i=(int)agemin; i <= (int)agemax+3; i++){ @@ -1702,10 +1706,12 @@ void concatwav(int wav[], int **dh, int if (j <= jmin) jmin=j; sum=sum+j; /*if (j<0) printf("j=%d num=%d \n",j,i); */ + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ } } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; @@ -1740,7 +1746,6 @@ void concatwav(int wav[], int **dh, int bh[mi][i]=ju; /* At least one step */ printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i); } - if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm); } } /* end if mle */ } /* end wave */ @@ -1844,7 +1849,7 @@ void evsij(char fileres[], double ***eij * This is mainly to measure the difference between two models: for example * if stepm=24 months pijx are given only every 2 years and by summing them * we are calculating an estimate of the Life Expectancy assuming a linear - * progression inbetween and thus overestimating or underestimating according + * progression in between and thus overestimating or underestimating according * to the curvature of the survival function. If, for the same date, we * estimate the model with stepm=1 month, we can keep estepm to 24 months * to compare the new estimate of Life expectancy with the same linear @@ -2034,7 +2039,7 @@ void varevsij(char optionfilefiname[], d } 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); - fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1*p1j+w2*p2j+... stand dev in()\n"); + 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); fprintf(ficresprobmorprev,"# Age cov=%-d",ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ fprintf(ficresprobmorprev," p.%-d SE",j); @@ -2056,8 +2061,8 @@ void varevsij(char optionfilefiname[], d exit(0); } else{ - fprintf(fichtm,"\n

  • Computing probabilities of dying as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); - fprintf(fichtm,"\n
    %s (à revoir)
    \n",digitp); + fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); + fprintf(fichtm,"\n
    %s
    \n",digitp); } varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); @@ -2091,7 +2096,7 @@ void varevsij(char optionfilefiname[], d and note for a fixed period like k years */ /* 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 only each two years of age and if + 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 results. So we changed our mind and took the option of the best precision. */ @@ -2107,7 +2112,7 @@ void varevsij(char optionfilefiname[], d for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++){ /* Computes gradient */ + for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ xp[i] = x[i] + (i==theta ?delti[theta]:0); } hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); @@ -2129,14 +2134,17 @@ void varevsij(char optionfilefiname[], d gp[h][j] += prlim[i][i]*p3mat[i][j][h]; } } - /* This for computing forces of mortality (h=1)as a weighted average */ - for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) + /* This for computing probability of death (h=1 means + computed over hstepm matrices product = hstepm*stepm months) + as a weighted average of prlim. + */ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gpp[j]=0.; i<= nlstate; i++) gpp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end force of mortality */ + /* end probability of death */ - for(i=1; i<=npar; i++) /* Computes gradient */ + for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); @@ -2157,17 +2165,21 @@ void varevsij(char optionfilefiname[], d gm[h][j] += prlim[i][i]*p3mat[i][j][h]; } } - /* This for computing force of mortality (h=1)as a weighted average */ - for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) - gmp[j] += prlim[i][i]*p3mat[i][j][1]; + /* This for computing probability of death (h=1 means + computed over hstepm matrices product = hstepm*stepm months) + as a weighted average of prlim. + */ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.; i<= nlstate; i++) + gmp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end force of mortality */ + /* end probability of death */ for(j=1; j<= nlstate; j++) /* vareij */ for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } @@ -2184,6 +2196,7 @@ void varevsij(char optionfilefiname[], d for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ for(theta=1; theta <=npar; theta++) trgradgp[j][theta]=gradgp[theta][j]; + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) @@ -2199,7 +2212,7 @@ void varevsij(char optionfilefiname[], d vareij[i][j][(int)age] += doldm[i][j]*hf*hf; } } - + /* pptj */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); @@ -2207,6 +2220,7 @@ void varevsij(char optionfilefiname[], d for(i=nlstate+1;i<=nlstate+ndeath;i++) varppt[j][i]=doldmp[j][i]; /* end ppptj */ + /* x centered again */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); @@ -2219,13 +2233,16 @@ void varevsij(char optionfilefiname[], d prlim[i][i]=mobaverage[(int)age][i][ij]; } } - - /* This for computing force of mortality (h=1)as a weighted average */ - for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) + + /* This for computing probability of death (h=1 means + computed over hstepm (estepm) matrices product = hstepm*stepm months) + as a weighted average of prlim. + */ + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.;i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end force of mortality */ + /* end probability of death */ fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ @@ -2255,11 +2272,14 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); - fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); +/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ + fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",fileresprobmorprev); fprintf(fichtm,"\n
    File (multiple files are possible if covariates are present): %s\n",fileresprobmorprev,fileresprobmorprev); - fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", stepm,digitp,digit); + fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", estepm,digitp,digit); /* fprintf(fichtm,"\n
    Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

    \n", stepm,YEARM,digitp,digit); */ fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); @@ -2274,7 +2294,7 @@ void varevsij(char optionfilefiname[], d fclose(ficresprobmorprev); fclose(ficgp); fclose(fichtm); -} +} /************ 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 ij) @@ -2420,10 +2440,11 @@ void varprob(char optionfilefiname[], do fprintf(ficresprobcov," p%1d-%1d ",i,j); fprintf(ficresprobcor," p%1d-%1d ",i,j); } - fprintf(ficresprob,"\n"); + /* fprintf(ficresprob,"\n"); fprintf(ficresprobcov,"\n"); fprintf(ficresprobcor,"\n"); - xp=vector(1,npar); + */ + xp=vector(1,npar); dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); @@ -2462,14 +2483,14 @@ void varprob(char optionfilefiname[], do if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresprob, "**********\n#"); + fprintf(ficresprob, "**********\n#\n"); fprintf(ficresprobcov, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresprobcov, "**********\n#"); + fprintf(ficresprobcov, "**********\n#\n"); fprintf(ficgp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficgp, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#\n"); fprintf(fichtm, "\n
    ********** Variable "); @@ -2478,7 +2499,7 @@ void varprob(char optionfilefiname[], do fprintf(ficresprobcor, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficgp, "**********\n#"); + fprintf(ficresprobcor, "**********\n#"); } for (age=bage; age<=fage; age ++){ @@ -2797,12 +2818,12 @@ m=pow(2,cptcoveff); if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"Stable 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+1.96*$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)"); } - fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1); + fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-1.96*$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)"); @@ -2993,21 +3014,23 @@ int movingaverage(double ***probs, doubl /************** Forecasting ******************/ -prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){ - - int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; +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, c, cptcod, i, h; int *popage; - double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; + double agec; /* generic age */ + double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; double *popeffectif,*popcount; double ***p3mat; double ***mobaverage; char fileresf[FILENAMELENGTH]; - agelim=AGESUP; - calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; - - prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); - + agelim=AGESUP; + prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(fileresf,"f"); strcat(fileresf,fileres); @@ -3031,11 +3054,10 @@ prevforecast(char fileres[], double anpr stepsize=(int) (stepm+YEARM-1)/YEARM; if (stepm<=12) stepsize=1; - agelim=AGESUP; - hstepm=1; hstepm=hstepm/stepm; - yp1=modf(dateintmean,&yp); + yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and + fractional in yp1 */ anprojmean=yp; yp2=modf((yp1*12),&yp); mprojmean=yp; @@ -3044,49 +3066,54 @@ prevforecast(char fileres[], double anpr if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; - fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); + fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); - for(cptcov=1;cptcov<=i2;cptcov++){ + fprintf(ficresf,"#****** Routine prevforecast **\n"); + for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; fprintf(ficresf,"\n#******"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); } fprintf(ficresf,"******\n"); - fprintf(ficresf,"# StartingAge FinalAge"); - for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j); - - - for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { + fprintf(ficresf,"# Covariate valuofcovar yearproj age"); + for(j=1; j<=nlstate+ndeath;j++){ + for(i=1; i<=nlstate;i++) + fprintf(ficresf," p%d%d",i,j); + fprintf(ficresf," p.%d",j); + } + for (yearp=0; yearp<=(anproj2-anproj1);yearp++) { fprintf(ficresf,"\n"); - fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); + fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + for (agec=fage; agec>=ageminpar; 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; - hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedate+YEARM*cpt)) { - fprintf(ficresf,"\n %.f %.f ",anproj1+cpt,agedeb+h*hstepm/YEARM*stepm); + if (h==(int) (YEARM*yearp)) { + fprintf(ficresf,"\n"); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; + ppij=0.; for(i=1; i<=nlstate;i++) { if (mobilav==1) - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; + ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec+1][i][cptcod]; else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; + ppij=ppij+p3mat[i][j][h]*probs[(int)(agec+1)][i][cptcod]; } - + if (h==(int)(YEARM*yearp)) + fprintf(ficresf," %.3f", p3mat[i][j][h]); } - if (h==(int)(calagedate+12*cpt)){ - fprintf(ficresf," %.3f", kk1); - + if (h==(int)(YEARM*yearp)){ + fprintf(ficresf," %.3f", ppij); } } } @@ -3100,12 +3127,13 @@ prevforecast(char fileres[], double anpr fclose(ficresf); } -/************** Forecasting ******************/ + +/************** Forecasting *****not tested NB*************/ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; int *popage; - double calagedate, agelim, kk1, kk2; + double calagedatem, agelim, kk1, kk2; double *popeffectif,*popcount; double ***p3mat,***tabpop,***tabpopprev; double ***mobaverage; @@ -3114,9 +3142,9 @@ populforecast(char fileres[], double anp tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); agelim=AGESUP; - calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; + calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; - prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(filerespop,"pop"); @@ -3162,7 +3190,7 @@ populforecast(char fileres[], double anp for (i=1; i=(ageminpar-((int)calagedate %12)/12.); agedeb--){ + for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; @@ -3186,7 +3214,7 @@ populforecast(char fileres[], double anp hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedate+YEARM*cpt)) { + if (h==(int) (calagedatem+YEARM*cpt)) { fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { @@ -3198,7 +3226,7 @@ populforecast(char fileres[], double anp kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; } } - if (h==(int)(calagedate+12*cpt)){ + if (h==(int)(calagedatem+12*cpt)){ tabpop[(int)(agedeb)][j][cptcod]=kk1; /*fprintf(ficrespop," %.3f", kk1); if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/ @@ -3209,10 +3237,10 @@ populforecast(char fileres[], double anp for(j=1; j<=nlstate;j++){ kk1= kk1+tabpop[(int)(agedeb)][j][cptcod]; } - tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)]; + tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)]; } - if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) + if (h==(int)(calagedatem+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]); } free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); @@ -3223,7 +3251,7 @@ populforecast(char fileres[], double anp for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ + for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; @@ -3231,7 +3259,7 @@ populforecast(char fileres[], double anp oldm=oldms;savm=savms; hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedate+YEARM*cpt)) { + if (h==(int) (calagedatem+YEARM*cpt)) { fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { @@ -3239,7 +3267,7 @@ populforecast(char fileres[], double anp for(i=1; i<=nlstate;i++) { kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod]; } - if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %15.2f", kk1); + if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1); } } free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); @@ -3286,9 +3314,10 @@ int main(int argc, char *argv[]) int ju,jl, mi; int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab; + int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; int hstepm, nhstepm; - double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate; + double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1; double bage, fage, age, agelim, agebase; double ftolpl=FTOL; @@ -3703,7 +3732,9 @@ int main(int argc, char *argv[]) } } } - else if(s[m][i] !=9){ /* Should no more exist */ + else if(s[m][i] !=9){ /* Standard case, age in fractional + years but with the precision of a + month */ agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); if(mint[m][i]==99 || anint[m][i]==9999) agev[m][i]=1; @@ -3922,6 +3953,8 @@ int main(int argc, char *argv[]) fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); @@ -3932,8 +3965,8 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); - dateprev1=anprev1+mprev1/12.+jprev1/365.; - dateprev2=anprev2+mprev2/12.+jprev2/365.; + dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.; + dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; fscanf(ficpar,"pop_based=%d\n",&popbased); fprintf(ficparo,"pop_based=%d\n",popbased); @@ -3947,10 +3980,12 @@ int main(int argc, char *argv[]) } ungetc(c,ficpar); - fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2); - fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); - fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); - + fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj); + fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + /* day and month of proj2 are not used but only year anproj2.*/ while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); @@ -4057,7 +4092,9 @@ Interval (in months) between two waves: for (age=agebase; age<=agelim; age++){ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); - fprintf(ficrespl,"%.0f",age ); + fprintf(ficrespl,"%.0f ",age ); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); for(i=1; i<=nlstate;i++) fprintf(ficrespl," %.5f", prlim[i][i]); fprintf(ficrespl,"\n"); @@ -4085,6 +4122,7 @@ Interval (in months) between two waves: /* hstepm=1; aff par mois*/ + fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; @@ -4102,13 +4140,13 @@ Interval (in months) between two waves: p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); oldm=oldms;savm=savms; hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); - fprintf(ficrespij,"# Age"); + fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespij," %1d-%1d",i,j); fprintf(ficrespij,"\n"); for (h=0; h<=nhstepm; h++){ - fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); + fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespij," %.5f", p3mat[i][j][h]); @@ -4126,14 +4164,17 @@ Interval (in months) between two waves: /*---------- Forecasting ------------------*/ - if((stepm == 1) && (strcmp(model,".")==0)){ - prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); - if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1); - } - else{ - erreur=108; - printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); - fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); + /*if((stepm == 1) && (strcmp(model,".")==0)){*/ + if(prevfcast==1){ + if(stepm ==1){ + prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); + if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1); + } + else{ + erreur=108; + printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); + fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); + } } @@ -4167,9 +4208,7 @@ 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); - calagedate=-1; - - prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);