function [G2, H2, G, H, x0] = rcsp()
%RCSP: Robust common spatial pattern analysis.
%
% 	[G2, H2] = rcsp(ex) produces matrix-valued functions G2 and H2
% 	for the second-order NEPv from robust common spatial pattern
% 	analysis for a random linear mixing model. 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].
%
%	[G2, H2, G, H, x0] = rcsp(...) 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:
%	[G2, H2, G, H, x0] = rcsp(); 
%	G2x0 = G2(x0); % size 10*10
%	H2x0 = H2(x0); % size 10*10
%
%
% -------
% 	NEPv features: Real symmetric G and H; Generalized eigenproblem;
% 	Single vector X; Unitarily invariant; Smallest positive eigenvalue;
%
% 	References:
% 	[1] M. Kawanabe, W. Samek, K.-R. Muller, and C. Vidaurre.
%		Robust common spatial filters with a maxmin approach, 
%		Neural Comput., 26 (2014), pp. 349–376.
%	[2]	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.
%

% Generate data from a random linear mixing model
rng(0); % fix random seed for reproducible testing
[V1, lam1, Sbar1, V2, lam2, Sbar2, Tst1, Tst2, Train1, Train2] = datagen();
delta = 6; % perturbation level, large value may lead to nonsmooth problems 

G2 = @(w) G2MAT(w, lam2, V2, delta, Sbar2);
H2 = @(w) H2MAT(w, lam1, V1, delta, Sbar1);

G = @(w) Sbar2 + build_optmat(w, lam2, V2, delta);
H = @(w) NegTrunc(Sbar1 - build_optmat(w, lam1, V1, delta));

[V, E] = eig(Sbar2, Sbar1);
[~, idx] = min(diag(E));
x0 = V(:,idx); x0 = x0/norm(x0);

return;
end % END OF RCSP



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

% ---------------------------------------------------- %
%   1. Function datagen generates random experiment data 
%	for a linear mixing model.
% ---------------------------------------------------- %

function [V1, lam1, Sbar1, V2, lam2, Sbar2, Tst1, Tst2, Train1, Train2] = datagen()

NUMTRIALS_TRAIN = 50;
TRAIN_SAMPLE_PERTRIAL = 200;
TST_SAMPLE_PERTRIAL = 200;

EPSLN_TR = 2;
EPSLN_TST = 30;

% generate rand rotation
[A,~] = qr(rand(10));

