--- imach/src/imach.c 2003/02/04 12:40:59 1.68 +++ imach/src/imach.c 2003/02/04 20:55:29 1.69 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.68 2003/02/04 12:40:59 lievre Exp $ +/* $Id: imach.c,v 1.69 2003/02/04 20:55:29 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -972,14 +972,6 @@ double func( double *x) * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ - /* if s2=-2 lli=out[1][1]+out[1][2];*/ - if (s2==-2) { - for (j=1,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= survp; - } - - else lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ /*if(lli ==000.0)*/ @@ -1424,7 +1416,7 @@ void freqsummary(char fileres[], int ag fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); + freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1539,8 +1531,8 @@ void freqsummary(char fileres[], int ag } } - for(jk=-2; jk <=nlstate+ndeath; jk++) - for(m=-2; m <=nlstate+ndeath; m++) + for(jk=-1; jk <=nlstate+ndeath; jk++) + for(m=-1; m <=nlstate+ndeath; m++) if(freq[jk][m][i] !=0 ) { if(first==1) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); @@ -1557,24 +1549,29 @@ void freqsummary(char fileres[], int ag dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); /* End of Freq */ } /************ 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, double calagedatem) +{ + /* 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); - freq=ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); + freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1584,12 +1581,12 @@ void prevalence(int agemin, float agemax for(i1=1; i1<=ncodemax[k1];i1++){ j1++; - for (i=-2; i<=nlstate+ndeath; i++) - for (jk=-2; jk<=nlstate+ndeath; jk++) + for (i=-1; i<=nlstate+ndeath; i++) + for (jk=-1; jk<=nlstate+ndeath; jk++) 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++) @@ -1597,20 +1594,20 @@ 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]; + if (calagedatem>0) /* We compute prevalence at exact age, agev in fractional years */ + freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedatem %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)(agemax+3)] += weight[i]; } } - } + } /* end selection of waves */ } } for(i=(int)agemin; i <= (int)agemax+3; i++){ @@ -1642,7 +1639,7 @@ void prevalence(int agemin, float agemax } /* end k1 */ - free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); } /* End of Freq */ @@ -1672,7 +1669,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1 || s[m][i]==-2) + if(s[m][i]>=1) mw[++mi][i]=m; if(m >=lastpass) break; @@ -2144,7 +2141,7 @@ void varevsij(char optionfilefiname[], d computed over hstepm matrices product = hstepm*stepm months) as a weighted average of prlim. */ - for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ + 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]; } @@ -2175,7 +2172,7 @@ void varevsij(char optionfilefiname[], d computed over hstepm matrices product = hstepm*stepm months) as a weighted average of prlim. */ - for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ + 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]; } @@ -2200,9 +2197,9 @@ void varevsij(char optionfilefiname[], d trgradg[h][j][theta]=gradg[h][theta][j]; for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ - for(theta=1; theta <=npar; theta++) { + 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++) @@ -2226,8 +2223,6 @@ void varevsij(char optionfilefiname[], d for(j=nlstate+1;j<=nlstate+ndeath;j++) 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); @@ -2449,10 +2444,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); @@ -2491,14 +2487,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 "); @@ -2507,7 +2503,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 ++){ @@ -2826,12 +2822,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)"); @@ -3023,20 +3019,28 @@ 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; + /* proj1, year, month, day of starting projection + agemin, agemax range of age + dateprev1 dateprev2 range of dates during which prevalence is computed + + */ + 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 calagedatem, agelim, kk1, kk2, 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; + calagedatem=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; /* offset + from the mean date of interviews between dateprev1 and dateprev2 + (in fractional years converted to fractional months) */ + /* Computing age-specific starting prevalence between exact + dateprev1 and dateprev2 (in float years) and + shifting ages from calagedatem. + */ + prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem); strcpy(fileresf,"f"); strcat(fileresf,fileres); @@ -3064,7 +3068,8 @@ prevforecast(char fileres[], double anpr 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; @@ -3073,9 +3078,9 @@ 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 ",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); - for(cptcov=1;cptcov<=i2;cptcov++){ + for(cptcov=1, k=0;cptcov<=i2;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; fprintf(ficresf,"\n#******"); @@ -3085,13 +3090,11 @@ prevforecast(char fileres[], double anpr 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++) { + for (yearp=0; yearp<=(anproj2-anproj1);yearp++) { fprintf(ficresf,"\n"); fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+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; @@ -3100,8 +3103,11 @@ prevforecast(char fileres[], double anpr hpxij(p3mat,nhstepm,agedeb,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) (calagedatem+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,agedeb+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { kk1=0.;kk2=0; @@ -3113,7 +3119,7 @@ prevforecast(char fileres[], double anpr } } - if (h==(int)(calagedate+12*cpt)){ + if (h==(int)(calagedatem+12*yearp)){ fprintf(ficresf," %.3f", kk1); } @@ -3134,7 +3140,7 @@ populforecast(char fileres[], double anp 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; @@ -3143,9 +3149,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, calagedatem); strcpy(filerespop,"pop"); @@ -3191,7 +3197,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; @@ -3215,7 +3221,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++) { @@ -3227,7 +3233,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]);*/ @@ -3238,10 +3244,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); @@ -3252,7 +3258,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; @@ -3260,7 +3266,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++) { @@ -3268,7 +3274,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); @@ -3315,9 +3321,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, calagedatem; double bage, fage, age, agelim, agebase; double ftolpl=FTOL; @@ -3718,10 +3725,11 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=1; (m<= maxwav); m++){ - if(s[m][i] >0 || s[m][i]==-2){ + if(s[m][i] >0){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) - if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i]; + if(moisdc[i]!=99 && andc[i]!=9999) + agev[m][i]=agedc[i]; /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ else { if (andc[i]!=9999){ @@ -3730,8 +3738,10 @@ int main(int argc, char *argv[]) agev[m][i]=-1; } } - } - 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; @@ -3783,8 +3793,7 @@ int main(int argc, char *argv[]) dh=imatrix(1,lastpass-firstpass+1,1,imx); bh=imatrix(1,lastpass-firstpass+1,1,imx); mw=imatrix(1,lastpass-firstpass+1,1,imx); - - + /* Concatenates waves */ concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); @@ -3951,6 +3960,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); @@ -3976,10 +3987,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); @@ -3993,7 +4006,6 @@ int main(int argc, char *argv[]) fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ @@ -4087,7 +4099,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"); @@ -4156,15 +4170,18 @@ 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, 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); } + } /*---------- Health expectancies and variances ------------*/ @@ -4188,7 +4205,6 @@ Interval (in months) between two waves: printf("Computing Health Expectancies: result on file '%s' \n", filerese); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); - strcpy(fileresv,"v"); strcat(fileresv,fileres); if((ficresvij=fopen(fileresv,"w"))==NULL) { @@ -4198,9 +4214,9 @@ 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; + calagedatem=-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, calagedatem); if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);