function eigvector = lde(X, Y, options)

%===== parameter check ================================================
if ~isfield(options,'ReducedDim')
    options.ReducedDim = 15;
end
if ~isfield(options,'t')
    options.t = 1;
end
if ~isfield(options,'intraK')
    options.intraK = 3;
end
if ~isfield(options,'intelK')
    options.intelK = 3;
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;
t = options.t;

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);
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;
% sig1 = 0.5;
sig1 = options.t;
W1 = exp(-S1/sig1);
W2 = exp(-S2/sig1);

% % 0-1 weight
% 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;
XLX1  = X' * L1 * X;
XLX2  = X' * L2 * X;

XLX1  = 0.5*(XLX1 + XLX1');
XLX2  = 0.5*(XLX2 + XLX2');

% XLX1 = XLX1 + 1e-5 * eye(size(XLX1));

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


function [V, D] = eig_large_g(A, B, ell)
[V, D] = eig(A, B);
[dummy, D_idx] = sort(-real(diag(D)));
V = V(:, D_idx(1:ell));
D = diag(D);
D = D(D_idx(1:ell));
return