version 1.79, 2003/06/05 15:17:23
|
version 1.84, 2003/06/13 21:44:43
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
|
$State$ |
|
$Log$ |
|
Revision 1.84 2003/06/13 21:44:43 brouard |
|
* imach.c (Repository): Replace "freqsummary" at a correct |
|
place. It differs from routine "prevalence" which may be called |
|
many times. Probs is memory consuming and must be used with |
|
parcimony. |
|
Version 0.95a2 (should output exactly the same maximization than 0.8a2) |
|
|
|
Revision 1.83 2003/06/10 13:39:11 lievre |
|
*** empty log message *** |
|
|
|
Revision 1.82 2003/06/05 15:57:20 brouard |
|
Add log in imach.c and fullversion number is now printed. |
|
|
|
*/ |
|
/* |
Interpolated Markov Chain |
Interpolated Markov Chain |
|
|
Short summary of the programme: |
Short summary of the programme: |
Line 58
|
Line 75
|
read parameterfile |
read parameterfile |
read datafile |
read datafile |
concatwav |
concatwav |
|
freqsummary |
if (mle >= 1) |
if (mle >= 1) |
mlikeli |
mlikeli |
print results files |
print results files |
Line 120
|
Line 138
|
#endif |
#endif |
|
|
/* $Id$ */ |
/* $Id$ */ |
/* $Log$ |
/* $State$ */ |
* Revision 1.79 2003/06/05 15:17:23 brouard |
|
* *** empty log message *** |
char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES "; |
* */ |
char fullversion[]="$Revision$ $Date$"; |
char version[80]="Imach version 0.95a1, June 2003, INED-EUROREVES "; |
|
int erreur; /* Error number */ |
int erreur; /* Error number */ |
int nvar; |
int nvar; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
Line 1168 double func( double *x)
|
Line 1185 double func( double *x)
|
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
}else{ /* ml=4 no inter-extrapolation */ |
}else if (mle==4){ /* ml=4 no inter-extrapolation */ |
|
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
|
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
|
for(mi=1; mi<= wav[i]-1; mi++){ |
|
for (ii=1;ii<=nlstate+ndeath;ii++) |
|
for (j=1;j<=nlstate+ndeath;j++){ |
|
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
|
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
|
} |
|
for(d=0; d<dh[mi][i]; d++){ |
|
newm=savm; |
|
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
for (kk=1; kk<=cptcovage;kk++) { |
|
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
|
} |
|
|
|
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
|
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
|
savm=oldm; |
|
oldm=newm; |
|
} /* end mult */ |
|
|
|
s1=s[mw[mi][i]][i]; |
|
s2=s[mw[mi+1][i]][i]; |
|
if( s2 > nlstate){ |
|
lli=log(out[s1][s2] - savm[s1][s2]); |
|
}else{ |
|
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
|
} |
|
ipmx +=1; |
|
sw += weight[i]; |
|
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
|
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ |
|
} /* end of wave */ |
|
} /* end of individual */ |
|
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
Line 1190 double func( double *x)
|
Line 1242 double func( double *x)
|
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
|
s1=s[mw[mi][i]][i]; |
|
s2=s[mw[mi+1][i]][i]; |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
|
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
} /* End of if */ |
} /* End of if */ |
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 */ |
|
/*exit(0); */ |
return -l; |
return -l; |
} |
} |
|
|
Line 1495 void lubksb(double **a, int n, int *indx
|
Line 1551 void lubksb(double **a, int n, int *indx
|
} |
} |
|
|
/************ Frequencies ********************/ |
/************ Frequencies ********************/ |
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) |
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint) |
{ /* Some frequencies */ |
{ /* Some frequencies */ |
|
|
int i, m, jk, k1,i1, j1, bool, z1,z2,j; |
int i, m, jk, k1,i1, j1, bool, z1,z2,j; |
Line 1549 void freqsummary(char fileres[], int ia
|
Line 1605 void freqsummary(char fileres[], int ia
|
if (bool==1){ |
if (bool==1){ |
for(m=firstpass; m<=lastpass; m++){ |
for(m=firstpass; m<=lastpass; m++){ |
k2=anint[m][i]+(mint[m][i]/12.); |
k2=anint[m][i]+(mint[m][i]/12.); |
if ((k2>=dateprev1) && (k2<=dateprev2)) { |
/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ |
if(agev[m][i]==0) agev[m][i]=iagemax+1; |
if(agev[m][i]==0) agev[m][i]=iagemax+1; |
if(agev[m][i]==1) agev[m][i]=iagemax+2; |
if(agev[m][i]==1) agev[m][i]=iagemax+2; |
if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; |
if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; |
Line 1562 void freqsummary(char fileres[], int ia
|
Line 1618 void freqsummary(char fileres[], int ia
|
dateintsum=dateintsum+k2; |
dateintsum=dateintsum+k2; |
k2cpt++; |
k2cpt++; |
} |
} |
} |
/*}*/ |
} |
} |
} |
} |
} |
} |
|
|
fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
|
|
if (cptcovn>0) { |
if (cptcovn>0) { |
fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
Line 1628 void freqsummary(char fileres[], int ia
|
Line 1684 void freqsummary(char fileres[], int ia
|
if( i <= iagemax){ |
if( i <= iagemax){ |
if(pos>=1.e-5){ |
if(pos>=1.e-5){ |
fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); |
probs[i][jk][j1]= pp[jk]/pos; |
/*probs[i][jk][j1]= pp[jk]/pos;*/ |
/*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ |
/*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ |
} |
} |
else |
else |
Line 1661 void freqsummary(char fileres[], int ia
|
Line 1717 void freqsummary(char fileres[], int ia
|
} |
} |
|
|
/************ Prevalence ********************/ |
/************ Prevalence ********************/ |
void prevalence(double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
{ |
{ |
/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
in each health status at the date of interview (if between dateprev1 and dateprev2). |
in each health status at the date of interview (if between dateprev1 and dateprev2). |
Line 2395 void varevsij(char optionfilefiname[], d
|
Line 2451 void varevsij(char optionfilefiname[], d
|
fclose(ficresprobmorprev); |
fclose(ficresprobmorprev); |
fclose(ficgp); |
fclose(ficgp); |
fclose(fichtm); |
fclose(fichtm); |
} |
} /* end varevsij */ |
|
|
/************ Variance of prevlim ******************/ |
/************ Variance of prevlim ******************/ |
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) |
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) |
Line 2992 m=pow(2,cptcoveff);
|
Line 3048 m=pow(2,cptcoveff);
|
fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1); |
fprintf(ficgp,"\nset out \"p%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] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1); |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1); |
|
|
for (i=1; i<= nlstate ; i ++) |
for (i=1; i< nlstate ; i ++) |
fprintf(ficgp,"+$%d",k+i+1); |
fprintf(ficgp,"+$%d",k+i+1); |
fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1); |
fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1); |
|
|
Line 3131 prevforecast(char fileres[], double anpr
|
Line 3187 prevforecast(char fileres[], double anpr
|
char fileresf[FILENAMELENGTH]; |
char fileresf[FILENAMELENGTH]; |
|
|
agelim=AGESUP; |
agelim=AGESUP; |
prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
|
|
strcpy(fileresf,"f"); |
strcpy(fileresf,"f"); |
strcat(fileresf,fileres); |
strcat(fileresf,fileres); |
Line 3254 populforecast(char fileres[], double anp
|
Line 3310 populforecast(char fileres[], double anp
|
agelim=AGESUP; |
agelim=AGESUP; |
calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; |
calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; |
|
|
prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
|
|
|
|
strcpy(filerespop,"pop"); |
strcpy(filerespop,"pop"); |
Line 3396 populforecast(char fileres[], double anp
|
Line 3452 populforecast(char fileres[], double anp
|
free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
fclose(ficrespop); |
fclose(ficrespop); |
} |
} /* End of popforecast */ |
|
|
/***********************************************/ |
/***********************************************/ |
/**************** Main Program *****************/ |
/**************** Main Program *****************/ |
Line 3459 int main(int argc, char *argv[])
|
Line 3515 int main(int argc, char *argv[])
|
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
getcwd(pathcd, size); |
getcwd(pathcd, size); |
|
|
printf("\n%s",version); |
printf("\n%s\n%s",version,fullversion); |
if(argc <=1){ |
if(argc <=1){ |
printf("\nEnter the parameter file name: "); |
printf("\nEnter the parameter file name: "); |
scanf("%s",pathtot); |
scanf("%s",pathtot); |
Line 3837 int main(int argc, char *argv[])
|
Line 3893 int main(int argc, char *argv[])
|
if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ |
if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ |
printf("Error! Month of death of individual %d on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); |
printf("Error! Month of death of individual %d on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); |
fprintf(ficlog,"Error! Month of death of individual %d on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); |
fprintf(ficlog,"Error! Month of death of individual %d on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); |
s[m][i]=-1; |
s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ |
} |
} |
} |
} |
} |
} |
Line 3959 int main(int argc, char *argv[])
|
Line 4015 int main(int argc, char *argv[])
|
|
|
/* 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); |
|
|
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 */ |
Line 4002 int main(int argc, char *argv[])
|
Line 4059 int main(int argc, char *argv[])
|
} |
} |
} |
} |
} |
} |
if(mle==1){ |
if(mle!=0){ |
/* Computing hessian and covariance matrix */ |
/* Computing hessian and covariance matrix */ |
ftolhess=ftol; /* Usually correct */ |
ftolhess=ftol; /* Usually correct */ |
hesscov(matcov, p, npar, delti, ftolhess, func); |
hesscov(matcov, p, npar, delti, ftolhess, func); |
Line 4133 int main(int argc, char *argv[])
|
Line 4190 int main(int argc, char *argv[])
|
fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
|
|
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
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,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); |
/*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
|
|
/*------------ gnuplot -------------*/ |
/*------------ gnuplot -------------*/ |
strcpy(optionfilegnuplot,optionfilefiname); |
strcpy(optionfilegnuplot,optionfilefiname); |
Line 4298 Interval (in months) between two waves:
|
Line 4355 Interval (in months) between two waves:
|
|
|
fclose(ficrespij); |
fclose(ficrespij); |
|
|
|
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
|
|
/*---------- Forecasting ------------------*/ |
/*---------- Forecasting ------------------*/ |
/*if((stepm == 1) && (strcmp(model,".")==0)){*/ |
/*if((stepm == 1) && (strcmp(model,".")==0)){*/ |
Line 4345 Interval (in months) between two waves:
|
Line 4403 Interval (in months) between two waves:
|
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
|
|
/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ |
/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ |
prevalence(agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
/* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ |
/* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ |
ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
*/ |
*/ |