--- imach/src/imach-0.98r.c 2016/02/19 09:46:04 1.1 +++ imach/src/imach-0.98r.c 2019/05/21 15:55:06 1.2 @@ -1,6 +1,9 @@ -/* $Id: imach-0.98r.c,v 1.1 2016/02/19 09:46:04 brouard Exp $ +/* $Id: imach-0.98r.c,v 1.2 2019/05/21 15:55:06 brouard Exp $ $State: Exp $ $Log: imach-0.98r.c,v $ + Revision 1.2 2019/05/21 15:55:06 brouard + Summary: RPM and SRPMS + Revision 1.1 2016/02/19 09:46:04 brouard Summary: Kind of 0.98r? series, starting with r7 @@ -747,7 +750,9 @@ Back prevalence and projections: #define DEBUGHESSIJ /* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ #define POWELL /* Instead of NLOPT */ -#define POWELLF1F3 /* Skip test */ +/* #ifndef POWELLNOF3INFF1TEST */ +/* #define POWELLNOF3INFF1TEST /\* Skip test *\/ */ +/* #endif */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ @@ -837,12 +842,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach-0.98r.c,v 1.1 2016/02/19 09:46:04 brouard Exp $ */ +/* $Id: imach-0.98r.c,v 1.2 2019/05/21 15:55:06 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" -char version[]=__IMACH_VERSION__; -char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.1 $ $Date: 2016/02/19 09:46:04 $"; +char version[]=__IMACH_VERSION_98R__; +char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; +char fullversion[]="$Revision: 1.2 $ $Date: 2019/05/21 15:55:06 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1545,12 +1550,12 @@ double brent(double ax, double bx, doubl etemp=e; e=d; if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) - d=CGOLD*(e=(x >= xm ? a-x : b-x)); + d=CGOLD*(e=(x >= xm ? a-x : b-x)); else { - d=p/q; - u=x+d; - if (u-a < tol2 || b-u < tol2) - d=SIGN(tol1,xm-x); + d=p/q; + u=x+d; + if (u-a < tol2 || b-u < tol2) + d=SIGN(tol1,xm-x); } } else { d=CGOLD*(e=(x >= xm ? a-x : b-x)); @@ -1564,13 +1569,13 @@ double brent(double ax, double bx, doubl } else { if (u < x) a=u; else b=u; if (fu <= fw || w == x) { - v=w; - w=u; - fv=fw; - fw=fu; + v=w; + w=u; + fv=fw; + fw=fu; } else if (fu <= fv || v == x || v == w) { - v=u; - fv=fu; + v=u; + fv=fu; } } } @@ -1611,24 +1616,24 @@ values at the three points, fa, fb , and *cx=(*bx)+GOLD*(*bx-*ax); *fc=(*func)(*cx); #ifdef DEBUG - printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); - fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); + printf("mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc); + fprintf(ficlog,"mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc); #endif - while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */ + while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc. If fc=inf it exits*/ r=(*bx-*ax)*(*fb-*fc); - q=(*bx-*cx)*(*fb-*fa); + q=(*bx-*cx)*(*fb-*fa); /* What if fa=inf */ u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscissa of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */ ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */ - if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is between b and c */ + if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is between b and c */ fu=(*func)(u); #ifdef DEBUG /* f(x)=A(x-u)**2+f(u) */ double A, fparabu; A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u); fparabu= *fa - A*(*ax-u)*(*ax-u); - printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu); - fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu); + printf("\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r); + fprintf(ficlog,"\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r); /* And thus,it can be that fu > *fc even if fparabu < *fc */ /* mnbrak (*ax=7.666299858533, *fa=299039.693133272231), (*bx=8.595447774979, *fb=298976.598289369489), (*cx=10.098840694817, *fc=298946.631474258087), (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */ @@ -1661,9 +1666,12 @@ values at the three points, fa, fb , and /* fu = *fc; */ /* *fc =dum; */ /* } */ -#ifdef DEBUG - printf("mnbrak34 fu < or >= fc \n"); - fprintf(ficlog, "mnbrak34 fu < fc\n"); +#ifdef DEBUGMNBRAK + double A, fparabu; + A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u); + fparabu= *fa - A*(*ax-u)*(*ax-u); + printf("\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r); + fprintf(ficlog,"\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r); #endif dum=u; /* Shifting c and u */ u = *cx; @@ -1674,38 +1682,45 @@ values at the three points, fa, fb , and #endif } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */ #ifdef DEBUG - printf("mnbrak2 u after c but before ulim\n"); - fprintf(ficlog, "mnbrak2 u after c but before ulim\n"); + printf("\nmnbrak2 u=%lf after c=%lf but before ulim\n",u,*cx); + fprintf(ficlog,"\nmnbrak2 u=%lf after c=%lf but before ulim\n",u,*cx); #endif fu=(*func)(u); if (fu < *fc) { #ifdef DEBUG - printf("mnbrak2 u after c but before ulim AND fu < fc\n"); - fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu = 0.0) { /* u outside ulim (verifying that ulim is beyond c) */ #ifdef DEBUG - printf("mnbrak2 u outside ulim (verifying that ulim is beyond c)\n"); - fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n"); + printf("\nmnbrak2 u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx); + fprintf(ficlog,"\nmnbrak2 u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx); #endif u=ulim; fu=(*func)(u); } else { /* u could be left to b (if r > q parabola has a maximum) */ #ifdef DEBUG - printf("mnbrak2 u could be left to b (if r > q parabola has a maximum)\n"); - fprintf(ficlog, "mnbrak2 u could be left to b (if r > q parabola has a maximum)\n"); + printf("\nmnbrak2 u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q); + fprintf(ficlog,"\nmnbrak2 u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q); #endif u=(*cx)+GOLD*(*cx-*bx); fu=(*func)(u); +#ifdef DEBUG + printf("\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx); + fprintf(ficlog,"\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx); +#endif } /* end tests */ SHFT(*ax,*bx,*cx,u) SHFT(*fa,*fb,*fc,fu) #ifdef DEBUG - printf("mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu); - fprintf(ficlog, "mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu); + printf("\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc); + fprintf(ficlog, "\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc); #endif } /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */ } @@ -1720,7 +1735,11 @@ int ncom; double *pcom,*xicom; double (*nrfunc)(double []); +#ifdef LINMINORIGINAL void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [])) +#else +void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []), int *flat) +#endif { double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin); @@ -1764,28 +1783,40 @@ void linmin(double p[], double xi[], int #ifdef LINMINORIGINAL #else if (fx != fx){ - xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ - printf("|"); - fprintf(ficlog,"|"); + xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ + printf("|"); + fprintf(ficlog,"|"); #ifdef DEBUGLINMIN - printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); + printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); #endif } - }while(fx != fx); + }while(fx != fx && xxs > 1.e-5); #endif #ifdef DEBUGLINMIN printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); #endif +#ifdef LINMINORIGINAL +#else + if(fb == fx){ /* Flat function in the direction */ + xmin=xx; + *flat=1; + }else{ + *flat=0; +#endif *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/ /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */ /* fmin = f(p[j] + xmin * xi[j]) */ /* P+lambda n in that direction (lambdamin), with TOL between abscisses */ /* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */ #ifdef DEBUG - printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); - fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); + printf("retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin); + fprintf(ficlog,"retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin); +#endif +#ifdef LINMINORIGINAL +#else + } #endif #ifdef DEBUGLINMIN printf("linmin end "); @@ -1835,17 +1866,33 @@ such that failure to decrease by more th output, p is set to the best point found, xi is the then-current direction set, fret is the returned function value at p , and iter is the number of iterations taken. The routine linmin is used. */ +#ifdef LINMINORIGINAL +#else + int *flatdir; /* Function is vanishing in that direction */ + int flat=0; /* Function is vanishing in that direction */ +#endif void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, double (*func)(double [])) { - void linmin(double p[], double xi[], int n, double *fret, +#ifdef LINMINORIGINAL + void linmin(double p[], double xi[], int n, double *fret, double (*func)(double [])); +#else + void linmin(double p[], double xi[], int n, double *fret, + double (*func)(double []),int *flat); +#endif int i,ibig,j; double del,t,*pt,*ptt,*xit; double directest; double fp,fptt; double *xits; int niterf, itmp; +#ifdef LINMINORIGINAL +#else + + flatdir=ivector(1,n); + for (j=1;j<=n;j++) flatdir[j]=0; +#endif pt=vector(1,n); ptt=vector(1,n); @@ -1879,18 +1926,18 @@ void powell(double p[], double **xi, int rforecast_time=rcurr_time; itmp = strlen(strcurr); if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ - strcurr[itmp-1]='\0'; + strcurr[itmp-1]='\0'; printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); for(niterf=10;niterf<=30;niterf+=10){ - rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); - forecast_time = *localtime(&rforecast_time); - strcpy(strfor,asctime(&forecast_time)); - itmp = strlen(strfor); - if(strfor[itmp-1]=='\n') - strfor[itmp-1]='\0'; - printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); - fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); + rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); + forecast_time = *localtime(&rforecast_time); + strcpy(strfor,asctime(&forecast_time)); + itmp = strlen(strfor); + if(strfor[itmp-1]=='\n') + strfor[itmp-1]='\0'; + printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); + fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); } } for (i=1;i<=n;i++) { /* For each direction i */ @@ -1902,27 +1949,32 @@ void powell(double p[], double **xi, int #endif printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); +#ifdef LINMINORIGINAL linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ - /* Outputs are fret(new point p) p is updated and xit rescaled */ +#else + linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ + flatdir[i]=flat; /* Function is vanishing in that direction i */ +#endif + /* Outputs are fret(new point p) p is updated and xit rescaled */ if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */ - /* because that direction will be replaced unless the gain del is small */ - /* in comparison with the 'probable' gain, mu^2, with the last average direction. */ - /* Unless the n directions are conjugate some gain in the determinant may be obtained */ - /* with the new direction. */ - del=fabs(fptt-(*fret)); - ibig=i; + /* because that direction will be replaced unless the gain del is small */ + /* in comparison with the 'probable' gain, mu^2, with the last average direction. */ + /* Unless the n directions are conjugate some gain in the determinant may be obtained */ + /* with the new direction. */ + del=fabs(fptt-(*fret)); + ibig=i; } #ifdef DEBUG printf("%d %.12e",i,(*fret)); fprintf(ficlog,"%d %.12e",i,(*fret)); for (j=1;j<=n;j++) { - xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); - printf(" x(%d)=%.12e",j,xit[j]); - fprintf(ficlog," x(%d)=%.12e",j,xit[j]); + xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); + printf(" x(%d)=%.12e",j,xit[j]); + fprintf(ficlog," x(%d)=%.12e",j,xit[j]); } for(j=1;j<=n;j++) { - printf(" p(%d)=%.12e",j,p[j]); - fprintf(ficlog," p(%d)=%.12e",j,p[j]); + printf(" p(%d)=%lf ",j,p[j]); + fprintf(ficlog," p(%d)=%lf ",j,p[j]); } printf("\n"); fprintf(ficlog,"\n"); @@ -1931,6 +1983,13 @@ void powell(double p[], double **xi, int /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ /* New value of last point Pn is not computed, P(n-1) */ + for(j=1;j<=n;j++) { + printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); + fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); + } + printf("\n"); + fprintf(ficlog,"\n"); + if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ @@ -1939,7 +1998,7 @@ void powell(double p[], double **xi, int /* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */ /* By using V1+V2+V3, the gain should be 7.82, compared with basic 1+age. */ /* By adding 10 parameters more the gain should be 18.31 */ - + /* Starting the program with initial values given by a former maximization will simply change */ /* the scales of the directions and the directions, because the are reset to canonical directions */ /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */ @@ -1967,7 +2026,10 @@ void powell(double p[], double **xi, int } #endif - +#ifdef LINMINORIGINAL +#else + free_ivector(flatdir,1,n); +#endif free_vector(xit,1,n); free_vector(xits,1,n); free_vector(ptt,1,n); @@ -1981,7 +2043,7 @@ void powell(double p[], double **xi, int pt[j]=p[j]; } fptt=(*func)(ptt); /* f_3 */ -#ifdef POWELLF1F3 +#ifdef POWELLNOF3INFF1TEST /* skips test F3 0 */ + /* mu² and del² are equal when f3=f1 */ + /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */ + /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */ + /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */ + /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */ #ifdef NRCORIGINAL t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/ #else @@ -2013,54 +2083,71 @@ void powell(double p[], double **xi, int if (t < 0.0) { /* Then we use it for new direction */ #else if (directest*t < 0.0) { /* Contradiction between both tests */ - printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); + printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); - fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); + fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); } if (directest < 0.0) { /* Then we use it for new direction */ #endif #ifdef DEBUGLINMIN - printf("Before linmin in direction P%d-P0\n",n); - for (j=1;j<=n;j++) { - printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0){ - printf("\n"); - fprintf(ficlog,"\n"); - } - } + printf("Before linmin in direction P%d-P0\n",n); + for (j=1;j<=n;j++) { + printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } #endif - linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ -#ifdef DEBUGLINMIN - for (j=1;j<=n;j++) { - printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0){ - printf("\n"); - fprintf(ficlog,"\n"); - } - } +#ifdef LINMINORIGINAL + linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ +#else + linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ + flatdir[i]=flat; /* Function is vanishing in that direction i */ #endif - for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ - xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ - } - printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); +#ifdef DEBUGLINMIN + for (j=1;j<=n;j++) { + printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } +#endif + for (j=1;j<=n;j++) { + xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ + xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ + } +#ifdef LINMINORIGINAL +#else + printf("Flat directions\n"); + fprintf(ficlog,"Flat directions\n"); + for (j=1;j<=n;j++) { + printf("flatdir[%d]=%d ",j,flatdir[j]); + fprintf(ficlog,"flatdir[%d]=%d ",j,flatdir[j]); + } + printf("\n"); + fprintf(ficlog,"\n"); +#endif + printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + #ifdef DEBUG - printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - 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]); - fprintf(ficlog," %.12e",xit[j]); - } - printf("\n"); - fprintf(ficlog,"\n"); + printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); + 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(" %lf",xit[j]); + fprintf(ficlog," %lf",xit[j]); + } + printf("\n"); + fprintf(ficlog,"\n"); #endif } /* end of t or directest negative */ -#ifdef POWELLF1F3 +#ifdef POWELLNOF3INFF1TEST #else } /* end if (fptt < fp) */ #endif @@ -2337,65 +2424,65 @@ double **pmij(double **ps, double *cov, /*double t34;*/ int i,j, nc, ii, jj; - for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ - } - ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ - } - } + for(i=1; i<= nlstate; i++){ + for(j=1; ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ + } + ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ + } + } - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ - ps[i][i]=1./(s1+1.); - /* Computing other pijs */ - for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ + ps[i][i]=1./(s1+1.); + /* Computing other pijs */ + for(j=1; j(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ - } - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ - ipmx +=1; - sw += weight[i]; - ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - /* if (lli < log(mytinydouble)){ */ - /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ - /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ - /* } */ - } /* end of wave */ - } /* end of individual */ - } else if(mle==2){ - for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; - for(mi=1; mi<= wav[i]-1; mi++){ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - for(d=0; d<=dh[mi][i]; d++){ - newm=savm; - agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; - cov[2]=agexact; - if(nagesqr==1) - cov[3]= agexact*agexact; - for (kk=1; kk<=cptcovage;kk++) { - cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; - } - out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, - 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); - savm=oldm; - oldm=newm; - } /* end mult */ + } else if ( s2==-1 ) { /* alive */ + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + /*survp += out[s1][j]; */ + lli= log(survp); + } + else if (s2==-4) { + for (j=3,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + } + else if (s2==-5) { + for (j=1,survp=0. ; j<=2; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + } + else{ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ + } + /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ + /*if(lli ==000.0)*/ + /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /* if (lli < log(mytinydouble)){ */ + /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ + /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ + /* } */ + } /* end of wave */ + } /* end of individual */ + } else if(mle==2){ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d<=dh[mi][i]; d++){ + newm=savm; + agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; + cov[2]=agexact; + if(nagesqr==1) + cov[3]= agexact*agexact; + for (kk=1; kk<=cptcovage;kk++) { + cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; + } + out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, + 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); + savm=oldm; + oldm=newm; + } /* end mult */ - s1=s[mw[mi][i]][i]; - s2=s[mw[mi+1][i]][i]; - bbh=(double)bh[mi][i]/(double)stepm; - lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ - ipmx +=1; - sw += weight[i]; - ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - } /* end of wave */ - } /* end of individual */ - } else if(mle==3){ /* exponential inter-extrapolation */ - for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; - for(mi=1; mi<= wav[i]-1; mi++){ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - for(d=0; d(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + } /* end of wave */ + } /* end of individual */ + } else if(mle==3){ /* exponential inter-extrapolation */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ - ipmx +=1; - sw += weight[i]; - ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - } /* end of wave */ - } /* end of individual */ - }else if (mle==4){ /* ml=4 no inter-extrapolation */ - for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; - for(mi=1; mi<= wav[i]-1; mi++){ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - for(d=0; d1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + } /* end of wave */ + } /* end of individual */ + }else if (mle==4){ /* ml=4 no inter-extrapolation */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d nlstate){ - lli=log(out[s1][s2] - savm[s1][s2]); - } else if ( s2==-1 ) { /* alive */ - for (j=1,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= log(survp); - }else{ - lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ - } - ipmx +=1; - sw += weight[i]; - ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; + if( s2 > nlstate){ + lli=log(out[s1][s2] - savm[s1][s2]); + } else if ( s2==-1 ) { /* alive */ + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= log(survp); + }else{ + lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ + } + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ - } /* end of wave */ - } /* end of individual */ - }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ - for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; - for(mi=1; mi<= wav[i]-1; mi++){ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - for(d=0; d nlstate && (mle <5) ){ /* Jackson */ - lli=log(out[s1][s2] - savm[s1][s2]); + lli=log(out[s1][s2] - savm[s1][s2]); } else if ( s2==-1 ) { /* alive */ - for (j=1,survp=0. ; j<=nlstate; j++) - survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; - lli= log(survp); + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); }else if (mle==1){ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ } else if(mle==2){ - lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ + lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ } else if(mle==3){ /* exponential inter-extrapolation */ - lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ + lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ } else if (mle==4){ /* mle=4 no inter-extrapolation */ - lli=log(out[s1][s2]); /* Original formula */ + lli=log(out[s1][s2]); /* Original formula */ } else{ /* mle=0 back to 1 */ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ - /*lli=log(out[s1][s2]); */ /* Original formula */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /*lli=log(out[s1][s2]); */ /* Original formula */ } /* End of if */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ + fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, - 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); - for(k=1,llt=0.,l=0.; k<=nlstate; k++){ - llt +=ll[k]*gipmx/gsw; - fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); - } - fprintf(ficresilk," %10.6f\n", -llt); + num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, + 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); + for(k=1,llt=0.,l=0.; k<=nlstate; k++){ + llt +=ll[k]*gipmx/gsw; + fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); + } + fprintf(ficresilk," %10.6f\n", -llt); } } /* end of wave */ } /* end of individual */ @@ -3685,319 +3772,319 @@ void pstamp(FILE *fichier) int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies */ - int i, m, jk, j1, bool, z1,j; - int iind=0, iage=0; - int mi; /* Effective wave */ - int first; - double ***freq; /* Frequencies */ - double *pp, **prop, *posprop, *pospropt; - double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; - char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; - double agebegin, ageend; - - pp=vector(1,nlstate); - prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); - posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ - pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ - /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ - strcpy(fileresp,"P_"); - strcat(fileresp,fileresu); - /*strcat(fileresphtm,fileresu);*/ - if((ficresp=fopen(fileresp,"w"))==NULL) { - printf("Problem with prevalence resultfile: %s\n", fileresp); - fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); - exit(0); - } + int i, m, jk, j1, bool, z1,j; + int iind=0, iage=0; + int mi; /* Effective wave */ + int first; + double ***freq; /* Frequencies */ + double *pp, **prop, *posprop, *pospropt; + double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; + char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; + double agebegin, ageend; + + pp=vector(1,nlstate); + prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ + pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ + /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ + strcpy(fileresp,"P_"); + strcat(fileresp,fileresu); + /*strcat(fileresphtm,fileresu);*/ + if((ficresp=fopen(fileresp,"w"))==NULL) { + printf("Problem with prevalence resultfile: %s\n", fileresp); + fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); + exit(0); + } - strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); - if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { - printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); - fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); - fflush(ficlog); - exit(70); - } - else{ - fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \ + strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); + if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { + printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \
\n\ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ - fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); - } - fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition

