--- imach/src/imach.c 2002/05/21 18:44:41 1.42
+++ imach/src/imach.c 2002/06/20 14:03:39 1.49
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.42 2002/05/21 18:44:41 brouard Exp $
+/* $Id: imach.c,v 1.49 2002/06/20 14:03:39 lievre Exp $
Interpolated Markov Chain
Short summary of the programme:
@@ -75,11 +75,16 @@
#define YEARM 12. /* Number of months per year */
#define AGESUP 130
#define AGEBASE 40
+#ifdef windows
+#define DIRSEPARATOR '\\'
+#else
+#define DIRSEPARATOR '/'
+#endif
-
+char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES ";
int erreur; /* Error number */
int nvar;
-int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
+int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
@@ -96,13 +101,25 @@ double jmean; /* Mean space between 2 wa
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,*ficrespop;
-FILE *ficgp,*ficresprob,*ficpop;
+FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
+FILE *fichtm; /* Html File */
FILE *ficreseij;
- char filerese[FILENAMELENGTH];
- FILE *ficresvij;
- char fileresv[FILENAMELENGTH];
- FILE *ficresvpl;
- char fileresvpl[FILENAMELENGTH];
+char filerese[FILENAMELENGTH];
+FILE *ficresvij;
+char fileresv[FILENAMELENGTH];
+FILE *ficresvpl;
+char fileresvpl[FILENAMELENGTH];
+char title[MAXLINE];
+char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
+char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
+
+char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
+
+char filerest[FILENAMELENGTH];
+char fileregp[FILENAMELENGTH];
+char popfile[FILENAMELENGTH];
+
+char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
#define NR_END 1
#define FREE_ARG char*
@@ -161,11 +178,7 @@ static int split( char *path, char *dirc
l1 = strlen( path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
-#ifdef windows
- s = strrchr( path, '\\' ); /* find last / */
-#else
- s = strrchr( path, '/' ); /* find last / */
-#endif
+ s = strrchr( path, DIRSEPARATOR ); /* find last / */
if ( s == NULL ) { /* no directory, so use current */
#if defined(__bsd__) /* get current working directory */
extern char *getwd( );
@@ -1641,14 +1654,10 @@ void evsij(char fileres[], double ***eij
}
}
}
-
-
-
for(j=1; j<= nlstate*2; j++)
for(h=0; h<=nhstepm-1; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
}
-
}
/* End theta */
@@ -1658,15 +1667,15 @@ void evsij(char fileres[], double ***eij
for(h=0; h<=nhstepm-1; h++)
for(j=1; j<=nlstate*2;j++)
for(theta=1; theta <=npar; theta++)
- trgradg[h][j][theta]=gradg[h][theta][j];
-
+ trgradg[h][j][theta]=gradg[h][theta][j];
+
for(i=1;i<=nlstate*2;i++)
for(j=1;j<=nlstate*2;j++)
varhe[i][j][(int)age] =0.;
- printf("%d||",(int)age);fflush(stdout);
- for(h=0;h<=nhstepm-1;h++){
+ printf("%d|",(int)age);fflush(stdout);
+ for(h=0;h<=nhstepm-1;h++){
for(k=0;k<=nhstepm-1;k++){
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
@@ -1675,8 +1684,6 @@ void evsij(char fileres[], double ***eij
varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
}
}
-
-
/* Computing expectancies */
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
@@ -1702,6 +1709,8 @@ void evsij(char fileres[], double ***eij
free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
+ printf("\n");
+
free_vector(xp,1,npar);
free_matrix(dnewm,1,nlstate*2,1,npar);
free_matrix(doldm,1,nlstate*2,1,nlstate*2);
@@ -1724,7 +1733,7 @@ void varevsij(char fileres[], double ***
double age,agelim, hf;
int theta;
- fprintf(ficresvij,"# Covariances of life expectancies\n");
+ fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n");
fprintf(ficresvij,"# Age");
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
@@ -1859,7 +1868,7 @@ void varprevlim(char fileres[], double *
double age,agelim;
int theta;
- fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");
+ fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");
fprintf(ficresvpl,"# Age");
for(i=1; i<=nlstate;i++)
fprintf(ficresvpl," %1d-%1d",i,i);
@@ -1928,67 +1937,131 @@ void varprevlim(char fileres[], double *
}
/************ Variance of one-step probabilities ******************/
-void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
+void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
{
- int i, j, i1, k1, j1, z1;
- int k=0, cptcode;
+ int i, j=0, i1, k1, l1, t, tj;
+ int k2, l2, j1, z1;
+ int k=0,l, cptcode;
+ int first=1;
+ double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;
double **dnewm,**doldm;
double *xp;
double *gp, *gm;
double **gradg, **trgradg;
+ double **mu;
double age,agelim, cov[NCOVMAX];
+ double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
int theta;
char fileresprob[FILENAMELENGTH];
+ char fileresprobcov[FILENAMELENGTH];
+ char fileresprobcor[FILENAMELENGTH];
+
+ double ***varpij;
strcpy(fileresprob,"prob");
strcat(fileresprob,fileres);
if((ficresprob=fopen(fileresprob,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprob);
}
+ strcpy(fileresprobcov,"probcov");
+ strcat(fileresprobcov,fileres);
+ if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprobcov);
+ }
+ strcpy(fileresprobcor,"probcor");
+ strcat(fileresprobcor,fileres);
+ if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprobcor);
+ }
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
+ printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
+ printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
-fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
+ fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
fprintf(ficresprob,"# Age");
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=(nlstate+ndeath);j++)
- fprintf(ficresprob," p%1d-%1d (SE)",i,j);
+ fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
+ fprintf(ficresprobcov,"# Age");
+ fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
+ fprintf(ficresprobcov,"# Age");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=(nlstate+ndeath);j++){
+ fprintf(ficresprob," p%1d-%1d (SE)",i,j);
+ fprintf(ficresprobcov," p%1d-%1d ",i,j);
+ fprintf(ficresprobcor," p%1d-%1d ",i,j);
+ }
fprintf(ficresprob,"\n");
+ fprintf(ficresprobcov,"\n");
+ fprintf(ficresprobcor,"\n");
+ xp=vector(1,npar);
+ dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+ doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
+ mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
+ varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage);
+ first=1;
+ if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
+ printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
+ exit(0);
+ }
+ else{
+ fprintf(ficgp,"\n# Routine varprob");
+ }
+ if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
+ printf("Problem with html file: %s\n", optionfilehtm);
+ exit(0);
+ }
+ else{
+ fprintf(fichtm,"\n
Computing matrix of variance-covariance of step probabilities
\n");
+ fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
\n");
+ fprintf(fichtm,"\n
We have drawn x'cov-1x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis.
When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.
\n");
+ }
- xp=vector(1,npar);
- dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
- doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
-
+
cov[1]=1;
- j=cptcoveff;
- if (cptcovn<1) {j=1;ncodemax[1]=1;}
+ tj=cptcoveff;
+ if (cptcovn<1) {tj=1;ncodemax[1]=1;}
j1=0;
- for(k1=1; k1<=1;k1++){
- for(i1=1; i1<=ncodemax[k1];i1++){
- j1++;
-
- if (cptcovn>0) {
- fprintf(ficresprob, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
- fprintf(ficresprob, "**********\n#");
- }
-
+ for(t=1; t<=tj;t++){
+ for(i1=1; i1<=ncodemax[t];i1++){
+ j1++;
+
+ if (cptcovn>0) {
+ fprintf(ficresprob, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficresprob, "**********\n#");
+ fprintf(ficresprobcov, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficresprobcov, "**********\n#");
+
+ fprintf(ficgp, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficgp, "**********\n#");
+
+
+ fprintf(fichtm, "\n
********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(fichtm, "**********\n
");
+
+ fprintf(ficresprobcor, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+ fprintf(ficgp, "**********\n#");
+ }
+
for (age=bage; age<=fage; age ++){
cov[2]=age;
for (k=1; k<=cptcovn;k++) {
cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
-
}
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
for (k=1; k<=cptcovprod;k++)
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
- gradg=matrix(1,npar,1,9);
- trgradg=matrix(1,9,1,npar);
- gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
- gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
+ gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
+ trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+ gp=vector(1,(nlstate)*(nlstate+ndeath));
+ gm=vector(1,(nlstate)*(nlstate+ndeath));
for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++)
@@ -1997,7 +2070,7 @@ fprintf(ficresprob,"#One-step probabilit
pmij(pmmij,cov,ncovmodel,xp,nlstate);
k=0;
- for(i=1; i<= (nlstate+ndeath); i++){
+ for(i=1; i<= (nlstate); i++){
for(j=1; j<=(nlstate+ndeath);j++){
k=k+1;
gp[k]=pmmij[i][j];
@@ -2009,92 +2082,169 @@ fprintf(ficresprob,"#One-step probabilit
pmij(pmmij,cov,ncovmodel,xp,nlstate);
k=0;
- for(i=1; i<=(nlstate+ndeath); i++){
+ for(i=1; i<=(nlstate); i++){
for(j=1; j<=(nlstate+ndeath);j++){
k=k+1;
gm[k]=pmmij[i][j];
}
}
- for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
+ for(i=1; i<= (nlstate)*(nlstate+ndeath); i++)
gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];
}
- for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
+ for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
for(theta=1; theta <=npar; theta++)
trgradg[j][theta]=gradg[theta][j];
- matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
- matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
+ matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
+ matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
pmij(pmmij,cov,ncovmodel,x,nlstate);
k=0;
- for(i=1; i<=(nlstate+ndeath); i++){
+ for(i=1; i<=(nlstate); i++){
for(j=1; j<=(nlstate+ndeath);j++){
k=k+1;
- gm[k]=pmmij[i][j];
+ mu[k][(int) age]=pmmij[i][j];
}
}
-
- /*printf("\n%d ",(int)age);
- for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
+ for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
+ for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
+ varpij[i][j][(int)age] = doldm[i][j];
+
+ /*printf("\n%d ",(int)age);
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
}*/
fprintf(ficresprob,"\n%d ",(int)age);
+ fprintf(ficresprobcov,"\n%d ",(int)age);
+ fprintf(ficresprobcor,"\n%d ",(int)age);
- for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
- fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
-
- }
- }
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
+ fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
+ for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
+ fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
+ fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
+ }
+ i=0;
+ for (k=1; k<=(nlstate);k++){
+ for (l=1; l<=(nlstate+ndeath);l++){
+ i=i++;
+ fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
+ fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
+ for (j=1; j<=i;j++){
+ fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
+ fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
+ }
+ }
+ }/* end of loop for state */
+ } /* end of loop for age */
+ /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
+ for (k1=1; k1<=(nlstate);k1++){
+ for (l1=1; l1<=(nlstate+ndeath);l1++){
+ if(l1==k1) continue;
+ i=(k1-1)*(nlstate+ndeath)+l1;
+ for (k2=1; k2<=(nlstate);k2++){
+ for (l2=1; l2<=(nlstate+ndeath);l2++){
+ if(l2==k2) continue;
+ j=(k2-1)*(nlstate+ndeath)+l2;
+ if(j<=i) continue;
+ for (age=bage; age<=fage; age ++){
+ if ((int)age %5==0){
+ v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
+ v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM;
+ cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
+ mu1=mu[i][(int) age]/stepm*YEARM ;
+ mu2=mu[j][(int) age]/stepm*YEARM;
+ /* Computing eigen value of matrix of covariance */
+ lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
+ lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
+ printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);
+ /* Eigen vectors */
+ v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
+ v21=sqrt(1.-v11*v11);
+ v12=-v21;
+ v22=v11;
+ /*printf(fignu*/
+ /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
+ /* mu2+ v21*lc1*cost + v21*lc2*sin(t) */
+ if(first==1){
+ first=0;
+ fprintf(ficgp,"\nset parametric;set nolabel");
+ fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1);
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+ fprintf(fichtm,"\n
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1 :varpijgr%s%d%1d%1d-%1d%1d.png, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1);
+ fprintf(fichtm,"\n
",optionfilefiname, j1,k2,l2,k1,l1);
+ fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
+ fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\
+ mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
+ }else{
+ first=0;
+ fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
+ fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
+ fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\
+ mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
+ }/* if first */
+ } /* age mod 5 */
+ } /* end loop age */
+ fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1);
+ first=1;
+ } /*l12 */
+ } /* k12 */
+ } /*l1 */
+ }/* k1 */
+ } /* loop covariates */
+ free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+ free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
}
free_vector(xp,1,npar);
fclose(ficresprob);
-
+ fclose(ficresprobcov);
+ fclose(ficresprobcor);
+ fclose(ficgp);
+ fclose(fichtm);
}
+
/******************* 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[],char rfileres[], char optionfilegnuplot[],\
- char version[], int popforecast, int estepm ){
+ int lastpass, int stepm, int weightopt, char model[],\
+ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
+ int popforecast, int estepm ,\
+ double jprev1, double mprev1,double anprev1, \
+ double jprev2, double mprev2,double anprev2){
int jj1, k1, i1, cpt;
- FILE *fichtm;
/*char optionfilehtm[FILENAMELENGTH];*/
-
- strcpy(optionfilehtm,optionfile);
- strcat(optionfilehtm,".htm");
- if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
+ if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
printf("Problem with %s \n",optionfilehtm), exit(0);
}
- fprintf(fichtm," %s
\n
-Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n
-\n
-Total number of observations=%d
\n
-Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n
-
- - Outputs files
\n
- - Copy of the parameter file: o%s
\n
- - Gnuplot file name: %s
\n
- - Observed prevalence in each state: p%s
\n
- - Stationary prevalence in each state: pl%s
\n
- - Transition probabilities: pij%s
\n
- - Life expectancies by age and initial health status (estepm=%2d months): e%s
\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
-
- fprintf(fichtm,"\n
- - Parameter file with estimated parameters and the covariance matrix: %s
\n
- - Variance of one-step probabilities: prob%s
\n
- - Variances of life expectancies by age and initial health status (estepm=%d months): v%s
\n
- - Health expectancies with their variances: t%s
\n
- - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
+ fprintf(fichtm,"Result files (first order: no variance)
\n
+ - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): p%s
\n
+ - Estimated transition probabilities over %d (stepm) months: pij%s
\n
+ - Stable prevalence in each health state: pl%s
\n
+ - Life expectancies by age and initial health status (estepm=%2d months):
+ e%s
\n ", \
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
+
+ fprintf(fichtm,"\n
Result files (second order: variances)
\n
+ - Parameter file with estimated parameters and covariance matrix: %s
\n
+ - Variance of one-step probabilities: prob%s
\n
+ - Variance-covariance of one-step probabilities: probcov%s
\n
+ - Correlation matrix of one-step probabilities: probcor%s
\n
+ - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): v%s
\n
+ - Health expectancies with their variances (no covariance): t%s
\n
+ - Standard deviation of stable prevalences: vpl%s
\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
if(popforecast==1) fprintf(fichtm,"\n
- Prevalences forecasting: f%s
\n
@@ -2102,7 +2252,7 @@ Interval (in months) between two waves:
",fileres,fileres,fileres,fileres);
else
fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)
\n",popforecast, stepm, model);
-fprintf(fichtm," - Graphs
");
+fprintf(fichtm,"
- Graphs
");
m=cptcoveff;
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
@@ -2110,17 +2260,22 @@ fprintf(fichtm,"
- Graphs
");
jj1=0;
for(k1=1; k1<=m;k1++){
for(i1=1; i1<=ncodemax[k1];i1++){
- jj1++;
- if (cptcovn > 0) {
- fprintf(fichtm,"
************ 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
");
- }
- fprintf(fichtm,"
- Probabilities: pe%s%d.png
-",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
+ jj1++;
+ if (cptcovn > 0) {
+ fprintf(fichtm,"
************ 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
");
+ }
+ /* Pij */
+ fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png
+",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
+ /* Quasi-incidences */
+ fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png
+",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
+ /* Stable prevalence in each health state */
for(cpt=1; cpt- Prevalence of disability : p%s%d%d.png
+ fprintf(fichtm,"
- Stable prevalence in each health state : p%s%d%d.png
",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
}
for(cpt=1; cpt<=nlstate;cpt++) {
@@ -2135,20 +2290,17 @@ interval) in state (%d): v%s%d%d.png
- Total life expectancy by age and
health expectancies in states (1) and (2): e%s%d.png
",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
-fprintf(fichtm,"\n");
- }
}
+ }
fclose(fichtm);
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplot(char fileres[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
-
- strcpy(optionfilegnuplot,optionfilefiname);
- strcat(optionfilegnuplot,".gp.txt");
- if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+ int ng;
+ if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
printf("Problem with file %s",optionfilegnuplot);
}
@@ -2162,11 +2314,11 @@ m=pow(2,cptcoveff);
for (k1=1; k1<= m ; k1 ++) {
#ifdef windows
- fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+ fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
#endif
#ifdef unix
-fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
#endif
@@ -2193,7 +2345,7 @@ fprintf(ficgp,"\nset ter png small\nset
/*2 eme*/
for (k1=1; k1<= m ; k1 ++) {
- fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1);
+ fprintf(ficgp,"\nset out \"e%s%d.png\" \n",strtok(optionfile, "."),k1);
fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
for (i=1; i<= nlstate+1 ; i ++) {
@@ -2226,7 +2378,7 @@ fprintf(ficgp,"\nset ter png small\nset
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<= nlstate ; cpt ++) {
k=2+nlstate*(2*cpt-2);
- fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
+ fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
@@ -2247,7 +2399,7 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt %s
\n
+Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n
+\n
+Total number of observations=%d
\n
+Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n
+
+ Parameter files
\n
+ - Copy of the parameter file: o%s
\n
+ - Gnuplot file name: %s
\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);
+ fclose(fichtm);
+
+ printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
/*------------ free_vector -------------*/
chdir(path);
@@ -3264,11 +3440,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
fclose(ficparo);
fclose(ficres);
-/*--------- index.htm --------*/
- printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm);
-
-
/*--------------- Prevalence limit --------------*/
strcpy(filerespl,"pl");
@@ -3329,7 +3501,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){
agelim=AGESUP;
hstepm=stepsize*YEARM; /* Every year of age */
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
-
+
+ /* hstepm=1; aff par mois*/
+
k=0;
for(cptcov=1;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
@@ -3342,6 +3516,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+
+ /* nhstepm=nhstepm*YEARM; aff par mois*/
+
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);
@@ -3351,7 +3528,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
fprintf(ficrespij," %1d-%1d",i,j);
fprintf(ficrespij,"\n");
for (h=0; h<=nhstepm; h++){
- fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
+ fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)
fprintf(ficrespij," %.5f", p3mat[i][j][h]);
@@ -3363,7 +3540,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
}
}
- varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
+ varprob(optionfilefiname, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
fclose(ficrespij);
@@ -3519,6 +3696,11 @@ free_matrix(mint,1,maxwav,1,n);
free_matrix(agev,1,maxwav,1,imx);
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
+ fprintf(fichtm,"\n