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