samer@0: package samer.models;
samer@0: 
samer@0: import samer.maths.*;
samer@0: import samer.maths.opt.*;
samer@0: import samer.tools.*;
samer@0: import samer.core.*;
samer@0: import samer.core.types.*;
samer@0: import java.util.*;
samer@0: 
samer@0: /**
samer@0: 	Non-adaptive generalised exponential factorial prior: the pointy
samer@0: 	bit of the usual GeneralisedExponential has been smoothed out
samer@0: 	by blending with a quadratic.
samer@0:  */
samer@0: 
samer@0: public class SmoothGeneralisedExponential extends NamedTask implements Model, Observer {
samer@0: 	Vec		input;
samer@0: 	VVector	alpha, e, grad;
samer@0: 	int			N;
samer@0: 	VDouble	E;
samer@0: 	double[]	x, g, e0, _e, a, eps;
samer@0: 	double	last_eps;
samer@0: 
samer@0: 	public SmoothGeneralisedExponential(Vec x) { this(x.size()); setInput(x); }
samer@0: 	public SmoothGeneralisedExponential(int n)  {
samer@0: 		super("genExp");
samer@0: 		Shell.push(node);
samer@0: 		N=n;
samer@0: 		E=new VDouble("E");
samer@0: 		e=new VVector("e",N);
samer@0: 		grad=new VVector("phi",N);
samer@0: 		alpha=new VVector("alpha",N);
samer@0: 		eps=new double[N];
samer@0: 		alpha.addSaver();
samer@0: 		Shell.pop();
samer@0: 
samer@0: 		alpha.addObserver(this);
samer@0: 		g=grad.array(); // new double[N];
samer@0: 		e0=new double[N];
samer@0: 		a=alpha.array();
samer@0: 		_e=e.array();
samer@0: 		
samer@0: 		Mathx.setAll(a,1.0); 
samer@0: 		setEps(0.1);
samer@0: 	}
samer@0: 
samer@0: 	public void update(Observable obs, Object arg) {
samer@0: 		setEps(last_eps);
samer@0: 	}
samer@0: 	public void setEps(double e) {
samer@0: 		last_eps=e;
samer@0: 		for (int i=0; i<N; i++) { eps[i]=Math.pow(e,2-a[i]); }
samer@0: 	}
samer@0: 	public void setInput(Vec in) { input=in; x=input.array(); }
samer@0: 	public int getSize() { return N; }
samer@0: 	public void dispose() {
samer@0: 		alpha.dispose();
samer@0: 		grad.dispose();
samer@0: 		E.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public VVector	getEnergyVector() { return e; }
samer@0: 	public VDouble	getEnergySignal() { return E; }
samer@0: 	public double	getEnergy() { return E.value; }
samer@0: 	public double [] getGradient() { return g; 	}
samer@0: 	public VVector  getAlpha() { return alpha; }
samer@0: 	
samer@0: 	public void run() { compute(); }
samer@0: 	public void infer() {}
samer@0: 	public void compute() {
samer@0: 		// compute log likelihood
samer@0: 
samer@0: 		for (int i=0; i<N; i++) {
samer@0: 			double t=Math.abs(x[i]);
samer@0: 			_e[i] = t*t/(Math.pow(t,2-a[i]) + eps[i]);
samer@0: 		}
samer@0: 
samer@0: 		// compute gradient g_i = dL/dx_i
samer@0: 		for (int i=0; i<N; i++) {
samer@0: 			if (x[i]==0) g[i]=0;
samer@0: 			else g[i] = a[i]*(_e[i]/x[i]);
samer@0: 		}
samer@0: 
samer@0: 		e.changed();
samer@0: 		grad.changed();
samer@0: 		E.set(Mathx.sum(_e));
samer@0: 	}
samer@0: 
samer@0: 	public Functionx functionx() {
samer@0: 		return new Functionx() {
samer@0: 			double [] __e=new double[N];
samer@0: 			public void 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++) {
samer@0: 					if (x[i]==0) { g[i]=0; __e[i]=0; }
samer@0: 					else {
samer@0: 						double t=Math.abs(x[i]);
samer@0: 						__e[i] = t*t/(Math.pow(t,2-a[i]) + eps[i]);
samer@0: 						g[i] = a[i]*(__e[i]/x[i]);
samer@0: 					}
samer@0: 				}
samer@0: 				return  Mathx.sum(__e); // +L0;
samer@0: 			}
samer@0: 		};
samer@0: 	}
samer@0: }