view arrows/mscaler.m @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

% mscaler - Dynamic additive and multiplicative normalisation in Matlab
%
% mscaler ::
%    model(real) ~'model for observations',
%    options {
%       scale   :: nonneg /1  ~'initial scale factor';
%       offset  :: real   /0  ~'initial offset';
%       scale_rate  :: nonneg/0.02 ~'scale adaptation rate';
%       offset_rate :: nonneg/1e-7 ~'offset adaptation rate';
%       batch   :: natural/4  ~'batch size for updates'
%    }
% -> arrow({[[N]]},{[[N]]},mscaler_state).

function o=mscaler(model,varargin)
	opts=options('scale',1,'offset',0, ...
		'scale_rate',0.02,'offset_rate',1e-7,'batch',4, ...
		'nargout', 1, ...
		varargin{:});

	rates=[opts.scale_rate;opts.offset_rate];
	batch=opts.batch;
	score=scorefn(model);
	z3 = [0;0;0];

	if (opts.batch==1)
		o=loop(@update,@(s)[opts.scale;opts.offset]);
	else
		o=loop(@update_batched,@(s){z3,[opts.scale;opts.offset]});
	end

	function ss=stats(y,phi), 
		n=size(y,1); 
		ss=[n;sum(y.*phi)-n;sum(phi)]; 
	end

	function [y,phi]=infer(params,x)
		y = (x-params(2))/params(1);
		phi = score(y);
	end

	function params=updparams(params,stat1,stat2)
		params = [	params(1)*exp(rates(1)*stat1); 
						params(2)+rates(2)*params(1)*stat2 ];
	end

	function [y,state]=update(x,state)
		[y,phi] = infer(state,x);
		state = updparams(state, mean(y(:).*phi(:))-1, mean(phi(:)));
	end

	function [y,state]=update_batched(x,state)
		params=state{2};
		[y,phi] = infer(params,x);
		ss = state{1} + stats(y(:),phi(:));
		if (ss(1)<batch) state{1}=ss;
		else state = {z3, updparams(params,ss(2)/ss(1),ss(3)/ss(1))};
		end
	end
end