function [G2, H2, G, H, x0] = rgec(example)
%RGEC: Robust generalized eigenvalue classifier.
%
% 	[G2, H2] = rgec(example) produces matrix-valued functions G2 and H2
% 	for the second-order NEPv from robust generalized eigenvalue
% 	classification. 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).
%
%	example = 'rnd' (default) returns an NEPv of size 3 for randomly generated data;
%	example = 'pima' returns an NEPv of size 9 for the Pima Indians Diabetes data set [2].
%	Both are testing examples used in [3].
%
%	[G2, H2, G, H, x0] = rgec(...) also returns the first-order NEPv
%		G(x)*x = mu *H(x)*x, 
%	where G and H are symmetric matrices of size n. 
%
%	x0 is an approximate eigenvector to be used for iterative methods such as SCF.
%
%
% -------
%	Usage:
%	ex = 'rnd';
%	[G2, H2, G, H, x0] = rec(ex); 
%	G2x0 = G2(x0); % size 3*3
%	H2x0 = H2(x0); % size 3*3
%
%	ex = 'pima';
%	[G2, H2, G, H, x0] = rec(ex); 
%	G2x0 = G2(x0); % size 9*9
%	H2x0 = H2(x0); % size 9*9
%
% -------
% 	NEPv features: Real symmetric G and H; Generalized eigenproblem;
% 	Single vector X; Unitarily invariant; Smallest positive eigenvalue;
%
% 	References:
% 	[1] O. L. Mangasarian and E. W. Wild, Multisurface proximal support
% 		vector machine classification via generalized eigenvalues, IEEE
% 		Trans. Pattern Anal. Mach. Intell., 27 (2005), pp. 1–6. 
%	[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.
%

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

if strcmp(example, 'rng')
	load rng.mat; 

	G2 = @(x) G2MAT(x, A, SIGA);
	H2 = @(x) H2MAT(x, B, SIGB);

	G = @(x) GMAT(x, A, SIGA);
	H = @(x) HMAT(x, B, SIGB);

	AA= [A,-ones(length(A),1)]; AA = AA'*AA;
	BB= [B,-ones(length(B),1)]; BB = BB'*BB;
	[V,E] = eig( AA, BB );
	[~, idx] = min(diag(E));
	x0 = V(:,idx); x0 = x0/norm(x0);
	
elseif strcmp(example, 'pima')

	% pima is an n-by-9 matrix
	% n is the number of data entries
	% 9 is the number of attributes
	% Each Attribute: (all numeric-valued)
	%   1. Number of times pregnant
	%   2. Plasma glucose concentration a 2 hours in an oral glucose tolerance test
	%   3. Diastolic blood pressure (mm Hg)
	%   4. Triceps skin fold thickness (mm)
	%   5. 2-Hour serum insulin (mu U/ml)
	%   6. Body mass index (weight in kg/(height in m)^2)
	%   7. Diabetes pedigree function
	%   8. Age (years)
	%   9. Class variable (0 or 1)
	
	load pima.mat;
	n = length(pima);   % #of entries
	nf= size(pima,2)-1; % #of features
	e = mean(pima(:,1:nf)); % mean 
	nn = round(.7*n); % number of training instances, first 70% data

	% Define ellipsoids 
	SA = eye(nf)/diag(e.^2)/(.5^2); 
	SA(end,end) = SA(end,end)/(.002^2);    % the age is not perturbed
	SA(1,1) = SA(1,1)/(.002^2);            % the number of pregnant
	SB = eye(nf)/diag(e.^2)/(.5^2); 
	SB(end,end) = SB(end,end)/(.002^2);    % the age is not perturbed
	SB(1,1) = SB(1,1)/(.002^2);            % the number of pregnant

    %
    % Generate training (ratio=tstr) and testing data
    %
    % pima = pima(randperm(n),:);    % apply random row permulation to database
    T = pima(1:nn,:);      % training data
    C = pima(nn+1:end,:);  % testing data
    D = sortrows(T,nf+1);   % Sort the training data according to the class label
    nA = length(find(D(:,end)==0)); % #data entry in class A, labeled 0
    nB = length(D) - nA;    % #data entry in class B, labeled 1
    A = D(1:nA,:);
    B = D(nA+1:end,:);
    SIGA = kron(ones(1,nA),SA);
    SIGB = kron(ones(1,nB),SB);

	G2 = @(x) G2MAT(x, A(:,1:nf), SIGA);
	H2 = @(x) H2MAT(x, B(:,1:nf), SIGB);

	G = @(x) GMAT(x, A(:,1:nf), SIGA);
	H = @(x) HMAT(x, B(:,1:nf), SIGB);
	AA = [A(:,1:nf),-ones(nA,1)]; AA = AA'*AA;
	BB = [B(:,1:nf),-ones(nB,1)]; BB = BB'*BB;
	[V,E] = eig( AA, BB );
	[~, idx] = min(diag(E));
	x0 = V(:,idx); x0 = x0/norm(x0);

else
	disp("Error: Input example = 'rng' or 'pima'.");
	return
end

return;
end % END OF RGEC



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

% -------------------------------------------
%   1. Function GMAT evaluates the G(z1) matrix.
% -------------------------------------------
function [AA, TA] = TAMAT(z1, A, SIGA) 
n = size(A,2);
m = size(A,1);
for j = 1 : m	% Compute AA = G(z1)
   S = SIGA(:,(j-1)*n+1:j*n);
   tp = S\z1(1:n);
   TA(j,:) = sign( -z1(end) + A(j,:)*z1(1:n) ) * tp.' / sqrt(z1(1:n)'*tp);
