--- imach/src/imach.c 2014/12/23 13:26:59 1.171 +++ imach/src/imach.c 2015/08/18 13:32:00 1.194 @@ -1,6 +1,89 @@ -/* $Id: imach.c,v 1.171 2014/12/23 13:26:59 brouard Exp $ +/* $Id: imach.c,v 1.194 2015/08/18 13:32:00 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.194 2015/08/18 13:32:00 brouard + Summary: Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line. + + Revision 1.193 2015/08/04 07:17:42 brouard + Summary: 0.98q4 + + Revision 1.192 2015/07/16 16:49:02 brouard + Summary: Fixing some outputs + + Revision 1.191 2015/07/14 10:00:33 brouard + Summary: Some fixes + + Revision 1.190 2015/05/05 08:51:13 brouard + Summary: Adding digits in output parameters (7 digits instead of 6) + + Fix 1+age+. + + Revision 1.189 2015/04/30 14:45:16 brouard + Summary: 0.98q2 + + Revision 1.188 2015/04/30 08:27:53 brouard + *** empty log message *** + + Revision 1.187 2015/04/29 09:11:15 brouard + *** empty log message *** + + Revision 1.186 2015/04/23 12:01:52 brouard + Summary: V1*age is working now, version 0.98q1 + + Some codes had been disabled in order to simplify and Vn*age was + working in the optimization phase, ie, giving correct MLE parameters, + but, as usual, outputs were not correct and program core dumped. + + Revision 1.185 2015/03/11 13:26:42 brouard + Summary: Inclusion of compile and links command line for Intel Compiler + + Revision 1.184 2015/03/11 11:52:39 brouard + Summary: Back from Windows 8. Intel Compiler + + Revision 1.183 2015/03/10 20:34:32 brouard + Summary: 0.98q0, trying with directest, mnbrak fixed + + We use directest instead of original Powell test; probably no + incidence on the results, but better justifications; + We fixed Numerical Recipes mnbrak routine which was wrong and gave + wrong results. + + Revision 1.182 2015/02/12 08:19:57 brouard + Summary: Trying to keep directest which seems simpler and more general + Author: Nicolas Brouard + + Revision 1.181 2015/02/11 23:22:24 brouard + Summary: Comments on Powell added + + Author: + + Revision 1.180 2015/02/11 17:33:45 brouard + Summary: Finishing move from main to function (hpijx and prevalence_limit) + + Revision 1.179 2015/01/04 09:57:06 brouard + Summary: back to OS/X + + Revision 1.178 2015/01/04 09:35:48 brouard + *** empty log message *** + + Revision 1.177 2015/01/03 18:40:56 brouard + Summary: Still testing ilc32 on OSX + + Revision 1.176 2015/01/03 16:45:04 brouard + *** empty log message *** + + Revision 1.175 2015/01/03 16:33:42 brouard + *** empty log message *** + + Revision 1.174 2015/01/03 16:15:49 brouard + Summary: Still in cross-compilation + + Revision 1.173 2015/01/03 12:06:26 brouard + Summary: trying to detect cross-compilation + + Revision 1.172 2014/12/27 12:07:47 brouard + Summary: Back from Visual Studio and Intel, options for compiling for Windows XP + Revision 1.171 2014/12/23 13:26:59 brouard Summary: Back from Visual C @@ -528,7 +611,12 @@ end */ +/* #define DEBUG */ +/* #define DEBUGBRENT */ #define POWELL /* Instead of NLOPT */ +#define POWELLF1F3 /* Skip test */ +/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ +/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ #include #include @@ -537,6 +625,8 @@ #ifdef _WIN32 #include +#include +#include #else #include #endif @@ -598,6 +688,7 @@ typedef struct { #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 #define AGEBASE 40 +#define AGEOVERFLOW 1.e20 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ #ifdef _WIN32 #define DIRSEPARATOR '\\' @@ -609,15 +700,15 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.171 2014/12/23 13:26:59 brouard Exp $ */ +/* $Id: imach.c,v 1.194 2015/08/18 13:32:00 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98p, December 2014,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.171 $ $Date: 2014/12/23 13:26:59 $"; +char version[]="Imach version 0.98q5, August 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.194 $ $Date: 2015/08/18 13:32:00 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar=0, nforce=0; /* Number of variables, number of forces */ +int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ /* Number of covariates model=V2+V1+ V3*age+V2*V4 */ int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */ int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */ @@ -728,7 +819,7 @@ static double maxarg1,maxarg2; #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) #define rint(a) floor(a+0.5) /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */ -/* #define mytinydouble 1.0e-16 */ +#define mytinydouble 1.0e-16 /* #define DEQUAL(a,b) (fabs((a)-(b)) Tvar[1]= 2 */ int *Ndum; /** Freq of modality (tricode */ @@ -776,7 +873,7 @@ static int split( char *path, char *dirc the name of the file (name), its extension only (ext) and its first part of the name (finame) */ char *ss; /* pointer */ - int l1, l2; /* length counters */ + int l1=0, l2=0; /* length counters */ l1 = strlen(path ); /* length of path */ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); @@ -787,7 +884,11 @@ static int split( char *path, char *dirc printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ /* get current working directory */ /* extern char* getcwd ( char *buf , int len);*/ - if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { +#ifdef WIN32 + if (_getcwd( dirc, FILENAME_MAX ) == NULL ) { +#else + if (getcwd(dirc, FILENAME_MAX) == NULL) { +#endif return( GLOCK_ERROR_GETCWD ); } /* got dirc from getcwd*/ @@ -798,7 +899,7 @@ static int split( char *path, char *dirc if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); strcpy( name, ss ); /* save file name */ strncpy( dirc, path, l1 - l2 ); /* now the directory */ - dirc[l1-l2] = 0; /* add zero */ + dirc[l1-l2] = '\0'; /* add zero */ printf(" DIRC2 = %s \n",dirc); } /* We add a separator at the end of dirc if not exists */ @@ -850,11 +951,54 @@ char *trimbb(char *out, char *in) return s; } +/* char *substrchaine(char *out, char *in, char *chain) */ +/* { */ +/* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */ +/* char *s, *t; */ +/* t=in;s=out; */ +/* while ((*in != *chain) && (*in != '\0')){ */ +/* *out++ = *in++; */ +/* } */ + +/* /\* *in matches *chain *\/ */ +/* while ((*in++ == *chain++) && (*in != '\0')){ */ +/* printf("*in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* } */ +/* in--; chain--; */ +/* while ( (*in != '\0')){ */ +/* printf("Bef *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* *out++ = *in++; */ +/* printf("Aft *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ +/* } */ +/* *out='\0'; */ +/* out=s; */ +/* return out; */ +/* } */ +char *substrchaine(char *out, char *in, char *chain) +{ + /* Substract chain 'chain' from 'in', return and output 'out' */ + /* in="V1+V1*age+age*age+V2", chain="age*age" */ + + char *strloc; + + strcpy (out, in); + strloc = strstr(out, chain); /* strloc points to out at age*age+V2 */ + printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out); + if(strloc != NULL){ + /* will affect out */ /* strloc+strlenc(chain)=+V2 */ /* Will also work in Unicode */ + memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1); + /* strcpy (strloc, strloc +strlen(chain));*/ + } + printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out); + return out; +} + + char *cutl(char *blocc, char *alocc, char *in, char occ) { - /* cuts string in into blocc and alocc where blocc ends before first occurence of char 'occ' + /* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ' and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') - gives blocc="abcdef2ghi" and alocc="j". + gives blocc="abcdef" and alocc="ghi2j". If occ is not found blocc is null and alocc is equal to in. Returns blocc */ char *s, *t; @@ -880,7 +1024,7 @@ char *cutl(char *blocc, char *alocc, cha } char *cutv(char *blocc, char *alocc, char *in, char occ) { - /* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ' + /* cuts string in into blocc and alocc where blocc ends before LAST occurence of char 'occ' and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') gives blocc="abcdef2ghi" and alocc="j". If occ is not found blocc is null and alocc is equal to in. Returns alocc @@ -1191,7 +1335,13 @@ double f1dim(double x) /*****************brent *************************/ double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin) -{ +{ + /* Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is + * between ax and cx, and f(bx) is less than both f(ax) and f(cx) ), this routine isolates + * the minimum to a fractional precision of about tol using Brent’s method. The abscissa of + * the minimum is returned as xmin, and the minimum function value is returned as brent , the + * returned function value. + */ int iter; double a,b,d,etemp; double fu=0,fv,fw,fx; @@ -1244,19 +1394,19 @@ double brent(double ax, double bx, doubl if (fu <= fx) { if (u >= x) a=x; else b=x; SHFT(v,w,x,u) - SHFT(fv,fw,fx,fu) - } else { - if (u < x) a=u; else b=u; - if (fu <= fw || w == x) { - v=w; - w=u; - fv=fw; - fw=fu; - } else if (fu <= fv || v == x || v == w) { - v=u; - fv=fu; - } - } + SHFT(fv,fw,fx,fu) + } else { + if (u < x) a=u; else b=u; + if (fu <= fw || w == x) { + v=w; + w=u; + fv=fw; + fw=fu; + } else if (fu <= fv || v == x || v == w) { + v=u; + fv=fu; + } + } } nrerror("Too many iterations in brent"); *xmin=x; @@ -1267,25 +1417,44 @@ double brent(double ax, double bx, doubl void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, double (*func)(double)) -{ +{ /* Given a function func , and given distinct initial points ax and bx , this routine searches in +the downhill direction (defined by the function as evaluated at the initial points) and returns +new points ax , bx , cx that bracket a minimum of the function. Also returned are the function +values at the three points, fa, fb , and fc such that fa > fb and fb < fc. + */ double ulim,u,r,q, dum; double fu; - - *fa=(*func)(*ax); - *fb=(*func)(*bx); + + double scale=10.; + int iterscale=0; + + *fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ + *fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ + + + /* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */ + /* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */ + /* *bx = *ax - (*ax - *bx)/scale; */ + /* *fb=(*func)(*bx); /\* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) *\/ */ + /* } */ + if (*fb > *fa) { SHFT(dum,*ax,*bx,dum) - SHFT(dum,*fb,*fa,dum) - } + SHFT(dum,*fb,*fa,dum) + } *cx=(*bx)+GOLD*(*bx-*ax); *fc=(*func)(*cx); - while (*fb > *fc) { /* Declining fa, fb, fc */ +#ifdef DEBUG + printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); + fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); +#endif + while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */ r=(*bx-*ax)*(*fb-*fc); q=(*bx-*cx)*(*fb-*fa); u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ - (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */ - ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */ - if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */ + (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 */ fu=(*func)(u); #ifdef DEBUG /* f(x)=A(x-u)**2+f(u) */ @@ -1294,23 +1463,85 @@ void mnbrak(double *ax, double *bx, doub 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); + /* 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) */ + /* In that case, there is no bracket in the output! Routine is wrong with many consequences.*/ #endif +#ifdef MNBRAKORIGINAL +#else +/* if (fu > *fc) { */ +/* #ifdef DEBUG */ +/* printf("mnbrak4 fu > fc \n"); */ +/* fprintf(ficlog, "mnbrak4 fu > fc\n"); */ +/* #endif */ +/* /\* SHFT(u,*cx,*cx,u) /\\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\\/ *\/ */ +/* /\* SHFT(*fa,*fc,fu,*fc) /\\* (b, u, c) is a bracket while test fb > fc will be fu > fc will exit *\\/ *\/ */ +/* dum=u; /\* Shifting c and u *\/ */ +/* u = *cx; */ +/* *cx = dum; */ +/* dum = fu; */ +/* fu = *fc; */ +/* *fc =dum; */ +/* } else { /\* end *\/ */ +/* #ifdef DEBUG */ +/* printf("mnbrak3 fu < fc \n"); */ +/* fprintf(ficlog, "mnbrak3 fu < fc\n"); */ +/* #endif */ +/* dum=u; /\* Shifting c and u *\/ */ +/* u = *cx; */ +/* *cx = dum; */ +/* dum = fu; */ +/* fu = *fc; */ +/* *fc =dum; */ +/* } */ +#ifdef DEBUG + printf("mnbrak34 fu < or >= fc \n"); + fprintf(ficlog, "mnbrak34 fu < fc\n"); +#endif + dum=u; /* Shifting c and u */ + u = *cx; + *cx = dum; + dum = fu; + fu = *fc; + *fc =dum; +#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"); +#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"); +#endif u=ulim; fu=(*func)(u); - } else { + } 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"); +#endif u=(*cx)+GOLD*(*cx-*bx); fu=(*func)(u); - } + } /* end tests */ SHFT(*ax,*bx,*cx,u) - SHFT(*fa,*fb,*fc,fu) - } + 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); +#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) */ } /*************** linmin ************************/ @@ -1333,6 +1564,8 @@ void linmin(double p[], double xi[], int int j; double xx,xmin,bx,ax; double fx,fb,fa; + + double scale=10., axs, xxs, xxss; /* Scale added for infinity */ ncom=n; pcom=vector(1,n); @@ -1342,18 +1575,57 @@ void linmin(double p[], double xi[], int pcom[j]=p[j]; xicom[j]=xi[j]; } - ax=0.0; - xx=1.0; - mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */ - *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */ + + /* axs=0.0; */ + /* xxss=1; /\* 1 and using scale *\/ */ + xxs=1; + /* do{ */ + ax=0.; + xx= xxs; + mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */ + /* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */ + /* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */ + /* Outputs: fa=f(p(j)) and fx=f(p(j) + xxs * xi(j) ) and f(bx)= f(p(j)+ bx* xi(j)) */ + /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */ + /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */ + /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/ + /* if (fx != fx){ */ + /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */ + /* 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); */ + /* } */ + /* }while(fx != fx); */ + +#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); +#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); #endif +#ifdef DEBUGLINMIN + printf("linmin end "); +#endif for (j=1;j<=n;j++) { - xi[j] *= xmin; - p[j] += xi[j]; + /* printf(" before xi[%d]=%12.8f", j,xi[j]); */ + xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ + /* if(xxs <1.0) */ + /* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */ + p[j] += xi[j]; /* Parameters values are updated accordingly */ } + /* printf("\n"); */ +#ifdef DEBUGLINMIN + printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); + for (j=1;j<=n;j++) { + printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]); + if(j % ncovmodel == 0) + printf("\n"); + } +#endif free_vector(xicom,1,n); free_vector(pcom,1,n); } @@ -1375,6 +1647,7 @@ void powell(double p[], double **xi, int double (*func)(double [])); int i,ibig,j; double del,t,*pt,*ptt,*xit; + double directest; double fp,fptt; double *xits; int niterf, itmp; @@ -1387,7 +1660,7 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) pt[j]=p[j]; rcurr_time = time(NULL); for (*iter=1;;++(*iter)) { - fp=(*fret); + fp=(*fret); /* From former iteration or initial value */ ibig=0; del=0.0; rlast_time=rcurr_time; @@ -1397,7 +1670,7 @@ void powell(double p[], double **xi, int printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ - for (i=1;i<=n;i++) { + for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); @@ -1425,17 +1698,22 @@ void powell(double p[], double **xi, int 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 (j=1;j<=n;j++) xit[j]=xi[j][i]; + for (i=1;i<=n;i++) { /* For each direction i */ + for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */ fptt=(*fret); #ifdef DEBUG printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); #endif - printf("%d",i);fflush(stdout); + printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); - linmin(p,xit,n,fret,func); - if (fabs(fptt-(*fret)) > del) { + 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 */ + 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; } @@ -1454,8 +1732,23 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } /* end i */ - if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { + } /* end loop on each direction i */ + /* 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) */ + 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 */ + /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ + /* decreased of more than 3.84 */ + /* 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 */ + /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */ #ifdef DEBUG int k[2],l; k[0]=1; @@ -1485,31 +1778,35 @@ void powell(double p[], double **xi, int free_vector(ptt,1,n); free_vector(pt,1,n); return; - } + } /* enough precision */ if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { /* Computes an extrapolated point */ + for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j]; pt[j]=p[j]; } - fptt=(*func)(ptt); + fptt=(*func)(ptt); /* f_3 */ +#ifdef POWELLF1F3 +#else if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ +#endif /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ /* From x1 (P0) distance of x2 is at h and x3 is 2h */ /* Let f"(x2) be the 2nd derivative equal everywhere. */ /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */ /* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h**2 */ - /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */ - /* Thus we compare delta(2h) with observed f1-f3 */ - /* or best gain on one ancient line 'del' with total */ - /* gain f1-f2 = f1 - f2 - 'del' with del */ + /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */ /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */ - - t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); +#ifdef NRCORIGINAL + t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/ +#else + t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */ t= t- del*SQR(fp-fptt); - printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); - fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); +#endif + directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */ #ifdef DEBUG + printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); + fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), @@ -1517,14 +1814,39 @@ void powell(double p[], double **xi, int printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); #endif - if (t < 0.0) { /* Then we use it for last direction */ - linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/ +#ifdef POWELLORIGINAL + if (t < 0.0) { /* Then we use it for new direction */ +#else + if (directest*t < 0.0) { /* Contradiction between both tests */ + printf("directest= %.12lf, 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, 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]); + if(j % ncovmodel == 0) + printf("\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++) { - xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ - xi[j][n]=xit[j]; /* and nth direction by the extrapolated */ + printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0) + printf("\n"); } - printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); +#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 DEBUG printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); @@ -1536,9 +1858,12 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } /* end of t negative */ + } /* end of t or directest negative */ +#ifdef POWELLF1F3 +#else } /* end if (fptt < fp) */ - } +#endif + } /* loop iteration */ } /**** Prevalence limit (stable or period prevalence) ****************/ @@ -1567,14 +1892,16 @@ double **prevalim(double **prlim, int nl newm=savm; /* Covariates have to be included here again */ cov[2]=agefin; - + if(nagesqr==1) + cov[3]= agefin*agefin;; for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/ } - /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */ - /* cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; */ + /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2]; + for (k=1; k<=cptcovprod;k++) /* Useless */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ @@ -1728,6 +2055,7 @@ double ***hpxij(double ***po, int nhstep int i, j, d, h, k; double **out, cov[NCOVMAX+1]; double **newm; + double agexact; /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) @@ -1741,13 +2069,17 @@ double ***hpxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + cov[2]=agexact; + if(nagesqr==1) + cov[3]= agexact*agexact; for (k=1; k<=cptcovn;k++) - cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - for (k=1; k<=cptcovage;k++) - cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; + for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ + /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ @@ -1799,6 +2131,7 @@ double func( double *x) int s1, s2; double bbh, survp; long ipmx; + double agexact; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -1820,7 +2153,7 @@ double func( double *x) to be observed in j being in i according to the model. */ for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ - cov[2+k]=covar[Tvar[k]][i]; + cov[2+nagesqr+k]=covar[Tvar[k]][i]; } /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] @@ -1833,9 +2166,12 @@ double func( double *x) } for(d=0; d ncodemax[j]) break; - } /* end of loop on */ - } /* end of loop on modality */ + ij=0; /* ij is similar to i but can jumps over null modalities */ + for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/ + if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */ + break; + } + ij++; + nbcode[Tvar[j]][ij]=i; /* stores the original 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 */ + /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ + /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ + /* } */ + /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ + /* if (ij > ncodemax[j]) { */ + /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ + /* break; */ + /* } */ + /* } /\* end of loop on modality k *\/ */ } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */ + for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ /* 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]++; + Ndum[ij]++; /* Might be supersed V1 + V1*age */ } - ij=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) */ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ if((Ndum[i]!=0) && (i<=ncovcol)){ + ij++; /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ Tvaraff[ij]=i; /*For printing (unclear) */ - ij++; - }else - Tvaraff[ij]=0; + }else{ + /* Tvaraff[ij]=0; */ + } } - ij--; + /* ij--; */ cptcoveff=ij; /*Number of total covariates*/ } @@ -3872,17 +4268,20 @@ To be simple, these graphs help to under gm=vector(1,(nlstate)*(nlstate+ndeath)); for (age=bage; age<=fage; age ++){ cov[2]=age; + if(nagesqr==1) + cov[3]= age*age; for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4 + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4 * 1 1 1 1 1 * 2 2 1 1 1 * 3 1 2 1 1 */ /* nbcode[1][1]=0 nbcode[1][2]=1;*/ } - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; + /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; for(theta=1; theta <=npar; theta++){ @@ -4114,12 +4513,14 @@ fprintf(fichtm," \n
  • Graphs jj1=0; for(k1=1; k1<=m;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ jj1++; if (cptcovn > 0) { fprintf(fichtm,"
    ************ Results for covariates"); - for (cpt=1; cpt<=cptcoveff;cpt++) + for (cpt=1; cpt<=cptcoveff;cpt++){ fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); + printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout); + } fprintf(fichtm," ************\n
    "); } /* Pij */ @@ -4138,16 +4539,16 @@ fprintf(fichtm," \n
    • Graphs fprintf(fichtm,"\n
      - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : %s%d%d.png
      \ ",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); } - } /* end i1 */ + /* } /\* end i1 *\/ */ }/* End k1 */ fprintf(fichtm,"
    "); - fprintf(fichtm,"\ \n

  • Result files (second order: variances)

    \n\ - - Parameter file with estimated parameters and covariance matrix: %s
    \n", rfileres,rfileres); + - Parameter file with estimated parameters and covariance matrix: %s
    \ + - 95%% confidence intervals and T statistics are in the log file.
    \n", rfileres,rfileres); - fprintf(fichtm," - Variance of one-step probabilities: %s
    \n", + fprintf(fichtm," - Standard deviation of one-step probabilities: %s
    \n", subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); fprintf(fichtm,"\ - Variance-covariance of one-step probabilities: %s
    \n", @@ -4188,7 +4589,7 @@ fprintf(fichtm," \n
    • Graphs jj1=0; for(k1=1; k1<=m;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ + /* for(i1=1; i1<=ncodemax[k1];i1++){ */ jj1++; if (cptcovn > 0) { fprintf(fichtm,"
      ************ Results for covariates"); @@ -4207,7 +4608,7 @@ true period expectancies (those weighted drawn in addition to the population based expectancies computed using\ observed and cahotic prevalences: %s%d.png
      \ ",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); - } /* end i1 */ + /* } /\* end i1 *\/ */ }/* End k1 */ fprintf(fichtm,"
    "); fflush(fichtm); @@ -4342,20 +4743,42 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) } /* end covariate */ /* proba elementaires */ + fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); for(i=1,jk=1; i <=nlstate; i++){ + 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]); + fprintf(ficgp,"p%d=%f; ",jk,p[jk]); jk++; - fprintf(ficgp,"\n"); } + fprintf(ficgp,"\n"); } } } + fprintf(ficgp,"##############\n#\n"); + /*goto avoid;*/ + fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n"); + fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); + fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n"); + fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n"); + fprintf(ficgp,"# logi(p13/p11)=p6 +p7*age +p8*age*age+ p9*V1+ p10*V1*age\n"); + fprintf(ficgp,"# p12+p13+p14+p11=1=p11(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); + fprintf(ficgp,"# p11=1/(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); + fprintf(ficgp,"# p12=exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)/\n"); + fprintf(ficgp,"# (1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); + fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n"); + fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n"); + fprintf(ficgp,"#\n"); for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/ + fprintf(ficgp,"# ng=%d\n",ng); + fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m); for(jk=1; jk <=m; jk++) { + fprintf(ficgp,"# jk=%d\n",jk); fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); if (ng==2) fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n"); @@ -4368,30 +4791,40 @@ plot [%.f:%.f] ", ageminpar, agemaxpar) for(k=1; k<=(nlstate+ndeath); k++) { if (k != k2){ if(ng==2) - fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); + if(nagesqr==0) + fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); + else /* nagesqr =1 */ + fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr); else - fprintf(ficgp," exp(p%d+p%d*x",i,i+1); + if(nagesqr==0) + fprintf(ficgp," exp(p%d+p%d*x",i,i+1); + else /* nagesqr =1 */ + fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); ij=1;/* To be checked else nbcode[0][0] wrong */ - for(j=3; j <=ncovmodel; j++) { - /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */ - /* /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */ - /* ij++; */ - /* } */ - /* else */ - fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); + for(j=3; j <=ncovmodel-nagesqr; j++) { + if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */ + fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); + ij++; + } + else + fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); } fprintf(ficgp,")/(1"); - for(k1=1; k1 <=nlstate; k1++){ - fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); + for(k1=1; k1 <=nlstate; k1++){ + if(nagesqr==0) + fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); + else /* nagesqr =1 */ + fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr); + ij=1; - for(j=3; j <=ncovmodel; j++){ - /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */ - /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */ - /* ij++; */ - /* } */ - /* else */ - fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); + for(j=3; j <=ncovmodel-nagesqr; j++){ + if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { + fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); + ij++; + } + else + fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]); } fprintf(ficgp,")"); } @@ -5185,9 +5618,10 @@ void removespace(char *str) { } 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 - * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 - * - cptcovn or number of covariates k of the models excluding age*products =6 + * 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 @@ -5202,22 +5636,14 @@ int decodemodel ( char model[], int last int j1, k1, k2; char modelsav[80]; char stra[80], strb[80], strc[80], strd[80],stre[80]; + char *strpt; /*removespace(model);*/ if (strlen(model) >1){ /* If there is at least 1 covariate */ j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0; - j=nbocc(model,'+'); /**< j=Number of '+' */ - j1=nbocc(model,'*'); /**< j1=Number of '*' */ - cptcovs=j+1-j1; /**< Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2 */ - cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 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 */ - strcpy(modelsav,model); if (strstr(model,"AGE") !=0){ - printf("Error. AGE must be in lower case 'age' model=%s ",model); - fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog); + printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model); + fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog); return 1; } if (strstr(model,"v") !=0){ @@ -5225,117 +5651,153 @@ int decodemodel ( char model[], int last fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog); return 1; } - - /* Design - * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight - * < ncovcol=8 > - * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 - * 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 - * 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 - * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11 - * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product - * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8 - * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2]; - * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted - * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 - * < ncovcol=8 > - * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 - * k= 1 2 3 4 5 6 7 8 9 10 11 12 - * Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 - * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} - * p Tprod[1]@2={ 6, 5} - *p Tvard[1][1]@4= {7, 8, 5, 6} - * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 - * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; - *How to reorganize? - * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age - * Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} - * {2, 1, 4, 8, 5, 6, 3, 7} - * Struct [] - */ + strcpy(modelsav,model); + 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 \ + '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 \ + '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; + } - /* This loop fills the array Tvar from the string 'model'.*/ - /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ - /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */ - /* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */ - /* k=3 V4 Tvar[k=3]= 4 (from V4) */ - /* k=2 V1 Tvar[k=2]= 1 (from V1) */ - /* k=1 Tvar[1]=2 (from V2) */ - /* k=5 Tvar[5] */ - /* for (k=1; k<=cptcovn;k++) { */ - /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ - /* } */ - /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - /* - * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ - for(k=cptcovt; k>=1;k--) /**< Number of covariates */ + nagesqr=1; + if (strstr(model,"+age*age") !=0) + substrchaine(modelsav, model, "+age*age"); + else if (strstr(model,"age*age+") !=0) + substrchaine(modelsav, model, "age*age+"); + else + 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 */ + 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). */ + cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */ + cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ + + + /* Design + * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight + * < ncovcol=8 > + * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + * 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 + * 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 + * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11 + * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product + * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8 + * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2]; + * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted + * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 + * < ncovcol=8 > + * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 + * k= 1 2 3 4 5 6 7 8 9 10 11 12 + * Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 + * p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} + * p Tprod[1]@2={ 6, 5} + *p Tvard[1][1]@4= {7, 8, 5, 6} + * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 + * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; + *How to reorganize? + * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age + * Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} + * {2, 1, 4, 8, 5, 6, 3, 7} + * Struct [] + */ + + /* This loop fills the array Tvar from the string 'model'.*/ + /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ + /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */ + /* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */ + /* k=3 V4 Tvar[k=3]= 4 (from V4) */ + /* k=2 V1 Tvar[k=2]= 1 (from V1) */ + /* k=1 Tvar[1]=2 (from V2) */ + /* k=5 Tvar[5] */ + /* for (k=1; k<=cptcovn;k++) { */ + /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ + /* } */ + /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */ + /* + * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ + for(k=cptcovt; k>=1;k--) /**< Number of covariates */ 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 */ - cptcovage++; /* Sums the number of covariates which include age as a product */ - Tage[cptcovage]=k; /* Tage[1] = 4 */ - /*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 */ + 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++; - 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 + */ - } /* end model */ + 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) */ /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ @@ -5451,47 +5913,138 @@ int calandcheckages(int imx, int maxwav, return (1); } -void syscompilerinfo() +#if defined(_MSC_VER) +/*printf("Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ +/*fprintf(ficlog, "Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ +//#include "stdafx.h" +//#include +//#include +//#include +//#include +typedef BOOL(WINAPI *LPFN_ISWOW64PROCESS) (HANDLE, PBOOL); + +LPFN_ISWOW64PROCESS fnIsWow64Process; + +BOOL IsWow64() +{ + BOOL bIsWow64 = FALSE; + + //typedef BOOL (APIENTRY *LPFN_ISWOW64PROCESS) + // (HANDLE, PBOOL); + + //LPFN_ISWOW64PROCESS fnIsWow64Process; + + HMODULE module = GetModuleHandle(_T("kernel32")); + const char funcName[] = "IsWow64Process"; + fnIsWow64Process = (LPFN_ISWOW64PROCESS) + GetProcAddress(module, funcName); + + if (NULL != fnIsWow64Process) + { + if (!fnIsWow64Process(GetCurrentProcess(), + &bIsWow64)) + //throw std::exception("Unknown error"); + printf("Unknown error\n"); + } + return bIsWow64 != FALSE; +} +#endif + +void syscompilerinfo(int logged) { /* #include "syscompilerinfo.h"*/ - /* #include */ /* Only on gnu */ + /* command line Intel compiler 32bit windows, XP compatible:*/ + /* /GS /W3 /Gy + /Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D + "_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D + "UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo + /Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch" + */ + /* 64 bits */ + /* + /GS /W3 /Gy + /Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" /D "NDEBUG" + /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo /Zc:forScope + /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Fo"x64\Release\" /Qprof-dir + "x64\Release\" /Fp"x64\Release\IMaCh.pch" */ + /* Optimization are useless and O3 is slower than O2 */ + /* + /GS /W3 /Gy /Zc:wchar_t /Zi /O3 /Fd"x64\Release\vc120.pdb" /D "WIN32" + /D "NDEBUG" /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo + /Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel + /Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" + */ + /* Link is */ /* /OUT:"visual studio + 2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT + /PDB:"visual studio + 2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE + "kernel32.lib" "user32.lib" "gdi32.lib" "winspool.lib" + "comdlg32.lib" "advapi32.lib" "shell32.lib" "ole32.lib" + "oleaut32.lib" "uuid.lib" "odbc32.lib" "odbccp32.lib" + /MACHINE:X86 /OPT:REF /SAFESEH /INCREMENTAL:NO + /SUBSYSTEM:CONSOLE",5.01" /MANIFESTUAC:"level='asInvoker' + uiAccess='false'" + /ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF + /NOLOGO /TLBID:1 + */ +#if defined __INTEL_COMPILER +#if defined(__GNUC__) + struct utsname sysInfo; /* For Intel on Linux and OS/X */ +#endif +#elif defined(__GNUC__) +#ifndef __APPLE__ +#include /* Only on gnu */ +#endif + struct utsname sysInfo; + int cross = CROSS; + if (cross){ + printf("Cross-"); + if(logged) fprintf(ficlog, "Cross-"); + } +#endif + #include - printf("Compiled with:");fprintf(ficlog,"Compiled with:"); + + printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:"); #if defined(__clang__) - printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ + printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ #endif #if defined(__ICC) || defined(__INTEL_COMPILER) - printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */ + printf(" Intel ICC/ICPC");if(logged)fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */ #endif #if defined(__GNUC__) || defined(__GNUG__) - printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */ + printf(" GNU GCC/G++");if(logged)fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */ #endif #if defined(__HP_cc) || defined(__HP_aCC) - printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */ + printf(" Hewlett-Packard C/aC++");if(logged)fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */ #endif #if defined(__IBMC__) || defined(__IBMCPP__) - printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */ + printf(" IBM XL C/C++"); if(logged) fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */ #endif #if defined(_MSC_VER) - printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */ + printf(" Microsoft Visual Studio");if(logged)fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */ #endif #if defined(__PGI) - printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */ + printf(" Portland Group PGCC/PGCPP");if(logged) fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */ #endif #if defined(__SUNPRO_C) || defined(__SUNPRO_CC) - printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */ + printf(" Oracle Solaris Studio");if(logged)fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */ #endif - printf(". ");fprintf(ficlog,". "); + printf(" for "); if (logged) fprintf(ficlog, " for "); // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros #ifdef _WIN32 // note the underscore: without it, it's not msdn official! // Windows (x64 and x86) + printf("Windows (x64 and x86) ");if(logged) fprintf(ficlog,"Windows (x64 and x86) "); #elif __unix__ // all unices, not all compilers // Unix + printf("Unix ");if(logged) fprintf(ficlog,"Unix "); #elif __linux__ // linux + printf("linux ");if(logged) fprintf(ficlog,"linux "); #elif __APPLE__ - // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though... + // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though.. + printf("Mac OS ");if(logged) fprintf(ficlog,"Mac OS "); #endif /* __MINGW32__ */ @@ -5505,22 +6058,13 @@ void syscompilerinfo() /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */ #if UINTPTR_MAX == 0xffffffff - printf(" 32-bit.\n"); fprintf(ficlog," 32-bit.\n");/* 32-bit */ + printf(" 32-bit"); if(logged) fprintf(ficlog," 32-bit");/* 32-bit */ #elif UINTPTR_MAX == 0xffffffffffffffff - printf(" 64-bit.\n"); fprintf(ficlog," 64-bit.\n");/* 64-bit */ + printf(" 64-bit"); if(logged) fprintf(ficlog," 64-bit");/* 64-bit */ #else - printf(" wtf-bit.\n"); fprintf(ficlog," wtf-bit.\n");/* wtf */ + printf(" wtf-bit"); if(logged) fprintf(ficlog," wtf-bit");/* wtf */ #endif -/* struct utsname sysInfo; - - if (uname(&sysInfo) != -1) { - printf(" %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); - fprintf(ficlog," %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); - } - else - perror("uname() error"); - */ #if defined(__GNUC__) # if defined(__GNUC_PATCHLEVEL__) # define __GNUC_VERSION__ (__GNUC__ * 10000 \ @@ -5530,18 +6074,181 @@ void syscompilerinfo() # define __GNUC_VERSION__ (__GNUC__ * 10000 \ + __GNUC_MINOR__ * 100) # endif - printf("GNU C version %d.\n", __GNUC_VERSION__); - fprintf(ficlog, "GNU C version %d.\n", __GNUC_VERSION__); + printf(" using GNU C version %d.\n", __GNUC_VERSION__); + if(logged) fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__); + + if (uname(&sysInfo) != -1) { + printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); + if(logged) fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); + } + else + perror("uname() error"); + //#ifndef __INTEL_COMPILER +#if !defined (__INTEL_COMPILER) && !defined(__APPLE__) + printf("GNU libc version: %s\n", gnu_get_libc_version()); + if(logged) fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version()); #endif +#endif + + // void main() + // { #if defined(_MSC_VER) - /*printf("Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ - /*fprintf(ficlog, "Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ + if (IsWow64()){ + printf("\nThe program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); + if (logged) fprintf(ficlog, "\nThe program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); + } + else{ + printf("\nThe program is not running under WOW64 (i.e probably on a 64bit Windows).\n"); + if (logged) fprintf(ficlog, "\nThe programm is not running under WOW64 (i.e probably on a 64bit Windows).\n"); + } + // printf("\nPress Enter to continue..."); + // getchar(); + // } + #endif - /* printf("GNU libc version: %s\n", gnu_get_libc_version()); */ } +int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){ + /*--------------- Prevalence limit (period or stable prevalence) --------------*/ + int i, j, k, i1 ; + double ftolpl = 1.e-10; + double age, agebase, agelim; + + strcpy(filerespl,"pl"); + strcat(filerespl,fileres); + if((ficrespl=fopen(filerespl,"w"))==NULL) { + printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; + fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; + } + printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); + pstamp(ficrespl); + fprintf(ficrespl,"# Period (stable) prevalence \n"); + fprintf(ficrespl,"#Age "); + for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); + fprintf(ficrespl,"\n"); + + /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */ + + agebase=ageminpar; + agelim=agemaxpar; + + i1=pow(2,cptcoveff); + if (cptcovn < 1){i1=1;} + + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */ + //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + /* to clean */ + //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); + fprintf(ficrespl,"\n#******"); + printf("\n#******"); + fprintf(ficlog,"\n#******"); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + fprintf(ficrespl,"******\n"); + printf("******\n"); + fprintf(ficlog,"******\n"); + + fprintf(ficrespl,"#Age "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); + fprintf(ficrespl,"\n"); + + for (age=agebase; age<=agelim; age++){ + /* for (age=agebase; age<=agebase; age++){ */ + prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + fprintf(ficrespl,"%.0f ",age ); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + for(i=1; i<=nlstate;i++) + fprintf(ficrespl," %.5f", prlim[i][i]); + fprintf(ficrespl,"\n"); + } /* Age */ + /* was end of cptcod */ + } /* cptcov */ + return 0; +} + +int hPijx(double *p, int bage, int fage){ + /*------------- h Pij x at various ages ------------*/ + + int stepsize; + int agelim; + int hstepm; + int nhstepm; + int h, i, i1, j, k; + + double agedeb; + double ***p3mat; + + strcpy(filerespij,"pij"); strcat(filerespij,fileres); + if((ficrespij=fopen(filerespij,"w"))==NULL) { + printf("Problem with Pij resultfile: %s\n", filerespij); return 1; + fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1; + } + printf("Computing pij: result on file '%s' \n", filerespij); + fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); + + stepsize=(int) (stepm+YEARM-1)/YEARM; + /*if (stepm<=24) stepsize=2;*/ + + agelim=AGESUP; + hstepm=stepsize*YEARM; /* Every year of age */ + hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ + + /* hstepm=1; aff par mois*/ + pstamp(ficrespij); + fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); + i1= pow(2,cptcoveff); + /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ + /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */ + /* k=k+1; */ + for (k=1; k <= (int) pow(2,cptcoveff); k++){ + fprintf(ficrespij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespij,"******\n"); + + for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + + /* nhstepm=nhstepm*YEARM; aff par mois*/ + + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + oldm=oldms;savm=savms; + hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate+ndeath;j++) + fprintf(ficrespij," %1d-%1d",i,j); + fprintf(ficrespij,"\n"); + for (h=0; h<=nhstepm; h++){ + /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/ + fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate+ndeath;j++) + fprintf(ficrespij," %.5f", p3mat[i][j][h]); + fprintf(ficrespij,"\n"); + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + fprintf(ficrespij,"\n"); + } + /*}*/ + } + return 0; +} + + /***********************************************/ /**************** Main Program *****************/ /***********************************************/ @@ -5568,11 +6275,12 @@ int main(int argc, char *argv[]) /* FILE *fichtm; *//* Html File */ /* FILE *ficgp;*/ /*Gnuplot File */ struct stat info; - double agedeb; - double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; + double agedeb=0.; + + double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; double fret; - double dum; /* Dummy variable */ + double dum=0.; /* Dummy variable */ double ***p3mat; double ***mobaverage; @@ -5582,18 +6290,20 @@ int main(int argc, char *argv[]) char *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; int c, h , cpt; - int jl; - int i1, j1, jk, stepsize; + int jl=0; + int i1, j1, jk, stepsize=0; + int count=0; + int *tab; int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; - int hstepm, nhstepm; + int hstepm=0, nhstepm=0; int agemortsup; float sumlpop=0.; double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000; - double bage=0, fage=110, age, agelim, agebase; + double bage=0, fage=110., age, agelim=0., agebase=0.; double ftolpl=FTOL; double **prlim; double ***param; /* Matrix of parameters */ @@ -5649,8 +6359,12 @@ int main(int argc, char *argv[]) nberr=0; /* Number of errors and warnings */ nbwarn=0; +#ifdef WIN32 + _getcwd(pathcd, size); +#else getcwd(pathcd, size); - +#endif + syscompilerinfo(0); printf("\n%s\n%s",version,fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); @@ -5685,9 +6399,14 @@ int main(int argc, char *argv[]) /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); +#ifdef WIN32 + _chdir(path); /* Can be a relative path */ + if(_getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ +#else chdir(path); /* Can be a relative path */ - if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ - printf("Current directory %s!\n",pathcd); + if (getcwd(pathcd, MAXLINE) > 0) /* So pathcd is the full path */ +#endif + printf("Current directory %s!\n",pathcd); strcpy(command,"mkdir "); strcat(command,optionfilefiname); if((outcmd=system(command)) != 0){ @@ -5702,7 +6421,7 @@ int main(int argc, char *argv[]) /*-------- arguments in the command line --------*/ - /* Log file */ + /* Main Log file */ strcat(filelog, optionfilefiname); strcat(filelog,".log"); /* */ if((ficlog=fopen(filelog,"w"))==NULL) { @@ -5718,7 +6437,7 @@ int main(int argc, char *argv[]) optionfilext=%s\n\ optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); - syscompilerinfo(); + syscompilerinfo(0); printf("Local time (at start):%s",strstart); fprintf(ficlog,"Local time (at start): %s",strstart); @@ -5731,7 +6450,7 @@ int main(int argc, char *argv[]) strcat(fileres, optionfilefiname); strcat(fileres,".txt"); /* Other files have txt extension */ - /*---------arguments file --------*/ + /* Main ---------arguments file --------*/ if((ficpar=fopen(optionfile,"r"))==NULL) { printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); @@ -5764,12 +6483,24 @@ int main(int argc, char *argv[]) } ungetc(c,ficpar); - fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); - numlinepar++; - printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); - fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); - fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); + numlinepar=numlinepar+3; /* In general */ + printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); + if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */ + model[strlen(model)-1]='\0'; + fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); fflush(ficlog); + /* if(model[0]=='#'|| model[0]== '\0'){ */ + if(model[0]=='#'){ + printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \ + 'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \ + 'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \ + if(mle != -1){ + printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n"); + exit(1); + } + } while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); @@ -5788,10 +6519,9 @@ int main(int argc, char *argv[]) v1+v2*age+v2*v3 makes cptcovn = 3 */ if (strlen(model)>1) - ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ + ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7,age*age makes 3*/ else - ncovmodel=2; - nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ + ncovmodel=2; /* Constant and age */ nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ npar= nforce*ncovmodel; /* Number of parameters like aij*/ if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ @@ -5806,18 +6536,18 @@ int main(int argc, char *argv[]) /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ if(mle==-1){ /* Print a wizard for help writing covariance matrix */ prwizard(ncovmodel, nlstate, ndeath, model, ficparo); - printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); - fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); + printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); + fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); fclose (ficparo); fclose (ficlog); goto end; exit(0); } - else if(mle==-3) { + else if(mle==-3) { /* Main Wizard */ prwizard(ncovmodel, nlstate, ndeath, model, ficparo); - printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); - fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + fprintf(ficlog," You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso); param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); matcov=matrix(1,npar,1,npar); } @@ -5841,7 +6571,7 @@ int main(int argc, char *argv[]) if(jj==i) continue; j++; fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1 != i) && (j1 != j)){ + if ((i1 != i) || (j1 != jj)){ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ It might be a problem of design; if ncovcol and the model are correct\n \ run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1); @@ -5849,8 +6579,8 @@ run imach with mle=-1 to get a correct t } fprintf(ficparo,"%1d%1d",i1,j1); if(mle==1) - printf("%1d%1d",i,j); - fprintf(ficlog,"%1d%1d",i,j); + printf("%1d%1d",i,jj); + fprintf(ficlog,"%1d%1d",i,jj); for(k=1; k<=ncovmodel;k++){ fscanf(ficpar," %lf",¶m[i][j][k]); if(mle==1){ @@ -5931,12 +6661,22 @@ run imach with mle=-1 to get a correct t for(i=1; i <=npar; i++) for(j=1; j <=npar; j++) matcov[i][j]=0.; + /* Scans npar lines */ for(i=1; i <=npar; i++){ - fscanf(ficpar,"%s",str); + count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk); + if(count != 3){ + printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + exit(1); + }else if(mle==1) - printf("%s",str); - fprintf(ficlog,"%s",str); - fprintf(ficparo,"%s",str); + printf("%1d%1d%1d",i1,j1,jk); + fprintf(ficlog,"%1d%1d%1d",i1,j1,jk); + fprintf(ficparo,"%1d%1d%1d",i1,j1,jk); for(j=1; j <=i; j++){ fscanf(ficpar," %le",&matcov[i][j]); if(mle==1){ @@ -5952,6 +6692,7 @@ run imach with mle=-1 to get a correct t fprintf(ficlog,"\n"); fprintf(ficparo,"\n"); } + /* End of read covariance matrix npar lines */ for(i=1; i <=npar; i++) for(j=i+1;j<=npar;j++) matcov[i][j]=matcov[j][i]; @@ -5974,7 +6715,8 @@ run imach with mle=-1 to get a correct t fprintf(ficres,"#%s\n",version); } /* End of mle != -3 */ - + /* Main data + */ n= lastobs; num=lvector(1,n); moisnais=vector(1,n); @@ -5990,6 +6732,7 @@ run imach with mle=-1 to get a correct t s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ + ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ /* Reads data from file datafile */ if (readdata(datafile, firstobs, lastobs, &imx)==1) @@ -6025,6 +6768,9 @@ run imach with mle=-1 to get a correct t Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 */ +/* Main decodemodel */ + + if(decodemodel(model, lastobs) == 1) goto end; @@ -6068,11 +6814,15 @@ run imach with mle=-1 to get a correct t nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; Ndum =ivector(-1,NCOVMAX); - if (ncovmodel > 2) + if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ + /* Nbcode gives the value of the lth modality of jth covariate, in + V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ + /* 1 to ncodemax[j] is the maximum value of this jth covariate */ codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/ + /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ h=0; @@ -6089,7 +6839,9 @@ run imach with mle=-1 to get a correct t if (h>m) h=1; /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 - * h 1 2 3 4 + * For k=4 covariates, h goes from 1 to 2**k + * codtabm(h,k)= 1 & (h-1) >> (k-1) ; + * h\k 1 2 3 4 *______________________________ * 1 i=1 1 i=1 1 i=1 1 i=1 1 * 2 2 1 1 1 @@ -6109,6 +6861,7 @@ run imach with mle=-1 to get a correct t * 16 2 2 2 1 */ codtab[h][k]=j; + /* codtab[12][3]=1; */ /*codtab[h][Tvar[k]]=j;*/ printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); } @@ -6129,7 +6882,7 @@ run imach with mle=-1 to get a correct t - /*------------ gnuplot -------------*/ + /* Initialisation of ----------- gnuplot -------------*/ strcpy(optionfilegnuplot,optionfilefiname); if(mle==-3) strcat(optionfilegnuplot,"-mort"); @@ -6145,7 +6898,9 @@ run imach with mle=-1 to get a correct t fprintf(ficgp,"set datafile missing 'NaNq'\n"); } /* fclose(ficgp);*/ - /*--------- index.htm --------*/ + + + /* Initialisation of --------- index.htm --------*/ strcpy(optionfilehtm,optionfilefiname); /* Main html file */ if(mle==-3) @@ -6187,7 +6942,12 @@ Title=%s
    Datafile=%s Firstpass=%d La strcpy(pathr,path); strcat(pathr,optionfilefiname); +#ifdef WIN32 + _chdir(optionfilefiname); /* Move to directory named optionfile */ +#else chdir(optionfilefiname); /* Move to directory named optionfile */ +#endif + /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -6210,10 +6970,10 @@ Interval (in months) between two waves: p=param[1][1]; /* *(*(*(param +1)+1)+0) */ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ - + /* For mortality only */ if (mle==-3){ ximort=matrix(1,NDIM,1,NDIM); -/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ + /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ cens=ivector(1,n); ageexmed=vector(1,n); agecens=vector(1,n); @@ -6303,8 +7063,8 @@ Interval (in months) between two waves: /* Initialize method and iterate */ /* p[1]=0.0268; p[NDIM]=0.083; */ -/* gsl_vector_set(x, 0, 0.0268); */ -/* gsl_vector_set(x, 1, 0.083); */ + /* gsl_vector_set(x, 0, 0.0268); */ + /* gsl_vector_set(x, 1, 0.083); */ gsl_vector_set(x, 0, p[1]); gsl_vector_set(x, 1, p[2]); @@ -6372,9 +7132,10 @@ Interval (in months) between two waves: } printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); - for (i=1;i<=NDIM;i++) + for (i=1;i<=NDIM;i++) { printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); - + fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + } lsurv=vector(1,AGESUP); lpop=vector(1,AGESUP); tpop=vector(1,AGESUP); @@ -6406,8 +7167,15 @@ Interval (in months) between two waves: replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); - + if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + }else + printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ stepm, weightopt,\ model,imx,p,matcov,agemortsup); @@ -6421,27 +7189,28 @@ Interval (in months) between two waves: free_ivector(dcwave,1,n); free_matrix(ximort,1,NDIM,1,NDIM); #endif - } /* Endof if mle==-3 */ - + } /* Endof if mle==-3 mortality only */ + /* Standard maximisation */ else{ /* For mle >=1 */ globpr=0;/* debug */ + /* Computes likelihood for initial parameters */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); printf("\n"); - globpr=1; /* to print the contributions */ + globpr=1; /* again, to print the contributions */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); printf("\n"); - if(mle>=1){ /* Could be 1 or 2 */ + if(mle>=1){ /* Could be 1 or 2, Real Maximisation */ mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); } /*--------- results files --------------*/ - fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); + fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); @@ -6454,9 +7223,9 @@ Interval (in months) between two waves: fprintf(ficlog,"%d%d ",i,k); fprintf(ficres,"%1d%1d ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%lf ",p[jk]); - fprintf(ficlog,"%lf ",p[jk]); - fprintf(ficres,"%lf ",p[jk]); + printf("%12.7f ",p[jk]); + fprintf(ficlog,"%12.7f ",p[jk]); + fprintf(ficres,"%12.7f ",p[jk]); jk++; } printf("\n"); @@ -6470,6 +7239,24 @@ Interval (in months) between two waves: ftolhess=ftol; /* Usually correct */ hesscov(matcov, p, npar, delti, ftolhess, func); } + printf("Parameters and 95%% confidence intervals\n"); + fprintf(ficlog, "Parameters, T and confidence intervals\n"); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); + fprintf(ficlog,"%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + } + } + } + fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); printf("# Scales (for hessian or gradient estimation)\n"); fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); @@ -6597,7 +7384,8 @@ Interval (in months) between two waves: fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); - + + /* Other stuffs, more or less useful */ while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); @@ -6625,6 +7413,7 @@ Interval (in months) between two waves: dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; fscanf(ficpar,"pop_based=%d\n",&popbased); + fprintf(ficlog,"pop_based=%d\n",popbased); fprintf(ficparo,"pop_based=%d\n",popbased); fprintf(ficres,"pop_based=%d\n",popbased); @@ -6649,7 +7438,15 @@ Interval (in months) between two waves: /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + }else + printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ @@ -6670,8 +7467,13 @@ Interval (in months) between two waves: fclose(ficres); + /* Other results (useful)*/ + + /*--------------- Prevalence limit (period or stable prevalence) --------------*/ -#include "prevlim.h" /* Use ficrespl, ficlog */ + /*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ + prlim=matrix(1,nlstate,1,nlstate); + prevalence_limit(p, prlim, ageminpar, agemaxpar); fclose(ficrespl); #ifdef FREEEXIT2 @@ -6679,7 +7481,8 @@ Interval (in months) between two waves: #endif /*------------- h Pij x at various ages ------------*/ -#include "hpijx.h" + /*#include "hpijx.h"*/ + hPijx(p, bage, fage); fclose(ficrespij); /*-------------- Variance of one-step probabilities---*/ @@ -6706,7 +7509,8 @@ Interval (in months) between two waves: /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ /* } */ } - + + /* ------ Other prevalence ratios------------ */ /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ @@ -6937,6 +7741,7 @@ Interval (in months) between two waves: free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); free_ivector(ncodemax,1,NCOVMAX); + free_ivector(ncodemaxwundef,1,NCOVMAX); free_ivector(Tvar,1,NCOVMAX); free_ivector(Tprod,1,NCOVMAX); free_ivector(Tvaraff,1,NCOVMAX); @@ -6981,9 +7786,15 @@ Interval (in months) between two waves: printf("Before Current directory %s!\n",pathcd); +#ifdef WIN32 + if (_chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(_getcwd(pathcd,MAXLINE) > 0) +#else if(chdir(pathcd) != 0) - printf("Can't move to directory %s!\n",path); - if(getcwd(pathcd,MAXLINE) > 0) + printf("Can't move to directory %s!\n", path); + if (getcwd(pathcd, MAXLINE) > 0) +#endif printf("Current directory %s!\n",pathcd); /*strcat(plotcmd,CHARSEPARATOR);*/ sprintf(plotcmd,"gnuplot");