samer@0: /*
samer@0:  *	Histogram.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.units;
samer@0: 
samer@0: import samer.core.*;
samer@0: import samer.core.types.*;
samer@0: import samer.core.util.*;
samer@0: import samer.core.Agent.*;
samer@0: import samer.tools.*;
samer@0: import samer.maths.*;
samer@0: import java.util.*;
samer@0: 
samer@0: /**
samer@0: 	This actually builds several histograms in parallel: one
samer@0: 	for each element of a Vec or double array.  A IMap provides
samer@0: 	the mapping from element values to histogram bins.
samer@0: 	The result is a <i>Matrix</i>, each <i>row</i> of which
samer@0: 	is the histogram for the corresponding element of the
samer@0: 	input vector.
samer@0:   */
samer@0: 
samer@0: public class Histogram extends Viewable implements Task, Agent
samer@0: {
samer@0: 	private Vec			input;
samer@0: 	private Matrix			bins;
samer@0: 	private VMap			map;
samer@0: 	private double[][]	binArray;
samer@0: 	private double[]		lbw; // log of bin widths
samer@0: 	private VVector		L;
samer@0: 	private VDouble		sumL;
samer@0: 	private int				N, count;
samer@0: 	private double[]		x;
samer@0: 
samer@0: 	public Histogram( Vec input, int binCount)
samer@0: 	{
samer@0: 		super("histogram");
samer@0: 
samer@0: 		Shell.push(getNode());
samer@0: 
samer@0: 		N = input.size();
samer@0: 		Shell.setAutoRegister(false);
samer@0: 		bins = new Matrix("bins",N,binCount);
samer@0: 		map = new VMap(new LinearMap(binCount));
samer@0: 		L = new VVector("likelihood",N);
samer@0: 		sumL = new VDouble("sumL");
samer@0: 		Shell.setAutoRegister(true);
samer@0: 		count=0; x=input.array();
samer@0: 		lbw = new double[binCount];
samer@0: 		Shell.pop();
samer@0: 
samer@0: 		setAgent(this);
samer@0: 		binArray = bins.getArray();
samer@0: 		Shell.registerViewable(this);
samer@0: 	}
samer@0: 
samer@0: 	public Matrix getBinMatrix() { return bins; }
samer@0: 	public VVector getLikelihoodVector() { return L; }
samer@0: 	public VDouble getLikelihoodSignal() { return sumL; }
samer@0: 
samer@0: 	public void starting() {}
samer@0: 	public void stopping() {}
samer@0: 	public void run()
samer@0: 	{
samer@0: 		IMap	binmap=map.getMap();
samer@0: 		double [] l=L.array();
samer@0: 
samer@0: 		if (x!=null) {
samer@0: 			for(int k=0; k<N; k++) {
samer@0: 				int j=binmap.clipInt(x[k]);
samer@0: 				l[k] = -Math.log(++binArray[k][j]) + lbw[j];
samer@0: 			}
samer@0: 		} else {
samer@0: 			Vec.Iterator i=input.iterator();
samer@0: 			for(int k=0; i.more(); k++) {
samer@0: 				int j=binmap.clipInt(i.next());
samer@0: 				l[k] = -Math.log(++binArray[k][j]) + lbw[j];
samer@0: 			}
samer@0: 		}
samer@0: 
samer@0: 		count++;
samer@0: 		Mathx.add(l,Math.log(count));
samer@0: 		sumL.value=Mathx.sum(l);
samer@0: 
samer@0: 		L.changed();
samer@0: 		sumL.changed();
samer@0: 		bins.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void clear() { bins.zero(); bins.changed(); count=0; }
samer@0: 	public double[] normalise() {
samer@0: 		count=(int)Mathx.sum(binArray[0]);
samer@0: 		Shell.print("Histogram: count="+count);
samer@0: 		Shell.print("Computing bin widths...");
samer@0: 		IMap	binmap=map.getMap();
samer@0: 
samer@0: 		for (int i=0; i<lbw.length; i++) {
samer@0: 			lbw[i]=Math.log(binmap.inverseFromInt(i+1) - binmap.inverseFromInt(i));
samer@0: 		}
samer@0: 		return lbw;
samer@0: 	}
samer@0: 
samer@0: 	public void dispose() {
samer@0: 		bins.dispose();
samer@0: 		map.dispose();
samer@0: 		L.dispose();
samer@0: 		super.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public Viewer getViewer() {
samer@0: 		DefaultViewer vwr=new DefaultViewer(this);
samer@0: 		vwr.add(bins.viewable());
samer@0: 		vwr.add(L);
samer@0: 		vwr.add(sumL);
samer@0: 		return vwr;
samer@0: 	}
samer@0: 
samer@0: 	public void getCommands(Registry r) {
samer@0: 		r.add("clear").add("normalise");
samer@0: 		r.group(); map.getCommands(r);
samer@0: 	}
samer@0: 
samer@0: 	public void execute(String cmd, Environment env) throws Exception {
samer@0: 		if (cmd.equals("clear")) clear();
samer@0: 		else if (cmd.equals("normalise")) normalise();
samer@0: 		map.execute(cmd,env);
samer@0: 	}
samer@0: 
samer@0: 	public Equaliser getEqualiser() { return new Equaliser(); }
samer@0: 
samer@0: 	/** This is a function which uses a cumulative probability function
samer@0: 		derived from the histogram to compute a transformation that
samer@0: 		results in a uniformly distributed variable. It takes a SNAPSHOT
samer@0: 		of the current histogram map and bin counts. You must update()
samer@0: 		to update from current histogram.
samer@0: 		*/
samer@0: 	public class Equaliser extends VectorFunctionOfVector
samer@0: 	{
samer@0: 		IMap	map;
samer@0: 		int 	N, M;
samer@0: 		Matrix	C;
samer@0: 		double K;
samer@0: 		double [][]P;
samer@0: 
samer@0: 		public Equaliser() {
samer@0: 			N=bins.getRowDimension();
samer@0: 			M=bins.getColumnDimension();
samer@0: 			C=new Matrix("cumulative",N,M);
samer@0: 			P=C.getArray();
samer@0: 			update();
samer@0: 		}
samer@0: 
samer@0: 		public void dispose() { C.dispose(); }
samer@0: 		public void update() {
samer@0: 			map=Histogram.this.map.getMap();
samer@0: 			K=map.getIntRange();
samer@0: 			for (int i=0; i<N; i++) {
samer@0: 				double T=0;
samer@0: 				for (int j=0; j<M; j++) {
samer@0: 					P[i][j]=(T+=binArray[i][j]);
samer@0: 				}
samer@0: 				Mathx.mul(P[i],1/T);
samer@0: 			}
samer@0: 			C.changed();
samer@0: 		}
samer@0: 
samer@0: 		public void apply(double [] x) { apply(x,x); }
samer@0: 		public void apply(double[] x, double[] y) {
samer@0: 			for (int i=0; i<N; i++) {
samer@0: 				double z=K*map.map(x[i]);
samer@0: 				int j=(int)z;
samer@0: 				if (j<0) y[i]=0;
samer@0: 				else if (j>=M) y[i]=1;
samer@0: 				else {
samer@0: 					// linearly interpolate 
samer@0: 					double P1=P[i][j];
samer@0: 					double P0=(j==0)? 0 : P[i][j-1];
samer@0: 					double d=z-j;
samer@0: 					y[i]=(1-d)*P0 + d*P1;
samer@0: 				}
samer@0: 			}
samer@0: 		}
samer@0: 	}
samer@0: }
samer@0: