Ported to javascript by Philippe Rivière, from the C++ implementation found at https://github.com/yongyanghz/LAPJV-algorithm-c
See Fil/lap-jv for details, and a simpler example.
<!DOCTYPE html>
<head>
<meta charset="utf-8">
<script src="https://d3js.org/d3.v4.min.js"></script>
<style>
body { margin:0;position:fixed;top:0;right:0;bottom:0;left:0; }
</style>
</head>
<body>
<script>
// Feel free to change or delete any of the code you see in this editor!
var svg = d3.select("body").append("svg")
.attr("width", 960)
.attr("height", 500)
const m = 40, n = m * m, w = Math.ceil(480/m);
// const data = d3.range(n).map(k => [Math.floor(100 * m * (2 + Math.cos( k ))/4)/100, Math.floor(100 * m * (2 + Math.sin(k))/4)/100]);
const data = d3.range(n).map(k => [m * Math.random(), m * Math.random()]);
// const data = d3.range(n).map(k => [m * (1 + Math.cos(k))/2, m * (1 + Math.sin(k))/2]);
// const data = d3.range(n).map(k => { const i = k % m, j = (k-i)/m; return [i + Math.random(),j + Math.random()]; });
data.map(d => d.color = d3.rgb(Math.random()*255, Math.random()*255, Math.random()*255));
svg.selectAll('circle')
.data(data)
.enter()
.append('circle')
.attr('r', 3);
svg.selectAll('circle')
.attr('cx', d => w * d[0])
.attr('cy', d => w * d[1])
.attr('fill', d => d.color);
createWorker('worker.js', (worker) => {
const start = performance.now();
worker.postMessage({
m: m,
data: data,
});
worker.onmessage = function (e) {
if (e.data.assign) {
console.log(n + ': ' + Math.ceil(performance.now() - start) + 'ms')
draw(e.data.assign);
}
}
});
function draw(res) {
res.col.map((c, k) => {
const i = k % m, j = (k-i)/m;
data[c].i = i;
data[c].j = j;
});
svg.selectAll('line')
.data(data)
.enter()
.append('line')
.attr('x1', d => w * d[0])
.attr('y1', d => w * d[1])
.attr('x2', d => w * d[0])
.attr('y2', d => w * d[1])
.attr('stroke', d => d.color)
.attr('opacity', 0.8)
.transition()
.duration(1500)
.attr('x2', d => w/2 + w * d.i)
.attr('y2', d => w/2 + w * d.j)
;
setTimeout(function() {
svg.selectAll('circle')
.transition()
.duration(500)
.attr('cx', d => w/2 + w * d.i)
.attr('cy', d => w/2 + w * d.j);
if (false) svg.selectAll('line')
.transition()
.attr('x1', d => w/2 + w * d.i)
.attr('y1', d => w/2 + w * d.j)
.duration(500)
}, 1500);
}
function createWorker(s, cb) {
const worker = new Worker(window.URL.createObjectURL(new Blob([`
importScripts('https://raw.githack.com/Fil/lap-jv/master/lap.js');
self.onmessage = function(e){
const m = e.data.m,
n = m * m,
data = e.data.data;
const costs = function(a, k) {
const d = data[a],
i = k % m,
j = (k-i)/m;
const dx = d[0] - i - 0.5, dy = d[1] - j - 0.5;
return dx * dx + dy * dy;
};
self.postMessage({
assign: lap(n, costs)
});
};`], {
type: "text/javascript"
})));
if (typeof cb == 'function') {
cb(worker);
}
}
</script>
</body>
/************************************************************************
*
* lap.js -- ported to javascript from
lap.cpp
version 1.0 - 4 September 1996
author: Roy Jonker @ MagicLogic Optimization Inc.
e-mail: roy_jonker@magiclogic.com
Code for Linear Assignment Problem, according to
"A Shortest Augmenting Path Algorithm for Dense and Sparse Linear
Assignment Problems," Computing 38, 325-340, 1987
by
R. Jonker and A. Volgenant, University of Amsterdam.
*
PORTED TO JAVASCRIPT 2017-01-02 by Philippe Riviere(fil@rezo.net)
CHANGED 2016-05-13 by Yang Yong(yangyongeducation@163.com) in column reduction part according to
matlab version of LAPJV algorithm(Copyright (c) 2010, Yi Cao All rights reserved)--
https://www.mathworks.com/matlabcentral/fileexchange/26836-lapjv-jonker-volgenant-algorithm-for-linear-assignment-problem-v3-0:
*
*************************************************************************/
/* This function is the jv shortest augmenting path algorithm to solve the assignment problem */
function lap(dim, assigncost) {
// input:
// dim - problem size
// assigncost - cost matrix
// output:
// rowsol - column assigned to row in solution
// colsol - row assigned to column in solution
// u - dual variables, row reduction numbers
// v - dual variables, column reduction numbers
const sum = assigncost.map(d => d.reduce((a, b) => a + b, 0)).reduce((a, b) => a + b, 0);
const BIG = 10000 * (sum / dim);
const epsilon = (sum / dim) / 10000;
const rowsol = new Int32Array(dim),
colsol = new Int32Array(dim),
u = new Float64Array(dim),
v = new Float64Array(dim);
let unassignedfound;
/* row */
let i, imin, numfree = 0,
prvnumfree, f, i0, k, freerow; // *pred, *free
/* col */
let j, j1, j2, endofpath, last, low, up; // *collist, *matches
/* cost */
let min, h, umin, usubmin, v2; // *d
const free = new Int32Array(dim); // list of unassigned rows.
const collist = new Int32Array(dim); // list of columns to be scanned in various ways.
const matches = new Int32Array(dim); // counts how many times a row could be assigned.
const d = new Float64Array(dim); // 'cost-distance' in augmenting path calculation.
const pred = new Int32Array(dim); // row-predecessor of column in augmenting/alternating path.
// init how many times a row will be assigned in the column reduction.
for (i = 0; i < dim; i++)
matches[i] = 0;
// COLUMN REDUCTION
for (j = dim; j--;) // reverse order gives better results.
{
// find minimum cost over rows.
min = assigncost[0][j];
imin = 0;
for (i = 1; i < dim; i++)
if (assigncost[i][j] < min) {
min = assigncost[i][j];
imin = i;
}
v[j] = min;
if (++matches[imin] == 1) {
// init assignment if minimum row assigned for first time.
rowsol[imin] = j;
colsol[j] = imin;
} else if (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.
}
// REDUCTION TRANSFER
for (i = 0; i < dim; i++)
if (matches[i] == 0) // fill list of unassigned 'free' rows.
free[numfree++] = i;
else
if (matches[i] == 1) // transfer reduction from rows that are assigned once.
{
j1 = rowsol[i];
min = BIG;
for (j = 0; j < dim; j++)
if (j != j1)
if (assigncost[i][j] - v[j] < min + epsilon)
min = assigncost[i][j] - v[j];
v[j1] = v[j1] - min;
}
// AUGMENTING ROW REDUCTION
let loopcnt = 0; // do-loop to be done twice.
do {
loopcnt++;
// scan all free rows.
// in some cases, a free row may be replaced with another one to be scanned next.
k = 0;
prvnumfree = numfree;
numfree = 0; // start list of rows still free after augmenting row reduction.
while (k < prvnumfree) {
i = free[k];
k++;
// find minimum and second minimum reduced cost over columns.
umin = assigncost[i][0] - v[0];
j1 = 0;
usubmin = BIG;
for (j = 1; j < dim; j++) {
h = assigncost[i][j] - v[j];
if (h < usubmin)
if (h >= umin) {
usubmin = h;
j2 = j;
} else {
usubmin = umin;
umin = h;
j2 = j1;
j1 = j;
}
}
i0 = colsol[j1];
if (umin < usubmin + epsilon)
// change the reduction of the minimum column to increase the minimum
// reduced cost in the row to the subminimum.
v[j1] = v[j1] - (usubmin + epsilon - umin);
else // minimum and subminimum equal.
if (i0 > -1) // minimum column j1 is assigned.
{
// swap columns j1 and j2, as j2 may be unassigned.
j1 = j2;
i0 = colsol[j2];
}
// (re-)assign i to j1, possibly de-assigning an i0.
rowsol[i] = j1;
colsol[j1] = i;
if (i0 > -1) // minimum column j1 assigned earlier.
if (umin < usubmin)
// put in current k, and go back to that k.
// continue augmenting path i - j1 with i0.
free[--k] = i0;
else
// no further augmenting reduction possible.
// store i0 in list of free rows for next phase.
free[numfree++] = i0;
}
}
while (loopcnt < 2); // repeat once.
// AUGMENT SOLUTION for each free row.
for (f = 0; f < numfree; f++) {
freerow = free[f]; // start row of augmenting path.
// Dijkstra shortest path algorithm.
// runs until unassigned column added to shortest path tree.
for (j = dim; j--;) {
d[j] = assigncost[freerow][j] - v[j];
pred[j] = freerow;
collist[j] = j; // init column list.
}
low = 0; // columns in 0..low-1 are ready, now none.
up = 0; // columns in low..up-1 are to be scanned for current minimum, now none.
// columns in up..dim-1 are to be considered later to find new minimum,
// at this stage the list simply contains all columns
unassignedfound = false;
do {
if (up == low) // no more columns to be scanned for current minimum.
{
last = low - 1;
// scan columns for up..dim-1 to find all indices for which new minimum occurs.
// store these indices between low..up-1 (increasing up).
min = d[collist[up++]];
for (k = up; k < dim; k++) {
j = collist[k];
h = d[j];
if (h <= min) {
if (h < min) // new minimum.
{
up = low; // restart list at index low.
min = h;
}
// new index with same minimum, put on undex up, and extend list.
collist[k] = collist[up];
collist[up++] = j;
}
}
// check if any of the minimum columns happens to be unassigned.
// if so, we have an augmenting path right away.
for (k = low; k < up; k++)
if (colsol[collist[k]] < 0) {
endofpath = collist[k];
unassignedfound = true;
break;
}
}
if (!unassignedfound) {
// update 'distances' between freerow and all unscanned columns, via next scanned column.
j1 = collist[low];
low++;
i = colsol[j1];
h = assigncost[i][j1] - v[j1] - min;
for (k = up; k < dim; k++) {
j = collist[k];
v2 = assigncost[i][j] - v[j] - h;
if (v2 < d[j]) {
pred[j] = i;
if (v2 == min) // new column found at same minimum value
if (colsol[j] < 0) {
// if unassigned, shortest augmenting path is complete.
endofpath = j;
unassignedfound = true;
break;
}
// else add to list to be scanned right away.
else {
collist[k] = collist[up];
collist[up++] = j;
}
d[j] = v2;
}
}
}
}
while (!unassignedfound);
// update column prices.
for (k = last + 1; k--;) {
j1 = collist[k];
v[j1] = v[j1] + d[j1] - min;
}
// reset row and column assignments along the alternating path.
do {
i = pred[endofpath];
colsol[endofpath] = i;
j1 = endofpath;
endofpath = rowsol[i];
rowsol[i] = j1;
}
while (i != freerow);
}
// calculate optimal cost.
let lapcost = 0;
for (i = dim; i--;) {
j = rowsol[i];
u[i] = assigncost[i][j] - v[j];
lapcost = lapcost + assigncost[i][j];
}
return {
cost: lapcost,
row: rowsol,
col: colsol,
u: u,
v: v
};
}