% [a,t,rank] = recall1(X,Y,g,r) Recall@R
%
% This implements the recall@R measure in information retrieval: the average
% rate of queries for which the 1-nearest neighbor is ranked in the top R
% positions.
% We have a (large) dataset X of D-dim vectors where we want to search for
% queries, and a collection of D-dim query vectors Y. For each query y(m) we
% are given the index g(m) in X of its true nearest neighbor (ground-truth).
% We want to know whether x(g(m)) is among the R nearest neighbors of y(m) in
% X.
% This is useful when the vectors X and Y are an approximate version of the
% original vectors. Typically, the original vectors are high-dimensional and X
% and Y are low-dimensional and binary (as for binary hashing), so that the
% nearest neighbors in X are an approximation to the true nearest neighbors in
% the high-dimensional dataset.
% This can be used with the datasets for approximate nearest neighbor search
% in http://corpus-texmex.irisa.fr (such as SIFT1M or SIFT1B).
% 
% Notes:
% - We use the squared Euclidean distance. This equals the Hamming distance
%   if the vectors are binary, but can also be used when the vectors are
%   real-valued.
% - The algorithm below is somewhat inefficient. A more efficient version
%   would prune the search as follows. For each query y, we scan X, computing
%   the distance to y and keeping track of its current rank, that is,
%   incrementing the rank every time dist(x,y) < d (where d is the
%   ground-truth distance for y). As soon as the rank exceeds max(r) we can
%   stop the search. For queries for which the ground truth is ranked poorly
%   this will speed up the search considerably. If the distances in (X,Y)
%   approximate well the ground-truth distances then there will be little
%   speed up over the worst case, which is O(D.N.M). The vectorized code
%   below runs the worst case all the time but this is efficient in Matlab.
%   (It does not search for queries for which the ground-truth distance is 0,
%   since the rank for these is 1.)
% - If your dataset X doesn't fit in RAM memory, call recall1 on subsets and
%   add the ranks. For example:
%     d = sum((Y-X(g,:)).^2,2); % ground-truth distances
%     load X1; [a,t,rank] = recall1(X1,Y,d,r); R = rank;
%     load X2; [a,t,rank] = recall1(X2,Y,d,r); R = R + rank - 1;
%     load X3; [a,t,rank] = recall1(X3,Y,d,r); R = R + rank - 1;
%     t = bsxfun(@le,R,r); a = mean(t,1);
%   You should modify recall1.m to use the ground-truth distance d directly
%   rather than the index g.
%
% In:
%   X: NxD matrix of N row D-dim vectors, the base vectors (where we search).
%   Y: MxD matrix of M row D-dim vectors, the query vectors.
%   g: Mx1 list of indices, the ground-truth nearest neighbors of Y in X.
%   r: 1xR list of R values to try.
% Out:
%   a: 1xR list where a(i) is the average recall@r(i) in [0,1] over all the
%      M queries.
%   t: MxR logical matrix where t(m,i) is 1 iff x(g(m)) (the ground-truth
%      nearest neighbor of y(m)) is among the r(i) nearest neighbors of y(m)
%      in X.
%   rank: Mx1 list where rank(m) is the rank (order statistic) of x(g(m)) in
%      sorted distance of y(m) to every vector in X.

% Copyright (c) 2015 by Miguel A. Carreira-Perpinan

function [a,t,rank] = recall1(X,Y,g,r)

N = size(X,1); M = size(Y,1);

% Block size (see below), select as large as your memory allows
m = 1; B = floor((m*1024^3)/(4*N*8)/2);	% This should fit in m GB RAM

% Compute ground-truth distances
d = sum((Y-X(g,:)).^2,2);

% If the ground-truth distance is 0 then the rank is 1, no need to search
nz = find(d>0);
rank = ones(M,1);

if ~isempty(nz)
  % Compute distances and ranks
  %rank = sum(bsxfun(@lt,sqdist(X,Y),d'),1)+1; % Single pass, memory-consuming
  % Process X by blocks to save memory
  i1 = 1; i2 = min(N,B);
  x = sum(X.^2,2); y = sum(Y(nz,:).^2,2)';
  while i1 <= N
    rank(nz) = rank(nz) +...
        sum(bsxfun(@lt,bsxfun(@plus,x(i1:i2),...
            bsxfun(@minus,y,2*X(i1:i2,:)*Y(nz,:)')),d(nz)'),1)';
    i1 = i1 + B; i2 = min(N,i1+B-1);
  end
end

% Average recall
t = bsxfun(@le,rank,r);
a = mean(t,1);

