version 1.86, 2003/06/17 20:04:08
|
version 1.87, 2003/06/18 12:26:01
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.87 2003/06/18 12:26:01 brouard |
|
Version 0.96 |
|
|
Revision 1.86 2003/06/17 20:04:08 brouard |
Revision 1.86 2003/06/17 20:04:08 brouard |
(Module): Change position of html and gnuplot routines and added |
(Module): Change position of html and gnuplot routines and added |
routine fileappend. |
routine fileappend. |
Line 161
|
Line 164
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
|
char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES "; |
char version[]="Imach version 0.96, June 2003, INED-EUROREVES "; |
char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
int erreur; /* Error number */ |
int erreur; /* Error number */ |
int nvar; |
int nvar; |
Line 175 int popbased=0;
|
Line 178 int popbased=0;
|
int *wav; /* Number of waves for this individuual 0 is possible */ |
int *wav; /* Number of waves for this individuual 0 is possible */ |
int maxwav; /* Maxim number of waves */ |
int maxwav; /* Maxim number of waves */ |
int jmin, jmax; /* min, max spacing between 2 waves */ |
int jmin, jmax; /* min, max spacing between 2 waves */ |
|
int gipmx, gsw; /* Global variables on the number of contributions |
|
to the likelihood and the sum of weights (done by funcone)*/ |
int mle, weightopt; |
int mle, weightopt; |
int **mw; /* mw[mi][i] is number of the mi wave for this individual */ |
int **mw; /* mw[mi][i] is number of the mi wave for this individual */ |
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ |
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ |
Line 1290 double func( double *x)
|
Line 1295 double func( double *x)
|
/*************** log-likelihood *************/ |
/*************** log-likelihood *************/ |
double funcone( double *x) |
double funcone( double *x) |
{ |
{ |
|
/* Same as likeli but slower because of a lot of printf and if */ |
int i, ii, j, k, mi, d, kk; |
int i, ii, j, k, mi, d, kk; |
double l, ll[NLSTATEMAX], cov[NCOVMAX]; |
double l, ll[NLSTATEMAX], cov[NCOVMAX]; |
double **out; |
double **out; |
double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
|
double llt; |
int s1, s2; |
int s1, s2; |
double bbh, survp; |
double bbh, survp; |
/*extern weight */ |
/*extern weight */ |
Line 1354 double funcone( double *x)
|
Line 1361 double funcone( double *x)
|
%10.6f %10.6f %10.6f ", \ |
%10.6f %10.6f %10.6f ", \ |
num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], |
num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], |
2*weight[i]*lli,out[s1][s2],savm[s1][s2]); |
2*weight[i]*lli,out[s1][s2],savm[s1][s2]); |
for(k=1,l=0.; k<=nlstate; k++) |
for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
fprintf(ficresilk," %10.6f",ll[k]); |
llt +=ll[k]*gipmx/gsw; |
fprintf(ficresilk,"\n"); |
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); |
|
} |
|
fprintf(ficresilk," %10.6f\n", -llt); |
} |
} |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
|
if(globpr==0){ /* First time we count the contributions and weights */ |
|
gipmx=ipmx; |
|
gsw=sw; |
|
} |
return -l; |
return -l; |
} |
} |
|
|
|
|
void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpr, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) |
void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) |
{ |
{ |
/* This routine should help understanding what is done with the selection of individuals/waves and |
/* This routine should help understanding what is done with |
|
the selection of individuals/waves and |
to check the exact contribution to the likelihood. |
to check the exact contribution to the likelihood. |
Plotting could be done. |
Plotting could be done. |
*/ |
*/ |
int k; |
int k; |
if(globpr !=0){ /* Just counts and sums no printings */ |
|
|
if(*globpri !=0){ /* Just counts and sums no printings */ |
strcpy(fileresilk,"ilk"); |
strcpy(fileresilk,"ilk"); |
strcat(fileresilk,fileres); |
strcat(fileresilk,fileres); |
if((ficresilk=fopen(fileresilk,"w"))==NULL) { |
if((ficresilk=fopen(fileresilk,"w"))==NULL) { |
printf("Problem with resultfile: %s\n", fileresilk); |
printf("Problem with resultfile: %s\n", fileresilk); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); |
} |
} |
fprintf(ficresilk, "#individual(line's record) s1 s2 wave# effective_wave# number_of_product_matrix pij weight 2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state\n"); |
fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); |
fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight out sav "); |
fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight out sav "); |
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ |
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ |
for(k=1; k<=nlstate; k++) |
for(k=1; k<=nlstate; k++) |
fprintf(ficresilk," ll[%d]",k); |
fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); |
fprintf(ficresilk,"\n"); |
fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); |
} |
} |
|
|
*fretone=(*funcone)(p); |
*fretone=(*funcone)(p); |
if(globpr !=0){ |
if(*globpri !=0){ |
fclose(ficresilk); |
fclose(ficresilk); |
if((fichtm=fopen(optionfilehtm,"a"))==NULL) { |
fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk,fileresilk); |
printf("Problem with html file: %s\n", optionfilehtm); |
fflush(fichtm); |
fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); |
} |
exit(0); |
|
} |
|
else{ |
|
fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",fileresilk); |
|
fclose(fichtm); |
|
} |
|
} |
|
return; |
return; |
} |
} |
|
|
Line 2368 void varevsij(char optionfilefiname[], d
|
Line 2376 void varevsij(char optionfilefiname[], d
|
else{ |
else{ |
fprintf(ficgp,"\n# Routine varevsij"); |
fprintf(ficgp,"\n# Routine varevsij"); |
} |
} |
if((fichtm=fopen(optionfilehtm,"a"))==NULL) { |
/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ |
printf("Problem with html file: %s\n", optionfilehtm); |
/* printf("Problem with html file: %s\n", optionfilehtm); */ |
fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); |
/* fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */ |
exit(0); |
/* exit(0); */ |
} |
/* } */ |
else{ |
/* else{ */ |
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
} |
/* } */ |
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
|
|
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,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); |
Line 2606 void varevsij(char optionfilefiname[], d
|
Line 2614 void varevsij(char optionfilefiname[], d
|
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
fclose(ficresprobmorprev); |
fclose(ficresprobmorprev); |
fclose(ficgp); |
fclose(ficgp); |
fclose(fichtm); |
/* fclose(fichtm); */ |
} /* end varevsij */ |
} /* end varevsij */ |
|
|
/************ Variance of prevlim ******************/ |
/************ Variance of prevlim ******************/ |
Line 2771 void varprob(char optionfilefiname[], do
|
Line 2779 void varprob(char optionfilefiname[], do
|
else{ |
else{ |
fprintf(ficgp,"\n# Routine varprob"); |
fprintf(ficgp,"\n# Routine varprob"); |
} |
} |
if((fichtm=fopen(optionfilehtm,"a"))==NULL) { |
/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ |
printf("Problem with html file: %s\n", optionfilehtm); |
/* printf("Problem with html file: %s\n", optionfilehtm); */ |
fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); |
/* fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm); */ |
exit(0); |
/* exit(0); */ |
} |
/* } */ |
else{ |
/* else{ */ |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
|
|
Line 2784 void varprob(char optionfilefiname[], do
|
Line 2792 void varprob(char optionfilefiname[], do
|
fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 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. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n"); |
fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 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. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n"); |
|
|
} |
/* } */ |
|
|
cov[1]=1; |
cov[1]=1; |
tj=cptcoveff; |
tj=cptcoveff; |
Line 3001 void varprob(char optionfilefiname[], do
|
Line 3009 void varprob(char optionfilefiname[], do
|
fclose(ficresprobcov); |
fclose(ficresprobcov); |
fclose(ficresprobcor); |
fclose(ficresprobcor); |
fclose(ficgp); |
fclose(ficgp); |
fclose(fichtm); |
|
} |
} |
|
|
|
|
Line 3014 void printinghtml(char fileres[], char t
|
Line 3021 void printinghtml(char fileres[], char t
|
double jprev2, double mprev2,double anprev2){ |
double jprev2, double mprev2,double anprev2){ |
int jj1, k1, i1, cpt; |
int jj1, k1, i1, cpt; |
/*char optionfilehtm[FILENAMELENGTH];*/ |
/*char optionfilehtm[FILENAMELENGTH];*/ |
if((fichtm=fopen(optionfilehtm,"a"))==NULL) { |
/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ |
printf("Problem with %s \n",optionfilehtm), exit(0); |
/* printf("Problem with %s \n",optionfilehtm), exit(0); */ |
fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); |
/* fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */ |
} |
/* } */ |
|
|
fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \ |
fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \ |
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \ |
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n \ |
Line 3104 interval) in state (%d): v%s%d%d.png <br
|
Line 3111 interval) in state (%d): v%s%d%d.png <br
|
} /* end i1 */ |
} /* end i1 */ |
}/* End k1 */ |
}/* End k1 */ |
fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
fclose(fichtm); |
fflush(fichtm); |
} |
} |
|
|
/******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
Line 3611 populforecast(char fileres[], double anp
|
Line 3618 populforecast(char fileres[], double anp
|
fclose(ficrespop); |
fclose(ficrespop); |
} /* End of popforecast */ |
} /* End of popforecast */ |
|
|
int fileappend(FILE *fichier, char *optionfile) |
int fileappend(FILE *fichier, char *optionfich) |
{ |
{ |
if((fichier=fopen(optionfile,"a"))==NULL) { |
if((fichier=fopen(optionfich,"a"))==NULL) { |
printf("Problem with file: %s\n", optionfile); |
printf("Problem with file: %s\n", optionfich); |
fprintf(ficlog,"Problem with file: %s\n", optionfile); |
fprintf(ficlog,"Problem with file: %s\n", optionfich); |
return (1); |
return (0); |
} |
} |
|
fflush(fichier); |
|
return (1); |
} |
} |
/***********************************************/ |
/***********************************************/ |
/**************** Main Program *****************/ |
/**************** Main Program *****************/ |
Line 3630 int main(int argc, char *argv[])
|
Line 3638 int main(int argc, char *argv[])
|
int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; |
int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; |
int jj; |
int jj; |
int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
|
/* FILE *fichtm; *//* Html File */ |
|
/* FILE *ficgp;*/ /*Gnuplot File */ |
double agedeb, agefin,hf; |
double agedeb, agefin,hf; |
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
|
|
Line 4268 int main(int argc, char *argv[])
|
Line 4278 int main(int argc, char *argv[])
|
fclose(ficgp); |
fclose(ficgp); |
/*--------- index.htm --------*/ |
/*--------- index.htm --------*/ |
|
|
strcpy(optionfilehtm,optionfile); |
strcpy(optionfilehtm,optionfilefiname); |
strcat(optionfilehtm,".htm"); |
strcat(optionfilehtm,".htm"); |
if((fichtm=fopen(optionfilehtm,"w"))==NULL) { |
if((fichtm=fopen(optionfilehtm,"w"))==NULL) { |
printf("Problem with %s \n",optionfilehtm), exit(0); |
printf("Problem with %s \n",optionfilehtm), exit(0); |
} |
} |
|
|
fprintf(fichtm,"<body> <font size=\"2\">%s <br> %s</font> \ |
fprintf(fichtm,"<body>\n<title>IMaCh %s</title>\n <font size=\"2\">%s <br> %s</font> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\ |
\n\ |
\n\ |
Line 4282 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4292 Title=%s <br>Datafile=%s Firstpass=%d La
|
<ul><li><h4>Parameter files</h4>\n\ |
<ul><li><h4>Parameter files</h4>\n\ |
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\ |
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\ |
- Log file of the run: <a href=\"%s\">%s</a><br>\n\ |
- Log file of the run: <a href=\"%s\">%s</a><br>\n\ |
- Gnuplot file name: <a href=\"%s\">%s</a>\n\ |
- Gnuplot file name: <a href=\"%s\">%s</a><br>\n\ |
- Date and time at start: %s</ul>\n",\ |
- Date and time at start: %s</ul>\n",\ |
version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\ |
fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\ |
model,fileres,fileres,\ |
model,fileres,fileres,\ |
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strt); |
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strt); |
fclose(fichtm); |
/*fclose(fichtm);*/ |
|
fflush(fichtm); |
|
|
/* Calculates basic frequencies. Computes observed prevalence at single age |
/* Calculates basic frequencies. Computes observed prevalence at single age |
and prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); |
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); |
|
|
if(fileappend(fichtm, optionfilehtm)){ |
/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ |
|
/* printf("Problem with file: %s\n", optionfilehtm); */ |
|
/* fprintf(ficlog,"Problem with file: %s\n", optionfilehtm); */ |
|
/* } */ |
|
|
|
|
|
/* if(fileappend(fichtm, optionfilehtm)){ */ |
|
fprintf(fichtm,"\n"); |
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
imx,agemin,agemax,jmin,jmax,jmean); |
imx,agemin,agemax,jmin,jmax,jmean); |
fclose(fichtm); |
/* fclose(fichtm); */ |
} |
/* } */ |
|
|
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
Line 4837 ageminpar, agemax, s[lastpass][imx], age
|
Line 4856 ageminpar, agemax, s[lastpass][imx], age
|
printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); |
printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); |
fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); |
fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); |
/* printf("Total time was %d uSec.\n", total_usecs);*/ |
/* printf("Total time was %d uSec.\n", total_usecs);*/ |
if(fileappend(fichtm,optionfilehtm)){ |
/* if(fileappend(fichtm,optionfilehtm)){ */ |
fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend); |
fprintf(fichtm,"<br>Localtime at start %s and at end=%s<br>",strt, strtend); |
fclose(fichtm); |
fclose(fichtm); |
} |
|
/*------ End -----------*/ |
/*------ End -----------*/ |
|
|
end: |
end: |