function A = RR_SVD(X,E,sigsq,runctrl);
%% RR_SVD: the RR-SVD algo of Zhang et. al. for regularized LDA %FOLDUP
% 
% A = RR_SVD(X,E,sigsq,runctrl);
% A = RR_SVD(X,classnum,sigsq,runctrl);
%
% does the LDA algorithm 1. of Zhang et. al.
%
% input:
%  X                      n x p matrix; each row an independent observation;
%                           this should not have nans in it!
%  E                      n x c 0/1 indicator matrix; rows correspond to X;
%                           sum of each row is 1; 
%             * or *
%  classnum               n x 1 vector of integers corresponding to the
%                           classes;
%  sigsq                  the regularization term; defaults to 1;
%  runctrl                a struct with fields; none yet...
%   .scale_by_sv            whether to scale by the sqrt singular values at the
%                           end; defaults to true;
% output:
%  A                      p x l matrix transforming into a space
%                           of 'good' separation of X; that is
%                           for new vector v, comparing v'A to XA is
%                           a good separator;
%
% global:
% env var:
% nb: 
% See also: Z. Zhang, G. Dai, M. I. Jordan, 'A Flexible and Efficient Algorithm for 
% Regularized Fisher Discriminant Analysis,' 2009.
% www.cs.berkeley.edu/~jordan/papers/zhang-dai-jordan-ecml09.pdf
% todo:
% changelog: 
%
% Created: 2009.08.04
% Copyright: Cerebellum Capital, 2009-2009
% Author: Steven E. Pav
% Comments: Steven E. Pav
%% %UNFOLD

%#ok<*WNTAG,*CTCH,*NOSEM,*LERR> 			%for matlab mlint; look away.

%% check input %FOLDUP
error(nargchk(2,4,nargin));

[n,p] 			= size(X);

%there are 2 modes of input E: the classnums, or a characteristic matrix;
%we guess which one we have here:

if (isvector(E))
	%these are the class numbers; convert them to a characteristic matrix:
	[clasvals,clasrep,whichclas] 	= unique(E(:));
	E 				= zeros(n,numel(clasvals));
	E(sub2ind(size(E),(1:n)',whichclas)) 	= 1;
end

%we now should have a characteristic function;
my_assert(size(X,1) == size(E,1),'X, E same # of rows');

c 					= size(E,2);

if (islogical(E))
	E 			= double(E);
end

my_assert(all(ismember(E(:),[0,1]) | isnan(E(:))),'E is 0/1/nan');
my_assert(all(sum(E,2) == 1),'E is a characteristic matrix');

%fill in value for sigsq
if (~exist('sigsq','var') || isempty(sigsq))
	sigsq	= 1;
end

%fill in value for runctrl
if (~exist('runctrl','var') || ~isstruct(runctrl))
	runctrl		= struct();
end

runctrl			= setDefaultValues(runctrl,...
	'scale_by_sv',true);
%UNFOLD

%the number of elements in each class:
PiDiag 			= nan_sum(E,1);

%note that we do not get any help from SMW for (13) & (14) b/c they are already
%in the most rank reduced form. that is, for p < n, using (13), we have to
%invert a p x p matrix, which is better than inverting an n x n; SMW is not
%much help here, b/c we have a rank-n update to a pxp invertible identity. 

%note that centering is idempotent & H is symmetric, so
%H'HX = HHX = HX,
%with the last equality b/c you are re-centering a centered matrix;
%so X'HHX = X'HX;
%thus it suffices to center the X now and be done with it;
X 		= bsxfun(@minus,X,mean(X,1));   
%we never had to construct H, rather we compute HX implicitly;

%similarly it suffices to normalize E
E 		= bsxfun(@rdivide,E,sqrt(PiDiag));

%this can be done faster w/ accumarray, I believe, for sufficiently large n,p,c;
XE 		= X'*E;

if (n < p)
	%use (14)
	invP 	= (X*X' + sigsq * eye(n,n)) \ E;
	G 		= X' * invP;
	R 		= XE' * G;
	%R 		= (XE'*X') * invP; 			%this might be faster; needs profiling;
else
	%use (13)
	G 		= (X'*X + sigsq * eye(p,p)) \ XE;
	R 		= XE' * G;
end

%thin SVD for a symmetric matrix:
%symmetrize sneakily, compute the eigenvalues;
[V,gam] 	= eig(0.5 * (R+R'));

%the first eigenvalue in gam is essentially zero; we drop it;
A 			= G * V(:,2:end);

if (runctrl.scale_by_sv)
	scalby 	= sqrt(diag(gam));
	scalby 	= scalby(2:end);
	A 		= bsxfun(@rdivide,A,scalby(:)');
end

end %function
%%% helper functions:
function my_assert(check,msg);
%% my_assert: assertion checking 

if (~check)
	fprintf(2,'Assertion Failure: %s\n',msg);
	throwAsCaller(MException('MATLAB:Utilities:AssertionFailure','Assertion FAILS: %s',msg));
end

end %function
function params = setDefaultValues(params,varargin);
%% setDefaultValues: if struct doesn't have a value, set it;

for i=1:floor(nargin/2)
	if (~isfield(params,varargin{2*i-1}))
		params.(varargin{2*i-1}) = varargin{2*i};
	end
end

end %function