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  java.util.*;
samer@0: 
samer@0: /**
samer@0: 		Constraints and line minimisations for functions with gradient
samer@0: 		discontinuity at co-ordinate zeros (hyperplanes). 
samer@0: 		The ZeroCrossingSparsity class handles this by checking for
samer@0: 		zero crossings during the line search. If a gradient 
samer@0: 		discontinuity causes a change of sign in the gradient,
samer@0: 		that dimension is constrained and becomes inactive (drops
samer@0: 		out of minimisation). It can later be reactivated if
samer@0: 		the algorithm determines that the gradient no longer 
samer@0: 		changes sign at that point.
samer@0:   */
samer@0: 
samer@0: public class ZeroCrossingSparsity extends Constraints
samer@0: {
samer@0: 	public double	alphas[];
samer@0: 	public int		crossings[];
samer@0: 	public int		numcrossings;
samer@0: 	public int		clipped=0;
samer@0: 	private int		released, releasable;
samer@0: 	private double	XEPS=1e-6, GEPS=1e-4;
samer@0: 	private double	GTHRESH=1;
samer@0: 
samer@0: 	public static Factory factory() {
samer@0: 		return new Factory() {
samer@0: 			public Constraints build(MinimiserBase minimiser) { 
samer@0: 				return new ZeroCrossingSparsity(minimiser);
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	public ZeroCrossingSparsity(MinimiserBase minimiser)
samer@0: 	{
samer@0: 		super(minimiser);
samer@0: 		crossings=new int[n];
samer@0: 		alphas=new double[n];
samer@0: 
samer@0: 		minimiser.add(new VParameter("XEPS", new DoubleModel() {
samer@0: 				public void set(double t) { XEPS=t; }
samer@0: 				public double get() { return XEPS; }
samer@0: 			} )
samer@0: 		);
samer@0: 
samer@0: 		minimiser.add(new VParameter("GEPS", new DoubleModel() {
samer@0: 				public void set(double t) { GEPS=t; }
samer@0: 				public double get() { return GEPS; }
samer@0: 			} )
samer@0: 		);
samer@0: 
samer@0: 		try {
samer@0: 			final Object jump=Shell.get("ZeroCrossingSparsity.jump");
samer@0: 			if (jump instanceof Double) {
samer@0: 				GTHRESH=((Double)jump).doubleValue();
samer@0: 				Shell.print("Gradient jump="+GTHRESH);
samer@0: 			} else if (jump instanceof VDouble) {
samer@0: 				GTHRESH=((VDouble)jump).value;
samer@0: 				((VDouble)jump).addObserver( new Observer() {
samer@0: 					public void update(Observable o, Object arg) {
samer@0: 						GTHRESH=((VDouble)jump).value;
samer@0: 					}
samer@0: 				} );
samer@0: 				Shell.print("Gradient jump in VDouble="+GTHRESH);
samer@0: 			}
samer@0: 		} catch (Exception ex) {
samer@0: 			Shell.print("*** no jump object found! ***");
samer@0: 			Shell.print(ex.toString());
samer@0: 			ex.printStackTrace();
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	public void setGThreshold(double gt) { GTHRESH=gt; }
samer@0: 
samer@0: 	/** constrain all dimensions which are currently set
samer@0: 	    to zero.
samer@0: 	  */
samer@0: 
samer@0: 	public void initialise()
samer@0: 	{
samer@0: 		int p=0, q=0;
samer@0: 
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			if (S.P1.x[i]==0) inactive[p++]=i;
samer@0: 			else active[q++]=i;
samer@0: 		}
samer@0: 		m=q;
samer@0: 	}
samer@0: 
samer@0: 
samer@0: 	/** find all the zero crossings */
samer@0: 		
samer@0: 	private boolean check() 
samer@0: 	{
samer@0: 		// got to find sign swappers
samer@0: 		numcrossings=0;
samer@0: 		for (int j=0; j<m; j++) {
samer@0: 			int i=active[j];
samer@0: 			if ( (S.P1.x[i]>XEPS && S.P2.x[i]<XEPS)
samer@0: 			  || (S.P2.x[i]>-XEPS && S.P1.x[i]<-XEPS) ) {
samer@0: 				crossings[numcrossings]=i;
samer@0: 				alphas[numcrossings]=-S.P1.x[i]/S.h[i];
samer@0: 				numcrossings++;
samer@0: 			}
samer@0: 		}
samer@0: 		return numcrossings>0;
samer@0: 	}
samer@0: 
samer@0: 	public void report() 
samer@0: 	{
samer@0: 		if (numcrossings>0) {
samer@0: 			Shell.trace("crossings: "+numcrossings);
samer@0: 			for (int i=0; i<numcrossings; i++) {
samer@0: 				Shell.trace("dim: "+crossings[i]+" alpha: "+alphas[i]);
samer@0: 			}
samer@0: 		}
samer@0: 		super.report();
samer@0: 	}
samer@0: 	
samer@0: 
samer@0: 	/**
samer@0: 		This finds the constrained dimension with the
samer@0: 		largest gradient, and if that absolute value of
samer@0: 		that gradient is bigger than 1, then releases 
samer@0: 		that dimension so that it can take part in the 
samer@0: 		optimization
samer@0: 	 */
samer@0: 
samer@0: 	public boolean release(Datum P)
samer@0: 	{
samer@0: 		if (m==n) return false; // all active
samer@0: 
samer@0: 		int imax=-1;
samer@0: 		double grad=GTHRESH+GEPS;
samer@0: 
samer@0: 		releasable=0;
samer@0: 		for (int k=0; k<(n-m); k++) {
samer@0: 			int j=inactive[k];
samer@0: 			if (Math.abs(P.g[j])>grad) {
samer@0: 				releasable++;
samer@0: 				grad=Math.abs(P.g[j]);
samer@0: 				imax=j;
samer@0: 			}
samer@0: 		}
samer@0: 		
samer@0: 		if (imax!=-1) {	
samer@0: 			activate(imax); released=imax; 
samer@0: 			return true; 
samer@0: 		}
samer@0: 		else return false;
samer@0: 	}
samer@0: 
samer@0: 	public int getReleasedDimension() { return released; }
samer@0: 	public int getReleasableDimensions() { return releasable; }
samer@0: 
samer@0: 	public void lineSearch(Condition convergence, CubicLineSearch ls, double beta0)
samer@0: 	{
samer@0: 		boolean zcdone=false;
samer@0: 
samer@0: 		step(beta0);
samer@0: 		eval2();
samer@0: 
samer@0: 		while (!convergence.test()) {
samer@0: 
samer@0: 			// formulate any other direction change conditions
samer@0: 
samer@0: 			if (S.P2.f<S.P1.f && S.P2.s<0) {
samer@0: 				double beta=ls.extrapolate();
samer@0: 				S.move();
samer@0: 				step(beta);
samer@0: 				eval2();
samer@0: 				zcdone=false;
samer@0: 			} else { 
samer@0: 				/* couple of options here:
samer@0: 					we could check all the zero crossings in one go.
samer@0: 					if we find one that needs constraining, then
samer@0: 					we're done.
samer@0: 					otherwise, we don't need to check zero crossings
samer@0: 					again unless we extrapolate
samer@0: 
samer@0: 					other option is to backtrack one zero crossing
samer@0: 					at a time
samer@0: 
samer@0: 					OR find zero crossing just after (or before) 
samer@0: 					expected minimum
samer@0: 					
samer@0: 				 */
samer@0: 
samer@0: 				// get cubic step
samer@0: 				double interp = ls.interpolate();
samer@0: 
samer@0: 				if (!zcdone) { check(); zcdone=true; }
samer@0: 
samer@0: 				// try to pick a good step length
samer@0: 				// bearing in mind the zero crossings we have
samer@0: 				if (numcrossings>0) {
samer@0: 
samer@0: 					// we are going to look at the alphas, and
samer@0: 					// find the smallest one that is between 
samer@0: 					// interp and current alpha
samer@0: 
samer@0: 					int best=-1;
samer@0: 					double bestAlpha=S.alpha;
samer@0: 
samer@0: 					for (int k=0; k<numcrossings; k++) {
samer@0: 						double beta=alphas[k];
samer@0: 						if (beta>=interp && beta<bestAlpha) {
samer@0: 							best=k; bestAlpha=beta;
samer@0: 						}
samer@0: 					}
samer@0: 	
samer@0: 					if (best!=-1) {
samer@0: 						int		j=crossings[best];
samer@0: 
samer@0: 						// kth zero crossing of jth dimension
samer@0: 						step(bestAlpha);
samer@0: 						S.P2.x[j]=0;
samer@0: 						eval2();
samer@0: 						if (Math.abs(S.P2.g[j])<=GTHRESH-GEPS) {
samer@0: 							inactivate(j);
samer@0: 							// Shell.trace("- inactivating: "+j);
samer@0: 							return;
samer@0: 
samer@0: 							// what if SEVERAL dimensions
samer@0: 							// end up close to zero:
samer@0: 							// should we constrain them all?
samer@0: 						} else {
samer@0: 							// Shell.trace("truncating");
samer@0: 							// what?
samer@0: 							// should we stick with this point, or
samer@0: 							// move to interp?
samer@0: 							continue;
samer@0: 						}
samer@0: 					} // else fall through to ordinary step
samer@0: 				} 
samer@0: 
samer@0: 				// no crossings, take ordinary step
samer@0: 				step(interp);
samer@0: 				eval2();
samer@0: 			}
samer@0: 		}
samer@0: 	}
samer@0: }
samer@0: