function [eigvector, eigvalue, elapse] = SDA(gnd,fea,LabelIdx,UnlabelIdx,options)
% SDA: Semi-supervised Discriminant Analysis
%
%       [eigvector, eigvalue, elapse] = SDA(gnd,feaLabel,feaUnlabel,options)
% 
%             Input:
%               gnd     - Label vector.  
%               fea     - data matrix. Each row vector of fea is a data point. 
%
%              LabelIdx - fea(LabelIdx,:) is the labeled data matrix. 
%            UnlabelIdx - fea(UnlabelIdx,:) is the unlabeled data matrix. 
%
%               options - Struct value in Matlab. The fields in options
%                         that can be set:
%
%                     WOptions     Please see ConstructW.m for detailed options.  
%                       or
%                      W           You can construct the W outside.
%
%                     beta         Paramter to tune the weight between
%                                  supervised info and local info
%                                    Default 0.1. 
%                                       beta*L+\tilde{I} 
%
%                         Please see LGE.m for other options.
%
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           data point (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The eigvalue of SDA eigen-problem. sorted from
%                           smallest to largest. 
%               elapse    - Time spent on different steps 
% 
%
%    Examples:
%
%       
% 
%
% See also LPP, LGE
%
%Reference:
%
%   Deng Cai, Xiaofei He and Jiawei Han, "Semi-Supervised Discriminant
%   Analysis ", IEEE International Conference on Computer Vision (ICCV),
%   Rio de Janeiro, Brazil, Oct. 2007.   
%
%   version 2.0 --July/2007 
%   version 1.0 --May/2006 
%
%   Written by Deng Cai (dengcai2 AT cs.uiuc.edu)



if ~isfield(options,'ReguType')
    options.ReguType = 'Ridge';
end
if ~isfield(options,'ReguAlpha')
    options.ReguAlpha = 0.1;
end


[nSmp,nFea] = size(fea);
nSmpLabel = length(LabelIdx);
nSmpUnlabel = length(UnlabelIdx);


if nSmpLabel+nSmpUnlabel ~= nSmp
    error('input error!');
end

gnd = gnd(LabelIdx);

classLabel = unique(gnd);
nClass = length(classLabel);
Dim = nClass;


if ~isfield(options,'W') 
    [W, timeW] = constructW(fea,options.WOptions);

    if isfield(options.WOptions,'bSemiSupervised') & options.WOptions.bSemiSupervised
        if ~isfield(options.WOptions,'SameCategoryWeight')
            options.WOptions.SameCategoryWeight = 1;
        end

        G2 = zeros(nSmpLabel,nSmpLabel);
        Label = unique(gnd);
        nLabel = length(Label);
        for idx=1:nLabel
            classIdx = find(gnd==Label(idx));
            G2(classIdx,classIdx) = options.WOptions.SameCategoryWeight;
        end
        W(LabelIdx,LabelIdx) = G2;
    end
else
    W = options.W;
    timeW = 0;
end



tmp_T = cputime;

D = full(sum(W,2));
W = -W;
for i=1:size(W,1)
    W(i,i) = W(i,i) + D(i);
end

beta = 0.1;
if isfield(options,'beta') & (options.beta > 0) 
    beta = options.beta;
end

D = W*beta;
for i=1:nSmpLabel
    D(LabelIdx(i),LabelIdx(i)) = D(LabelIdx(i),LabelIdx(i)) + 1;
end

elapse.timeW = timeW + cputime - tmp_T;


tmp_T = cputime;

%==========================
% If data is too large, the following centering codes can be commented
%==========================
if isfield(options,'keepMean') & options.keepMean
    ;
else
    if issparse(fea)
        fea = full(fea);
    end
    sampleMean = mean(fea,1);
    fea = (fea - repmat(sampleMean,nSmp,1));
end
%==========================


DPrime = fea'*D*fea;

switch lower(options.ReguType)
    case {lower('Ridge')}
        for i=1:size(DPrime,1)
            DPrime(i,i) = DPrime(i,i) + options.ReguAlpha;
        end
    case {lower('Tensor')}
        DPrime = DPrime + options.ReguAlpha*options.regularizerR;
    case {lower('Custom')}
        DPrime = DPrime + options.ReguAlpha*options.regularizerR;
    otherwise
        error('ReguType does not exist!');
end

DPrime = max(DPrime,DPrime');



feaLabel = fea(LabelIdx,:);

Hb = zeros(nClass,nFea);
for i = 1:nClass,
    index = find(gnd==classLabel(i));
    classMean = mean(feaLabel(index,:),1);
    Hb (i,:) = sqrt(length(index))*classMean;
end
WPrime = Hb'*Hb;
WPrime = max(WPrime,WPrime');


elapse.timePCA = cputime - tmp_T;
tmp_T = cputime;


dimMatrix = size(WPrime,2);

if Dim > dimMatrix
    Dim = dimMatrix; 
end


if isfield(options,'bEigs')
    if options.bEigs
        bEigs = 1;
    else
        bEigs = 0;
    end
else
    if (dimMatrix > 1000 & Dim < dimMatrix/10) | (dimMatrix > 500 & Dim < dimMatrix/20) | (dimMatrix > 250 & Dim < dimMatrix/30) 
        bEigs = 1;
    else
        bEigs = 0;
    end
end


if bEigs
    %disp('use eigs to speed up!');
    option = struct('disp',0);
    [eigvector, eigvalue] = eigs(WPrime,DPrime,Dim,'la',option);
    eigvalue = diag(eigvalue);
else
    [eigvector, eigvalue] = eig(WPrime,DPrime);
    eigvalue = diag(eigvalue);
    
    [junk, index] = sort(-eigvalue);
    eigvalue = eigvalue(index);
    eigvector = eigvector(:,index);

    if Dim < size(eigvector,2)
        eigvector = eigvector(:, 1:Dim);
        eigvalue = eigvalue(1:Dim);
    end
end

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


elapse.timeMethod = cputime - tmp_T; 
elapse.timeAll =  elapse.timeW + elapse.timePCA + elapse.timeMethod;