function A = ERLDA1(X,E,o)
%{
ERLDA = Efficient Regularized Linear Discriminant Analysis
    A = answer
    X = Data observations
    E = Class assignments
    o = sigma = regularization weight
%}
[n p] = size(X); % n = # of observations % p = # of attributes
H_n = eye(n) - (1/n)*ones(n,1)*ones(1,n); % Centering matrix %sec 2.1 %p.3)
H_p = eye(p) - (1/p)*ones(p,1)*ones(1,p);

c = size(E,2); % Number of classes
if c == 1
 % Allows class specification as 1D vector rather than indicator matrix
    Classes = unique(E);
    c = numel(Classes);
    if c < 2
        error('You must have at least 2 classes');
    end
    E = ones(n,1) * Classes' == E * ones(1,c);
end
E = double(E); % Don't want it to be logical

pi = sum(E,1)';
PI = diag(pi); % number of observations in each class (algorithm parameter in the paper)
H_pi = eye(c) - (1/n)*sqrt(pi)*sqrt(pi)';   %p.4

M = inv(PI)*E'*X;   %(3)

%    -----
% Algorithm 1
%    -----

% STEP 2
if p < n % inverting matrix is O(n^3) so we take the smaller one
    %G = inv(X'*H_n*X + (o^2)*eye(p))*X'*H_n*E*PI^(-1/2);      %(13)
    G = inv(X'*H_n*X + (o^2)*eye(p))*M'*sqrt(PI)*H_pi;     %(13 line 2)
else 
    G = X'*H_n*inv(H_n*X*X'*H_n + (o^2)*eye(n))*E*PI^(-1/2);     %(14)
end
R = H_pi*PI^(1/2)*M*G;     %(15)

% STEP 3 
[V_R Gamma_R V_R_transpose] = svd(R); %%% "Thin" SVD?

% STEP 4
%A = G*V_R;
A = G*V_R*Gamma_R^(-1/2);