--- imach096d/src/imach.c 2002/05/24 16:34:18 1.45 +++ imach096d/src/imach.c 2002/06/10 13:12:49 1.48 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.45 2002/05/24 16:34:18 lievre Exp $ +/* $Id: imach.c,v 1.48 2002/06/10 13:12:49 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -75,8 +75,13 @@ #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; @@ -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++){ + 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); @@ -1928,39 +1937,86 @@ 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, i1, k1, l1; + 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+ndeath)*(nlstate+ndeath),1,npar); - doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); - + 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,"\n
We 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"); + + } cov[1]=1; j=cptcoveff; if (cptcovn<1) {j=1;ncodemax[1]=1;} @@ -1971,24 +2027,35 @@ fprintf(ficresprob,"#One-step probabilit if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); + fprintf(ficresprobcov, "\n#********** Variable "); + fprintf(ficgp, "\n#********** Variable "); + fprintf(fichtm, "\n

********** Variable

\n "); + fprintf(ficresprobcor, "\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 (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcov, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(fichtm, "**********\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 +2064,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,85 +2076,153 @@ 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%1d%1d-%1d%1d.png, ",k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1); + fprintf(fichtm,"\n
, ",optionfilefiname,k2,l2,k1,l1); + fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\"",optionfilefiname,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%1d%1d-%1d%1d.png\";replot;",optionfilefiname,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 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 -
- \n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot); - fprintf(fichtm,"