--- imach/src/imach.c 2002/02/26 17:11:54 1.25 +++ imach/src/imach.c 2002/02/27 15:42:00 1.26 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.25 2002/02/26 17:11:54 lievre Exp $ +/* $Id: imach.c,v 1.26 2002/02/27 15:42:00 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -1176,7 +1176,7 @@ void lubksb(double **a, int n, int *indx } /************ Frequencies ********************/ -void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2) +void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) { /* Some frequencies */ int i, m, jk, k1,i1, j1, bool, z1,z2,j; @@ -1236,6 +1236,9 @@ void freqsummary(char fileres[], int ag } } } + + fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -1620,7 +1623,7 @@ void varevsij(char fileres[], double *** for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; } - + for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) @@ -1632,7 +1635,7 @@ void varevsij(char fileres[], double *** 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); - + if (popbased==1) { for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; @@ -1684,7 +1687,7 @@ void varevsij(char fileres[], double *** free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } /* End age */ - + free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,npar); free_matrix(dnewm,1,nlstate,1,nlstate); @@ -1876,7 +1879,7 @@ if (i== 4) fprintf(ficresprob,"%.3e %.3e } free_vector(xp,1,npar); fclose(ficresprob); - exit(0); + } /******************* Printing html file ***********/ @@ -2763,11 +2766,11 @@ printf("Total number of individuals= %d, fputs(line,ficparo); } ungetc(c,ficpar); - fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mob_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav); -fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); -fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); + fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav); +fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); +fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p); @@ -2903,23 +2906,20 @@ fprintf(ficres,"popforecast=%d popfile=% /*---------- Forecasting ------------------*/ calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; - /*printf("calage= %f", calagedate);*/ + prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); - + strcpy(fileresf,"f"); strcat(fileresf,fileres); if((ficresf=fopen(fileresf,"w"))==NULL) { printf("Problem with forecast resultfile: %s\n", fileresf);goto end; } printf("Computing forecasting: result on file '%s' \n", fileresf); - + free_matrix(mint,1,maxwav,1,n); free_matrix(anint,1,maxwav,1,n); - free_matrix(agev,1,maxwav,1,imx); - - /* Mobile average */ if (cptcoveff==0) ncodemax[cptcoveff]=1; @@ -2930,11 +2930,11 @@ fprintf(ficres,"popforecast=%d popfile=% stepsize=(int) (stepm+YEARM-1)/YEARM; if (stepm<=12) stepsize=1; - + agelim=AGESUP; - /*hstepm=stepsize*YEARM; *//* Every year of age */ + hstepm=1; - hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ + hstepm=hstepm/stepm; yp1=modf(dateintmean,&yp); anprojmean=yp; yp2=modf((yp1*12),&yp); @@ -2943,9 +2943,9 @@ fprintf(ficres,"popforecast=%d popfile=% jprojmean=yp; if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; - + fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); - + if (popforecast==1) { if((ficpop=fopen(popfile,"r"))==NULL) { printf("Problem with population file : %s\n",popfile);goto end; @@ -2953,7 +2953,7 @@ fprintf(ficres,"popforecast=%d popfile=% popage=ivector(0,AGESUP); popeffectif=vector(0,AGESUP); popcount=vector(0,AGESUP); - + i=1; while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) { @@ -2963,7 +2963,7 @@ fprintf(ficres,"popforecast=%d popfile=% for (i=1; i=(agemin-((int)calagedate %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*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); + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ + nhstepm=(int) rint((agelim-agedeb)*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); - for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedate+YEARM*cpt)) { - fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); - } - for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; - for(i=1; i<=nlstate;i++) { - if (mobilav==1) - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; - else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; - + for (h=0; h<=nhstepm; h++){ + if (h==(int) (calagedate+YEARM*cpt)) { + fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); + } + for(j=1; j<=nlstate+ndeath;j++) { + kk1=0.;kk2=0; + for(i=1; i<=nlstate;i++) { + if (mobilav==1) + kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; + else { + kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; + } + if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb]; } - - if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb]; - } - - if (h==(int)(calagedate+12*cpt)){ - fprintf(ficresf," %.3f", kk1); + if (h==(int)(calagedate+12*cpt)){ + fprintf(ficresf," %.3f", kk1); - if (popforecast==1) fprintf(ficresf," [%.f]", kk2); + if (popforecast==1) fprintf(ficresf," [%.f]", kk2); + } } } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } - } + } } } + if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); if (popforecast==1) { @@ -3024,16 +3022,18 @@ fprintf(ficres,"popforecast=%d popfile=% free_vector(popeffectif,0,AGESUP); free_vector(popcount,0,AGESUP); } + free_imatrix(s,1,maxwav+1,1,n); free_vector(weight,1,n); fclose(ficresf); - }/* End forecasting */ + } + + /* End forecasting */ else{ erreur=108; printf("Error %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm); - } - + /*---------- Health expectancies and variances ------------*/ @@ -3083,12 +3083,14 @@ fprintf(ficres,"popforecast=%d popfile=% evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); - + varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); + + + fprintf(ficrest,"#Total LEs with variances: e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); - + hf=1; if (stepm >= YEARM) hf=stepm/YEARM; epj=vector(1,nlstate+1); @@ -3127,7 +3129,7 @@ fprintf(ficres,"popforecast=%d popfile=% fclose(ficpar); free_vector(epj,1,nlstate+1); /* scanf("%d ",i); */ - + /*------- Variance limit prevalence------*/ strcpy(fileresvpl,"vpl"); @@ -3169,7 +3171,7 @@ strcpy(fileresvpl,"vpl"); free_matrix(matcov,1,npar,1,npar); free_vector(delti,1,npar); - + free_matrix(agev,1,maxwav,1,imx); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); if(erreur >0)