function [rowsol,cost,v,u,costMat] = lapjv(costMat,resolution)
% LAPJV Jonker-Volgenant Algorithm for Linear Assignment Problem.
%
% [ROWSOL,COST,v,u,rMat] = LAPJV(COSTMAT, resolution) returns the optimal column indices,
% ROWSOL, assigned to row in solution, and the minimum COST based on the
% assignment problem represented by the COSTMAT, where the (i,j)th element
% represents the cost to assign the jth job to the ith worker.
% The second optional input can be used to define data resolution to
% accelerate speed.
% Other output arguments are:
% v: dual variables, column reduction numbers.
% u: dual variables, row reduction numbers.
% rMat: the reduced cost matrix.
%
% For a rectangular (nonsquare) costMat, rowsol is the index vector of the
% larger dimension assigned to the smaller dimension.
%
% [ROWSOL,COST,v,u,rMat] = LAPJV(COSTMAT,resolution) accepts the second
% input argument as the minimum resolution to differentiate costs between
% assignments. The default is eps.
%
% Known problems: The original algorithm was developed for integer costs.
% When it is used for real (floating point) costs, sometime the algorithm
% will take an extreamly long time. In this case, using a reasonable large
% resolution as the second arguments can significantly increase the
% solution speed.
%
% See also munkres, Hungarian
% version 1.0 by Yi Cao at Cranfield University on 3rd March 2010
% version 1.1 by Yi Cao at Cranfield University on 19th July 2010
% version 1.2 by Yi Cao at Cranfield University on 22nd July 2010
% version 2.0 by Yi Cao at Cranfield University on 28th July 2010
% version 2.1 by Yi Cao at Cranfield University on 13th August 2010
% version 2.2 by Yi Cao at Cranfield University on 17th August 2010
% version 3.0 by Yi Cao at Cranfield University on 10th April 2013
% This Matlab version is developed based on the orginal C++ version coded
% by Roy Jonker @ MagicLogic Optimization Inc on 4 September 1996.
% Reference:
% R. Jonker and A. Volgenant, "A shortest augmenting path algorithm for
% dense and spare linear assignment problems", Computing, Vol. 38, pp.
% 325-340, 1987.
%
% Examples
% Example 1: a 5 x 5 example
%{
[rowsol,cost] = lapjv(magic(5));
disp(rowsol); % 3 2 1 5 4
disp(cost); %15
%}
% Example 2: 1000 x 1000 random data
%{
n=1000;
A=randn(n)./rand(n);
tic
[a,b]=lapjv(A);
toc % about 0.5 seconds
%}
% Example 3: nonsquare test
%{
n=100;
A=1./randn(n);
tic
[a,b]=lapjv(A);
toc % about 0.2 sec
A1=[A zeros(n,1)+max(max(A))];
tic
[a1,b1]=lapjv(A1);
toc % about 0.01 sec. The nonsquare one can be done faster!
%check results
disp(norm(a-a1))
disp(b-b)
%}
if nargin<2
maxcost=min(1e16,max(max(costMat)));
resolution=eps(maxcost);
end
% Prepare working data
[rdim,cdim] = size(costMat);
M=min(min(costMat));
if rdim>cdim
costMat = costMat';
[rdim,cdim] = size(costMat);
swapf=true;
else
swapf=false;
end
dim=cdim;
costMat = [costMat;2*M+zeros(cdim-rdim,cdim)];
costMat(costMat~=costMat)=Inf;
maxcost=max(costMat(costMat<Inf))*dim+1;
if isempty(maxcost)
maxcost = Inf;
end
costMat(costMat==Inf)=maxcost;
% free = zeros(dim,1); % list of unssigned rows
% colist = 1:dim; % list of columns to be scaed in various ways
% d = zeros(1,dim); % 'cost-distance' in augmenting path calculation.
% pred = zeros(dim,1); % row-predecessor of column in augumenting/alternating path.
v = zeros(1,dim); % dual variables, column reduction numbers.
rowsol = zeros(1,dim)-1; % column assigned to row in solution
colsol = zeros(dim,1)-1; % row assigned to column in solution
numfree=0;
free = zeros(dim,1); % list of unssigned rows
matches = zeros(dim,1); % counts how many times a row could be assigned.
% The Initilization Phase
% column reduction
for j=dim:-1:1 % reverse order gives better results
% find minimum cost over rows
[v(j), imin] = min(costMat(:,j));
if ~matches(imin)
% init assignement if minimum row assigned for first time
rowsol(imin)=j;
colsol(j)=imin;
elseif v(j)<v(rowsol(imin))
j1=rowsol(imin);
rowsol(imin)=j;
colsol(j)=imin;
colsol(j1)=-1;
else
colsol(j)=-1; % row already assigned, column not assigned.
end
matches(imin)=matches(imin)+1;
end
% Reduction transfer from unassigned to assigned rows
for i=1:dim
if ~matches(i) % fill list of unaasigned 'free' rows.
numfree=numfree+1;
free(numfree)=i;
else
if matches(i) == 1 % transfer reduction from rows that are assigned once.
j1 = rowsol(i);
x = costMat(i,:)-v;
x(j1) = maxcost;
v(j1) = v(j1) - min(x);
end
end
end
% Augmenting reduction of unassigned rows
loopcnt = 0;
while loopcnt < 2
loopcnt = loopcnt + 1;
% scan all free rows
% in some cases, a free row may be replaced with another one to be scaed next
k = 0;
prvnumfree = numfree;
numfree = 0; % start list of rows still free after augmenting row reduction.
while k < prvnumfree
k = k+1;
i = free(k);
% find minimum and second minimum reduced cost over columns
x = costMat(i,:) - v;
[umin, j1] = min(x);
x(j1) = maxcost;
[usubmin, j2] = min(x);
i0 = colsol(j1);
if usubmin - umin > resolution
% change the reduction of the minmum column to increase the
% minimum reduced cost in the row to the subminimum.
v(j1) = v(j1) - (usubmin - umin);
else % minimum and subminimum equal.
if i0 > 0 % minimum column j1 is assigned.
% swap columns j1 and j2, as j2 may be unassigned.
j1 = j2;
i0 = colsol(j2);
end
end
% reassign i to j1, possibly de-assigning an i0.
rowsol(i) = j1;
colsol(j1) = i;
if i0 > 0 % ,inimum column j1 assigned easier
if usubmin - umin > resolution
% put in current k, and go back to that k.
% continue augmenting path i - j1 with i0.
free(k)=i0;
k=k-1;
else
% no further augmenting reduction possible
% store i0 in list of free rows for next phase.
numfree = numfree + 1;
free(numfree) = i0;
end
end
end
end
% Augmentation Phase
% augment solution for each free rows
for f=1:numfree
freerow = free(f); % start row of augmenting path
% Dijkstra shortest path algorithm.
% runs until unassigned column added to shortest path tree.
d = costMat(freerow,:) - v;
pred = freerow(1,ones(1,dim));
collist = 1:dim;
low = 1; % columns in 1...low-1 are ready, now none.
up = 1; % columns in low...up-1 are to be scaed for current minimum, now none.
% columns in up+1...dim are to be considered later to find new minimum,
% at this stage the list simply contains all columns.
unassignedfound = false;
while ~unassignedfound
if up == low % no more columns to be scaned for current minimum.
last = low-1;
% scan columns for up...dim to find all indices for which new minimum occurs.
% store these indices between low+1...up (increasing up).
minh = d(collist(up));
up = up + 1;
for k=up:dim
j = collist(k);
h = d(j);
if h<=minh
if h<minh
up = low;
minh = h;
end
% new index with same minimum, put on index up, and extend list.
collist(k) = collist(up);
collist(up) = j;
up = up +1;
end
end
% check if any of the minimum columns happens to be unassigned.
% if so, we have an augmenting path right away.
for k=low:up-1
if colsol(collist(k)) < 0
endofpath = collist(k);
unassignedfound = true;
break
end
end
end
if ~unassignedfound
% update 'distances' between freerow and all unscanned columns,
% via next scanned column.
j1 = collist(low);
low=low+1;
i = colsol(j1); %line 215
x = costMat(i,:)-v;
h = x(j1) - minh;
xh = x-h;
k=up:dim;
j=collist(k);
vf0 = xh<d;
vf = vf0(j);
vj = j(vf);
vk = k(vf);
pred(vj)=i;
v2 = xh(vj);
d(vj)=v2;
vf = v2 == minh; % new column found at same minimum value
j2 = vj(vf);
k2 = vk(vf);
cf = colsol(j2)<0;
if any(cf) % unassigned, shortest augmenting path is complete.
i2 = find(cf,1);
endofpath = j2(i2);
unassignedfound = true;
else
i2 = numel(cf)+1;
end
% add to list to be scaned right away
for k=1:i2-1
collist(k2(k)) = collist(up);
collist(up) = j2(k);
up = up + 1;
end
end
end
% update column prices
j1=collist(1:last+1);
v(j1) = v(j1) + d(j1) - minh;
% reset row and column assignments along the alternating path
while 1
i=pred(endofpath);
colsol(endofpath)=i;
j1=endofpath;
endofpath=rowsol(i);
rowsol(i)=j1;
if (i==freerow)
break
end
end
end
rowsol = rowsol(1:rdim);
u=diag(costMat(:,rowsol))-v(rowsol)';
u=u(1:rdim);
v=v(1:cdim);
cost = sum(u)+sum(v(rowsol));
costMat=costMat(1:rdim,1:cdim);
costMat = costMat - u(:,ones(1,cdim)) - v(ones(rdim,1),:);
if swapf
costMat = costMat';
t=u';
u=v';
v=t;
end
if cost>maxcost
cost=Inf;
end