function [G2, H2, G, Hf, x0] = rlda(example, alpha)
%RLDA: Robust Fisher’s linear discriminant analysis.
%
% 	[G2, H2] = rlda(example, alpha) produces matrix-valued functions G2 and H2
% 	for the second-order NEPv from robust LDA analysis. The NEPv is 
%	in the form of 
%		G2(x)x = lam*H2(x)x,
%	where x is a real vector, G2 and H2 are symmetric matrices of size
%	n, and lam is a smallest positive eigenvalue of the pencil G2(x) - mu*H2(x).
%	This is a test example used in [2].
%
%	INPUT (optional): 
%	example --- {sonar}, iono % two testing data sets 
%	alpha --- sampling rate {0.4}; a number between 0 and 1, indicating
%	the ratio of data points sampled for testing. 
%
%	[G2, H2, G, Hf, x0] = rlda(...) also returns the first-order NEPv
%		G(x)*x = mu *H(x)*x, 
%	where G and H(x) = Hf(x)*Hf(x)' are symmetric matrices of size n. 
%
%	x0 is an approximate eigenvector to be used for iterative methods such as SCF.
%
%
% -------
%	Usage:
%	ex = 'sonar';
%	[G2, H2, G, H, x0] = rec(ex); 
%	G2x0 = G2(x0); % size 60*60
%	H2x0 = H2(x0); % size 60*60
%
%	ex = 'iono';
%	[G2, H2, G, H, x0] = rec(ex); 
%	G2x0 = G2(x0); % size 34*34
%	H2x0 = H2(x0); % size 34*34
%
%
% -------
% 	NEPv features: Real symmetric G and H; Generalized eigenproblem;
% 	Single vector X; Unitarily invariant; Smallest positive eigenvalue;
%
% 	References:
% 	[1] S.-J. Kim, A. Magnani, and S. Boyd, Robust Fisher discriminant analysis, 
%		in Advances in Neural Information Processing Systems 18 (NIPS 2005), 
%		MIT Press, 2006, pp. 659–666.
%	[2] M. Lichman, UCI Machine Learning Repository, 2013,
%		http://archive.ics.uci.edu/ml.
%	[3]	Z. Bai, D. Lu, and B. Vandereycken.  Robust Rayleigh quotient
%		minimization and nonlinear eigenvalue problems,  SIAM J. Matrix
%		Anal. Appl., 40 (2018), pp.A3495-A3522.
%
%	Ding.Lu@uky.edu. Version date July-20-2022.
%

rng(0); % fix random seed for reproducible testing

if nargin < 1, example = 'sonar'; end

if strcmp(example, 'sonar')
    load sonar.mat;
    sonar_label = (sonar_label == 77); % turn to 0-1 label
    DATA = [sonar_data, sonar_label]; % merge data and label for simplicity 
    %ALPHA = [0.3, 0.4, 0.5, 0.6, 0.7, 0.8];
	if nargin < 2, alpha = 0.4; end;

elseif strcmp(example, 'iono')
    load ionosphere.mat;
    ionosphere_label = (ionosphere_label == 103); % turn to 0-1 label
    DATA = [ionosphere_data, ionosphere_label]; % merge data and label for simplicity 
    %ALPHA = [0.10, 0.15, 0.2, 0.4, 0.6, 0.7];
	if nargin < 2, alpha = 0.4; end;
else
	disp("Error: Input example = 'sonar' or 'iono'.");
	return
end

n = size(DATA, 2) - 1;
N_SAMPLE = size(DATA, 1);
N_TRAIN = floor(N_SAMPLE*alpha);
N_RESAMPLETIME = 100;

% Uncomment the following if you want to randomly shuffle the data points
RDPM = randperm(N_SAMPLE);
DATA_TRAIN = DATA(RDPM(1:N_TRAIN),:);

% --------------------------------------------------
% Loop to obtain parameters for the uncertainty set 
% -------------------------------------------------

