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

/**
		Constraints and line minimisations for functions with gradient
		discontinuity at co-ordinate zeros (hyperplanes). 
		The ZeroCrossingSparsity class handles this by checking for
		zero crossings during the line search. If a gradient 
		discontinuity causes a change of sign in the gradient,
		that dimension is constrained and becomes inactive (drops
		out of minimisation). It can later be reactivated if
		the algorithm determines that the gradient no longer 
		changes sign at that point.
  */

public class ZeroCrossingSparsity extends Constraints
{
	public double	alphas[];
	public int		crossings[];
	public int		numcrossings;
	public int		clipped=0;
	private int		released, releasable;
	private double	XEPS=1e-6, GEPS=1e-4;
	private double	GTHRESH=1;

	public static Factory factory() {
		return new Factory() {
			public Constraints build(MinimiserBase minimiser) { 
				return new ZeroCrossingSparsity(minimiser);
			}
		};
	}

	public ZeroCrossingSparsity(MinimiserBase minimiser)
	{
		super(minimiser);
		crossings=new int[n];
		alphas=new double[n];

		minimiser.add(new VParameter("XEPS", new DoubleModel() {
				public void set(double t) { XEPS=t; }
				public double get() { return XEPS; }
			} )
		);

		minimiser.add(new VParameter("GEPS", new DoubleModel() {
				public void set(double t) { GEPS=t; }
				public double get() { return GEPS; }
			} )
		);

		try {
			final Object jump=Shell.get("ZeroCrossingSparsity.jump");
			if (jump instanceof Double) {
				GTHRESH=((Double)jump).doubleValue();
				Shell.print("Gradient jump="+GTHRESH);
			} else if (jump instanceof VDouble) {
				GTHRESH=((VDouble)jump).value;
				((VDouble)jump).addObserver( new Observer() {
					public void update(Observable o, Object arg) {
						GTHRESH=((VDouble)jump).value;
					}
				} );
				Shell.print("Gradient jump in VDouble="+GTHRESH);
			}
		} catch (Exception ex) {
			Shell.print("*** no jump object found! ***");
			Shell.print(ex.toString());
			ex.printStackTrace();
		}
	}

	public void setGThreshold(double gt) { GTHRESH=gt; }

	/** constrain all dimensions which are currently set
	    to zero.
	  */

	public void initialise()
	{
		int p=0, q=0;

		for (int i=0; i<n; i++) {
			if (S.P1.x[i]==0) inactive[p++]=i;
			else active[q++]=i;
		}
		m=q;
	}


	/** find all the zero crossings */
		
	private boolean check() 
	{
		// got to find sign swappers
		numcrossings=0;
		for (int j=0; j<m; j++) {
			int i=active[j];
			if ( (S.P1.x[i]>XEPS && S.P2.x[i]<XEPS)
			  || (S.P2.x[i]>-XEPS && S.P1.x[i]<-XEPS) ) {
				crossings[numcrossings]=i;
				alphas[numcrossings]=-S.P1.x[i]/S.h[i];
				numcrossings++;
			}
		}
		return numcrossings>0;
	}

	public void report() 
	{
		if (numcrossings>0) {
			Shell.trace("crossings: "+numcrossings);
			for (int i=0; i<numcrossings; i++) {
				Shell.trace("dim: "+crossings[i]+" alpha: "+alphas[i]);
			}
		}
		super.report();
	}
	

	/**
		This finds the constrained dimension with the
		largest gradient, and if that absolute value of
		that gradient is bigger than 1, then releases 
		that dimension so that it can take part in the 
		optimization
	 */

	public boolean release(Datum P)
	{
		if (m==n) return false; // all active

		int imax=-1;
		double grad=GTHRESH+GEPS;

		releasable=0;
		for (int k=0; k<(n-m); k++) {
			int j=inactive[k];
			if (Math.abs(P.g[j])>grad) {
				releasable++;
				grad=Math.abs(P.g[j]);
				imax=j;
			}
		}
		
		if (imax!=-1) {	
			activate(imax); released=imax; 
			return true; 
		}
		else return false;
	}

	public int getReleasedDimension() { return released; }
	public int getReleasableDimensions() { return releasable; }

	public void lineSearch(Condition convergence, CubicLineSearch ls, double beta0)
	{
		boolean zcdone=false;

		step(beta0);
		eval2();

		while (!convergence.test()) {

			// formulate any other direction change conditions

			if (S.P2.f<S.P1.f && S.P2.s<0) {
				double beta=ls.extrapolate();
				S.move();
				step(beta);
				eval2();
				zcdone=false;
			} else { 
				/* couple of options here:
					we could check all the zero crossings in one go.
					if we find one that needs constraining, then
					we're done.
					otherwise, we don't need to check zero crossings
					again unless we extrapolate

					other option is to backtrack one zero crossing
					at a time

					OR find zero crossing just after (or before) 
					expected minimum
					
				 */

				// get cubic step
				double interp = ls.interpolate();

				if (!zcdone) { check(); zcdone=true; }

				// try to pick a good step length
				// bearing in mind the zero crossings we have
				if (numcrossings>0) {

					// we are going to look at the alphas, and
					// find the smallest one that is between 
					// interp and current alpha

					int best=-1;
					double bestAlpha=S.alpha;

					for (int k=0; k<numcrossings; k++) {
						double beta=alphas[k];
						if (beta>=interp && beta<bestAlpha) {
							best=k; bestAlpha=beta;
						}
					}
	
					if (best!=-1) {
						int		j=crossings[best];

						// kth zero crossing of jth dimension
						step(bestAlpha);
						S.P2.x[j]=0;
						eval2();
						if (Math.abs(S.P2.g[j])<=GTHRESH-GEPS) {
							inactivate(j);
							// Shell.trace("- inactivating: "+j);
							return;

							// what if SEVERAL dimensions
							// end up close to zero:
							// should we constrain them all?
						} else {
							// Shell.trace("truncating");
							// what?
							// should we stick with this point, or
							// move to interp?
							continue;
						}
					} // else fall through to ordinary step
				} 

				// no crossings, take ordinary step
				step(interp);
				eval2();
			}
		}
	}
}