auction_match(distMatrix)

function [assignment,r,p,u,v,value] = auction_match(distMatrix)
%       auction_match:   Compute optimal (maximal) weighted assignment
%       and the corresponding "lattice of dual prices" supporting the
%       optimal assignment.
%		auction_match(disMatrix) computes the optimal assignment for the
%		given rectangular value matrix, for example the assignment
%		of bidders (in rows) to objects (in columns) and vice versa.

%		[assignment,r,p,u,v,value] = ASSIGNMENTOPTIMAL(DISTMATRIX) returns the assignment
%		vector (in assignment) and the overall value (in value) and
%		v: surplus of columns if columns were bidding for rows.
%		u: the corresponding prices of rows.
%		p: prices for columns if rows were bidding for columns
%		r: the corresponding surplus of rows.
%

%       Note that (p,-r) forms the lower corner and (v,-u) forms the
%       upper corner in the lattice of optimal dual vector supporting 
%       the optimal assignment thus giving the complete lattice.
%       Ref. the survey "From the Assignment Model to Combinatorial Auctions" 
%       by S. Bikhchandani and J. Ostroy


%		This is update of the assignment code by Markus Buehren which used Munkres
%       Algorithm for MINIMAL weighted matching. A description of Munkres algorithm
%      (also called Hungarian algorithm) can easily be found on the web.


%% Author : Anuj Kumar www.columbia.edu/~ak2108


% save original distMatrix for value computation
originalDistMatrix    = distMatrix;
% get matrix dimensions
[nOfRows, nOfColumns] = size(distMatrix);
r=zeros(nOfRows,1);
p=zeros(nOfColumns,1);
u=r; v=p;
% memory allocation
coveredColumns = zeros(1,nOfColumns);
coveredRows    = zeros(nOfRows,1);
starMatrix     = zeros(nOfRows, nOfColumns);
primeMatrix    = zeros(nOfRows, nOfColumns);

%%%preliminary steps
if nOfRows < nOfColumns
    error('nOfRows < nOfColumns .. Please enter the disMatrix transpose, you will get what you are looking for')
else % nOfRows > nOfColumns
    minDim = nOfColumns;

    % find the smallest element of each column
    maxVector = max(distMatrix);
    v=v+maxVector';

    % subtract the smallest element of each column from the column
    distMatrix =  distMatrix - repmat(maxVector, nOfRows, 1);

    % Steps 1 and 2

    for col = 1:nOfColumns
        for row = find(distMatrix(:,col)==0)'
            if ~coveredRows(row)
                starMatrix(row, col) = 1;
                coveredColumns(col)  = 1;
                coveredRows(row)     = 1;
                break
            end
        end
    end
    coveredRows(:) = 0; % was used auxiliary above
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Compute the (p,-r) supporting this initial assignment by computing
    % the shortest path in the residual network corresponding to this
    % initial assignment.
    for i=1:nOfRows+nOfColumns+1
        distance(i) = - inf;
    end
    distance(1) = 0;
    for loop1=1:nOfRows+nOfColumns+1
        for i=1:nOfRows
            if (0 < distance(1+i))
                distance(1+i) = 0;
            end
        end
        for i=1:nOfRows
            for j=1:nOfColumns
                if starMatrix(i,j)==0
                    if (distance(i+1)~= -inf)
                        new_distance = distance(i+1) + originalDistMatrix(i,j);
                        if (new_distance > distance(1+nOfRows+j))
                            distance(1+nOfRows+j) = new_distance;
                        end
                    end
                elseif starMatrix(i,j)==1
                    if (distance(1+nOfRows+j)~=-inf)
                        new_distance = distance(1+nOfRows+j) - originalDistMatrix(i,j);
                        if (new_distance > distance(1+i))
                            distance(1+i) = new_distance;
                        end
                    end
                end
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    p(1:nOfColumns)=distance(1+nOfRows+1:1+nOfRows+nOfColumns);
    for i=1:nOfRows
        if sum(starMatrix(i,:))==0
            r(i)=0;
        else
            r(i)=  originalDistMatrix(i,find(starMatrix(i,:))) - p(find(starMatrix(i,:)));
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


if sum(coveredColumns) == minDim
    % algorithm finished
    assignment = buildassignmentvector(starMatrix);
else
    % move to step 3
    [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v);
end

% compute value and remove invalid assignments
[assignment, value] = computeassignmentcost(assignment, originalDistMatrix, nOfRows);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function assignment = buildassignmentvector(starMatrix)

[maxValue, assignment] = max(starMatrix, [], 2);
assignment(maxValue == 0) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, value] = computeassignmentcost(assignment, distMatrix, nOfRows)

rowIndex   = find(assignment);
costVector = distMatrix(rowIndex + nOfRows * (assignment(rowIndex)-1));
value = sum(costVector);

% Step 2: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step2(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v)

% cover every column containing a starred zero
maxValue = max(starMatrix);
coveredColumns(maxValue == 1) = 1;

if sum(coveredColumns) == minDim
    % algorithm finished
    assignment = buildassignmentvector(starMatrix);
else
    % move to step 3
    [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v);
end

% Step 3: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v)

zerosFound = 1;
while zerosFound

    zerosFound = 0;
    for col = find(~coveredColumns)
        for row = find(~coveredRows')
            if distMatrix(row,col) == 0

                primeMatrix(row, col) = 1;
                starCol = find(starMatrix(row,:));
                if isempty(starCol)
                    % move to step 4
                    [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step4(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim,r,p,u,v);
                    return
                else
                    coveredRows(row)        = 1;
                    coveredColumns(starCol) = 0;
                    zerosFound              = 1;
                    break % go on in next column
                end
            end
        end
    end
end

% move to step 5
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step5(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v);

% Step 4: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step4(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim,r,p,u,v)

newStarMatrix          = starMatrix;
newStarMatrix(row,col) = 1;

starCol = col;
starRow = find(starMatrix(:, starCol));

while ~isempty(starRow)

    % unstar the starred zero
    newStarMatrix(starRow, starCol) = 0;

    % find primed zero in row
    primeRow = starRow;
    primeCol = find(primeMatrix(primeRow, :));

    % star the primed zero
    newStarMatrix(primeRow, primeCol) = 1;

    % find starred zero in column
    starCol = primeCol;
    starRow = find(starMatrix(:, starCol));

end
starMatrix = newStarMatrix;

primeMatrix(:) = 0;
coveredRows(:) = 0;

% move to step 2
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step2(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v);


% Step 5: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step5(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v)

% find smallest uncovered element
uncoveredRowsIndex    = find(~coveredRows');
uncoveredColumnsIndex = find(~coveredColumns);
[s, index1] = max(distMatrix(uncoveredRowsIndex,uncoveredColumnsIndex));
[s, index2] = max(s);
h = distMatrix(uncoveredRowsIndex(index1(index2)), uncoveredColumnsIndex(index2));

% add h to each covered row
index = find(coveredRows);
distMatrix(index, :) = distMatrix(index, :) + h;
u(index) = u(index)-h;
r(index) = r(index)-h;

% subtract h from each uncovered column
distMatrix(:, uncoveredColumnsIndex) = distMatrix(:, uncoveredColumnsIndex) - h;
v(uncoveredColumnsIndex) = v(uncoveredColumnsIndex)+h;
p(uncoveredColumnsIndex) = p(uncoveredColumnsIndex)+h;
% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows,r,p,u,v] = step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim,r,p,u,v);