function a = topNrankagree(refsim, candsim, N, alpha_r, alpha_c)
% a = topNrankagree(refsim, candsim, N, alpha_r, alpha_c)
%     Calculate the "top-N rank agreement score" as defined in 
%     Berenzweig, Logan, Ellis, Whitman @ ISMIR2003
%     refsim is the square reference similarity matrix;
%     candsim is the square candidata similarity matrix;
%     Calculate their normalized similarity between 0 and 1.
%     Defaults: N=10, alpha_r = 0.5^(1/3), alpha_c = 0.5^(2/3)
% 2003-08-12 dpwe@ee.columbia.edu

% Constants
if nargin < 3
  % How far down the ranking to look?
  N = 10;
end
if nargin < 4
  % Decay constant for reference
  alpha_r = 0.5^(1/3);
end
if nargin < 5
  % Decay constant for candidate
  alpha_c = 0.5^(2/3);
end

nr = size(refsim,1);

agree = 0;

for i = 1:nr
  
  % Calculate rankings under each sim
  refsims = refsim(i,[1:nr] ~= i);
  cndsims = candsim(i,[1:nr] ~= i);

%  [dummy, reforder] = sort(-refsims);
%  [dummy, cndorder] = sort(-cndsims);
  % Randomly permute them before sorting to jumble ties
  rp1 = randperm(nr-1);
  rp2 = randperm(nr-1);
  [dummy, refix] = sort(-refsims(rp1));
  [dummy, cndix] = sort(-cndsims(rp2));
  reforder = rp1(refix);
  cndorder = rp2(cndix);

  % Where in the candidate ranking does the artist ranked i by the reference 
  % appear?
  [dummy, cndinv] = sort(cndorder);
  ref2cnd = cndinv(reforder);
  % Calculate the agreement for this row
  agree = agree + sum( (alpha_r.^(1:N)) .* (alpha_c.^ref2cnd(1:N)) );
  
end

a = agree/nr/sum((alpha_r*alpha_c).^(1:N));