--- imach/src/imach.c 2015/01/03 16:33:42 1.175 +++ imach/src/imach.c 2015/10/27 17:36:57 1.207 @@ -1,6 +1,130 @@ -/* $Id: imach.c,v 1.175 2015/01/03 16:33:42 brouard Exp $ +/* $Id: imach.c,v 1.207 2015/10/27 17:36:57 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.207 2015/10/27 17:36:57 brouard + *** empty log message *** + + Revision 1.206 2015/10/24 07:14:11 brouard + *** empty log message *** + + Revision 1.205 2015/10/23 15:50:53 brouard + Summary: 0.98r3 some clarification for graphs on likelihood contributions + + Revision 1.204 2015/10/01 16:20:26 brouard + Summary: Some new graphs of contribution to likelihood + + Revision 1.203 2015/09/30 17:45:14 brouard + Summary: looking at better estimation of the hessian + + Also a better criteria for convergence to the period prevalence And + therefore adding the number of years needed to converge. (The + prevalence in any alive state shold sum to one + + Revision 1.202 2015/09/22 19:45:16 brouard + Summary: Adding some overall graph on contribution to likelihood. Might change + + Revision 1.201 2015/09/15 17:34:58 brouard + Summary: 0.98r0 + + - Some new graphs like suvival functions + - Some bugs fixed like model=1+age+V2. + + Revision 1.200 2015/09/09 16:53:55 brouard + Summary: Big bug thanks to Flavia + + Even model=1+age+V2. did not work anymore + + Revision 1.199 2015/09/07 14:09:23 brouard + Summary: 0.98q6 changing default small png format for graph to vectorized svg. + + Revision 1.198 2015/09/03 07:14:39 brouard + Summary: 0.98q5 Flavia + + Revision 1.197 2015/09/01 18:24:39 brouard + *** empty log message *** + + Revision 1.196 2015/08/18 23:17:52 brouard + Summary: 0.98q5 + + Revision 1.195 2015/08/18 16:28:39 brouard + Summary: Adding a hack for testing purpose + + After reading the title, ftol and model lines, if the comment line has + a q, starting with #q, the answer at the end of the run is quit. It + permits to run test files in batch with ctest. The former workaround was + $ echo q | imach foo.imach + + 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 *** @@ -540,7 +664,16 @@ end */ +/* #define DEBUG */ +/* #define DEBUGBRENT */ +/* #define DEBUGLINMIN */ +/* #define DEBUGHESS */ +#define DEBUGHESSIJ +/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ #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 @@ -607,11 +740,12 @@ typedef struct { #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ -#define codtabm(h,k) 1 & (h-1) >> (k-1) ; +#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 #define MAXN 20000 #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 '\\' @@ -623,15 +757,16 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.175 2015/01/03 16:33:42 brouard Exp $ */ +/* $Id: imach.c,v 1.207 2015/10/27 17:36:57 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.175 $ $Date: 2015/01/03 16:33:42 $"; +#include "version.h" +char version[]=__IMACH_VERSION__; +char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; +char fullversion[]="$Revision: 1.207 $ $Date: 2015/10/27 17:36:57 $"; 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) */ @@ -693,7 +828,7 @@ char command[FILENAMELENGTH]; int outcmd=0; char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; - +char fileresu[FILENAMELENGTH]; /* fileres without r in front */ char filelog[FILENAMELENGTH]; /* Log file */ char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; @@ -742,7 +877,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 *Tage; int *Ndum; /** Freq of modality (tricode */ -int **codtab; /**< codtab=imatrix(1,100,1,10); */ +/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ int **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; @@ -790,7 +932,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 ); @@ -801,18 +943,22 @@ 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*/ printf(" DIRC = %s \n",dirc); - } else { /* strip direcotry from path */ + } else { /* strip directory from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ 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 */ @@ -864,11 +1010,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; @@ -894,7 +1083,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 @@ -1205,7 +1394,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; @@ -1258,19 +1453,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; @@ -1281,25 +1476,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) */ @@ -1308,23 +1522,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 ************************/ @@ -1347,27 +1623,97 @@ void linmin(double p[], double xi[], int int j; double xx,xmin,bx,ax; double fx,fb,fa; - + +#ifdef LINMINORIGINAL +#else + double scale=10., axs, xxs; /* Scale added for infinity */ +#endif + ncom=n; pcom=vector(1,n); xicom=vector(1,n); nrfunc=func; for (j=1;j<=n;j++) { pcom[j]=p[j]; - xicom[j]=xi[j]; + xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */ } - 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 */ + +#ifdef LINMINORIGINAL + xx=1.; +#else + axs=0.0; + xxs=1.; + do{ + xx= xxs; +#endif + ax=0.; + 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]]*/ +#ifdef LINMINORIGINAL +#else + if (fx != fx){ + xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ + printf("|"); + fprintf(ficlog,"|"); +#ifdef DEBUGLINMIN + printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); +#endif + } + }while(fx != fx); +#endif + +#ifdef DEBUGLINMIN + printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); + fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); +#endif + *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 "); + fprintf(ficlog,"linmin end "); +#endif for (j=1;j<=n;j++) { +#ifdef LINMINORIGINAL xi[j] *= xmin; - p[j] += xi[j]; +#else +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" before xi[%d]=%12.8f", j,xi[j]); +#endif + xi[j] *= xmin*xxs; /* xi rescaled by xmin and number of loops: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ +#ifdef DEBUGLINMIN + 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 ); +#endif +#endif + p[j] += xi[j]; /* Parameters values are updated accordingly */ } +#ifdef DEBUGLINMIN + printf("\n"); + printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); + fprintf(ficlog,"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]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]); + fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } +#else +#endif free_vector(xicom,1,n); free_vector(pcom,1,n); } @@ -1389,6 +1735,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; @@ -1399,9 +1746,9 @@ void powell(double p[], double **xi, int xits=vector(1,n); *fret=(*func)(p); for (j=1;j<=n;j++) pt[j]=p[j]; - rcurr_time = time(NULL); + 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; @@ -1411,7 +1758,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]); @@ -1439,17 +1786,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); + 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; } @@ -1468,8 +1820,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; @@ -1499,31 +1866,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 delta 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), @@ -1531,14 +1902,45 @@ 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 (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); + printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); + fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); + fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); + } + if (directest < 0.0) { /* Then we use it for new direction */ +#endif +#ifdef DEBUGLINMIN + printf("Before linmin in direction P%d-P0\n",n); + for (j=1;j<=n;j++) { + printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } + } +#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]); + fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ + printf("\n"); + fprintf(ficlog,"\n"); + } } - 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); @@ -1550,24 +1952,45 @@ 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) ****************/ -double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) +double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij) { /* Computes the prevalence limit in each live state at age x by left multiplying the unit - matrix by transitions matrix until convergence is reached */ - + matrix by transitions matrix until convergence is reached with precision ftolpl */ + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + int i, ii,j,k; double min, max, maxmin, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **pmij(); double **newm; - double agefin, delaymax=50 ; /* Max number of years to converge */ + double agefin, delaymax=100 ; /* Max number of years to converge */ + int ncvloop=0; for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ @@ -1577,18 +2000,25 @@ double **prevalim(double **prlim, int nl cov[1]=1.; /* Even if hstepm = 1, at least one multiplication by the unit matrix */ + /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ + ncvloop++; 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]]]; - /*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]]]; */ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ + } + /*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]][codtabm(ij,Tvar[k])]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; + for (k=1; k<=cptcovprod;k++) /* Useless */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; /*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]);*/ @@ -1607,17 +2037,23 @@ double **prevalim(double **prlim, int nl sumnew=0; for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; prlim[i][j]= newm[i][j]/(1-sumnew); - /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/ max=FMAX(max,prlim[i][j]); min=FMIN(min,prlim[i][j]); + printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); } - maxmin=max-min; + maxmin=(max-min)/(max+min)*2; maxmax=FMAX(maxmax,maxmin); } /* j loop */ + *ncvyear= (int)age- (int)agefin; + printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); if(maxmax < ftolpl){ + /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ return prlim; } } /* age loop */ + printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\ +Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); +/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ return prlim; /* should not reach here */ } @@ -1742,6 +2178,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++) @@ -1755,13 +2192,20 @@ 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]][codtabm(ij,k)]; + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(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]]][codtabm(ij,k)]*cov[2]; + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(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]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ @@ -1813,6 +2257,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]);*/ @@ -1834,7 +2279,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] @@ -1847,9 +2292,12 @@ double func( double *x) } for(d=0; dFile of contributions to the likelihood: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fflush(fichtm); - } + if (mle ==0) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); + else if(mle >=1) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); + fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); + + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fflush(fichtm); + + for (k=1; k<= nlstate ; k++) { + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j %s-p%dj.png
\ +",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); + } + } return; } @@ -2225,7 +2727,7 @@ void mlikeli(FILE *ficres,double p[], in for (j=1;j<=npar;j++) xi[i][j]=(i==j ? 1.0 : 0.0); printf("Powell\n"); fprintf(ficlog,"Powell\n"); - strcpy(filerespow,"pow"); + strcpy(filerespow,"POW_"); strcat(filerespow,fileres); if((ficrespow=fopen(filerespow,"w"))==NULL) { printf("Problem with resultfile: %s\n", filerespow); @@ -2270,30 +2772,30 @@ void mlikeli(FILE *ficres,double p[], in fclose(ficrespow); printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); - fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); } /**** Computes Hessian and covariance matrix ***/ -void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) +void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) { double **a,**y,*x,pd; - double **hess; + /* double **hess; */ int i, j; int *indx; double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); - double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar); + double hessij(double p[], double **hess, double delti[], int i, int j,double (*func)(double []),int npar); void lubksb(double **a, int npar, int *indx, double b[]) ; void ludcmp(double **a, int npar, int *indx, double *d) ; double gompertz(double p[]); - hess=matrix(1,npar,1,npar); + /* hess=matrix(1,npar,1,npar); */ printf("\nCalculation of the hessian matrix. Wait...\n"); fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n"); for (i=1;i<=npar;i++){ - printf("%d",i);fflush(stdout); - fprintf(ficlog,"%d",i);fflush(ficlog); + printf("%d-",i);fflush(stdout); + fprintf(ficlog,"%d-",i);fflush(ficlog); hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); @@ -2304,9 +2806,9 @@ void hesscov(double **matcov, double p[] for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { if (j>i) { - printf(".%d%d",i,j);fflush(stdout); - fprintf(ficlog,".%d%d",i,j);fflush(ficlog); - hess[i][j]=hessij(p,delti,i,j,func,npar); + printf(".%d-%d",i,j);fflush(stdout); + fprintf(ficlog,".%d-%d",i,j);fflush(ficlog); + hess[i][j]=hessij(p,hess, delti,i,j,func,npar); hess[j][i]=hess[i][j]; /*printf(" %lf ",hess[i][j]);*/ @@ -2340,53 +2842,78 @@ void hesscov(double **matcov, double p[] fprintf(ficlog,"\n#Hessian matrix#\n"); for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { - printf("%.3e ",hess[i][j]); - fprintf(ficlog,"%.3e ",hess[i][j]); + printf("%.6e ",hess[i][j]); + fprintf(ficlog,"%.6e ",hess[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } + /* printf("\n#Covariance matrix#\n"); */ + /* fprintf(ficlog,"\n#Covariance matrix#\n"); */ + /* for (i=1;i<=npar;i++) { */ + /* for (j=1;j<=npar;j++) { */ + /* printf("%.6e ",matcov[i][j]); */ + /* fprintf(ficlog,"%.6e ",matcov[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + /* Recompute Inverse */ - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; - ludcmp(a,npar,indx,&pd); + /* for (i=1;i<=npar;i++) */ + /* for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */ + /* ludcmp(a,npar,indx,&pd); */ + + /* printf("\n#Hessian matrix recomputed#\n"); */ + + /* for (j=1;j<=npar;j++) { */ + /* for (i=1;i<=npar;i++) x[i]=0; */ + /* x[j]=1; */ + /* lubksb(a,npar,indx,x); */ + /* for (i=1;i<=npar;i++){ */ + /* y[i][j]=x[i]; */ + /* printf("%.3e ",y[i][j]); */ + /* fprintf(ficlog,"%.3e ",y[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + + /* Verifying the inverse matrix */ +#ifdef DEBUGHESS + y=matprod2(y,hess,1,npar,1,npar,1,npar,matcov); - /* printf("\n#Hessian matrix recomputed#\n"); + printf("\n#Verification: multiplying the matrix of covariance by the Hessian matrix, should be unity:#\n"); + fprintf(ficlog,"\n#Verification: multiplying the matrix of covariance by the Hessian matrix. Should be unity:#\n"); for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); for (i=1;i<=npar;i++){ - y[i][j]=x[i]; - printf("%.3e ",y[i][j]); - fprintf(ficlog,"%.3e ",y[i][j]); + printf("%.2f ",y[i][j]); + fprintf(ficlog,"%.2f ",y[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } - */ +#endif free_matrix(a,1,npar,1,npar); free_matrix(y,1,npar,1,npar); free_vector(x,1,npar); free_ivector(indx,1,npar); - free_matrix(hess,1,npar,1,npar); + /* free_matrix(hess,1,npar,1,npar); */ } /*************** hessian matrix ****************/ double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) -{ +{ /* Around values of x, computes the function func and returns the scales delti and hessian */ int i; int l=1, lmax=20; - double k1,k2; + double k1,k2, res, fx; double p2[MAXPARM+1]; /* identical to x */ - double res; double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; - double fx; int k=0,kmax=10; double l1; @@ -2402,9 +2929,9 @@ double hessii(double x[], double delta, p2[theta]=x[theta]-delt; k2=func(p2)-fx; /*res= (k1-2.0*fx+k2)/delt/delt; */ - res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ + res= (k1+k2)/delt/delt/2.; /* Divided by 2 because L and not 2*L */ -#ifdef DEBUGHESS +#ifdef DEBUGHESSII printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); #endif @@ -2418,48 +2945,122 @@ double hessii(double x[], double delta, else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ delts=delt; } - } + } /* End loop k */ } delti[theta]=delts; return res; } -double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) +double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; int l=1, lmax=20; double k1,k2,k3,k4,res,fx; double p2[MAXPARM+1]; - int k; - + int k, kmax=1; + double v1, v2, cv12, lc1, lc2; + fx=func(x); - for (k=1; k<=2; k++) { + for (k=1; k<=kmax; k=k+10) { for (i=1;i<=npar;i++) p2[i]=x[i]; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k1=func(p2)-fx; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k2=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k3=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k4=func(p2)-fx; - res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */ -#ifdef DEBUG - printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); - fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ + if(k1*k2*k3*k4 <0.){ + kmax=kmax+10; + if(kmax >=10){ + printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } + } +#ifdef DEBUGHESSIJ + v1=hess[thetai][thetai]; + v2=hess[thetaj][thetaj]; + cv12=res; + /* Computing eigen value of Hessian matrix */ + lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } #endif } return res; } + /* Not done yet: Was supposed to fix if not exactly at the maximum */ +/* double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) */ +/* { */ +/* int i; */ +/* int l=1, lmax=20; */ +/* double k1,k2,k3,k4,res,fx; */ +/* double p2[MAXPARM+1]; */ +/* double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; */ +/* int k=0,kmax=10; */ +/* double l1; */ + +/* fx=func(x); */ +/* for(l=0 ; l <=lmax; l++){ /\* Enlarging the zone around the Maximum *\/ */ +/* l1=pow(10,l); */ +/* delts=delt; */ +/* for(k=1 ; k khi/nkhif) || (k2 >khi/nkhif) || (k4 >khi/nkhif) || (k4 >khi/nkhif)){ /\* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. *\/ */ +/* k=kmax; l=lmax*10; */ +/* } */ +/* else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ */ +/* delts=delt; */ +/* } */ +/* } /\* End loop k *\/ */ +/* } */ +/* delti[theta]=delts; */ +/* return res; */ +/* } */ + + /************** Inverse of matrix **************/ void ludcmp(double **a, int n, int *indx, double *d) { @@ -2536,7 +3137,7 @@ void lubksb(double **a, int n, int *indx void pstamp(FILE *fichier) { - fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart); + fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); } /************ Frequencies ********************/ @@ -2552,8 +3153,8 @@ void freqsummary(char fileres[], int ia pp=vector(1,nlstate); prop=matrix(1,nlstate,iagemin,iagemax+3); - strcpy(fileresp,"p"); - strcat(fileresp,fileres); + strcpy(fileresp,"P_"); + strcat(fileresp,fileresu); if((ficresp=fopen(fileresp,"w"))==NULL) { printf("Problem with prevalence resultfile: %s\n", fileresp); fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); @@ -2588,13 +3189,13 @@ void freqsummary(char fileres[], int ia bool=1; if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){ + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* Tests if the value of each of the covariates of i is equal to filter j1 */ bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1], - j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/ + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), + j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ + /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } } @@ -2623,10 +3224,10 @@ void freqsummary(char fileres[], int ia pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresp, "**********\n#"); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficlog, "**********\n#"); } for(i=1; i<=nlstate;i++) @@ -2754,7 +3355,7 @@ void prevalence(double ***probs, double bool=1; if (cptcovn>0) { for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) bool=0; } if (bool==1) { @@ -2959,12 +3560,12 @@ void tricode(int *Tvar, int **nbcode, in cptcoveff=0; - for (k=-1; k < maxncov; k++) Ndum[k]=0; for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ /* Loop on covariates without age and products */ - for (j=1; j<=(cptcovs); j++) { /* model V1 + V2*age+ V3 + V3*V4 : V1 + V3 = 2 only */ - for (i=1; i<=imx; i++) { /* Lopp on individuals: reads the data file to get the maximum value of the + for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */ + for (k=-1; k < maxncov; k++) Ndum[k]=0; + for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i * If product of Vn*Vm, still boolean *: @@ -2986,66 +3587,87 @@ void tricode(int *Tvar, int **nbcode, in /* getting the maximum value of the modality of the covariate (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and female is 1, then modmaxcovj=1.*/ - } + } /* end for loop on individuals i */ printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); + fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); cptcode=modmaxcovj; /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ /*for (i=0; i<=cptcode; i++) {*/ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */ - printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]); - if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */ - ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */ + for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */ + printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); + fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); + if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */ + if( k != -1){ + ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered excluding + undefined. Usually 2: 0 and 1. */ + } + ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th + covariate for which somebody answered including + undefined. Usually 3: -1, 0 and 1. */ } /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ } /* Ndum[-1] number of undefined modalities */ /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ - /* For covariate j, modalities could be 1, 2, 3, 4. If Ndum[2]=0 ncodemax[j] is not 4 but 3 */ - /* If Ndum[3}= 635; Ndum[4]=0; Ndum[5]=0; Ndum[6]=27; Ndum[7]=125; + /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. + If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; modmincovj=3; modmaxcovj = 7; - There are only 3 modalities non empty (or 2 if 27 is too few) : ncodemax[j]=3; - which will be coded 0, 1, 2 which in binary on 3-1 digits are 0=00 1=01, 2=10; defining two dummy - variables V1_1 and V1_2. + There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; + which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; + defining two dummy variables: variables V1_1 and V1_2. nbcode[Tvar[j]][ij]=k; nbcode[Tvar[j]][1]=0; nbcode[Tvar[j]][2]=1; nbcode[Tvar[j]][3]=2; + To be continued (not working yet). */ - ij=1; /* 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 */ - for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */ - /*recode from 0 */ - if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ - nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode. - 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; */ - ij++; - } - if (ij > ncodemax[j]) break; - } /* end of loop on */ - } /* end of loop on modality */ + ij=0; /* ij is similar to i but can jump over null modalities */ + for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + if (Ndum[i] == 0) { /* If nobody responded to this modality k */ + break; + } + ij++; + nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ + cptcode = ij; /* New max modality for covar j */ + } /* 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*/ } @@ -3371,7 +3993,7 @@ void cvevsij(double ***eij, double x[], } /************ Variance ******************/ -void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) + void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -3400,11 +4022,11 @@ void varevsij(char optionfilefiname[], d if(popbased==1){ if(mobilav!=0) - strcpy(digitp,"-populbased-mobilav-"); - else strcpy(digitp,"-populbased-nomobil-"); + strcpy(digitp,"-POPULBASED-MOBILAV_"); + else strcpy(digitp,"-POPULBASED-NOMOBIL_"); } else - strcpy(digitp,"-stablbased-"); + strcpy(digitp,"-STABLBASED_"); if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); @@ -3414,12 +4036,12 @@ void varevsij(char optionfilefiname[], d } } - strcpy(fileresprobmorprev,"prmorprev"); + strcpy(fileresprobmorprev,"PRMORPREV-"); sprintf(digit,"%-d",ij); /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/ strcat(fileresprobmorprev,digit); /* Tvar to be done */ strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */ - strcat(fileresprobmorprev,fileres); + strcat(fileresprobmorprev,fileresu); if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobmorprev); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); @@ -3437,7 +4059,8 @@ void varevsij(char optionfilefiname[], d } fprintf(ficresprobmorprev,"\n"); fprintf(ficgp,"\n# Routine varevsij"); - /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ + fprintf(ficgp,"\nunset title \n"); +/* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); /* } */ @@ -3496,7 +4119,7 @@ void varevsij(char optionfilefiname[], d xp[i] = x[i] + (i==theta ?delti[theta]:0); } hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); if (popbased==1) { if(mobilav ==0){ @@ -3527,7 +4150,7 @@ void varevsij(char optionfilefiname[], d for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij); if (popbased==1) { if(mobilav ==0){ @@ -3602,7 +4225,7 @@ void varevsij(char optionfilefiname[], d /* end ppptj */ /* x centered again */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij); if (popbased==1) { if(mobilav ==0){ @@ -3649,9 +4272,11 @@ void varevsij(char optionfilefiname[], d free_vector(gmp,nlstate+1,nlstate+ndeath); free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); + /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */ + fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ @@ -3659,11 +4284,11 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev)); fprintf(fichtm,"\n
    File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); - fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); - /* fprintf(fichtm,"\n
    Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

    \n", stepm,YEARM,digitp,digit); + fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); + /* fprintf(fichtm,"\n
    Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

    \n", stepm,YEARM,digitp,digit); */ -/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */ - fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); +/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */ + fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,nlstate); @@ -3678,9 +4303,9 @@ void varevsij(char optionfilefiname[], d } /* end varevsij */ /************ Variance of prevlim ******************/ -void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[]) + void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[]) { - /* Variance of prevalence limit */ + /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ double **dnewm,**doldm; @@ -3717,13 +4342,13 @@ void varprevlim(char fileres[], double * for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); for(i=1;i<=nlstate;i++) gp[i] = prlim[i][i]; for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); for(i=1;i<=nlstate;i++) gm[i] = prlim[i][i]; @@ -3739,8 +4364,13 @@ void varprevlim(char fileres[], double * for(i=1;i<=nlstate;i++) varpl[i][(int)age] =0.; + if((int)age==67 ||(int)age== 66 ){ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + }else{ + matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ @@ -3781,20 +4411,20 @@ void varprob(char optionfilefiname[], do char fileresprobcor[FILENAMELENGTH]; double ***varpij; - strcpy(fileresprob,"prob"); + strcpy(fileresprob,"PROB_"); strcat(fileresprob,fileres); if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); } - strcpy(fileresprobcov,"probcov"); - strcat(fileresprobcov,fileres); + strcpy(fileresprobcov,"PROBCOV_"); + strcat(fileresprobcov,fileresu); if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobcov); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); } - strcpy(fileresprobcor,"probcor"); - strcat(fileresprobcor,fileres); + strcpy(fileresprobcor,"PROBCOR_"); + strcat(fileresprobcor,fileresu); if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobcor); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); @@ -3836,10 +4466,9 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of pairs of step probabilities (drawings)

    • \n",optionfilehtmcov); - fprintf(fichtmcov,"\n

    Matrix of variance-covariance of pairs of step probabilities

    \n\ - file %s
    \n",optionfilehtmcov); - fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated\ + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.
    • \n",optionfilehtmcov); + fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ They are expressed in year-1 in order to be less dependent of stepm.
    \n"); fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ @@ -3860,23 +4489,23 @@ To be simple, these graphs help to under /*j1++;*/ if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprob, "**********\n#\n"); fprintf(ficresprobcov, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprobcov, "**********\n#\n"); fprintf(ficgp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficgp, "**********\n#\n"); fprintf(fichtmcov, "\n
    ********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(fichtmcov, "**********\n
    "); fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprobcor, "**********\n#"); } @@ -3886,17 +4515,21 @@ 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]][codtabm(j1,k)]; + /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(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]]][codtabm(ij,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]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; for(theta=1; theta <=npar; theta++){ @@ -4044,17 +4677,18 @@ To be simple, these graphs help to under /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ if(first==1){ first=0; + fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); fprintf(ficgp,"\nset parametric;unset label"); fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); - fprintf(ficgp,"\nset ter png small size 320, 240"); + fprintf(ficgp,"\nset ter svg size 640, 480"); fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ - :\ -%s%d%1d%1d-%1d%1d.png, ",k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + :\ +%s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\ + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,\ + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); fprintf(fichtmcov,"\n
    Correlation at age %d (%.3f),",(int) age, c12); - fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ @@ -4071,7 +4705,7 @@ To be simple, these graphs help to under }/* if first */ } /* age mod 5 */ } /* end loop age */ - fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); first=1; } /*l12 */ } /* k12 */ @@ -4093,7 +4727,7 @@ To be simple, these graphs help to under /******************* Printing html file ***********/ -void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ +void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ int popforecast, int estepm ,\ @@ -4106,20 +4740,20 @@ void printinghtml(char fileres[], char t "); fprintf(fichtm,"