view src/samer/maths/opt/GillMurray.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.*;

public class GillMurray
{
	public Matrix		hessin;
	Jama.Matrix	H0;
	double[][]	H;				// approximate inverse Hessian
	double[]		u, v, Hv;	// used in updating hessian
	double		uv;
	State			S;
	boolean		uvc;

	public GillMurray(State s)
	{
		S = s;
		hessin = new Matrix("H",s.n,s.n);
		H = hessin.getArray();
		u = new double[s.n];
		v = new double[s.n];
		Hv= new double[s.n];
		uv= 0;

		H0=Matrix.identity(new Jama.Matrix(s.n,s.n));
		resetHessian();
	}

	public void dispose() { 
		H=null;
		u=v=Hv=null;
		hessin.dispose();
		S=null;
	}

	public void setInitialHessian(Jama.Matrix H0) { this.H0=H0; }
	public void resetHessian() { hessin.assign(H0); }
	public void init() { computeHg(S.P1.g); }

	public void steepest() {
		Mathx.negate(S.P1.g,S.h);
		S.normh = Util.maxabs(S.h);
	}

	public void computeHg(double [] g) {
		Mathx.mul(S.h,H,g);
		Mathx.negate(S.h);
		S.normh = Util.maxabs(S.h);
	}

	protected boolean uvCheck( double uv, double uu, double vv) {
		return uvc = (uv>1e-6*Math.sqrt(uu*vv));
	}

	public void update()
	{
		// calculate new hessian
		Mathx.sub(v,S.P2.g,S.P1.g);
 		Mathx.sub(u,S.P2.x,S.P1.x);
		double uv=Mathx.dot(u,v);

		if (uvCheck(uv,Mathx.norm2(u),Mathx.norm2(v))) {

			Mathx.mul(Hv,H,v);
			double a = 1/uv;
			double b = a*(1+a*Mathx.dot(v,Hv));

			for (int i=0; i<S.n; i++) {
				for (int j=0; j<i; j++) {
					H[i][j] = H[i][j] - a*(Hv[i]*u[j]+u[i]*Hv[j]) + b*u[i]*u[j]; 
					H[j][i] = H[i][j];
				}
				H[i][i] = H[i][i] - a*(2*Hv[i]*u[i]) + b*u[i]*u[i];
			}
		}
		computeHg(S.P2.g);
		S.setSlope2(Mathx.dot(S.P2.g,S.h));
	}
}