function eigvector = SELDLPP(YL, XL, XU, options)
% Semisupervised LDA + LPP

%===== parameter check ================================================
if ~isfield(options,'ReducedDim')
    options.ReducedDim = 15;
end
if ~isfield(options,'t')
    options.t = 1;  % heat kernel parameter in LPP
end
if ~isfield(options,'k')
    options.k = 5;  % knn parameter in LPP
end
Dim = options.ReducedDim;
k = options.k;
t = options.t;

[nSmpLabel,   nFea] = size(XL);
[nSmpUnlabel, nFea] = size(XU);
nSmp = nSmpLabel + nSmpUnlabel;
X = [XL; XU];

% Compute between-class scatter matrix
classLabel = unique(YL);  nClass = length(classLabel);
Hb = zeros(nClass,nFea);
for i = 1 : nClass
    index = find(YL==classLabel(i));
    classMean = mean(XL(index,:),1);
    Hb(i,:) = sqrt(length(index))*classMean;
end
Sb = Hb'*Hb;
Sb = max(Sb,Sb'); 

% Compute within-class scatter matrix
St = XL'*XL;  St = max(St,St');  % total scatter matrix St
Sw = St - Sb;   Sw = max(Sw,Sw');

% using unlabeled samples to construct the graph D and L
% Construct neighborhood graph
options = [];  options.Metric = 'Euclidean';
options.NeighborMode = 'KNN';       options.k = k;
options.WeightMode = 'HeatKernel';  options.t = t;
W = constructW(XU, options);

% Construct diagonal matrix D
D = diag(sum(W, 2));

% Compute Laplacian
L = D - W;
L(isnan(L)) = 0; D(isnan(D)) = 0;

% Compute SELDup and SELDdn
SELDup = Sb + XU'*D*XU;
SELDdn = Sw + XU'*L*XU;

SELDup = max(SELDup, SELDup');
SELDdn = max(SELDdn, SELDdn');

option = struct('disp', 0);
[eigvector_seldlpp, Ei] = eigs(SELDup, SELDdn, Dim, 'la', option);
eigvector = eigvector_seldlpp;

