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: import  samer.core.*;
samer@0: import  samer.core.Agent.*;
samer@0: import  samer.core.types.*;
samer@0: import  samer.tools.*;
samer@0: import  samer.maths.*;
samer@0: 
samer@0: public class ICAWithScaler extends ICA
samer@0: {
samer@0: 	VVector  		k;			// state, scaling vector
samer@0: 	double[]		__k, tmp;
samer@0: 
samer@0: 	public ICAWithScaler(Vec input) { this(input.size()); setInput(input); }
samer@0: 	public ICAWithScaler(int N)
samer@0: 	{
samer@0: 		super(N);
samer@0: 		Shell.push(node);
samer@0: 		k = new VVector("k",n);
samer@0: 		Shell.pop();
samer@0: 		__k=k.array();
samer@0: 		tmp=new double[n];
samer@0: 
samer@0: 		for (int i=0; i<n; i++) __k[i]=1;
samer@0: 		k.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void dispose() { k.dispose(); super.dispose(); }
samer@0: 	public void infer() {	// this overrides ICA.infer
samer@0: 		infer.run(); // compute s=Wx
samer@0: 		Mathx.mul(s.array(),__k);
samer@0: 		s.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void compute() {
samer@0: 		Mathx.mul(tmp,sourceModel.getGradient(),__k);
samer@0: 		grad.apply(tmp,_g);
samer@0: 	}
samer@0: 
samer@0: 	public samer.maths.opt.Functionx functionx() { return null; }
samer@0: 
samer@0: 	public void fold() {
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			Mathx.mul(W.getArray()[i],__k[i]);
samer@0: 			__k[i]=1;
samer@0: 		}
samer@0: 		k.changed();
samer@0: 		W.changed();
samer@0: 	}
samer@0: 
samer@0: 	public Trainer getDiffTrainer() { return new DifferentialTrainer(); }
samer@0: 	public Trainer getScaleTrainer() { return new ScalerTrainer(); }
samer@0: 
samer@0: 
samer@0: 	public class DifferentialTrainer extends ON3Trainer {
samer@0: 		VDouble		scaleRate,stretchRate;
samer@0: 		double		lastflush;
samer@0: 
samer@0: 		public DifferentialTrainer() {
samer@0: 			scaleRate=new VDouble("alpha",0.1);
samer@0: 			stretchRate=new VDouble("beta",0.05);
samer@0: 		}
samer@0: 
samer@0: 		public void reset() { super.reset(); lastflush=0; }
samer@0: 		public void flush() { diffFlush(); fold(); super.flush(); lastflush=0; }
samer@0: 		public void diffFlush() // flush multipliers to k instead of W
samer@0: 		{
samer@0: 			double batchlet=batch-lastflush;
samer@0: 			if (batchlet==0) return;
samer@0: 
samer@0: 			// do differential learning on trace & diagonal of G
samer@0: 			double alpha=scaleRate.value/batchlet;
samer@0: 			double beta=stretchRate.value/batchlet;
samer@0: 
samer@0: 			// compute factors and scale each row of W
samer@0: 			double mu=G.trace()/n,dl;
samer@0: 			for (int i=0; i<_n; i++) {
samer@0: 				dl=alpha*mu+beta*(_G[i][i]-mu);
samer@0: 				double tmp=Math.exp(-dl);
samer@0: 				// if (Double.isNaN(tmp)) throw new Error("alt: NaN"+i);
samer@0: 				__k[i]*=tmp; // instead of Mathx.mul(_W[i],tmp);
samer@0: 				_G[i][i]=0;
samer@0: 			}
samer@0: 			k.changed();
samer@0: 			lastflush=batch;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 	/** This one trains ONLY the scaler part, not the ICA part, so is a lot faster
samer@0: 		than using the differential trainer with a zero learning rate. */
samer@0: 	
samer@0: 	public class ScalerTrainer extends AnonymousTask implements Model.Trainer
samer@0: 	{
samer@0: 		VVector		G;
samer@0: 		double[]		_G;
samer@0: 		double[]     _g,_s;
samer@0: 		VDouble		scaleRate,stretchRate;
samer@0: 		int				_n;
samer@0: 		double		batch,thresh;
samer@0: 		
samer@0: 		public ScalerTrainer()
samer@0: 		{
samer@0: 			_n=n;
samer@0: 			G=new VVector("G",n);
samer@0: 			thresh=Shell.getDouble("anomaly",20*n);
samer@0: 			scaleRate=new VDouble("alpha",0.02);
samer@0: 			stretchRate=new VDouble("beta",0.002);
samer@0: 			batch=0;
samer@0: 
samer@0: 			_s=s.array();
samer@0: 			_G=G.array();
samer@0: 		}
samer@0: 
samer@0: 		public void starting() { reset(); }
samer@0: 		public void run() { accumulate(); }
samer@0: 
samer@0: 		public void dispose() { G.dispose(); scaleRate.dispose(); stretchRate.dispose(); super.dispose(); }
samer@0: 		public void oneshot() { accumulate(); flush(); }
samer@0: 		public void reset() { Mathx.zero(_G); batch=0; }
samer@0: 		public void accumulate() { accumulate(1); }
samer@0: 		public void accumulate(double w) {
samer@0: 			// HACK!
samer@0: 			if (sourceModel.getEnergy()>thresh) return; 
samer@0: 			batch+=w;
samer@0: 
samer@0: 			double[] phi=sourceModel.getGradient();
samer@0: 			for (int i=0; i<_n; i++) _G[i] += w*(phi[i]*_s[i] - 1);
samer@0: 		}
samer@0: 
samer@0: 		public void flush()
samer@0: 		{
samer@0: 			if (batch==0) return;
samer@0: 
samer@0: 			G.changed();
samer@0: 			
samer@0: 			// do differential learning on trace & diagonal of G
samer@0: 			double alpha=scaleRate.value/batch;
samer@0: 			double beta=stretchRate.value/batch;
samer@0: 			
samer@0: 			// compute factors and scale each row of W
samer@0: 			double mu=Mathx.sum(_G)/n, dl;
samer@0: 			for (int i=0; i<_n; i++) {
samer@0: 				dl=alpha*mu+beta*(_G[i]-mu);
samer@0: 				double tmp=Math.exp(-dl);
samer@0: 				if (Double.isNaN(tmp)) throw new Error("alt: NaN"+i);
samer@0: 				__k[i]*=tmp; // instead of Mathx.mul(_W[i],tmp);
samer@0: 			}
samer@0: 			k.changed();
samer@0: 		
samer@0: 			reset(); // ready for next batch
samer@0: 		}
samer@0: 	}
samer@0: }
samer@0: 
samer@0: