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: public class GillMurray
samer@0: {
samer@0: 	public Matrix		hessin;
samer@0: 	Jama.Matrix	H0;
samer@0: 	double[][]	H;				// approximate inverse Hessian
samer@0: 	double[]		u, v, Hv;	// used in updating hessian
samer@0: 	double		uv;
samer@0: 	State			S;
samer@0: 	boolean		uvc;
samer@0: 
samer@0: 	public GillMurray(State s)
samer@0: 	{
samer@0: 		S = s;
samer@0: 		hessin = new Matrix("H",s.n,s.n);
samer@0: 		H = hessin.getArray();
samer@0: 		u = new double[s.n];
samer@0: 		v = new double[s.n];
samer@0: 		Hv= new double[s.n];
samer@0: 		uv= 0;
samer@0: 
samer@0: 		H0=Matrix.identity(new Jama.Matrix(s.n,s.n));
samer@0: 		resetHessian();
samer@0: 	}
samer@0: 
samer@0: 	public void dispose() { 
samer@0: 		H=null;
samer@0: 		u=v=Hv=null;
samer@0: 		hessin.dispose();
samer@0: 		S=null;
samer@0: 	}
samer@0: 
samer@0: 	public void setInitialHessian(Jama.Matrix H0) { this.H0=H0; }
samer@0: 	public void resetHessian() { hessin.assign(H0); }
samer@0: 	public void init() { computeHg(S.P1.g); }
samer@0: 
samer@0: 	public void steepest() {
samer@0: 		Mathx.negate(S.P1.g,S.h);
samer@0: 		S.normh = Util.maxabs(S.h);
samer@0: 	}
samer@0: 
samer@0: 	public void computeHg(double [] g) {
samer@0: 		Mathx.mul(S.h,H,g);
samer@0: 		Mathx.negate(S.h);
samer@0: 		S.normh = Util.maxabs(S.h);
samer@0: 	}
samer@0: 
samer@0: 	protected boolean uvCheck( double uv, double uu, double vv) {
samer@0: 		return uvc = (uv>1e-6*Math.sqrt(uu*vv));
samer@0: 	}
samer@0: 
samer@0: 	public void update()
samer@0: 	{
samer@0: 		// calculate new hessian
samer@0: 		Mathx.sub(v,S.P2.g,S.P1.g);
samer@0:  		Mathx.sub(u,S.P2.x,S.P1.x);
samer@0: 		double uv=Mathx.dot(u,v);
samer@0: 
samer@0: 		if (uvCheck(uv,Mathx.norm2(u),Mathx.norm2(v))) {
samer@0: 
samer@0: 			Mathx.mul(Hv,H,v);
samer@0: 			double a = 1/uv;
samer@0: 			double b = a*(1+a*Mathx.dot(v,Hv));
samer@0: 
samer@0: 			for (int i=0; i<S.n; i++) {
samer@0: 				for (int j=0; j<i; j++) {
samer@0: 					H[i][j] = H[i][j] - a*(Hv[i]*u[j]+u[i]*Hv[j]) + b*u[i]*u[j]; 
samer@0: 					H[j][i] = H[i][j];
samer@0: 				}
samer@0: 				H[i][i] = H[i][i] - a*(2*Hv[i]*u[i]) + b*u[i]*u[i];
samer@0: 			}
samer@0: 		}
samer@0: 		computeHg(S.P2.g);
samer@0: 		S.setSlope2(Mathx.dot(S.P2.g,S.h));
samer@0: 	}
samer@0: }