% generate sample data
for i = 1:NUMTRIALS_TRAIN
	% sample covariance in trial i
	X1 = []; X2 = [];
	Y1 = []; Y2 = [];
	for j = 1:TRAIN_SAMPLE_PERTRIAL
		s1 = [randn*sqrt(1.8); randn*sqrt(0.6); randn(8,1)];
		dice = rand;
		XX = A*s1 + randn*sqrt(EPSLN_TR);
		X1 = [X1, XX];

		s1 = [randn*sqrt(1.8); randn*sqrt(0.6); randn(8,1)];
		dice = rand;
		XX = A*s1 + randn*sqrt(EPSLN_TST);
		Y1 = [Y1, XX];

		s2 = [randn*sqrt(0.2); randn*sqrt(1.4); randn(8,1)];
		dice = rand;
		XX = A*s2 + randn*sqrt(EPSLN_TR);
		X2 = [X2, XX];

		s2 = [randn*sqrt(0.2); randn*sqrt(1.4); randn(8,1)];
		dice = rand;
		XX = A*s2 + randn*sqrt(EPSLN_TST);
		Y2 = [Y2, XX];
	end
	S1{i} = cov(X1');%X1 * X1' / (SAMPLE_PERTRIAL - 1);
	S2{i} = cov(X2');%X2 * X2' / (SAMPLE_PERTRIAL - 1);

	Train1{i} = X1; Train2{i} = X2;
	Tst1{i} = Y1; Tst2{i} = Y2;
end

[V1, Sbar1, lam1] = PCA_V(S1, 10);
[V2, Sbar2, lam2] = PCA_V(S2, 10);
return
end %END OF datagen


% ---------------------------------------------------- %
%   2.  Function PCA_V applies PCA for the V matrices 
% ---------------------------------------------------- %

function [V, Sbar1, lam] = PCA_V(S1, n)
%  truncate the tailing n-k1 singular values.
%  you can also set tol to truncate singular values less than |S1|*tol.

% compute Vs by SVD
S1_mat = cell2mat(S1);
S1_vec = reshape(S1_mat, n*n, []);
Sbar1 = mean(S1_vec')'; 
Sbar1 = reshape(Sbar1, n,n);

CovM = cov(S1_vec');
[VV, E] = eig(CovM);
lam1 = diag(E);

k1 = 10; % fixed k1, may also use the following
% tol = eps;
% k1 = nnz(lam1 > tol*lam1(end));

lam = mat2cell(lam1(end-k1+1:end)', 1, ones(1,k1));
VV = VV(:, end-k1+1:end);
V1 = reshape(VV, n, []);

% symmetrize (can also use tril part to represent matrices)
V1T = V1'; V1T = reshape(V1T, n, k1, n); 
V1T = permute(V1T, [1, 3, 2]); 
V1T = reshape(V1T, n, n*k1);
V1 = (V1+V1T)/2;
V = mat2cell(V1, n, n*ones(1,k1));
return
end %END OF PCA_V


% ---------------------------------------------------- %
%   3. Function G2MAT compute G and tilde G
% ---------------------------------------------------- %

function G2 = G2MAT(w, lam2, V2, delta, Sbar1)
[G, DG] = build_optmat(w, lam2, V2, delta);
G = Sbar1 + G;
G2 = G + DG;
end %END OF G2MAT


% ---------------------------------------------------- %
%   4. Function H2MAT compute H and tilde H
% ---------------------------------------------------- %

function H2 = H2MAT(w, lam1, V1, delta, Sbar1)
[H, DH] = build_optmat(w, lam1, V1, delta);
%H = Sbar1 - H;
[H, info] = NegTrunc(Sbar1 - H); % suggested in [Kawanabe, et al., 2014]
if ~info, disp('nonsmooth point in H matrix'), end;
DH = -DH;
H2 = H + DH;
end %END OF H2MAT


% ---------------------------------------------------- %
%   5. Function NegTrunc truncate the negative eig mode in H
% ---------------------------------------------------- %

function [H issmth] = NegTrunc(H)
[V, D] = eig(H);
H = V*max(D,zeros(size(D)))*V';
issmth = (min(diag(D)) >= 0);
end %END OF NegTrunc


% ---------------------------------------------------- %
%   6. Function build_optmat
% ---------------------------------------------------- %

function [P, dP] = build_optmat(w, lam, V, delta)
% build_optmat computes the optimal perturbation matrix, for
%   generating the matrices
%       G(w) = Sigma + build_optmat(w, lam, V, delta)
%       H(w) = Sigma - build_optmat(w, lam, V, delta)
%
%   Input: lam{i} = lambda_i, V{i} = V_i.

P = zeros(size(V{1}));

% evaluate wtVw{i} = w'*V{i}*w
wtVw_fun = @(Vi) w'*Vi*w;
wtVw = cellfun(wtVw_fun, V);
nn = sqrt( wtVw.^2*cell2mat(lam)' );

for i = 1:length(V);
    alpha = delta*(lam{i}*wtVw(i))/nn;
    P = P+alpha*V{i};
end

if nargout < 2
    return;
end

% evaluate dP
S1 = 0; S2 = 0; S3 = 0;
for i = 1:length(V)
    S1 = S1 + (V{i}'*w)*(w'*V{i})*lam{i};
    S2 = S2 + wtVw(i)*V{i}'*w*lam{i};
end

dP = 2*delta*(S1/nn-S2*S2'/nn^3);

return
end %END OF BUILD_OPTMAT
