wolffd@0: function [Y, Loss] = separationOracleNDCG(q, D, pos, neg, k)
wolffd@0: %
wolffd@0: %   [Y,Loss]  = separationOracleNDCG(q, D, pos, neg, k)
wolffd@0: %
wolffd@0: %   q   = index of the query point
wolffd@0: %   D   = the current distance matrix
wolffd@0: %   pos = indices of relevant results for q
wolffd@0: %   neg = indices of irrelevant results for q
wolffd@0: %   k   = length of the list to consider 
wolffd@0: %
wolffd@0: %   Y is a permutation 1:n corresponding to the maximally
wolffd@0: %   violated constraint
wolffd@0: %
wolffd@0: %   Loss is the loss for Y, in this case, 1-NDCG(Y)
wolffd@0: 
wolffd@0: 
wolffd@0:     % First, sort the documents in descending order of W'Phi(q,x)
wolffd@0:     % Phi = - (X(q) - X(x)) * (X(q) - X(x))'
wolffd@0:     
wolffd@0:     % Sort the positive documents
wolffd@0:     ScorePos        = - D(pos, q);
wolffd@0:     [Vpos, Ipos]    = sort(full(ScorePos'), 'descend');
wolffd@0:     Ipos            = pos(Ipos);
wolffd@0:     
wolffd@0:     % Sort the negative documents
wolffd@0:     ScoreNeg        = - D(neg, q);
wolffd@0:     [Vneg, Ineg]    = sort(full(ScoreNeg'), 'descend');
wolffd@0:     Ineg            = neg(Ineg);
wolffd@0: 
wolffd@0:     % Now, solve the DP for the interleaving
wolffd@0: 
wolffd@0:     numPos  = length(pos);
wolffd@0:     numNeg  = length(neg);
wolffd@0:     n       = numPos + numNeg;
wolffd@0: 
wolffd@0:     % From Chakrabarti (KDD08)
wolffd@0:     k       = min(k, numPos);
wolffd@0: 
wolffd@0:     cVneg   = cumsum(Vneg);
wolffd@0: 
wolffd@0:     Discount = zeros(k, 1);
wolffd@0:     Discount(1:2) = 1;
wolffd@0:     Discount(3:k) = 1./ log2(3:k);
wolffd@0: 
wolffd@0:     DCGstar         = sum(Discount);
wolffd@0: 
wolffd@0: 
wolffd@0:     % Pre-compute the loss table
wolffd@0:     LossTab     = padarray(  hankel(- Discount / DCGstar), ...
wolffd@0:                             max(0, [numNeg numPos] - k), 0, 'post');
wolffd@0:     if sum(size(LossTab) > [numNeg, numPos])
wolffd@0:         LossTab     = LossTab(1:numNeg, 1:numPos);
wolffd@0:     end
wolffd@0: 
wolffd@0:     % 2010-01-17 09:13:41 by Brian McFee <bmcfee@cs.ucsd.edu>
wolffd@0:     % initialize the score table
wolffd@0: 
wolffd@0:     pcVneg = [0 cVneg];
wolffd@0:     % Pre-compute cellScore
wolffd@0:     cellValue = bsxfun(@times, Vpos / (numPos * numNeg), numNeg - 2 * ((1:numNeg)-1)');
wolffd@0:     cellValue = bsxfun(@plus, (2 * pcVneg(1:numNeg) - cVneg(end))' / (numPos * numNeg), cellValue);
wolffd@0:     cellValue = cellValue + LossTab;
wolffd@0: 
wolffd@0:     S       = zeros(numNeg, numPos);
wolffd@0:     P       = zeros(numNeg, numPos);
wolffd@0:     
wolffd@0:     % Initialize first column
wolffd@0:     P(:,1) = 1;
wolffd@0:     S(:,1) = cellValue(:,1);
wolffd@0: 
wolffd@0:     % Initialize first row
wolffd@0:     P(1,:) = 1;
wolffd@0:     S(1,:) = cumsum(cellValue(1,:));
wolffd@0:     
wolffd@0:     % For the rest, use the recurrence
wolffd@0: 
wolffd@0:     for g = 2:numPos
wolffd@0:         [m, pointer]    = cummax(S(:,g-1));
wolffd@0:         P(:,g)          = pointer;
wolffd@0:         S(:,g)          = m' + cellValue(:,g);
wolffd@0:     end
wolffd@0: 
wolffd@0:     % Now reconstruct the permutation from the DP table
wolffd@0:     Y           = nan * ones(n,1);
wolffd@0:     [m,p]       = max(S(:,numPos));
wolffd@0:     
wolffd@0:     Loss        = 1 + LossTab(p,numPos);
wolffd@0:     
wolffd@0:     NegsBefore      = zeros(numPos,1);
wolffd@0:     NegsBefore(numPos)  = p-1;
wolffd@0: 
wolffd@0:     for a = numPos:-1:2
wolffd@0:         p               = P(p,a);
wolffd@0:         NegsBefore(a-1) = p-1;
wolffd@0:         Loss            = Loss + LossTab(p,a-1);
wolffd@0:     end
wolffd@0:     Y((1:numPos)' + NegsBefore)     = Ipos;
wolffd@0:     Y(isnan(Y))                     = Ineg;
wolffd@0: 
wolffd@0: end