annotate src/samer/models/DiffScaler.java @ 0:bf79fb79ee13

Initial Mercurial check in.

/*
 *	Copyright (c) 2002, 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.models;

import samer.core.*;
import samer.core.types.*;
import samer.maths.*;
import samer.maths.opt.*;
import samer.tools.*;

/**
	Differential scaler: scales and offsets each element of a vector
	independently aiming to match a given prior model. This is
	like ICA using a diagonal weight matrix, with a 'zero-mean'
	normalisation built in, though it doesn't actually the mean to
	do this, but some statistic based on the prior model.
*/

public class DiffScaler extends AnonymousTask implements Model
{
	private Model 		M; // models P(s)
	private int			n;
	private Vec			x;
	private VVector	s;
	private VVector	w;		// vector of multipliers
	private VVector	mu;	// vector of offsets
	private VDouble	logA;

	double []		_x, _s, _g, _w, _m, phi;

	public DiffScaler( Vec input, Model M) { this(input); setOutputModel(M); M.setInput(s); }
	public DiffScaler( Vec input) { this(input.size()); setInput(input); }
	public DiffScaler( int N)
	{
		n = N;

		x = null;
		mu = new VVector("mu",n);	mu.addSaver();
		w = new VVector("w",n);		w.addSaver();
		s = new VVector("output",n);
		logA = new VDouble("log|A|");

		_s = s.array();
		_w = w.array();
		_m = mu.array();
		_g = new double[n];
		phi = null;
		reset();
	}

	public int getSize() { return n; }
	public VVector output() { return s; }
	public VVector weights() { return w; }
	public VVector offsets() { return mu; }
	public Model getOutputModel() { return M;	}
	public void setOutputModel(Model m) { M=m;	}
	public void setInput(Vec in) { x=in; _x=x.array(); }
	public void reset() { reset(1.0); }
	public void reset(double k) {
		Mathx.setAll(_w,k); w.changed();
		Mathx.setAll(_m,0); mu.changed();
		logA.set(-sumlog(_w));
	}

	public String toString() { return "DiffScaler("+x+")"; }
	private static double sumlog(double [] x) {
		double S=0;
		for (int i=0; i<x.length; i++) S += Math.log(x[i]);
		return S;
	}

	public void dispose()
	{
		logA.dispose();
		mu.dispose();
		w.dispose();
		s.dispose();
		super.dispose();
	}

	public void infer() {
		for (int i=0; i<n; i++) _s[i] = _w[i]*(_x[i]-_m[i]);
		s.changed();
	}

	public void compute() {
		Mathx.mul(_g,M.getGradient(),_w);
	}

	public double	getEnergy() { return M.getEnergy() + logA.value; }
	public double [] getGradient() { return _g; }

	public Functionx functionx() {
		return new Functionx() {
			Functionx fM=M.functionx();
			double [] s=new double[n];

			public void dispose() { fM.dispose(); }
			public void evaluate(Datum P) { P.f=evaluate(P.x,P.g); }
			public double evaluate(double [] x, double [] g) {
				for (int i=0; i<n; i++) s[i] = _w[i]*(x[i]-_m[i]);
				double E=fM.evaluate(s,g);
				Mathx.mul(g,_w);
				return E+logA.value;
			}
		};
	}

	public void starting() { logA.set(-sumlog(_w)); }
	public void stopping() {}
	public void run() { infer(); }

	public Model.Trainer getTrainer() { return new Trainer(); }
	public Model.Trainer getOffsetTrainer() { return new OffsetTrainer(); }
	public Model.Trainer getScaleTrainer() { return new ScaleTrainer(); }
	public Model.Trainer getTensionedTrainer() { return  new TensionedTrainer(); }


	public class Trainer extends AnonymousTask implements Model.Trainer
	{
		VDouble		rate1=new VDouble("scaleRate",0.0001);
		VDouble		rate2=new VDouble("offsetRate",0.000001);
		double[]		G,H;
		double []	_s;
		double		count; // n


		public Trainer() {
//			n=DiffScaler.this.n;
			_s = s.array();
			G = new double[n];
			H = new double[n];
		}

		public String toString() { return "Trainer:"+DiffScaler.this; }

		public VDouble getScaleRate() { return rate1; }
		public VDouble getOffsetRate() { return rate2; }

		public void reset() { count=0; Mathx.zero(G); Mathx.zero(H); }
		public void accumulate() { accumulate(1); }
		public void accumulate(double w) {
			double [] phi=M.getGradient();
			for (int i=0; i<n; i++) {
				G[i] += w*(phi[i]*_s[i] - 1);
				H[i] += w*phi[i];
			}
			count+=w;
		}

