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: 
samer@0: public abstract class Constraints
samer@0: {
samer@0: 	public State		S;
samer@0: 	public int[]		active;			// active dimensions
samer@0: 	public int[]		inactive;		// inactive dimensions
samer@0: 	public int			n,m;			// number of active dimensions
samer@0: 
samer@0: 	public Constraints(State s)
samer@0: 	{
samer@0: 		S=s; n=S.n;
samer@0: 		active = new int[n];
samer@0: 		inactive = new int[n];
samer@0: 		activateAll();
samer@0: 	}
samer@0: 
samer@0: 	public final int activeCount() { return m; }
samer@0: 
samer@0: 	public void report() 
samer@0: 	{
samer@0: 		Shell.print("inactive dimensions:");
samer@0: 		for (int i=0; i<(n-m); i++) {
samer@0: 			Shell.print(" - "+inactive[i]);
samer@0: 		}
samer@0: 		Shell.print("active dimensions:");
samer@0: 		for (int i=0; i<m; i++) {
samer@0: 			Shell.print(" + "+active[i]);
samer@0: 		}
samer@0: 		Shell.print("  active: "+m);
samer@0: 		Shell.print("inactive: "+(n-m));
samer@0: 	}
samer@0: 
samer@0: 	public void activateAll() {
samer@0: 		for (int i=0; i<n; i++) active[i]=i;
samer@0: 		m=n;
samer@0: 	}
samer@0: 
samer@0: 	public void deactivateAll() { 
samer@0: 		for (int i=0; i<n; i++) inactive[i]=i;
samer@0: 		m=0; 
samer@0: 	}
samer@0: 
samer@0: 	public void zeroInactive(double [] z) {
samer@0: 		for (int i=0; i<(n-m); i++) z[inactive[i]]=0;
samer@0: 	}
samer@0: 
samer@0: 	/**
samer@0: 			Inactivate the i'th dimension
samer@0: 			Presupposes that it is not already inactive
samer@0: 	  */	
samer@0: 	public final void inactivate(int i) 
samer@0: 	{
samer@0: 		inactive[n-m]=i;
samer@0: 		int j=Util.find(i,active,m);
samer@0: 		if (j<(m-1)) {
samer@0: 			System.arraycopy(active,j+1,active,j,m-j-1); 
samer@0: 		}
samer@0: 		m--;
samer@0: 	}
samer@0: 
samer@0: 	/**
samer@0: 			activate the i'th dimension
samer@0: 			Presupposes that it is not already inactive
samer@0: 	  */	
samer@0: 	public final void activate(int i) 
samer@0: 	{
samer@0: 		active[m]=i;
samer@0: 		int j=Util.find(i,inactive,n-m);
samer@0: 		if (j<(n-m-1)) {
samer@0: 			System.arraycopy(inactive,j+1,inactive,j,n-m-j-1); 
samer@0: 		}
samer@0: 		m++;
samer@0: 	}
samer@0: 
samer@0: 	/**
samer@0: 			take a step of lambda*h, set alpha, but don't
samer@0: 			evaluate function at new point
samer@0: 	  */
samer@0: 	public void step(double lambda)
samer@0: 	{
samer@0: 		for (int i=0; i<m; i++) {
samer@0: 			int j=active[i];
samer@0: 			S.P2.x[j] = S.P1.x[j] + lambda*S.h[j];
samer@0: 		}
samer@0: 		S.alpha=lambda;
samer@0: 	}
samer@0: 
samer@0: 	/**		evaluate function and gradient at P2 */
samer@0: 	public void eval2()
samer@0: 	{
samer@0: 		S.func.evaluate(S.P2);
samer@0: 		S.P2.s = dot(S.P2.g,S.h);
samer@0: 	}
samer@0: 
samer@0: 	/**		sparse matrix mutiplication */
samer@0: 	public final void mul( double [] y, double [][]A, double [] x)
samer@0: 	{
samer@0: 		double [] arow;
samer@0: 
samer@0: 		for (int k=0; k<m; k++) {
samer@0: 			int i=active[k];
samer@0: 			double t=0;
samer@0: 			arow=A[i];
samer@0: 			for (int l=0; l<m; l++) { int j=active[l]; t+=arow[j]*x[j]; }
samer@0: 			y[i]=t;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	/**		sparse vector dot product */
samer@0: 	public final double dot( double [] a, double [] b)
samer@0: 	{
samer@0: 		double t=0;
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; t+=a[i]*b[i]; }
samer@0: 		return t;
samer@0: 	}
samer@0: 
samer@0: 	/**		sparse vector time scalar */
samer@0: 	public final void mul( double [] a, double b)
samer@0: 	{
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; a[i]*=b; }
samer@0: 	}
samer@0: 
samer@0: 	/**		sparse vector subtraction */
samer@0: 	public final void sub( double [] out, double [] a, double [] b)
samer@0: 	{
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; out[i]=a[i]-b[i]; }
samer@0: 	}
samer@0: 
samer@0: 	/**		sparse vector subtraction in place */
samer@0: 	public final void sub( double [] a, double [] b)
samer@0: 	{
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; a[i]-=b[i]; }
samer@0: 	}
samer@0: 
samer@0: 
samer@0: 	/**		sparse vector negation in place */
samer@0: 	public final void negate(double [] x) 
samer@0: 	{
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; x[i]=-x[i]; }
samer@0: 	}
samer@0: 	
samer@0: 	/**		sparse vector negation in place */
samer@0: 	public final void negate(double [] x, double [] y) 
samer@0: 	{
samer@0: 		for (int k=0; k<m; k++) { int i=active[k]; y[i]=-x[i]; }
samer@0: 	}
samer@0: 	
samer@0: 	/**		sparse max(abs(x)) */
samer@0: 	public final double maxabs(double [] x)
samer@0: 	{
samer@0: 		double max=0;
samer@0: 		for (int k=0; k<m; k++) {
samer@0: 			double xk = Math.abs(x[active[k]]);
samer@0: 			if (xk>max) max=xk;
samer@0: 		}
samer@0: 		return max;
samer@0: 	}
samer@0: 
samer@0: 	// ugh
samer@0: 	public boolean release(Datum P) { return false; }
samer@0: 	public int getReleasedDimension() { return -1; }
samer@0: 	public int getReleasableDimensions() { return 0; }
samer@0: 	public abstract void lineSearch(Condition c, CubicLineSearch ls, double step);
samer@0: 	public abstract void initialise();
samer@0: 
samer@0: 	public interface Factory {
samer@0: 		public Constraints build(MinimiserBase min);
samer@0: 	}
samer@0: }