--- imach/src/imach.c 2002/03/10 13:43:02 1.31
+++ imach/src/imach.c 2002/03/13 17:19:16 1.34
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.31 2002/03/10 13:43:02 brouard Exp $
+/* $Id: imach.c,v 1.34 2002/03/13 17:19:16 brouard Exp $
Interpolated Markov Chain
Short summary of the programme:
@@ -82,7 +82,7 @@ int cptcovn, cptcovage=0, cptcoveff=0,cp
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
-int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
int popbased=0;
int *wav; /* Number of waves for this individuual 0 is possible */
@@ -1513,7 +1513,7 @@ void tricode(int *Tvar, int **nbcode, in
ij=1;
for (i=1; i<=10; i++) {
- if((Ndum[i]!=0) && (i<=ncov)){
+ if((Ndum[i]!=0) && (i<=ncovcol)){
Tvaraff[ij]=i;
ij++;
}
@@ -1540,14 +1540,14 @@ void evsij(char fileres[], double ***eij
k=1; /* For example stepm=6 months */
hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */
- hstepm=1; /* or (b) We decided to compute the life expectancy with the smallest unit */
+ hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */
/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
nhstepm is the number of hstepm from age to agelim
nstepm is the number of stepm from age to agelin.
Look at hpijx to understand the reason of that which relies in memory size
and note for a fixed period like k years */
/* We decided (b) to get a life expectancy respecting the most precise curvature of the
- survival function given par stepm (the optimization length). Unfortunately it
+ survival function given by stepm (the optimization length). Unfortunately it
means that if the survival funtion is printed only each two years of age and if
you sum them up and add 1 year (area under the trapezoids) you won't get the same
results. So we changed our mind and took the option of the best precision.
@@ -1559,13 +1559,13 @@ void evsij(char fileres[], double ***eij
/* nhstepm age range expressed in number of stepm */
nstepm=(int) rint((agelim-age)*YEARM/stepm);
/* Typically if 20 years nstepm = 20*12/6=40 stepm */
- if (stepm >= YEARM) hstepm=1;
+ /* if (stepm >= YEARM) hstepm=1;*/
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
/* Computed by stepm unit matrices, product of hstepm matrices, stored
in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);
- hf=hstepm/YEARM; /* Duration of hstepm expressed in year unit. */
+ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
@@ -1903,7 +1903,7 @@ void printinghtml(char fileres[], char t
printf("Problem with %s \n",optionfilehtm), exit(0);
}
- fprintf(fichtm,"
Imach, Version 0.71a
+ fprintf(fichtm," Imach, Version 0.8
Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
Total number of observations=%d
@@ -2466,7 +2466,7 @@ int main(int argc, char *argv[])
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
- char version[80]="Imach version 0.71a, March 2002, INED-EUROREVES ";
+ char version[80]="Imach version 0.8, March 2002, INED-EUROREVES ";
char *alph[]={"a","a","b","c","d","e"}, str[4];
@@ -2527,9 +2527,9 @@ int main(int argc, char *argv[])
}
ungetc(c,ficpar);
- fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
- printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
- fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
+ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+ fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
@@ -2686,7 +2686,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
- for (j=ncov;j>=1;j--){
+ for (j=ncovcol;j>=1;j--){
cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
}
num[i]=atol(stra);
@@ -2755,7 +2755,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
}
else {
cutv(strb,stre,strc,'V');
- Tvar[i]=ncov+k1;
+ Tvar[i]=ncovcol+k1;
cutv(strb,strc,strd,'V');
Tprod[k1]=i;
Tvard[k1][1]=atoi(strc);
@@ -2763,7 +2763,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
Tvar[cptcovn+k2]=Tvard[k1][1];
Tvar[cptcovn+k2+1]=Tvard[k1][2];
for (k=1; k<=lastobs;k++)
- covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
+ covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
k1++;
k2=k2+2;
}
@@ -2917,12 +2917,12 @@ printf("Total number of individuals= %d,
}
/*--------- results files --------------*/
- fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, weightopt,model);
+ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
jk=1;
- fprintf(ficres,"# Parameters\n");
- printf("# Parameters\n");
+ fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
+ printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i)
@@ -2944,8 +2944,8 @@ printf("Total number of individuals= %d,
ftolhess=ftol; /* Usually correct */
hesscov(matcov, p, npar, delti, ftolhess, func);
}
- fprintf(ficres,"# Scales\n");
- printf("# Scales\n");
+ fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
+ printf("# Scales (for hessian or gradient estimation)\n");
for(i=1,jk=1; i <=nlstate; i++){
for(j=1; j <=nlstate+ndeath; j++){
if (j!=i) {
@@ -2963,8 +2963,8 @@ printf("Total number of individuals= %d,
}
k=1;
- fprintf(ficres,"# Covariance\n");
- printf("# Covariance\n");
+ fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
+ printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
for(i=1;i<=npar;i++){
/* if (k>nlstate) k=1;
i1=(i-1)/(ncovmodel*nlstate)+1;
@@ -3175,15 +3175,15 @@ while((c=getc(ficpar))=='#' && c!= EOF){
/*---------- Forecasting ------------------*/
- if((stepm == 1) && (model==".")){
+ if((stepm == 1) && (strcmp(model,".")==0)){
prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
-if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
+ if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
free_matrix(mint,1,maxwav,1,n);
free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
free_vector(weight,1,n);}
else{
erreur=108;
- printf("Error %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);
+ printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
}
@@ -3243,8 +3243,6 @@ if (popforecast==1) populforecast(filere
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
fprintf(ficrest,"\n");
- hf=1;
- if (stepm >= YEARM) hf=stepm/YEARM;
epj=vector(1,nlstate+1);
for(age=bage; age <=fage ;age++){
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
@@ -3253,19 +3251,19 @@ if (popforecast==1) populforecast(filere
prlim[i][i]=probs[(int)age][i][k];
}
- fprintf(ficrest," %.0f",age);
+ fprintf(ficrest," %4.0f",age);
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*hf*eij[i][j][(int)age];
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
}
epj[nlstate+1] +=epj[j];
}
for(i=1, vepp=0.;i <=nlstate;i++)
for(j=1;j <=nlstate;j++)
vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %.2f (%.2f)", epj[nlstate+1],hf*sqrt(vepp));
+ fprintf(ficrest," %7.2f (%7.2f)", epj[nlstate+1],sqrt(vepp));
for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %.2f (%.2f)", epj[j],hf*sqrt(vareij[j][j][(int)age]));
+ fprintf(ficrest," %7.2f (%7.2f)", epj[j],sqrt(vareij[j][j][(int)age]));
}
fprintf(ficrest,"\n");
}
@@ -3323,7 +3321,7 @@ if (popforecast==1) populforecast(filere
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
if(erreur >0)
- printf("End of Imach with error %d\n",erreur);
+ printf("End of Imach with error or warning %d\n",erreur);
else printf("End of Imach\n");
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */