wolffd@0: function [Y, Loss] = separationOracleMAP(q, D, pos, neg, k)
wolffd@0: %
wolffd@0: %   [Y,Loss]  = separationOracleMAP(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 (unused in MAP)
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-AP(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:     
wolffd@0:     % Pre-generate the precision scores
wolffd@0: %     H       = triu(1./bsxfun(@minus, (0:(numPos-1))', 1:n));
wolffd@0:     H       = tril(1./bsxfun(@minus, 0:(numPos-1), (1:n)'));
wolffd@0: 
wolffd@0:     % Padded cumulative Vneg
wolffd@0:     pcVneg  = cumsum([0  Vneg]);
wolffd@0: 
wolffd@0:     % Generate the discriminant scores
wolffd@0:     H       = H + scoreChangeMatrix(Vpos, Vneg, n, pcVneg);
wolffd@0: 
wolffd@0:     % Cost of inserting the first + at position b
wolffd@0:     P       = zeros(size(H));
wolffd@0: 
wolffd@0:     % Now recurse
wolffd@0:     for a = 2:numPos
wolffd@0: 
wolffd@0:         % Fill in the back-pointers
wolffd@0:         [m,p]           = cummax(H(:,a-1));
wolffd@0:         % The best point is the previous row, up to b-1
wolffd@0:         H(a:n,a)        = H(a:n,a) + (a-1)/a .* m(a-1:n-1)';
wolffd@0:         P(a+1:n,a)      = p(a:n-1);
wolffd@0:         P(a,a)          = a-1;
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(H(:,numPos));
wolffd@0:     Y(p)        = Ipos(numPos);
wolffd@0: 
wolffd@0:     for a = numPos:-1:2
wolffd@0:         p       = P(p,a);
wolffd@0:         Y(p)    = Ipos(a-1);
wolffd@0:     end
wolffd@0:     Y(isnan(Y)) = Ineg;
wolffd@0: 
wolffd@0:     % Compute loss for this list
wolffd@0:     Loss        = 1 - AP(Y, pos);
wolffd@0: end
wolffd@0: 
wolffd@0: function C = scoreChangeMatrix(Vpos, Vneg, n, pcVneg)
wolffd@0:     numNeg  = length(Vneg);
wolffd@0:     numPos  = length(Vpos);
wolffd@0: 
wolffd@0:     % Inserting the a'th relevant document at position b
wolffd@0:     % There are (b - (a - 1)) negative docs before a
wolffd@0:     % And (numNeg - (b - (a - 1))) negative docs after
wolffd@0:     %
wolffd@0:     % The change in score is proportional to: 
wolffd@0:     %
wolffd@0:     %   sum_{negative j}  (Vpos(a) - Vneg(j)) * y_{aj}
wolffd@0:     %
wolffd@0:     %   = (numNeg - (b - (a - 1))) * Vpos(a)            # Negatives after a
wolffd@0:     %   - (cVneg(end) - cVneg(b - (a - 1)))             Weight of negs after a
wolffd@0:     %   - (b - (a - 1)) * Vpos(a)                       # Negatives before a
wolffd@0:     %   + cVneg(b - (a - 1))                            Weight of negs before a
wolffd@0:     %
wolffd@0:     %   Rearrange:
wolffd@0:     %
wolffd@0:     %   (numNeg - 2 * (b - a + 1)) * Vpos(a)
wolffd@0:     %   - cVneg(end) + 2 * cVneg(b - a + 1)
wolffd@0:     %
wolffd@0:     % Valid range of a:  1:numPos
wolffd@0:     % Valid range of b:  a:n
wolffd@0: 
wolffd@0:     D   = bsxfun(@plus, 1-(1:numPos), (1:n)');
wolffd@0:     C   = numNeg - 2 * D;
wolffd@0:     C   = bsxfun(@times, Vpos, C);
wolffd@0: 
wolffd@0:     D(D < 1)                = 1;
wolffd@0:     D(D > length(pcVneg))   = length(pcVneg);
wolffd@0: 
wolffd@0:     %     FIXME:  2011-01-28 21:13:37 by Brian McFee <bmcfee@cs.ucsd.edu>
wolffd@0:     % brutal hack to get around matlab's screwy matrix reshaping 
wolffd@0:     if numPos == 1
wolffd@0:         pcVneg = pcVneg';
wolffd@0:     end
wolffd@0: 
wolffd@0:     C   = C + 2 * pcVneg(D) - pcVneg(end);
wolffd@0: 
wolffd@0:     % Normalize
wolffd@0:     C   = bsxfun(@ldivide, (1:numPos) * numNeg, C);
wolffd@0: 
wolffd@0:     % -Inf out the infeasible regions
wolffd@0:     C   = C - triu(Inf * bsxfun(@gt, (1:numPos), (1:n)'),1);
wolffd@0: 
wolffd@0: 
wolffd@0: end
wolffd@0: 
wolffd@0: function x = AP(Y, pos)
wolffd@0:     % Indicator for relevant documents
wolffd@0:     rel     = ismember(Y, pos);
wolffd@0: 
wolffd@0:     % Prec@k for all k
wolffd@0:     Prec    = cumsum(rel)' ./ (1:length(Y));
wolffd@0: 
wolffd@0:     % Prec@k averaged over relevant positions
wolffd@0:     x       = mean(Prec(rel));
wolffd@0: end