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: import  samer.core.types.*;
samer@0: import  samer.core.util.heavy.*;
samer@0: 
samer@0: /**
samer@0: 		Constrained minimiser.
samer@0: 
samer@0: 		-	ConjugateGradient 
samer@0: 		-	OR Quasi-newton (using GillMurray uydates)
samer@0: 		-	Safeguarded cubic interpolation line search using gradients
samer@0: 
samer@0: 		We expect a Class object to be in object Space
samer@0: 		to tell us what kind of constraints to create
samer@0:   */
samer@0: 
samer@0: 
samer@0: public class ConstrainedMinimiser extends MinimiserBase implements Agent
samer@0: {
samer@0: 	ConstrainedGillMurray	dir;
samer@0: //	ConstrainedConjGrad		dir;
samer@0: 	Constraints					cons;	// handles line search
samer@0: 	VDouble						percentage;
samer@0: 	boolean toggle=false;
samer@0: 
samer@0: 	/** expected stack top: Vec, Functionx, Class */
samer@0: 	public ConstrainedMinimiser(Vec v, Functionx f, Constraints.Factory constraintFactory) 
samer@0: 	{
samer@0: 		super(v,f);
samer@0: 		cons = constraintFactory.build(this);
samer@0: 		init();
samer@0: 	}
samer@0: 
samer@0: 	public ConstrainedMinimiser(Vec v, Functionx f, Class conscl)
samer@0: 	{
samer@0: 		super(v,f);
samer@0: 
samer@0: 		if (conscl==ZeroCrossingSparsity.class) {
samer@0: 			cons = new ZeroCrossingSparsity(this);
samer@0: 		} else if (conscl==Positivity.class) {
samer@0: 			cons = new Positivity(this);
samer@0: 		} else {
samer@0: 			throw new Error("Unrecognised constrain class: "+conscl);
samer@0: 		}
samer@0: 		init();
samer@0: 	}
samer@0: 
samer@0: 	private void init() {
samer@0: 		dir = new ConstrainedGillMurray(this,cons);
samer@0: //		dir = new ConstrainedConjGrad(this,cons);
samer@0: 		percentage=new VDouble("releaseFactor",0.01);
samer@0: 
samer@0: 		// add( dir.hessin.viewable());
samer@0: 		add( percentage);
samer@0: 		cons.activateAll();
samer@0: 	}
samer@0: 
samer@0: 	public void setHessian(Jama.Matrix H) { dir.setInitialHessian(H); }
samer@0: 
samer@0: 	public void getCommands(Agent.Registry r) {
samer@0: 		super.getCommands(r);
samer@0: 		r.add("all").add("run");
samer@0: 		r.add("report").add("release").add("resetHess");
samer@0: 	}
samer@0: 
samer@0: 	public void execute(String cmd, Environment env) throws Exception
samer@0: 	{
samer@0: 		if (cmd.equals("run"))  			run();
samer@0: 		else if (cmd.equals("report")) 	cons.report();
samer@0: 		else if (cmd.equals("release")) 	releaseDimensions(P1);
samer@0: 		else if (cmd.equals("all")) {
samer@0: 			toggle = X._bool(env.datum(),!toggle);
samer@0: 			if (toggle) cons.activateAll();
samer@0: 			else cons.deactivateAll();
samer@0: 		} else if (cmd.equals("resetHess")) {
samer@0: 			dir.resetHessian();
samer@0: 		} else super.execute(cmd,env);
samer@0: 	}
samer@0: 
samer@0: 	private boolean releaseDimensions(Datum P)
samer@0: 	{
samer@0: 		if (cons.release(P)) {
samer@0: 			dir.resetDimension(cons.getReleasedDimension());
samer@0: 
samer@0: 			int p=cons.getReleasableDimensions()-1;
samer@0: 
samer@0: 			// going to release a percentage of releasable dims?
samer@0: 			// or absolute number?
samer@0: 			p=(int)(percentage.value*p);
samer@0: 			// if (p>1) p=1;
samer@0: 			for (; p>0; p--) {
samer@0: 				cons.release(P);
samer@0: 				dir.resetDimension(cons.getReleasedDimension());
samer@0: 			}
samer@0: 
samer@0:  			return true;
samer@0: 		}
samer@0: 		return false;
samer@0: 	}
samer@0: 
samer@0: 	public void run()
samer@0: 	{
samer@0: 		double	beta;
samer@0: 		int		i, maxiter=getMaxiter();
samer@0: 		boolean	triedSteepest=false;
samer@0: 
samer@0: 		// we start from the current state instead.
samer@0: 		// cons.initialise should set up the constraints
samer@0: 		// to reflect the current state.
samer@0: 		cons.initialise();
samer@0: 
samer@0: 		evaluate();
samer@0: 		releaseDimensions(P1);
samer@0: 
samer@0: 		if (cons.activeCount()==0) {
samer@0: 			iters.next(0);
samer@0: 			sig3.off();
samer@0: 			return;
samer@0: 		}
samer@0: 
samer@0: 
samer@0: 		sig1.off();
samer@0: 
samer@0: 		// must reset hessian for active dimensions
samer@0: 		cons.zeroInactive(P2.x);
samer@0: 		dir.resetHessian();
samer@0: 		dir.hessin.changed();
samer@0: 		dir.init();
samer@0: 		setSlope();
samer@0: 		beta = initialStep();
samer@0: 		// beta = beta.value; ??
samer@0: 
samer@0: 		for (i=0; i<maxiter; i++) {
samer@0: 
samer@0: 			// constrained line search with initial step beta0
samer@0: 			lstest.init();
samer@0: 			cons.lineSearch(lstest,ls,beta);
samer@0: 			lsiters.next(lstest.count);
samer@0: 			steplength.next(alpha);
samer@0: 
samer@0: 			if (lstest.tiny && (P2.f>=P1.f)) {
samer@0: 				//  tiny step was no good
samer@0: 				sig1.on();
samer@0: 				sig2.off();
samer@0: 				if (gconv.isSatisfied(P1.g,cons)) break;
samer@0: 				if (triedSteepest) break;
samer@0: 				sig2.on();
samer@0: 
samer@0: 				beta = this.beta.value;
samer@0: 				dir.update();			// ??
samer@0: 				if (dir.uvc) dir.init();
samer@0: 				else		 dir.steepest();
samer@0: 				setSlope();
samer@0: 				triedSteepest=true;
samer@0: 				continue;
samer@0: 			}
samer@0: 
samer@0: 
samer@0: 			// check gradients at P2
samer@0: 			if (!releaseDimensions(P2)) {
samer@0: 				// don't quit if any dimensions released
samer@0: 				if (cons.activeCount()==0) break;
samer@0: 				if (xfconv.isSatisfied(this)) {
samer@0: 					if (gconv.isSatisfied(P2.g,cons)) break;
samer@0: 				}
samer@0: 			}
samer@0: 
samer@0: 			// beta = beta.value;
samer@0: 			beta = nextStep();
samer@0: 			dir.update();
samer@0: 			move();
samer@0: 			perIteration();
samer@0: 		}
samer@0: 
samer@0: 		perOptimisation(i);
samer@0: 		dir.hessin.changed();
samer@0: 	}
samer@0: }
samer@0: 
samer@0: 
samer@0: 
samer@0: