function [eigvector L1 L2 XLX1 XLX2] = RLDE(X, Y, options)
% RLDE: Supervised LDE + PCA
% LDE: (X'*L2*X) * w = lambda * (X'*L1*X) * w 
% LDEPCA: [ X'*((1-a)*L2+a*I)*X ] * w = lambda * [ (1-a)*(X'*L1*X)+a*diag(X'*L1*X) ] * w 
% a = 0, LDE; a = 1, similar to PCA

%===== parameter check ================================================
if ~isfield(options,'ReducedDim')
    options.ReducedDim = 15;
end
if ~isfield(options,'intraK')
    options.intraK = 5;
end
if ~isfield(options,'intelK')
    options.intelK = 5;
end
if ~isfield(options,'alpha')
    options.alpha = 0.5;
end
if ~isfield(options,'PCARatio') || options.PCARatio == 1
    bPCA = 0;
else
    bPCA = 1;
    PCARatio = options.PCARatio;
end

Dim = options.ReducedDim;
K1 = options.intraK;
K2 = options.intelK;
alpha = options.alpha;

nSmp = size(X,1);
if bPCA
    Xm = X - ones(nSmp,1) * mean(X,1);
    [L, S, V] = svd(Xm);
    if PCARatio > 1
        eigvector_PCA = V(:, 1:PCARatio);
    else
        S = sum(S, 1).^2;
        S_sum = cumsum(S);
        thres = S_sum(end)*PCARatio;
        for i = 2 : length(S)
            if S_sum(i) > thres
                break;
            end
        end
        eigvector_PCA = V(:, 1:i-1);
    end
    X = X*eigvector_PCA;
end

S1 = repmat(inf, nSmp, nSmp);  
S2 = repmat(inf, nSmp, nSmp);
S1 = sparse(S1);  S2 = sparse(S2);
for i = 1 : nSmp
    for j = i : nSmp
        if Y(i) == Y(j)
            S1(i, j) = sum((X(i,:) - X(j,:)).^2);
        else
            S2(i, j) = sum((X(i,:) - X(j,:)).^2);
        end
    end
end
S1 = min(S1,S1');
S2 = min(S2,S2');

[dummy, ind] = sort(S1);
for i = 1 : nSmp
    S1(i, ind((2+K1):end,i)) = inf;
end
S1 = min(S1,S1');

[dummy, ind] = sort(S2);
for i = 1 : nSmp
    S2(i, ind((1+K2):end,i)) = inf;
end
S2 = min(S2,S2');


% sig1 = median(S1((S1(:)~=inf & S1(:)~=0)))*30;  % empirical setting
% sig1 = 0.5;
sig1 = options.t;
W1 = exp(-S1/sig1);
W2 = exp(-S2/sig1);

% % 0-1 weight, not used
% W1 = S1; W1(W1==inf) = 0;  W1(W1~=0) = 1;
% W2 = S2; W2(W2==inf) = 0;  W2(W2~=0) = 1;

D1 = diag(sum(W1, 2));
D2 = diag(sum(W2, 2));
L1 = D1 - W1;
L2 = D2 - W2;

% % LDE + PCA
% XLX1  = (1-alpha)*(X'*L1*X) +  alpha*eye(size(X,2));
% XLX2 = X' * ( (1-alpha)*L2+ alpha*eye(size(L2)) )  * X;

XLX1 = (1-alpha)*(X'*L1*X) + alpha*diag(diag(X'*L1*X));
XLX2 = X' * ( (1-alpha)*L2 + alpha*eye(size(L2)) )  * X;


XLX1 = max(XLX1, XLX1');
XLX2 = max(XLX2, XLX2');

if bPCA
    if size(eigvector_PCA,2) > Dim
        option = struct('disp', 0);
        [eigvector_lde, Ei] = eigs(XLX2, XLX1, Dim, 'la', option);
        eigvector = eigvector_PCA * eigvector_lde;
    else
        eigvector = eigvector_PCA;
    end
else
    option = struct('disp', 0);
    [eigvector_lde, Ei] = eigs(XLX2, XLX1, Dim, 'la', option);
    eigvector = eigvector_lde;
end

% for i = 1:size(eigvector,2)
%     eigvector(:,i) = eigvector(:,i)./norm(eigvector(:,i));
% end



