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

Initial Mercurial check in.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/samer/maths/opt/Positivity.java	Tue Jan 17 17:50:20 2012 +0000
@@ -0,0 +1,179 @@
+/*
+ *	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.*;
+
+/**
+	Cubic line search with positivity constraints
+	dimensions become inactive when they go to
+	zero and the gradient > GEPS
+ */
+ 
+public class Positivity extends Constraints
+{
+	double		XEPS=1e-6;
+	double		GEPS=1e-4;
+	int []		zerod;
+	int			released;
+
+	public Positivity(MinimiserBase minimiser)
+	{ 
+		super(minimiser);
+		zerod=new int[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; }
+			} )
+		);
+	}
+	
+	/** 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; S.P1.x[i]=0; }
+			else active[q++]=i;
+		}
+		m=q;
+	}
+
+	/**
+		This finds the constrained dimension with the
+		largest positive gradient, and if 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;
+
+		int imax=-1;
+		double grad=-GEPS; // tolerance
+
+		for (int k=0; k<(n-m); k++) {
+			int j=inactive[k];
+			if (P.g[j]<grad) {
+				grad=P.g[j];
+				imax=j;
+			}
+		}
+		
+		if (imax!=-1) {	activate(imax); released=imax; return true; }
+		else return false;
+	}
+
+	public int getReleasedDimension() { return released; }
+
+	/**
+			Take a step but make sure none of the coordinates
+			becomes negative. This function returns true if
+			any dimensions were inactivated - implying that
+			we should finish the line search. The alternative
+			is that even the clipped step was too long and
+			we need to interpolate some more.
+	  */
+
+	public boolean clipStep(double beta)
+	{
+
+		step(beta);
+
+		
+		// got to find sign swappers
+		int imin=-1;
+		double alphamin=beta; // ??
+			
+		for (int k=0; k<m; k++) {
+			int i=active[k];
+			if (S.P2.x[i]<=XEPS) {
+				double alpha=-S.P1.x[i]/S.h[i];
+				if (imin==-1 || alpha<alphamin) {
+					imin=i;
+					alphamin=alpha;
+				}
+			}
+		}
+
+		// if no negative values, evaluate and return		
+		if (imin==-1) { S.eval2(); return false; } 
+		
+		// atlease one dimension has become 
+		// zero or negative. imin was the first
+		// what if alphamin=0?
+		step(alphamin);
+			
+		// now we should check for ANY unconstrained
+		// x[i] being less than XEPS, and set it to
+		// exactly zero, including, obviously x[j]
+		int numzerod=0;
+		for (int k=0; k<m; k++) {
+			int j=active[k];
+			if (S.P2.x[j]<=XEPS) {
+				S.P2.x[j]=0;
+				zerod[numzerod++]=j;
+			}
+		}
+
+		eval2();
+
+		// next, check zeroed dimensions
+		// for constraint eligibilty
+		boolean flag=false;
+		
+		for (int k=0; k<numzerod; k++) {
+			int j=zerod[k];
+			if (S.P2.g[j]>=-GEPS) {
+				inactivate(j);
+				flag=true;
+			}
+		}
+
+		return flag;
+	}
+
+
+	public void lineSearch(Condition convergence, CubicLineSearch ls, double beta0)
+	{
+		if (clipStep(beta0)) return;
+
+		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();
+				if (clipStep(beta)) break;
+			} else { 
+				// since both P1 and P2 must satisfy positivity
+				// constraints, all points between them do alse,
+				// so no worries about interpolation
+				step(ls.interpolate());
+				eval2();
+			}
+		}
+	}
+}
+