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 ConstrainedGillMurray extends GillMurray
samer@0: {
samer@0: 	State		S;
samer@0: 	Constraints C;
samer@0: 	
samer@0: 	public ConstrainedGillMurray(State s, Constraints c) {
samer@0: 		super(s); S=s; C=c;
samer@0: 	}
samer@0: 
samer@0: 	public void computeHg(double [] g)
samer@0: 	{
samer@0: 		C.zeroInactive(S.h);
samer@0: 		C.mul(S.h,H,g);
samer@0: 		C.negate(S.h);
samer@0: 		S.normh = Util.maxabs(S.h);
samer@0: 	}
samer@0: 
samer@0: 	public void steepest() {
samer@0: 		C.negate(S.P1.g,S.h);
samer@0: 		C.zeroInactive(S.h);
samer@0: 		S.normh = Util.maxabs(S.h);
samer@0: 	}
samer@0: 
samer@0: 	public void resetDimension(int k) {
samer@0: 		double [] hk=H0.getArray()[k];
samer@0: 		for (int i=0; i<S.n; i++) H[i][k]=H[k][i]=hk[i];
samer@0: 	}
samer@0: 
samer@0: 	/* calculate new hessian and search direction */
samer@0: 	public void update()
samer@0: 	{
samer@0: 		// sparse subtraction and dot product
samer@0: 		C.sub(u,S.P2.x,S.P1.x);
samer@0: 		C.sub(v,S.P2.g,S.P1.g);
samer@0: 		double uv = C.dot(u,v);
samer@0: 
samer@0: 		if (uvCheck(uv,C.dot(u,u),C.dot(v,v))) {
samer@0: 			C.mul(Hv,H,v);
samer@0: 
samer@0: 			double a = 1/uv;
samer@0: 			double b = a*(1+a*C.dot(v,Hv));
samer@0: 
samer@0: 			for (int k=0; k<C.m; k++) {
samer@0: 				int i=C.active[k];
samer@0: 				for (int l=0; l<k; l++) {
samer@0: 					int j=C.active[l];
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(C.dot(S.P2.g,S.h));
samer@0: 	}
samer@0: }