view src/samer/maths/opt/PolynomialLineSearch.java @ 8:5e3cbbf173aa tip

Reorganise some more

/*
 *	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;
		}
	}
}