\n",fileresphtm, fileresphtm); + fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition

\n",fileresphtm, fileresphtm); - strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); - if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { - printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); - fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); - fflush(ficlog); - exit(70); - } - else{ - fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \ + strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); + if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { + printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \
\n\ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ - fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); - } - fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions by age at begin of transition

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions by age at begin of transition

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); - freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); - j1=0; + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); + j1=0; - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} + j=cptcoveff; + if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; + first=1; - /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: - reference=low_education V1=0,V2=0 - med_educ V1=1 V2=0, - high_educ V1=0 V2=1 - Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff - */ + /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: + reference=low_education V1=0,V2=0 + med_educ V1=1 V2=0, + high_educ V1=0 V2=1 + Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff + */ - for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ - posproptt=0.; - /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); - scanf("%d", i);*/ - for (i=-5; i<=nlstate+ndeath; i++) - for (jk=-5; jk<=nlstate+ndeath; jk++) - for(m=iagemin; m <= iagemax+3; m++) - freq[i][jk][m]=0; + for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ + posproptt=0.; + /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); + scanf("%d", i);*/ + for (i=-5; i<=nlstate+ndeath; i++) + for (jk=-5; jk<=nlstate+ndeath; jk++) + for(m=iagemin; m <= iagemax+3; m++) + freq[i][jk][m]=0; - for (i=1; i<=nlstate; i++) { - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0; - posprop[i]=0; - pospropt[i]=0; - } + for (i=1; i<=nlstate; i++) { + for(m=iagemin; m <= iagemax+3; m++) + prop[i][m]=0; + posprop[i]=0; + pospropt[i]=0; + } - dateintsum=0; - k2cpt=0; - - for (iind=1; iind<=imx; iind++) { /* For each individual iind */ - bool=1; - if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - for (z1=1; z1<=cptcoveff; z1++) { - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ - /* Tests if the value of each of the covariates of i is equal to filter j1 */ - bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + dateintsum=0; + k2cpt=0; + /* For that comination of covariate j1, we count and print the frequencies */ + for (iind=1; iind<=imx; iind++) { /* For each individual iind */ + bool=1; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) { + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ + /* Tests if the value of each of the covariates of i is equal to filter j1 */ + bool=0; + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ - } - } /* end z1 */ - } /* cptcovn > 0 */ - - if (bool==1){ - /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; mi=firstpass && m <=lastpass){ - k2=anint[m][iind]+(mint[m][iind]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ - if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ - if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ - prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ - if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { - dateintsum=dateintsum+k2; - k2cpt++; - /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ - } - /*}*/ - } /* end m */ - } /* end bool */ - } /* end iind = 1 to imx */ + /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ + } + } /* end z1 */ + } /* cptcovn > 0 */ + + if (bool==1){ /* We selected an individual iin satisfying combination j1 */ + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][iind]+(mint[m][iind]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + if (m1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + /*}*/ + } /* end m */ + } /* end bool */ + } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as - freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ + freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ - /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - pstamp(ficresp); - if (cptcovn>0) { - fprintf(ficresp, "\n#********** Variable "); - fprintf(ficresphtm, "\n

********** Variable "); - fprintf(ficresphtmfr, "\n

********** Variable "); - for (z1=1; z1<=cptcoveff; z1++){ - fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - } - fprintf(ficresp, "**********\n#"); - fprintf(ficresphtm, "**********

\n"); - fprintf(ficresphtmfr, "**********\n"); - fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "**********\n"); - } - fprintf(ficresphtm,""); - for(i=1; i<=nlstate;i++) { - fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); - fprintf(ficresphtm, "",i,i); - } - fprintf(ficresp, "\n"); - fprintf(ficresphtm, "\n"); + /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ + pstamp(ficresp); + if (cptcovn>0) { + fprintf(ficresp, "\n#********** Variable "); + fprintf(ficresphtm, "\n

********** Variable "); + fprintf(ficresphtmfr, "\n

********** Variable "); + for (z1=1; z1<=cptcoveff; z1++){ + fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } + fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "**********

\n"); + fprintf(ficresphtmfr, "**********\n"); + fprintf(ficlog, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficlog, "**********\n"); + } + fprintf(ficresphtm,"
AgePrev(%d)N(%d)N
"); + for(i=1; i<=nlstate;i++) { + fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); + fprintf(ficresphtm, "",i,i); + } + fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); - /* Header of frequency table by age */ - fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); - fprintf(ficresphtmfr," "); - for(jk=-1; jk <=nlstate+ndeath; jk++){ - for(m=-1; m <=nlstate+ndeath; m++){ - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",jk,m); - } - } - fprintf(ficresphtmfr, "\n"); + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
Age%d%d
"); + fprintf(ficresphtmfr," "); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",jk,m); + } + } + fprintf(ficresphtmfr, "\n"); - /* For each age */ - for(iage=iagemin; iage <= iagemax+3; iage++){ - fprintf(ficresphtm,""); - if(iage==iagemax+1){ - fprintf(ficlog,"1"); - fprintf(ficresphtmfr," "); - }else if(iage==iagemax+2){ - fprintf(ficlog,"0"); - fprintf(ficresphtmfr," "); - }else if(iage==iagemax+3){ - fprintf(ficlog,"Total"); - fprintf(ficresphtmfr," "); - }else{ - if(first==1){ - first=0; - printf("See log file for details...\n"); - } - fprintf(ficresphtmfr," ",iage); - fprintf(ficlog,"Age %d", iage); - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][iage]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][iage]; - if(pp[jk]>=1.e-10){ - if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[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++){ - /* posprop[jk]=0; */ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ - pp[jk] += freq[jk][m][iage]; - } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ - - for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ - pos += pp[jk]; /* pos is the total number of transitions until this age */ - posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ - pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state - from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ - } - for(jk=1; jk <=nlstate ; jk++){ - if(pos>=1.e-5){ - if(first==1) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - 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( iage <= iagemax){ - if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); - /*probs[iage][jk][j1]= pp[jk]/pos;*/ - /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ - } - else{ - fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); - fprintf(ficresphtm,"",iage, prop[jk][iage],pospropta); - } - } - pospropt[jk] +=posprop[jk]; - } /* end loop jk */ - /* pospropt=0.; */ - for(jk=-1; jk <=nlstate+ndeath; jk++){ - for(m=-1; m <=nlstate+ndeath; m++){ - if(freq[jk][m][iage] !=0 ) { /* minimizing output */ - if(first==1){ - printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); - } - if(jk!=0 && m!=0) - fprintf(ficresphtmfr," ",freq[jk][m][iage]); - } - } /* end loop jk */ - posproptt=0.; - for(jk=1; jk <=nlstate; jk++){ - posproptt += pospropt[jk]; - } - fprintf(ficresphtmfr,"\n "); - if(iage <= iagemax){ - fprintf(ficresp,"\n"); - fprintf(ficresphtm,"\n"); - } - if(first==1) - printf("Others in log...\n"); - fprintf(ficlog,"\n"); - } /* end loop age iage */ - fprintf(ficresphtm,""); - for(jk=1; jk <=nlstate ; jk++){ - if(posproptt < 1.e-5){ - fprintf(ficresphtm,"",pospropt[jk],posproptt); - }else{ - fprintf(ficresphtm,"",pospropt[jk]/posproptt,pospropt[jk],posproptt); - } - } - fprintf(ficresphtm,"\n"); - fprintf(ficresphtm,"
Age%d%d
0
Unknown
Total
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f
TotNanq%.0f%.0f%.5f%.0f%.0f
\n"); - fprintf(ficresphtmfr,"\n"); - if(posproptt < 1.e-5){ - fprintf(ficresphtm,"\n

This combination (%d) is not valid and no result will be produced

",j1); - fprintf(ficresphtmfr,"\n

This combination (%d) is not valid and no result will be produced

",j1); - fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); - invalidvarcomb[j1]=1; - }else{ - fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); - invalidvarcomb[j1]=0; - } - fprintf(ficresphtmfr,"\n"); - } /* end selected combination of covariate j1 */ - dateintmean=dateintsum/k2cpt; + /* For each age */ + for(iage=iagemin; iage <= iagemax+3; iage++){ + fprintf(ficresphtm,""); + if(iage==iagemax+1){ + fprintf(ficlog,"1"); + fprintf(ficresphtmfr,"0 "); + }else if(iage==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr,"Unknown "); + }else if(iage==iagemax+3){ + fprintf(ficlog,"Total"); + fprintf(ficresphtmfr,"Total "); + }else{ + if(first==1){ + first=0; + printf("See log file for details...\n"); + } + fprintf(ficresphtmfr,"%d ",iage); + fprintf(ficlog,"Age %d", iage); + } + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) + pp[jk] += freq[jk][m][iage]; + } + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pos=0; m <=0 ; m++) + pos += freq[jk][m][iage]; + if(pp[jk]>=1.e-10){ + if(first==1){ + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[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++){ + /* posprop[jk]=0; */ + for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ + pp[jk] += freq[jk][m][iage]; + } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ + + for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ + pos += pp[jk]; /* pos is the total number of transitions until this age */ + posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state + from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ + } + for(jk=1; jk <=nlstate ; jk++){ + if(pos>=1.e-5){ + if(first==1) + printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + 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( iage <= iagemax){ + if(pos>=1.e-5){ + fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + fprintf(ficresphtm,"%d%.5f%.0f%.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); + /*probs[iage][jk][j1]= pp[jk]/pos;*/ + /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ + } + else{ + fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); + fprintf(ficresphtm,"%dNaNq%.0f%.0f",iage, prop[jk][iage],pospropta); + } + } + pospropt[jk] +=posprop[jk]; + } /* end loop jk */ + /* pospropt=0.; */ + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][iage] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); + } + if(jk!=0 && m!=0) + fprintf(ficresphtmfr,"%.0f ",freq[jk][m][iage]); + } + } /* end loop jk */ + posproptt=0.; + for(jk=1; jk <=nlstate; jk++){ + posproptt += pospropt[jk]; + } + fprintf(ficresphtmfr,"\n "); + if(iage <= iagemax){ + fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); + } + if(first==1) + printf("Others in log...\n"); + fprintf(ficlog,"\n"); + } /* end loop age iage */ + fprintf(ficresphtm,"Tot"); + for(jk=1; jk <=nlstate ; jk++){ + if(posproptt < 1.e-5){ + fprintf(ficresphtm,"Nanq%.0f%.0f",pospropt[jk],posproptt); + }else{ + fprintf(ficresphtm,"%.5f%.0f%.0f",pospropt[jk]/posproptt,pospropt[jk],posproptt); + } + } + fprintf(ficresphtm,"\n"); + fprintf(ficresphtm,"\n"); + fprintf(ficresphtmfr,"\n"); + if(posproptt < 1.e-5){ + fprintf(ficresphtm,"\n

This combination (%d) is not valid and no result will be produced

",j1); + fprintf(ficresphtmfr,"\n

This combination (%d) is not valid and no result will be produced

",j1); + fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); + invalidvarcomb[j1]=1; + }else{ + fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); + invalidvarcomb[j1]=0; + } + fprintf(ficresphtmfr,"\n"); + } /* end selected combination of covariate j1 */ + dateintmean=dateintsum/k2cpt; - fclose(ficresp); - fclose(ficresphtm); - fclose(ficresphtmfr); - free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); - free_vector(pospropt,1,nlstate); - free_vector(posprop,1,nlstate); - free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); - free_vector(pp,1,nlstate); - /* End of Freq */ + fclose(ficresp); + fclose(ficresphtm); + fclose(ficresphtmfr); + free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); + free_vector(pospropt,1,nlstate); + free_vector(posprop,1,nlstate); + free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); + free_vector(pp,1,nlstate); + /* End of Freq */ } /************ Prevalence ********************/ @@ -4104,173 +4191,199 @@ void concatwav(int wav[], int **dh, int and mw[mi+1][i]. dh depends on stepm. */ - int i, mi, m; + int i=0, mi=0, m=0, mli=0; /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; double sum=0., jmean=0.;*/ - int first, firstwo, firsthree, firstfour; + int first=0, firstwo=0, firsthree=0, firstfour=0, firstfiv=0; int j, k=0,jk, ju, jl; double sum=0.; - first=0; - firstwo=0; - firsthree=0; - firstfour=0; jmin=100000; jmax=-1; jmean=0.; + +/* Treating live states */ for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ - mi=0; + mi=0; /* First valid wave */ + mli=0; /* Last valid wave? */ m=firstpass; while(s[m][i] <= nlstate){ /* a live state */ - if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ - mw[++mi][i]=m; - } - if(m >=lastpass){ - if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ - if(firsthree == 0){ - printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); - firsthree=1; - } - fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); - mw[++mi][i]=m; - } - if(s[m][i]==-2){ /* Vital status is really unknown */ - nbwarn++; - if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ - printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); - fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); - } - break; - } - break; + if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ + mli=m-1;/* mw[++mi][i]=m-1; */ + }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ + mw[++mi][i]=m; + mli=m; + } /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ + if(m < lastpass){ /* m < lastpass, standard case */ + m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ } - else - m++; + else{ /* m >= lastpass, eventual special issue with warning */ +#ifdef UNKNOWNSTATUSNOTCONTRIBUTING + break; +#else + if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ + if(firsthree == 0){ + printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + firsthree=1; + } + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); + mw[++mi][i]=m; + mli=m; + } + if(s[m][i]==-2){ /* Vital status is really unknown */ + nbwarn++; + if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ + printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); + } + break; + } + break; +#endif + }/* End m >= lastpass */ }/* end while */ - + + /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */ /* After last pass */ +/* Treating death states */ if (s[m][i] > nlstate){ /* In a death state */ + /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */ + /* } */ mi++; /* Death is another wave */ /* if(mi==0) never been interviewed correctly before death */ - /* Only death is a correct wave */ + /* Only death is a correct wave */ mw[mi][i]=m; - }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */ + } +#ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE + else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */ /* m++; */ /* mi++; */ /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */ /* mw[mi][i]=m; */ - nberr++; if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ - if(firstwo==0){ - printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); - firstwo=1; - } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ + nbwarn++; + if(firstfiv==0){ + printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + firstfiv=1; + }else{ + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + } + }else{ /* Death occured afer last wave potential bias */ + nberr++; + if(firstwo==0){ + printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + firstwo=1; + } + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + } }else{ /* end date of interview is known */ - /* death is known but not confirmed by death status at any wave */ - if(firstfour==0){ - printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); - firstfour=1; - } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + /* death is known but not confirmed by death status at any wave */ + if(firstfour==0){ + printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + firstfour=1; + } + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); } - } - wav[i]=mi; + } /* end if date of death is known */ +#endif + wav[i]=mi; /* mi should be the last effective wave (or mli) */ + /* wav[i]=mw[mi][i]; */ if(mi==0){ nbwarn++; if(first==0){ - printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); - first=1; + printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); + first=1; } if(first==1){ - fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); + fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); } } /* end mi==0 */ } /* End individuals */ /* wav and mw are no more changed */ - + for(i=1; i<=imx; i++){ for(mi=1; mi nlstate) { /* A death */ - if (agedc[i] < 2*AGESUP) { - j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); - if(j==0) j=1; /* Survives at least one month after exam */ - else if(j<0){ - nberr++; - printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - j=1; /* Temporary Dangerous patch */ - printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); - fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); - } - k=k+1; - if (j >= jmax){ - jmax=j; - ijmax=i; - } - if (j <= jmin){ - jmin=j; - ijmin=i; - } - sum=sum+j; - /*if (j<0) printf("j=%d num=%d \n",j,i);*/ - /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ - } - } - else{ - j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + if (s[mw[mi+1][i]][i] > nlstate) { /* A death */ + if (agedc[i] < 2*AGESUP) { + j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); + if(j==0) j=1; /* Survives at least one month after exam */ + else if(j<0){ + nberr++; + printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + j=1; /* Temporary Dangerous patch */ + printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); + fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); + } + k=k+1; + if (j >= jmax){ + jmax=j; + ijmax=i; + } + if (j <= jmin){ + jmin=j; + ijmin=i; + } + sum=sum+j; + /*if (j<0) printf("j=%d num=%d \n",j,i);*/ + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ + } + } + else{ + j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); /* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ - - k=k+1; - if (j >= jmax) { - jmax=j; - ijmax=i; - } - else if (j <= jmin){ - jmin=j; - ijmin=i; - } - /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ - /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ - if(j<0){ - nberr++; - printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - } - sum=sum+j; - } - jk= j/stepm; - jl= j -jk*stepm; - ju= j -(jk+1)*stepm; - if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ - if(jl==0){ - dh[mi][i]=jk; - bh[mi][i]=0; - }else{ /* We want a negative bias in order to only have interpolation ie - * to avoid the price of an extra matrix product in likelihood */ - dh[mi][i]=jk+1; - bh[mi][i]=ju; - } - }else{ - if(jl <= -ju){ - dh[mi][i]=jk; - bh[mi][i]=jl; /* bias is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ - } - else{ - dh[mi][i]=jk+1; - bh[mi][i]=ju; - } - if(dh[mi][i]==0){ - dh[mi][i]=1; /* At least one step */ - bh[mi][i]=ju; /* At least one step */ - /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ - } - } /* end if mle */ + + k=k+1; + if (j >= jmax) { + jmax=j; + ijmax=i; + } + else if (j <= jmin){ + jmin=j; + ijmin=i; + } + /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ + /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ + if(j<0){ + nberr++; + printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + } + sum=sum+j; + } + jk= j/stepm; + jl= j -jk*stepm; + ju= j -(jk+1)*stepm; + if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ + if(jl==0){ + dh[mi][i]=jk; + bh[mi][i]=0; + }else{ /* We want a negative bias in order to only have interpolation ie + * to avoid the price of an extra matrix product in likelihood */ + dh[mi][i]=jk+1; + bh[mi][i]=ju; + } + }else{ + if(jl <= -ju){ + dh[mi][i]=jk; + bh[mi][i]=jl; /* bias is positive if real duration + * is higher than the multiple of stepm and negative otherwise. + */ + } + else{ + dh[mi][i]=jk+1; + bh[mi][i]=ju; + } + if(dh[mi][i]==0){ + dh[mi][i]=1; /* At least one step */ + bh[mi][i]=ju; /* At least one step */ + /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ + } + } /* end if mle */ } } /* end wave */ } @@ -4345,7 +4458,7 @@ void concatwav(int wav[], int **dh, int undefined. Usually 3: -1, 0 and 1. */ } /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for - historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ + * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ } /* Ndum[-1] number of undefined modalities */ /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ @@ -4366,11 +4479,11 @@ void concatwav(int wav[], int **dh, int if (Ndum[i] == 0) { /* If nobody responded to this modality k */ break; } - ij++; - nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ - cptcode = ij; /* New max modality for covar j */ + ij++; + nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ + cptcode = ij; /* New max modality for covar j */ } /* end of loop on modality i=-1 to 1 or more */ - + /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ /* /\*recode from 0 *\/ */ /* k is a modality. If we have model=V1+V1*sex */ @@ -4392,7 +4505,7 @@ void concatwav(int wav[], int **dh, int /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ Ndum[ij]++; /* Might be supersed V1 + V1*age */ - } + } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */ ij=0; for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ @@ -4407,8 +4520,9 @@ void concatwav(int wav[], int **dh, int } /* ij--; */ /* cptcoveff=ij; /\*Number of total covariates*\/ */ - *cptcov=ij; /*Number of total covariates*/ - + *cptcov=ij; /*Number of total real effective covariates: effective + * because they can be excluded from the model and real + * if in the model but excluded because missing values*/ } @@ -5259,7 +5373,7 @@ To be simple, these graphs help to under tj = (int) pow(2,cptcoveff); if (cptcovn<1) {tj=1;ncodemax[1]=1;} j1=0; - for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */ + for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates or only once*/ if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); @@ -5700,7 +5814,7 @@ true period expectancies (those weighted } /******************* Gnuplot file **************/ - void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){ +void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){ char dirfileres[132],optfileres[132]; char gplotcondition[132]; @@ -5963,7 +6077,7 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) /* Survival functions (period) from state i in state j by final state j */ for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ - + fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ @@ -5975,8 +6089,8 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; } fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); @@ -5986,24 +6100,24 @@ unset log y\n plot [%.f:%.f] ", ageminpar, agemaxpar); k=3; for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ - if(j==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); - else - fprintf(ficgp,", '' "); - l=(nlstate+ndeath)*(cpt-1) +j; - fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l); - /* for (i=2; i<= nlstate+ndeath ; i ++) */ - /* fprintf(ficgp,"+$%d",k+l+i-1); */ - fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j); + if(j==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); + else + fprintf(ficgp,", '' "); + l=(nlstate+ndeath)*(cpt-1) +j; + fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l); + /* for (i=2; i<= nlstate+ndeath ; i ++) */ + /* fprintf(ficgp,"+$%d",k+l+i-1); */ + fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j); } /* nlstate */ fprintf(ficgp,", '' "); fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1); for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ - l=(nlstate+ndeath)*(cpt-1) +j; - if(j < nlstate) - fprintf(ficgp,"$%d +",k+l); - else - fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt); + l=(nlstate+ndeath)*(cpt-1) +j; + if(j < nlstate) + fprintf(ficgp,"$%d +",k+l); + else + fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt); } fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ @@ -6013,7 +6127,7 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) /* CV preval stable (period) for each covariate */ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ - + fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ @@ -6025,32 +6139,32 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } fprintf(ficgp,"\n#\n"); if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; } - + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ -set ter svg size 640, 480\n\ -unset log y\n\ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); k=3; /* Offset */ for (i=1; i<= nlstate ; i ++){ - if(i==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); - else - fprintf(ficgp,", '' "); - l=(nlstate+ndeath)*(i-1)+1; - fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); - for (j=2; j<= nlstate ; j ++) - fprintf(ficgp,"+$%d",k+l+j-1); - fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); + if(i==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); + else + fprintf(ficgp,", '' "); + l=(nlstate+ndeath)*(i-1)+1; + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); + for (j=2; j<= nlstate ; j ++) + fprintf(ficgp,"+$%d",k+l+j-1); + fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); } /* nlstate */ fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ - - + + /* 7eme */ if(backcast == 1){ /* CV back preval stable (period) for each covariate */ @@ -6062,37 +6176,37 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ - vlv= nbcode[Tvaraff[k]][lv]; - fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); - } - fprintf(ficgp,"\n#\n"); - if(invalidvarcomb[k1]){ - fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); - continue; - } + vlv= nbcode[Tvaraff[k]][lv]; + fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); + } + fprintf(ficgp,"\n#\n"); + if(invalidvarcomb[k1]){ + fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); + continue; + } - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter svg size 640, 480\n \ unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); - k=3; /* Offset */ - for (i=1; i<= nlstate ; i ++){ - if(i==1) - fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); - else - fprintf(ficgp,", '' "); - /* l=(nlstate+ndeath)*(i-1)+1; */ - l=(nlstate+ndeath)*(cpt-1)+1; - /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ - /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ - fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */ - /* for (j=2; j<= nlstate ; j ++) */ - /* fprintf(ficgp,"+$%d",k+l+j-1); */ - /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */ - fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt); - } /* nlstate */ - fprintf(ficgp,"\nset out\n"); + k=3; /* Offset */ + for (i=1; i<= nlstate ; i ++){ + if(i==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_")); + else + fprintf(ficgp,", '' "); + /* l=(nlstate+ndeath)*(i-1)+1; */ + l=(nlstate+ndeath)*(cpt-1)+1; + /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ + /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */ + /* for (j=2; j<= nlstate ; j ++) */ + /* fprintf(ficgp,"+$%d",k+l+j-1); */ + /* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */ + fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt); + } /* nlstate */ + fprintf(ficgp,"\nset out\n"); } /* end cpt state*/ } /* end covariate */ } /* End if backcast */ @@ -6118,11 +6232,11 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) continue; } - fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); - fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ -set ter svg size 640, 480\n \ -unset log y\n \ + fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); + fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ +set ter svg size 640, 480\n \ +unset log y\n \ plot [%.f:%.f] ", ageminpar, agemaxpar); for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ @@ -6196,17 +6310,17 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) fprintf(ficgp,"# initial state %d\n",i); for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { - fprintf(ficgp,"# current state %d\n",k); - for(j=1; j <=ncovmodel; j++){ - fprintf(ficgp,"p%d=%f; ",jk,p[jk]); - jk++; - } - fprintf(ficgp,"\n"); + fprintf(ficgp,"# current state %d\n",k); + for(j=1; j <=ncovmodel; j++){ + fprintf(ficgp,"p%d=%f; ",jk,p[jk]); + jk++; + } + fprintf(ficgp,"\n"); } } } fprintf(ficgp,"##############\n#\n"); - + /*goto avoid;*/ fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n"); fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); @@ -7196,8 +7310,8 @@ int readdata(char datafile[], int firsto } trimbb(linetmp,line); /* Trims multiple blanks in line */ strcpy(line, linetmp); - -/* Loops on waves */ + + /* Loops on waves */ for (j=maxwav;j>=1;j--){ cutv(stra, strb, line, ' '); if(strb[0]=='.') { /* Missing status */ @@ -7212,19 +7326,21 @@ int readdata(char datafile[], int firsto return 1; } } + s[j][i]=lval; - + + /* Date of Interview */ strcpy(line,stra); cutv(stra, strb,line,' '); if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ } else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ - month=99; - year=9999; + month=99; + year=9999; }else{ - printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); - return 1; + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); + return 1; } anint[j][i]= (double) year; mint[j][i]= (double)month; @@ -7240,8 +7356,8 @@ int readdata(char datafile[], int firsto year=9999; }else{ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); - return 1; + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; } andc[i]=(double) year; moisdc[i]=(double) month; @@ -7257,12 +7373,12 @@ int readdata(char datafile[], int firsto }else{ printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); - return 1; + return 1; } if (year==9999) { printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line); fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); - return 1; + return 1; } annais[i]=(double)(year); @@ -7286,18 +7402,18 @@ int readdata(char datafile[], int firsto for (j=ncovcol;j>=1;j--){ cutv(stra, strb,line,' '); if(strb[0]=='.') { /* Missing covariate value */ - lval=-1; + lval=-1; }else{ - errno=0; - lval=strtol(strb,&endptr,10); - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); - fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); - return 1; - } + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); + return 1; + } } if(lval <-1 || lval >1){ - printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ For example, for multinomial values like 1, 2 and 3,\n \ @@ -7306,7 +7422,7 @@ int readdata(char datafile[], int firsto and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j); - fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ For example, for multinomial values like 1, 2 and 3,\n \ @@ -7315,7 +7431,7 @@ int readdata(char datafile[], int firsto and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j);fflush(ficlog); - return 1; + return 1; } covar[j][i]=(double)(lval); strcpy(line,stra); @@ -7339,13 +7455,11 @@ int readdata(char datafile[], int firsto return (0); /* endread: */ - printf("Exiting readdata: "); - fclose(fic); - return (1); - - - + printf("Exiting readdata: "); + fclose(fic); + return (1); } + void removespace(char *str) { char *p1 = str, *p2 = str; do @@ -7354,20 +7468,21 @@ void removespace(char *str) { while (*p1++ == *p2++); } -int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns: - * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age - * - nagesqr = 1 if age*age in the model, otherwise 0. - * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age - * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age - * - cptcovage number of covariates with age*products =2 - * - cptcovs number of simple covariates - * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 - * which is a new column after the 9 (ncovcol) variables. - * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual - * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage - * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6. - * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 . - */ +int decodemodel ( char model[], int lastobs) + /**< This routine decode the model and returns: + * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age + * - nagesqr = 1 if age*age in the model, otherwise 0. + * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age + * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age + * - cptcovage number of covariates with age*products =2 + * - cptcovs number of simple covariates + * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 + * which is a new column after the 9 (ncovcol) variables. + * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual + * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage + * Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6. + * - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 . + */ { int i, j, k, ks; int j1, k1, k2; @@ -7392,34 +7507,34 @@ int decodemodel ( char model[], int last if ((strpt=strstr(model,"age*age")) !=0){ printf(" strpt=%s, model=%s\n",strpt, model); if(strpt != model){ - printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); - fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ + fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ 'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ corresponding column of parameters.\n",model); fflush(ficlog); - return 1; - } + return 1; + } nagesqr=1; if (strstr(model,"+age*age") !=0) - substrchaine(modelsav, model, "+age*age"); + substrchaine(modelsav, model, "+age*age"); else if (strstr(model,"age*age+") !=0) - substrchaine(modelsav, model, "age*age+"); + substrchaine(modelsav, model, "age*age+"); else - substrchaine(modelsav, model, "age*age"); + substrchaine(modelsav, model, "age*age"); }else nagesqr=0; if (strlen(modelsav) >1){ j=nbocc(modelsav,'+'); /**< j=Number of '+' */ j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */ - cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =2 */ + cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =5-3=2 */ cptcovt= j+1; /* Number of total covariates in the model, not including - * cst, age and age*age - * V1+V1*age+ V3 + V3*V4+age*age=> 4*/ - /* including age products which are counted in cptcovage. - * but the covariates which are products must be treated - * separately: ncovn=4- 2=2 (V1+V3). */ + * cst, age and age*age + * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/ + /* including age products which are counted in cptcovage. + * but the covariates which are products must be treated + * separately: ncovn=4- 2=2 (V1+V3). */ cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */ cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ @@ -7431,8 +7546,8 @@ int decodemodel ( char model[], int last * k= 1 2 3 4 5 6 7 8 * cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8 * covar[k,i], value of kth covariate if not including age for individual i: - * covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8) - * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8 + * covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8) + * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[2]=1 Tvar[4]=3 Tvar[8]=8 * if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and * Tage[++cptcovage]=k * if products, new covar are created after ncovcol with k1 @@ -7476,62 +7591,62 @@ int decodemodel ( char model[], int last Tvar[k]=0; cptcovage=0; for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ - cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' - modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ - if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ - /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ - /*scanf("%d",i);*/ - if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ - cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ - if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ - /* covar is not filled and then is empty */ - cptcovprod--; - cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ - Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ - cptcovage++; /* Sums the number of covariates which include age as a product */ - Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ - /*printf("stre=%s ", stre);*/ - } else if (strcmp(strd,"age")==0) { /* or age*Vn */ - cptcovprod--; - cutl(stre,strb,strc,'V'); - Tvar[k]=atoi(stre); - cptcovage++; - Tage[cptcovage]=k; - } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ - /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ - cptcovn++; - cptcovprodnoage++;k1++; - cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ - Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but - because this model-covariate is a construction we invent a new column - ncovcol + k1 - If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 - Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ - cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ - Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ - Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ - Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ - k2=k2+2; - Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ - Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ - for (i=1; i<=lastobs;i++){ - /* Computes the new covariate which is a product of - covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ - covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; - } - } /* End age is not in the model */ - } /* End if model includes a product */ - else { /* no more sum */ - /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ - /* scanf("%d",i);*/ - cutl(strd,strc,strb,'V'); - ks++; /**< Number of simple covariates */ - cptcovn++; - Tvar[k]=atoi(strd); - } - strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ - /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); - scanf("%d",i);*/ + cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' + modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ + if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ + /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ + /*scanf("%d",i);*/ + if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ + cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ + if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ + /* covar is not filled and then is empty */ + cptcovprod--; + cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ + Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ + cptcovage++; /* Sums the number of covariates which include age as a product */ + Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ + /*printf("stre=%s ", stre);*/ + } else if (strcmp(strd,"age")==0) { /* or age*Vn */ + cptcovprod--; + cutl(stre,strb,strc,'V'); + Tvar[k]=atoi(stre); + cptcovage++; + Tage[cptcovage]=k; + } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ + /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ + cptcovn++; + cptcovprodnoage++;k1++; + cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ + Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but + because this model-covariate is a construction we invent a new column + ncovcol + k1 + If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 + Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ + cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ + Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ + Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ + Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ + k2=k2+2; + Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ + Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ + for (i=1; i<=lastobs;i++){ + /* Computes the new covariate which is a product of + covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ + covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; + } + } /* End age is not in the model */ + } /* End if model includes a product */ + else { /* no more sum */ + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutl(strd,strc,strb,'V'); + ks++; /**< Number of simple covariates */ + cptcovn++; + Tvar[k]=atoi(strd); + } + strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ + /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); + scanf("%d",i);*/ } /* end of loop + on total covariates */ } /* end if strlen(modelsave == 0) age*age might exist */ } /* end if strlen(model == 0) */ @@ -7540,16 +7655,17 @@ int decodemodel ( char model[], int last If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); - printf("cptcovprod=%d ", cptcovprod); - fprintf(ficlog,"cptcovprod=%d ", cptcovprod); + printf("cptcovprod=%d ", cptcovprod); + fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - scanf("%d ",i);*/ + scanf("%d ",i);*/ +/* Dispatching in quantitative and time varying covariates */ return (0); /* with covar[new additional covariate if product] and Tage if age */ /*endread:*/ - printf("Exiting decodemodel: "); - return (1); + printf("Exiting decodemodel: "); + return (1); } int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) @@ -8742,7 +8858,7 @@ Please run with mle=-1 to get a correct printf("Problem writing new parameter file: %s\n", rfileres);goto end; fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end; } - fprintf(ficres,"#%s\n",version); + fprintf(ficres,"#IMaCh version %s\n",version); } /* End of mle != -3 */ /* Main data @@ -8810,7 +8926,7 @@ Please run with mle=-1 to get a correct /* Main decodemodel */ - if(decodemodel(model, lastobs) == 1) + if(decodemodel(model, lastobs) == 1) /* In order to get Tvar[k] V4+V3+V5 p Tvar[1]@3 = {4, 3, 5}*/ goto end; if((double)(lastobs-imx)/(double)imx > 1.10){ @@ -9067,7 +9183,7 @@ Interval (in months) between two waves: ageexmed=vector(1,n); agecens=vector(1,n); dcwave=ivector(1,n); - + for (i=1; i<=imx; i++){ dcwave[i]=-1; for (m=firstpass; m<=lastpass; m++) @@ -9567,9 +9683,9 @@ Please run with mle=-1 to get a correct ungetc(c,ficpar); fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj); - fprintf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); - fprintf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); - fprintf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); + fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); + fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); + fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); /* day and month of proj2 are not used but only year anproj2.*/