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: public class Mixture extends NamedTask implements Model
samer@0: {
samer@0: 	private Model 		M[]; 	// models
samer@0: 	private int			n, m; // size of vector, num models
samer@0: 	private Vec			x;		// input
samer@0: 	private VVector	w;		// prior weights
samer@0: 	private VVector	s;		// posterior
samer@0: 	private int			k; 	// MAP estimate
samer@0: 	private VDouble	Z;		// Parition function, ie p(x)
samer@0: 	private double[]  _x,_s,_w,_g;
samer@0: 
samer@0: 	public Mixture( Vec input, int m) { this(input.size(), m); setInput(input); }
samer@0: 	public Mixture( int N, int L)
samer@0: 	{
samer@0: 		super("mixture");
samer@0: 		Shell.push(node);
samer@0: 
samer@0: 		n = N;
samer@0: 		m = L;
samer@0: 
samer@0: 		x = null;
samer@0: 		w = new VVector("prior",m);
samer@0: 		s = new VVector("posterior",m);
samer@0: 		Z = new VDouble("Z");
samer@0: 		M = new Model[m];
samer@0: 		Shell.pop();
samer@0: 	
samer@0: 		_s=s.array();
samer@0: 		_w=w.array();
samer@0: 		_g=new double[n];
samer@0: 		Mathx.set(_w,new Constant(1.0/L));
samer@0: 	}
samer@0: 
samer@0: 	public VVector prior() 		{ return w; }
samer@0: 	public VVector posterior() { return s; }
samer@0: 	public void setModel(int i, Model m) { M[i]=m;	}
samer@0: 	public void setInput(Vec in) { x=in; _x=x.array(); }
samer@0: 	public int getSize() { return n; }
samer@0: 
samer@0: 	public void dispose()
samer@0: 	{
samer@0: 		s.dispose();
samer@0: 		w.dispose();
samer@0: 		Z.dispose();
samer@0: 		for (int i=0; i<m; i++) M[i].dispose();
samer@0: 		super.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public void infer() {
samer@0: 		// get models to compute energies.
samer@0: 		// for (int i=0; i<m; i++) { M[i].infer(); M[i].compute(); }
samer@0: 
samer@0: 		// compute relative posterior
samer@0: 		for (int i=0; i<m; i++) _s[i] = M[i].getEnergy(); // collect energies
samer@0: 		double Emin=Mathx.min(_s);
samer@0: 		for (int i=0; i<m; i++) _s[i] = _w[i]*Math.exp(Emin-_s[i]);
samer@0: 
samer@0: 		// compute partition function, normalise posterior
samer@0: 		Z.set(Mathx.sum(_s));	// compute parition function
samer@0: 		k=Mathx.argmax(_s);	// get MAP model
samer@0: 		Mathx.mul(_s,1/Z.value);	// normalise posterior
samer@0: 		s.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void compute()
samer@0: 	{
samer@0: 		/* compute gradients weighted by posterior */
samer@0: 		Mathx.zero(_g);
samer@0: 		for (int i=0; i<m; i++) {
samer@0: 			double [] phi = M[i].getGradient();
samer@0: 			for (int j=0; j<n; j++) _g[j] += _s[i]*phi[j];
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	public double	getEnergy() { return -Math.log(Z.value); }
samer@0: 	public double [] getGradient() { return _g; }
samer@0: 
samer@0: 	public Functionx functionx() { return null; }
samer@0: 
samer@0: 	public void run() { infer(); }
samer@0: 
samer@0: 	public Trainer getTrainer() { return new Trainer(); }
samer@0: 
samer@0: 	public class Trainer implements Model.Trainer
samer@0: 	{
samer@0: 		Model.Trainer	T[];
samer@0: 		VDouble		rate;
samer@0: 		VVector		dw;
samer@0: 		double		batch, _dw[];
samer@0: 
samer@0: 		public Trainer() {
samer@0: 			T=new Model.Trainer[m]; // should all be null
samer@0: 			rate=new VDouble("rate",0.001);
samer@0: 			dw=new VVector("dw",m);
samer@0: 			_dw=dw.array();
samer@0: 		}
samer@0: 
samer@0: 		public void setTrainer(int i,Model.Trainer t) { T[i]=t; }
samer@0: 		public void dispose() { rate.dispose(); dw.dispose(); }
samer@0: 
samer@0: 		public void accumulate() { accumulate(1.0); }
samer@0: 		public void accumulate(double w) {
samer@0: 			batch+=w;
samer@0: 			for (int i=0;i<m; i++) {
samer@0: 				if (T[i]!=null) T[i].accumulate(w*_s[i]); // sweet
samer@0: 			}
samer@0: 
samer@0: 			// now accumulate info about priors
samer@0: 			Mathx.add(_dw,_s);
samer@0: 		}
samer@0: 
samer@0: 		public void oneshot() { accumulate(1.0); flush(); }
samer@0: 		public void flush() {
samer@0: 			for (int i=0; i<m; i++) if (T[i]!=null) T[i].flush();
samer@0:          double lambda=Mathx.dot(_w,_dw)/Mathx.dot(_w,_w);
samer@0: 			double nu=rate.value/batch;
samer@0: 
samer@0: 			dw.changed();
samer@0: 			for (int i=0; i<m; i++) {
samer@0: 				_w[i] *= Math.exp(nu*(_dw[i]-lambda*_w[i])); // update w
samer@0: 			}
samer@0: 			Mathx.zero(_dw); batch=0;
samer@0: 
samer@0: 			// normalise
samer@0: 			Mathx.mul(_w,1/Mathx.sum(_w));
samer@0: 			w.changed();
samer@0: 		}
samer@0: 		public void reset() { 
samer@0: 			for (int i=0; i<m; i++) if (T[i]!=null) T[i].reset();
samer@0: 			Mathx.zero(_dw);
samer@0: 			batch=0;
samer@0: 		}
samer@0: 	}
samer@0: }
samer@0: