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.tools.*;
samer@0: import	java.util.*;
samer@0: 
samer@0: public class MOGVector extends NamedTask
samer@0: {
samer@0: 	Vec		input;
samer@0: 
samer@0: 	int		N, M;					// num states, number of inputs
samer@0: 	VVector	s;						// states
samer@0: 	VVector	l;						// likelihoods
samer@0: 	VDouble	L;						// total likelihood
samer@0: 		
samer@0: 	Matrix	t;						// deviations: t=(x-mu)/sig
samer@0: 	Matrix	post;					// posterior over states
samer@0: 	Matrix	w, mu, sig;			// weights, means, stddevs
samer@0: 	Matrix	dw, dmu, dsig;		// accumulated stats for learning
samer@0: 
samer@0: 	VDouble	nuw, numu, nusig;	// learning rates
samer@0: 
samer@0: 
samer@0: 	// ------------ private bits ------------------------
samer@0: 
samer@0: 	private final static double Q=0.5*Math.log(2*Math.PI);
samer@0: 	private double[] tmp;
samer@0: 	private double [] _s;
samer@0: 	private double [] _l;
samer@0: 	private double [][] _t;
samer@0: 	private double [][] _w;
samer@0: 	private double [][] _mu;
samer@0: 	private double [][] _sig;
samer@0: 	private double [][] _post;
samer@0: 	// ---------------------------------------------------
samer@0: 
samer@0: 	public MOGVector(Vec in, int n)
samer@0: 	{
samer@0: 		super("mog");
samer@0: 		input=in;
samer@0: 
samer@0: 		Shell.push(node);
samer@0: 		N=n; M=in.size();
samer@0: 
samer@0: 		s = new VVector("state",M);
samer@0: 		l = new VVector("-log p(s)",M);
samer@0: 		L = new VDouble("likelihood");
samer@0: 
samer@0: 		t		= new Matrix("t",N,M);
samer@0: 		post	= new Matrix("p(s|x)",N,M);
samer@0: 		mu		= new Matrix("means",N,M);
samer@0: 		sig	= new Matrix("sigmas",N,M);
samer@0: 		w		= new Matrix("weight",N,M);
samer@0: 		dmu	= new Matrix("dmu",N,M);
samer@0: 		dsig	= new Matrix("dsig",N,M);
samer@0: 		dw		= new Matrix("dw",N,M);
samer@0: 
samer@0: 		tmp=new double[M];
samer@0: 
samer@0: 		Shell.pop();
samer@0: 
samer@0: 		// initialise parameters
samer@0: 
samer@0: 		// weights...
samer@0: 		w.set(new Constant(1.0/N));
samer@0: 		w.changed();
samer@0: 
samer@0: 		// means...
samer@0: 		for (int i=0; i<N; i++) {
samer@0: 			mu.setMatrix(i,i,0,M-1,new Jama.Matrix(1,M,(double)i));
samer@0: 		}
samer@0: 		mu.changed();
samer@0: 
samer@0: 		// sigmas
samer@0: 		sig.set(new Constant(1));
samer@0: 		sig.changed();
samer@0: 
samer@0: 
samer@0: 		// ----- initialise private bits ----------
samer@0: 		_s=s.array();
samer@0: 		_l=l.array();
samer@0: 		_t=t.getArray();
samer@0: 		_w=w.getArray();
samer@0: 		_mu=mu.getArray();
samer@0: 		_sig=sig.getArray();
samer@0: 		_post=post.getArray();
samer@0: 	}
samer@0: 
samer@0: 	// normalise column sum of w, return array of sums
samer@0: 	private void normalise(double [][] w, double[] sum) {
samer@0: 		// make sure weights are normalised properly
samer@0: 		Mathx.zero(sum);
samer@0: 		for (int i=0; i<N; i++) Mathx.add(sum,w[i]);
samer@0: 		for (int i=0; i<N; i++) Mathx.div(w[i],sum);
samer@0: 	}
samer@0: 
samer@0: 	public VFunction getPDF() { 
samer@0: 		PDF		 fn=new PDF();
samer@0: 		VFunction vfn=new VFunction("pdf",fn);
samer@0: 		fn.setupObservers(vfn);
samer@0: 		return vfn; 
samer@0: 	}
samer@0: 
samer@0: 	public void run() 
samer@0: 	{
samer@0: 		// matlab equivalent code:
samer@0: 		// t = (x - mu)./sig
samer@0: 		// post = (w./sig).*exp(-0.5*t.^2)
samer@0: 		// [dummy, s] = max(post,2)
samer@0: 		// Z = sum(post,2)
samer@0: 		// post = post./(Z*ones(1,N))
samer@0: 		// Z=Z*Q;
samer@0: 
samer@0: 		{
samer@0: 			// Vec.Iterator it=input.iterator();
samer@0: 
samer@0: 			double [] x=input.array();
samer@0: 
samer@0: 			for (int i=0; i<N; i++) {					
samer@0: 				for (int j=0; j<M; j++) {
samer@0: 					_t[i][j] = (x[j]-_mu[i][j])/_sig[i][j];
samer@0: 					_post[i][j] = (_w[i][j]/_sig[i][j])*Math.exp(-0.5*_t[i][j]*_t[i][j]);
samer@0: 				}
samer@0: 			}
samer@0: 
samer@0: 			// this computes partition function 
samer@0: 			// and normalises posterior in one go
samer@0: 			normalise(_post,_l);
samer@0: 
samer@0: 			// get MAP state
samer@0: 			for (int j=0; j<M; j++) {
samer@0: 				int		state=0;	
samer@0: 				double	pmax=_post[0][j];
samer@0: 				
samer@0: 				for (int i=1; i<N; i++) {
samer@0: 					if (_post[i][j]>pmax) { state=i; pmax=_post[i][j]; }
samer@0: 				}
samer@0: 				_s[j]=state;
samer@0: 			}
samer@0: 
samer@0: 			for (int j=0; j<M; j++) _l[j]=-Math.log(_l[j]);
samer@0: 
samer@0: 			L.value = Mathx.sum(_l);
samer@0: 		}
samer@0: 
samer@0: 		L.changed();
samer@0: 		t.changed();
samer@0: 		s.changed();
samer@0: 		l.changed();
samer@0: 		post.changed();
samer@0: 	}
samer@0: 
samer@0: 	public Task learnTask() {
samer@0: 		return new AnonymousTask() {
samer@0: 			// buffer changes to parameters
samer@0: 			double [][] _dw=dw.getArray();
samer@0: 			double [][] _dmu=dmu.getArray();
samer@0: 			double [][] _dsig=dsig.getArray();
samer@0: 
samer@0: 			public void starting() { 
samer@0: 				dw.zero();	
samer@0: 				dmu.zero();	
samer@0: 				dsig.zero();
samer@0: 				normalise(w.getArray(),tmp);
samer@0: 			}
samer@0: 
samer@0: 			public void run() {
samer@0: 				for (int i=0; i<N; i++) {
samer@0: 					for (int j=0; j<M; j++) {
samer@0: 						double pp=_post[i][j], tt=_t[i][j];
samer@0: 						_dw[i][j]  += pp;
samer@0: 						_dmu[i][j] += pp*tt;
samer@0: 						_dsig[i][j]+= pp*(tt*tt-1);
samer@0: 					}
samer@0: 				}
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: 
samer@0: 	public Task flushTask() { 
samer@0: 		Shell.push(node);
samer@0: 
samer@0: 		try {
samer@0: 			return new AnonymousTask() {
samer@0: 
samer@0: 				VDouble nuw=new VDouble("weights.learn.rate",.001);
samer@0: 				VDouble numu=new VDouble("means.learn.rate",.001);
samer@0: 				VDouble nusig=new VDouble("sigmas.learn.rate",.001);
samer@0: 
samer@0: 				double [][] _w=w.getArray();
samer@0: 				double [][] _dw=dw.getArray();
samer@0: 
samer@0: 				public void run()
samer@0: 				{
samer@0: 					// mu += numu*dmu.*sig
samer@0: 					// sig += nusig*dsig.*sig;
samer@0: 					// lambda = sum(w.*dw,2)./sum(w.*w,2)
samer@0: 					//	dw -= (lambda*ones(1,N)) .* w
samer@0: 					// w *= exp(nu*dw)
samer@0: 
samer@0: 					dmu.arrayTimesEquals(sig);
samer@0: 					dmu.timesEquals(nusig.value);
samer@0: 					mu.plusEquals(dmu);
samer@0: 					mu.changed();
samer@0: 					dmu.zero();
samer@0: 
samer@0: 					dsig.arrayTimesEquals(sig);
samer@0: 					dsig.timesEquals(nusig.value);
samer@0: 					sig.plusEquals(dsig);
samer@0: 					sig.changed();
samer@0: 					dsig.zero();
samer@0: 
samer@0: 					{
samer@0: 						// the effect of this is to project
samer@0: 						// dw away from w. the resulting vector
samer@0: 						// is then added to the log of w
samer@0: 
samer@0: 						double nu=nuw.value;
samer@0: 
samer@0: 						for (int j=0; j<M; j++) {
samer@0: 							double w2=0, lambda=0;
samer@0: 							for (int i=0; i<N; i++) {
samer@0: 								lambda += _w[i][j]*_dw[i][j];
samer@0: 								w2     += _w[i][j]*_w[i][j];
samer@0: 							}
samer@0: 							tmp[j]=lambda/w2;
samer@0: 						}
samer@0: 
samer@0: 						for (int i=0; i<N; i++) {
samer@0: 							for (int j=0; j<M; j++) {
samer@0: 								_w[i][j]  *= Math.exp(
samer@0: 									nu*(_dw[i][j] - tmp[j]*_w[i][j])
samer@0: 								);  // update w
samer@0: 							}
samer@0: 						}
samer@0: 					}
samer@0: 
samer@0: 					normalise(w.getArray(),tmp);
samer@0: 					w.changed();
samer@0: 					dw.zero();	
samer@0: 				}
samer@0: 			};
samer@0: 		} finally { Shell.pop(); }
samer@0: 	}
samer@0: 
samer@0: 	class PDF extends Function implements Observer {
samer@0: 		VInteger		index;
samer@0: 		Viewable		vbl=null;
samer@0: 
samer@0: 		public PDF() { 
samer@0: 			index=new VInteger("index",M/2); 
samer@0: 			index.setRange(0,M-1); 
samer@0: 			index.addObserver(this);
samer@0: 		}
samer@0: 
samer@0: 		public void dispose() { index.dispose(); }
samer@0: 
samer@0: 		public double apply(double x) {
samer@0: 			double Z=0;
samer@0: 			for (int i=0; i<N; i++) {
samer@0: 				int	 j=index.value;
samer@0: 				double t=(x-_mu[i][j])/_sig[i][j];
samer@0: 				Z+=(_w[i][j]/_sig[i][j])*Math.exp(-0.5*t*t);
samer@0: 			}
samer@0: 			return Z;
samer@0: 		}
samer@0: 		public String format(String x) {
samer@0: 			return "mogpdf("+x+")";
samer@0: 		}
samer@0: 
samer@0: 		public void setupObservers(Viewable v) {
samer@0: 			w.addObserver(this);
samer@0: 			mu.addObserver(this);
samer@0: 			sig.addObserver(this);
samer@0: 			vbl=v;
samer@0: 		}
samer@0: 
samer@0: 		public void update(Observable o, Object arg) {
samer@0: 			if (arg!=Viewable.DISPOSING) vbl.changed();
samer@0: 		}
samer@0: 	}
samer@0: }