end
TA = A + TA;

eA = -ones(m,1);
AA= [TA,eA]'*[TA,eA];

end % END OF TAMAT

% -------------------------------------------
%   2. Function GMAT evaluates the G(z1) matrix.
% -------------------------------------------
function AA = GMAT(z1, A, SIGA) 
	[AA, ~] = TAMAT(z1, A, SIGA);
end % END OF GMAT

% -------------------------------------------
%   3. Function G2MAT evaluates tilde G(z1) matrix.
% -------------------------------------------
function G2 = G2MAT(z1, A, SIGA)
[AA, TA] = TAMAT(z1, A, SIGA);
n = size(A,2);
m = size(A,1);
eA = -ones(m,1);
Y = [TA,eA]*z1; 	% Compute G2
G2=[];
for jj = 1 : n
    ej = zeros(n,1); ej(jj) = 1;
    for j = 1 : m
        S = SIGA(:,(j-1)*n+1:j*n);
        tp = S\z1(1:n);
        TTA(j,:) = sign( -z1(end) + A(j,:)*z1(1:n) )/ sqrt(z1(1:n)'*tp)...
            * ( (S\ej)' - tp'*ej/(z1(1:n)'*tp)*tp' );
    end
    G2 = [G2,[TTA,zeros(m,1)]'*Y];
end
G2 = [G2,zeros(n+1,1)];
G2 = 0.5*(G2+G2') + AA;
end % END OF G2MAT

% -------------------------------------------
%   4. Function TBMAT evaluates the TB(z1) matrix.
% -------------------------------------------
function [BB, TB] = TBMAT(z1, B, SIGB)
n = size(B,2);
p = size(B,1);
RB = []; TB = []; 	% Compute BB = H(z1)
for j = 1 : p
   S = SIGB(:,(j-1)*n+1:j*n);
   tp = S\z1(1:n);
   % check if S and hyperplane connected
   % i.e., sqrt(w' S^-1 w) < gamma + w^T a_i
   if ( sqrt(z1(1:n)'*tp) > abs(-z1(end)+B(j,:)*z1(1:n)) )
      % RB = [RB; B(j,:)];
      % TB = [TB; ( -B(j,:)*z1(1:n)+z1(end) ) * tp.' / (z1(1:n)'*tp)];
      continue;
   else
       RB = [RB; B(j,:)];
       TB = [TB; sign( -B(j,:)*z1(1:n)+z1(end) ) * tp.' / sqrt(z1(1:n)'*tp)];
   end
end
TB = RB + TB;
eB = -ones(size(TB,1),1); 
BB= [TB,eB]'*[TB,eB];
if length(TB) == 0
   error('Error: All class B ellipsoids intersects plane H(z)=0!')
end 

end % END OF TBMAT

% -------------------------------------------
%   5. Function HMAT evaluates the H(z1) matrix.
% -------------------------------------------
function BB = HMAT(z1, B, SIGB)
 [BB, ~] = TBMAT(z1, B, SIGB);
end % END OF HMAT

% -------------------------------------------
%   6. Function H2MAT evaluates tilde H(z1) matrix.
% -------------------------------------------
function H2 = H2MAT(z1, B, SIGB) 
[BB, TB] = TBMAT(z1, B, SIGB);
n = size(B,2);
p = size(B,1);
eB = -ones(size(TB,1),1);
Y = [TB,eB]*z1; 	% Compute H2
H2=[];
for jj = 1 : n
    TTB=[];
    ej = zeros(n,1); ej(jj) = 1;
    for j = 1 : p
        S = SIGB(:,(j-1)*n+1:j*n);
        tp = S\z1(1:n);
        if ( sqrt(z1(1:n)'*tp) > abs(-z1(end)+B(j,:)*z1(1:n)) )
            % delete the column
            continue;
        end
        tt = sign( -B(j,:)*z1(1:n)+z1(end) )  / sqrt(z1(1:n)'*tp) ...
            * ( (S\ej)' - tp'*ej/(z1(1:n)'*tp)*tp' );
        TTB = [TTB; tt];
    end
    H2 = [H2,[TTB,zeros(size(TTB,1),1)]'*Y];
end
H2 = [H2,zeros(n+1,1)]; 
H2 = 0.5*(H2+H2') + BB;
end % END OF H2MAT
