samer@0: /*
samer@0:  *	Copyright (c) 2002, 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: 
samer@0: import samer.core.*;
samer@0: import samer.core.types.*;
samer@0: import samer.maths.*;
samer@0: import samer.maths.opt.*;
samer@0: import samer.tools.*;
samer@0: 
samer@0: /**
samer@0: 	Differential scaler: scales and offsets each element of a vector
samer@0: 	independently aiming to match a given prior model. This is
samer@0: 	like ICA using a diagonal weight matrix, with a 'zero-mean'
samer@0: 	normalisation built in, though it doesn't actually the mean to
samer@0: 	do this, but some statistic based on the prior model.
samer@0: */
samer@0: 
samer@0: public class DiffScaler extends AnonymousTask implements Model
samer@0: {
samer@0: 	private Model 		M; // models P(s)
samer@0: 	private int			n;
samer@0: 	private Vec			x;
samer@0: 	private VVector	s;
samer@0: 	private VVector	w;		// vector of multipliers
samer@0: 	private VVector	mu;	// vector of offsets
samer@0: 	private VDouble	logA;
samer@0: 
samer@0: 	double []		_x, _s, _g, _w, _m, phi;
samer@0: 
samer@0: 	public DiffScaler( Vec input, Model M) { this(input); setOutputModel(M); M.setInput(s); }
samer@0: 	public DiffScaler( Vec input) { this(input.size()); setInput(input); }
samer@0: 	public DiffScaler( int N)
samer@0: 	{
samer@0: 		n = N;
samer@0: 
samer@0: 		x = null;
samer@0: 		mu = new VVector("mu",n);	mu.addSaver();
samer@0: 		w = new VVector("w",n);		w.addSaver();
samer@0: 		s = new VVector("output",n);
samer@0: 		logA = new VDouble("log|A|");
samer@0: 
samer@0: 		_s = s.array();
samer@0: 		_w = w.array();
samer@0: 		_m = mu.array();
samer@0: 		_g = new double[n];
samer@0: 		phi = null;
samer@0: 		reset();
samer@0: 	}
samer@0: 
samer@0: 	public int getSize() { return n; }
samer@0: 	public VVector output() { return s; }
samer@0: 	public VVector weights() { return w; }
samer@0: 	public VVector offsets() { return mu; }
samer@0: 	public Model getOutputModel() { return M;	}
samer@0: 	public void setOutputModel(Model m) { M=m;	}
samer@0: 	public void setInput(Vec in) { x=in; _x=x.array(); }
samer@0: 	public void reset() { reset(1.0); }
samer@0: 	public void reset(double k) {
samer@0: 		Mathx.setAll(_w,k); w.changed();
samer@0: 		Mathx.setAll(_m,0); mu.changed();
samer@0: 		logA.set(-sumlog(_w));
samer@0: 	}
samer@0: 
samer@0: 	public String toString() { return "DiffScaler("+x+")"; }
samer@0: 	private static double sumlog(double [] x) {
samer@0: 		double S=0;
samer@0: 		for (int i=0; i<x.length; i++) S += Math.log(x[i]);
samer@0: 		return S;
samer@0: 	}
samer@0: 
samer@0: 	public void dispose()
samer@0: 	{
samer@0: 		logA.dispose();
samer@0: 		mu.dispose();
samer@0: 		w.dispose();
samer@0: 		s.dispose();
samer@0: 		super.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public void infer() {
samer@0: 		for (int i=0; i<n; i++) _s[i] = _w[i]*(_x[i]-_m[i]);
samer@0: 		s.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void compute() {
samer@0: 		Mathx.mul(_g,M.getGradient(),_w);
samer@0: 	}
samer@0: 
samer@0: 	public double	getEnergy() { return M.getEnergy() + logA.value; }
samer@0: 	public double [] getGradient() { return _g; }
samer@0: 
samer@0: 	public Functionx functionx() {
samer@0: 		return new Functionx() {
samer@0: 			Functionx fM=M.functionx();
samer@0: 			double [] s=new double[n];
samer@0: 
samer@0: 			public void dispose() { fM.dispose(); }
samer@0: 			public void evaluate(Datum P) { P.f=evaluate(P.x,P.g); }
samer@0: 			public double evaluate(double [] x, double [] g) {
samer@0: 				for (int i=0; i<n; i++) s[i] = _w[i]*(x[i]-_m[i]);
samer@0: 				double E=fM.evaluate(s,g);
samer@0: 				Mathx.mul(g,_w);
samer@0: 				return E+logA.value;
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	public void starting() { logA.set(-sumlog(_w)); }
samer@0: 	public void stopping() {}
samer@0: 	public void run() { infer(); }
samer@0: 
samer@0: 	public Model.Trainer getTrainer() { return new Trainer(); }
samer@0: 	public Model.Trainer getOffsetTrainer() { return new OffsetTrainer(); }
samer@0: 	public Model.Trainer getScaleTrainer() { return new ScaleTrainer(); }
samer@0: 	public Model.Trainer getTensionedTrainer() { return  new TensionedTrainer(); }
samer@0: 
samer@0: 
samer@0: 	public class Trainer extends AnonymousTask implements Model.Trainer
samer@0: 	{
samer@0: 		VDouble		rate1=new VDouble("scaleRate",0.0001);
samer@0: 		VDouble		rate2=new VDouble("offsetRate",0.000001);
samer@0: 		double[]		G,H;
samer@0: 		double []	_s;
samer@0: 		double		count; // n
samer@0: 
samer@0: 
samer@0: 		public Trainer() {
samer@0: //			n=DiffScaler.this.n;
samer@0: 			_s = s.array();
samer@0: 			G = new double[n];
samer@0: 			H = new double[n];
samer@0: 		}
samer@0: 
samer@0: 		public String toString() { return "Trainer:"+DiffScaler.this; }
samer@0: 
samer@0: 		public VDouble getScaleRate() { return rate1; }
samer@0: 		public VDouble getOffsetRate() { return rate2; }
samer@0: 
samer@0: 		public void reset() { count=0; Mathx.zero(G); Mathx.zero(H); }
samer@0: 		public void accumulate() { accumulate(1); }
samer@0: 		public void accumulate(double w) {
samer@0: 			double [] phi=M.getGradient();
samer@0: 			for (int i=0; i<n; i++) {
samer@0: 				G[i] += w*(phi[i]*_s[i] - 1);
samer@0: 				H[i] += w*phi[i];
samer@0: 			}
samer@0: 			count+=w;
samer@0: 		}
samer@0: 
samer@0: 		public void flush() {
samer@0: 			if (count==0) return; // nothing to do
samer@0: 
samer@0: 			double eta1 = rate1.value/count;
samer@0: 			double eta2 = rate2.value/count;
samer@0: 
samer@0: 			for (int i=0; i<n; i++) {
samer@0: 				_m[i] += eta2*H[i]/_w[i];
samer@0: 				_w[i] *= Math.exp(-eta1*G[i]);
samer@0: 			}
samer@0: 			logA.value+=eta1*Mathx.sum(G);
samer@0: 			logA.changed();
samer@0: 			mu.changed();
samer@0: 			w.changed();
samer@0: 			reset();
samer@0: 		}
samer@0: 
samer@0: 		public void oneshot() { reset(); accumulate(); flush(); }
samer@0: 		public void dispose() { rate1.dispose(); rate2.dispose(); }
samer@0: 		public void starting() { reset(); }
samer@0: 		public void run() { accumulate(); flush(); }
samer@0: 	}
samer@0: 
samer@0: 	public class ScaleTrainer extends AnonymousTask implements Model.Trainer
samer@0: 	{
samer@0: 		VDouble	scale=new VDouble("scale",0.001);
samer@0: 		VDouble	stretch=new VDouble("stretch",0.001/n);
samer@0: 		double		thresh=Shell.getDouble("anomaly",20*n);
samer@0: 		double[]		G, H, _s;
samer@0: 		double		count; // n
samer@0: 
samer@0: 
samer@0: 		public ScaleTrainer() {
samer@0: 			_s = s.array();
samer@0: 			G = new double[n];
samer@0: 			// H = new double[n];
samer@0: 		}
samer@0: 
samer@0: 		public String toString() { return "ScaleTrainer:"+DiffScaler.this; }
samer@0: 
samer@0: 		public VDouble getScaleRate() { return scale; }
samer@0: 		public VDouble getStretchRate() { return stretch; }
samer@0: 
samer@0: 		public void reset() { count=0; Mathx.zero(G); }
samer@0: 		public void accumulate() { accumulate(1); }
samer@0: 		public void accumulate(double w) {
samer@0: 			if (M.getEnergy()>thresh) return;
samer@0: 			double [] phi=M.getGradient();
samer@0: 			for (int i=0; i<n; i++) {
samer@0: 				G[i] += w*(phi[i]*_s[i] - 1);
samer@0: 				// H[i] += w*phi[i];
samer@0: 			}
samer@0: 			count+=w;
samer@0: 		}
samer@0: 
samer@0: 		public void flush() {
samer@0: 			if (count==0) return; // nothing to do
samer@0: 
samer@0: 
samer@0: 			{	// filter elements of G
samer@0: 				double beta = stretch.value/scale.value;
samer@0: 				double meanG = Mathx.sum(G)/n;
samer@0: 				for (int i=0; i<n; i++) {
samer@0: 					G[i]=meanG+beta*(G[i]-meanG);
samer@0: 				}
samer@0: 			}
samer@0: 
samer@0: 			double alpha = scale.value/count;
samer@0: 			for (int i=0; i<n; i++) {
samer@0: 				double tmp=Math.exp(-alpha*G[i]);
samer@0: 				if (Double.isNaN(tmp)) throw new Error("alt: NaN"+i);
samer@0: 				_w[i] *= tmp;
samer@0: 			}
samer@0: 			logA.value+=alpha*Mathx.sum(G);
samer@0: 			logA.changed();
samer@0: 			w.changed();
samer@0: 			reset();
samer@0: 		}
samer@0: 
samer@0: 		public void oneshot() { reset(); accumulate(); flush(); }
samer@0: 		public void dispose() { scale.dispose(); stretch.dispose(); }
samer@0: 		public void starting() { reset(); }
samer@0: 		public void run() { accumulate(); flush(); }
samer@0: 	}
samer@0: 
samer@0: 	public class OffsetTrainer extends AnonymousTask implements Model.Trainer
samer@0: 	{
samer@0: 		VDouble		rate2=new VDouble("offsetRate",0.000001);
samer@0: 		double[]		H;
samer@0: 		double		count; // n
samer@0: 
samer@0: 
samer@0: 		public OffsetTrainer() {
samer@0: 	//		n=DiffScaler.this.n;
samer@0: 			H = new double[n];
samer@0: 		}
samer@0: 
samer@0: 		public String toString() { return "OffsetTrainer:"+DiffScaler.this; }
samer@0: 
samer@0: 		public VDouble getOffsetRate() { return rate2; }
samer@0: 
samer@0: 		public void reset() { count=0; Mathx.zero(H); }
samer@0: 		public void accumulate() { accumulate(1); }
samer@0: 		public void accumulate(double w) {
samer@0: 			double [] phi=M.getGradient();
samer@0: 			for (int i=0; i<n; i++) H[i] += w*phi[i];
samer@0: 			count+=w;
samer@0: 		}
samer@0: 
samer@0: 		public void flush() {
samer@0: 			if (count==0) return; // nothing to do
samer@0: 
samer@0: 			double eta2 = rate2.value/count;
samer@0: 
samer@0: 			for (int i=0; i<n; i++) {
samer@0: 				_m[i] += eta2*H[i]/_w[i];
samer@0: 			}
samer@0: 			mu.changed();
samer@0: 			reset();
samer@0: 		}
samer@0: 
samer@0: 		public void oneshot() { reset(); accumulate(1); flush(); }
samer@0: 		public void dispose() { rate2.dispose(); }
samer@0: 		public void starting() { reset(); }
samer@0: 		public void run() { accumulate(1); flush(); }
samer@0: 	}
samer@0: 	
samer@0: 	public class TensionedTrainer extends Trainer
samer@0: 	{
samer@0: 		VDouble		tension=new VDouble("tension",0.01);
samer@0: 		double			lw[]=new double[n];
samer@0: 
samer@0: 		public TensionedTrainer() {}
samer@0: 
samer@0: 		public String toString() { return "TensionedTrainer:"+DiffScaler.this; }
samer@0: 
samer@0: 		public VDouble getTension() { return tension; }
samer@0: 
samer@0: 		public void flush() {
samer@0: 			double T=tension.value;
samer@0: 			int i;
samer@0: 
samer@0: 			// go through and modify G
samer@0: 			for (i=0; i<n; i++) { lw[i]=Math.log(_w[i]); }
samer@0: 			G[0] -= T*(lw[1]-lw[0]);	i=n-1;
samer@0: 			G[i] -= T*(lw[i-1]-lw[i]); 	i--;
samer@0: 			for (; i>0; i--) {
samer@0: 				G[i] -= T*(lw[i-1] -2*lw[i]+lw[i+1]);
samer@0: 			}
samer@0: 
samer@0: 			// do the same for H?
samer@0: 			// may have to modify T.
samer@0: 			H[0] -= T*(_m[1]-_m[0]);	i=n-1;
samer@0: 			H[i] -= T*(_m[i-1]-_m[i]); 	i--;
samer@0: 			for (; i>0; i--) {
samer@0: 				H[i] -= T*(_m[i-1] -2*_m[i]+_m[i+1]);
samer@0: 			}
samer@0: 
samer@0: 			super.flush();
samer@0: 		}
samer@0: 
samer@0: 		public void dispose() { tension.dispose(); super.dispose(); }
samer@0: 	}
samer@0: }
samer@0: