--- imach/src/imach.c 2003/02/04 12:40:59 1.68 +++ imach/src/imach.c 2003/03/28 13:32:54 1.71 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.68 2003/02/04 12:40:59 lievre Exp $ +/* $Id: imach.c,v 1.71 2003/03/28 13:32:54 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.91, 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; @@ -972,19 +972,38 @@ 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 */ + if( s2 > nlstate){ + /* i.e. if s2 is a death state and if the date of death is known then the contribution + to the likelihood is the probability to die between last step unit time and current + step unit time, which is also the differences between probability to die before dh + and probability to die before dh-stepm . + In version up to 0.92 likelihood was computed + as if date of death was unknown. Death was treated as any other + health state: the date of the interview describes the actual state + and not the date of a change in health state. The former idea was + to consider that at each interview the state was recorded + (healthy, disable or death) and IMaCh was corrected; but when we + introduced the exact date of death then we should have modified + the contribution of an exact death to the likelihood. This new + contribution is smaller and very dependent of the step unit + stepm. It is no more the probability to die between last interview + and month of death but the probability to survive from last + interview up to one month before death multiplied by the + probability to die within a month. Thanks to Chris + Jackson for correcting this bug. Former versions increased + mortality artificially. The bad side is that we add another loop + which slows down the processing. The difference can be up to 10% + lower mortality. + */ + lli=log(out[s1][s2] - savm[s1][s2]); + }else{ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /* 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)*/ /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ - ipmx +=1; + ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; } /* end of wave */ @@ -1424,7 +1443,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 +1558,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 +1576,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, 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); - 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 +1608,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 +1621,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++){ @@ -1642,7 +1663,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 +1693,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; @@ -1713,6 +1734,7 @@ void concatwav(int wav[], int **dh, int 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);*/ + /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ } } else{ @@ -1750,7 +1772,7 @@ void concatwav(int wav[], int **dh, int if(dh[mi][i]==0){ dh[mi][i]=1; /* At least one step */ 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); + /* 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);*/ } } } /* end if mle */ @@ -2144,7 +2166,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 +2197,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 +2222,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++) @@ -2218,16 +2240,13 @@ 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); - - for(j=nlstate+1;j<=nlstate+ndeath;j++) - for(i=nlstate+1;i<=nlstate+ndeath;i++){ + 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); @@ -2242,7 +2261,7 @@ void varevsij(char optionfilefiname[], d prlim[i][i]=mobaverage[(int)age][i][ij]; } } - + /* This for computing probability of death (h=1 means computed over hstepm (estepm) matrices product = hstepm*stepm months) as a weighted average of prlim. @@ -2288,10 +2307,10 @@ void varevsij(char optionfilefiname[], d 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", estepm,digitp,digit); + fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,digitp,optionfilefiname,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); + fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,nlstate); @@ -2303,7 +2322,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) @@ -2449,10 +2468,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 +2511,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 +2527,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 +2846,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)"); @@ -3022,21 +3042,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); @@ -3060,11 +3082,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; @@ -3073,49 +3094,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-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; - 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; - for(i=1; i<=nlstate;i++) { + 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][i][cptcod]; else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; + ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][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); } } } @@ -3129,12 +3155,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; @@ -3143,9 +3170,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"); @@ -3191,7 +3218,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 +3242,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 +3254,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 +3265,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 +3279,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 +3287,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 +3295,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 +3342,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; @@ -3615,7 +3643,11 @@ int main(int argc, char *argv[]) if (s[4][i]==9) s[4][i]=-1; printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/ + for (i=1; i<=imx; i++) + /*if ((s[3][i]==3) || (s[4][i]==3)) weight[i]=0.08; + else weight[i]=1;*/ + /* Calculation of the number of parameter from char model*/ Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ Tprod=ivector(1,15); @@ -3717,11 +3749,12 @@ 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){ + for(m=firstpass; (m<= lastpass); m++){ + 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){ @@ -3731,7 +3764,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; @@ -3757,7 +3792,7 @@ int main(int argc, char *argv[]) } for (i=1; i<=imx; i++) { - for(m=1; (m<= maxwav); m++){ + for(m=firstpass; (m<=lastpass); m++){ if (s[m][i] > (nlstate+ndeath)) { printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); @@ -3766,6 +3801,13 @@ int main(int argc, char *argv[]) } } + /*for (i=1; i<=imx; i++){ + for (m=firstpass; (m