version 1.48, 2002/06/10 13:12:49
|
version 1.50, 2002/06/26 23:25:02
|
Line 77
|
Line 77
|
#define AGEBASE 40
|
#define AGEBASE 40
|
#ifdef windows
|
#ifdef windows
|
#define DIRSEPARATOR '\\'
|
#define DIRSEPARATOR '\\'
|
|
#define ODIRSEPARATOR '/'
|
#else
|
#else
|
#define DIRSEPARATOR '/'
|
#define DIRSEPARATOR '/'
|
|
#define ODIRSEPARATOR '\\'
|
#endif
|
#endif
|
|
|
char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES ";
|
char version[80]="Imach version 0.8i, June 2002, INED-EUROREVES ";
|
int erreur; /* Error number */
|
int erreur; /* Error number */
|
int nvar;
|
int nvar;
|
int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
|
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
|
int npar=NPARMAX;
|
int npar=NPARMAX;
|
int nlstate=2; /* Number of live states */
|
int nlstate=2; /* Number of live states */
|
int ndeath=1; /* Number of dead states */
|
int ndeath=1; /* Number of dead states */
|
Line 101 double jmean; /* Mean space between 2 wa
|
Line 103 double jmean; /* Mean space between 2 wa
|
double **oldm, **newm, **savm; /* Working pointers to matrices */
|
double **oldm, **newm, **savm; /* Working pointers to matrices */
|
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
|
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
|
|
FILE *ficlog;
|
FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
|
FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
|
|
FILE *ficresprobmorprev;
|
FILE *fichtm; /* Html File */
|
FILE *fichtm; /* Html File */
|
FILE *ficreseij;
|
FILE *ficreseij;
|
char filerese[FILENAMELENGTH];
|
char filerese[FILENAMELENGTH];
|
Line 114 char optionfile[FILENAMELENGTH], datafil
|
Line 118 char optionfile[FILENAMELENGTH], datafil
|
char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
|
char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
|
|
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
|
|
char filelog[FILENAMELENGTH]; /* Log file */
|
char filerest[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char popfile[FILENAMELENGTH];
|
char popfile[FILENAMELENGTH];
|
Line 178 static int split( char *path, char *dirc
|
Line 182 static int split( char *path, char *dirc
|
|
|
l1 = strlen( path ); /* length of path */
|
l1 = strlen( path ); /* length of path */
|
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
|
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
|
s = strrchr( path, DIRSEPARATOR ); /* find last / */
|
s= strrchr( path, DIRSEPARATOR ); /* find last / */
|
if ( s == NULL ) { /* no directory, so use current */
|
if ( s == NULL ) { /* no directory, so use current */
|
|
/*if(strrchr(path, ODIRSEPARATOR )==NULL)
|
|
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
|
#if defined(__bsd__) /* get current working directory */
|
#if defined(__bsd__) /* get current working directory */
|
extern char *getwd( );
|
extern char *getwd( );
|
|
|
Line 245 int nbocc(char *s, char occ)
|
Line 251 int nbocc(char *s, char occ)
|
|
|
void cutv(char *u,char *v, char*t, char occ)
|
void cutv(char *u,char *v, char*t, char occ)
|
{
|
{
|
|
/* cuts string t into u and v where u is ended by char occ excluding it
|
|
and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)
|
|
gives u="abcedf" and v="ghi2j" */
|
int i,lg,j,p=0;
|
int i,lg,j,p=0;
|
i=0;
|
i=0;
|
for(j=0; j<=strlen(t)-1; j++) {
|
for(j=0; j<=strlen(t)-1; j++) {
|
Line 441 double brent(double ax, double bx, doubl
|
Line 450 double brent(double ax, double bx, doubl
|
tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
|
tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
|
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
|
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
|
printf(".");fflush(stdout);
|
printf(".");fflush(stdout);
|
|
fprintf(ficlog,".");fflush(ficlog);
|
#ifdef DEBUG
|
#ifdef DEBUG
|
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);
|
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);
|
|
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);
|
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
|
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
|
#endif
|
#endif
|
if (fabs(x-xm) <= (tol2-0.5*(b-a))){
|
if (fabs(x-xm) <= (tol2-0.5*(b-a))){
|
Line 567 void linmin(double p[], double xi[], int
|
Line 578 void linmin(double p[], double xi[], int
|
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
|
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
|
#ifdef DEBUG
|
#ifdef DEBUG
|
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
|
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
|
|
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
|
#endif
|
#endif
|
for (j=1;j<=n;j++) {
|
for (j=1;j<=n;j++) {
|
xi[j] *= xmin;
|
xi[j] *= xmin;
|
Line 597 void powell(double p[], double **xi, int
|
Line 609 void powell(double p[], double **xi, int
|
ibig=0;
|
ibig=0;
|
del=0.0;
|
del=0.0;
|
printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
|
printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
|
|
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
|
for (i=1;i<=n;i++)
|
for (i=1;i<=n;i++)
|
printf(" %d %.12f",i, p[i]);
|
printf(" %d %.12f",i, p[i]);
|
|
fprintf(ficlog," %d %.12f",i, p[i]);
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
for (i=1;i<=n;i++) {
|
for (i=1;i<=n;i++) {
|
for (j=1;j<=n;j++) xit[j]=xi[j][i];
|
for (j=1;j<=n;j++) xit[j]=xi[j][i];
|
fptt=(*fret);
|
fptt=(*fret);
|
#ifdef DEBUG
|
#ifdef DEBUG
|
printf("fret=%lf \n",*fret);
|
printf("fret=%lf \n",*fret);
|
|
fprintf(ficlog,"fret=%lf \n",*fret);
|
#endif
|
#endif
|
printf("%d",i);fflush(stdout);
|
printf("%d",i);fflush(stdout);
|
|
fprintf(ficlog,"%d",i);fflush(ficlog);
|
linmin(p,xit,n,fret,func);
|
linmin(p,xit,n,fret,func);
|
if (fabs(fptt-(*fret)) > del) {
|
if (fabs(fptt-(*fret)) > del) {
|
del=fabs(fptt-(*fret));
|
del=fabs(fptt-(*fret));
|
Line 614 void powell(double p[], double **xi, int
|
Line 631 void powell(double p[], double **xi, int
|
}
|
}
|
#ifdef DEBUG
|
#ifdef DEBUG
|
printf("%d %.12e",i,(*fret));
|
printf("%d %.12e",i,(*fret));
|
|
fprintf(ficlog,"%d %.12e",i,(*fret));
|
for (j=1;j<=n;j++) {
|
for (j=1;j<=n;j++) {
|
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
|
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
|
printf(" x(%d)=%.12e",j,xit[j]);
|
printf(" x(%d)=%.12e",j,xit[j]);
|
|
fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
|
}
|
}
|
for(j=1;j<=n;j++)
|
for(j=1;j<=n;j++) {
|
printf(" p=%.12e",p[j]);
|
printf(" p=%.12e",p[j]);
|
|
fprintf(ficlog," p=%.12e",p[j]);
|
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
#endif
|
#endif
|
}
|
}
|
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
|
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
|
Line 629 void powell(double p[], double **xi, int
|
Line 651 void powell(double p[], double **xi, int
|
k[0]=1;
|
k[0]=1;
|
k[1]=-1;
|
k[1]=-1;
|
printf("Max: %.12e",(*func)(p));
|
printf("Max: %.12e",(*func)(p));
|
for (j=1;j<=n;j++)
|
fprintf(ficlog,"Max: %.12e",(*func)(p));
|
|
for (j=1;j<=n;j++) {
|
printf(" %.12e",p[j]);
|
printf(" %.12e",p[j]);
|
|
fprintf(ficlog," %.12e",p[j]);
|
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
for(l=0;l<=1;l++) {
|
for(l=0;l<=1;l++) {
|
for (j=1;j<=n;j++) {
|
for (j=1;j<=n;j++) {
|
ptt[j]=p[j]+(p[j]-pt[j])*k[l];
|
ptt[j]=p[j]+(p[j]-pt[j])*k[l];
|
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
|
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
|
|
fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
|
}
|
}
|
printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
|
printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
|
|
fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
|
}
|
}
|
#endif
|
#endif
|
|
|
Line 665 void powell(double p[], double **xi, int
|
Line 693 void powell(double p[], double **xi, int
|
}
|
}
|
#ifdef DEBUG
|
#ifdef DEBUG
|
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
|
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
|
for(j=1;j<=n;j++)
|
fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
|
|
for(j=1;j<=n;j++){
|
printf(" %.12e",xit[j]);
|
printf(" %.12e",xit[j]);
|
|
fprintf(ficlog," %.12e",xit[j]);
|
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
#endif
|
#endif
|
}
|
}
|
}
|
}
|
Line 938 void mlikeli(FILE *ficres,double p[], in
|
Line 970 void mlikeli(FILE *ficres,double p[], in
|
for (i=1;i<=npar;i++)
|
for (i=1;i<=npar;i++)
|
for (j=1;j<=npar;j++)
|
for (j=1;j<=npar;j++)
|
xi[i][j]=(i==j ? 1.0 : 0.0);
|
xi[i][j]=(i==j ? 1.0 : 0.0);
|
printf("Powell\n");
|
printf("Powell\n"); fprintf(ficlog,"Powell\n");
|
powell(p,xi,npar,ftol,&iter,&fret,func);
|
powell(p,xi,npar,ftol,&iter,&fret,func);
|
|
|
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
|
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
|
|
fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
|
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
|
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
|
|
|
}
|
}
|
Line 962 void hesscov(double **matcov, double p[]
|
Line 995 void hesscov(double **matcov, double p[]
|
hess=matrix(1,npar,1,npar);
|
hess=matrix(1,npar,1,npar);
|
|
|
printf("\nCalculation of the hessian matrix. Wait...\n");
|
printf("\nCalculation of the hessian matrix. Wait...\n");
|
|
fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n");
|
for (i=1;i<=npar;i++){
|
for (i=1;i<=npar;i++){
|
printf("%d",i);fflush(stdout);
|
printf("%d",i);fflush(stdout);
|
|
fprintf(ficlog,"%d",i);fflush(ficlog);
|
hess[i][i]=hessii(p,ftolhess,i,delti);
|
hess[i][i]=hessii(p,ftolhess,i,delti);
|
/*printf(" %f ",p[i]);*/
|
/*printf(" %f ",p[i]);*/
|
/*printf(" %lf ",hess[i][i]);*/
|
/*printf(" %lf ",hess[i][i]);*/
|
Line 973 void hesscov(double **matcov, double p[]
|
Line 1008 void hesscov(double **matcov, double p[]
|
for (j=1;j<=npar;j++) {
|
for (j=1;j<=npar;j++) {
|
if (j>i) {
|
if (j>i) {
|
printf(".%d%d",i,j);fflush(stdout);
|
printf(".%d%d",i,j);fflush(stdout);
|
|
fprintf(ficlog,".%d%d",i,j);fflush(ficlog);
|
hess[i][j]=hessij(p,delti,i,j);
|
hess[i][j]=hessij(p,delti,i,j);
|
hess[j][i]=hess[i][j];
|
hess[j][i]=hess[i][j];
|
/*printf(" %lf ",hess[i][j]);*/
|
/*printf(" %lf ",hess[i][j]);*/
|
Line 980 void hesscov(double **matcov, double p[]
|
Line 1016 void hesscov(double **matcov, double p[]
|
}
|
}
|
}
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
|
|
printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
|
printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
|
|
fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");
|
|
|
a=matrix(1,npar,1,npar);
|
a=matrix(1,npar,1,npar);
|
y=matrix(1,npar,1,npar);
|
y=matrix(1,npar,1,npar);
|
Line 1001 void hesscov(double **matcov, double p[]
|
Line 1039 void hesscov(double **matcov, double p[]
|
}
|
}
|
|
|
printf("\n#Hessian matrix#\n");
|
printf("\n#Hessian matrix#\n");
|
|
fprintf(ficlog,"\n#Hessian matrix#\n");
|
for (i=1;i<=npar;i++) {
|
for (i=1;i<=npar;i++) {
|
for (j=1;j<=npar;j++) {
|
for (j=1;j<=npar;j++) {
|
printf("%.3e ",hess[i][j]);
|
printf("%.3e ",hess[i][j]);
|
|
fprintf(ficlog,"%.3e ",hess[i][j]);
|
}
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
}
|
}
|
|
|
/* Recompute Inverse */
|
/* Recompute Inverse */
|
Line 1022 void hesscov(double **matcov, double p[]
|
Line 1063 void hesscov(double **matcov, double p[]
|
for (i=1;i<=npar;i++){
|
for (i=1;i<=npar;i++){
|
y[i][j]=x[i];
|
y[i][j]=x[i];
|
printf("%.3e ",y[i][j]);
|
printf("%.3e ",y[i][j]);
|
|
fprintf(ficlog,"%.3e ",y[i][j]);
|
}
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
}
|
}
|
*/
|
*/
|
|
|
Line 1065 double hessii( double x[], double delta,
|
Line 1108 double hessii( double x[], double delta,
|
|
|
#ifdef DEBUG
|
#ifdef DEBUG
|
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);
|
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);
|
|
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);
|
#endif
|
#endif
|
/*if(fabs(k1-2.0*fx+k2) <1.e-13){ */
|
/*if(fabs(k1-2.0*fx+k2) <1.e-13){ */
|
if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)){
|
if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)){
|
Line 1112 double hessij( double x[], double delti[
|
Line 1156 double hessij( double x[], double delti[
|
res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */
|
res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */
|
#ifdef DEBUG
|
#ifdef DEBUG
|
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);
|
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);
|
|
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);
|
#endif
|
#endif
|
}
|
}
|
return res;
|
return res;
|
Line 1196 void freqsummary(char fileres[], int ag
|
Line 1241 void freqsummary(char fileres[], int ag
|
{ /* 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;
|
|
int first;
|
double ***freq; /* Frequencies */
|
double ***freq; /* Frequencies */
|
double *pp;
|
double *pp;
|
double pos, k2, dateintsum=0,k2cpt=0;
|
double pos, k2, dateintsum=0,k2cpt=0;
|
Line 1208 void freqsummary(char fileres[], int ag
|
Line 1254 void freqsummary(char fileres[], int ag
|
strcat(fileresp,fileres);
|
strcat(fileresp,fileres);
|
if((ficresp=fopen(fileresp,"w"))==NULL) {
|
if((ficresp=fopen(fileresp,"w"))==NULL) {
|
printf("Problem with prevalence resultfile: %s\n", fileresp);
|
printf("Problem with prevalence resultfile: %s\n", fileresp);
|
|
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
|
exit(0);
|
exit(0);
|
}
|
}
|
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
|
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
|
Line 1215 void freqsummary(char fileres[], int ag
|
Line 1262 void freqsummary(char fileres[], int ag
|
|
|
j=cptcoveff;
|
j=cptcoveff;
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
|
|
|
first=1;
|
|
|
for(k1=1; k1<=j;k1++){
|
for(k1=1; k1<=j;k1++){
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
j1++;
|
j1++;
|
Line 1267 void freqsummary(char fileres[], int ag
|
Line 1316 void freqsummary(char fileres[], int ag
|
fprintf(ficresp, "\n");
|
fprintf(ficresp, "\n");
|
|
|
for(i=(int)agemin; i <= (int)agemax+3; i++){
|
for(i=(int)agemin; i <= (int)agemax+3; i++){
|
if(i==(int)agemax+3)
|
if(i==(int)agemax+3){
|
printf("Total");
|
fprintf(ficlog,"Total");
|
else
|
}else{
|
printf("Age %d", i);
|
if(first==1){
|
|
first=0;
|
|
printf("See log file for details...\n");
|
|
}
|
|
fprintf(ficlog,"Age %d", i);
|
|
}
|
for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
pp[jk] += freq[jk][m][i];
|
pp[jk] += freq[jk][m][i];
|
Line 1278 void freqsummary(char fileres[], int ag
|
Line 1332 void freqsummary(char fileres[], int ag
|
for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
for(m=-1, pos=0; m <=0 ; m++)
|
for(m=-1, pos=0; m <=0 ; m++)
|
pos += freq[jk][m][i];
|
pos += freq[jk][m][i];
|
if(pp[jk]>=1.e-10)
|
if(pp[jk]>=1.e-10){
|
|
if(first==1){
|
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
|
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
|
else
|
}
|
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
|
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
|
|
}else{
|
|
if(first==1)
|
|
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
|
|
fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
|
|
}
|
}
|
}
|
|
|
for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
Line 1292 void freqsummary(char fileres[], int ag
|
Line 1352 void freqsummary(char fileres[], int ag
|
for(jk=1,pos=0; jk <=nlstate ; jk++)
|
for(jk=1,pos=0; jk <=nlstate ; jk++)
|
pos += pp[jk];
|
pos += pp[jk];
|
for(jk=1; jk <=nlstate ; jk++){
|
for(jk=1; jk <=nlstate ; jk++){
|
if(pos>=1.e-5)
|
if(pos>=1.e-5){
|
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
|
if(first==1)
|
else
|
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
|
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
|
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
|
|
}else{
|
|
if(first==1)
|
|
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
|
|
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
|
|
}
|
if( i <= (int) agemax){
|
if( i <= (int) agemax){
|
if(pos>=1.e-5){
|
if(pos>=1.e-5){
|
fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
|
fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
|
Line 1309 void freqsummary(char fileres[], int ag
|
Line 1374 void freqsummary(char fileres[], int ag
|
|
|
for(jk=-1; jk <=nlstate+ndeath; jk++)
|
for(jk=-1; jk <=nlstate+ndeath; jk++)
|
for(m=-1; m <=nlstate+ndeath; m++)
|
for(m=-1; m <=nlstate+ndeath; m++)
|
if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
|
if(freq[jk][m][i] !=0 ) {
|
|
if(first==1)
|
|
printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
|
|
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
|
|
}
|
if(i <= (int) agemax)
|
if(i <= (int) agemax)
|
fprintf(ficresp,"\n");
|
fprintf(ficresp,"\n");
|
printf("\n");
|
if(first==1)
|
|
printf("Others in log...\n");
|
|
fprintf(ficlog,"\n");
|
}
|
}
|
}
|
}
|
}
|
}
|
Line 1399 void prevalence(int agemin, float agemax
|
Line 1470 void prevalence(int agemin, float agemax
|
probs[i][jk][j1]= pp[jk]/pos;
|
probs[i][jk][j1]= pp[jk]/pos;
|
}
|
}
|
}
|
}
|
}
|
}/* end jk */
|
|
}/* end i */
|
}
|
} /* end i1 */
|
}
|
} /* end k1 */
|
}
|
|
|
|
|
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
Line 1425 void concatwav(int wav[], int **dh, int
|
Line 1495 void concatwav(int wav[], int **dh, int
|
int i, mi, m;
|
int i, mi, m;
|
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
|
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
|
double sum=0., jmean=0.;*/
|
double sum=0., jmean=0.;*/
|
|
int first;
|
int j, k=0,jk, ju, jl;
|
int j, k=0,jk, ju, jl;
|
double sum=0.;
|
double sum=0.;
|
|
first=0;
|
jmin=1e+5;
|
jmin=1e+5;
|
jmax=-1;
|
jmax=-1;
|
jmean=0.;
|
jmean=0.;
|
Line 1450 void concatwav(int wav[], int **dh, int
|
Line 1521 void concatwav(int wav[], int **dh, int
|
}
|
}
|
|
|
wav[i]=mi;
|
wav[i]=mi;
|
if(mi==0)
|
if(mi==0){
|
printf("Warning, no any valid information for:%d line=%d\n",num[i],i);
|
if(first==0){
|
|
printf("Warning, no any valid information for:%d line=%d and may be others, see log file\n",num[i],i);
|
|
first=1;
|
|
}
|
|
if(first==1){
|
|
fprintf(ficlog,"Warning, no any valid information for:%d line=%d\n",num[i],i);
|
|
}
|
|
} /* end mi==0 */
|
}
|
}
|
|
|
for(i=1; i<=imx; i++){
|
for(i=1; i<=imx; i++){
|
Line 1492 void concatwav(int wav[], int **dh, int
|
Line 1570 void concatwav(int wav[], int **dh, int
|
}
|
}
|
jmean=sum/k;
|
jmean=sum/k;
|
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
|
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
|
|
fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
|
}
|
}
|
|
|
/*********** Tricode ****************************/
|
/*********** Tricode ****************************/
|
void tricode(int *Tvar, int **nbcode, int imx)
|
void tricode(int *Tvar, int **nbcode, int imx)
|
{
|
{
|
Line 1532 void tricode(int *Tvar, int **nbcode, in
|
Line 1612 void tricode(int *Tvar, int **nbcode, in
|
for (k=0; k<19; k++) Ndum[k]=0;
|
for (k=0; k<19; k++) Ndum[k]=0;
|
|
|
for (i=1; i<=ncovmodel-2; i++) {
|
for (i=1; i<=ncovmodel-2; i++) {
|
ij=Tvar[i];
|
ij=Tvar[i];
|
Ndum[ij]++;
|
Ndum[ij]++;
|
}
|
}
|
|
|
ij=1;
|
ij=1;
|
for (i=1; i<=10; i++) {
|
for (i=1; i<=10; i++) {
|
Line 1544 void tricode(int *Tvar, int **nbcode, in
|
Line 1624 void tricode(int *Tvar, int **nbcode, in
|
}
|
}
|
}
|
}
|
|
|
cptcoveff=ij-1;
|
cptcoveff=ij-1;
|
}
|
}
|
|
|
/*********** Health Expectancies ****************/
|
/*********** Health Expectancies ****************/
|
Line 1675 void evsij(char fileres[], double ***eij
|
Line 1755 void evsij(char fileres[], double ***eij
|
varhe[i][j][(int)age] =0.;
|
varhe[i][j][(int)age] =0.;
|
|
|
printf("%d|",(int)age);fflush(stdout);
|
printf("%d|",(int)age);fflush(stdout);
|
|
fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
|
for(h=0;h<=nhstepm-1;h++){
|
for(h=0;h<=nhstepm-1;h++){
|
for(k=0;k<=nhstepm-1;k++){
|
for(k=0;k<=nhstepm-1;k++){
|
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
|
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
|
Line 1710 void evsij(char fileres[], double ***eij
|
Line 1791 void evsij(char fileres[], double ***eij
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
}
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
|
|
free_vector(xp,1,npar);
|
free_vector(xp,1,npar);
|
free_matrix(dnewm,1,nlstate*2,1,npar);
|
free_matrix(dnewm,1,nlstate*2,1,npar);
|
Line 1718 void evsij(char fileres[], double ***eij
|
Line 1800 void evsij(char fileres[], double ***eij
|
}
|
}
|
|
|
/************ Variance ******************/
|
/************ Variance ******************/
|
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)
|
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)
|
{
|
{
|
/* Variance of health expectancies */
|
/* Variance of health expectancies */
|
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
|
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
|
double **newm;
|
/* double **newm;*/
|
double **dnewm,**doldm;
|
double **dnewm,**doldm;
|
|
double **dnewmp,**doldmp;
|
int i, j, nhstepm, hstepm, h, nstepm ;
|
int i, j, nhstepm, hstepm, h, nstepm ;
|
int k, cptcode;
|
int k, cptcode;
|
double *xp;
|
double *xp;
|
double **gp, **gm;
|
double **gp, **gm; /* for var eij */
|
double ***gradg, ***trgradg;
|
double ***gradg, ***trgradg; /*for var eij */
|
|
double **gradgp, **trgradgp; /* for var p point j */
|
|
double *gpp, *gmp; /* for var p point j */
|
|
double **varppt; /* for var p point j nlstate to nlstate+ndeath */
|
double ***p3mat;
|
double ***p3mat;
|
double age,agelim, hf;
|
double age,agelim, hf;
|
int theta;
|
int theta;
|
|
char digit[4];
|
|
char digitp[16];
|
|
|
|
char fileresprobmorprev[FILENAMELENGTH];
|
|
|
|
if(popbased==1)
|
|
strcpy(digitp,"-populbased-");
|
|
else
|
|
strcpy(digitp,"-stablbased-");
|
|
|
|
strcpy(fileresprobmorprev,"prmorprev");
|
|
sprintf(digit,"%-d",ij);
|
|
/*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
|
|
strcat(fileresprobmorprev,digit); /* Tvar to be done */
|
|
strcat(fileresprobmorprev,digitp); /* Popbased or not */
|
|
strcat(fileresprobmorprev,fileres);
|
|
if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
|
|
printf("Problem with resultfile: %s\n", fileresprobmorprev);
|
|
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
|
|
}
|
|
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
|
|
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
|
|
fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1*p1j+w2*p2j+... stand dev in()\n");
|
|
fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
|
|
for(j=nlstate+1; j<=(nlstate+ndeath);j++){
|
|
fprintf(ficresprobmorprev," p.%-d SE",j);
|
|
for(i=1; i<=nlstate;i++)
|
|
fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
|
|
}
|
|
fprintf(ficresprobmorprev,"\n");
|
|
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
|
|
printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
|
|
fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);
|
|
exit(0);
|
|
}
|
|
else{
|
|
fprintf(ficgp,"\n# Routine varevsij");
|
|
}
|
|
if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
|
|
printf("Problem with html file: %s\n", optionfilehtm);
|
|
fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
|
|
exit(0);
|
|
}
|
|
else{
|
|
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");
|
|
}
|
|
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");
|
fprintf(ficresvij,"# Age");
|
fprintf(ficresvij,"# Age");
|
Line 1743 void varevsij(char fileres[], double ***
|
Line 1876 void varevsij(char fileres[], double ***
|
xp=vector(1,npar);
|
xp=vector(1,npar);
|
dnewm=matrix(1,nlstate,1,npar);
|
dnewm=matrix(1,nlstate,1,npar);
|
doldm=matrix(1,nlstate,1,nlstate);
|
doldm=matrix(1,nlstate,1,nlstate);
|
|
dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
|
|
doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
|
|
|
|
gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
|
|
gpp=vector(nlstate+1,nlstate+ndeath);
|
|
gmp=vector(nlstate+1,nlstate+ndeath);
|
|
trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
|
|
|
if(estepm < stepm){
|
if(estepm < stepm){
|
printf ("Problem %d lower than %d\n",estepm, stepm);
|
printf ("Problem %d lower than %d\n",estepm, stepm);
|
Line 1770 void varevsij(char fileres[], double ***
|
Line 1910 void varevsij(char fileres[], double ***
|
gp=matrix(0,nhstepm,1,nlstate);
|
gp=matrix(0,nhstepm,1,nlstate);
|
gm=matrix(0,nhstepm,1,nlstate);
|
gm=matrix(0,nhstepm,1,nlstate);
|
|
|
|
|
for(theta=1; theta <=npar; theta++){
|
for(theta=1; theta <=npar; theta++){
|
for(i=1; i<=npar; i++){ /* Computes gradient */
|
for(i=1; i<=npar; i++){ /* Computes gradient */
|
xp[i] = x[i] + (i==theta ?delti[theta]:0);
|
xp[i] = x[i] + (i==theta ?delti[theta]:0);
|
Line 1788 void varevsij(char fileres[], double ***
|
Line 1929 void varevsij(char fileres[], double ***
|
gp[h][j] += prlim[i][i]*p3mat[i][j][h];
|
gp[h][j] += prlim[i][i]*p3mat[i][j][h];
|
}
|
}
|
}
|
}
|
|
/* This for computing forces of mortality (h=1)as a weighted average */
|
|
for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){
|
|
for(i=1; i<= nlstate; i++)
|
|
gpp[j] += prlim[i][i]*p3mat[i][j][1];
|
|
}
|
|
/* end force of mortality */
|
|
|
for(i=1; i<=npar; i++) /* Computes gradient */
|
for(i=1; i<=npar; i++) /* Computes gradient */
|
xp[i] = x[i] - (i==theta ?delti[theta]:0);
|
xp[i] = x[i] - (i==theta ?delti[theta]:0);
|
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
|
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
|
Line 1805 void varevsij(char fileres[], double ***
|
Line 1952 void varevsij(char fileres[], double ***
|
gm[h][j] += prlim[i][i]*p3mat[i][j][h];
|
gm[h][j] += prlim[i][i]*p3mat[i][j][h];
|
}
|
}
|
}
|
}
|
|
/* This for computing force of mortality (h=1)as a weighted average */
|
|
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
|
|
for(i=1; i<= nlstate; i++)
|
|
gmp[j] += prlim[i][i]*p3mat[i][j][1];
|
|
}
|
|
/* end force of mortality */
|
|
|
for(j=1; j<= nlstate; j++)
|
for(j=1; j<= nlstate; j++) /* vareij */
|
for(h=0; h<=nhstepm; h++){
|
for(h=0; h<=nhstepm; h++){
|
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
|
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
|
}
|
}
|
|
for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
|
|
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
|
|
}
|
|
|
} /* End theta */
|
} /* End theta */
|
|
|
trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
|
trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
|
|
|
for(h=0; h<=nhstepm; h++)
|
for(h=0; h<=nhstepm; h++) /* veij */
|
for(j=1; j<=nlstate;j++)
|
for(j=1; j<=nlstate;j++)
|
for(theta=1; theta <=npar; theta++)
|
for(theta=1; theta <=npar; theta++)
|
trgradg[h][j][theta]=gradg[h][theta][j];
|
trgradg[h][j][theta]=gradg[h][theta][j];
|
|
|
|
for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
|
|
for(theta=1; theta <=npar; theta++)
|
|
trgradgp[j][theta]=gradgp[theta][j];
|
|
|
hf=hstepm*stepm/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(i=1;i<=nlstate;i++)
|
for(j=1;j<=nlstate;j++)
|
for(j=1;j<=nlstate;j++)
|
Line 1834 void varevsij(char fileres[], double ***
|
Line 1995 void varevsij(char fileres[], double ***
|
}
|
}
|
}
|
}
|
|
|
|
/* pptj */
|
|
matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
|
|
matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
|
|
for(j=nlstate+1;j<=nlstate+ndeath;j++)
|
|
for(i=nlstate+1;i<=nlstate+ndeath;i++)
|
|
varppt[j][i]=doldmp[j][i];
|
|
/* end ppptj */
|
|
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
|
|
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
|
|
|
|
if (popbased==1) {
|
|
for(i=1; i<=nlstate;i++)
|
|
prlim[i][i]=probs[(int)age][i][ij];
|
|
}
|
|
|
|
/* This for computing force of mortality (h=1)as a weighted average */
|
|
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
|
|
for(i=1; i<= nlstate; i++)
|
|
gmp[j] += prlim[i][i]*p3mat[i][j][1];
|
|
}
|
|
/* end force of mortality */
|
|
|
|
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
|
|
for(j=nlstate+1; j<=(nlstate+ndeath);j++){
|
|
fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
|
|
for(i=1; i<=nlstate;i++){
|
|
fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
|
|
}
|
|
}
|
|
fprintf(ficresprobmorprev,"\n");
|
|
|
fprintf(ficresvij,"%.0f ",age );
|
fprintf(ficresvij,"%.0f ",age );
|
for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
for(j=1; j<=nlstate;j++){
|
for(j=1; j<=nlstate;j++){
|
Line 1846 void varevsij(char fileres[], double ***
|
Line 2038 void varevsij(char fileres[], double ***
|
free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
|
free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
} /* End age */
|
} /* End age */
|
|
free_vector(gpp,nlstate+1,nlstate+ndeath);
|
|
free_vector(gmp,nlstate+1,nlstate+ndeath);
|
|
free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
|
|
free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
|
|
fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");
|
|
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
|
|
fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
|
|
fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm);
|
|
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm);
|
|
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm);
|
|
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);
|
|
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);
|
|
fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);
|
|
|
free_vector(xp,1,npar);
|
free_vector(xp,1,npar);
|
free_matrix(doldm,1,nlstate,1,npar);
|
free_matrix(doldm,1,nlstate,1,nlstate);
|
free_matrix(dnewm,1,nlstate,1,nlstate);
|
free_matrix(dnewm,1,nlstate,1,npar);
|
|
free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
|
|
free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
|
|
free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
|
|
fclose(ficresprobmorprev);
|
|
fclose(ficgp);
|
|
fclose(fichtm);
|
|
|
}
|
}
|
|
|
Line 1939 void varprevlim(char fileres[], double *
|
Line 2150 void varprevlim(char fileres[], double *
|
/************ Variance of one-step probabilities ******************/
|
/************ Variance of one-step probabilities ******************/
|
void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
|
void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
|
{
|
{
|
int i, j, i1, k1, l1;
|
int i, j=0, i1, k1, l1, t, tj;
|
int k2, l2, j1, z1;
|
int k2, l2, j1, z1;
|
int k=0,l, cptcode;
|
int k=0,l, cptcode;
|
int first=1;
|
int first=1, first1;
|
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;
|
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;
|
double **dnewm,**doldm;
|
double **dnewm,**doldm;
|
double *xp;
|
double *xp;
|
Line 1962 void varprob(char optionfilefiname[], do
|
Line 2173 void varprob(char optionfilefiname[], do
|
strcat(fileresprob,fileres);
|
strcat(fileresprob,fileres);
|
if((ficresprob=fopen(fileresprob,"w"))==NULL) {
|
if((ficresprob=fopen(fileresprob,"w"))==NULL) {
|
printf("Problem with resultfile: %s\n", fileresprob);
|
printf("Problem with resultfile: %s\n", fileresprob);
|
|
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
|
}
|
}
|
strcpy(fileresprobcov,"probcov");
|
strcpy(fileresprobcov,"probcov");
|
strcat(fileresprobcov,fileres);
|
strcat(fileresprobcov,fileres);
|
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
|
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
|
printf("Problem with resultfile: %s\n", fileresprobcov);
|
printf("Problem with resultfile: %s\n", fileresprobcov);
|
|
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
|
}
|
}
|
strcpy(fileresprobcor,"probcor");
|
strcpy(fileresprobcor,"probcor");
|
strcat(fileresprobcor,fileres);
|
strcat(fileresprobcor,fileres);
|
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
|
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
|
printf("Problem with resultfile: %s\n", fileresprobcor);
|
printf("Problem with resultfile: %s\n", fileresprobcor);
|
|
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
|
}
|
}
|
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
|
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
|
|
fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
|
printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
|
printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
|
|
fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
|
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
|
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
|
|
fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
|
|
|
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
|
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
|
fprintf(ficresprob,"# Age");
|
fprintf(ficresprob,"# Age");
|
Line 2002 void varprob(char optionfilefiname[], do
|
Line 2219 void varprob(char optionfilefiname[], do
|
first=1;
|
first=1;
|
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
|
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
|
printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
|
printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
|
|
fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);
|
exit(0);
|
exit(0);
|
}
|
}
|
else{
|
else{
|
Line 2009 void varprob(char optionfilefiname[], do
|
Line 2227 void varprob(char optionfilefiname[], do
|
}
|
}
|
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);
|
exit(0);
|
exit(0);
|
}
|
}
|
else{
|
else{
|
fprintf(fichtm,"\n<H2> Computing matrix of variance-covariance of step probabilities</h2>\n");
|
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
|
fprintf(fichtm,"\n<br> We 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");
|
|
|
|
fprintf(fichtm,"\n<li><h4> Computing matrix of variance-covariance of step probabilities</h4></li>\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;
|
j=cptcoveff;
|
tj=cptcoveff;
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
if (cptcovn<1) {tj=1;ncodemax[1]=1;}
|
j1=0;
|
j1=0;
|
for(k1=1; k1<=1;k1++){
|
for(t=1; t<=tj;t++){
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
for(i1=1; i1<=ncodemax[t];i1++){
|
j1++;
|
j1++;
|
|
|
if (cptcovn>0) {
|
if (cptcovn>0) {
|
fprintf(ficresprob, "\n#********** Variable ");
|
fprintf(ficresprob, "\n#********** Variable ");
|
fprintf(ficresprobcov, "\n#********** Variable ");
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
fprintf(ficgp, "\n#********** Variable ");
|
fprintf(ficresprob, "**********\n#");
|
fprintf(fichtm, "\n<h4>********** Variable</h4>\n ");
|
fprintf(ficresprobcov, "\n#********** Variable ");
|
fprintf(ficresprobcor, "\n#********** Variable ");
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
fprintf(ficresprobcov, "**********\n#");
|
fprintf(ficresprob, "**********\n#");
|
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
fprintf(ficgp, "\n#********** Variable ");
|
fprintf(ficresprobcov, "**********\n#");
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
fprintf(ficgp, "**********\n#");
|
fprintf(ficgp, "**********\n#");
|
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
|
fprintf(ficgp, "**********\n#");
|
fprintf(fichtm, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
fprintf(fichtm, "**********\n#");
|
fprintf(fichtm, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
|
}
|
|
|
fprintf(ficresprobcor, "\n#********** Variable ");
|
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
|
|
fprintf(ficgp, "**********\n#");
|
|
}
|
|
|
for (age=bage; age<=fage; age ++){
|
for (age=bage; age<=fage; age ++){
|
cov[2]=age;
|
cov[2]=age;
|
for (k=1; k<=cptcovn;k++) {
|
for (k=1; k<=cptcovn;k++) {
|
Line 2110 void varprob(char optionfilefiname[], do
|
Line 2338 void varprob(char optionfilefiname[], do
|
/*printf("\n%d ",(int)age);
|
/*printf("\n%d ",(int)age);
|
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
|
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
|
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
|
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
|
|
fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
|
}*/
|
}*/
|
|
|
fprintf(ficresprob,"\n%d ",(int)age);
|
fprintf(ficresprob,"\n%d ",(int)age);
|
Line 2135 void varprob(char optionfilefiname[], do
|
Line 2364 void varprob(char optionfilefiname[], do
|
}
|
}
|
}/* end of loop for state */
|
}/* end of loop for state */
|
} /* end of loop for age */
|
} /* end of loop for age */
|
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
|
|
|
/* Confidence intervalle of pij */
|
|
/*
|
|
fprintf(ficgp,"\nset noparametric;unset label");
|
|
fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
|
|
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
|
|
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);
|
|
fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
|
|
fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
|
|
fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
|
|
*/
|
|
|
|
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
|
|
first1=1;
|
for (k1=1; k1<=(nlstate);k1++){
|
for (k1=1; k1<=(nlstate);k1++){
|
for (l1=1; l1<=(nlstate+ndeath);l1++){
|
for (l1=1; l1<=(nlstate+ndeath);l1++){
|
if(l1==k1) continue;
|
if(l1==k1) continue;
|
Line 2155 void varprob(char optionfilefiname[], do
|
Line 2397 void varprob(char optionfilefiname[], do
|
/* Computing eigen value of matrix of covariance */
|
/* Computing eigen value of matrix of covariance */
|
lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
|
lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
|
lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
|
lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
|
printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);
|
if(first1==1){
|
|
first1=0;
|
|
printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\nOthers in log...\n",v1,v2,cv12,lc1,lc2);
|
|
}
|
|
fprintf(ficlog,"Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);
|
/* Eigen vectors */
|
/* Eigen vectors */
|
v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
|
v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
|
v21=sqrt(1.-v11*v11);
|
v21=sqrt(1.-v11*v11);
|
Line 2166 void varprob(char optionfilefiname[], do
|
Line 2412 void varprob(char optionfilefiname[], do
|
/* mu2+ v21*lc1*cost + v21*lc2*sin(t) */
|
/* mu2+ v21*lc1*cost + v21*lc2*sin(t) */
|
if(first==1){
|
if(first==1){
|
first=0;
|
first=0;
|
fprintf(ficgp,"\nset parametric;set nolabel");
|
fprintf(ficgp,"\nset parametric;set nolabel");
|
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);
|
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);
|
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
|
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
|
fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%1d%1d-%1d%1d.png\">varpijgr%s%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1);
|
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);
|
fprintf(fichtm,"\n<br><img src=\"varpijgr%s%1d%1d-%1d%1d.png\">, ",optionfilefiname,k2,l2,k1,l1);
|
fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k2,l2,k1,l1);
|
fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\"",optionfilefiname,k2,l2,k1,l1);
|
fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);
|
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
|
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
|
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
|
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
|
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\"",\
|
/* 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\"",\
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
|
|
*/
|
|
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",\
|
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));
|
}else{
|
}else{
|
first=0;
|
first=0;
|
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
|
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
|
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
|
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
|
|
/*
|
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\"",\
|
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\"",\
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
|
|
*/
|
|
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",\
|
|
mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
|
|
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));
|
}/* if first */
|
}/* if first */
|
} /* age mod 5 */
|
} /* age mod 5 */
|
} /* end loop age */
|
} /* end loop age */
|
fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\";replot;",optionfilefiname,k2,l2,k1,l1);
|
fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1);
|
first=1;
|
first=1;
|
} /*l12 */
|
} /*l12 */
|
} /* k12 */
|
} /* k12 */
|
Line 2221 void printinghtml(char fileres[], char t
|
Line 2476 void printinghtml(char fileres[], char t
|
/*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(fichtm,"<ul><li>Result files (first order: no variance)<br>\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
|
- Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
|
- Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
|
- Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
|
- Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
|
Line 2231 void printinghtml(char fileres[], char t
|
Line 2487 void printinghtml(char fileres[], char t
|
<a href=\"e%s\">e%s</a> <br>\n</li>", \
|
<a href=\"e%s\">e%s</a> <br>\n</li>", \
|
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
|
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
|
|
|
fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n
|
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
|
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
|
|
- Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
|
|
- Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
|
|
- Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n
|
|
- Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
|
|
- Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
|
|
- Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
|
|
|
|
if(popforecast==1) fprintf(fichtm,"\n
|
|
- Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
|
|
- Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
|
|
<br>",fileres,fileres,fileres,fileres);
|
|
else
|
|
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);
|
|
fprintf(fichtm," <li>Graphs</li><p>");
|
|
|
|
m=cptcoveff;
|
m=cptcoveff;
|
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
|
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
|
Line 2272 fprintf(fichtm," <li>Graphs</li><p>");
|
Line 2513 fprintf(fichtm," <li>Graphs</li><p>");
|
fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
|
fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
|
<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
}
|
}
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
|
interval) in state (%d): v%s%d%d.png <br>
|
|
<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
|
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
|
<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
Line 2284 interval) in state (%d): v%s%d%d.png <br
|
Line 2520 interval) in state (%d): v%s%d%d.png <br
|
fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
fprintf(fichtm,"\n<br>- Total life expectancy by age and
|
health expectancies in states (1) and (2): e%s%d.png<br>
|
health expectancies in states (1) and (2): e%s%d.png<br>
|
<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
|
}
|
} /* end i1 */
|
}
|
}/* End k1 */
|
|
fprintf(fichtm,"</ul>");
|
|
|
|
|
|
fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n
|
|
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
|
|
- Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
|
|
- Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
|
|
- Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n
|
|
- Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
|
|
- Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
|
|
- Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
|
|
|
|
if(popforecast==1) fprintf(fichtm,"\n
|
|
- Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
|
|
- Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
|
|
<br>",fileres,fileres,fileres,fileres);
|
|
else
|
|
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);
|
|
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
|
|
|
|
m=cptcoveff;
|
|
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
|
|
|
|
jj1=0;
|
|
for(k1=1; k1<=m;k1++){
|
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
|
jj1++;
|
|
if (cptcovn > 0) {
|
|
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
|
|
for (cpt=1; cpt<=cptcoveff;cpt++)
|
|
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
|
|
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
|
|
}
|
|
for(cpt=1; cpt<=nlstate;cpt++) {
|
|
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
|
|
interval) in state (%d): v%s%d%d.png <br>
|
|
<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
|
|
}
|
|
} /* end i1 */
|
|
}/* End k1 */
|
|
fprintf(fichtm,"</ul>");
|
fclose(fichtm);
|
fclose(fichtm);
|
}
|
}
|
|
|
Line 2296 void printinggnuplot(char fileres[], dou
|
Line 2573 void printinggnuplot(char fileres[], dou
|
int ng;
|
int ng;
|
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
|
if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
|
printf("Problem with file %s",optionfilegnuplot);
|
printf("Problem with file %s",optionfilegnuplot);
|
|
fprintf(ficlog,"Problem with file %s",optionfilegnuplot);
|
}
|
}
|
|
|
#ifdef windows
|
#ifdef windows
|
Line 2415 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
Line 2693 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
for(k=1; k <=(nlstate+ndeath); k++){
|
for(k=1; k <=(nlstate+ndeath); k++){
|
if (k != i) {
|
if (k != i) {
|
for(j=1; j <=ncovmodel; j++){
|
for(j=1; j <=ncovmodel; j++){
|
|
|
fprintf(ficgp,"p%d=%f ",jk,p[jk]);
|
fprintf(ficgp,"p%d=%f ",jk,p[jk]);
|
jk++;
|
jk++;
|
fprintf(ficgp,"\n");
|
fprintf(ficgp,"\n");
|
Line 2469 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
Line 2746 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
|
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
|
i=i+ncovmodel;
|
i=i+ncovmodel;
|
}
|
}
|
}
|
} /* end k */
|
}
|
} /* end k2 */
|
}
|
} /* end jk */
|
}
|
} /* end ng */
|
fclose(ficgp);
|
fclose(ficgp);
|
} /* end gnuplot */
|
} /* end gnuplot */
|
|
|
Line 2520 calagedate=(anproj1+mproj1/12.+jproj1/36
|
Line 2797 calagedate=(anproj1+mproj1/12.+jproj1/36
|
strcat(fileresf,fileres);
|
strcat(fileresf,fileres);
|
if((ficresf=fopen(fileresf,"w"))==NULL) {
|
if((ficresf=fopen(fileresf,"w"))==NULL) {
|
printf("Problem with forecast resultfile: %s\n", fileresf);
|
printf("Problem with forecast resultfile: %s\n", fileresf);
|
|
fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
|
}
|
}
|
printf("Computing forecasting: result on file '%s' \n", fileresf);
|
printf("Computing forecasting: result on file '%s' \n", fileresf);
|
|
fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);
|
|
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
|
|
Line 2624 populforecast(char fileres[], double anp
|
Line 2903 populforecast(char fileres[], double anp
|
strcat(filerespop,fileres);
|
strcat(filerespop,fileres);
|
if((ficrespop=fopen(filerespop,"w"))==NULL) {
|
if((ficrespop=fopen(filerespop,"w"))==NULL) {
|
printf("Problem with forecast resultfile: %s\n", filerespop);
|
printf("Problem with forecast resultfile: %s\n", filerespop);
|
|
fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
|
}
|
}
|
printf("Computing forecasting: result on file '%s' \n", filerespop);
|
printf("Computing forecasting: result on file '%s' \n", filerespop);
|
|
fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);
|
|
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
|
|
Line 2645 populforecast(char fileres[], double anp
|
Line 2926 populforecast(char fileres[], double anp
|
if (popforecast==1) {
|
if (popforecast==1) {
|
if((ficpop=fopen(popfile,"r"))==NULL) {
|
if((ficpop=fopen(popfile,"r"))==NULL) {
|
printf("Problem with population file : %s\n",popfile);exit(0);
|
printf("Problem with population file : %s\n",popfile);exit(0);
|
|
fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0);
|
}
|
}
|
popage=ivector(0,AGESUP);
|
popage=ivector(0,AGESUP);
|
popeffectif=vector(0,AGESUP);
|
popeffectif=vector(0,AGESUP);
|
Line 2773 int main(int argc, char *argv[])
|
Line 3055 int main(int argc, char *argv[])
|
double ***p3mat;
|
double ***p3mat;
|
int *indx;
|
int *indx;
|
char line[MAXLINE], linepar[MAXLINE];
|
char line[MAXLINE], linepar[MAXLINE];
|
char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
|
char path[80],pathc[80],pathcd[80],pathtot[80],model[80];
|
int firstobs=1, lastobs=10;
|
int firstobs=1, lastobs=10;
|
int sdeb, sfin; /* Status at beginning and end */
|
int sdeb, sfin; /* Status at beginning and end */
|
int c, h , cpt,l;
|
int c, h , cpt,l;
|
Line 2834 int main(int argc, char *argv[])
|
Line 3116 int main(int argc, char *argv[])
|
|
|
/*-------- arguments in the command line --------*/
|
/*-------- arguments in the command line --------*/
|
|
|
|
/* Log file */
|
|
strcat(filelog, optionfilefiname);
|
|
strcat(filelog,".log"); /* */
|
|
if((ficlog=fopen(filelog,"w"))==NULL) {
|
|
printf("Problem with logfile %s\n",filelog);
|
|
goto end;
|
|
}
|
|
fprintf(ficlog,"Log filename:%s\n",filelog);
|
|
fprintf(ficlog,"\n%s",version);
|
|
fprintf(ficlog,"\nEnter the parameter file name: ");
|
|
fprintf(ficlog,"pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
|
|
fflush(ficlog);
|
|
|
|
/* */
|
strcpy(fileres,"r");
|
strcpy(fileres,"r");
|
strcat(fileres, optionfilefiname);
|
strcat(fileres, optionfilefiname);
|
strcat(fileres,".txt"); /* Other files have txt extension */
|
strcat(fileres,".txt"); /* Other files have txt extension */
|
Line 2842 int main(int argc, char *argv[])
|
Line 3138 int main(int argc, char *argv[])
|
|
|
if((ficpar=fopen(optionfile,"r"))==NULL) {
|
if((ficpar=fopen(optionfile,"r"))==NULL) {
|
printf("Problem with optionfile %s\n",optionfile);
|
printf("Problem with optionfile %s\n",optionfile);
|
|
fprintf(ficlog,"Problem with optionfile %s\n",optionfile);
|
goto end;
|
goto end;
|
}
|
}
|
|
|
strcpy(filereso,"o");
|
strcpy(filereso,"o");
|
strcat(filereso,fileres);
|
strcat(filereso,fileres);
|
if((ficparo=fopen(filereso,"w"))==NULL) {
|
if((ficparo=fopen(filereso,"w"))==NULL) {
|
printf("Problem with Output resultfile: %s\n", filereso);goto end;
|
printf("Problem with Output resultfile: %s\n", filereso);
|
|
fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
|
|
goto end;
|
}
|
}
|
|
|
/* Reads comments: lines beginning with '#' */
|
/* Reads comments: lines beginning with '#' */
|
Line 2860 int main(int argc, char *argv[])
|
Line 3159 int main(int argc, char *argv[])
|
}
|
}
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
|
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);
|
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);
|
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);
|
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);
|
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){
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2894 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3193 while((c=getc(ficpar))=='#' && c!= EOF){
|
for(j=1; j <=nlstate+ndeath-1; j++){
|
for(j=1; j <=nlstate+ndeath-1; j++){
|
fscanf(ficpar,"%1d%1d",&i1,&j1);
|
fscanf(ficpar,"%1d%1d",&i1,&j1);
|
fprintf(ficparo,"%1d%1d",i1,j1);
|
fprintf(ficparo,"%1d%1d",i1,j1);
|
printf("%1d%1d",i,j);
|
if(mle==1)
|
|
printf("%1d%1d",i,j);
|
|
fprintf(ficlog,"%1d%1d",i,j);
|
for(k=1; k<=ncovmodel;k++){
|
for(k=1; k<=ncovmodel;k++){
|
fscanf(ficpar," %lf",¶m[i][j][k]);
|
fscanf(ficpar," %lf",¶m[i][j][k]);
|
printf(" %lf",param[i][j][k]);
|
if(mle==1){
|
|
printf(" %lf",param[i][j][k]);
|
|
fprintf(ficlog," %lf",param[i][j][k]);
|
|
}
|
|
else
|
|
fprintf(ficlog," %lf",param[i][j][k]);
|
fprintf(ficparo," %lf",param[i][j][k]);
|
fprintf(ficparo," %lf",param[i][j][k]);
|
}
|
}
|
fscanf(ficpar,"\n");
|
fscanf(ficpar,"\n");
|
printf("\n");
|
if(mle==1)
|
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
fprintf(ficparo,"\n");
|
fprintf(ficparo,"\n");
|
}
|
}
|
|
|
Line 2949 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3257 while((c=getc(ficpar))=='#' && c!= EOF){
|
matcov=matrix(1,npar,1,npar);
|
matcov=matrix(1,npar,1,npar);
|
for(i=1; i <=npar; i++){
|
for(i=1; i <=npar; i++){
|
fscanf(ficpar,"%s",&str);
|
fscanf(ficpar,"%s",&str);
|
printf("%s",str);
|
if(mle==1)
|
|
printf("%s",str);
|
|
fprintf(ficlog,"%s",str);
|
fprintf(ficparo,"%s",str);
|
fprintf(ficparo,"%s",str);
|
for(j=1; j <=i; j++){
|
for(j=1; j <=i; j++){
|
fscanf(ficpar," %le",&matcov[i][j]);
|
fscanf(ficpar," %le",&matcov[i][j]);
|
printf(" %.5le",matcov[i][j]);
|
if(mle==1){
|
|
printf(" %.5le",matcov[i][j]);
|
|
fprintf(ficlog," %.5le",matcov[i][j]);
|
|
}
|
|
else
|
|
fprintf(ficlog," %.5le",matcov[i][j]);
|
fprintf(ficparo," %.5le",matcov[i][j]);
|
fprintf(ficparo," %.5le",matcov[i][j]);
|
}
|
}
|
fscanf(ficpar,"\n");
|
fscanf(ficpar,"\n");
|
printf("\n");
|
if(mle==1)
|
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
fprintf(ficparo,"\n");
|
fprintf(ficparo,"\n");
|
}
|
}
|
for(i=1; i <=npar; i++)
|
for(i=1; i <=npar; i++)
|
for(j=i+1;j<=npar;j++)
|
for(j=i+1;j<=npar;j++)
|
matcov[i][j]=matcov[j][i];
|
matcov[i][j]=matcov[j][i];
|
|
|
printf("\n");
|
if(mle==1)
|
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
|
|
|
|
/*-------- Rewriting paramater file ----------*/
|
/*-------- Rewriting paramater file ----------*/
|
Line 2974 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3293 while((c=getc(ficpar))=='#' && c!= EOF){
|
strcat(rfileres,optionfilext); /* Other files have txt extension */
|
strcat(rfileres,optionfilext); /* Other files have txt extension */
|
if((ficres =fopen(rfileres,"w"))==NULL) {
|
if((ficres =fopen(rfileres,"w"))==NULL) {
|
printf("Problem writing new parameter file: %s\n", fileres);goto end;
|
printf("Problem writing new parameter file: %s\n", fileres);goto end;
|
|
fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
|
}
|
}
|
fprintf(ficres,"#%s\n",version);
|
fprintf(ficres,"#%s\n",version);
|
|
|
/*-------- data file ----------*/
|
/*-------- data file ----------*/
|
if((fic=fopen(datafile,"r"))==NULL) {
|
if((fic=fopen(datafile,"r"))==NULL) {
|
printf("Problem with datafile: %s\n", datafile);goto end;
|
printf("Problem with datafile: %s\n", datafile);goto end;
|
|
fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;
|
}
|
}
|
|
|
n= lastobs;
|
n= lastobs;
|
Line 3045 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3366 while((c=getc(ficpar))=='#' && c!= EOF){
|
|
|
|
|
/* Calculation of the number of parameter from char model*/
|
/* Calculation of the number of parameter from char model*/
|
Tvar=ivector(1,15);
|
Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
|
Tprod=ivector(1,15);
|
Tprod=ivector(1,15);
|
Tvaraff=ivector(1,15);
|
Tvaraff=ivector(1,15);
|
Tvard=imatrix(1,15,1,2);
|
Tvard=imatrix(1,15,1,2);
|
Line 3061 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3382 while((c=getc(ficpar))=='#' && c!= EOF){
|
strcpy(modelsav,model);
|
strcpy(modelsav,model);
|
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
|
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
|
printf("Error. Non available option model=%s ",model);
|
printf("Error. Non available option model=%s ",model);
|
|
fprintf(ficlog,"Error. Non available option model=%s ",model);
|
goto end;
|
goto end;
|
}
|
}
|
|
|
for(i=(j+1); i>=1;i--){
|
for(i=(j+1); i>=1;i--){
|
cutv(stra,strb,modelsav,'+');
|
cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */
|
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);
|
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */
|
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
|
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
|
/*scanf("%d",i);*/
|
/*scanf("%d",i);*/
|
if (strchr(strb,'*')) {
|
if (strchr(strb,'*')) { /* Model includes a product */
|
cutv(strd,strc,strb,'*');
|
cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *age)*/
|
if (strcmp(strc,"age")==0) {
|
if (strcmp(strc,"age")==0) { /* Vn*age */
|
cptcovprod--;
|
cptcovprod--;
|
cutv(strb,stre,strd,'V');
|
cutv(strb,stre,strd,'V');
|
Tvar[i]=atoi(stre);
|
Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/
|
cptcovage++;
|
cptcovage++;
|
Tage[cptcovage]=i;
|
Tage[cptcovage]=i;
|
/*printf("stre=%s ", stre);*/
|
/*printf("stre=%s ", stre);*/
|
}
|
}
|
else if (strcmp(strd,"age")==0) {
|
else if (strcmp(strd,"age")==0) { /* or age*Vn */
|
cptcovprod--;
|
cptcovprod--;
|
cutv(strb,stre,strc,'V');
|
cutv(strb,stre,strc,'V');
|
Tvar[i]=atoi(stre);
|
Tvar[i]=atoi(stre);
|
cptcovage++;
|
cptcovage++;
|
Tage[cptcovage]=i;
|
Tage[cptcovage]=i;
|
}
|
}
|
else {
|
else { /* Age is not in the model */
|
cutv(strb,stre,strc,'V');
|
cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/
|
Tvar[i]=ncovcol+k1;
|
Tvar[i]=ncovcol+k1;
|
cutv(strb,strc,strd,'V');
|
cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
|
Tprod[k1]=i;
|
Tprod[k1]=i;
|
Tvard[k1][1]=atoi(strc);
|
Tvard[k1][1]=atoi(strc); /* m*/
|
Tvard[k1][2]=atoi(stre);
|
Tvard[k1][2]=atoi(stre); /* n */
|
Tvar[cptcovn+k2]=Tvard[k1][1];
|
Tvar[cptcovn+k2]=Tvard[k1][1];
|
Tvar[cptcovn+k2+1]=Tvard[k1][2];
|
Tvar[cptcovn+k2+1]=Tvard[k1][2];
|
for (k=1; k<=lastobs;k++)
|
for (k=1; k<=lastobs;k++)
|
Line 3101 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3423 while((c=getc(ficpar))=='#' && c!= EOF){
|
k2=k2+2;
|
k2=k2+2;
|
}
|
}
|
}
|
}
|
else {
|
else { /* no more sum */
|
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
|
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
|
/* scanf("%d",i);*/
|
/* scanf("%d",i);*/
|
cutv(strd,strc,strb,'V');
|
cutv(strd,strc,strb,'V');
|
Line 3110 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3432 while((c=getc(ficpar))=='#' && c!= EOF){
|
strcpy(modelsav,stra);
|
strcpy(modelsav,stra);
|
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
|
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
|
scanf("%d",i);*/
|
scanf("%d",i);*/
|
}
|
} /* end of loop + */
|
}
|
} /* end model */
|
|
|
/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
|
/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
|
printf("cptcovprod=%d ", cptcovprod);
|
printf("cptcovprod=%d ", cptcovprod);
|
|
fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
|
scanf("%d ",i);*/
|
scanf("%d ",i);*/
|
fclose(fic);
|
fclose(fic);
|
|
|
Line 3147 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3470 while((c=getc(ficpar))=='#' && c!= EOF){
|
else {
|
else {
|
if (andc[i]!=9999){
|
if (andc[i]!=9999){
|
printf("Warning negative age at death: %d line:%d\n",num[i],i);
|
printf("Warning negative age at death: %d line:%d\n",num[i],i);
|
|
fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
|
agev[m][i]=-1;
|
agev[m][i]=-1;
|
}
|
}
|
}
|
}
|
Line 3179 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3503 while((c=getc(ficpar))=='#' && c!= EOF){
|
for (i=1; i<=imx; i++) {
|
for (i=1; i<=imx; i++) {
|
for(m=1; (m<= maxwav); m++){
|
for(m=1; (m<= maxwav); m++){
|
if (s[m][i] > (nlstate+ndeath)) {
|
if (s[m][i] > (nlstate+ndeath)) {
|
printf("Error: Wrong value in nlstate or ndeath\n");
|
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);
|
|
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);
|
goto end;
|
goto end;
|
}
|
}
|
}
|
}
|
}
|
}
|
|
|
printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
|
printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
|
|
fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
|
|
|
free_vector(severity,1,maxwav);
|
free_vector(severity,1,maxwav);
|
free_imatrix(outcome,1,maxwav+1,1,n);
|
free_imatrix(outcome,1,maxwav+1,1,n);
|
Line 3261 printf("Total number of individuals= %d,
|
Line 3587 printf("Total number of individuals= %d,
|
jk=1;
|
jk=1;
|
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
|
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
|
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
|
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
|
|
fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
|
for(i=1,jk=1; i <=nlstate; i++){
|
for(i=1,jk=1; i <=nlstate; i++){
|
for(k=1; k <=(nlstate+ndeath); k++){
|
for(k=1; k <=(nlstate+ndeath); k++){
|
if (k != i)
|
if (k != i)
|
{
|
{
|
printf("%d%d ",i,k);
|
printf("%d%d ",i,k);
|
|
fprintf(ficlog,"%d%d ",i,k);
|
fprintf(ficres,"%1d%1d ",i,k);
|
fprintf(ficres,"%1d%1d ",i,k);
|
for(j=1; j <=ncovmodel; j++){
|
for(j=1; j <=ncovmodel; j++){
|
printf("%f ",p[jk]);
|
printf("%f ",p[jk]);
|
|
fprintf(ficlog,"%f ",p[jk]);
|
fprintf(ficres,"%f ",p[jk]);
|
fprintf(ficres,"%f ",p[jk]);
|
jk++;
|
jk++;
|
}
|
}
|
printf("\n");
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
fprintf(ficres,"\n");
|
fprintf(ficres,"\n");
|
}
|
}
|
}
|
}
|
}
|
}
|
if(mle==1){
|
if(mle==1){
|
/* 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);
|
}
|
}
|
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
|
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
|
printf("# Scales (for hessian or gradient estimation)\n");
|
printf("# Scales (for hessian or gradient estimation)\n");
|
for(i=1,jk=1; i <=nlstate; i++){
|
fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
|
for(j=1; j <=nlstate+ndeath; j++){
|
for(i=1,jk=1; i <=nlstate; i++){
|
if (j!=i) {
|
for(j=1; j <=nlstate+ndeath; j++){
|
fprintf(ficres,"%1d%1d",i,j);
|
if (j!=i) {
|
printf("%1d%1d",i,j);
|
fprintf(ficres,"%1d%1d",i,j);
|
for(k=1; k<=ncovmodel;k++){
|
printf("%1d%1d",i,j);
|
printf(" %.5e",delti[jk]);
|
fprintf(ficlog,"%1d%1d",i,j);
|
fprintf(ficres," %.5e",delti[jk]);
|
for(k=1; k<=ncovmodel;k++){
|
jk++;
|
printf(" %.5e",delti[jk]);
|
}
|
fprintf(ficlog," %.5e",delti[jk]);
|
printf("\n");
|
fprintf(ficres," %.5e",delti[jk]);
|
fprintf(ficres,"\n");
|
jk++;
|
}
|
}
|
}
|
printf("\n");
|
|
fprintf(ficlog,"\n");
|
|
fprintf(ficres,"\n");
|
|
}
|
}
|
}
|
|
}
|
k=1;
|
|
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");
|
k=1;
|
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");
|
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");
|
for(i=1;i<=npar;i++){
|
if(mle==1)
|
/* if (k>nlstate) k=1;
|
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");
|
i1=(i-1)/(ncovmodel*nlstate)+1;
|
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");
|
fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
|
for(i=1;i<=npar;i++){
|
printf("%s%d%d",alph[k],i1,tab[i]);*/
|
/* if (k>nlstate) k=1;
|
fprintf(ficres,"%3d",i);
|
i1=(i-1)/(ncovmodel*nlstate)+1;
|
printf("%3d",i);
|
fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
|
for(j=1; j<=i;j++){
|
printf("%s%d%d",alph[k],i1,tab[i]);*/
|
fprintf(ficres," %.5e",matcov[i][j]);
|
fprintf(ficres,"%3d",i);
|
printf(" %.5e",matcov[i][j]);
|
if(mle==1)
|
}
|
printf("%3d",i);
|
fprintf(ficres,"\n");
|
fprintf(ficlog,"%3d",i);
|
printf("\n");
|
for(j=1; j<=i;j++){
|
k++;
|
fprintf(ficres," %.5e",matcov[i][j]);
|
}
|
if(mle==1)
|
|
printf(" %.5e",matcov[i][j]);
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
fprintf(ficlog," %.5e",matcov[i][j]);
|
ungetc(c,ficpar);
|
}
|
fgets(line, MAXLINE, ficpar);
|
fprintf(ficres,"\n");
|
puts(line);
|
if(mle==1)
|
fputs(line,ficparo);
|
printf("\n");
|
}
|
fprintf(ficlog,"\n");
|
ungetc(c,ficpar);
|
k++;
|
estepm=0;
|
}
|
fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
|
|
if (estepm==0 || estepm < stepm) estepm=stepm;
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
if (fage <= 2) {
|
ungetc(c,ficpar);
|
bage = ageminpar;
|
fgets(line, MAXLINE, ficpar);
|
fage = agemaxpar;
|
puts(line);
|
}
|
fputs(line,ficparo);
|
|
}
|
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
|
ungetc(c,ficpar);
|
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
estepm=0;
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
|
|
if (estepm==0 || estepm < stepm) estepm=stepm;
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
if (fage <= 2) {
|
ungetc(c,ficpar);
|
bage = ageminpar;
|
fgets(line, MAXLINE, ficpar);
|
fage = agemaxpar;
|
puts(line);
|
}
|
fputs(line,ficparo);
|
|
}
|
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
|
ungetc(c,ficpar);
|
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
|
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
|
ungetc(c,ficpar);
|
|
fgets(line, MAXLINE, ficpar);
|
|
puts(line);
|
|
fputs(line,ficparo);
|
|
}
|
|
ungetc(c,ficpar);
|
|
|
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
|
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
|
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
fgets(line, MAXLINE, ficpar);
|
fgets(line, MAXLINE, ficpar);
|
puts(line);
|
puts(line);
|
fputs(line,ficparo);
|
fputs(line,ficparo);
|
}
|
}
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
|
|
|
dateprev1=anprev1+mprev1/12.+jprev1/365.;
|
dateprev1=anprev1+mprev1/12.+jprev1/365.;
|
Line 3415 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 3757 Title=%s <br>Datafile=%s Firstpass=%d La
|
Total number of observations=%d <br>\n
|
Total number of observations=%d <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
|
<hr size=\"2\" color=\"#EC5E5E\">
|
<hr size=\"2\" color=\"#EC5E5E\">
|
<ul><li>Parameter files<br>\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
|
- Gnuplot file name: <a href=\"%s\">%s</a><br></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);
|
- Log file of the run: <a href=\"%s\">%s</a><br>\n
|
|
- 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);
|
fclose(fichtm);
|
fclose(fichtm);
|
|
|
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
|
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
|
Line 3441 Interval (in months) between two waves:
|
Line 3784 Interval (in months) between two waves:
|
strcat(filerespl,fileres);
|
strcat(filerespl,fileres);
|
if((ficrespl=fopen(filerespl,"w"))==NULL) {
|
if((ficrespl=fopen(filerespl,"w"))==NULL) {
|
printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;
|
printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;
|
|
fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;
|
}
|
}
|
printf("Computing prevalence limit: result on file '%s' \n", filerespl);
|
printf("Computing prevalence limit: result on file '%s' \n", filerespl);
|
|
fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);
|
fprintf(ficrespl,"#Prevalence limit\n");
|
fprintf(ficrespl,"#Prevalence limit\n");
|
fprintf(ficrespl,"#Age ");
|
fprintf(ficrespl,"#Age ");
|
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
|
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
|
Line 3466 Interval (in months) between two waves:
|
Line 3811 Interval (in months) between two waves:
|
k=k+1;
|
k=k+1;
|
/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
|
/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
|
fprintf(ficrespl,"\n#******");
|
fprintf(ficrespl,"\n#******");
|
for(j=1;j<=cptcoveff;j++)
|
printf("\n#******");
|
|
fprintf(ficlog,"\n#******");
|
|
for(j=1;j<=cptcoveff;j++) {
|
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
|
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
|
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
|
}
|
fprintf(ficrespl,"******\n");
|
fprintf(ficrespl,"******\n");
|
|
printf("******\n");
|
|
fprintf(ficlog,"******\n");
|
|
|
for (age=agebase; age<=agelim; age++){
|
for (age=agebase; age<=agelim; age++){
|
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
|
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
|
Line 3486 Interval (in months) between two waves:
|
Line 3838 Interval (in months) between two waves:
|
strcpy(filerespij,"pij"); strcat(filerespij,fileres);
|
strcpy(filerespij,"pij"); strcat(filerespij,fileres);
|
if((ficrespij=fopen(filerespij,"w"))==NULL) {
|
if((ficrespij=fopen(filerespij,"w"))==NULL) {
|
printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
|
printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
|
|
fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;
|
}
|
}
|
printf("Computing pij: result on file '%s' \n", filerespij);
|
printf("Computing pij: result on file '%s' \n", filerespij);
|
|
fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
|
|
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
/*if (stepm<=24) stepsize=2;*/
|
/*if (stepm<=24) stepsize=2;*/
|
Line 3495 Interval (in months) between two waves:
|
Line 3849 Interval (in months) between two waves:
|
agelim=AGESUP;
|
agelim=AGESUP;
|
hstepm=stepsize*YEARM; /* Every year of age */
|
hstepm=stepsize*YEARM; /* Every year of age */
|
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
|
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
|
|
|
|
/* hstepm=1; aff par mois*/
|
|
|
k=0;
|
k=0;
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
Line 3508 Interval (in months) between two waves:
|
Line 3864 Interval (in months) between two waves:
|
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
|
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
|
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
|
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
|
|
|
|
/* nhstepm=nhstepm*YEARM; aff par mois*/
|
|
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
oldm=oldms;savm=savms;
|
oldm=oldms;savm=savms;
|
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
|
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
|
Line 3517 Interval (in months) between two waves:
|
Line 3876 Interval (in months) between two waves:
|
fprintf(ficrespij," %1d-%1d",i,j);
|
fprintf(ficrespij," %1d-%1d",i,j);
|
fprintf(ficrespij,"\n");
|
fprintf(ficrespij,"\n");
|
for (h=0; h<=nhstepm; h++){
|
for (h=0; h<=nhstepm; h++){
|
fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
|
fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
|
for(i=1; i<=nlstate;i++)
|
for(i=1; i<=nlstate;i++)
|
for(j=1; j<=nlstate+ndeath;j++)
|
for(j=1; j<=nlstate+ndeath;j++)
|
fprintf(ficrespij," %.5f", p3mat[i][j][h]);
|
fprintf(ficrespij," %.5f", p3mat[i][j][h]);
|
Line 3542 Interval (in months) between two waves:
|
Line 3901 Interval (in months) between two waves:
|
else{
|
else{
|
erreur=108;
|
erreur=108;
|
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);
|
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);
|
|
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);
|
}
|
}
|
|
|
|
|
Line 3551 Interval (in months) between two waves:
|
Line 3911 Interval (in months) between two waves:
|
strcat(filerest,fileres);
|
strcat(filerest,fileres);
|
if((ficrest=fopen(filerest,"w"))==NULL) {
|
if((ficrest=fopen(filerest,"w"))==NULL) {
|
printf("Problem with total LE resultfile: %s\n", filerest);goto end;
|
printf("Problem with total LE resultfile: %s\n", filerest);goto end;
|
|
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
|
}
|
}
|
printf("Computing Total LEs with variances: file '%s' \n", filerest);
|
printf("Computing Total LEs with variances: file '%s' \n", filerest);
|
|
fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest);
|
|
|
|
|
strcpy(filerese,"e");
|
strcpy(filerese,"e");
|
strcat(filerese,fileres);
|
strcat(filerese,fileres);
|
if((ficreseij=fopen(filerese,"w"))==NULL) {
|
if((ficreseij=fopen(filerese,"w"))==NULL) {
|
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
|
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
|
|
fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
|
}
|
}
|
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
|
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
|
|
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
|
|
|
strcpy(fileresv,"v");
|
strcpy(fileresv,"v");
|
strcat(fileresv,fileres);
|
strcat(fileresv,fileres);
|
if((ficresvij=fopen(fileresv,"w"))==NULL) {
|
if((ficresvij=fopen(fileresv,"w"))==NULL) {
|
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
|
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
|
|
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
|
}
|
}
|
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
|
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
|
|
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
|
calagedate=-1;
|
calagedate=-1;
|
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
k=0;
|
k=0;
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
|
k=k+1;
|
k=k+1;
|
fprintf(ficrest,"\n#****** ");
|
fprintf(ficrest,"\n#****** ");
|
for(j=1;j<=cptcoveff;j++)
|
for(j=1;j<=cptcoveff;j++)
|
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
Line 3596 prevalence(ageminpar, agemax, s, agev, n
|
Line 3962 prevalence(ageminpar, agemax, s, agev, n
|
|
|
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
|
oldm=oldms;savm=savms;
|
oldm=oldms;savm=savms;
|
varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
|
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0);
|
|
if(popbased==1){
|
|
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased);
|
|
}
|
|
|
|
|
fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
|
fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
|
Line 3690 free_matrix(mint,1,maxwav,1,n);
|
Line 4058 free_matrix(mint,1,maxwav,1,n);
|
fclose(ficgp);
|
fclose(ficgp);
|
|
|
|
|
if(erreur >0)
|
if(erreur >0){
|
printf("End of Imach with error or warning %d\n",erreur);
|
printf("End of Imach with error or warning %d\n",erreur);
|
else printf("End of Imach\n");
|
fprintf(ficlog,"End of Imach with error or warning %d\n",erreur);
|
|
}else{
|
|
printf("End of Imach\n");
|
|
fprintf(ficlog,"End of Imach\n");
|
|
}
|
|
printf("See log file on %s\n",filelog);
|
|
fclose(ficlog);
|
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */
|
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */
|
|
|
/* 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);*/
|