version 1.22, 2002/02/22 17:54:20
|
version 1.26, 2002/02/27 15:42:00
|
Line 1
|
Line 1
|
/* $Id$
|
/* $Id$
|
Interpolate Markov Chain
|
Interpolated Markov Chain
|
|
|
Short summary of the programme:
|
Short summary of the programme:
|
|
|
Line 95 double jmean; /* Mean space between 2 wa
|
Line 95 double jmean; /* Mean space between 2 wa
|
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;
|
FILE *ficgp, *fichtm,*ficresprob,*ficpop;
|
FILE *ficgp,*ficresprob,*ficpop;
|
FILE *ficreseij;
|
FILE *ficreseij;
|
char filerese[FILENAMELENGTH];
|
char filerese[FILENAMELENGTH];
|
FILE *ficresvij;
|
FILE *ficresvij;
|
Line 123 FILE *ficreseij;
|
Line 123 FILE *ficreseij;
|
static double maxarg1,maxarg2;
|
static double maxarg1,maxarg2;
|
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))
|
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))
|
#define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))
|
#define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))
|
|
|
#define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
|
#define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
|
#define rint(a) floor(a+0.5)
|
#define rint(a) floor(a+0.5)
|
|
|
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 ***********/
|
|
void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[] ){
|
|
int jj1, k1, i1, cpt;
|
|
FILE *fichtm;
|
|
/*char optionfilehtm[FILENAMELENGTH];*/
|
|
|
|
strcpy(optionfilehtm,optionfile);
|
|
strcat(optionfilehtm,".htm");
|
|
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
|
|
printf("Problem with %s \n",optionfilehtm), exit(0);
|
|
}
|
|
|
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">
|
|
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
|
|
|
Total number of observations=%d <br>
|
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
|
<hr size=\"2\" color=\"#EC5E5E\">
|
|
<li>Outputs files<br><br>\n
|
|
- Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
|
|
- Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
|
|
- Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
|
|
- Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
|
|
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
|
|
- Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
|
|
- Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
|
|
- Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
|
|
- Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
|
|
- Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>
|
|
<br>",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);
|
|
|
|
fprintf(fichtm," <li>Graphs</li><p>");
|
|
|
|
m=cptcoveff;
|
|
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
|
|
|
|
jj1=0;
|
|
for(k1=1; k1<=m;k1++){
|
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
|
jj1++;
|
|
if (cptcovn > 0) {
|
|
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
|
|
for (cpt=1; cpt<=cptcoveff;cpt++)
|
|
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
|
|
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
|
|
}
|
|
fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
|
|
<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
|
for(cpt=1; cpt<nlstate;cpt++){
|
|
fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
|
|
<img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
|
interval) in state (%d): v%s%d%d.gif <br>
|
|
<img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
|
|
<img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
|
}
|
|
fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
|
health expectancies in states (1) and (2): e%s%d.gif<br>
|
|
<img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
|
fprintf(fichtm,"\n</body>");
|
|
}
|
|
}
|
|
fclose(fichtm);
|
|
}
|
|
|
|
/******************* Gnuplot file **************/
|
|
void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){
|
|
|
|
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
|
|
|
|
strcpy(optionfilegnuplot,optionfilefiname);
|
|
strcat(optionfilegnuplot,".plt");
|
|
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
|
|
printf("Problem with file %s",optionfilegnuplot);
|
|
}
|
|
|
|
#ifdef windows
|
|
fprintf(ficgp,"cd \"%s\" \n",pathc);
|
|
#endif
|
|
m=pow(2,cptcoveff);
|
|
|
|
/* 1eme*/
|
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
|
|
#ifdef windows
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);
|
|
#endif
|
|
#ifdef unix
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);
|
|
#endif
|
|
|
|
for (i=1; i<= nlstate ; i ++) {
|
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$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);
|
|
for (i=1; i<= nlstate ; i ++) {
|
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
|
|
#ifdef unix
|
|
fprintf(ficgp,"\nset ter gif small size 400,300");
|
|
#endif
|
|
fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
/*2 eme*/
|
|
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
|
|
|
|
for (i=1; i<= nlstate+1 ; i ++) {
|
|
k=2*i;
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
|
|
else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"\" w l 0,");
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
|
|
else fprintf(ficgp,"\" t\"\" w l 0,");
|
|
}
|
|
fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
|
|
}
|
|
|
|
/*3eme*/
|
|
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
|
k=2+nlstate*(cpt-1);
|
|
fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);
|
|
for (i=1; i< nlstate ; i ++) {
|
|
fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);
|
|
}
|
|
fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
|
|
/* CV preval stat */
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
for (cpt=1; cpt<nlstate ; cpt ++) {
|
|
k=3;
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);
|
|
|
|
for (i=1; i< nlstate ; i ++)
|
|
fprintf(ficgp,"+$%d",k+i+1);
|
|
fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
|
|
|
|
l=3+(nlstate+ndeath)*cpt;
|
|
fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);
|
|
for (i=1; i< nlstate ; i ++) {
|
|
l=3+(nlstate+ndeath)*cpt;
|
|
fprintf(ficgp,"+$%d",l+i+1);
|
|
}
|
|
fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
|
|
fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
|
|
/* proba elementaires */
|
|
for(i=1,jk=1; i <=nlstate; i++){
|
|
for(k=1; k <=(nlstate+ndeath); k++){
|
|
if (k != i) {
|
|
for(j=1; j <=ncovmodel; j++){
|
|
|
|
fprintf(ficgp,"p%d=%f ",jk,p[jk]);
|
|
jk++;
|
|
fprintf(ficgp,"\n");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for(jk=1; jk <=m; jk++) {
|
|
fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemax);
|
|
i=1;
|
|
for(k2=1; k2<=nlstate; k2++) {
|
|
k3=i;
|
|
for(k=1; k<=(nlstate+ndeath); k++) {
|
|
if (k != k2){
|
|
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
|
|
ij=1;
|
|
for(j=3; j <=ncovmodel; j++) {
|
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
|
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
|
ij++;
|
|
}
|
|
else
|
|
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
|
}
|
|
fprintf(ficgp,")/(1");
|
|
|
|
for(k1=1; k1 <=nlstate; k1++){
|
|
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
|
|
ij=1;
|
|
for(j=3; j <=ncovmodel; j++){
|
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
|
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
|
ij++;
|
|
}
|
|
else
|
|
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
|
}
|
|
fprintf(ficgp,")");
|
|
}
|
|
fprintf(ficgp,") t \"p%d%d\" ", k2,k);
|
|
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
|
i=i+ncovmodel;
|
|
}
|
|
}
|
|
}
|
|
fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
|
|
}
|
|
|
|
fclose(ficgp);
|
|
} /* end gnuplot */
|
|
|
|
|
|
/*************** Moving average **************/
|
|
void movingaverage(double agedeb, double fage,double agemin, double ***mobaverage){
|
|
|
|
int i, cpt, cptcod;
|
|
for (agedeb=agemin; agedeb<=fage; agedeb++)
|
|
for (i=1; i<=nlstate;i++)
|
|
for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
|
|
mobaverage[(int)agedeb][i][cptcod]=0.;
|
|
|
|
for (agedeb=agemin+4; agedeb<=fage; agedeb++){
|
|
for (i=1; i<=nlstate;i++){
|
|
for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
|
for (cpt=0;cpt<=4;cpt++){
|
|
mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
|
|
}
|
|
mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
|
|
}
|
|
}
|
|
}
|
|
|
}
|
}
|
|
|
/***********************************************/
|
/***********************************************/
|
Line 1898 int main(int argc, char *argv[])
|
Line 2165 int main(int argc, char *argv[])
|
int *indx;
|
int *indx;
|
char line[MAXLINE], linepar[MAXLINE];
|
char line[MAXLINE], linepar[MAXLINE];
|
char title[MAXLINE];
|
char title[MAXLINE];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
|
char optionfilext[10], optionfilefiname[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], fileresf[FILENAMELENGTH];
|
|
|
char filerest[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
Line 2472 printf("Total number of individuals= %d,
|
Line 2739 printf("Total number of individuals= %d,
|
}
|
}
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
|
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
|
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
|
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
|
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
|
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
Line 2499 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\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
|
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\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
|
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\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
|
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);
|
/*chdir(pathcd);*/
|
|
strcpy(optionfilegnuplot,optionfilefiname);
|
/*------------ free_vector -------------*/
|
strcat(optionfilegnuplot,".plt");
|
chdir(path);
|
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
|
|
printf("Problem with file %s",optionfilegnuplot);goto end;
|
free_ivector(wav,1,imx);
|
}
|
free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
|
#ifdef windows
|
free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
|
fprintf(ficgp,"cd \"%s\" \n",pathc);
|
free_ivector(num,1,n);
|
#endif
|
free_vector(agedc,1,n);
|
m=pow(2,cptcoveff);
|
/*free_matrix(covar,1,NCOVMAX,1,n);*/
|
|
fclose(ficparo);
|
/* 1eme*/
|
fclose(ficres);
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
|
|
#ifdef windows
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);
|
|
#endif
|
|
#ifdef unix
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);
|
|
#endif
|
|
|
|
for (i=1; i<= nlstate ; i ++) {
|
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$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);
|
|
for (i=1; i<= nlstate ; i ++) {
|
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
|
|
#ifdef unix
|
|
fprintf(ficgp,"\nset ter gif small size 400,300");
|
|
#endif
|
|
fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
/*2 eme*/
|
|
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
|
|
|
|
for (i=1; i<= nlstate+1 ; i ++) {
|
|
k=2*i;
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
|
|
else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
fprintf(ficgp,"\" t\"\" w l 0,");
|
|
fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
|
|
for (j=1; j<= nlstate+1 ; j ++) {
|
|
if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
|
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
|
}
|
|
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
|
|
else fprintf(ficgp,"\" t\"\" w l 0,");
|
|
}
|
|
fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
|
|
}
|
|
|
|
/*3eme*/
|
|
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
|
k=2+nlstate*(cpt-1);
|
|
fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);
|
|
for (i=1; i< nlstate ; i ++) {
|
|
fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);
|
|
}
|
|
fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
|
|
/* CV preval stat */
|
|
for (k1=1; k1<= m ; k1 ++) {
|
|
for (cpt=1; cpt<nlstate ; cpt ++) {
|
|
k=3;
|
|
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);
|
|
for (i=1; i< nlstate ; i ++)
|
|
fprintf(ficgp,"+$%d",k+i+1);
|
|
fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
|
|
|
|
l=3+(nlstate+ndeath)*cpt;
|
|
fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);
|
|
for (i=1; i< nlstate ; i ++) {
|
|
l=3+(nlstate+ndeath)*cpt;
|
|
fprintf(ficgp,"+$%d",l+i+1);
|
|
}
|
|
fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
|
|
fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
|
}
|
|
}
|
|
|
|
/* proba elementaires */
|
|
for(i=1,jk=1; i <=nlstate; i++){
|
|
for(k=1; k <=(nlstate+ndeath); k++){
|
|
if (k != i) {
|
|
for(j=1; j <=ncovmodel; j++){
|
|
/*fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/
|
|
/*fprintf(ficgp,"%s",alph[1]);*/
|
|
fprintf(ficgp,"p%d=%f ",jk,p[jk]);
|
|
jk++;
|
|
fprintf(ficgp,"\n");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for(jk=1; jk <=m; jk++) {
|
|
fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemax);
|
|
i=1;
|
|
for(k2=1; k2<=nlstate; k2++) {
|
|
k3=i;
|
|
for(k=1; k<=(nlstate+ndeath); k++) {
|
|
if (k != k2){
|
|
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
|
|
ij=1;
|
|
for(j=3; j <=ncovmodel; j++) {
|
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
|
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
|
ij++;
|
|
}
|
|
else
|
|
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
|
}
|
|
fprintf(ficgp,")/(1");
|
|
|
|
for(k1=1; k1 <=nlstate; k1++){
|
|
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
|
|
ij=1;
|
|
for(j=3; j <=ncovmodel; j++){
|
|
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
|
|
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
|
|
ij++;
|
|
}
|
|
else
|
|
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
|
|
}
|
|
fprintf(ficgp,")");
|
|
}
|
|
fprintf(ficgp,") t \"p%d%d\" ", k2,k);
|
|
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
|
i=i+ncovmodel;
|
|
}
|
|
}
|
|
}
|
|
fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
|
|
}
|
|
|
|
fclose(ficgp);
|
|
/* end gnuplot */
|
|
|
|
chdir(path);
|
|
|
|
free_ivector(wav,1,imx);
|
|
free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
|
|
free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
|
|
free_ivector(num,1,n);
|
|
free_vector(agedc,1,n);
|
|
/*free_matrix(covar,1,NCOVMAX,1,n);*/
|
|
fclose(ficparo);
|
|
fclose(ficres);
|
|
/* }*/
|
|
|
|
/*________fin mle=1_________*/
|
|
|
|
|
|
|
|
/* No more information from the sample is required now */
|
|
/* Reads comments: lines beginning with '#' */
|
/* Reads comments: lines beginning with '#' */
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
Line 2704 chdir(path);
|
Line 2801 chdir(path);
|
fprintf(ficparo,"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 --------*/
|
/*--------- index.htm --------*/
|
|
|
strcpy(optionfilehtm,optionfile);
|
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);
|
strcat(optionfilehtm,".htm");
|
|
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
|
|
printf("Problem with %s \n",optionfilehtm);goto end;
|
|
}
|
|
|
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">
|
|
Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
|
Total number of observations=%d <br>
|
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
|
<hr size=\"2\" color=\"#EC5E5E\">
|
|
<li>Outputs files<br><br>\n
|
|
- Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
|
|
- Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
|
|
- Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
|
|
- Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
|
|
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
|
|
- Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
|
|
- Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
|
|
- Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
|
|
- Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
|
|
- Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>
|
|
<br>",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);
|
|
|
|
fprintf(fichtm," <li>Graphs</li><p>");
|
|
|
|
m=cptcoveff;
|
|
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
|
|
|
|
j1=0;
|
|
for(k1=1; k1<=m;k1++){
|
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
|
j1++;
|
|
if (cptcovn > 0) {
|
|
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
|
|
for (cpt=1; cpt<=cptcoveff;cpt++)
|
|
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);
|
|
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
|
|
}
|
|
fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
|
|
<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);
|
|
for(cpt=1; cpt<nlstate;cpt++){
|
|
fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
|
|
<img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
|
interval) in state (%d): v%s%d%d.gif <br>
|
|
<img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
|
|
<img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);
|
|
}
|
|
fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
|
health expectancies in states (1) and (2): e%s%d.gif<br>
|
|
<img src=\"e%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);
|
|
fprintf(fichtm,"\n</body>");
|
|
}
|
|
}
|
|
fclose(fichtm);
|
|
|
|
|
|
/*--------------- Prevalence limit --------------*/
|
/*--------------- Prevalence limit --------------*/
|
|
|
strcpy(filerespl,"pl");
|
strcpy(filerespl,"pl");
|
Line 2868 fclose(fichtm);
|
Line 2906 fclose(fichtm);
|
/*---------- 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;
|
|
|
if (mobilav==1) {
|
if (mobilav==1) {
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
|
movingaverage(agedeb, fage, agemin, mobaverage);
|
for (i=1; i<=nlstate;i++)
|
|
for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
|
|
mobaverage[(int)agedeb][i][cptcod]=0.;
|
|
|
|
for (agedeb=bage+4; agedeb<=fage; agedeb++){
|
|
for (i=1; i<=nlstate;i++){
|
|
for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
|
for (cpt=0;cpt<=4;cpt++){
|
|
mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
|
|
}
|
|
mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
|
|
}
|
|
}
|
|
}
|
|
}
|
}
|
|
|
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 2921 fclose(fichtm);
|
Line 2943 fclose(fichtm);
|
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 2931 fclose(fichtm);
|
Line 2953 fclose(fichtm);
|
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 2941 fclose(fichtm);
|
Line 2963 fclose(fichtm);
|
|
|
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 2953 fclose(fichtm);
|
Line 2975 fclose(fichtm);
|
fprintf(ficresf,"# StartingAge FinalAge");
|
fprintf(ficresf,"# StartingAge FinalAge");
|
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<4;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>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */
|
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
|
nhstepm = nhstepm/hstepm;
|
|
/*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/
|
|
|
|
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 (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
|
if (h==(int) (calagedate+YEARM*cpt)) {
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
|
nhstepm = nhstepm/hstepm;
|
}
|
|
for(j=1; j<=nlstate+ndeath;j++) {
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
kk1=0.;kk2=0;
|
oldm=oldms;savm=savms;
|
for(i=1; i<=nlstate;i++) {
|
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
|
if (mobilav==1)
|
|
kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
|
for (h=0; h<=nhstepm; h++){
|
else {
|
if (h==(int) (calagedate+YEARM*cpt)) {
|
kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
|
fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
|
/* fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h], probs[(int)(agedeb)+1][i][cptcod]);*/
|
}
|
|
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)){
|
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) {
|
free_ivector(popage,0,AGESUP);
|
free_ivector(popage,0,AGESUP);
|
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 ------------*/
|
|
|
strcpy(filerest,"t");
|
strcpy(filerest,"t");
|
Line 3059 fclose(fichtm);
|
Line 3083 fclose(fichtm);
|
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 3103 fclose(fichtm);
|
Line 3129 fclose(fichtm);
|
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 3145 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)
|