Daniel@0: function [Y, Loss] = separationOracleMAP(q, D, pos, neg, k)
Daniel@0: %
Daniel@0: %   [Y,Loss]  = separationOracleMAP(q, D, pos, neg, k)
Daniel@0: %
Daniel@0: %   q   = index of the query point
Daniel@0: %   D   = the current distance matrix
Daniel@0: %   pos = indices of relevant results for q
Daniel@0: %   neg = indices of irrelevant results for q
Daniel@0: %   k   = length of the list to consider (unused in MAP)
Daniel@0: %
Daniel@0: %   Y is a permutation 1:n corresponding to the maximally
Daniel@0: %   violated constraint
Daniel@0: %
Daniel@0: %   Loss is the loss for Y, in this case, 1-AP(Y)
Daniel@0: 
Daniel@0: 
Daniel@0:     % First, sort the documents in descending order of W'Phi(q,x)
Daniel@0:     % Phi = - (X(q) - X(x)) * (X(q) - X(x))'
Daniel@0:     
Daniel@0:     % Sort the positive documents
Daniel@0:     ScorePos        = - D(pos,q);
Daniel@0:     [Vpos, Ipos]    = sort(full(ScorePos'), 'descend');
Daniel@0:     Ipos            = pos(Ipos);
Daniel@0:     
Daniel@0:     % Sort the negative documents
Daniel@0:     ScoreNeg        = - D(neg,q);
Daniel@0:     [Vneg, Ineg]    = sort(full(ScoreNeg'), 'descend');
Daniel@0:     Ineg            = neg(Ineg);
Daniel@0: 
Daniel@0:     % Now, solve the DP for the interleaving
Daniel@0: 
Daniel@0:     numPos  = length(pos);
Daniel@0:     numNeg  = length(neg);
Daniel@0:     n       = numPos + numNeg;
Daniel@0: 
Daniel@0:     
Daniel@0:     % Pre-generate the precision scores
Daniel@0: %     H       = triu(1./bsxfun(@minus, (0:(numPos-1))', 1:n));
Daniel@0:     H       = tril(1./bsxfun(@minus, 0:(numPos-1), (1:n)'));
Daniel@0: 
Daniel@0:     % Padded cumulative Vneg
Daniel@0:     pcVneg  = cumsum([0  Vneg]);
Daniel@0: 
Daniel@0:     % Generate the discriminant scores
Daniel@0:     H       = H + scoreChangeMatrix(Vpos, Vneg, n, pcVneg);
Daniel@0: 
Daniel@0:     % Cost of inserting the first + at position b
Daniel@0:     P       = zeros(size(H));
Daniel@0: 
Daniel@0:     % Now recurse
Daniel@0:     for a = 2:numPos
Daniel@0: 
Daniel@0:         % Fill in the back-pointers
Daniel@0:         [m,p]           = cummax(H(:,a-1));
Daniel@0:         % The best point is the previous row, up to b-1
Daniel@0:         H(a:n,a)        = H(a:n,a) + (a-1)/a .* m(a-1:n-1)';
Daniel@0:         P(a+1:n,a)      = p(a:n-1);
Daniel@0:         P(a,a)          = a-1;
Daniel@0:     end
Daniel@0: 
Daniel@0:     % Now reconstruct the permutation from the DP table
Daniel@0:     Y           = nan * ones(n,1);
Daniel@0:     [m,p]       = max(H(:,numPos));
Daniel@0:     Y(p)        = Ipos(numPos);
Daniel@0: 
Daniel@0:     for a = numPos:-1:2
Daniel@0:         p       = P(p,a);
Daniel@0:         Y(p)    = Ipos(a-1);
Daniel@0:     end
Daniel@0:     Y(isnan(Y)) = Ineg;
Daniel@0: 
Daniel@0:     % Compute loss for this list
Daniel@0:     Loss        = 1 - AP(Y, pos);
Daniel@0: end
Daniel@0: 
Daniel@0: function C = scoreChangeMatrix(Vpos, Vneg, n, pcVneg)
Daniel@0:     numNeg  = length(Vneg);
Daniel@0:     numPos  = length(Vpos);
Daniel@0: 
Daniel@0:     % Inserting the a'th relevant document at position b
Daniel@0:     % There are (b - (a - 1)) negative docs before a
Daniel@0:     % And (numNeg - (b - (a - 1))) negative docs after
Daniel@0:     %
Daniel@0:     % The change in score is proportional to: 
Daniel@0:     %
Daniel@0:     %   sum_{negative j}  (Vpos(a) - Vneg(j)) * y_{aj}
Daniel@0:     %
Daniel@0:     %   = (numNeg - (b - (a - 1))) * Vpos(a)            # Negatives after a
Daniel@0:     %   - (cVneg(end) - cVneg(b - (a - 1)))             Weight of negs after a
Daniel@0:     %   - (b - (a - 1)) * Vpos(a)                       # Negatives before a
Daniel@0:     %   + cVneg(b - (a - 1))                            Weight of negs before a
Daniel@0:     %
Daniel@0:     %   Rearrange:
Daniel@0:     %
Daniel@0:     %   (numNeg - 2 * (b - a + 1)) * Vpos(a)
Daniel@0:     %   - cVneg(end) + 2 * cVneg(b - a + 1)
Daniel@0:     %
Daniel@0:     % Valid range of a:  1:numPos
Daniel@0:     % Valid range of b:  a:n
Daniel@0: 
Daniel@0:     D   = bsxfun(@plus, 1-(1:numPos), (1:n)');
Daniel@0:     C   = numNeg - 2 * D;
Daniel@0:     C   = bsxfun(@times, Vpos, C);
Daniel@0: 
Daniel@0:     D(D < 1)                = 1;
Daniel@0:     D(D > length(pcVneg))   = length(pcVneg);
Daniel@0: 
Daniel@0:     %     FIXME:  2011-01-28 21:13:37 by Brian McFee <bmcfee@cs.ucsd.edu>
Daniel@0:     % brutal hack to get around matlab's screwy matrix reshaping 
Daniel@0:     if numPos == 1
Daniel@0:         pcVneg = pcVneg';
Daniel@0:     end
Daniel@0: 
Daniel@0:     C   = C + 2 * pcVneg(D) - pcVneg(end);
Daniel@0: 
Daniel@0:     % Normalize
Daniel@0:     C   = bsxfun(@ldivide, (1:numPos) * numNeg, C);
Daniel@0: 
Daniel@0:     % -Inf out the infeasible regions
Daniel@0:     C   = C - triu(Inf * bsxfun(@gt, (1:numPos), (1:n)'),1);
Daniel@0: 
Daniel@0: 
Daniel@0: end
Daniel@0: 
Daniel@0: function x = AP(Y, pos)
Daniel@0:     % Indicator for relevant documents
Daniel@0:     rel     = ismember(Y, pos);
Daniel@0: 
Daniel@0:     % Prec@k for all k
Daniel@0:     Prec    = cumsum(rel)' ./ (1:length(Y));
Daniel@0: 
Daniel@0:     % Prec@k averaged over relevant positions
Daniel@0:     x       = mean(Prec(rel));
Daniel@0: end