view src/samer/maths/opt/Constraints.java @ 0:bf79fb79ee13

Initial Mercurial check in.
author samer
date Tue, 17 Jan 2012 17:50:20 +0000
parents
children
line wrap: on
line source
/*
 *	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.*;

public abstract class Constraints
{
	public State		S;
	public int[]		active;			// active dimensions
	public int[]		inactive;		// inactive dimensions
	public int			n,m;			// number of active dimensions

	public Constraints(State s)
	{
		S=s; n=S.n;
		active = new int[n];
		inactive = new int[n];
		activateAll();
	}

	public final int activeCount() { return m; }

	public void report() 
	{
		Shell.print("inactive dimensions:");
		for (int i=0; i<(n-m); i++) {
			Shell.print(" - "+inactive[i]);
		}
		Shell.print("active dimensions:");
		for (int i=0; i<m; i++) {
			Shell.print(" + "+active[i]);
		}
		Shell.print("  active: "+m);
		Shell.print("inactive: "+(n-m));
	}

	public void activateAll() {
		for (int i=0; i<n; i++) active[i]=i;
		m=n;
	}

	public void deactivateAll() { 
		for (int i=0; i<n; i++) inactive[i]=i;
		m=0; 
	}

	public void zeroInactive(double [] z) {
		for (int i=0; i<(n-m); i++) z[inactive[i]]=0;
	}

	/**
			Inactivate the i'th dimension
			Presupposes that it is not already inactive
	  */	
	public final void inactivate(int i) 
	{
		inactive[n-m]=i;
		int j=Util.find(i,active,m);
		if (j<(m-1)) {
			System.arraycopy(active,j+1,active,j,m-j-1); 
		}
		m--;
	}

	/**
			activate the i'th dimension
			Presupposes that it is not already inactive
	  */	
	public final void activate(int i) 
	{
		active[m]=i;
		int j=Util.find(i,inactive,n-m);
		if (j<(n-m-1)) {
			System.arraycopy(inactive,j+1,inactive,j,n-m-j-1); 
		}
		m++;
	}

	/**
			take a step of lambda*h, set alpha, but don't
			evaluate function at new point
	  */
	public void step(double lambda)
	{
		for (int i=0; i<m; i++) {
			int j=active[i];
			S.P2.x[j] = S.P1.x[j] + lambda*S.h[j];
		}
		S.alpha=lambda;
	}

	/**		evaluate function and gradient at P2 */
	public void eval2()
	{
		S.func.evaluate(S.P2);
		S.P2.s = dot(S.P2.g,S.h);
	}

	/**		sparse matrix mutiplication */
	public final void mul( double [] y, double [][]A, double [] x)
	{
		double [] arow;

		for (int k=0; k<m; k++) {
			int i=active[k];
			double t=0;
			arow=A[i];
			for (int l=0; l<m; l++) { int j=active[l]; t+=arow[j]*x[j]; }
			y[i]=t;
		}
	}

	/**		sparse vector dot product */
	public final double dot( double [] a, double [] b)
	{
		double t=0;
		for (int k=0; k<m; k++) { int i=active[k]; t+=a[i]*b[i]; }
		return t;
	}

	/**		sparse vector time scalar */
	public final void mul( double [] a, double b)
	{
		for (int k=0; k<m; k++) { int i=active[k]; a[i]*=b; }
	}

	/**		sparse vector subtraction */
	public final void sub( double [] out, double [] a, double [] b)
	{
		for (int k=0; k<m; k++) { int i=active[k]; out[i]=a[i]-b[i]; }
	}

	/**		sparse vector subtraction in place */
	public final void sub( double [] a, double [] b)
	{
		for (int k=0; k<m; k++) { int i=active[k]; a[i]-=b[i]; }
	}


	/**		sparse vector negation in place */
	public final void negate(double [] x) 
	{
		for (int k=0; k<m; k++) { int i=active[k]; x[i]=-x[i]; }
	}
	
	/**		sparse vector negation in place */
	public final void negate(double [] x, double [] y) 
	{
		for (int k=0; k<m; k++) { int i=active[k]; y[i]=-x[i]; }
	}
	
	/**		sparse max(abs(x)) */
	public final double maxabs(double [] x)
	{
		double max=0;
		for (int k=0; k<m; k++) {
			double xk = Math.abs(x[active[k]]);
			if (xk>max) max=xk;
		}
		return max;
	}

	// ugh
	public boolean release(Datum P) { return false; }
	public int getReleasedDimension() { return -1; }
	public int getReleasableDimensions() { return 0; }
	public abstract void lineSearch(Condition c, CubicLineSearch ls, double step);
	public abstract void initialise();

	public interface Factory {
		public Constraints build(MinimiserBase min);
	}
}