samer@0: /*
samer@0:  *	Matrix.java	
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;
samer@0: import	samer.core.*;
samer@0: import   samer.core.util.*;
samer@0: import	java.util.*;
samer@0: import	java.io.*;
samer@0: 
samer@0: /**
samer@0: 	Represents a bundling of a Jama.Matrix with Viewable bits.
samer@0: 	Create one of these and you can either keep it invisible,
samer@0: 	expose its UI yourself, or let the Viewable framework
samer@0: 	expose a default viewer automatically
samer@0:  */
samer@0: 
samer@0: public class Matrix extends Jama.Matrix implements Mat, Saveable
samer@0: {
samer@0: 	Viewable	viewable;
samer@0: 
samer@0: 	// ............... Viewable bits .................
samer@0: 
samer@0: 	private class Vbl extends Viewable {
samer@0: 		Vbl(String nm) { super(nm); }
samer@0: 		Vbl() { super("matrix"); }
samer@0: 
samer@0: 		/** returns current Agent if there is one, creates a MatrixAgent if not */
samer@0: 		public Agent getAgent() {
samer@0: 			if (super.getAgent()==null) setAgent(new MatrixAgent(Matrix.this));
samer@0: 			return super.getAgent();
samer@0: 		}
samer@0: 
samer@0: 		public Matrix getMatrix() { return Matrix.this; }
samer@0: 		public Viewer getViewer()
samer@0: 		{
samer@0: 			BaseViewer		vwr = new BaseViewer(this);
samer@0: 			StringBuffer	buf = new StringBuffer();
samer@0: 
samer@0: 			buf.append(getLabel())
samer@0: 				.append(" (").append(getRowDimension())
samer@0: 				.append(" by ").append(getColumnDimension())
samer@0: 				.append(")");
samer@0: 			vwr.setText(buf.toString());
samer@0: 
samer@0: 			return vwr;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	public Observable	observable() { return viewable; }
samer@0: 	public Viewable	viewable() { return viewable; }
samer@0: 	public boolean		equals(Observable o) { return viewable==o; }
samer@0: 	public Node		getNode() { return viewable.getNode(); }
samer@0: 
samer@0: 	public final void	changed(Object o)	{ viewable.changed(o); }
samer@0: 	public final void	changed() { viewable.changed(); }
samer@0: 	public void	addObserver(Observer o) { viewable.addObserver(o); }
samer@0: 	public void	deleteObserver(Observer o) { viewable.deleteObserver(o); }
samer@0: 	public void dispose() { viewable.dispose(); }
samer@0: 
samer@0: 	public void copyFrom(double[][] src) {
samer@0: 		double[][] dest=getArray();
samer@0: 		for (int i=0; i<src.length; i++) {
samer@0: 			Mathx.copy(src[i],dest[i]);
samer@0: 			changed();
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	// ................................................
samer@0: 
samer@0: 	public Matrix( String name, int m, int n) { this(name,m,n,1); }
samer@0: 	public Matrix( String name, int m, int n, int flags) {
samer@0: 		super(m,n); viewable=new Vbl(name);
samer@0: 		if (flags==1) Shell.registerViewable(viewable);
samer@0: 	}
samer@0: 
samer@0: 
samer@0: 	public Matrix( int m, int n) { super(m,n); viewable=new Vbl(); Shell.registerViewable(viewable); }
samer@0: 	public Matrix( double array[][]) { super(array); viewable=new Vbl(); Shell.registerViewable(viewable); }
samer@0: 
samer@0: 	public String toString() {
samer@0: 		StringBuffer	buf = new StringBuffer();
samer@0: 
samer@0: 		buf.append(viewable.getLabel())
samer@0: 			.append(" (").append(getRowDimension())
samer@0: 			.append(" by ").append(getColumnDimension())
samer@0: 			.append(")");
samer@0: 		return buf.toString();
samer@0: 	}
samer@0: 
samer@0: 	// implement Mat interface
samer@0: 	public int width() { return getColumnDimension(); }
samer@0: 	public int height() { return getRowDimension(); }
samer@0: 
samer@0: 	public Agent getAgent() { return viewable.getAgent(); }
samer@0: 
samer@0: 
samer@0: 	// ................................................
samer@0: 	// some useful additions to Matrix functionality
samer@0: 
samer@0: 	public void set( Generator g) {
samer@0: 		for (int i=0; i<m; i++) g.next(A[i]);
samer@0: 	}
samer@0: 
samer@0: 	public void apply( Function f) {
samer@0: 		for (int i=0; i<m; i++) f.apply(A[i]);
samer@0: 	}
samer@0: 
samer@0: 	public void add( Generator g)
samer@0: 	{
samer@0: 		for (int i=0; i<m; i++) {
samer@0: 			for (int j=0; j<n; j++) {
samer@0: 				A[i][j]+=g.next();
samer@0: 			}
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 
samer@0:    /** Set to identity matrix */
samer@0: 
samer@0:    public void identity()
samer@0:    {
samer@0:       for (int i = 0; i < m; i++) {
samer@0:          for (int j = 0; j < n; j++) {
samer@0:             A[i][j] = (i == j ? 1.0 : 0.0);
samer@0:          }
samer@0:       }
samer@0:    }
samer@0: 
samer@0: 	public static Jama.Matrix identity(Jama.Matrix a)
samer@0: 	{
samer@0: 		double [][] A=a.getArray();
samer@0: 		int m=a.getRowDimension();
samer@0: 		int n=a.getColumnDimension();
samer@0: 		for (int i = 0; i < m; i++) {
samer@0: 			for (int j = 0; j < n; j++) {
samer@0: 				A[i][j] = (i == j ? 1.0 : 0.0);
samer@0: 			}
samer@0: 		}
samer@0: 		return a;
samer@0: 	}
samer@0: 	
samer@0:    /** Set to zero matrix */
samer@0: 
samer@0:    public void zero() {
samer@0:       for (int i = 0; i < m; i++) { Mathx.zero(A[i]); }
samer@0:    }
samer@0: 
samer@0:    /** Set to given matrix
samer@0:      * (need to make sure both matrices are the same size)
samer@0: 		@param b the matrix to be copied
samer@0:      */
samer@0: 
samer@0:    public void assign( Jama.Matrix b)
samer@0:    {
samer@0:       double[][] B=b.getArray();
samer@0: 	  int mmax = B.length>m ? m : B.length;
samer@0: 	  int nmax = B[0].length>n ? n : B[0].length;
samer@0: 
samer@0:       for (int i = 0; i < mmax; i++) {
samer@0: 		System.arraycopy(B[i],0,A[i],0,nmax);
samer@0:       }
samer@0:    }
samer@0: 
samer@0: 	// Saveable implementation
samer@0: 
samer@0: 	public void load(InputStream in) throws Exception {
samer@0: 		ObjectInputStream	o=new ObjectInputStream(in);
samer@0: 		assign((Jama.Matrix)o.readObject());
samer@0: 		changed();
samer@0: 	}
samer@0: 
samer@0: 	public void save(OutputStream out) throws Exception {
samer@0: 		ObjectOutputStream	o=new ObjectOutputStream(out);
samer@0: 		o.writeObject(new Jama.Matrix(getArray()));
samer@0: 	}
samer@0: 
samer@0: 	public void write(Writer wr) throws Exception {
samer@0: 		print(new PrintWriter(wr,true), new DoubleFormat(16),4);
samer@0: 	}
samer@0: 
samer@0: 	public void read(Reader rdr) throws Exception {
samer@0: 		Jama.Matrix A=Jama.Matrix.read(new BufferedReader(rdr));
samer@0: 		assign(A);	changed();
samer@0: 	}
samer@0: 
samer@0: /** Symmetrise in place preserving trace.
samer@0: 		Will fail if matrix is not square. */
samer@0: 	public void symmetrise()
samer@0: 	{
samer@0: 		double a, _G[][]=getArray();
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			for (int j=i+1; j<n; j++) {
samer@0: 				a=(_G[i][j]+_G[j][i])/2;
samer@0: 				_G[i][j]=_G[j][i]=a;
samer@0: 			}
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	/** Antisymmetrise in place.
samer@0: 		Will fail if matrix is not square. */
samer@0: 	public void antisymmetrise()
samer@0: 	{
samer@0: 		double a, _G[][]=getArray();
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			for (int j=i+1; j<n; j++) {
samer@0: 				a=(_G[i][j]-_G[j][i])/2;
samer@0: 				_G[i][j]=a;	_G[j][i]=-a;
samer@0: 			}
samer@0: 			_G[i][i]=0;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	/** return log(abs(det(this))) (must be square matrix) */
samer@0: 	
samer@0: 	public double logdet() {
samer@0: 		int n=getRowDimension();
samer@0: 		double [][] x=lu().getU().getArray();
samer@0: 		double L=0;
samer@0: 		for (int i=0; i<n; i++) L+= Math.log(Math.abs(x[i][i]));
samer@0: 		return L;
samer@0: 	}
samer@0: 	// ................................................
samer@0: 
samer@0: 	public Vec getRow(int k) { return new Row(k); }
samer@0: 	public Vec getColumn(int k) { return new Column(k); }
samer@0: 
samer@0: 	public static Vec getDiagonal(final Jama.Matrix A) {
samer@0: 		return new Vec() {
samer@0: 			double[][] x=A.getArray();
samer@0: 
samer@0: 			public int size() { return x.length; }
samer@0: 			public double[] array() { return null; }
samer@0: 			public Mat mat() {
samer@0: 				return new Mat() {
samer@0: 					// Mat implementation
samer@0: 					public int width() { return x.length; }
samer@0: 					public int height() { return 1; }
samer@0: 					public double get( int i, int j) { return x[j][j]; }
samer@0: 					public void set( int i, int j, double t) { x[j][j]=t; }
samer@0: 				};
samer@0: 			}
samer@0: 			public Iterator iterator() {
samer@0: 				return new Iterator() {
samer@0: 					int i=0, n=x.length;
samer@0: 					public void rewind() { i=0; }
samer@0: 					public boolean more() { return i<n; }
samer@0: 					public double next() {
samer@0: 						double y=x[i][i]; i++; return y;
samer@0: 					}
samer@0: 				};
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	private abstract class Slice implements Vec, Vec.Iterator
samer@0: 	{
samer@0: 		protected int i, j, n;
samer@0: 
samer@0: 		Slice(int n, int j) { this.n=n; this.j=j; i=0; }
samer@0: 
samer@0: 		// Vec implementation
samer@0: 		public final Vec.Iterator	iterator() { i=0; return this; }
samer@0: 		public final int	size() { return n; }
samer@0: 		public final void	rewind() { i=0; }
samer@0: 		public final boolean more() { return i<n; }
samer@0: 	}
samer@0: 
samer@0: 	private class Row extends Slice {
samer@0: 		double []	x;
samer@0: 
samer@0: 		public Row(int rownum) {
samer@0: 			super(Matrix.this.getColumnDimension(),rownum);
samer@0: 			x=getArray()[rownum];
samer@0: 		}
samer@0: 
samer@0: 		public final double		next() { return x[i++]; }
samer@0: 		public final double[]	array() { return x; }
samer@0: 		public final Mat		mat() {
samer@0: 			return new Mat() {
samer@0: 				// Mat implementation
samer@0: 				public int width() { return Row.this.n; }
samer@0: 				public int height() { return 1; }
samer@0: 				public double get( int i, int j) { return x[j]; }
samer@0: 				public void set( int i, int j, double t) { x[j]=t; }
samer@0: 			};
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	private class Column extends Slice {
samer@0: 		double [][]	x;
samer@0: 		int			j;
samer@0: 
samer@0: 		public Column(int colnum) {
samer@0: 			super(Matrix.this.getRowDimension(),colnum);
samer@0: 			x=getArray();
samer@0: 			j=colnum;
samer@0: 		}
samer@0: 
samer@0: 		public final double next() { return x[i++][j]; }
samer@0: 		public final double[] array() { return null; }
samer@0: 		public final Mat		mat() {
samer@0: 			return new Mat() {
samer@0: 				public final int width() { return Column.this.n; }
samer@0: 				public final int height() { return 1; }
samer@0: 				public final double get( int i, int k) { return x[k][j]; }
samer@0: 				public final void set( int i, int k, double t) { x[k][j]=t; }
samer@0: 			};
samer@0: 		}
samer@0: 	}
samer@0: }