%****************************************************************************
% This repeats the third experiment (SSL) on USPS in the paper
% "Analyzing the Harmonic Structure in Graph-Based Learning"
% by Xiao-Ming Wu et al. in NIPS 2013
% Please report problems to xmwu@ee.columbia.edu 
%
% We use all 9268 instances of USPS and randomly sample 20 instances as 
% labeled data and make sure there is at least one label for each class. 
% We compare 3 methods: 1) Hitting time (many to one / from the labeled
% points to hit the unlabeled point) 2) Hitting time of the opposite direction 
% (one to many / from the unlabeled point to hit the unlabeled points) 
% 3) the Moore-Penrose pseudoinverse of Graph Laplacian.
% Last Updated: Jan 30, 2014
%*****************************************************************************


clc,clear,close all;
addpath 'util'

extAllusps;
npt = length(labels);         % number of points
unil = unique(labels);        % unique labels
k = length(unil);             % number of classes
nl = 20;                      % number of labeled samples
nr = 100;                     % number of runs

dm = squareform(pdist(X));    % compute the distance matrix
KNN = 20;                       
W = graph_knn(dm, KNN);       % build the unweighted K-NN graph
if(~connectedG(W))
    error('Graph not connected. Please use more neighbors.');
end
clear dm;

d = sum(W,2);
L = diag(d) - W;                   % Laplacian matrix
pL = pinv(full(L));                % Moore-Penrose pseudoinverse of Laplacian matrix
clear L;

nd = d.^(-1/2);
nW = W.*(nd*nd');
nL = speye(npt) - nW;              % normalized Laplacian matrix
clear W nW;
pnL = pinv(full(nL));              % Moore-Penrose pseudoinverse of the normalized Laplacian matrix
clear nL;

acc_hit1 = zeros(nr,1);            % accuracy of method 1
acc_hit2 = zeros(nr,1);            % accuracy of method 2
acc_L = zeros(nr,1);               % accuracy of method 3

for n = 1:nr
    
    Y = labeling_rd(labels,nl);                       % randomly select nl labeled points
    lv = lm2lv(Y);                                    % convert the label matrix to a label vector
    idx = find(lv);
    lll = labels(idx);                                % labels of the labeled instances
    
    %**********************************************************************
    % The hitting time formula:
    % H_{ij} \propto \langle \frac{1}{\sqrt{d_j}}e_j, L_{sym}^\dag (\frac{1}{\sqrt{d_j}}e_j - \frac{1}{\sqrt{d_i}}e_i) \rangle
    % where L_{sym}^\dag = d^{-1/2} L d^{-1/2} is the normalized Laplacian matrix
    
    % 1) Compare hitting times (from all the labeled points to hit one unlabeled point)
    F = -bsxfun(@rdivide,pnL(:,idx),sqrt(d(idx)'));   % note we only compute the part that matters for comparison
    [~,idxx] = min(F,[],2);
    results = lll(idxx);                              % classification results
    acc_hit1(n) = sum(results==labels)/npt;           % accuracy
    clear F idxx results
    
    % 2) Compare hitting times (from one unlabeled point to hit all the labeled points)
    F = bsxfun(@rdivide,pnL(idx,:),sqrt(d(idx)));
    F = bsxfun(@rdivide,F,sqrt(d)');
    lidx = sub2ind(size(pnL),idx,idx);                % convert to linear index
    tmp = pnL(lidx)./d(idx);
    F = bsxfun(@plus,-F,tmp);
    F = F';
    [~,idxx] = min(F,[],2);
    results = lll(idxx);
    acc_hit2(n) = sum(results==labels)/npt;  
    clear F lidx idxx results
    
    % 3) the pseudo-inverse of Graph Laplacian.
    F = pL(:,idx);
    [~,idxx] = max(F,[],2);   
    results = lll(idxx);
    acc_L(n) = sum(results==labels)/npt;
    clear F idxx results idx
end

A(1) = mean(acc_hit1);
A(2) = mean(acc_hit2);
A(3) = mean(acc_L),


