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.models;
samer@0: import  samer.core.*;
samer@0: import  samer.core.types.*;
samer@0: import  samer.tools.*;
samer@0: import  samer.maths.*;
samer@0: import  samer.maths.*;
samer@0: import  samer.maths.opt.*;
samer@0: 
samer@0: public class NoisyICA extends NamedTask implements Model
samer@0: {
samer@0: 	Vec				x;				// data (input)
samer@0: 	int				n, m;			// sizes (data, sources)
samer@0: 	Matrix		A;				// basis matrix
samer@0: 	Model		Ms, Me;	// source and noise models
samer@0: 	VVector	s, z, e;		// sources, reconstruction, error
samer@0: 	VDouble		E;				// energy
samer@0: 
samer@0: 	Task			inference=new NullTask();
samer@0: 
samer@0: 
samer@0: 	// ----- variables used in computations -----------------
samer@0: 
samer@0: 	private	double []	e_;
samer@0: 	private	double []	x_;
samer@0: 	private	double []	z_;
samer@0: 	private	double []	s_;
samer@0: 	private	VectorFunctionOfVector tA, tAt;
samer@0: 
samer@0: 
samer@0: 	public NoisyICA(Vec in,int outs) { this(new Node("noisyica"),in.size(),outs); setInput(in); }
samer@0: 	public NoisyICA(int ins,int outs) { this(new Node("noisyica"),ins,outs); }
samer@0: 	public NoisyICA(Node node, int inputs, int outputs)
samer@0: 	{
samer@0: 		super(node);
samer@0: 		Shell.push(node);
samer@0: 
samer@0: 		n = inputs;
samer@0: 		m = outputs;
samer@0: 
samer@0: 		s = new VVector("s",m);
samer@0: 		z = new VVector("z",n);
samer@0: 		e = new VVector("e",n);
samer@0: 		E = new VDouble("E");
samer@0: 		A = new Matrix("A",n,m);
samer@0: 		A.identity();
samer@0: 
samer@0: 		e_ = e.array();
samer@0: 		z_ = z.array();
samer@0: 		s_ = s.array();
samer@0: 		tA = new MatrixTimesVector(A);
samer@0: 		tAt= new MatrixTransposeTimesVector(A);
samer@0: 
samer@0: 		Shell.pop();
samer@0: 	}
samer@0: 
samer@0: 	public Model getSourceModel() { return Ms; }
samer@0: 	public Model getNoiseModel() { return Me; }
samer@0: 	public void  setSourceModel(Model m) { Ms=m; Ms.setInput(s); }
samer@0: 	public void  setNoiseModel(Model m) { Me=m; Me.setInput(e); }
samer@0: 	public void  setInput(Vec in) { x=in; x_ = x.array(); }
samer@0: 	public Matrix basisMatrix() { return A; }
samer@0: 	public VVector output() { return s; }
samer@0: 	public VVector error() { return e; }
samer@0: 	public VVector reconstruction() { return z; }
samer@0: 
samer@0: 	public int  getSize() { return n; }
samer@0: 
samer@0: 	public void setInferenceTask(Task t) { inference=t; }
samer@0: 
samer@0: 	public void infer() {
samer@0: 		try {	inference.run(); }
samer@0: 		catch (Exception ex) { 
samer@0: 			Shell.trace("error: "+ex);
samer@0: 			ex.printStackTrace();
samer@0: 			throw new Error("inference failed: "+ex); }
samer@0: 		tA.apply(s_,z_); Mathx.sub(e_,x_,z_);
samer@0: 		e.changed(); z.changed(); s.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void compute() {
samer@0: 		E.set(Me.getEnergy() + Ms.getEnergy());
samer@0: 		// what about dE/dx?
samer@0: 	}
samer@0: 	public double getEnergy() { return E.value; }
samer@0: 	public double [] getGradient() { return null; } // this is wrong
samer@0: 
samer@0: 	/** get basis vector norms into given array */
samer@0: 	public void norms(double [] na)
samer@0: 	{
samer@0: 		double [][] M=A.getArray();
samer@0: 
samer@0: 		Mathx.zero(na);
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double [] Mi=M[i];
samer@0: 			for (int j=0; j<m; j++) {
samer@0: 				na[j] += Mi[j]*Mi[j];
samer@0: 			}
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	public void run() { infer(); }
samer@0: 	public void dispose()
samer@0: 	{
samer@0: 		s.dispose();
samer@0: 		A.dispose();
samer@0: 		e.dispose();
samer@0: 		z.dispose();
samer@0: 		E.dispose();
samer@0: 		tA.dispose();
samer@0: 		tAt.dispose();
samer@0: 
samer@0: 		super.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public Functionx functionx() { return null; }
samer@0: 
samer@0: 	public Functionx posterior()	{
samer@0: 		// returns Functionx which evaluates E and dE/ds at current x
samer@0: 		return new Functionx() {
samer@0: 			Functionx fMs=Ms.functionx();
samer@0: 			Functionx fMe=Me.functionx();
samer@0: 			double [] e=new double[n];
samer@0: 			double [] ge=new double[n];
samer@0: 			double [] gs=new double[m];
samer@0: 
samer@0: 			public void dispose() { fMs.dispose(); fMe.dispose(); }
samer@0: 			public void evaluate( Datum P) { P.f=evaluate(P.x,P.g); }
samer@0: 			public double evaluate( double [] s, double [] g)
samer@0: 			{
samer@0: 				tA.apply(s,z_); Mathx.sub(e,x_,z_);
samer@0: 				double Ee=fMe.evaluate(e,ge); tAt.apply(ge,gs); // gs=A'*gamma(e)
samer@0: 				double Es=fMs.evaluate(s,g);	Mathx.sub(g,gs); // g=gamma(s)-ge
samer@0: 				return Es+Ee;
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	public Trainer learnHebbian() 	{ return new MatrixTrainer(e_,A,s_); }
samer@0: 
samer@0: 	public Trainer learnLewickiSejnowski()
samer@0: 	{
samer@0: 		final double [] h =new double[n];
samer@0: 		final double [] f =new double[m];
samer@0: 		return new MatrixTrainer(h,A,s_) {
samer@0: 
samer@0: 			public void accumulate(double w) {
samer@0: 				// tAt.apply(Me.getGradient(),f); tA.apply(f,h);
samer@0: 				tA.apply(Ms.getGradient(),h);
samer@0: 				super.accumulate(w);
samer@0: 			}
samer@0: 			public void flush() { // flush with decay
samer@0: 				for (int j=0; j<this.m; j++)
samer@0: 					for (int i=0; i<this.n; i++)
samer@0: 						_T[i][j] -= count*_A[i][j];
samer@0: 				super.flush();
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	public Trainer learnDecayWhenActive()
samer@0: 	{
samer@0: 		final double [] h =new double[n];
samer@0: 		return new MatrixTrainer(h,A,s_) {
samer@0: 			VDouble				th=new VDouble("threshold",0.01);
samer@0: 			public void accumulate(double w) {
samer@0: 				// perhaps would like to get info out of optimiser here
samer@0: 				tA.apply(Ms.getGradient(),h);
samer@0: 				super.accumulate(w);
samer@0: 
samer@0: 				// decay when active part
samer@0: 				double thresh=th.value;
samer@0: 				for (int j=0; j<this.m; j++)
samer@0: 					if (isActive(s_[j],thresh))
samer@0: 						for (int i=0; i<this.n; i++)
samer@0: 							_T[i][j] -= _A[i][j];
samer@0: 			}
samer@0: 
samer@0: 			public VDouble getThreshold() { return th; }
samer@0: 		};
samer@0: 	}
samer@0: 	static boolean isActive(double s, double t) { return s>t || s<-t; }
samer@0: }