From: N. Brouard Date: Wed, 24 Apr 2024 21:21:17 +0000 (+0000) Subject: Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegen... X-Git-Tag: imach-099s7~28 X-Git-Url: https://henry.ined.fr/git/?a=commitdiff_plain;h=30f643bce1dcbf6de5ff3b84e4109866411c1a66;p=.git Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes --- diff --git a/src/imach.c b/src/imach.c index ca7e4e9..2ee93a2 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,22 +1,20 @@ /* $Id$ $State$ $Log$ - Revision 1.357 2023/06/14 14:55:52 brouard - * imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE + Revision 1.6 2024/04/24 21:10:29 brouard + Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes - Revision 1.356 2023/05/23 12:08:43 brouard - Summary: 0.99r46 + Revision 1.5 2023/10/09 09:10:01 brouard + Summary: trying to reconsider - * imach.c (Module): Fixed PROB_r + Revision 1.4 2023/06/22 12:50:51 brouard + Summary: stil on going - Revision 1.355 2023/05/22 17:03:18 brouard - Summary: 0.99r46 - - * imach.c (Module): In the ILK....txt file, the number of columns - before the covariates values is dependent of the number of states (16+nlstate): 0.99r46 + Revision 1.3 2023/06/22 11:28:07 brouard + *** empty log message *** - Revision 1.354 2023/05/21 05:05:17 brouard - Summary: Temporary change for imachprax + Revision 1.2 2023/06/22 11:22:40 brouard + Summary: with svd but not working yet Revision 1.353 2023/05/08 18:48:22 brouard *** empty log message *** @@ -1295,7 +1293,8 @@ Important routines /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ /* #define FLATSUP *//* Suppresses directions where likelihood is flat */ -#define POWELLORIGINCONJUGATE /* Don't use conjugate but biggest decrease if valuable */ +/* #define POWELLORIGINCONJUGATE /\* Don't use conjugate but biggest decrease if valuable *\/ */ +/* #define NOTMINFIT */ #include #include @@ -1392,7 +1391,7 @@ double gnuplotversion=GNUPLOTVERSION; /* $State$ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="Testing conjugate April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; +char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; char fullversion[]="$Revision$ $Date$"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; @@ -1440,7 +1439,8 @@ int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav=0; /* Maxim number of waves */ int jmin=0, jmax=0; /* min, max spacing between 2 waves */ int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */ -int gipmx=0, gsw=0; /* Global variables on the number of contributions +int gipmx = 0; +double gsw = 0; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ int mle=1, weightopt=0; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ @@ -2627,6 +2627,1512 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [ free_vector(pcom,1,n); } +/**** praxis gegen ****/ + +/* This has been tested by Visual C from Microsoft and works */ +/* meaning tha valgrind could be wrong */ +/*********************************************************************/ +/* f u n c t i o n p r a x i s */ +/* */ +/* praxis is a general purpose routine for the minimization of a */ +/* function in several variables. the algorithm used is a modifi- */ +/* cation of conjugate gradient search method by powell. the changes */ +/* are due to r.p. brent, who gives an algol-w program, which served */ +/* as a basis for this function. */ +/* */ +/* references: */ +/* - powell, m.j.d., 1964. an efficient method for finding */ +/* the minimum of a function in several variables without */ +/* calculating derivatives, computer journal, 7, 155-162 */ +/* - brent, r.p., 1973. algorithms for minimization without */ +/* derivatives, prentice hall, englewood cliffs. */ +/* */ +/* problems, suggestions or improvements are always wellcome */ +/* karl gegenfurtner 07/08/87 */ +/* c - version */ +/*********************************************************************/ +/* */ +/* usage: min = praxis(tol, macheps, h, n, prin, x, func) */ +/* macheps has been suppressed because it is replaced by DBL_EPSILON */ +/* and if it was an argument of praxis (as it is in original brent) */ +/* it should be declared external */ +/* usage: min = praxis(tol, h, n, prin, x, func) */ +/* was min = praxis(fun, x, n); */ +/* */ +/* fun the function to be minimized. fun is called from */ +/* praxis with x and n as arguments */ +/* x a double array containing the initial guesses for */ +/* the minimum, which will contain the solution on */ +/* return */ +/* n an integer specifying the number of unknown */ +/* parameters */ +/* min praxis returns the least calculated value of fun */ +/* */ +/* some additional global variables control some more aspects of */ +/* the inner workings of praxis. setting them is optional, they */ +/* are all set to some reasonable default values given below. */ +/* */ +/* prin controls the printed output from the routine. */ +/* 0 -> no output */ +/* 1 -> print only starting and final values */ +/* 2 -> detailed map of the minimization process */ +/* 3 -> print also eigenvalues and vectors of the */ +/* search directions */ +/* the default value is 1 */ +/* tol is the tolerance allowed for the precision of the */ +/* solution. praxis returns if the criterion */ +/* 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol */ +/* is fulfilled more than ktm times. */ +/* the default value depends on the machine precision */ +/* ktm see just above. default is 1, and a value of 4 leads */ +/* to a very(!) cautious stopping criterion. */ +/* h0 or step is a steplength parameter and should be set equal */ +/* to the expected distance from the solution. */ +/* exceptionally small or large values of step lead to */ +/* slower convergence on the first few iterations */ +/* the default value for step is 1.0 */ +/* scbd is a scaling parameter. 1.0 is the default and */ +/* indicates no scaling. if the scales for the different */ +/* parameters are very different, scbd should be set to */ +/* a value of about 10.0. */ +/* illc should be set to true (1) if the problem is known to */ +/* be ill-conditioned. the default is false (0). this */ +/* variable is automatically set, when praxis finds */ +/* the problem to be ill-conditioned during iterations. */ +/* maxfun is the maximum number of calls to fun allowed. praxis */ +/* will return after maxfun calls to fun even when the */ +/* minimum is not yet found. the default value of 0 */ +/* indicates no limit on the number of calls. */ +/* this return condition is only checked every n */ +/* iterations. */ +/* */ +/*********************************************************************/ + +#include +#include +#include +#include /* for DBL_EPSILON */ +/* #include "machine.h" */ + + +/* extern void minfit(int n, double eps, double tol, double **ab, double q[]); */ +/* extern void minfit(int n, double eps, double tol, double ab[N][N], double q[]); */ +/* control parameters */ +/* control parameters */ +#define SQREPSILON 1.0e-19 +/* #define EPSILON 1.0e-8 */ /* in main */ + +double tol = SQREPSILON, + scbd = 1.0, + step = 1.0; +int ktm = 1, + /* prin = 2, */ + maxfun = 0, + illc = 0; + +/* some global variables */ +static int i, j, k, k2, nl, nf, kl, kt; +/* static double s; */ +double sl, dn, dmin, + fx, f1, lds, ldt, sf, df, + qf1, qd0, qd1, qa, qb, qc, + m2, m4, small_windows, vsmall, large, + vlarge, ldfac, t2; +/* static double d[N], y[N], z[N], */ +/* q0[N], q1[N], v[N][N]; */ + +static double *d, *y, *z; +static double *q0, *q1, **v; +double *tflin; /* used in flin: return (*fun)(tflin, n); */ +double *e; /* used in minfit, don't konw how to free memory and thus made global */ +/* static double s, sl, dn, dmin, */ +/* fx, f1, lds, ldt, sf, df, */ +/* qf1, qd0, qd1, qa, qb, qc, */ +/* m2, m4, small, vsmall, large, */ +/* vlarge, ldfac, t2; */ +/* static double d[N], y[N], z[N], */ +/* q0[N], q1[N], v[N][N]; */ + +/* these will be set by praxis to point to it's arguments */ +static int prin; /* added */ +static int n; +static double *x; +static double (*fun)(); +/* static double (*fun)(double *x, int n); */ + +/* these will be set by praxis to the global control parameters */ +/* static double h, macheps, t; */ +extern double macheps; +static double h; +static double t; + +static double +drandom() /* return random no between 0 and 1 */ +{ + return (double)(rand()%(8192*2))/(double)(8192*2); +} + +static void sort() /* d and v in descending order */ +{ + int k, i, j; + double s; + + for (i=1; i<=n-1; i++) { + k = i; s = d[i]; + for (j=i+1; j<=n; j++) { + if (d[j] > s) { + k = j; + s = d[j]; + } + } + if (k > i) { + d[k] = d[i]; + d[i] = s; + for (j=1; j<=n; j++) { + s = v[j][i]; + v[j][i] = v[j][k]; + v[j][k] = s; + } + } + } +} + +double randbrent ( int *naught ) +{ + double ran1, ran3[127], half; + int ran2, q, r, i, j; + int init=0; /* false */ + double rr; + /* REAL*8 RAN1,RAN3(127),HALF */ + + /* INTEGER RAN2,Q,R */ + /* LOGICAL INIT */ + /* DATA INIT/.FALSE./ */ + /* IF (INIT) GO TO 3 */ + if(!init){ +/* R = MOD(NAUGHT,8190) + 1 *//* 1804289383 rand () */ + r = *naught % 8190 + 1;/* printf(" naught r %d %d",*naught,r); */ + ran2=127; + for(i=ran2; i>0; i--){ +/* RAN2 = 128 */ +/* DO 2 I=1,127 */ + ran2 = ran2-1; +/* RAN2 = RAN2 - 1 */ + ran1 = -pow(2.0,55); +/* RAN1 = -2.D0**55 */ +/* DO 1 J=1,7 */ + for(j=1; j<=7;j++){ +/* R = MOD(1756*R,8191) */ + r = (1756*r) % 8191;/* printf(" i=%d (1756*r)%8191=%d",j,r); */ + q=r/32; +/* Q = R/32 */ +/* 1 RAN1 = (RAN1 + Q)*(1.0D0/256) */ + ran1 =(ran1+q)*(1.0/256); + } +/* 2 RAN3(RAN2) = RAN1 */ + ran3[ran2] = ran1; /* printf(" ran2=%d ran1=%.7g \n",ran2,ran1); */ + } +/* INIT = .TRUE. */ + init=1; +/* 3 IF (RAN2.EQ.1) RAN2 = 128 */ + } + if(ran2 == 0) ran2 = 126; + else ran2 = ran2 -1; + /* RAN2 = RAN2 - 1 */ + /* RAN1 = RAN1 + RAN3(RAN2) */ + ran1 = ran1 + ran3[ran2];/* printf("BIS ran2=%d ran1=%.7g \n",ran2,ran1); */ + half= 0.5; + /* HALF = .5D0 */ + /* IF (RAN1.GE.0.D0) HALF = -HALF */ + if(ran1 >= 0.) half =-half; + ran1 = ran1 +half; + ran3[ran2] = ran1; + rr= ran1+0.5; + /* RAN1 = RAN1 + HALF */ + /* RAN3(RAN2) = RAN1 */ + /* RANDOM = RAN1 + .5D0 */ +/* r = ( ( double ) ( *seed ) ) * 4.656612875E-10; */ + return rr; +} +static void matprint(char *s, double **v, int m, int n) +/* char *s; */ +/* double v[N][N]; */ +{ +#define INCX 8 + int i; + + int i2hi; + int ihi; + int ilo; + int i2lo; + int jlo=1; + int j; + int j2hi; + int jhi; + int j2lo; + ilo=1; + ihi=n; + jlo=1; + jhi=n; + + printf ("\n" ); + printf ("%s\n", s ); + for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) + { + j2hi = j2lo + INCX - 1; + if ( n < j2hi ) + { + j2hi = n; + } + if ( jhi < j2hi ) + { + j2hi = jhi; + } + + /* fprintf ( ficlog, "\n" ); */ + printf ("\n" ); +/* + For each column J in the current range... + + Write the header. +*/ + /* fprintf ( ficlog, " Col: "); */ + printf ("Col:"); + for ( j = j2lo; j <= j2hi; j++ ) + { + /* fprintf ( ficlog, " %7d ", j - 1 ); */ + /* printf (" %9d ", j - 1 ); */ + printf (" %9d ", j ); + } + /* fprintf ( ficlog, "\n" ); */ + /* fprintf ( ficlog, " Row\n" ); */ + /* fprintf ( ficlog, "\n" ); */ + printf ("\n" ); + printf (" Row\n" ); + printf ("\n" ); +/* + Determine the range of the rows in this strip. +*/ + if ( 1 < ilo ){ + i2lo = ilo; + }else{ + i2lo = 1; + } + if ( m < ihi ){ + i2hi = m; + }else{ + i2hi = ihi; + } + + for ( i = i2lo; i <= i2hi; i++ ){ +/* + Print out (up to) 5 entries in row I, that lie in the current strip. +*/ + /* fprintf ( ficlog, "%5d:", i - 1 ); */ + /* printf ("%5d:", i - 1 ); */ + printf ("%5d:", i ); + for ( j = j2lo; j <= j2hi; j++ ) + { + /* fprintf ( ficlog, " %14g", a[i-1+(j-1)*m] ); */ + /* printf ("%14.7g ", a[i-1+(j-1)*m] ); */ + /* printf("%14.7f ", v[i-1][j-1]); */ + printf("%14.7f ", v[i][j]); + /* fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); */ + } + /* fprintf ( ficlog, "\n" ); */ + printf ("\n" ); + } + } + + /* printf("%s\n", s); */ + /* for (k=0; k 0) { /* linear search */ + /* for (i=0; i F1. NITS CONTROLS THE NUMBER OF TIMES */ + /* AN ATTEMPT IS MADE TO HALVE THE INTERVAL. */ + /* SIDE EFFECTS: USES AND ALTERS X, FX, NF, NL. */ + /* IF J < 1 USES VARIABLES Q... . */ + /* USES H, N, T, M2, M4, LDT, DMIN, MACHEPS; */ + int k, i, dz; + double x2, xm, f0, f2, fm, d1, t2, sf1, sx1; + double s; + double macheps; + macheps=pow(16.0,-13.0); + sf1 = f1; sx1 = *x1; + k = 0; xm = 0.0; fm = f0 = fx; dz = *d2 < macheps; + /* h=1.0;*/ /* To be revised */ +#ifdef DEBUGPRAX + /* printf("min macheps=%14g h=%14g step=%14g t=%14g fx=%14g\n",macheps,h, step,t, fx); */ + /* Where is fx coming from */ + printf(" min macheps=%14g h=%14g t=%14g fx=%.9lf dirj=%d\n",macheps, h, t, fx, j); + matprint(" min vectors:",v,n,n); +#endif + /* find step size */ + s = 0.; + /* for (i=0; i s) t2 = s; + if (t2 < small_windows) t2 = small_windows; + if (t2 > 0.01*h) t2 = 0.01 * h; + if (fk && f1 <= fm) { + xm = *x1; + fm = f1; + } +#ifdef DEBUGPRAX + printf(" additional flin X1=%14.7f t2=%14.7f *f1=%14.7f fm=%14.7f fk=%d\n",*x1,t2,f1,fm,fk); +#endif + if (!fk || fabs(*x1) < t2) { + *x1 = (*x1 >= 0 ? t2 : -t2); + /* *x1 = (*x1 > 0 ? t2 : -t2); */ /* kind of error */ +#ifdef DEBUGPRAX + printf(" additional flin X1=%16.10e dirj=%d fk=%d\n",*x1, j, fk); +#endif + f1 = flin(*x1, j); +#ifdef DEBUGPRAX + printf(" after flin f1=%18.12e dirj=%d fk=%d\n",f1, j,fk); +#endif + } + if (f1 <= fm) { + xm = *x1; + fm = f1; + } +L0: /*L0 loop or next */ +/* + Evaluate FLIN at another point and estimate the second derivative. +*/ + if (dz) { + x2 = (f0 < f1 ? -(*x1) : 2*(*x1)); +#ifdef DEBUGPRAX + printf(" additional second flin x2=%14.8e x1=%14.8e f0=%14.8e f1=%18.12e dirj=%d\n",x2,*x1,f0,f1,j); +#endif + f2 = flin(x2, j); +#ifdef DEBUGPRAX + printf(" additional second flin x2=%16.10e x1=%16.10e f1=%18.12e f0=%18.10e f2=%18.10e fm=%18.10e\n",x2, *x1, f1,f0,f2,fm); +#endif + if (f2 <= fm) { + xm = x2; + fm = f2; + } + /* d2 is the curvature or double difference f1 doesn't seem to be accurately computed */ + *d2 = (x2*(f1-f0) - (*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2)); +#ifdef DEBUGPRAX + double d11,d12; + d11=(f1-f0)/(*x1);d12=(f2-f0)/x2; + printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/(x2-(*x1))=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/(x2-(*x1))); + printf(" original computing f1=%18.12e *d2=%16.10e f0=%18.12e f1-f0=%16.10e f2-f0=%16.10e\n",f1,*d2,f0,f1-f0, f2-f0); + double ff1=7.783920622852e+04; + double f1mf0=9.0344736236e-05; + *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); + /* *d2 = (ff1-f0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */ + printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0); + *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2); + printf(" overlifi computing *d2=%16.10e\n",*d2); +#endif + *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2); + } +#ifdef DEBUGPRAX + printf(" additional second flin xm=%14.8e fm=%14.8e *d2=%14.8e\n",xm, fm,*d2); +#endif + /* + Estimate the first derivative at 0. + */ + d1 = (f1-f0)/(*x1) - *x1**d2; dz = 1; + /* + Predict the minimum. + */ + if (*d2 <= small_windows) { + x2 = (d1 < 0 ? h : -h); + } + else { + x2 = - 0.5*d1/(*d2); + } +#ifdef DEBUGPRAX + printf(" AT d1=%14.8e d2=%14.8e small=%14.8e dz=%d x1=%14.8e x2=%14.8e\n",d1,*d2,small_windows,dz,*x1,x2); +#endif + if (fabs(x2) > h) + x2 = (x2 > 0 ? h : -h); +L1: /* L1 or try loop */ +#ifdef DEBUGPRAX + printf(" AT predicted minimum flin x2=%14.8e x1=%14.8e K=%14d NITS=%14d dirj=%d\n",x2,*x1,k,nits,j); +#endif + f2 = flin(x2, j); /* x[i]+x2*v[i][j] */ +#ifdef DEBUGPRAX + printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e\n",f0,f1,f2, fm); +#endif + if ((k < nits) && (f2 > f0)) { +#ifdef DEBUGPRAX + printf(" NO SUCCESS SO TRY AGAIN;\n"); +#endif + k++; + if ((f0 < f1) && (*x1*x2 > 0.0)) + goto L0; /* or next */ + x2 *= 0.5; + goto L1; + } + nl++; +#ifdef DEBUGPRAX + printf(" bebeBE end of min x1=%14.8e x2=%14.8e f1=%14.8e f2=%14.8e f0=%14.8e fm=%14.8e d2=%14.8e\n",*x1, x2, f1, f2, f0, fm, *d2); +#endif + if (f2 > fm) x2 = xm; else fm = f2; + if (fabs(x2*(x2-*x1)) > small_windows) { + *d2 = (x2*(f1-f0) - *x1*(fm-f0))/(*x1*x2*(*x1-x2)); + } + else { + if (k > 0) *d2 = 0; + } +#ifdef DEBUGPRAX + printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2); +#endif + if (*d2 <= small_windows) *d2 = small_windows; + *x1 = x2; fx = fm; + if (sf1 < fx) { + fx = sf1; + *x1 = sx1; + } + /* + Update X for linear search. + */ +#ifdef DEBUGPRAX + printf(" end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2); +#endif + + /* if (j != -1) */ + /* for (i=0; i 0) + for (i=1; i<=n; i++) + x[i] += (*x1)*v[i][j]; +} + +void quad() /* look for a minimum along the curve q0, q1, q2 */ +{ + int i; + double l, s; + + s = fx; fx = qf1; qf1 = s; qd1 = 0.0; + /* for (i=0; i0.0 && qd1>0.0 &&nl>=3*n*n) { +#ifdef DEBUGPRAX + printf(" QUAD before min value=%14.8e \n",qf1); +#endif + /* min(-1, 2, &s, &l, qf1, 1); */ + minny(0, 2, &s, &l, qf1, 1); + qa = l*(l-qd1)/(qd0*(qd0+qd1)); + qb = (l+qd0)*(qd1-l)/(qd0*qd1); + qc = l*(l+qd0)/(qd1*(qd0+qd1)); + } + else { + fx = qf1; qa = qb = 0.0; qc = 1.0; + } +#ifdef DEBUGPRAX + printf("after eventual min qd0=%14.8e qd1=%14.8e nl=%d\n",qd0, qd1,nl); +#endif + qd0 = qd1; + /* for (i=0; i x) x = y; +#ifdef DEBUGPRAX + printf(" I Y=%d %.7g",i,y); +#endif +#ifdef DEBUGPRAX + printf(" i=%d e(i) %.7g",i,e[i]); +#endif + } /* end i */ + /* + Accumulation of right hand transformations */ + /* for (i=n-1; i >= 0; i--) { */ /* FOR I := N STEP -1 UNTIL 1 DO */ + /* We should avoid the overflow in Golub */ + /* ab[n-1][n-1] = 1.0; */ + /* g = e[n-1]; */ + ab[n][n] = 1.0; + g = e[n]; + l = n; + + /* for (i=n; i >= 1; i--) { */ + for (i=n-1; i >= 1; i--) { /* n-1 loops, different from brent and golub*/ + if (g != 0.0) { + /* h = ab[i-1][i]*g; */ + h = ab[i][i+1]*g; + for (j=l; j<=n; j++) ab[j][i] = ab[i][j] / h; + for (j=l; j<=n; j++) { + /* h = ab[i][i+1]*g; */ + /* for (j=l; j= 0; k--) { */ + for (k=n; k>= 1; k--) { + kt = 0; +TestFsplitting: +#ifdef DEBUGPRAX + printf(" TestFsplitting: k=%d kt=%d\n",k,kt); + /* for(i=1;i<=n;i++)printf(" e(%d)=%.14f",i,e[i]);printf("\n"); */ +#endif + kt = kt+1; +/* TestFsplitting: */ + /* if (++kt > 30) { */ + if (kt > 30) { + e[k] = 0.0; + fprintf(stderr, "\n+++ MINFIT - Fatal error\n"); + fprintf ( stderr, " The QR algorithm failed to converge.\n" ); + } + /* for (l2=k; l2>=0; l2--) { */ + for (l2=k; l2>=1; l2--) { + l = l2; +#ifdef DEBUGPRAX + printf(" l e(l)< eps %d %.7g %.7g ",l,e[l], eps); +#endif + /* if (fabs(e[l]) <= eps) */ + if (fabs(e[l]) <= eps) + goto TestFconvergence; + /* if (fabs(q[l-1]) <= eps)*/ /* missing if ( 1 < l ){ *//* printf(" q(l-1)< eps %d %.7g %.7g ",l-1,q[l-2], eps); */ + if (fabs(q[l-1]) <= eps) + break; /* goto Cancellation; */ + } + Cancellation: +#ifdef DEBUGPRAX + printf(" Cancellation:\n"); +#endif + c = 0.0; s = 1.0; + for (i=l; i<=k; i++) { + f = s * e[i]; e[i] *= c; + /* f = s * e[i]; e[i] *= c; */ + if (fabs(f) <= eps) + goto TestFconvergence; + /* g = q[i]; */ + g = q[i]; + if (fabs(f) < fabs(g)) { + double fg = f/g; + h = fabs(g)*sqrt(1.0+fg*fg); + } + else { + double gf = g/f; + h = (f!=0.0 ? fabs(f)*sqrt(1.0+gf*gf) : 0.0); + } + /* COMMENT: THE ABOVE REPLACES Q(I):=H:=LONGSQRT(G*G+F*F) */ + /* WHICH MAY GIVE INCORRECT RESULTS IF THE */ + /* SQUARES UNDERFLOW OR IF F = G = 0; */ + + /* q[i] = h; */ + q[i] = h; + if (h == 0.0) { h = 1.0; g = 1.0; } + c = g/h; s = -f/h; + } +TestFconvergence: + #ifdef DEBUGPRAX + printf(" TestFconvergence: l=%d k=%d\n",l,k); +#endif + /* z = q[k]; */ + z = q[k]; + if (l == k) + goto Convergence; + /* shift from bottom 2x2 minor */ + /* x = q[l]; y = q[k-l]; g = e[k-1]; h = e[k]; */ /* Error */ + x = q[l]; y = q[k-1]; g = e[k-1]; h = e[k]; + f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y); + g = sqrt(f*f+1.0); + if (f <= 0.0) + f = ((x-z)*(x+z) + h*(y/(f-g)-h))/x; + else + f = ((x-z)*(x+z) + h*(y/(f+g)-h))/x; + /* next qr transformation */ + s = c = 1.0; + for (i=l+1; i<=k; i++) { +#ifdef DEBUGPRAXQR + printf(" Before Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]); +#endif + /* g = e[i]; y = q[i]; h = s*g; g *= c; */ + g = e[i]; y = q[i]; h = s*g; g *= c; + if (fabs(f) < fabs(h)) { + double fh = f/h; + z = fabs(h) * sqrt(1.0 + fh*fh); + } + else { + double hf = h/f; + z = (f!=0.0 ? fabs(f)*sqrt(1.0+hf*hf) : 0.0); + } + /* e[i-1] = z; */ + e[i-1] = z; +#ifdef DEBUGPRAXQR + printf(" Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]); +#endif + if (z == 0.0) + f = z = 1.0; + c = f/z; s = h/z; + f = x*c + g*s; g = - x*s + g*c; h = y*s; + y *= c; + /* for (j=0; j 1) { + printf("\n------------- enter function praxis -----------\n"); + printf("... current parameter settings ...\n"); + printf("... scaling ... %20.10e\n", scbd); + printf("... tol ... %20.10e\n", t); + printf("... maxstep ... %20.10e\n", h); + printf("... illc ... %20u\n", illc); + printf("... ktm ... %20u\n", ktm); + printf("... maxfun ... %20u\n", maxfun); + } + if (prin) print2(); + +mloop: + biter++; /* Added to count the loops */ + /* sf = d[0]; */ + /* s = d[0] = 0.0; */ + printf("\n Big iteration %d \n",biter); + fprintf(ficlog,"\n Big iteration %d \n",biter); + sf = d[1]; + s = d[1] = 0.0; + + /* minimize along first direction V(*,1) */ +#ifdef DEBUGPRAX + printf(" Minimize along the first direction V(*,1). illc=%d\n",illc); + /* fprintf(ficlog," Minimize along the first direction V(*,1).\n"); */ +#endif +#ifdef DEBUGPRAX2 + printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); +#endif + /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */ + minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */ +#ifdef DEBUGPRAX + printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); +#endif + if (s <= 0.0) + /* for (i=0; i < n; i++) */ + for (i=1; i <= n; i++) + v[i][1] = -v[i][1]; + /* if ((sf <= (0.9 * d[0])) || ((0.9 * sf) >= d[0])) */ + if ((sf <= (0.9 * d[1])) || ((0.9 * sf) >= d[1])) + /* for (i=1; i 0); +next: + kl = k; + df = 0.0; + if (illc) { /* random step to get off resolution valley */ +#ifdef DEBUGPRAX + printf(" A random step follows, to avoid resolution valleys.\n"); + matprint(" before rand, vectors:",v,n,n); +#endif + for (i=1; i<=n; i++) { +#ifdef NOBRENTRAND + r = drandom(); +#else + seed=i; + /* seed=i+1; */ +#ifdef DEBUGRAND + printf(" Random seed=%d, brent i=%d",seed,i); /* YYYY i=5 j=1 vji= -0.0001170073 */ +#endif + r = randbrent ( &seed ); +#endif +#ifdef DEBUGRAND + printf(" Random r=%.7g \n",r); +#endif + z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (r - 0.5); + /* z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (drandom() - 0.5); */ + + s = z[i]; + for (j=1; j <= n; j++) + x[j] += s * v[j][i]; + } +#ifdef DEBUGRAND + matprint(" after rand, vectors:",v,n,n); +#endif +#ifdef NR_SHIFT + fx = (*fun)((x-1), n); +#else + fx = (*fun)(x, n); +#endif + /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ + nf++; + } + /* minimize along non-conjugate directions */ +#ifdef DEBUGPRAX + printf(" Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n); + /* fprintf(ficlog," Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n); */ +#endif + /* for (k2=k; k2 df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc); +#endif + illc = 1; + goto next; + } +#ifdef DEBUGPRAX + printf("\n SUCCESS, BREAKS inner loop K(=%d) because DF is big, fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e <= df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc); +#endif + + /* if ((k == 1) && (prin > 1)){ /\* be careful k=2 *\/ */ + if ((k == 2) && (prin > 1)){ /* be careful k=2 */ +#ifdef DEBUGPRAX + printf(" NEW D The second difference array d:\n" ); + /* fprintf(ficlog, " NEW D The second difference array d:\n" ); */ +#endif + vecprint(" NEW D The second difference array d:",d,n); + } + /* minimize along conjugate directions */ + /* + Minimize along the "conjugate" directions V(*,1),...,V(*,K-1). + */ +#ifdef DEBUGPRAX + printf("Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); + /* fprintf(ficlog,"Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); */ +#endif + /* for (k2=0; k2<=k-1; k2++) { */ + for (k2=1; k2<=k-1; k2++) { + s = 0.0; + /* min(k2-1, 2, &d[k2-1], &s, fx, 0); */ + minny(k2, 2, &d[k2], &s, fx, 0); + } + f1 = fx; + fx = sf; + lds = 0.0; + /* for (i=0; i small_windows) { +#ifdef DEBUGPRAX + printf("lds big enough to throw direction V(*,kl=%d). If no random step was taken, V(*,KL) is the 'non-conjugate' direction along which the greatest improvement was made.\n",kl); + matprint(" before shift new conjugate vectors:",v,n,n); +#endif + for (i=kl-1; i>=k; i--) { + /* for (j=0; j < n; j++) */ + for (j=1; j <= n; j++) + /* v[j][i+1] = v[j][i]; */ /* This is v[j][i+1]=v[j][i] i=kl-1 to k */ + v[j][i+1] = v[j][i]; /* This is v[j][i+1]=v[j][i] i=kl-1 to k */ + /* v[j][i+1] = v[j][i]; */ + /* d[i+1] = d[i];*/ /* last is d[k+1]= d[k] */ + d[i+1] = d[i]; /* last is d[k]= d[k-1] */ + } +#ifdef DEBUGPRAX + matprint(" after shift new conjugate vectors:",v,n,n); +#endif /* d[k] = 0.0; */ + d[k] = 0.0; + for (i=1; i <= n; i++) + v[i][k] = y[i] / lds; + /* v[i][k] = y[i] / lds; */ +#ifdef DEBUGPRAX + printf("Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x). d2=%14.7g lds=%.10f\n",k,d[k],lds); + /* fprintf(ficlog,"Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x).\n",k); */ + matprint(" before min new conjugate vectors:",v,n,n); +#endif + /* min(k-1, 4, &d[k-1], &lds, f1, 1); */ + minny(k, 4, &d[k], &lds, f1, 1); +#ifdef DEBUGPRAX + printf(" after min d(k)=%d %.7g lds=%14f\n",k,d[k],lds); + matprint(" after min vectors:",v,n,n); +#endif + if (lds <= 0.0) { + lds = -lds; +#ifdef DEBUGPRAX + printf(" lds changed sign lds=%.14f k=%d\n",lds,k); +#endif + /* for (i=0; i 0){ +#ifdef DEBUGPRAX + printf(" k=%d",k); + /* fprintf(ficlog," k=%d",k); */ +#endif + print2();/* n, x, prin, fx, nf, nl ); */ + } + t2 = 0.0; + /* for (i=0; i (0.5 * t2)) + kt = 0; + else + kt++; +#ifdef DEBUGPRAX + printf("if kt=%d >? ktm=%d gotoL2 loop\n",kt,ktm); +#endif + if (kt > ktm){ + if ( 0 < prin ){ + /* printf("\nr8vec_print\n X:\n"); */ + /* fprintf(ficlog,"\nr8vec_print\n X:\n"); */ + vecprint ("END X:", x, n ); + } + goto fret; + } +#ifdef DEBUGPRAX + matprint(" end of L2 loop vectors:",v,n,n); +#endif + + } + /* printf("The inner loop ends here.\n"); */ + /* fprintf(ficlog,"The inner loop ends here.\n"); */ + /* + The inner loop ends here. + + Try quadratic extrapolation in case we are in a curved valley. + */ +#ifdef DEBUGPRAX + printf("Try QUAD ratic extrapolation in case we are in a curved valley.\n"); +#endif + /* try quadratic extrapolation in case */ + /* we are stuck in a curved valley */ + quad(); + dn = 0.0; + /* for (i=0; i 2) + matprint(" NEW DIRECTIONS vectors:",v,n,n); + /* for (j=0; j 1.0) { /* scale axis to reduce condition number */ +#ifdef DEBUGPRAX + printf("Scale the axes to try to reduce the condition number.\n"); +#endif + /* fprintf(ficlog,"Scale the axes to try to reduce the condition number.\n"); */ + s = vlarge; + /* for (i=0; i z[i]) + s = z[i]; + } + /* for (i=0; i scbd) { + sl = 1.0 / scbd; + z[i] = scbd; + } + } + } + for (i=1; i<=n; i++) + /* for (j=0; j<=i-1; j++) { */ + /* for (j=1; j<=i; j++) { */ + for (j=1; j<=i-1; j++) { + s = v[i][j]; + v[i][j] = v[j][i]; + v[j][i] = s; + } +#ifdef DEBUGPRAX + printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n"); +#endif + /* + MINFIT finds the singular value decomposition of V. + + This gives the principal values and principal directions of the + approximating quadratic form without squaring the condition number. + */ + #ifdef DEBUGPRAX + printf(" MINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n approximating quadratic form without squaring the condition number...\n"); +#endif + + minfit(n, macheps, vsmall, v, d); + /* for(i=0; i 1.0) { +#ifdef DEBUGPRAX + printf(" Unscale the axes.\n"); +#endif + /* for (i=0; i large) + d[i] = vsmall; + else if ((dn * d[i]) < small_windows) + d[i] = vlarge; + else + d[i] = 1.0 / dni / dni; /* added for compatibility with buckhardt but not brent */ + /* d[i] = pow(dn * d[i],-2.0); */ + } +#ifdef DEBUGPRAX + vecprint ("\n Before sort Eigenvalues of a:",d,n ); +#endif + + sort(); /* the new eigenvalues and eigenvectors */ +#ifdef DEBUGPRAX + vecprint( " After sort the eigenvalues ....\n", d, n); + matprint( " After sort the eigenvectors....\n", v, n,n); +#endif +#ifdef DEBUGPRAX + printf(" Determine the smallest eigenvalue.\n"); +#endif + /* dmin = d[n-1]; */ + dmin = d[n]; + if (dmin < small_windows) + dmin = small_windows; + /* + The ratio of the smallest to largest eigenvalue determines whether + the system is ill conditioned. + */ + + /* illc = (m2 * d[0]) > dmin; */ + illc = (m2 * d[1]) > dmin; +#ifdef DEBUGPRAX + printf(" The ratio of the smallest to largest eigenvalue determines whether\n the system is ill conditioned=%d . dmin=%.10lf < m2=%.10lf * d[1]=%.10lf \n",illc, dmin,m2, d[1]); +#endif + + if ((prin > 2) && (scbd > 1.0)) + vecprint("\n The scale factors:",z,n); + if (prin > 2) + vecprint(" Principal values (EIGEN VALUES OF A) of the quadratic form:",d,n); + if (prin > 2) + matprint(" The principal axes (EIGEN VECTORS OF A:",v,n, n); + + if ((maxfun > 0) && (nf > maxfun)) { + if (prin) + printf("\n... maximum number of function calls reached ...\n"); + goto fret; + } +#ifdef DEBUGPRAX + printf("Goto main loop\n"); +#endif + goto mloop; /* back to main loop */ + +fret: + if (prin > 0) { + vecprint("\n X:", x, n); + /* printf("\n... ChiSq reduced to %20.10e ...\n", fx); */ + /* printf("... after %20u function calls.\n", nf); */ + } + free_vector(d, 1, n); + free_vector(y, 1, n); + free_vector(z, 1, n); + free_vector(q0, 1, n); + free_vector(q1, 1, n); + free_matrix(v, 1, n, 1, n); + /* double *d, *y, *z, */ + /* *q0, *q1, **v; */ + free_vector(tflin, 1, n); + /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */ + free_vector(e, 1, n); + /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */ + + return(fx); +} + +/* end praxis gegen */ /*************** powell ************************/ /* @@ -2678,7 +4184,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret curr_time = *localtime(&rcurr_time); /* printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); */ /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */ - Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */ + /* Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /\* Big iteration, i.e on ncovmodel cycle *\/ */ + Bigter=(*iter - (*iter-1) % n)/n +1; /* Big iteration, i.e on ncovmodel cycle */ printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec); @@ -2749,9 +4256,10 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret fprintf(ficlog," - if your program needs %d BIG iterations (%d iterations) to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); } } - 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); + for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */ + for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales. xi is not changed but one dim xit */ + + fptt=(*fret); /* Computes likelihood for parameters xit */ #ifdef DEBUG printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); @@ -2759,33 +4267,33 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); #ifdef LINMINORIGINAL - linmin(p,xit,n,fret,func); /* New point i minimizing in direction i has coordinates p[j].*/ + linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit, i has coordinates p[j].*/ /* xit[j] gives the n coordinates of direction i as input.*/ /* *fret gives the maximum value on direction xit */ #else linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ - flatdir[i]=flat; /* Function is vanishing in that direction i */ + flatdir[i]=flat; /* Function is vanishing in that direction i */ #endif - /* Outputs are fret(new point p) p is updated and xit rescaled */ + /* Outputs are fret(new point p) p is updated and xit rescaled */ if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */ - /* because that direction will be replaced unless the gain del is small */ - /* in comparison with the 'probable' gain, mu^2, with the last average direction. */ - /* Unless the n directions are conjugate some gain in the determinant may be obtained */ - /* with the new direction. */ - del=fabs(fptt-(*fret)); - ibig=i; + /* because that direction will be replaced unless the gain del is small */ + /* in comparison with the 'probable' gain, mu^2, with the last average direction. */ + /* Unless the n directions are conjugate some gain in the determinant may be obtained */ + /* with the new direction. */ + del=fabs(fptt-(*fret)); + ibig=i; } #ifdef DEBUG printf("%d %.12e",i,(*fret)); fprintf(ficlog,"%d %.12e",i,(*fret)); for (j=1;j<=n;j++) { - xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); - printf(" x(%d)=%.12e",j,xit[j]); - fprintf(ficlog," x(%d)=%.12e",j,xit[j]); + xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); + printf(" x(%d)=%.12e",j,xit[j]); + fprintf(ficlog," x(%d)=%.12e",j,xit[j]); } for(j=1;j<=n;j++) { - printf(" p(%d)=%.12e",j,p[j]); - fprintf(ficlog," p(%d)=%.12e",j,p[j]); + printf(" p(%d)=%.12e",j,p[j]); + fprintf(ficlog," p(%d)=%.12e",j,p[j]); } printf("\n"); fprintf(ficlog,"\n"); @@ -2793,6 +4301,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret } /* end loop on each direction i */ /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */ /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ + /* New value of last point Pn is not computed, P(n-1) */ for(j=1;j<=n;j++) { if(flatdir[j] >0){ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); @@ -2845,13 +4354,19 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret return; } /* enough precision */ if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ + for (j=1;j<=n;j++) { /* Computes the extrapolated point and value f3, P_0 + 2 (P_n-P_0)=2Pn-P0 and xit is direction Pn-P0 */ ptt[j]=2.0*p[j]-pt[j]; - xit[j]=p[j]-pt[j]; - pt[j]=p[j]; - } + xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-P_0 */ +#ifdef DEBUG + printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]); +#endif + pt[j]=p[j]; /* New P0 is Pn */ + } +#ifdef DEBUG + printf("\n"); +#endif fptt=(*func)(ptt); /* f_3 */ -#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */ +#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in directions until some iterations are done */ if (*iter <=4) { #else #endif @@ -2870,10 +4385,10 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */ /* Even if f3 0 */ /* mu² and del² are equal when f3=f1 */ - /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */ - /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */ - /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */ - /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */ + /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */ + /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */ + /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */ + /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */ #ifdef NRCORIGINAL t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/ #else @@ -2895,7 +4410,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret if (t < 0.0) { /* Then we use it for new direction */ #else if (directest*t < 0.0) { /* Contradiction between both tests */ - printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); + printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); @@ -2930,17 +4445,10 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret } } #endif -#ifdef POWELLORIGINCONJUGATE 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 */ } -#else - for (j=1;j<=n-1;j++) { - xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */ - xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ - } -#endif #ifdef LINMINORIGINAL #else for (j=1, flatd=0;j<=n;j++) { @@ -2969,15 +4477,6 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret #endif printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */ - /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */ - /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */ - - /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */ - /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$. As $P^{(2)}_1$ and */ - /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */ - /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate. } */ - #ifdef DEBUG printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); @@ -2990,6 +4489,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret fprintf(ficlog,"\n"); #endif } /* end of t or directest negative */ + printf(" Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n); + fprintf(ficlog," Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n); #ifdef POWELLNOF3INFF1TEST #else } /* end if (fptt < fp) */ @@ -3198,7 +4699,9 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret first++; } - /* 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); */ + /* 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); */ free_vector(min,1,nlstate); free_vector(max,1,nlstate); free_vector(meandiff,1,nlstate); @@ -3233,7 +4736,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* 0.51326036147820708, 0.48673963852179264} */ /* If we start from prlim again, prlim tends to a constant matrix */ - int i, ii,j,k, k1; + int i, ii,j, k1; int first=0; double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ @@ -3500,9 +5003,9 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) /* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too. * Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. */ - int i, ii, j,k; + int ii, j; - double **out, **pmij(); + double **pmij(); double sumnew=0.; double agefin; double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ @@ -3715,11 +5218,11 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in */ - int i, j, d, h, k, k1; + int i, j, d, h, k1; double **out, cov[NCOVMAX+1]; double **newm; double agexact; - double agebegin, ageend; + /*double agebegin, ageend;*/ /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) @@ -3896,11 +5399,11 @@ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, do The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output */ - int i, j, d, h, k, k1; + int i, j, d, h, k1; double **out, cov[NCOVMAX+1], **bmij(); double **newm, ***newmm; double agexact; - double agebegin, ageend; + /*double agebegin, ageend;*/ double **oldm, **savm; newmm=po; /* To be saved */ @@ -4419,7 +5922,7 @@ double func( double *x) double funcone( double *x) { /* Same as func but slower because of a lot of printf and if */ - int i, ii, j, k, mi, d, kk, kv=0, kf=0; + int i, ii, j, k, mi, d, kv=0, kf=0; int ioffset=0; int ipos=0,iposold=0,ncovv=0; @@ -4969,13 +6472,13 @@ void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, lon void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double [])) { - int i,j,k, jk, jkk=0, iter=0; + int i,j, jkk=0, iter=0; double **xi; - double fret; - double fretone; /* Only one call to likelihood */ + /*double fret;*/ + /*double fretone;*/ /* Only one call to likelihood */ /* char filerespow[FILENAMELENGTH];*/ - double * p1; /* Shifted parameters from 0 instead of 1 */ + /*double * p1;*/ /* Shifted parameters from 0 instead of 1 */ #ifdef NLOPT int creturn; nlopt_opt opt; @@ -4991,7 +6494,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub for (i=1;i<=npar;i++) /* Starting with canonical directions j=1,n xi[i=1,n][j] */ for (j=1;j<=npar;j++) xi[i][j]=(i==j ? 1.0 : 0.0); - printf("Powell\n"); fprintf(ficlog,"Powell\n"); + printf("Powell-prax\n"); fprintf(ficlog,"Powell-prax\n"); strcpy(filerespow,"POW_"); strcat(filerespow,fileres); if((ficrespow=fopen(filerespow,"w"))==NULL) { @@ -5055,7 +6558,23 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub } powell(p,xi,npar,ftol,&iter,&fret,flatdir,func); #else /* FLATSUP */ - powell(p,xi,npar,ftol,&iter,&fret,func); +/* powell(p,xi,npar,ftol,&iter,&fret,func);*/ +/* praxis ( t0, h0, n, prin, x, beale_f ); */ + int prin=1; + double h0=0.25; + double macheps; + double fmin; + macheps=pow(16.0,-13.0); +/* #include "praxis.h" */ + /* Be careful that praxis start at x[0] and powell start at p[1] */ + /* praxis ( ftol, h0, npar, prin, p, func ); */ +/* p1= (p+1); */ /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */ +printf("Praxis Gegenfurtner \n"); +fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog); +/* praxis ( ftol, h0, npar, prin, p1, func ); */ + /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */ + fmin = praxis(ftol,macheps, h0, npar, prin, p, func); +printf("End Praxis\n"); #endif /* FLATSUP */ #ifdef LINMINORIGINAL @@ -6186,7 +7705,7 @@ void prevalence(double ***probs, double agemin, double agemax, int **s, double * int i, m, jk, j1, bool, z1,j, iv; int mi; /* Effective wave */ int iage; - double agebegin, ageend; + double agebegin; /*, ageend;*/ double **prop; double posprop; @@ -6425,10 +7944,10 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc if(j==0) j=1; /* Survives at least one month after exam */ else if(j<0){ nberr++; - printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); j=1; /* Temporary Dangerous patch */ printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); - fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); } k=k+1; @@ -6462,8 +7981,8 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ if(j<0){ nberr++; - printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + printf("Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); } sum=sum+j; } @@ -8073,7 +9592,7 @@ void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \ double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){ - int jj1, k1, i1, cpt, k4, nres; + int jj1, k1, cpt, nres; /* In fact some results are already printed in fichtm which is open */ fprintf(fichtm,"
  • Result files (first order: no variance)\n \
  • Result files (second order (variance)\n \ @@ -8210,21 +9729,21 @@ divided by h: hPij/h : %s_%d-3- ",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
    \n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages.     %s_%d-%d-%d.svg
    ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); + fprintf(fichtm,"
    \n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. %s_%d-%d-%d.svg
    ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
    ",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
    \n- Survival functions in state %d and in any other live state (total).\ - And probability to be observed in various states (up to %d) being in state %d at different ages. \ + And probability to be observed in various states (up to %d) being in state %d at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \ %s_%d-%d-%d.svg
    ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
    ",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } /* Period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
    \n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
    ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); + fprintf(fichtm,"
    \n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be alive in state %d some years after. %s_%d-%d-%d.svg
    ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
    ",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } @@ -8249,8 +9768,8 @@ divided by h: hPij/h : %s_%d-3- /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
    \n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ - from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \ - account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \ + from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \ + account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \ with weights corresponding to observed prevalence at different ages.     %s_%d-%d-%d.svg", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
    ",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_")); fprintf(fichtm," ", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); @@ -9937,10 +11456,10 @@ void prevforecast(char fileres[], double dateintmean, double dateprojd, double d */ /* double anprojd, mprojd, jprojd; */ /* double anprojf, mprojf, jprojf; */ - int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; + int yearp, stepsize, hstepm, nhstepm, j, k, i, h, nres=0; double agec; /* generic age */ - double agelim, ppij, yp,yp1,yp2; - double *popeffectif,*popcount; + double agelim, ppij; + /*double *popcount;*/ double ***p3mat; /* double ***mobaverage; */ char fileresf[FILENAMELENGTH]; @@ -10088,10 +11607,10 @@ void prevforecast(char fileres[], double dateintmean, double dateprojd, double d anback2 year of end of backprojection (same day and month as back1). prevacurrent and prev are prevalences. */ - int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0; + int yearp, stepsize, hstepm, nhstepm, j, k, i, h, nres=0; double agec; /* generic age */ - double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/ - double *popeffectif,*popcount; + double agelim, ppij, ppi; /* ,jintmean,mintmean,aintmean;*/ + /*double *popcount;*/ double ***p3mat; /* double ***mobaverage; */ char fileresfb[FILENAMELENGTH]; @@ -10774,7 +12293,7 @@ void printinggnuplotmort(char fileresu[], char optionfilefiname[], double agemin char dirfileres[132],optfileres[132]; - int ng; + /*int ng;*/ /*#ifdef windows */ @@ -10798,7 +12317,7 @@ int readdata(char datafile[], int firstobs, int lastobs, int *imax) /*-------- data file ----------*/ FILE *fic; char dummy[]=" "; - int i=0, j=0, n=0, iv=0, v; + int i = 0, j = 0, n = 0, iv = 0;/* , v;*/ int lstra; int linei, month, year,iout; int noffset=0; /* This is the offset if BOM data file */ @@ -11415,7 +12934,7 @@ int decodemodel( char model[], int lastobs) */ /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ { - int i, j, k, ks, v; + int i, j, k, ks;/* , v;*/ int n,m; int j1, k1, k11, k12, k2, k3, k4; char modelsav[300]; @@ -12843,7 +14362,7 @@ int hPijx(double *p, int bage, int fage){ int agelim; int hstepm; int nhstepm; - int h, i, i1, j, k, k4, nres=0; + int h, i, i1, j, k, nres=0; double agedeb; double ***p3mat; @@ -13049,7 +14568,7 @@ int main(int argc, char *argv[]) double fret; double dum=0.; /* Dummy variable */ - double ***p3mat; + /* double*** p3mat;*/ /* double ***mobaverage; */ double wald; @@ -13062,7 +14581,7 @@ int main(int argc, char *argv[]) char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */ - int c, h , cpt, c2; + int c, h; /* c2; */ int jl=0; int i1, j1, jk, stepsize=0; int count=0; @@ -13097,7 +14616,7 @@ int main(int argc, char *argv[]) double ***delti3; /* Scale */ double *delti; /* Scale */ double ***eij, ***vareij; - double **varpl; /* Variances of prevalence limits by age */ + //double **varpl; /* Variances of prevalence limits by age */ double *epj, vepp; @@ -13155,7 +14674,7 @@ int main(int argc, char *argv[]) getcwd(pathcd, size); #endif syscompilerinfo(0); - printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); + printf("\nIMaCh prax version %s, %s\n%s",version, copyright, fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); if(!fgets(pathr,FILENAMELENGTH,stdin)){ @@ -14057,7 +15576,7 @@ This file: %s
    Title=%s
    Datafile=%s /* Calculates basic frequencies. Computes observed prevalence at single age and for any valid combination of covariates and prints on file fileres'p'. */ - freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ + freqsummary(fileres, p, pstart, (double)agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); @@ -14148,7 +15667,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
    \n",\ #ifdef GSL printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); #else - printf("Powell\n"); fprintf(ficlog,"Powell\n"); + printf("Powell-mort\n"); fprintf(ficlog,"Powell-mort\n"); #endif strcpy(filerespow,"POW-MORT_"); strcat(filerespow,fileresu); @@ -14251,7 +15770,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
    \n",\ for(i=1; i <=NDIM; i++) for(j=i+1;j<=NDIM;j++) - matcov[i][j]=matcov[j][i]; + matcov[i][j]=matcov[j][i]; printf("\nCovariance matrix\n "); fprintf(ficlog,"\nCovariance matrix\n "); @@ -14705,7 +16224,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa } /* Results */ - /* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */ + /* Value of covariate in each resultine will be computed (if product) and sorted according to model rank */ /* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */ precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1); endishere=0; @@ -15108,9 +16627,9 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa /* */ if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */ continue; - printf("\n# model %s \n#****** Result for:", model); /* HERE model is empty */ - fprintf(ficrest,"\n# model %s \n#****** Result for:", model); - fprintf(ficlog,"\n# model %s \n#****** Result for:", model); + printf("\n# model=1+age+%s \n#****** Result for:", model); /* HERE model is empty */ + fprintf(ficrest,"\n# model=1+age+%s \n#****** Result for:", model); + fprintf(ficlog,"\n# model=1+age+%s \n#****** Result for:", model); /* It might not be a good idea to mix dummies and quantitative */ /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */ for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */ @@ -15305,6 +16824,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); } /* mle==-3 arrives here for freeing */ /* endfree:*/ + if(mle!=-3) free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/ free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -15366,7 +16886,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa free_ivector(TmodelInvind,1,NCOVMAX); free_ivector(TmodelInvQind,1,NCOVMAX); - free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/ + /* free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /\* Could be elsewhere ?*\/ */ free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); /* free_imatrix(codtab,1,100,1,10); */