version 1.25, 2002/02/26 17:11:54
|
version 1.27, 2002/02/28 17:49:07
|
Line 94 int **dh; /* dh[mi][i] is number of step
|
Line 94 int **dh; /* dh[mi][i] is number of step
|
double jmean; /* Mean space between 2 waves */
|
double jmean; /* Mean space between 2 waves */
|
double **oldm, **newm, **savm; /* Working pointers to matrices */
|
double **oldm, **newm, **savm; /* Working pointers to matrices */
|
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
|
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf;
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
|
FILE *ficgp,*ficresprob,*ficpop;
|
FILE *ficgp,*ficresprob,*ficpop;
|
FILE *ficreseij;
|
FILE *ficreseij;
|
char filerese[FILENAMELENGTH];
|
char filerese[FILENAMELENGTH];
|
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 2143 void movingaverage(double agedeb, double
|
Line 2146 void movingaverage(double agedeb, double
|
|
|
}
|
}
|
|
|
|
|
|
/************** Forecasting ******************/
|
|
prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, 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;
|
|
int *popage;
|
|
double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
|
|
double *popeffectif,*popcount;
|
|
double ***p3mat;
|
|
char fileresf[FILENAMELENGTH];
|
|
|
|
agelim=AGESUP;
|
|
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
|
|
|
|
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);
|
|
}
|
|
printf("Computing forecasting: result on file '%s' \n", fileresf);
|
|
|
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
|
|
|
if (mobilav==1) {
|
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
movingaverage(agedeb, fage, agemin, mobaverage);
|
|
}
|
|
|
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
|
if (stepm<=12) stepsize=1;
|
|
|
|
agelim=AGESUP;
|
|
|
|
hstepm=1;
|
|
hstepm=hstepm/stepm;
|
|
yp1=modf(dateintmean,&yp);
|
|
anprojmean=yp;
|
|
yp2=modf((yp1*12),&yp);
|
|
mprojmean=yp;
|
|
yp1=modf((yp2*30.5),&yp);
|
|
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);
|
|
|
|
for(cptcov=1;cptcov<=i2;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,"******\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,"\n");
|
|
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);
|
|
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];
|
|
}
|
|
|
|
}
|
|
if (h==(int)(calagedate+12*cpt)){
|
|
fprintf(ficresf," %.3f", kk1);
|
|
|
|
}
|
|
}
|
|
}
|
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
|
|
fclose(ficresf);
|
|
}
|
|
/************** Forecasting ******************/
|
|
populforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, 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;
|
|
int *popage;
|
|
double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
|
|
double *popeffectif,*popcount;
|
|
double ***p3mat,***tabpop,***tabpopprev;
|
|
char filerespop[FILENAMELENGTH];
|
|
|
|
tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
agelim=AGESUP;
|
|
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
|
|
|
|
prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
|
|
|
strcpy(filerespop,"pop");
|
|
strcat(filerespop,fileres);
|
|
if((ficrespop=fopen(filerespop,"w"))==NULL) {
|
|
printf("Problem with forecast resultfile: %s\n", filerespop);
|
|
}
|
|
printf("Computing forecasting: result on file '%s' \n", filerespop);
|
|
|
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
|
|
|
if (mobilav==1) {
|
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
movingaverage(agedeb, fage, agemin, mobaverage);
|
|
}
|
|
|
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
|
if (stepm<=12) stepsize=1;
|
|
|
|
agelim=AGESUP;
|
|
|
|
hstepm=1;
|
|
hstepm=hstepm/stepm;
|
|
|
|
if (popforecast==1) {
|
|
if((ficpop=fopen(popfile,"r"))==NULL) {
|
|
printf("Problem with population file : %s\n",popfile);exit(0);
|
|
}
|
|
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) i=i+1;
|
|
|
|
imx=i;
|
|
for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
|
|
}
|
|
|
|
for(cptcov=1;cptcov<=i2;cptcov++){
|
|
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
|
k=k+1;
|
|
fprintf(ficrespop,"\n#******");
|
|
for(j=1;j<=cptcoveff;j++) {
|
|
fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
|
}
|
|
fprintf(ficrespop,"******\n");
|
|
fprintf(ficrespop,"# StartingAge FinalAge");
|
|
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);
|
|
if (popforecast==1) fprintf(ficrespop," [Population]");
|
|
|
|
for (cpt=0; cpt<=0;cpt++) {
|
|
fprintf(ficrespop,"\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);
|
|
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(ficrespop,"\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 (h==(int)(calagedate+12*cpt)){
|
|
tabpop[(int)(agedeb)][j][cptcod]=kk1;
|
|
/*fprintf(ficrespop," %.3f", kk1);
|
|
if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
|
|
}
|
|
}
|
|
for(i=1; i<=nlstate;i++){
|
|
kk1=0.;
|
|
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)];
|
|
}
|
|
|
|
if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %.3f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
|
|
}
|
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
|
}
|
|
}
|
|
|
|
/******/
|
|
|
|
for (cpt=1; cpt<=4;cpt++) {
|
|
fprintf(ficrespop,"\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);
|
|
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++) {
|
|
kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];
|
|
}
|
|
if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %.3f", kk1);
|
|
}
|
|
}
|
|
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) {
|
|
free_ivector(popage,0,AGESUP);
|
|
free_vector(popeffectif,0,AGESUP);
|
|
free_vector(popcount,0,AGESUP);
|
|
}
|
|
free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
fclose(ficrespop);
|
|
}
|
|
|
/***********************************************/
|
/***********************************************/
|
/**************** Main Program *****************/
|
/**************** Main Program *****************/
|
/***********************************************/
|
/***********************************************/
|
Line 2165 int main(int argc, char *argv[])
|
Line 2423 int main(int argc, char *argv[])
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
|
char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
|
char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
|
|
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
|
|
|
char filerest[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
Line 2179 int main(int argc, char *argv[])
|
Line 2437 int main(int argc, char *argv[])
|
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
int mobilav=0,popforecast=0;
|
int mobilav=0,popforecast=0;
|
int hstepm, nhstepm;
|
int hstepm, nhstepm;
|
int *popage;/*boolprev=0 if date and zero if wave*/
|
|
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
|
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
|
|
|
double bage, fage, age, agelim, agebase;
|
double bage, fage, age, agelim, agebase;
|
Line 2195 int main(int argc, char *argv[])
|
Line 2452 int main(int argc, char *argv[])
|
double **varpl; /* Variances of prevalence limits by age */
|
double **varpl; /* Variances of prevalence limits by age */
|
double *epj, vepp;
|
double *epj, vepp;
|
double kk1, kk2;
|
double kk1, kk2;
|
double *popeffectif,*popcount;
|
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
|
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
|
|
double yp,yp1,yp2;
|
|
|
|
char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
|
char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
Line 2763 printf("Total number of individuals= %d,
|
Line 3019 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 2899 fprintf(ficres,"popforecast=%d popfile=%
|
Line 3155 fprintf(ficres,"popforecast=%d popfile=%
|
|
|
fclose(ficrespij);
|
fclose(ficrespij);
|
|
|
if(stepm == 1) {
|
|
/*---------- 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;
|
|
|
|
if (mobilav==1) {
|
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
movingaverage(agedeb, fage, agemin, mobaverage);
|
|
}
|
|
|
|
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 */
|
|
yp1=modf(dateintmean,&yp);
|
|
anprojmean=yp;
|
|
yp2=modf((yp1*12),&yp);
|
|
mprojmean=yp;
|
|
yp1=modf((yp2*30.5),&yp);
|
|
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;
|
|
}
|
|
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)
|
|
{
|
|
i=i+1;
|
|
}
|
|
imx=i;
|
|
|
|
for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
|
|
}
|
|
|
|
for(cptcov=1;cptcov<=i1;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,"******\n");
|
|
fprintf(ficresf,"# StartingAge FinalAge");
|
|
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
|
|
if (popforecast==1) fprintf(ficresf," [Population]");
|
|
|
|
for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {
|
|
fprintf(ficresf,"\n");
|
|
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);
|
|
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];
|
|
|
|
}
|
|
|
|
if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
|
|
}
|
|
|
|
if (h==(int)(calagedate+12*cpt)){
|
|
fprintf(ficresf," %.3f", kk1);
|
|
|
|
if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
|
|
}
|
|
}
|
|
}
|
|
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) {
|
/*---------- Forecasting ------------------*/
|
free_ivector(popage,0,AGESUP);
|
if(stepm == 1) {
|
free_vector(popeffectif,0,AGESUP);
|
prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
|
free_vector(popcount,0,AGESUP);
|
populforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
|
}
|
free_matrix(mint,1,maxwav,1,n);
|
free_imatrix(s,1,maxwav+1,1,n);
|
free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
|
free_vector(weight,1,n);
|
free_vector(weight,1,n);}
|
fclose(ficresf);
|
|
}/* 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 ------------*/
|
|
|
strcpy(filerest,"t");
|
strcpy(filerest,"t");
|
Line 3083 fprintf(ficres,"popforecast=%d popfile=%
|
Line 3217 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 3117 fprintf(ficres,"popforecast=%d popfile=%
|
Line 3253 fprintf(ficres,"popforecast=%d popfile=%
|
}
|
}
|
}
|
}
|
}
|
}
|
|
|
|
|
|
|
|
|
fclose(ficreseij);
|
fclose(ficreseij);
|
fclose(ficresvij);
|
fclose(ficresvij);
|
fclose(ficrest);
|
fclose(ficrest);
|
fclose(ficpar);
|
fclose(ficpar);
|
free_vector(epj,1,nlstate+1);
|
free_vector(epj,1,nlstate+1);
|
/* scanf("%d ",i); */
|
|
|
|
/*------- Variance limit prevalence------*/
|
/*------- Variance limit prevalence------*/
|
|
|
strcpy(fileresvpl,"vpl");
|
strcpy(fileresvpl,"vpl");
|
strcat(fileresvpl,fileres);
|
strcat(fileresvpl,fileres);
|
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
|
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
|
printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
|
printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
|
Line 3138 strcpy(fileresvpl,"vpl");
|
Line 3270 strcpy(fileresvpl,"vpl");
|
}
|
}
|
printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
|
printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
|
|
|
k=0;
|
k=0;
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
k=k+1;
|
k=k+1;
|
fprintf(ficresvpl,"\n#****** ");
|
fprintf(ficresvpl,"\n#****** ");
|
for(j=1;j<=cptcoveff;j++)
|
for(j=1;j<=cptcoveff;j++)
|
fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
fprintf(ficresvpl,"******\n");
|
fprintf(ficresvpl,"******\n");
|
|
|
varpl=matrix(1,nlstate,(int) bage, (int) fage);
|
varpl=matrix(1,nlstate,(int) bage, (int) fage);
|
oldm=oldms;savm=savms;
|
oldm=oldms;savm=savms;
|
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
|
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
|
}
|
}
|
}
|
}
|
|
|
fclose(ficresvpl);
|
fclose(ficresvpl);
|
Line 3169 strcpy(fileresvpl,"vpl");
|
Line 3301 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)
|