--- imach/src/imach.c 2002/02/26 17:11:54 1.25
+++ imach/src/imach.c 2002/03/01 17:56:23 1.28
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.25 2002/02/26 17:11:54 lievre Exp $
+/* $Id: imach.c,v 1.28 2002/03/01 17:56:23 lievre Exp $
Interpolated Markov Chain
Short summary of the programme:
@@ -94,7 +94,7 @@ int **dh; /* dh[mi][i] is number of step
double jmean; /* Mean space between 2 waves */
double **oldm, **newm, **savm; /* 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 *ficreseij;
char filerese[FILENAMELENGTH];
@@ -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;
@@ -1185,7 +1185,7 @@ void freqsummary(char fileres[], int ag
double pos, k2, dateintsum=0,k2cpt=0;
FILE *ficresp;
char fileresp[FILENAMELENGTH];
-
+
pp=vector(1,nlstate);
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
strcpy(fileresp,"p");
@@ -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]]);
@@ -1304,7 +1307,7 @@ void freqsummary(char fileres[], int ag
}
/************ Prevalence ********************/
-void prevalence(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 calagedate)
+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 */
int i, m, jk, k1, i1, j1, bool, z1,z2,j;
@@ -1329,14 +1332,14 @@ void prevalence(int agemin, int agemax,
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++) {
bool=1;
if (cptcovn>0) {
for (z1=1; z1<=cptcoveff; z1++)
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
bool=0;
- }
+ }
if (bool==1) {
for(m=firstpass; m<=lastpass; m++){
k2=anint[m][i]+(mint[m][i]/12.);
@@ -1344,12 +1347,11 @@ void prevalence(int agemin, int agemax,
if(agev[m][i]==0) agev[m][i]=agemax+1;
if(agev[m][i]==1) agev[m][i]=agemax+2;
freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
- freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];
+ /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; */
}
}
}
}
-
for(i=(int)agemin; i <= (int)agemax+3; i++){
for(jk=1; jk <=nlstate ; jk++){
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
@@ -1620,7 +1622,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 +1634,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 +1686,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 +1878,7 @@ if (i== 4) fprintf(ficresprob,"%.3e %.3e
}
free_vector(xp,1,npar);
fclose(ficresprob);
- exit(0);
+
}
/******************* Printing html file ***********/
@@ -1907,8 +1909,9 @@ Interval (in months) between two waves:
- Variances of life expectancies by age and initial health status: v%s
- Health expectancies with their variances: t%s
- Standard deviation of stationary prevalences: vpl%s
- - Prevalences and population forecasting: f%s
-
",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
+ - Prevalences forecasting: f%s
+ - Population forecasting (if popforecast=1): pop%s
+
",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
fprintf(fichtm," Graphs");
@@ -1950,7 +1953,7 @@ fclose(fichtm);
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){
+void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemaxpar, double fage , char pathc[], double p[]){
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
@@ -2045,7 +2048,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt=(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 anpyram,double mpyram,double jpyram,double agemin, 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, 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=(anpyram+mpyram/12.+jpyram/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=(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 %3.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," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
+ }
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ }
+ }
+
+ /******/
+
+ 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>=(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 %3.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," %15.2f", 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 *****************/
/***********************************************/
@@ -2152,7 +2411,7 @@ int main(int argc, char *argv[])
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
double agedeb, agefin,hf;
- double agemin=1.e20, agemax=-1.e20;
+ double agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
double fret;
double **xi,tmp,delta;
@@ -2165,7 +2424,7 @@ int main(int argc, char *argv[])
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[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 fileregp[FILENAMELENGTH];
@@ -2179,8 +2438,7 @@ int main(int argc, char *argv[])
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
int mobilav=0,popforecast=0;
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,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;
double bage, fage, age, agelim, agebase;
double ftolpl=FTOL;
@@ -2195,9 +2453,8 @@ int main(int argc, char *argv[])
double **varpl; /* Variances of prevalence limits by age */
double *epj, vepp;
double kk1, kk2;
- double *popeffectif,*popcount;
- double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
- double yp,yp1,yp2;
+ double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
+
char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
char *alph[]={"a","a","b","c","d","e"}, str[4];
@@ -2717,16 +2974,16 @@ printf("Total number of individuals= %d,
}
ungetc(c,ficpar);
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);
+ fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemaxpar, &bage, &fage);
if (fage <= 2) {
bage = agemin;
- fage = agemax;
+ fage = agemaxpar;
}
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
- fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
+ fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
@@ -2753,9 +3010,9 @@ printf("Total number of individuals= %d,
dateprev2=anprev2+mprev2/12.+jprev2/365.;
fscanf(ficpar,"pop_based=%d\n",&popbased);
- fprintf(ficparo,"pop_based=%d\n",popbased);
- fprintf(ficres,"pop_based=%d\n",popbased);
-
+ fprintf(ficparo,"pop_based=%d\n",popbased);
+ fprintf(ficres,"pop_based=%d\n",popbased);
+
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
@@ -2763,14 +3020,28 @@ 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);
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
+ fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
+fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
+fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
+
+
+while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ puts(line);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
+
+ fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1);
+ 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,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
/*------------ gnuplot -------------*/
- printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);
+ printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemaxpar,fage, pathc,p);
/*------------ free_vector -------------*/
chdir(path);
@@ -2783,19 +3054,7 @@ fprintf(ficres,"popforecast=%d popfile=%
/*free_matrix(covar,1,NCOVMAX,1,n);*/
fclose(ficparo);
fclose(ficres);
-
- /* Reads comments: lines beginning with '#' */
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
- puts(line);
- fputs(line,ficparo);
- }
- ungetc(c,ficpar);
-
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);
- printf("agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax, bage, fage);
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
+
/*--------- index.htm --------*/
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);
@@ -2822,7 +3081,7 @@ fprintf(ficres,"popforecast=%d popfile=%
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
k=0;
agebase=agemin;
- agelim=agemax;
+ agelim=agemaxpar;
ftolpl=1.e-10;
i1=cptcoveff;
if (cptcovn < 1){i1=1;}
@@ -2899,142 +3158,20 @@ fprintf(ficres,"popforecast=%d popfile=%
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=(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) {
- free_ivector(popage,0,AGESUP);
- 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 */
+ /*---------- Forecasting ------------------*/
+ if(stepm == 1) {
+ 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);
+ free_matrix(mint,1,maxwav,1,n);
+ free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
+ free_vector(weight,1,n);}
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 ------------*/
strcpy(filerest,"t");
@@ -3083,12 +3220,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);
@@ -3117,20 +3256,16 @@ fprintf(ficres,"popforecast=%d popfile=%
}
}
}
-
-
-
- fclose(ficreseij);
- fclose(ficresvij);
+ fclose(ficreseij);
+ fclose(ficresvij);
fclose(ficrest);
fclose(ficpar);
free_vector(epj,1,nlstate+1);
- /* scanf("%d ",i); */
-
+
/*------- Variance limit prevalence------*/
-strcpy(fileresvpl,"vpl");
+ strcpy(fileresvpl,"vpl");
strcat(fileresvpl,fileres);
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
@@ -3138,19 +3273,19 @@ strcpy(fileresvpl,"vpl");
}
printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
- k=0;
- for(cptcov=1;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- fprintf(ficresvpl,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficresvpl,"******\n");
-
- varpl=matrix(1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
+ k=0;
+ for(cptcov=1;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+ k=k+1;
+ fprintf(ficresvpl,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficresvpl,"******\n");
+
+ varpl=matrix(1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
- }
+ }
}
fclose(ficresvpl);
@@ -3169,7 +3304,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)