#define AGEBASE 40\r
#ifdef windows\r
#define DIRSEPARATOR '\\'\r
+#define ODIRSEPARATOR '/'\r
#else\r
#define DIRSEPARATOR '/'\r
+#define ODIRSEPARATOR '\\'\r
#endif\r
\r
-char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES ";\r
+char version[80]="Imach version 0.8i, June 2002, INED-EUROREVES ";\r
int erreur; /* Error number */\r
int nvar;\r
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;\r
double **oldm, **newm, **savm; /* Working pointers to matrices */\r
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */\r
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;\r
+FILE *ficlog;\r
FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;\r
+FILE *ficresprobmorprev;\r
FILE *fichtm; /* Html File */\r
FILE *ficreseij;\r
char filerese[FILENAMELENGTH];\r
char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];\r
\r
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];\r
-\r
+char filelog[FILENAMELENGTH]; /* Log file */\r
char filerest[FILENAMELENGTH];\r
char fileregp[FILENAMELENGTH];\r
char popfile[FILENAMELENGTH];\r
\r
l1 = strlen( path ); /* length of path */\r
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );\r
- s = strrchr( path, DIRSEPARATOR ); /* find last / */\r
+ s= strrchr( path, DIRSEPARATOR ); /* find last / */\r
if ( s == NULL ) { /* no directory, so use current */\r
+ /*if(strrchr(path, ODIRSEPARATOR )==NULL)\r
+ printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/\r
#if defined(__bsd__) /* get current working directory */\r
extern char *getwd( );\r
\r
\r
void cutv(char *u,char *v, char*t, char occ)\r
{\r
+ /* cuts string t into u and v where u is ended by char occ excluding it\r
+ and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)\r
+ gives u="abcedf" and v="ghi2j" */\r
int i,lg,j,p=0;\r
i=0;\r
for(j=0; j<=strlen(t)-1; j++) {\r
tol2=2.0*(tol1=tol*fabs(x)+ZEPS); \r
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/\r
printf(".");fflush(stdout);\r
+ fprintf(ficlog,".");fflush(ficlog);\r
#ifdef DEBUG\r
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);\r
+ fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);\r
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */\r
#endif\r
if (fabs(x-xm) <= (tol2-0.5*(b-a))){ \r
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); \r
#ifdef DEBUG\r
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);\r
+ fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);\r
#endif\r
for (j=1;j<=n;j++) { \r
xi[j] *= xmin; \r
ibig=0; \r
del=0.0; \r
printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);\r
+ fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);\r
for (i=1;i<=n;i++) \r
printf(" %d %.12f",i, p[i]);\r
+ fprintf(ficlog," %d %.12f",i, p[i]);\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
for (i=1;i<=n;i++) { \r
for (j=1;j<=n;j++) xit[j]=xi[j][i]; \r
fptt=(*fret); \r
#ifdef DEBUG\r
printf("fret=%lf \n",*fret);\r
+ fprintf(ficlog,"fret=%lf \n",*fret);\r
#endif\r
printf("%d",i);fflush(stdout);\r
+ fprintf(ficlog,"%d",i);fflush(ficlog);\r
linmin(p,xit,n,fret,func); \r
if (fabs(fptt-(*fret)) > del) { \r
del=fabs(fptt-(*fret)); \r
} \r
#ifdef DEBUG\r
printf("%d %.12e",i,(*fret));\r
+ fprintf(ficlog,"%d %.12e",i,(*fret));\r
for (j=1;j<=n;j++) {\r
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);\r
printf(" x(%d)=%.12e",j,xit[j]);\r
+ fprintf(ficlog," x(%d)=%.12e",j,xit[j]);\r
}\r
- for(j=1;j<=n;j++) \r
+ for(j=1;j<=n;j++) {\r
printf(" p=%.12e",p[j]);\r
+ fprintf(ficlog," p=%.12e",p[j]);\r
+ }\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
#endif\r
} \r
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {\r
k[0]=1;\r
k[1]=-1;\r
printf("Max: %.12e",(*func)(p));\r
- for (j=1;j<=n;j++) \r
+ fprintf(ficlog,"Max: %.12e",(*func)(p));\r
+ for (j=1;j<=n;j++) {\r
printf(" %.12e",p[j]);\r
+ fprintf(ficlog," %.12e",p[j]);\r
+ }\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
for(l=0;l<=1;l++) {\r
for (j=1;j<=n;j++) {\r
ptt[j]=p[j]+(p[j]-pt[j])*k[l];\r
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);\r
+ fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);\r
}\r
printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));\r
+ fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));\r
}\r
#endif\r
\r
}\r
#ifdef DEBUG\r
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);\r
- for(j=1;j<=n;j++)\r
+ fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);\r
+ for(j=1;j<=n;j++){\r
printf(" %.12e",xit[j]);\r
+ fprintf(ficlog," %.12e",xit[j]);\r
+ }\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
#endif\r
} \r
} \r
for (i=1;i<=npar;i++)\r
for (j=1;j<=npar;j++)\r
xi[i][j]=(i==j ? 1.0 : 0.0);\r
- printf("Powell\n");\r
+ printf("Powell\n"); fprintf(ficlog,"Powell\n");\r
powell(p,xi,npar,ftol,&iter,&fret,func);\r
\r
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));\r
+ fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));\r
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));\r
\r
}\r
hess=matrix(1,npar,1,npar);\r
\r
printf("\nCalculation of the hessian matrix. Wait...\n");\r
+ fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n");\r
for (i=1;i<=npar;i++){\r
printf("%d",i);fflush(stdout);\r
+ fprintf(ficlog,"%d",i);fflush(ficlog);\r
hess[i][i]=hessii(p,ftolhess,i,delti);\r
/*printf(" %f ",p[i]);*/\r
/*printf(" %lf ",hess[i][i]);*/\r
for (j=1;j<=npar;j++) {\r
if (j>i) { \r
printf(".%d%d",i,j);fflush(stdout);\r
+ fprintf(ficlog,".%d%d",i,j);fflush(ficlog);\r
hess[i][j]=hessij(p,delti,i,j);\r
hess[j][i]=hess[i][j]; \r
/*printf(" %lf ",hess[i][j]);*/\r
}\r
}\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
\r
printf("\nInverting the hessian to get the covariance matrix. Wait...\n");\r
+ fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");\r
\r
a=matrix(1,npar,1,npar);\r
y=matrix(1,npar,1,npar);\r
}\r
\r
printf("\n#Hessian matrix#\n");\r
+ fprintf(ficlog,"\n#Hessian matrix#\n");\r
for (i=1;i<=npar;i++) { \r
for (j=1;j<=npar;j++) { \r
printf("%.3e ",hess[i][j]);\r
+ fprintf(ficlog,"%.3e ",hess[i][j]);\r
}\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
}\r
\r
/* Recompute Inverse */\r
for (i=1;i<=npar;i++){ \r
y[i][j]=x[i];\r
printf("%.3e ",y[i][j]);\r
+ fprintf(ficlog,"%.3e ",y[i][j]);\r
}\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
}\r
*/\r
\r
\r
#ifdef DEBUG\r
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);\r
+ fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);\r
#endif\r
/*if(fabs(k1-2.0*fx+k2) <1.e-13){ */\r
if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)){\r
res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */\r
#ifdef DEBUG\r
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);\r
+ fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);\r
#endif\r
}\r
return res;\r
{ /* Some frequencies */\r
\r
int i, m, jk, k1,i1, j1, bool, z1,z2,j;\r
+ int first;\r
double ***freq; /* Frequencies */\r
double *pp;\r
double pos, k2, dateintsum=0,k2cpt=0;\r
strcat(fileresp,fileres);\r
if((ficresp=fopen(fileresp,"w"))==NULL) {\r
printf("Problem with prevalence resultfile: %s\n", fileresp);\r
+ fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);\r
exit(0);\r
}\r
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);\r
\r
j=cptcoveff;\r
if (cptcovn<1) {j=1;ncodemax[1]=1;}\r
- \r
+\r
+ first=1;\r
+\r
for(k1=1; k1<=j;k1++){\r
for(i1=1; i1<=ncodemax[k1];i1++){\r
j1++;\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
+ if(i==(int)agemax+3){\r
+ fprintf(ficlog,"Total");\r
+ }else{\r
+ if(first==1){\r
+ first=0;\r
+ printf("See log file for details...\n");\r
+ }\r
+ fprintf(ficlog,"Age %d", i);\r
+ }\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
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
+ if(pp[jk]>=1.e-10){\r
+ if(first==1){\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
+ fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);\r
+ }else{\r
+ if(first==1)\r
+ printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);\r
+ fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);\r
+ }\r
}\r
\r
for(jk=1; jk <=nlstate ; jk++){\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(pos>=1.e-5){\r
+ if(first==1)\r
+ printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);\r
+ fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);\r
+ }else{\r
+ if(first==1)\r
+ printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);\r
+ fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);\r
+ }\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
\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(freq[jk][m][i] !=0 ) {\r
+ if(first==1)\r
+ printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);\r
+ fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);\r
+ }\r
if(i <= (int) agemax)\r
fprintf(ficresp,"\n");\r
- printf("\n");\r
+ if(first==1)\r
+ printf("Others in log...\n");\r
+ fprintf(ficlog,"\n");\r
}\r
}\r
}\r
probs[i][jk][j1]= pp[jk]/pos;\r
}\r
}\r
- }\r
- \r
- }\r
- }\r
- }\r
+ }/* end jk */\r
+ }/* end i */\r
+ } /* end i1 */\r
+ } /* end k1 */\r
\r
\r
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);\r
int i, mi, m;\r
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;\r
double sum=0., jmean=0.;*/\r
-\r
+ int first;\r
int j, k=0,jk, ju, jl;\r
double sum=0.;\r
+ first=0;\r
jmin=1e+5;\r
jmax=-1;\r
jmean=0.;\r
}\r
\r
wav[i]=mi;\r
- if(mi==0)\r
- printf("Warning, no any valid information for:%d line=%d\n",num[i],i);\r
+ if(mi==0){\r
+ if(first==0){\r
+ printf("Warning, no any valid information for:%d line=%d and may be others, see log file\n",num[i],i);\r
+ first=1;\r
+ }\r
+ if(first==1){\r
+ fprintf(ficlog,"Warning, no any valid information for:%d line=%d\n",num[i],i);\r
+ }\r
+ } /* end mi==0 */\r
}\r
\r
for(i=1; i<=imx; i++){\r
}\r
jmean=sum/k;\r
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);\r
+ fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);\r
}\r
+\r
/*********** Tricode ****************************/\r
void tricode(int *Tvar, int **nbcode, int imx)\r
{\r
for (k=0; k<19; k++) Ndum[k]=0;\r
\r
for (i=1; i<=ncovmodel-2; i++) {\r
- ij=Tvar[i];\r
- Ndum[ij]++; \r
- }\r
+ ij=Tvar[i];\r
+ Ndum[ij]++; \r
+ }\r
\r
ij=1;\r
for (i=1; i<=10; i++) {\r
}\r
}\r
\r
- cptcoveff=ij-1;\r
+ cptcoveff=ij-1;\r
}\r
\r
/*********** Health Expectancies ****************/\r
varhe[i][j][(int)age] =0.;\r
\r
printf("%d|",(int)age);fflush(stdout);\r
+ fprintf(ficlog,"%d|",(int)age);fflush(ficlog);\r
for(h=0;h<=nhstepm-1;h++){\r
for(k=0;k<=nhstepm-1;k++){\r
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);\r
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
}\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
\r
free_vector(xp,1,npar);\r
free_matrix(dnewm,1,nlstate*2,1,npar);\r
}\r
\r
/************ Variance ******************/\r
-void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm)\r
+void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased)\r
{\r
/* Variance of health expectancies */\r
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/\r
- double **newm;\r
+ /* double **newm;*/\r
double **dnewm,**doldm;\r
+ double **dnewmp,**doldmp;\r
int i, j, nhstepm, hstepm, h, nstepm ;\r
int k, cptcode;\r
double *xp;\r
- double **gp, **gm;\r
- double ***gradg, ***trgradg;\r
+ double **gp, **gm; /* for var eij */\r
+ double ***gradg, ***trgradg; /*for var eij */\r
+ double **gradgp, **trgradgp; /* for var p point j */\r
+ double *gpp, *gmp; /* for var p point j */\r
+ double **varppt; /* for var p point j nlstate to nlstate+ndeath */\r
double ***p3mat;\r
double age,agelim, hf;\r
int theta;\r
+ char digit[4];\r
+ char digitp[16];\r
+\r
+ char fileresprobmorprev[FILENAMELENGTH];\r
+\r
+ if(popbased==1)\r
+ strcpy(digitp,"-populbased-");\r
+ else\r
+ strcpy(digitp,"-stablbased-");\r
+\r
+ strcpy(fileresprobmorprev,"prmorprev"); \r
+ sprintf(digit,"%-d",ij);\r
+ /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/\r
+ strcat(fileresprobmorprev,digit); /* Tvar to be done */\r
+ strcat(fileresprobmorprev,digitp); /* Popbased or not */\r
+ strcat(fileresprobmorprev,fileres);\r
+ if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {\r
+ printf("Problem with resultfile: %s\n", fileresprobmorprev);\r
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);\r
+ }\r
+ printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);\r
+ fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);\r
+ fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1*p1j+w2*p2j+... stand dev in()\n");\r
+ fprintf(ficresprobmorprev,"# Age cov=%-d",ij);\r
+ for(j=nlstate+1; j<=(nlstate+ndeath);j++){\r
+ fprintf(ficresprobmorprev," p.%-d SE",j);\r
+ for(i=1; i<=nlstate;i++)\r
+ fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);\r
+ } \r
+ fprintf(ficresprobmorprev,"\n");\r
+ if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {\r
+ printf("Problem with gnuplot file: %s\n", optionfilegnuplot);\r
+ fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);\r
+ exit(0);\r
+ }\r
+ else{\r
+ fprintf(ficgp,"\n# Routine varevsij");\r
+ }\r
+ if((fichtm=fopen(optionfilehtm,"a"))==NULL) {\r
+ printf("Problem with html file: %s\n", optionfilehtm);\r
+ fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);\r
+ exit(0);\r
+ }\r
+ else{\r
+ fprintf(fichtm,"\n<li><h4> Computing step probabilities of dying and weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");\r
+ }\r
+ varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);\r
\r
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");\r
fprintf(ficresvij,"# Age");\r
xp=vector(1,npar);\r
dnewm=matrix(1,nlstate,1,npar);\r
doldm=matrix(1,nlstate,1,nlstate);\r
+ dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);\r
+ doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);\r
+\r
+ gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);\r
+ gpp=vector(nlstate+1,nlstate+ndeath);\r
+ gmp=vector(nlstate+1,nlstate+ndeath);\r
+ trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/\r
\r
if(estepm < stepm){\r
printf ("Problem %d lower than %d\n",estepm, stepm);\r
gp=matrix(0,nhstepm,1,nlstate);\r
gm=matrix(0,nhstepm,1,nlstate);\r
\r
+\r
for(theta=1; theta <=npar; theta++){\r
for(i=1; i<=npar; i++){ /* Computes gradient */\r
xp[i] = x[i] + (i==theta ?delti[theta]:0);\r
gp[h][j] += prlim[i][i]*p3mat[i][j][h];\r
}\r
}\r
- \r
+ /* This for computing forces of mortality (h=1)as a weighted average */\r
+ for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){\r
+ for(i=1; i<= nlstate; i++)\r
+ gpp[j] += prlim[i][i]*p3mat[i][j][1];\r
+ } \r
+ /* end force of mortality */\r
+\r
for(i=1; i<=npar; i++) /* Computes gradient */\r
xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); \r
gm[h][j] += prlim[i][i]*p3mat[i][j][h];\r
}\r
}\r
-\r
- for(j=1; j<= nlstate; j++)\r
+ /* This for computing force of mortality (h=1)as a weighted average */\r
+ for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){\r
+ for(i=1; i<= nlstate; i++)\r
+ gmp[j] += prlim[i][i]*p3mat[i][j][1];\r
+ } \r
+ /* end force of mortality */\r
+\r
+ for(j=1; j<= nlstate; j++) /* vareij */\r
for(h=0; h<=nhstepm; h++){\r
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];\r
}\r
+ for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */\r
+ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];\r
+ }\r
+\r
} /* End theta */\r
\r
- trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);\r
+ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */\r
\r
- for(h=0; h<=nhstepm; h++)\r
+ for(h=0; h<=nhstepm; h++) /* veij */\r
for(j=1; j<=nlstate;j++)\r
for(theta=1; theta <=npar; theta++)\r
trgradg[h][j][theta]=gradg[h][theta][j];\r
\r
+ for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */\r
+ for(theta=1; theta <=npar; theta++)\r
+ trgradgp[j][theta]=gradgp[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
}\r
}\r
\r
+ /* pptj */\r
+ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);\r
+ matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);\r
+ for(j=nlstate+1;j<=nlstate+ndeath;j++)\r
+ for(i=nlstate+1;i<=nlstate+ndeath;i++)\r
+ varppt[j][i]=doldmp[j][i];\r
+ /* end ppptj */\r
+ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); \r
+ prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);\r
+ \r
+ if (popbased==1) {\r
+ for(i=1; i<=nlstate;i++)\r
+ prlim[i][i]=probs[(int)age][i][ij];\r
+ }\r
+ \r
+ /* This for computing force of mortality (h=1)as a weighted average */\r
+ for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){\r
+ for(i=1; i<= nlstate; i++)\r
+ gmp[j] += prlim[i][i]*p3mat[i][j][1];\r
+ } \r
+ /* end force of mortality */\r
+\r
+ fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);\r
+ for(j=nlstate+1; j<=(nlstate+ndeath);j++){\r
+ fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));\r
+ for(i=1; i<=nlstate;i++){\r
+ fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);\r
+ }\r
+ } \r
+ fprintf(ficresprobmorprev,"\n");\r
+\r
fprintf(ficresvij,"%.0f ",age );\r
for(i=1; i<=nlstate;i++)\r
for(j=1; j<=nlstate;j++){\r
free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);\r
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
} /* End age */\r
- \r
+ free_vector(gpp,nlstate+1,nlstate+ndeath);\r
+ free_vector(gmp,nlstate+1,nlstate+ndeath);\r
+ free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);\r
+ free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/\r
+ fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");\r
+ /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */\r
+ fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");\r
+ fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm);\r
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm);\r
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm);\r
+ fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);\r
+ fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);\r
+ fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);\r
+\r
free_vector(xp,1,npar);\r
- free_matrix(doldm,1,nlstate,1,npar);\r
- free_matrix(dnewm,1,nlstate,1,nlstate);\r
+ free_matrix(doldm,1,nlstate,1,nlstate);\r
+ free_matrix(dnewm,1,nlstate,1,npar);\r
+ free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);\r
+ free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);\r
+ free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);\r
+ fclose(ficresprobmorprev);\r
+ fclose(ficgp);\r
+ fclose(fichtm);\r
\r
}\r
\r
int i, j=0, i1, k1, l1, t, tj;\r
int k2, l2, j1, z1;\r
int k=0,l, cptcode;\r
- int first=1;\r
+ int first=1, first1;\r
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;\r
double **dnewm,**doldm;\r
double *xp;\r
strcat(fileresprob,fileres);\r
if((ficresprob=fopen(fileresprob,"w"))==NULL) {\r
printf("Problem with resultfile: %s\n", fileresprob);\r
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);\r
}\r
strcpy(fileresprobcov,"probcov"); \r
strcat(fileresprobcov,fileres);\r
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {\r
printf("Problem with resultfile: %s\n", fileresprobcov);\r
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);\r
}\r
strcpy(fileresprobcor,"probcor"); \r
strcat(fileresprobcor,fileres);\r
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {\r
printf("Problem with resultfile: %s\n", fileresprobcor);\r
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);\r
}\r
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);\r
+ fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);\r
printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);\r
+ fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);\r
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);\r
+ fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);\r
\r
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");\r
fprintf(ficresprob,"# Age");\r
first=1;\r
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {\r
printf("Problem with gnuplot file: %s\n", optionfilegnuplot);\r
+ fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);\r
exit(0);\r
}\r
else{\r
}\r
if((fichtm=fopen(optionfilehtm,"a"))==NULL) {\r
printf("Problem with html file: %s\n", optionfilehtm);\r
+ fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);\r
exit(0);\r
}\r
else{\r
+ fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");\r
+ fprintf(fichtm,"\n");\r
+\r
fprintf(fichtm,"\n<li><h4> Computing matrix of variance-covariance of step probabilities</h4></li>\n");\r
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");\r
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");\r
/*printf("\n%d ",(int)age);\r
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){\r
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));\r
+ fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));\r
}*/\r
\r
fprintf(ficresprob,"\n%d ",(int)age);\r
}\r
}/* end of loop for state */\r
} /* end of loop for age */\r
- /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/\r
+\r
+ /* Confidence intervalle of pij */\r
+ /*\r
+ fprintf(ficgp,"\nset noparametric;unset label");\r
+ fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");\r
+ fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");\r
+ fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);\r
+ fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);\r
+ fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);\r
+ fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);\r
+ */\r
+\r
+ /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/\r
+ first1=1;\r
for (k1=1; k1<=(nlstate);k1++){\r
for (l1=1; l1<=(nlstate+ndeath);l1++){ \r
if(l1==k1) continue;\r
/* Computing eigen value of matrix of covariance */\r
lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));\r
lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));\r
- printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);\r
+ if(first1==1){\r
+ first1=0;\r
+ printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\nOthers in log...\n",v1,v2,cv12,lc1,lc2);\r
+ }\r
+ fprintf(ficlog,"Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);\r
/* Eigen vectors */\r
v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));\r
v21=sqrt(1.-v11*v11);\r
/* mu2+ v21*lc1*cost + v21*lc2*sin(t) */\r
if(first==1){\r
first=0;\r
- fprintf(ficgp,"\nset parametric;set nolabel");\r
+ fprintf(ficgp,"\nset parametric;set nolabel");\r
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);\r
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");\r
fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1);\r
fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);\r
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);\r
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);\r
- 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\"",\\r
+ /* 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\"",\\r
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);\r
+ */\r
+ 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)) not",\\r
+ mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));\r
}else{\r
first=0;\r
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);\r
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);\r
+ /*\r
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\"",\\r
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);\r
+ */\r
+ 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)) not",\\r
+ mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
+ mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));\r
}/* if first */\r
} /* age mod 5 */\r
} /* end loop age */\r
/*char optionfilehtm[FILENAMELENGTH];*/\r
if((fichtm=fopen(optionfilehtm,"a"))==NULL) {\r
printf("Problem with %s \n",optionfilehtm), exit(0);\r
+ fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);\r
}\r
\r
fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n\r
<a href=\"e%s\">e%s</a> <br>\n</li>", \\r
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);\r
\r
+fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");\r
+\r
+ m=cptcoveff;\r
+ if (cptcovn < 1) {m=1;ncodemax[1]=1;}\r
+\r
+ jj1=0;\r
+ for(k1=1; k1<=m;k1++){\r
+ for(i1=1; i1<=ncodemax[k1];i1++){\r
+ jj1++;\r
+ if (cptcovn > 0) {\r
+ fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");\r
+ for (cpt=1; cpt<=cptcoveff;cpt++) \r
+ fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);\r
+ fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");\r
+ }\r
+ /* Pij */\r
+ fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>\r
+<img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); \r
+ /* Quasi-incidences */\r
+ fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>\r
+<img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); \r
+ /* Stable prevalence in each health state */\r
+ for(cpt=1; cpt<nlstate;cpt++){\r
+ fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>\r
+<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);\r
+ }\r
+ for(cpt=1; cpt<=nlstate;cpt++) {\r
+ fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>\r
+<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);\r
+ }\r
+ fprintf(fichtm,"\n<br>- Total life expectancy by age and\r
+health expectancies in states (1) and (2): e%s%d.png<br>\r
+<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);\r
+ } /* end i1 */\r
+ }/* End k1 */\r
+ fprintf(fichtm,"</ul>");\r
+\r
+\r
fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n\r
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n\r
- Variance of one-step probabilities: <a href=\"prob%s\">prob%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><b>Graphs</b></li><p>");\r
+fprintf(fichtm," <
ul><li><b>Graphs</b></li><p>");
\r \r
m=cptcoveff;\r
if (cptcovn < 1) {m=1;ncodemax[1]=1;}\r
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);\r
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");\r
}\r
- /* Pij */\r
- fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>\r
-<img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); \r
- /* Quasi-incidences */\r
- fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>\r
-<img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); \r
- /* Stable prevalence in each health state */\r
- for(cpt=1; cpt<nlstate;cpt++){\r
- fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>\r
-<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);\r
- }\r
- for(cpt=1; cpt<=nlstate;cpt++) {\r
+ for(cpt=1; cpt<=nlstate;cpt++) {\r
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident\r
interval) in state (%d): v%s%d%d.png <br>\r
<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); \r
}\r
- for(cpt=1; cpt<=nlstate;cpt++) {\r
- fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>\r
-<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);\r
- }\r
- fprintf(fichtm,"\n<br>- Total life expectancy by age and\r
-health expectancies in states (1) and (2): e%s%d.png<br>\r
-<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);\r
- }\r
- }\r
+ } /* end i1 */\r
+ }/* End k1 */\r
+ fprintf(fichtm,"</ul>");\r
fclose(fichtm);\r
}\r
\r
int ng;\r
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {\r
printf("Problem with file %s",optionfilegnuplot);\r
+ fprintf(ficlog,"Problem with file %s",optionfilegnuplot);\r
}\r
\r
#ifdef windows\r
for(k=1; k <=(nlstate+ndeath); k++){\r
if (k != i) {\r
for(j=1; j <=ncovmodel; j++){\r
- \r
fprintf(ficgp,"p%d=%f ",jk,p[jk]);\r
jk++; \r
fprintf(ficgp,"\n");\r
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");\r
i=i+ncovmodel;\r
}\r
- }\r
- }\r
- }\r
- }\r
+ } /* end k */\r
+ } /* end k2 */\r
+ } /* end jk */\r
+ } /* end ng */\r
fclose(ficgp); \r
} /* end gnuplot */\r
\r
strcat(fileresf,fileres);\r
if((ficresf=fopen(fileresf,"w"))==NULL) {\r
printf("Problem with forecast resultfile: %s\n", fileresf);\r
+ fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);\r
}\r
printf("Computing forecasting: result on file '%s' \n", fileresf);\r
+ fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);\r
\r
if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
\r
strcat(filerespop,fileres);\r
if((ficrespop=fopen(filerespop,"w"))==NULL) {\r
printf("Problem with forecast resultfile: %s\n", filerespop);\r
+ fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);\r
}\r
printf("Computing forecasting: result on file '%s' \n", filerespop);\r
+ fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);\r
\r
if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
\r
if (popforecast==1) {\r
if((ficpop=fopen(popfile,"r"))==NULL) {\r
printf("Problem with population file : %s\n",popfile);exit(0);\r
+ fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0);\r
} \r
popage=ivector(0,AGESUP);\r
popeffectif=vector(0,AGESUP);\r
double ***p3mat;\r
int *indx;\r
char line[MAXLINE], linepar[MAXLINE];\r
- char path[80],pathc[80],pathcd[80],pathtot[80],model[20];\r
+ char path[80],pathc[80],pathcd[80],pathtot[80],model[80];\r
int firstobs=1, lastobs=10;\r
int sdeb, sfin; /* Status at beginning and end */\r
int c, h , cpt,l;\r
\r
/*-------- arguments in the command line --------*/\r
\r
+ /* Log file */\r
+ strcat(filelog, optionfilefiname);\r
+ strcat(filelog,".log"); /* */\r
+ if((ficlog=fopen(filelog,"w"))==NULL) {\r
+ printf("Problem with logfile %s\n",filelog);\r
+ goto end;\r
+ }\r
+ fprintf(ficlog,"Log filename:%s\n",filelog);\r
+ fprintf(ficlog,"\n%s",version);\r
+ fprintf(ficlog,"\nEnter the parameter file name: ");\r
+ fprintf(ficlog,"pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);\r
+ fflush(ficlog);\r
+\r
+ /* */\r
strcpy(fileres,"r");\r
strcat(fileres, optionfilefiname);\r
strcat(fileres,".txt"); /* Other files have txt extension */\r
\r
if((ficpar=fopen(optionfile,"r"))==NULL) {\r
printf("Problem with optionfile %s\n",optionfile);\r
+ fprintf(ficlog,"Problem with optionfile %s\n",optionfile);\r
goto end;\r
}\r
\r
strcpy(filereso,"o");\r
strcat(filereso,fileres);\r
if((ficparo=fopen(filereso,"w"))==NULL) {\r
- printf("Problem with Output resultfile: %s\n", filereso);goto end;\r
+ printf("Problem with Output resultfile: %s\n", filereso);\r
+ fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);\r
+ goto end;\r
}\r
\r
/* Reads comments: lines beginning with '#' */\r
}\r
ungetc(c,ficpar);\r
\r
- 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);\r
+ 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 model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);\r
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);\r
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);\r
while((c=getc(ficpar))=='#' && c!= EOF){\r
for(j=1; j <=nlstate+ndeath-1; j++){\r
fscanf(ficpar,"%1d%1d",&i1,&j1);\r
fprintf(ficparo,"%1d%1d",i1,j1);\r
- printf("%1d%1d",i,j);\r
+ if(mle==1)\r
+ printf("%1d%1d",i,j);\r
+ fprintf(ficlog,"%1d%1d",i,j);\r
for(k=1; k<=ncovmodel;k++){\r
fscanf(ficpar," %lf",¶m[i][j][k]);\r
- printf(" %lf",param[i][j][k]);\r
+ if(mle==1){\r
+ printf(" %lf",param[i][j][k]);\r
+ fprintf(ficlog," %lf",param[i][j][k]);\r
+ }\r
+ else\r
+ fprintf(ficlog," %lf",param[i][j][k]);\r
fprintf(ficparo," %lf",param[i][j][k]);\r
}\r
fscanf(ficpar,"\n");\r
- printf("\n");\r
+ if(mle==1)\r
+ printf("\n");\r
+ fprintf(ficlog,"\n");\r
fprintf(ficparo,"\n");\r
}\r
\r
matcov=matrix(1,npar,1,npar);\r
for(i=1; i <=npar; i++){\r
fscanf(ficpar,"%s",&str);\r
- printf("%s",str);\r
+ if(mle==1)\r
+ printf("%s",str);\r
+ fprintf(ficlog,"%s",str);\r
fprintf(ficparo,"%s",str);\r
for(j=1; j <=i; j++){\r
fscanf(ficpar," %le",&matcov[i][j]);\r
- printf(" %.5le",matcov[i][j]);\r
+ if(mle==1){\r
+ printf(" %.5le",matcov[i][j]);\r
+ fprintf(ficlog," %.5le",matcov[i][j]);\r
+ }\r
+ else\r
+ fprintf(ficlog," %.5le",matcov[i][j]);\r
fprintf(ficparo," %.5le",matcov[i][j]);\r
}\r
fscanf(ficpar,"\n");\r
- printf("\n");\r
+ if(mle==1)\r
+ printf("\n");\r
+ fprintf(ficlog,"\n");\r
fprintf(ficparo,"\n");\r
}\r
for(i=1; i <=npar; i++)\r
for(j=i+1;j<=npar;j++)\r
matcov[i][j]=matcov[j][i];\r
\r
- printf("\n");\r
+ if(mle==1)\r
+ printf("\n");\r
+ fprintf(ficlog,"\n");\r
\r
\r
/*-------- Rewriting paramater file ----------*/\r
strcat(rfileres,optionfilext); /* Other files have txt extension */\r
if((ficres =fopen(rfileres,"w"))==NULL) {\r
printf("Problem writing new parameter file: %s\n", fileres);goto end;\r
+ fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;\r
}\r
fprintf(ficres,"#%s\n",version);\r
\r
/*-------- data file ----------*/\r
if((fic=fopen(datafile,"r"))==NULL) {\r
printf("Problem with datafile: %s\n", datafile);goto end;\r
+ fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;\r
}\r
\r
n= lastobs;\r
\r
\r
/* Calculation of the number of parameter from char model*/\r
- Tvar=ivector(1,15); \r
+ Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */\r
Tprod=ivector(1,15); \r
Tvaraff=ivector(1,15); \r
Tvard=imatrix(1,15,1,2);\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
+ fprintf(ficlog,"Error. Non available option model=%s ",model);\r
goto end;\r
}\r
\r
for(i=(j+1); i>=1;i--){\r
- cutv(stra,strb,modelsav,'+');\r
- if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); \r
+ cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ \r
+ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */\r
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/\r
/*scanf("%d",i);*/\r
- if (strchr(strb,'*')) {\r
- cutv(strd,strc,strb,'*');\r
- if (strcmp(strc,"age")==0) {\r
+ if (strchr(strb,'*')) { /* Model includes a product */\r
+ cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *age)*/\r
+ if (strcmp(strc,"age")==0) { /* Vn*age */\r
cptcovprod--;\r
cutv(strb,stre,strd,'V');\r
- Tvar[i]=atoi(stre);\r
+ Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/\r
cptcovage++;\r
Tage[cptcovage]=i;\r
/*printf("stre=%s ", stre);*/\r
}\r
- else if (strcmp(strd,"age")==0) {\r
+ else if (strcmp(strd,"age")==0) { /* or age*Vn */\r
cptcovprod--;\r
cutv(strb,stre,strc,'V');\r
Tvar[i]=atoi(stre);\r
cptcovage++;\r
Tage[cptcovage]=i;\r
}\r
- else {\r
- cutv(strb,stre,strc,'V');\r
+ else { /* Age is not in the model */\r
+ cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/\r
Tvar[i]=ncovcol+k1;\r
- cutv(strb,strc,strd,'V'); \r
+ cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */\r
Tprod[k1]=i;\r
- Tvard[k1][1]=atoi(strc);\r
- Tvard[k1][2]=atoi(stre);\r
+ Tvard[k1][1]=atoi(strc); /* m*/\r
+ Tvard[k1][2]=atoi(stre); /* n */\r
Tvar[cptcovn+k2]=Tvard[k1][1];\r
Tvar[cptcovn+k2+1]=Tvard[k1][2]; \r
for (k=1; k<=lastobs;k++) \r
k2=k2+2;\r
}\r
}\r
- else {\r
+ else { /* no more sum */\r
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/\r
/* scanf("%d",i);*/\r
cutv(strd,strc,strb,'V');\r
strcpy(modelsav,stra); \r
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);\r
scanf("%d",i);*/\r
- }\r
-}\r
+ } /* end of loop + */\r
+ } /* end model */\r
\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
+ fprintf(ficlog,"cptcovprod=%d ", cptcovprod);\r
scanf("%d ",i);*/\r
fclose(fic);\r
\r
else {\r
if (andc[i]!=9999){\r
printf("Warning negative age at death: %d line:%d\n",num[i],i);\r
+ fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);\r
agev[m][i]=-1;\r
}\r
}\r
for (i=1; i<=imx; i++) {\r
for(m=1; (m<= maxwav); m++){\r
if (s[m][i] > (nlstate+ndeath)) {\r
- printf("Error: Wrong value in nlstate or ndeath\n"); \r
+ printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); \r
+ fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); \r
goto end;\r
}\r
}\r
}\r
\r
printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);\r
+ fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); \r
\r
free_vector(severity,1,maxwav);\r
free_imatrix(outcome,1,maxwav+1,1,n);\r
jk=1;\r
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");\r
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");\r
+ fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");\r
for(i=1,jk=1; i <=nlstate; i++){\r
for(k=1; k <=(nlstate+ndeath); k++){\r
if (k != i) \r
{\r
printf("%d%d ",i,k);\r
+ fprintf(ficlog,"%d%d ",i,k);\r
fprintf(ficres,"%1d%1d ",i,k);\r
for(j=1; j <=ncovmodel; j++){\r
printf("%f ",p[jk]);\r
+ fprintf(ficlog,"%f ",p[jk]);\r
fprintf(ficres,"%f ",p[jk]);\r
jk++; \r
}\r
printf("\n");\r
+ fprintf(ficlog,"\n");\r
fprintf(ficres,"\n");\r
}\r
}\r
}\r
- if(mle==1){\r
- /* Computing hessian and covariance matrix */\r
- ftolhess=ftol; /* Usually correct */\r
- hesscov(matcov, p, npar, delti, ftolhess, func);\r
- }\r
- fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");\r
- printf("# Scales (for hessian or gradient estimation)\n");\r
- for(i=1,jk=1; i <=nlstate; i++){\r
- for(j=1; j <=nlstate+ndeath; j++){\r
- if (j!=i) {\r
- fprintf(ficres,"%1d%1d",i,j);\r
- printf("%1d%1d",i,j);\r
- for(k=1; k<=ncovmodel;k++){\r
- printf(" %.5e",delti[jk]);\r
- fprintf(ficres," %.5e",delti[jk]);\r
- jk++;\r
- }\r
- printf("\n");\r
- fprintf(ficres,"\n");\r
- }\r
- }\r
+ if(mle==1){\r
+ /* Computing hessian and covariance matrix */\r
+ ftolhess=ftol; /* Usually correct */\r
+ hesscov(matcov, p, npar, delti, ftolhess, func);\r
+ }\r
+ fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");\r
+ printf("# Scales (for hessian or gradient estimation)\n");\r
+ fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");\r
+ for(i=1,jk=1; i <=nlstate; i++){\r
+ for(j=1; j <=nlstate+ndeath; j++){\r
+ if (j!=i) {\r
+ fprintf(ficres,"%1d%1d",i,j);\r
+ printf("%1d%1d",i,j);\r
+ fprintf(ficlog,"%1d%1d",i,j);\r
+ for(k=1; k<=ncovmodel;k++){\r
+ printf(" %.5e",delti[jk]);\r
+ fprintf(ficlog," %.5e",delti[jk]);\r
+ fprintf(ficres," %.5e",delti[jk]);\r
+ jk++;\r
+ }\r
+ printf("\n");\r
+ fprintf(ficlog,"\n");\r
+ fprintf(ficres,"\n");\r
+ }\r
}\r
- \r
- k=1;\r
- 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");\r
- 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");\r
- for(i=1;i<=npar;i++){\r
- /* if (k>nlstate) k=1;\r
- i1=(i-1)/(ncovmodel*nlstate)+1; \r
- fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);\r
- printf("%s%d%d",alph[k],i1,tab[i]);*/\r
- fprintf(ficres,"%3d",i);\r
- printf("%3d",i);\r
- for(j=1; j<=i;j++){\r
- fprintf(ficres," %.5e",matcov[i][j]);\r
- printf(" %.5e",matcov[i][j]);\r
- }\r
- fprintf(ficres,"\n");\r
- printf("\n");\r
- k++;\r
- }\r
- \r
- while((c=getc(ficpar))=='#' && c!= EOF){\r
- ungetc(c,ficpar);\r
- fgets(line, MAXLINE, ficpar);\r
- puts(line);\r
- fputs(line,ficparo);\r
- }\r
- ungetc(c,ficpar);\r
- estepm=0;\r
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);\r
- if (estepm==0 || estepm < stepm) estepm=stepm;\r
- if (fage <= 2) {\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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);\r
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);\r
- \r
- while((c=getc(ficpar))=='#' && c!= EOF){\r
- ungetc(c,ficpar);\r
- fgets(line, MAXLINE, ficpar);\r
- puts(line);\r
- fputs(line,ficparo);\r
- }\r
- ungetc(c,ficpar);\r
+ }\r
+ \r
+ k=1;\r
+ 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");\r
+ if(mle==1)\r
+ 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");\r
+ fprintf(ficlog,"# 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");\r
+ for(i=1;i<=npar;i++){\r
+ /* if (k>nlstate) k=1;\r
+ i1=(i-1)/(ncovmodel*nlstate)+1; \r
+ fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);\r
+ printf("%s%d%d",alph[k],i1,tab[i]);*/\r
+ fprintf(ficres,"%3d",i);\r
+ if(mle==1)\r
+ printf("%3d",i);\r
+ fprintf(ficlog,"%3d",i);\r
+ for(j=1; j<=i;j++){\r
+ fprintf(ficres," %.5e",matcov[i][j]);\r
+ if(mle==1)\r
+ printf(" %.5e",matcov[i][j]);\r
+ fprintf(ficlog," %.5e",matcov[i][j]);\r
+ }\r
+ fprintf(ficres,"\n");\r
+ if(mle==1)\r
+ printf("\n");\r
+ fprintf(ficlog,"\n");\r
+ k++;\r
+ }\r
+ \r
+ while((c=getc(ficpar))=='#' && c!= EOF){\r
+ ungetc(c,ficpar);\r
+ fgets(line, MAXLINE, ficpar);\r
+ puts(line);\r
+ fputs(line,ficparo);\r
+ }\r
+ ungetc(c,ficpar);\r
+ estepm=0;\r
+ fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);\r
+ if (estepm==0 || estepm < stepm) estepm=stepm;\r
+ if (fage <= 2) {\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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);\r
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);\r
+ \r
+ while((c=getc(ficpar))=='#' && c!= EOF){\r
+ ungetc(c,ficpar);\r
+ fgets(line, MAXLINE, ficpar);\r
+ puts(line);\r
+ fputs(line,ficparo);\r
+ }\r
+ ungetc(c,ficpar);\r
\r
- fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);\r
- fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
- fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
- \r
- while((c=getc(ficpar))=='#' && c!= EOF){\r
- ungetc(c,ficpar);\r
- fgets(line, MAXLINE, ficpar);\r
- puts(line);\r
- fputs(line,ficparo);\r
- }\r
- ungetc(c,ficpar);\r
+ fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);\r
+ fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
+ fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);\r
+ \r
+ while((c=getc(ficpar))=='#' && c!= EOF){\r
+ ungetc(c,ficpar);\r
+ fgets(line, MAXLINE, ficpar);\r
+ puts(line);\r
+ fputs(line,ficparo);\r
+ }\r
+ ungetc(c,ficpar);\r
\r
\r
dateprev1=anprev1+mprev1/12.+jprev1/365.;\r
<hr size=\"2\" color=\"#EC5E5E\">\r
<ul><li><h4>Parameter files</h4>\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></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);\r
+ - Log file of the run: <a href=\"%s\">%s</a><br>\n\r
+ - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);\r
fclose(fichtm);\r
\r
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);\r
strcat(filerespl,fileres);\r
if((ficrespl=fopen(filerespl,"w"))==NULL) {\r
printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;\r
+ fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;\r
}\r
printf("Computing prevalence limit: result on file '%s' \n", filerespl);\r
+ fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);\r
fprintf(ficrespl,"#Prevalence limit\n");\r
fprintf(ficrespl,"#Age ");\r
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);\r
k=k+1;\r
/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/\r
fprintf(ficrespl,"\n#******");\r
- for(j=1;j<=cptcoveff;j++) \r
+ printf("\n#******");\r
+ fprintf(ficlog,"\n#******");\r
+ for(j=1;j<=cptcoveff;j++) {\r
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+ printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+ fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+ }\r
fprintf(ficrespl,"******\n");\r
+ printf("******\n");\r
+ fprintf(ficlog,"******\n");\r
\r
for (age=agebase; age<=agelim; age++){\r
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);\r
strcpy(filerespij,"pij"); strcat(filerespij,fileres);\r
if((ficrespij=fopen(filerespij,"w"))==NULL) {\r
printf("Problem with Pij resultfile: %s\n", filerespij);goto end;\r
+ fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;\r
}\r
printf("Computing pij: result on file '%s' \n", filerespij);\r
+ fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);\r
\r
stepsize=(int) (stepm+YEARM-1)/YEARM;\r
/*if (stepm<=24) stepsize=2;*/\r
else{\r
erreur=108;\r
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);\r
+ fprintf(ficlog,"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);\r
}\r
\r
\r
strcat(filerest,fileres);\r
if((ficrest=fopen(filerest,"w"))==NULL) {\r
printf("Problem with total LE resultfile: %s\n", filerest);goto end;\r
+ fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;\r
}\r
printf("Computing Total LEs with variances: file '%s' \n", filerest); \r
+ fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); \r
\r
\r
strcpy(filerese,"e");\r
strcat(filerese,fileres);\r
if((ficreseij=fopen(filerese,"w"))==NULL) {\r
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);\r
+ fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);\r
}\r
printf("Computing Health Expectancies: result on file '%s' \n", filerese);\r
+ fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);\r
\r
- strcpy(fileresv,"v");\r
+ strcpy(fileresv,"v");\r
strcat(fileresv,fileres);\r
if((ficresvij=fopen(fileresv,"w"))==NULL) {\r
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);\r
+ fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);\r
}\r
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);\r
+ fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);\r
calagedate=-1;\r
-prevalence(ageminpar, 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
k=0;\r
for(cptcov=1;cptcov<=i1;cptcov++){\r
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){\r
- k=k+1;\r
+ k=k+1; \r
fprintf(ficrest,"\n#****** ");\r
for(j=1;j<=cptcoveff;j++) \r
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
\r
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);\r
oldm=oldms;savm=savms;\r
- varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);\r
- \r
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0);\r
+ if(popbased==1){\r
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased);\r
+ }\r
\r
\r
fprintf(ficrest,"#Total LEs with variances: e.. (std) ");\r
fclose(ficgp);\r
\r
\r
- if(erreur >0)\r
+ if(erreur >0){\r
printf("End of Imach with error or warning %d\n",erreur);\r
- else printf("End of Imach\n");\r
+ fprintf(ficlog,"End of Imach with error or warning %d\n",erreur);\r
+ }else{\r
+ printf("End of Imach\n");\r
+ fprintf(ficlog,"End of Imach\n");\r
+ }\r
+ printf("See log file on %s\n",filelog);\r
+ fclose(ficlog);\r
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */\r
\r
/* 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);*/\r
git://henry.ined.fr / .git / commitdiff