UX = []; UY = [];
SX = []; SY = [];
for i = 1 : N_RESAMPLETIME
    % uniformly distributed resample
    RESAMPLE_IDX= randi(N_TRAIN, N_TRAIN, 1);
    TRAIN_TMP = DATA_TRAIN(RESAMPLE_IDX,:);
    TRLABEL = TRAIN_TMP(:, end);
    TRAIN_TMP_c0 = TRAIN_TMP(TRLABEL==0,1:n);
    TRAIN_TMP_c1 = TRAIN_TMP(TRLABEL==1,1:n);
    while size(TRAIN_TMP_c0,1) <= 1 || size(TRAIN_TMP_c1,1) <= 1
        % redo if too few sample
        RESAMPLE_IDX= randi(N_TRAIN, N_TRAIN, 1);
        TRAIN_TMP = DATA_TRAIN(RESAMPLE_IDX,:);
        TRLABEL = TRAIN_TMP(:, end);
        TRAIN_TMP_c0 = TRAIN_TMP(TRLABEL==0,1:n);
        TRAIN_TMP_c1 = TRAIN_TMP(TRLABEL==1,1:n);
    end
    ux0 = mean(TRAIN_TMP_c0); % mean of class 0
    sx0 = cov(TRAIN_TMP_c0); % cov of class 0
    uy0 = mean(TRAIN_TMP_c1); % mean  of class 1
    sy0 = cov(TRAIN_TMP_c1); % cov  of class 1

    UX = [UX; ux0];
    UY = [UY; uy0];
    SX = [SX; sx0(:)']; % vectorize cov
    SY = [SY; sy0(:)']; % vectorize cov

end

ux = mean(UX)'; 
uy = mean(UY)'; 
[Px, IPx] = trinv(cov(UX)); %+ 1.0E-8*eye(n);
Px = Px/n; IPx = IPx*n;
[Py, IPy] = trinv(cov(UY)); %+ 1.0E-8*eye(n);
Py = Py/n; IPy = IPy*n;

msx = mean(SX);
dfsx = SX - ones(N_RESAMPLETIME,1)*msx;
dx = sqrt(max(sum(dfsx.*dfsx, 2)));
Sx = reshape(msx, n,n); % + 1.0E-8*eye(n);

msy = mean(SY);
dfsy = SY - ones(N_RESAMPLETIME,1)*msy;
dy = sqrt(max(sum(dfsy.*dfsy, 2)));
Sy = reshape(msy, n,n); %+ 1.0E-8*eye(n);


w = pinv(Sx+Sy)*(ux-uy);
P0 = Sx+Sy+(dx+dy)*eye(n);


% Generate matrices
G2 = @(z) TGMAT(z, P0);
H2 = @(z) THMAT(z, ux, uy, IPx, IPy);

G = @(z) GMAT(z, P0);
Hf =  @(z) FW(z, ux, uy, IPx, IPy);

x0 = pinv(Sx+Sy)*(ux-uy);

% check feasibility of the initial x0
if norm(Hf(x0)) == 0
	disp('Warning: Computed initial vector non-feasible!')
end

return;

end % END OF RCSP



% *******************************************
% Local functions
% *******************************************

% ---------------------------------------------------- %
%   1. Function GMAT compute G 
% ---------------------------------------------------- %
function G = GMAT(w, M)
    G = M;
end

% ---------------------------------------------------- %
%   2. Function TGMAT compute TG 
% ---------------------------------------------------- %
function TG = TGMAT(w, M)
    TG = M;
end


% ---------------------------------------------------- %
%   3. Function HMAT compute H 
% ---------------------------------------------------- %
function H = HMAT(w, ux, uy, IPx, IPy)
	fw = FW(w, ux, uy, IPx, IPy)
	if norm(fw) == 0
	   	H = zeros(length(fw));
	else 
		H = fw*fw';
	end
end

% ---------------------------------------------------- %
%   4. Function THMAT compute tilde H
% ---------------------------------------------------- %
function TH = THMAT(w, ux, uy, IPx, IPy)
    n = length(ux);
    dxy = ux-uy;
    ipxw = IPx*w; wpxw = sqrt(w'*ipxw);
    ipyw = IPy*w; wpyw = sqrt(w'*ipyw);
    if abs(w'*dxy) <= (wpxw + wpyw )
        H = zeros(length(dxy));
        TH = zeros(length(dxy));
        return
    end
    fw = dxy - sign(w'*dxy)*(ipxw/wpxw + ipyw/wpyw);
    H = fw*fw'; H = (H+H')/2;
    DH = -sign(w'*dxy)*( (IPx - ipxw*ipxw'/ (wpxw)^2)/wpxw + (IPy - ipyw*ipyw'/ (wpyw)^2)/wpyw);
    DH = (DH + DH')/2;
    DH = (w'*fw) * DH;
	TH = H + DH;
end


% ---------------------------------------------------- %
%   5. Function FW: compute fw(z), s.t., H = f*f'
% ---------------------------------------------------- %
function fw = FW(w, ux, uy, IPx, IPy)
    n = length(ux);
    dxy = ux-uy;
    ipxw = IPx*w; wpxw = sqrt(w'*ipxw);
    ipyw = IPy*w; wpyw = sqrt(w'*ipyw);
    if abs(w'*dxy) <= (wpxw + wpyw )
        fw = zeros(length(dxy),1);
        return
    end
    fw = dxy - sign(w'*dxy)*(ipxw/wpxw + ipyw/wpyw);
end
 

% ---------------------------------------------------- %
%   6. Function TRINV resets small singular values of M 
%	and return pseudo inverse 
% ---------------------------------------------------- %
function [M, IM] = trinv(M)
tol = 1.0E-8;
[V,E] = eig(M); %M = V*E*V'
e = diag(E);

% set small e to max(e)*tol
e(e<max(e)*tol) = max(e)*tol;
IM = V*diag(e)*V'; IM = (IM+IM')/2;
E = pinv(diag(e));
M = V*E*V'; M = (M+M')/2;
end
