#include <unistd.h>\r
\r
#define MAXLINE 256\r
-#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"\r
+#define GNUPLOTPROGRAM "wgnuplot"\r
+/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/\r
#define FILENAMELENGTH 80\r
/*#define DEBUG*/\r
#define windows\r
\r
for (k=1; k<=cptcovn;k++) {\r
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];\r
- /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/\r
+ /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/\r
}\r
- for (k=1; k<=cptcovage;k++)\r
- cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];\r
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];\r
for (k=1; k<=cptcovprod;k++)\r
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];\r
\r
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/\r
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/\r
-\r
+ /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/\r
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);\r
\r
savm=oldm;\r
/************ Frequencies ********************/\r
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)\r
{ /* Some frequencies */\r
- \r
+ \r
int i, m, jk, k1,i1, j1, bool, z1,z2,j;\r
double ***freq; /* Frequencies */\r
double *pp;\r
double pos, k2, dateintsum=0,k2cpt=0;\r
FILE *ficresp;\r
char fileresp[FILENAMELENGTH];\r
- \r
+ \r
pp=vector(1,nlstate);\r
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
strcpy(fileresp,"p");\r
}\r
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);\r
j1=0;\r
-\r
+ \r
j=cptcoveff;\r
if (cptcovn<1) {j=1;ncodemax[1]=1;}\r
-\r
+ \r
for(k1=1; k1<=j;k1++){\r
- for(i1=1; i1<=ncodemax[k1];i1++){\r
- j1++;\r
- /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);\r
- scanf("%d", i);*/\r
- for (i=-1; i<=nlstate+ndeath; i++) \r
- for (jk=-1; jk<=nlstate+ndeath; jk++) \r
- for(m=agemin; m <= agemax+3; m++)\r
- freq[i][jk][m]=0;\r
-\r
- dateintsum=0;\r
- k2cpt=0;\r
- for (i=1; i<=imx; i++) {\r
- bool=1;\r
- if (cptcovn>0) {\r
- for (z1=1; z1<=cptcoveff; z1++) \r
- if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) \r
- bool=0;\r
- }\r
- if (bool==1) {\r
- for(m=firstpass; m<=lastpass; m++){\r
- k2=anint[m][i]+(mint[m][i]/12.);\r
- if ((k2>=dateprev1) && (k2<=dateprev2)) {\r
- if(agev[m][i]==0) agev[m][i]=agemax+1;\r
- if(agev[m][i]==1) agev[m][i]=agemax+2;\r
- freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
- freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];\r
- if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {\r
- dateintsum=dateintsum+k2;\r
- k2cpt++;\r
- }\r
-\r
- }\r
- }\r
- }\r
- }\r
+ for(i1=1; i1<=ncodemax[k1];i1++){\r
+ j1++;\r
+ /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);\r
+ scanf("%d", i);*/\r
+ for (i=-1; i<=nlstate+ndeath; i++) \r
+ for (jk=-1; jk<=nlstate+ndeath; jk++) \r
+ for(m=agemin; m <= agemax+3; m++)\r
+ freq[i][jk][m]=0;\r
+ \r
+ dateintsum=0;\r
+ k2cpt=0;\r
+ for (i=1; i<=imx; i++) {\r
+ bool=1;\r
+ if (cptcovn>0) {\r
+ for (z1=1; z1<=cptcoveff; z1++) \r
+ if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) \r
+ bool=0;\r
+ }\r
+ if (bool==1) {\r
+ for(m=firstpass; m<=lastpass; m++){\r
+ k2=anint[m][i]+(mint[m][i]/12.);\r
+ if ((k2>=dateprev1) && (k2<=dateprev2)) {\r
+ if(agev[m][i]==0) agev[m][i]=agemax+1;\r
+ if(agev[m][i]==1) agev[m][i]=agemax+2;\r
+ if (m<lastpass) {\r
+ freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
+ freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];\r
+ }\r
+ \r
+ if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {\r
+ dateintsum=dateintsum+k2;\r
+ k2cpt++;\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
\r
- fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
+ fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
\r
- if (cptcovn>0) {\r
- fprintf(ficresp, "\n#********** Variable "); \r
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
- fprintf(ficresp, "**********\n#");\r
+ if (cptcovn>0) {\r
+ fprintf(ficresp, "\n#********** Variable "); \r
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);\r
+ fprintf(ficresp, "**********\n#");\r
+ }\r
+ for(i=1; i<=nlstate;i++) \r
+ fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);\r
+ fprintf(ficresp, "\n");\r
+ \r
+ for(i=(int)agemin; i <= (int)agemax+3; i++){\r
+ if(i==(int)agemax+3)\r
+ printf("Total");\r
+ else\r
+ printf("Age %d", i);\r
+ for(jk=1; jk <=nlstate ; jk++){\r
+ for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)\r
+ pp[jk] += freq[jk][m][i]; \r
+ }\r
+ for(jk=1; jk <=nlstate ; jk++){\r
+ for(m=-1, pos=0; m <=0 ; m++)\r
+ pos += freq[jk][m][i];\r
+ if(pp[jk]>=1.e-10)\r
+ printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);\r
+ else\r
+ printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);\r
}\r
- for(i=1; i<=nlstate;i++) \r
- fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);\r
- fprintf(ficresp, "\n");\r
- \r
- for(i=(int)agemin; i <= (int)agemax+3; i++){\r
- if(i==(int)agemax+3)\r
- printf("Total");\r
- else\r
- printf("Age %d", i);\r
- for(jk=1; jk <=nlstate ; jk++){\r
- for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)\r
- pp[jk] += freq[jk][m][i]; \r
- }\r
- for(jk=1; jk <=nlstate ; jk++){\r
- for(m=-1, pos=0; m <=0 ; m++)\r
- pos += freq[jk][m][i];\r
- if(pp[jk]>=1.e-10)\r
- printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);\r
- else\r
- printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);\r
- }\r
\r
- for(jk=1; jk <=nlstate ; jk++){\r
- for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)\r
- pp[jk] += freq[jk][m][i];\r
- }\r
+ for(jk=1; jk <=nlstate ; jk++){\r
+ for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)\r
+ pp[jk] += freq[jk][m][i];\r
+ }\r
\r
- for(jk=1,pos=0; jk <=nlstate ; jk++)\r
- pos += pp[jk];\r
- for(jk=1; jk <=nlstate ; jk++){\r
- if(pos>=1.e-5)\r
- printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);\r
- else\r
- printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);\r
- if( i <= (int) agemax){\r
- if(pos>=1.e-5){\r
- fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);\r
- probs[i][jk][j1]= pp[jk]/pos;\r
- /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/\r
+ for(jk=1,pos=0; jk <=nlstate ; jk++)\r
+ pos += pp[jk];\r
+ for(jk=1; jk <=nlstate ; jk++){\r
+ if(pos>=1.e-5)\r
+ printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);\r
+ else\r
+ printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);\r
+ if( i <= (int) agemax){\r
+ if(pos>=1.e-5){\r
+ fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);\r
+ probs[i][jk][j1]= pp[jk]/pos;\r
+ /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/\r
+ }\r
+ else\r
+ fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);\r
+ }\r
}\r
- else\r
- fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);\r
+ \r
+ for(jk=-1; jk <=nlstate+ndeath; jk++)\r
+ for(m=-1; m <=nlstate+ndeath; m++)\r
+ if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);\r
+ if(i <= (int) agemax)\r
+ fprintf(ficresp,"\n");\r
+ printf("\n");\r
}\r
}\r
- for(jk=-1; jk <=nlstate+ndeath; jk++)\r
- for(m=-1; m <=nlstate+ndeath; m++)\r
- if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);\r
- if(i <= (int) agemax)\r
- fprintf(ficresp,"\n");\r
- printf("\n");\r
- }\r
- }\r
- }\r
+ }\r
dateintmean=dateintsum/k2cpt; \r
\r
fclose(ficresp);\r
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);\r
free_vector(pp,1,nlstate);\r
-\r
+ \r
/* End of Freq */\r
}\r
\r
if ((k2>=dateprev1) && (k2<=dateprev2)) {\r
if(agev[m][i]==0) agev[m][i]=agemax+1;\r
if(agev[m][i]==1) agev[m][i]=agemax+2;\r
- freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];\r
+ if (m<lastpass) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];\r
/* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; */ \r
}\r
}\r
for (k=0; k<=19; k++) {\r
if (Ndum[k] != 0) {\r
nbcode[Tvar[j]][ij]=k; \r
+ /* printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/\r
ij++;\r
}\r
if (ij > ncodemax[j]) break; \r
{\r
/* Health expectancies */\r
int i, j, nhstepm, hstepm, h, nstepm, k;\r
- double age, agelim,hf;\r
+ double age, agelim, hf;\r
double ***p3mat;\r
\r
fprintf(ficreseij,"# Health expectancies\n");\r
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/\r
double **newm;\r
double **dnewm,**doldm;\r
- int i, j, nhstepm, hstepm, h;\r
+ int i, j, nhstepm, hstepm, h, nstepm, kk;\r
int k, cptcode;\r
double *xp;\r
double **gp, **gm;\r
double ***gradg, ***trgradg;\r
double ***p3mat;\r
- double age,agelim;\r
+ double age,agelim, hf;\r
int theta;\r
\r
fprintf(ficresvij,"# Covariances of life expectancies\n");\r
dnewm=matrix(1,nlstate,1,npar);\r
doldm=matrix(1,nlstate,1,nlstate);\r
\r
- hstepm=1*YEARM; /* Every year of age */\r
- hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ \r
+ kk=1; /* For example stepm=6 months */\r
+ hstepm=kk*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */\r
+ hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */\r
+ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. \r
+ nhstepm is the number of hstepm from age to agelim \r
+ nstepm is the number of stepm from age to agelin. \r
+ Look at hpijx to understand the reason of that which relies in memory size\r
+ and note for a fixed period like k years */\r
+ /* We decided (b) to get a life expectancy respecting the most precise curvature of the\r
+ survival function given by stepm (the optimization length). Unfortunately it\r
+ means that if the survival funtion is printed only each two years of age and if\r
+ you sum them up and add 1 year (area under the trapezoids) you won't get the same \r
+ results. So we changed our mind and took the option of the best precision.\r
+ */\r
+ hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ \r
agelim = AGESUP;\r
for (age=bage; age<=fage; age ++){ /* If stepm=6 months */\r
- nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ \r
- if (stepm >= YEARM) hstepm=1;\r
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */\r
+ nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ \r
+ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */\r
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
gradg=ma3x(0,nhstepm,1,npar,1,nlstate);\r
gp=matrix(0,nhstepm,1,nlstate);\r
for(theta=1; theta <=npar; theta++)\r
trgradg[h][j][theta]=gradg[h][theta][j];\r
\r
+ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */\r
for(i=1;i<=nlstate;i++)\r
for(j=1;j<=nlstate;j++)\r
vareij[i][j][(int)age] =0.;\r
+\r
for(h=0;h<=nhstepm;h++){\r
for(k=0;k<=nhstepm;k++){\r
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);\r
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);\r
for(i=1;i<=nlstate;i++)\r
for(j=1;j<=nlstate;j++)\r
- vareij[i][j][(int)age] += doldm[i][j];\r
+ vareij[i][j][(int)age] += doldm[i][j]*hf*hf;\r
}\r
}\r
- h=1;\r
- if (stepm >= YEARM) h=stepm/YEARM;\r
+\r
fprintf(ficresvij,"%.0f ",age );\r
for(i=1; i<=nlstate;i++)\r
for(j=1; j<=nlstate;j++){\r
- fprintf(ficresvij," %.4f", h*vareij[i][j][(int)age]);\r
+ fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);\r
}\r
fprintf(ficresvij,"\n");\r
free_matrix(gp,0,nhstepm,1,nlstate);\r
}\r
\r
/******************* Printing html file ***********/\r
-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[] ){\r
+void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \\r
+ int lastpass, int stepm, int weightopt, char model[],\\r
+ int imx,int jmin, int jmax, double jmeanint,char optionfile[], \\r
+ char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\\r
+ char version[], int popforecast ){\r
int jj1, k1, i1, cpt;\r
FILE *fichtm;\r
/*char optionfilehtm[FILENAMELENGTH];*/\r
printf("Problem with %s \n",optionfilehtm), exit(0);\r
}\r
\r
- fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.8 </font> <hr size=\"2\" color=\"#EC5E5E\"> \r
-Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\r
-\r
-Total number of observations=%d <br>\r
-Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\r
+ fprintf(fichtm,"<body> <font size=\"2\">Imach, Version %s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n\r
+Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\r
+\n\r
+Total number of observations=%d <br>\n\r
+Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n\r
<hr size=\"2\" color=\"#EC5E5E\"> \r
-<li>Outputs files<br><br>\n\r
- - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n\r
-- Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\r
- - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>\r
- - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>\r
- - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\r
- - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>\r
- - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>\r
- - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\r
- - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\r
- - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\r
- - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\r
- <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);\r
- \r
+ <ul><li>Outputs files<br>\n\r
+ - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\r
+ - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\r
+ - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n\r
+ - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>\n\r
+ - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>\n\r
+ - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);\r
+\r
+ fprintf(fichtm,"\n\r
+ - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n\r
+ - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>\n\r
+ - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n\r
+ - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres);\r
+\r
+ if(popforecast==1) fprintf(fichtm,"\n\r
+ - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n\r
+ - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n\r
+ <br>",fileres,fileres,fileres,fileres);\r
+ else \r
+ fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);\r
fprintf(fichtm," <li>Graphs</li><p>");\r
\r
m=cptcoveff;\r
}\r
\r
/******************* Gnuplot file **************/\r
-void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemaxpar, double fage , char pathc[], double p[]){\r
+void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){\r
\r
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;\r
\r
strcpy(optionfilegnuplot,optionfilefiname);\r
- strcat(optionfilegnuplot,".plt");\r
+ strcat(optionfilegnuplot,".gp.txt");\r
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {\r
printf("Problem with file %s",optionfilegnuplot);\r
}\r
for (k1=1; k1<= m ; k1 ++) {\r
\r
#ifdef windows\r
- 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);\r
+ 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",ageminpar,fage,fileres,k1-1,k1-1);\r
#endif\r
#ifdef unix\r
-fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);\r
+fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);\r
#endif\r
\r
for (i=1; i<= nlstate ; i ++) {\r
/*2 eme*/\r
\r
for (k1=1; k1<= m ; k1 ++) { \r
- fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);\r
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);\r
\r
for (i=1; i<= nlstate+1 ; i ++) {\r
k=2*i;\r
for (k1=1; k1<= m ; k1 ++) { \r
for (cpt=1; cpt<= nlstate ; cpt ++) {\r
k=2+nlstate*(cpt-1);\r
- 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);\r
+ 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",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);\r
for (i=1; i< nlstate ; i ++) {\r
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);\r
} \r
for (k1=1; k1<= m ; k1 ++) { \r
for (cpt=1; cpt<nlstate ; cpt ++) {\r
k=3;\r
- 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,agemaxpar,fileres,k1,k+cpt+1,k+1);\r
+ 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",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);\r
\r
for (i=1; i< nlstate ; i ++)\r
fprintf(ficgp,"+$%d",k+i+1);\r
}\r
\r
for(jk=1; jk <=m; jk++) {\r
- fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemaxpar);\r
+ fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);\r
i=1;\r
for(k2=1; k2<=nlstate; k2++) {\r
k3=i;\r
\r
\r
/*************** Moving average **************/\r
-void movingaverage(double agedeb, double fage,double agemin, double ***mobaverage){\r
+void movingaverage(double agedeb, double fage,double ageminpar, double ***mobaverage){\r
\r
int i, cpt, cptcod;\r
- for (agedeb=agemin; agedeb<=fage; agedeb++)\r
+ for (agedeb=ageminpar; agedeb<=fage; agedeb++)\r
for (i=1; i<=nlstate;i++)\r
for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)\r
mobaverage[(int)agedeb][i][cptcod]=0.;\r
\r
- for (agedeb=agemin+4; agedeb<=fage; agedeb++){\r
+ for (agedeb=ageminpar+4; agedeb<=fage; agedeb++){\r
for (i=1; i<=nlstate;i++){\r
for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
for (cpt=0;cpt<=4;cpt++){\r
\r
\r
/************** Forecasting ******************/\r
-prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){\r
+prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){\r
\r
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;\r
int *popage;\r
agelim=AGESUP;\r
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;\r
\r
- prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
+ prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
\r
\r
strcpy(fileresf,"f"); \r
\r
if (mobilav==1) {\r
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
- movingaverage(agedeb, fage, agemin, mobaverage);\r
+ movingaverage(agedeb, fage, ageminpar, mobaverage);\r
}\r
\r
stepsize=(int) (stepm+YEARM-1)/YEARM;\r
fprintf(ficresf,"\n");\r
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); \r
\r
- for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+ for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ \r
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
nhstepm = nhstepm/hstepm; \r
\r
\r
for (h=0; h<=nhstepm; h++){\r
if (h==(int) (calagedate+YEARM*cpt)) {\r
- fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);\r
+ fprintf(ficresf,"\n %.f %.f ",anproj1+cpt,agedeb+h*hstepm/YEARM*stepm);\r
} \r
for(j=1; j<=nlstate+ndeath;j++) {\r
kk1=0.;kk2=0;\r
fclose(ficresf);\r
}\r
/************** Forecasting ******************/\r
-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){\r
+populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){\r
\r
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;\r
int *popage;\r
agelim=AGESUP;\r
calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;\r
\r
- prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
+ prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
\r
\r
strcpy(filerespop,"pop"); \r
\r
if (mobilav==1) {\r
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
- movingaverage(agedeb, fage, agemin, mobaverage);\r
+ movingaverage(agedeb, fage, ageminpar, mobaverage);\r
}\r
\r
stepsize=(int) (stepm+YEARM-1)/YEARM;\r
for (cpt=0; cpt<=0;cpt++) { \r
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); \r
\r
- for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+ for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ \r
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
nhstepm = nhstepm/hstepm; \r
\r
\r
for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { \r
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); \r
- for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+ for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ \r
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
nhstepm = nhstepm/hstepm; \r
\r
\r
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;\r
double agedeb, agefin,hf;\r
- double agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;\r
+ double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;\r
\r
double fret;\r
double **xi,tmp,delta;\r
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;\r
\r
\r
- char version[80]="Imach version 0.8, March 2002, INED-EUROREVES ";\r
+ char version[80]="Imach version 0.8a, March 2002, INED-EUROREVES ";\r
char *alph[]={"a","a","b","c","d","e"}, str[4];\r
\r
\r
struct timeval start_time, end_time;\r
\r
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */\r
-\r
+ getcwd(pathcd, size);\r
\r
printf("\n%s",version);\r
if(argc <=1){\r
if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;\r
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;\r
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;\r
- }\r
-\r
- for (i=1; i<=imx; i++)\r
- if (covar[1][i]==0) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));*/\r
-\r
+ }*/\r
+ \r
+ /* for (i=1; i<=imx; i++){\r
+ if (s[4][i]==9) s[4][i]=-1; \r
+ printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}\r
+ */\r
+ \r
/* Calculation of the number of parameter from char model*/\r
Tvar=ivector(1,15); \r
Tprod=ivector(1,15); \r
cptcovn=j+1;\r
cptcovprod=j1;\r
\r
- \r
strcpy(modelsav,model); \r
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){\r
printf("Error. Non available option model=%s ",model);\r
}\r
}\r
\r
- /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);\r
+ /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);\r
printf("cptcovprod=%d ", cptcovprod);\r
scanf("%d ",i);*/\r
fclose(fic);\r
/*-calculation of age at interview from date of interview and age at death -*/\r
agev=matrix(1,maxwav,1,imx);\r
\r
- for (i=1; i<=imx; i++) \r
- for(m=2; (m<= maxwav); m++)\r
+ for (i=1; i<=imx; i++) {\r
+ for(m=2; (m<= maxwav); m++) {\r
if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){\r
anint[m][i]=9999;\r
s[m][i]=-1;\r
}\r
- \r
+ if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;\r
+ }\r
+ }\r
+\r
for (i=1; i<=imx; i++) {\r
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);\r
for(m=1; (m<= maxwav); m++){\r
if(s[m][i] >0){\r
- if (s[m][i] == nlstate+1) {\r
+ if (s[m][i] >= nlstate+1) {\r
if(agedc[i]>0)\r
if(moisdc[i]!=99 && andc[i]!=9999)\r
- agev[m][i]=agedc[i];\r
- else {\r
+ agev[m][i]=agedc[i];\r
+ /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/\r
+ else {\r
if (andc[i]!=9999){\r
printf("Warning negative age at death: %d line:%d\n",num[i],i);\r
agev[m][i]=-1;\r
for(j=1; j <= ncodemax[k]; j++){\r
for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){\r
h++;\r
- if (h>m) h=1;codtab[h][k]=j;\r
+ if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;\r
+ /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/\r
} \r
}\r
}\r
} \r
-\r
-\r
- /*for(i=1; i <=m ;i++){ \r
- for(k=1; k <=cptcovn; k++){\r
- printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);\r
- }\r
- printf("\n");\r
- }\r
- scanf("%d",i);*/\r
+ /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); \r
+ codtab[1][2]=1;codtab[2][2]=2; */\r
+ /* for(i=1; i <=m ;i++){ \r
+ for(k=1; k <=cptcovn; k++){\r
+ printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);\r
+ }\r
+ printf("\n");\r
+ }\r
+ scanf("%d",i);*/\r
\r
/* Calculates basic frequencies. Computes observed prevalence at single age\r
and prints on file fileres'p'. */\r
}\r
ungetc(c,ficpar);\r
\r
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemaxpar, &bage, &fage);\r
+ fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&ageminpar,&agemaxpar, &bage, &fage);\r
\r
if (fage <= 2) {\r
- bage = agemin;\r
+ bage = ageminpar;\r
fage = agemaxpar;\r
}\r
\r
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");\r
- fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);\r
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);\r
+ fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage);\r
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage);\r
\r
while((c=getc(ficpar))=='#' && c!= EOF){\r
ungetc(c,ficpar);\r
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
\r
/*------------ gnuplot -------------*/\r
- printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemaxpar,fage, pathc,p);\r
+ printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,fage, pathc,p);\r
\r
/*------------ free_vector -------------*/\r
chdir(path);\r
\r
/*--------- index.htm --------*/\r
\r
- printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres);\r
+ printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast);\r
\r
\r
/*--------------- Prevalence limit --------------*/\r
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */\r
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */\r
k=0;\r
- agebase=agemin;\r
+ agebase=ageminpar;\r
agelim=agemaxpar;\r
ftolpl=1.e-10;\r
i1=cptcoveff;\r
\r
fprintf(ficreseij,"\n#****** ");\r
for(j=1;j<=cptcoveff;j++) \r
- fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);\r
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
fprintf(ficreseij,"******\n");\r
\r
fprintf(ficresvij,"\n#****** ");\r
for(j=1;j<=cptcoveff;j++) \r
- fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);\r
+ fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
fprintf(ficresvij,"******\n");\r
\r
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);\r
\r
#ifdef windows\r
while (z[0] != 'q') {\r
- chdir(path); \r
- printf("\nType e to edit output files, c to start again, and q for exiting: ");\r
+ /* chdir(path); */\r
+ printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");\r
scanf("%s",z);\r
if (z[0] == 'c') system("./imach");\r
- else if (z[0] == 'e') {\r
- chdir(path);\r
- system(optionfilehtm);\r
- }\r
+ else if (z[0] == 'e') system(optionfilehtm);\r
+ else if (z[0] == 'g') system(plotcmd);\r
else if (z[0] == 'q') exit(0);\r
}\r
#endif \r
git://henry.ined.fr / .git / commitdiff