New upstream version 2.0.1+repackupstream/2.0.1+repackupstream

1 files changed, 563 insertions, 149 deletions

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-<u)\n");
+	}
+	LDBG((" linmin: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf))
+	/* Got bracketed minimum between a -> 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-<u)\n");
 	}
-	LDBG(("linmin: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf))
+	LDBG((" linmind: Got bracket a:%f:%f x:%f:%f b:%f:%f\n",ax,af,xx,xf,bx,bf))
 	/* Got bracketed minimum between a -> 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
 
-/**************************************************/