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