%****************************************************************************
% This repeats the first 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 5 methods: 1) ARW-N-1NN (the proposed), 2) ARW-1NN, 3) ARW-CMN, 
% 4) Manifold Ranking (LGC), and 5) PARW with \Lambda = \alpha I.
% Last Updated: Jan 30, 2014
%*****************************************************************************

clc,clear,close all;
addpath 'util'

extAllusps;

alpha_I = 1e-6;                  % parameter for PARW with \Lambda = \alpha I
alpha_LGC = 0.9;                 % parameter for LGC

npt = length(labels);            % number of points
unil = unique(labels);           % unique labels
nl = 20;                         % number of labels
nr = 100;                         % number of runs

dm = squareform(pdist(X));       % compute the distance matrix
sdm = dm.*dm;
clear dm;

sig = 0.2*averagekmin_dm(sdm,20);
W = exp(-sdm/sig);
W(1:size(W,1)+1:end) = 0;

KNN = 20;                       
G = graph_knn(sdm, KNN);        % build the K-NN graph
if(~connectedG(G))
    error('Graph not connected. Please use more neighbors.');
end
W = W.*G;
clear G sdm;

acc_ARW_N_1NN = zeros(nr,1);      % accuracy of ARW-N-1NN proposed in the paper
acc_ARW_1NN = zeros(nr,1);        % accuracy of ARW-1NN
acc_ARW_CMN = zeros(nr,1);        % accuracy of ARW-CMN
acc_PARW = zeros(nr,1);           % accuracy of PARW with \Lambda = \alpha I
acc_LGC = zeros(nr,1);            % accuracy of Manifold Ranking (LGC)

for n = 1:nr
    
    Y = labeling_rd(labels,nl);               % randomly sample labeled points
    
    idxA = [1:npt];                           % all the index
    idxL = find(sum(Y,2)==1);                 % index of labeled points
    idxA(idxL) = [];
    relabels = [labels(idxL);labels(idxA)];   % reorder labels, put labeled points at the top
    lll = relabels(1:nl);                     % labels of labeled points
    
    X2 = [X(idxL,:);X(idxA,:)];               % reorder X
    [W2,Y2,xx] = reorder_hf(W,Y);             % reorder W, form Y2
    
    d = sum(W2,2);
    L = diag(d) - W2;                         % Laplacian matrix
    L = sparse(L);
    
    nd = d.^(-1/2);
    nW = W2.*(nd*nd');
    mL = eye(npt)- alpha_LGC*nW;              % LGC matrix
    mL = sparse(mL);
    s = [eye(nl); zeros(npt-nl,nl)];          % label matrix (idicator of starting vertices)
    
    %%%%------------------- Aborbing Random Walks -------------------
    
    [ARW_CMN, ARW_1NN, ARW_N_1NN] = harmonic_function(W2, Y2);
    
    rank = [eye(nl); ARW_N_1NN];
    [~, idx] = max(rank,[],2);
    results = lll(idx);
    acc_ARW_N_1NN(n)  = sum(results==relabels)/npt;                
    clear rank idx results;
    
    rank = [eye(nl); ARW_1NN];
    [~, idx] = max(rank,[],2);
    results = lll(idx);
    acc_ARW_1NN(n)  = sum(results==relabels)/npt;                
    clear rank idx results;
    
    rank = [Y2; ARW_CMN];
    [~, idx] = max(rank,[],2);
    results = unil(idx);
    acc_ARW_CMN(n)  = sum(results==relabels)/npt;                
    clear rank idx results;
   
    clear ARW_CMN ARW_1NN ARW_N_1NN ;
    
    %%%%------------------- Manifold Ranking (LGC) -------------------------------------------
        
    rank = mL\s;
    [~, idx] = max(rank*Y2,[],2);
    results = unil(idx);
    acc_LGC(n) = sum(results==relabels)/npt;                
    clear rank idx results;
    
    %%%%------------------- PARW with \Lambda = \alpha I --------------------------------------
    
    %  Let A be the npt*npt absorption probability matrix, 
    %  then A = (\Lambda + L)^{-1} \Lambda .
    %  Let s be the npt*nq indicator matrix,
    %  then A' * s = rank -> s = (A')^(-1)*rank -> rank = ((A')^(-1))\s 
    
    lmbd = ones(npt,1)* alpha_I;                % set the \Lambda 
    rank = (sparse(1:npt,1:npt,lmbd)+L)\s;      % compute the rank w.r.t. labels (first nl rows of A)
    [~, idx] = max(rank,[],2);                  % select the largest 
    results = lll(idx);                        
    acc_PARW(n)= sum(results==relabels)/npt; % classification accuracy               
    clear lmbd rank idx results;
    
    clear W2 Y2 nW L mL;
    
end

A(1) = mean(acc_ARW_N_1NN);
A(2) = mean(acc_ARW_1NN);   
A(3) = mean(acc_ARW_CMN);     
A(4) = mean(acc_LGC);      
A(5) = mean(acc_PARW),  



