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

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

[nSmpLabel,   nFea] = size(XL);
[nSmpUnlabel, nFea] = size(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 NPE M
% Compute the weight in NPE
tol = 1e-12;
Distance = EuDist2(XU, [], 0);
[sorted, index] = sort(Distance, 2);
neighborhood    = index(:, 2:(1+k) );
W = zeros(k, nSmpUnlabel);
for ii  = 1 : nSmpUnlabel
    z = XU(neighborhood(ii,:),:)-repmat(XU(ii,:),k,1); % shift ith pt to origin
    C = z*z';                                         % local covariance
    C = C + eye(size(C))*tol*trace(C);                % regularlization
    W(:,ii) = C\ones(k,1);                            % solve Cw=1
    W(:,ii) = W(:,ii)/sum(W(:,ii));                   % enforce sum(w)=1
end

M = sparse(1:nSmpUnlabel, 1:nSmpUnlabel, ones(1,nSmpUnlabel), nSmpUnlabel, nSmpUnlabel, 4*k*nSmpUnlabel);
for ii = 1 : nSmpUnlabel
    w = W(:,ii);
    jj = neighborhood(ii,:)';
    M(ii,jj) = M(ii,jj) - w';
    M(jj,ii) = M(jj,ii) - w;
    M(jj,jj) = M(jj,jj) + w*w';
end
M = max(M,M');  M = sparse(M);

% Compute SELDup and SELDdn
SELDup = Sb + XU'*eye(nSmpUnlabel)*XU;
SELDdn = Sw + XU'*M*XU;

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

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

