view src/samer/maths/opt/ConstrainedMinimiser.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.*;
import  samer.core.types.*;
import  samer.core.util.heavy.*;

/**
		Constrained minimiser.

		-	ConjugateGradient 
		-	OR Quasi-newton (using GillMurray uydates)
		-	Safeguarded cubic interpolation line search using gradients

		We expect a Class object to be in object Space
		to tell us what kind of constraints to create
  */


public class ConstrainedMinimiser extends MinimiserBase implements Agent
{
	ConstrainedGillMurray	dir;
//	ConstrainedConjGrad		dir;
	Constraints					cons;	// handles line search
	VDouble						percentage;
	boolean toggle=false;

	/** expected stack top: Vec, Functionx, Class */
	public ConstrainedMinimiser(Vec v, Functionx f, Constraints.Factory constraintFactory) 
	{
		super(v,f);
		cons = constraintFactory.build(this);
		init();
	}

	public ConstrainedMinimiser(Vec v, Functionx f, Class conscl)
	{
		super(v,f);

		if (conscl==ZeroCrossingSparsity.class) {
			cons = new ZeroCrossingSparsity(this);
		} else if (conscl==Positivity.class) {
			cons = new Positivity(this);
		} else {
			throw new Error("Unrecognised constrain class: "+conscl);
		}
		init();
	}

	private void init() {
		dir = new ConstrainedGillMurray(this,cons);
//		dir = new ConstrainedConjGrad(this,cons);
		percentage=new VDouble("releaseFactor",0.01);

		// add( dir.hessin.viewable());
		add( percentage);
		cons.activateAll();
	}

	public void setHessian(Jama.Matrix H) { dir.setInitialHessian(H); }

	public void getCommands(Agent.Registry r) {
		super.getCommands(r);
		r.add("all").add("run");
		r.add("report").add("release").add("resetHess");
	}

	public void execute(String cmd, Environment env) throws Exception
	{
		if (cmd.equals("run"))  			run();
		else if (cmd.equals("report")) 	cons.report();
		else if (cmd.equals("release")) 	releaseDimensions(P1);
		else if (cmd.equals("all")) {
			toggle = X._bool(env.datum(),!toggle);
			if (toggle) cons.activateAll();
			else cons.deactivateAll();
		} else if (cmd.equals("resetHess")) {
			dir.resetHessian();
		} else super.execute(cmd,env);
	}

	private boolean releaseDimensions(Datum P)
	{
		if (cons.release(P)) {
			dir.resetDimension(cons.getReleasedDimension());

			int p=cons.getReleasableDimensions()-1;

			// going to release a percentage of releasable dims?
			// or absolute number?
			p=(int)(percentage.value*p);
			// if (p>1) p=1;
			for (; p>0; p--) {
				cons.release(P);
				dir.resetDimension(cons.getReleasedDimension());
			}

 			return true;
		}
		return false;
	}

	public void run()
	{
		double	beta;
		int		i, maxiter=getMaxiter();
		boolean	triedSteepest=false;

		// we start from the current state instead.
		// cons.initialise should set up the constraints
		// to reflect the current state.
		cons.initialise();

		evaluate();
		releaseDimensions(P1);

		if (cons.activeCount()==0) {
			iters.next(0);
			sig3.off();
			return;
		}


		sig1.off();

		// must reset hessian for active dimensions
		cons.zeroInactive(P2.x);
		dir.resetHessian();
		dir.hessin.changed();
		dir.init();
		setSlope();
		beta = initialStep();
		// beta = beta.value; ??

		for (i=0; i<maxiter; i++) {

			// constrained line search with initial step beta0
			lstest.init();
			cons.lineSearch(lstest,ls,beta);
			lsiters.next(lstest.count);
			steplength.next(alpha);

			if (lstest.tiny && (P2.f>=P1.f)) {
				//  tiny step was no good
				sig1.on();
				sig2.off();
				if (gconv.isSatisfied(P1.g,cons)) break;
				if (triedSteepest) break;
				sig2.on();

				beta = this.beta.value;
				dir.update();			// ??
				if (dir.uvc) dir.init();
				else		 dir.steepest();
				setSlope();
				triedSteepest=true;
				continue;
			}


			// check gradients at P2
			if (!releaseDimensions(P2)) {
				// don't quit if any dimensions released
				if (cons.activeCount()==0) break;
				if (xfconv.isSatisfied(this)) {
					if (gconv.isSatisfied(P2.g,cons)) break;
				}
			}

			// beta = beta.value;
			beta = nextStep();
			dir.update();
			move();
			perIteration();
		}

		perOptimisation(i);
		dir.hessin.changed();
	}
}