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