From a0442ed58dee48a521ea053083ea967894507898 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Wed, 11 Jul 2018 22:19:56 +0200 Subject: New upstream version 2.0.1+repack --- numlib/powell.c | 712 ++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 563 insertions(+), 149 deletions(-) (limited to 'numlib/powell.c') diff --git a/numlib/powell.c b/numlib/powell.c index 47acc15..7fb57a7 100755 --- a/numlib/powell.c +++ b/numlib/powell.c @@ -40,8 +40,10 @@ #include "numsup.h" #include "powell.h" -#undef SLOPE_SANITY_CHECK /* exermental */ -#undef ABSTOL /* Make tollerance absolute */ +#undef SLOPE_SANITY_CHECK /* [und] expermental */ +#undef ABSTOL /* [und] Make tollerance absolute */ +#undef USE_LINMIND /* [und] Use limind for conjgrad (typically slower) */ + /* Some debugging printfs (not comprehensive): */ #undef PDEBUG /* Powell debug */ #undef CDEBUG /* Conjgrad debug */ @@ -71,6 +73,7 @@ # define LDBG(xxx) #endif +/* -------------------------------------- */ /* Standard interface for powell function */ /* return 0 on sucess, 1 on failure due to excessive itterions */ /* Result will be in cp */ @@ -91,10 +94,10 @@ void (*prog)(void *pdata, int perc), /* Optional progress percentage callback * void *pdata /* Opaque data needed by prog() */ ) { int i; - double **dmtx; /* Direction vector */ - double *spt; /* Sarting point before exploring all the directions */ - double *xpt; /* Extrapolated point */ - double *svec; /* Search vector */ + double **dmtx, *_dmtx[10], __dmtx[10 * 10] = { 0.0 }; /* Direction vector */ + double *spt, _spt[10]; /* Sarting point before exploring all the directions */ + double *xpt, _xpt[10]; /* Extrapolated point */ + double *svec, _svec[10]; /* Search vector */ int iter; double retv; /* Returned function value at p */ double stopth; /* Current stop threshold */ @@ -102,10 +105,20 @@ void *pdata /* Opaque data needed by prog() */ double curdel; /* Current change in function value */ int pc = 0; /* Percentage complete */ - dmtx = dmatrixz(0, di-1, 0, di-1); /* Zero filled */ - spt = dvector(0,di-1); - xpt = dvector(0,di-1); - svec = dvector(0,di-1); + if (di <= 10) { + int j; + for (j = i = 0; i < di; i++, j += di) + _dmtx[i] = __dmtx + j; + dmtx = _dmtx; + spt = _spt; + xpt = _xpt; + svec = _svec; + } else { + dmtx = dmatrixz(0, di-1, 0, di-1); /* Zero filled */ + spt = dvector(0, di-1); + xpt = dvector(0, di-1); + svec = dvector(0, di-1); + } /* Create initial direction matrix by */ /* placing search start on diagonal */ @@ -213,11 +226,14 @@ void *pdata /* Opaque data needed by prog() */ } //printf("~1 iters = %d\n",iter); + /* Free up all the temporary vectors and matrix */ - free_dvector(svec,0,di-1); - free_dvector(xpt,0,di-1); - free_dvector(spt,0,di-1); - free_dmatrix(dmtx, 0, di-1, 0, di-1); + if (di > 10) { + free_dvector(svec, 0, di-1); + free_dvector(xpt, 0, di-1); + free_dvector(spt, 0, di-1); + free_dmatrix(dmtx, 0, di-1, 0, di-1); + } if (prog != NULL) /* Report final progress */ prog(pdata, 100); @@ -232,12 +248,323 @@ void *pdata /* Opaque data needed by prog() */ return 1; /* Failed due to execessive itterations */ } +/* - - - - - - - - - - - - - - - - - */ + +#define POWELL_GOLD 1.618034 +#define POWELL_CGOLD 0.3819660 +#define POWELL_MAXIT 100 + +/* Line bracketing and minimisation routine. */ +/* Return value at minimum. */ +double linmin( +double cp[], /* Start point, and returned value */ +double xi[], /* Search vector */ +int di, /* Dimensionality */ +#ifdef ABSTOL +double ftol, /* Absolute tolerance to stop on */ +#else +double ftol, /* Relative tolerance to stop on */ +#endif +double (*func)(void *fdata, double tp[]), /* Error function to evaluate */ +void *fdata) /* Opaque data for func() */ +{ + int i; + double ax, xx, bx; /* Search vector multipliers */ + double af, xf, bf; /* Function values at those points */ + double *xt, _xt[10]; /* Trial point */ + + if (di <= 10) + xt = _xt; + else + xt = dvector(0, di-1); /* Vector for trial point */ + + /* -------------------------- */ + /* First bracket the solution */ + + LDBG((" linmin: Bracketing solution\n")) + + /* The line is measured as startpoint + offset * search vector. */ + /* (Search isn't symetric, but it seems to depend on cp being */ + /* best current solution ?) */ + ax = 0.0; + for (i = 0; i < di; i++) + xt[i] = cp[i] + ax * xi[i]; + af = (*func)(fdata, xt); + + /* xx being vector offset 0.618 */ + xx = 1.0/POWELL_GOLD; + for (i = 0; i < di; i++) + xt[i] = cp[i] + xx * xi[i]; + xf = (*func)(fdata, xt); + + LDBG((" linmin: Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) + + /* Fix it so that we are decreasing from point a -> x */ + if (xf > af) { + double tt; + tt = ax; ax = xx; xx = tt; + tt = af; af = xf; xf = tt; + } + LDBG((" linmin: Ordered Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) + + bx = xx + POWELL_GOLD * (xx-ax); /* Guess b beyond a -> x */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + bx * xi[i]; + bf = (*func)(fdata, xt); + + LDBG((" linmin: Initial bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf)) + +#ifdef SLOPE_SANITY_CHECK + /* If we're not seeing a slope indicitive of progress */ + /* of order ftol, give up straight away */ + if (2000.0 * fabs(xf - bf) <= ftol * (fabs(xf) + fabs(bf)) + && 2000.0 * fabs(af - xf) <= ftol * (fabs(af) + fabs(xf))) { + LDBG((" linmin: giving up because slope is too shallow\n")) + if (xt != _xt) + free_dvector(xt,0,di-1); + + if (bf < xf) { + xf = bf; + xx = bx; + } + goto done; + } +#endif /* SLOPE_SANITY_CHECK */ + + /* While not bracketed */ + while (xf > bf) { + double ulim, ux, uf; + double tt, r, q; + + LDBG((" linmin: Not bracketed because xf %f > bf %f\n",xf, bf)) + LDBG((" ax = %f, xx = %f, bx = %f\n",ax,xx,bx)) + + /* Compute ux by parabolic interpolation from a, x & b */ + q = (xx - bx) * (xf - af); + r = (xx - ax) * (xf - bf); + tt = q - r; + if (tt >= 0.0 && tt < 1e-20) /* If +ve too small */ + tt = 1e-20; + else if (tt <= 0.0 && tt > -1e-20) /* If -ve too small */ + tt = -1e-20; + ux = xx - ((xx - bx) * q - (xx - ax) * r) / (2.0 * tt); + ulim = xx + 100.0 * (bx - xx); /* Extrapolation limit */ + +//printf("~1 ux = %f, ulim = %f\n",ux,ulim); + if ((xx - ux) * (ux - bx) > 0.0) { /* u is between x and b */ + + for (i = 0; i < di; i++) /* Evaluate u */ + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); + +//printf("~1 u is between x and b, uf = %f\n",uf); + + if (uf < bf) { /* Minimum is between x and b */ +//printf("~1 min is between x and b\n"); + ax = xx; af = xf; + xx = ux; xf = uf; + break; + } else if (uf > xf) { /* Minimum is between a and u */ +//printf("~1 min is between a and u\n"); + bx = ux; bf = uf; + break; + } + + /* Parabolic fit didn't work, look further out in direction of b */ + ux = bx + POWELL_GOLD * (bx-xx); +//printf("~1 parabolic fit didn't work,look further in direction of b (%f)\n",ux); + + } else if ((bx - ux) * (ux - ulim) > 0.0) { /* u is between b and limit */ + for (i = 0; i < di; i++) /* Evaluate u */ + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); + +//printf("~1 u is between b and limit uf = %f\n",uf); + if (uf > bf) { /* Minimum is between x and u */ +//printf("~1 min is between x and uf\n"); + ax = xx; af = xf; + xx = bx; xf = bf; + bx = ux; bf = uf; + break; + } + xx = bx; xf = bf; /* Continue looking */ + bx = ux; bf = uf; + ux = bx + POWELL_GOLD * (bx - xx); /* Test beyond b */ +//printf("~1 continue looking beyond b (%f)\n",ux); + + } else if ((ux - ulim) * (ulim - bx) >= 0.0) { /* u is beyond limit */ + ux = ulim; +//printf("~1 use limit\n"); + } else { /* u is to left side of x ? */ + ux = bx + POWELL_GOLD * (bx-xx); +//printf("~1 look gold beyond b (%f)\n",ux); + } + /* Evaluate u, and move into place at b */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); +//printf("~1 lookup ux %f value uf = %f\n",ux,uf); + ax = xx; af = xf; + xx = bx; xf = bf; + bx = ux; bf = uf; +//printf("~1 move along to the right (a<-x, x<-b, b- x -> b */ +//printf("~1 got bracketed minimum at %f (%f), %f (%f), %f (%f)\n",ax,af,xx,xf,bx,bf); + + /* --------------------------------------- */ + /* Now use brent minimiser bewteen a and b */ + { + /* a and b bracket solution */ + /* x is best function value so far */ + /* w is second best function value so far */ + /* v is previous second best, or third best */ + /* u is most recently tested point */ + double wx, vx, ux; /* Search vector multipliers */ + double wf, vf = 0.0, uf; /* Function values at those points */ + int iter; + double de = 0.0; /* Distance moved on previous step */ + double e = 0.0; /* Distance moved on 2nd previous step */ + + /* Make sure a and b are in ascending order */ + if (ax > bx) { + double tt; + tt = ax; ax = bx; bx = tt; + tt = af; af = bf; bf = tt; + } + + wx = vx = xx; /* Initial values of other center points */ + wf = xf = xf; + + for (iter = 1; iter <= POWELL_MAXIT; iter++) { + double mx = 0.5 * (ax + bx); /* m is center of bracket values */ +#ifdef ABSTOL + double tol1 = ftol; /* Absolute tollerance */ +#else + double tol1 = ftol * fabs(xx) + 1e-10; +#endif + double tol2 = 2.0 * tol1; + + LDBG((" linmin it %d: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",iter,ax,af,xx,xf,bx,bf)) + + /* See if we're done */ +//printf("~1 linmin check %f <= %f\n",fabs(xx - mx), tol2 - 0.5 * (bx - ax)); + if (fabs(xx - mx) <= (tol2 - 0.5 * (bx - ax))) { + LDBG((" linmin: We're done because %e <= %e\n",fabs(xx - mx), tol2 - 0.5 * (bx - ax))) + break; + } + + LDBG((" linmin: e %e tol2 %e\n",e,tol1)) + + if (fabs(e) > tol1) { /* Do a trial parabolic fit */ + double te, p, q, r; + r = (xx-wx) * (xf-vf); + q = (xx-vx) * (xf-wf); + p = (xx-vx) * q - (xx-wx) * r; + q = 2.0 * (q - r); + if (q > 0.0) + p = -p; + else + q = -q; + te = e; /* Save previous e value */ + e = de; /* Previous steps distance moved */ + + LDBG((" linmin: Trial parabolic fit\n" )) + + if (fabs(p) >= fabs(0.5 * q * te) || p <= q * (ax-xx) || p >= q * (bx-xx)) { + /* Give up on the parabolic fit, and use the golden section search */ + e = ((xx >= mx) ? ax-xx : bx-xx); /* Override previous distance moved */ + de = POWELL_CGOLD * e; + LDBG((" linmin: Moving to golden section search\n" )) + } else { /* Use parabolic fit */ + de = p/q; /* Change in xb */ + ux = xx + de; /* Trial point according to parabolic fit */ + if ((ux - ax) < tol2 || (bx - ux) < tol2) { + if ((mx - xx) > 0.0) /* Don't use parabolic, use tol1 */ + de = tol1; /* tol1 is +ve */ + else + de = -tol1; + } + LDBG((" linmin: Using parabolic fit\n" )) + } + } else { /* Keep using the golden section search */ + e = ((xx >= mx) ? ax-xx : bx-xx); /* Override previous distance moved */ + de = POWELL_CGOLD * e; + LDBG((" linmin: Continuing golden section search\n" )) + } + + if (fabs(de) >= tol1) { /* If de moves as much as tol1 would */ + ux = xx + de; /* use it */ + LDBG((" linmin: ux = %f = xx %f + de %f\n",ux,xx,de)) + } else { /* else move by tol1 in direction de */ + if (de > 0.0) { + ux = xx + tol1; + LDBG((" linmin: ux = %f = xx %f + tol1 %e\n",ux,xx,tol1)) + } else { + ux = xx - tol1; + LDBG((" linmin: ux = %f = xx %f - tol1 %f\n",ux,xx,tol1)) + } + } + + /* Evaluate function */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); + + if (uf <= xf) { /* Found new best solution */ + LDBG((" linmin: found new best solution at %f val %f\n",ux,uf)) + if (ux >= xx) { + ax = xx; af = xf; /* New lower bracket */ + } else { + bx = xx; bf = xf; /* New upper bracket */ + } + vx = wx; vf = wf; /* New previous 2nd best solution */ + wx = xx; wf = xf; /* New 2nd best solution from previous best */ + xx = ux; xf = uf; /* New best solution from latest */ + } else { /* Found a worse solution */ + LDBG((" linmin: found new worse solution at %f val %f\n",ux,uf)) + LDBG((" linmin: current best at %f val %f\n",xx,xf)) + if (ux < xx) { + ax = ux; af = uf; /* New lower bracket */ + } else { + bx = ux; bf = uf; /* New upper bracket */ + } + if (uf <= wf || wx == xx) { /* New 2nd best solution, or equal best */ + vx = wx; vf = wf; /* New previous 2nd best solution */ + wx = ux; wf = uf; /* New 2nd best from latest */ + } else if (uf <= vf || vx == xx || vx == wx) { /* New 3rd best, or equal 1st & 2nd */ + vx = ux; vf = uf; /* New previous 2nd best from latest */ + } + } + } + /* !!! should do something if iter > POWELL_MAXIT !!!! */ + /* Solution is at xx, xf */ + } + + done:; + + /* Compute solution vector at xx */ + LDBG((" linmin: computing soln from best at %f val %f\n",xx,xf)) + for (i = 0; i < di; i++) + cp[i] += xx * xi[i]; + + if (xt != _xt) + free_dvector(xt,0,di-1); +//printf("~~~ line minimizer returning %e\n",xf); + return xf; +} + +#undef POWELL_GOLD +#undef POWELL_CGOLD +#undef POWELL_MAXIT + /* -------------------------------------- */ -/* Conjugate Gradient optimiser */ +/* Conjugate Gradient optimiser using partial derivatives. */ /* return 0 on sucess, 1 on failure due to excessive itterions */ /* Result will be in cp */ /* Note that we could use gradient in line minimiser, */ -/* but haven't bothered yet. */ +/* but this seems to be slower, so we don't use it. */ int conjgrad( double *rv, /* If not NULL, return the residual error */ int di, /* Dimentionality */ @@ -250,79 +577,93 @@ double ftol, /* Relative tollerance of error change to stop on */ #endif int maxit, /* Maximum iterations allowed */ double (*func)(void *fdata, double tp[]), /* Error function to evaluate */ -double (*dfunc)(void *fdata, double dp[], double tp[]), /* Gradient function to evaluate */ +double (*dfunc)(void *fdata, double dp[], double tp[]), /* Gradient & function to evaluate */ + /* dfunc() should return DFUNC_NRV if it doesn't return function value */ void *fdata, /* Opaque data needed by function */ void (*prog)(void *pdata, int perc), /* Optional progress percentage callback */ void *pdata /* Opaque data needed by prog() */ ) { int i, iter; - double *svec; /* Search vector */ - double *gvec; /* G direction vector */ - double *hvec; /* H direction vector */ + double *svec, _svec[10]; /* Search vector */ + double *ssvec, _ssvec[10]; /* s[] scaled search vector */ + double *gvec, _gvec[10]; /* G direction vector */ + double *hvec, _hvec[10]; /* H direction vector */ double retv; /* Returned function value at p */ double stopth; /* Current stop threshold */ double startdel = -1.0; /* Initial change in function value */ double curdel; /* Current change in function value */ - double svec_sca; /* initial svec scale factor */ + double brat; /* svec to s[] ratio */ + double svec_sca; /* svec scale factor */ int pc = 0; /* Percentage complete */ - svec = dvector(0,di-1); - gvec = dvector(0,di-1); - hvec = dvector(0,di-1); + if (di <= 10) { + svec = _svec; + ssvec = _ssvec; + gvec = _gvec; + hvec = _hvec; + } else { + svec = dvector(0, di-1); + ssvec = dvector(0, di-1); + gvec = dvector(0, di-1); + hvec = dvector(0, di-1); + } + + CDBG(("conjgrad with di %d\n", di)) if (prog != NULL) /* Report initial progress */ prog(pdata, pc); - /* Initial gradient function evaluation */ + /* Initial function and gradient evaluation */ retv = (*dfunc)(fdata, svec, cp); - - /* svec[] seems to be large after this. */ - /* Rescale it to conform to maximum of s[] */ - for (svec_sca = 0.0, i = 0; i < di; i++) { - if (fabs(svec[i]) > svec_sca) - svec_sca = fabs(svec[i]); + if (retv == DFUNC_NRV) + retv = (*func)(fdata, cp); + + /* svec[] seems to be large after this. Compute scaled version that */ + /* has maximum of s[] so that line search is guided by the search radius. */ + for (brat = 0.0, i = 0; i < di; i++) { + double rat = fabs(svec[i]) / fabs(s[i]); + if (rat > brat) + brat = rat; } - /* set scale so largest <= 1 */ - if (svec_sca < 1e-12) - svec_sca = 1.0; - else - svec_sca = 1.0/svec_sca; - - CDBG((" initial dir = %s\n", debPdv(di,svec))); - CDBG((" initial retv = %f\n",retv)); + svec_sca = 1.0/brat; /* Initial vector setup */ for (i = 0; i < di; i++) { - gvec[i] = hvec[i] = -svec[i]; /* Inverse gradient */ - svec[i] = s[i] * -svec[i] * svec_sca; /* Scale the search vector */ + svec[i] = gvec[i] = hvec[i] = -svec[i]; /* Inverse gradient */ + ssvec[i] = svec[i] * svec_sca; /* Scale the search vector to s[] size */ } - CDBG(("Initial svec = %s\n", debPdv(di,svec))); + + CDBG((" initial dir = %s\n", debPdv(di, ssvec))); + CDBG((" initial retv = %f\n",retv)); /* Itterate untill we converge on a solution, or give up. */ for (iter = 1; iter < maxit; iter++) { double gamden, gamnum, gam; double pretv; /* Previous function return value */ - CDBG(("conjrad: about to do linmin\n")) + CDBG(("conjrad it %d: about to do linmind\n",iter)) pretv = retv; - retv = linmin(cp, svec, di, ftol, func, fdata); +#ifdef USE_LINMIND + retv = linmind(cp, ssvec, di, ftol, func, dfunc, fdata); +#else + retv = linmin(cp, ssvec, di, ftol, func, fdata); +#endif #ifdef ABSTOL stopth = ftol; /* Absolute tollerance */ #else - stopth = ftol * 0.5 * (fabs(pretv) + fabs(retv) + DBL_EPSILON); // Old code + stopth = ftol * 0.5 * (fabs(pretv) + fabs(retv) + DBL_EPSILON); #endif curdel = fabs(pretv - retv); CDBG((" this retv = %f, pretv = %f, curdel = %f\n",retv,pretv,curdel)); if (startdel < 0.0) { startdel = curdel; - } else { + } else if (prog != NULL) { /* Update percentage */ int tt; tt = (int)(100.0 * pow((log(curdel) - log(startdel))/(log(stopth) - log(startdel)), 4.0) + 0.5); if (tt > pc && tt < 100) { pc = tt; - if (prog != NULL) /* Report initial progress */ - prog(pdata, pc); + prog(pdata, pc); /* Report initial progress */ } } @@ -336,7 +677,7 @@ void *pdata /* Opaque data needed by prog() */ CDBG(("conjrad: recomputing direction\n")) (*dfunc)(fdata, svec, cp); /* (Don't use retv as it wrecks stop test) */ - CDBG((" pderiv = %s\n", debPdv(di,svec))); + CDBG((" pderiv = %s\n", debPdv(di, svec))); /* Compute gamma */ for (gamnum = gamden = 0.0, i = 0; i < di; i++) { @@ -360,25 +701,26 @@ void *pdata /* Opaque data needed by prog() */ svec[i] = hvec[i] = gvec[i] + gam * hvec[i]; } - /* svec[] seems to be large after this. */ - /* Rescale it to conform to maximum of s[] */ - for (svec_sca = 0.0, i = 0; i < di; i++) { - if (fabs(svec[i]) > svec_sca) - svec_sca = fabs(svec[i]); + /* svec[] seems to be large after this. Compute scaled version that */ + /* has maximum of s[] so that line search is guided by the search radius. */ + for (brat = 0.0, i = 0; i < di; i++) { + double rat = fabs(svec[i]) / fabs(s[i]); + if (rat > brat) + brat = rat; } - /* set scale so largest <= 1 */ - if (svec_sca < 1e-12) - svec_sca = 1.0; - else - svec_sca = 1.0/svec_sca; + svec_sca = 1.0/brat; for (i = 0; i < di; i++) - svec[i] = svec[i] * s[i] * svec_sca; - CDBG((" svec = %s\n", debPdv(di,svec))); + ssvec[i] = svec[i] * svec_sca; + + CDBG((" ssvec = %s\n", debPdv(di,ssvec))); } /* Free up all the temporary vectors and matrix */ - free_dvector(hvec,0,di-1); - free_dvector(gvec,0,di-1); - free_dvector(svec,0,di-1); + if (di > 10) { + free_dvector(hvec, 0, di-1); + free_dvector(gvec, 0, di-1); + free_dvector(ssvec, 0, di-1); + free_dvector(svec, 0, di-1); + } if (prog != NULL) /* Report final progress */ prog(pdata, 100); @@ -394,14 +736,15 @@ void *pdata /* Opaque data needed by prog() */ return 1; /* Failed due to execessive itterations */ } -/*------------------------------*/ #define POWELL_GOLD 1.618034 -#define POWELL_CGOLD 0.3819660 #define POWELL_MAXIT 100 -/* Line bracketing and minimisation routine. */ +/* Line bracketing and minimisation routine using derivatives */ +/* This is not used, because it typically makes it slower */ +/* - it may take less itterations, but each itteration uses */ +/* a func() and dfunc() call, at least doubling itter overhead. */ /* Return value at minimum. */ -double linmin( +double linmind( double cp[], /* Start point, and returned value */ double xi[], /* Search vector */ int di, /* Dimensionality */ @@ -410,23 +753,29 @@ double ftol, /* Absolute tolerance to stop on */ #else double ftol, /* Relative tolerance to stop on */ #endif -double (*func)(void *fdata, double tp[]), /* Error function to evaluate */ +double (*func)(void *fdata, double tp[]), /* Error function to evaluate */ +double (*dfunc)(void *fdata, double dp[], double tp[]), /* Gradient function to evaluate */ + /* dfunc() should return DFUNC_NRV if it doesn't return function value */ void *fdata) /* Opaque data for func() */ { int i; double ax, xx, bx; /* Search vector multipliers */ double af, xf, bf; /* Function values at those points */ - double *xt, XT[10]; /* Trial point */ + double *xt, _xt[10]; /* Trial point */ + double *df, _df[10]; /* Derivative vector */ - if (di <= 10) - xt = XT; - else + if (di <= 10) { + xt = _xt; + df = _df; + } else { xt = dvector(0, di-1); /* Vector for trial point */ + df = dvector(0, di-1); /* Vector for trial point */ + } /* -------------------------- */ /* First bracket the solution */ - LDBG(("linmin: Bracketing solution\n")) + LDBG((" linmind: Bracketing solution\n")) /* The line is measured as startpoint + offset * search vector. */ /* (Search isn't symetric, but it seems to depend on cp being */ @@ -442,7 +791,7 @@ void *fdata) /* Opaque data for func() */ xt[i] = cp[i] + xx * xi[i]; xf = (*func)(fdata, xt); - LDBG(("linmin: Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) + LDBG((" linmind: Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) /* Fix it so that we are decreasing from point a -> x */ if (xf > af) { @@ -450,33 +799,31 @@ void *fdata) /* Opaque data for func() */ tt = ax; ax = xx; xx = tt; tt = af; af = xf; xf = tt; } - LDBG(("linmin: Ordered Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) + LDBG((" linmind: Ordered Initial points a:%f:%f -> b:%f:%f\n",ax,af,xx,xf)) bx = xx + POWELL_GOLD * (xx-ax); /* Guess b beyond a -> x */ for (i = 0; i < di; i++) xt[i] = cp[i] + bx * xi[i]; bf = (*func)(fdata, xt); - LDBG(("linmin: Initial bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf)) + LDBG((" linmind: Initial bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf)) #ifdef SLOPE_SANITY_CHECK /* If we're not seeing a slope indicitive of progress */ /* of order ftol, give up straight away */ if (2000.0 * fabs(xf - bf) <= ftol * (fabs(xf) + fabs(bf)) && 2000.0 * fabs(af - xf) <= ftol * (fabs(af) + fabs(xf))) { - LDBG(("linmin: giving up because slope is too shallow\n")) - if (xt != XT) - free_dvector(xt,0,di-1); + LDBG((" linmind: giving up because slope is too shallow\n")) + if (di > 10) { + free_dvector(df, 0, di-1); + free_dvector(xt, 0, di-1); + } if (bf < xf) { xf = bf; xx = bx; } - - /* Compute solution vector */ - for (i = 0; i < di; i++) - cp[i] += xx * xi[i]; - return xf; + goto done; } #endif /* SLOPE_SANITY_CHECK */ @@ -485,8 +832,7 @@ void *fdata) /* Opaque data for func() */ double ulim, ux, uf; double tt, r, q; -// LDBG(("linmin: Not bracketed a:%f:%f x:%f%f b:%f:%f\n",ax,af,xx,xf,bx,bf)) - LDBG(("linmin: Not bracketed because xf %f > bf %f\n",xf, bf)) + LDBG((" linmind: Not bracketed because xf %f > bf %f\n",xf, bf)) LDBG((" ax = %f, xx = %f, bx = %f\n",ax,xx,bx)) /* Compute ux by parabolic interpolation from a, x & b */ @@ -559,7 +905,7 @@ void *fdata) /* Opaque data for func() */ bx = ux; bf = uf; //printf("~1 move along to the right (a<-x, x<-b, b- x -> b */ //printf("~1 got bracketed minimum at %f (%f), %f (%f), %f (%f)\n",ax,af,xx,xf,bx,bf); @@ -571,8 +917,10 @@ void *fdata) /* Opaque data for func() */ /* w is second best function value so far */ /* v is previous second best, or third best */ /* u is most recently tested point */ + double wx, vx, ux; /* Search vector multipliers */ double wf, vf = 0.0, uf; /* Function values at those points */ + double xd, wd, vd, ud; /* Derivative values at those points */ int iter; double de = 0.0; /* Distance moved on previous step */ double e = 0.0; /* Distance moved on 2nd previous step */ @@ -587,6 +935,16 @@ void *fdata) /* Opaque data for func() */ wx = vx = xx; /* Initial values of other center points */ wf = xf = xf; + /* Lookup derivative at x (we already have xf from bracketing) */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + xx * xi[i]; + (*dfunc)(fdata, df, xt); + for (xd = 0.0, i = 0; i < di; i++) + xd += xi[i] * df[i]; + wd = ud = xd; + + LDBG((" linmind: xx %f, deriv. xd %f\n",xx,xd)) + for (iter = 1; iter <= POWELL_MAXIT; iter++) { double mx = 0.5 * (ax + bx); /* m is center of bracket values */ #ifdef ABSTOL @@ -596,113 +954,169 @@ void *fdata) /* Opaque data for func() */ #endif double tol2 = 2.0 * tol1; - LDBG(("linmin: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf)) + LDBG((" linmind it %d: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",iter, ax,af,xx,xf,bx,bf)) /* See if we're done */ -//printf("~1 linmin check %f <= %f\n",fabs(xx - mx), tol2 - 0.5 * (bx - ax)); if (fabs(xx - mx) <= (tol2 - 0.5 * (bx - ax))) { - LDBG(("linmin: We're done because %f <= %f\n",fabs(xx - mx), tol2 - 0.5 * (bx - ax))) + LDBG((" linmind: We're done because %e <= %e\n",fabs(xx - mx), tol2 - 0.5 * (bx - ax))) break; } - if (fabs(e) > tol1) { /* Do a trial parabolic fit */ - double te, p, q, r; - r = (xx-wx) * (xf-vf); - q = (xx-vx) * (xf-wf); - p = (xx-vx) * q - (xx-wx) * r; - q = 2.0 * (q - r); - if (q > 0.0) - p = -p; - else - q = -q; + LDBG((" linmind: e %f tol2 %f\n",e,tol1)) + + if (fabs(e) > tol1) { /* Do a trial secant fit */ + double te; + double dx1, dx2; /* Secant extrapolation points */ + double ux1, ux2; + int ch1, ch2; + + LDBG((" linmind: Doing trial secant fit\n" )) + + dx2 = dx1 = 2.0 * (bx - ax); /* Default to values out of the ax..bx bracket */ + + /* Extrapolated points from last two points (secant method) */ + if (wd != xd) + dx1 = (wx - xx) * xd/(xd - wd); + if (vd != xd) + dx2 = (vx - xx) * xd/(xd - vd); + + ux1 = xx + dx1; + ux2 = xx + dx2; + + /* Check which one is reasonable */ + ch1 = (ax - ux1) * (ux1 - bx) > 0.0 && xd * dx1 < 0.0; + ch2 = (ax - ux2) * (ux2 - bx) > 0.0 && xd * dx2 < 0.0; + + LDBG((" linmind: Doing dx1 %f dx2 %f ux1 %f ux2 %f ch1 %d ch2 %d\n",dx1,dx2,ux1,ux2,ch1,ch2)) + te = e; /* Save previous e value */ e = de; /* Previous steps distance moved */ - LDBG(("linmin: Trial parabolic fit\n" )) + if (!ch1 && !ch2) + goto bisect; - if (fabs(p) >= fabs(0.5 * q * te) || p <= q * (ax-xx) || p >= q * (bx-xx)) { - /* Give up on the parabolic fit, and use the golden section search */ - e = ((xx >= mx) ? ax-xx : bx-xx); /* Override previous distance moved */ - de = POWELL_CGOLD * e; - LDBG(("linmin: Moving to golden section search\n" )) - } else { /* Use parabolic fit */ - de = p/q; /* Change in xb */ - ux = xx + de; /* Trial point according to parabolic fit */ - if ((ux - ax) < tol2 || (bx - ux) < tol2) { - if ((mx - xx) > 0.0) /* Don't use parabolic, use tol1 */ - de = tol1; /* tol1 is +ve */ - else - de = -tol1; - } - LDBG(("linmin: Using parabolic fit\n" )) + /* Use smallest or the one that's valid */ + if (ch1 && ch2) + de = fabs(dx1) < fabs(dx2) ? dx1 : dx2; + if (ch1) + de = dx1; + else if (ch2) + de = dx2; + + LDBG((" linmind: set de %f\n",de)) + + if (fabs(de) > fabs(0.5 * te)) { + LDBG((" linmind: abs(de) %f > abs(te/2 = %f)\n",fabs(de),fabs(0.5 * te))) + goto bisect; } - } else { /* Keep using the golden section search */ - e = ((xx >= mx) ? ax-xx : bx-xx); /* Override previous distance moved */ - de = POWELL_CGOLD * e; - LDBG(("linmin: Continuing golden section search\n" )) + + ux = xx + de; + + if ((ux - ax) < tol2 || (bx - ux) < tol2) { + if ((mx - xx) < 0.0) + de = -fabs(tol1); + else + de = fabs(tol1); + LDBG((" linmind: Set de to tol1 %f\n",de)) + } +#ifdef LDEBUG + else { + LDBG((" linmind: Using secant fit de %f\n",de)) + } +#endif + + /* else bisect picking side using sign of derivative */ + } else { + bisect: + e = (xd >= 0.0 ? ax - xx : bx -xx); + de = 0.5 * e; + LDBG((" linmind: Continuing bisection search de %f\n",de)) } - if (fabs(de) >= tol1) { /* If de moves as much as tol1 would */ + if (fabs(de) >= tol1) { /* If de moves as much as tol1 or more */ ux = xx + de; /* use it */ - LDBG(("linmin: ux = %f = xx %f + de %f\n",ux,xx,de)) + + /* Evaluate function */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); + + LDBG((" linmind: ux = %f = xx %f + de %f, uf %f\n",ux,xx,de,uf)) + } else { /* else move by tol1 in direction de */ if (de > 0.0) { ux = xx + tol1; - LDBG(("linmin: ux = %f = xx %f + tol1 %f\n",ux,xx,tol1)) + LDBG((" linmind: ux = %f = xx %f + tol1 %f\n",ux,xx,tol1)) } else { ux = xx - tol1; - LDBG(("linmin: ux = %f = xx %f - tol1 %f\n",ux,xx,tol1)) + LDBG((" linmind: ux = %f = xx %f - tol1 %f\n",ux,xx,tol1)) + } + /* Evaluate function */ + for (i = 0; i < di; i++) + xt[i] = cp[i] + ux * xi[i]; + uf = (*func)(fdata, xt); + + LDBG((" linmind: uf %f\n",uf)) + + if (uf > xf) { /* If tol1 step downhill takes us uphill, we're done */ + goto done; } } - /* Evaluate function */ - for (i = 0; i < di; i++) - xt[i] = cp[i] + ux * xi[i]; - uf = (*func)(fdata, xt); + /* Evaluate derivative at trial point */ + (*dfunc)(fdata, df, xt); + for (ud = 0.0, i = 0; i < di; i++) + ud += xi[i] * df[i]; + + LDBG((" linmind: ux %f, deriv. ud %f\n",ux,ud)) + /* Houskeeping: */ if (uf <= xf) { /* Found new best solution */ - LDBG(("linmin: found new best solution at %f val %f\n",ux,uf)) + LDBG((" linmind: found new best solution at %f val %f dval %f\n",ux,uf,ud)) if (ux >= xx) { - ax = xx; af = xf; /* New lower bracket */ + ax = xx; af = xf; /* New lower bracket */ } else { bx = xx; bf = xf; /* New upper bracket */ } - vx = wx; vf = wf; /* New previous 2nd best solution */ - wx = xx; wf = xf; /* New 2nd best solution from previous best */ - xx = ux; xf = uf; /* New best solution from latest */ - } else { /* Found a worse solution */ - LDBG(("linmin: found new worse solution at %f val %f\n",ux,uf)) - LDBG(("linmin: current best at %f val %f\n",xx,xf)) + vx = wx; vf = wf; vd = wd; /* New previous 2nd best solution */ + wx = xx; wf = xf; wd = xd; /* New 2nd best solution from previous best */ + xx = ux; xf = uf; xd = ud; /* New best solution from latest */ + + } else { /* Found a worse solution */ + LDBG((" linmind: found new worse solution at %f val %f dval %f\n",ux,uf,ud)) + LDBG((" linmind: current best at %f val %f dval %f\n",xx,xf,xd)) if (ux < xx) { ax = ux; af = uf; /* New lower bracket */ } else { bx = ux; bf = uf; /* New upper bracket */ } if (uf <= wf || wx == xx) { /* New 2nd best solution, or equal best */ - vx = wx; vf = wf; /* New previous 2nd best solution */ - wx = ux; wf = uf; /* New 2nd best from latest */ + vx = wx; vf = wf; vd = wd; /* New previous 2nd best solution */ + wx = ux; wf = uf; wd = ud; /* New 2nd best from latest */ } else if (uf <= vf || vx == xx || vx == wx) { /* New 3rd best, or equal 1st & 2nd */ - vx = ux; vf = uf; /* New previous 2nd best from latest */ + vx = ux; vf = uf; vd = ud; /* New previous 2nd best from latest */ } } - } - /* !!! should do something if iter > POWELL_MAXIT !!!! */ + } /* Next itter */ + + /* !!! should do something if iter > POWELL_MAXIT ??? */ /* Solution is at xx, xf */ + done:; + if (di > 10) { + free_dvector(df, 0, di-1); + free_dvector(xt, 0, di-1); + } /* Compute solution vector */ - LDBG(("linmin: computing soln from best at %f val %f\n",xx,xf)) + LDBG((" linmind: computing soln from best at %f val %f dval %f\n",xx,xf,xd)) for (i = 0; i < di; i++) cp[i] += xx * xi[i]; - } - if (xt != XT) - free_dvector(xt,0,di-1); -//printf("~~~ line minimizer returning %e\n",xf); + } /* Minimizer context */ + return xf; } #undef POWELL_GOLD -#undef POWELL_CGOLD #undef POWELL_MAXIT -/**************************************************/ -- cgit v1.2.3