diff src/samer/maths/opt/ZeroCrossingSparsity.java @ 0:bf79fb79ee13

Initial Mercurial check in.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/samer/maths/opt/ZeroCrossingSparsity.java	Tue Jan 17 17:50:20 2012 +0000
@@ -0,0 +1,256 @@
+/*
+ *	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();
+			}
+		}
+	}
+}
+