%*************************************************************************
% This repeats the second experiment (ranking) on USPS in the paper
% "Learning with Partially Absorbing Random Walks" by Xiao-Ming Wu et al.
% Please report problems to xmwu@ee.columbia.edu 
%
% We use all 9268 samples of USPS and use each sample as a query. 
% We compare 3 methods: 1) Columns of PARW with \Lambda being random
% 2) Rows of the same PARW, and 3) PARW with \Lambda = \alpha I (symmetric).
% Last Updated: Jan 30, 2014
%*************************************************************************

clc,clear,close all;
addpath 'util'

alpha = 1e-6;                   % parameter for PARW 
MAP = [];                       % matrix for the MAP results

extAllusps;                     % Extract all USPS samples 
npt = length(labels);           % number of points
ncls = length(unique(labels));  % number of classes

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;

d = sum(W,2);               
L = diag(d) - W;                % Laplacian matrix
L = sparse(L);

rd = rand(npt,1);               % random positive 

for k = 1:ncls
    
    nq = ndgt(k);                        % number of queries
    idx_q = find(labels==k);             % index of the queries
    s = sparse(idx_q, 1:nq,1,npt,nq);
    
    %--------------------- Rows of PARW -----------------------------------
    
    %  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_row -> s = (A')^(-1)*rank_row -> rank_row = ((A')^(-1))\s 
    
    lmbd = alpha * rd;                               % set the \Lambda 
    
    % ranking by rows of PARW
    rank_row = ((sparse(1:npt,1:npt,lmbd)+L)*sparse(1:npt,1:npt,1./lmbd))\s;    
    
    % ranking by columns of the same PARW, note that (L + \Lambda)^{-1} is symmetric
    rank_col = rank_row./repmat(rd,[1,nq]);
    
    [~, idx] = sort(rank_row,'descend');
    out_labels = labels(idx);                       % labels of the ranking points
    [idx_k,~] = find(out_labels==k);                % index of the points of label k in the ranking
    idx_k = reshape(idx_k,nq,nq);
    prec = repmat((1:nq)',[1,nq])./idx_k;           % precision of ranking
    MAP(1,k) = mean(mean(prec),2);                  % mean average precision (MAP)
    clear lmbd rank_row idx idx_k out_labels prec ;
    
    %--------------------- Columns of PARW -----------------------
    
    [~, idx] = sort(rank_col,'descend');
    out_labels = labels(idx);                       % labels of the ranking points
    [idx_k,~] = find(out_labels==k);                % index of the points of label k in the ranking
    idx_k = reshape(idx_k,nq,nq);
    prec = repmat((1:nq)',[1,nq])./idx_k;           % precision of ranking
    MAP(2,k) = mean(mean(prec),2);                  % mean average precision (MAP)
    clear lmbd rank_col idx idx_k out_labels prec ;
    
    %--------------------- PARW with \Lambda = \alpha I (symmetric) -----------------------
    
    lmbd =  alpha * ones(npt,1);                    % set the \Lambda 
    rank = (sparse(1:npt,1:npt,lmbd)+L)\s;          % compute the ranking scores
    [~, idx] = sort(rank,'descend');
    out_labels = labels(idx);                       % labels of the ranking points
    [idx_k,~] = find(out_labels==k);                % index of the points of label k in the ranking
    idx_k = reshape(idx_k,nq,nq);
    prec = repmat((1:nq)',[1,nq])./idx_k;           % precision of ranking
    MAP(3,k) = mean(mean(prec),2);                  % mean average precision (MAP)
    clear lmbd rank idx idx_k out_labels prec ;
    
end

MAP(:,11)=mean(MAP,2),

