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

/**
	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();
			}
		}
	}
}