diff src/samer/maths/opt/PolynomialLineSearch.java @ 0:bf79fb79ee13

Initial Mercurial check in.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/samer/maths/opt/PolynomialLineSearch.java	Tue Jan 17 17:50:20 2012 +0000
@@ -0,0 +1,74 @@
+/*
+ *	Copyright (c) 2000, Samer Abdallah, King's College London.
+ *	All rights reserved.
+ *
+ *	This software is provided AS iS and WITHOUT ANY WARRANTY; 
+ *	without even the implied warranty of MERCHANTABILITY or 
+ *	FITNESS FOR A PARTICULAR PURPOSE.
+ */
+
+package samer.maths.opt;
+import  samer.maths.*;
+import  samer.core.*;
+
+/**
+		Line search from <b>Numerical Recipes in C</b>, pretty
+		much as supplied. 
+
+		I removed the max step and slope checking and made it the
+		responsibility of the calling code to check these.
+		Also, this works with ONE-DIMENSIONAL functions: the stuff
+		about vectors can be handled by a <code>LineFunction</code> object.
+ */
+
+public class PolynomialLineSearch
+{
+	State		S;
+
+	public PolynomialLineSearch(State s) { S=s; }
+
+	/**
+		This assumes that the initial step has already been taken
+		@param	alamin	minimum step size
+	  */
+
+	public void run(Condition termination)
+	{
+		double alam, alam2, f3;
+		double tmplam;
+
+		// if initial step is ok...
+		if (termination.test()) return;
+
+		// backtrack using quadratic interpolation
+		tmplam = -S.P1.s/(2*(S.P1.f-S.P2.f-S.P1.s));
+		alam = S.alpha;
+
+		for (;;) {
+
+			// remember failed parameters
+			alam2=alam;	f3=S.P2.f;
+			alam=Math.max(tmplam,0.1*alam);
+
+			S.step(alam);
+
+			// if step was ok...
+			if (termination.test()) return;
+
+			// backtrack using cubic interpolation
+			double rhs1 = S.P2.f - S.P1.f - alam*S.P1.s;
+			double rhs2 = f3     - S.P1.f - alam2*S.P1.s;
+			double a = (rhs1/(alam*alam) - rhs2/(alam2*alam2)) / (alam-alam2);
+			double b = (-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
+
+			if (a==0.0) tmplam = -S.P1.s/(2*b);
+			else {
+				double disc=b*b-2*a*S.P1.s;
+				if (disc<0) tmplam=alam/2;
+				else if (b<=0) tmplam=(-b+Math.sqrt(disc))/(3*a);
+				else tmplam=-S.P1.s/(b+Math.sqrt(disc));
+			}
+			if (tmplam>alam/2) tmplam=alam/2;
+		}
+	}
+}