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.
license: mit | |
border: no |
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 | |
}; | |
} |