		public void flush() {
			if (count==0) return; // nothing to do

			double eta1 = rate1.value/count;
			double eta2 = rate2.value/count;

			for (int i=0; i<n; i++) {
				_m[i] += eta2*H[i]/_w[i];
				_w[i] *= Math.exp(-eta1*G[i]);
			}
			logA.value+=eta1*Mathx.sum(G);
			logA.changed();
			mu.changed();
			w.changed();
			reset();
		}

		public void oneshot() { reset(); accumulate(); flush(); }
		public void dispose() { rate1.dispose(); rate2.dispose(); }
		public void starting() { reset(); }
		public void run() { accumulate(); flush(); }
	}

	public class ScaleTrainer extends AnonymousTask implements Model.Trainer
	{
		VDouble	scale=new VDouble("scale",0.001);
		VDouble	stretch=new VDouble("stretch",0.001/n);
		double		thresh=Shell.getDouble("anomaly",20*n);
		double[]		G, H, _s;
		double		count; // n


		public ScaleTrainer() {
			_s = s.array();
			G = new double[n];
			// H = new double[n];
		}

		public String toString() { return "ScaleTrainer:"+DiffScaler.this; }

		public VDouble getScaleRate() { return scale; }
		public VDouble getStretchRate() { return stretch; }

		public void reset() { count=0; Mathx.zero(G); }
		public void accumulate() { accumulate(1); }
		public void accumulate(double w) {
			if (M.getEnergy()>thresh) return;
			double [] phi=M.getGradient();
			for (int i=0; i<n; i++) {
				G[i] += w*(phi[i]*_s[i] - 1);
				// H[i] += w*phi[i];
			}
			count+=w;
		}

		public void flush() {
			if (count==0) return; // nothing to do


			{	// filter elements of G
				double beta = stretch.value/scale.value;
				double meanG = Mathx.sum(G)/n;
				for (int i=0; i<n; i++) {
					G[i]=meanG+beta*(G[i]-meanG);
				}
			}

			double alpha = scale.value/count;
			for (int i=0; i<n; i++) {
				double tmp=Math.exp(-alpha*G[i]);
				if (Double.isNaN(tmp)) throw new Error("alt: NaN"+i);
				_w[i] *= tmp;
			}
			logA.value+=alpha*Mathx.sum(G);
			logA.changed();
			w.changed();
			reset();
		}

		public void oneshot() { reset(); accumulate(); flush(); }
		public void dispose() { scale.dispose(); stretch.dispose(); }
		public void starting() { reset(); }
		public void run() { accumulate(); flush(); }
	}

	public class OffsetTrainer extends AnonymousTask implements Model.Trainer
	{
		VDouble		rate2=new VDouble("offsetRate",0.000001);
		double[]		H;
		double		count; // n


		public OffsetTrainer() {
	//		n=DiffScaler.this.n;
			H = new double[n];
		}

		public String toString() { return "OffsetTrainer:"+DiffScaler.this; }

		public VDouble getOffsetRate() { return rate2; }

		public void reset() { count=0; Mathx.zero(H); }
		public void accumulate() { accumulate(1); }
		public void accumulate(double w) {
			double [] phi=M.getGradient();
			for (int i=0; i<n; i++) H[i] += w*phi[i];
			count+=w;
		}

		public void flush() {
			if (count==0) return; // nothing to do

			double eta2 = rate2.value/count;

			for (int i=0; i<n; i++) {
				_m[i] += eta2*H[i]/_w[i];
			}
			mu.changed();
			reset();
		}

		public void oneshot() { reset(); accumulate(1); flush(); }
		public void dispose() { rate2.dispose(); }
		public void starting() { reset(); }
		public void run() { accumulate(1); flush(); }
	}
	
	public class TensionedTrainer extends Trainer
	{
		VDouble		tension=new VDouble("tension",0.01);
		double			lw[]=new double[n];

		public TensionedTrainer() {}

		public String toString() { return "TensionedTrainer:"+DiffScaler.this; }

		public VDouble getTension() { return tension; }

		public void flush() {
			double T=tension.value;
			int i;

			// go through and modify G
			for (i=0; i<n; i++) { lw[i]=Math.log(_w[i]); }
			G[0] -= T*(lw[1]-lw[0]);	i=n-1;
			G[i] -= T*(lw[i-1]-lw[i]); 	i--;
			for (; i>0; i--) {
				G[i] -= T*(lw[i-1] -2*lw[i]+lw[i+1]);
			}

			// do the same for H?
			// may have to modify T.
			H[0] -= T*(_m[1]-_m[0]);	i=n-1;
			H[i] -= T*(_m[i-1]-_m[i]); 	i--;
			for (; i>0; i--) {
				H[i] -= T*(_m[i-1] -2*_m[i]+_m[i+1]);
			}

			super.flush();
		}

		public void dispose() { tension.dispose(); super.dispose(); }
	}
}