I’ll work through some theory and then give some JavaScript code that meets 1–4 and 6, a modified version of 5, and additionally
- is robust to numerical issues
- has no external dependencies.
We can express 1–4 as a linear program. Your sample instance, for example, would translate as
minimize x4 - x1
subject to
(known minimum distances)
a: x3 - x1 >= 200
b: x4 - x1 >= 900
(fixed order)
c: x2 - x1 >= 0
d: x3 - x2 >= 0
e: x4 - x3 >= 0
x1, x2, x3, x4 unbounded
Every linear program has a dual program with the same optimum objective value (under conditions that are met here). The dual of the program above is below:
maximize 200 a + 900 b
subject to
x1: -a - b - c = -1
x2: c - d = 0
x3: a + d - e = 0
x4: b + e = 1
a, b, c, d, e >= 0
This is a max-cost flow problem with no capacity constraints. In consequence, some path is an optimal flow, so we can understand criteria 1–4 as a longest path problem.
Once we have the optimal objective, 900 here, we can shift criterion 4 to a constraint
x4 - x1 <= 900
and think about how to make the most pleasing layout. As written, we can achieve criterion 5 by solving a quadratic program:
minimize (x2 - x1)**2 + (x3 - x2)**2 + (x4 - x3)**2
subject to
(most efficient layout)
x4 - x1 <= 900
(known minimum distances)
x3 - x1 >= 200
x4 - x1 >= 900
(fixed order)
x2 - x1 >= 0
x3 - x2 >= 0
x4 - x3 >= 0
x1, x2, x3, x4 unbounded
The objective that I’ve written is not the standard deviation, but it doesn’t matter. Minimizing the standard deviation is equivalent to minimizing the variance. In turn, the variance of a random variable X is E[X**2] - E[X]**2, and E[X]**2 is a constant here (E[X] = 900/4 = 150 in the example). Finally, maximizing E[X**2] is equivalent to maximizing the sum of squares of the equally likely outcomes.
I would propose a different criterion however.
5′. Maximize (in lexicographic order) the list of gap lengths sorted small to large.
What this means is that we try as hard as possible to avoid small gaps. In an example from Roman Svistunov, suppose that we have constraints
x3 - x1 >= 2
x4 - x2 >= 5.5
x5 - x3 >= 8
The solution with minimum standard deviation sets
x1 = 0
x2 = 0.75
x3 = 2
x4 = 6.25
x5 = 10
whereas the solution that optimizes the replacement objective sets
x1 = 0
x2 = 1
x3 = 2
x4 = 6.5
x5 = 10
We’ve made x2 more centered at the cost of making x4 less centered. Since x2 is closer to its neighbors than x4, this seems like a reasonable trade to me, but you obviously have the final say.
To optimize 5′, we proceed as follows. Start by formulating a linear program like
maximize z
subject to
(most efficient layout)
x5 - x1 <= 10
(known minimum distances)
x3 - x1 >= 2
x4 - x2 >= 5.5
x5 - x3 >= 8
(fixed order)
x2 - x1 >= z
x3 - x2 >= z
x4 - x3 >= z
x5 - x4 >= z
x1, x2, x3, x4, x5 unbounded
We solve for the optimum, z = 1 here, and also record which gaps prevented it from being larger. We replace those z variables and solve again, repeating until there is no occurrence of z.
maximize z
subject to
(most efficient layout)
x5 - x1 <= 10
(known minimum distances)
x3 - x1 >= 2
x4 - x2 >= 5.5
x5 - x3 >= 8
(fixed order)
x2 - x1 = 1
x3 - x2 = 1
x4 - x3 >= z
x5 - x4 >= z
x1, x2, x3, x4, x5 unbounded
The optimum is z = 3.5.
maximize z
subject to
(most efficient layout)
x5 - x1 <= 10
(known minimum distances)
x3 - x1 >= 2
x4 - x2 >= 5.5
x5 - x3 >= 8
(fixed order)
x2 - x1 = 1
x3 - x2 = 1
x4 - x3 >= z
x5 - x4 = 3.5
x1, x2, x3, x4, x5 unbounded
The optimum is z = 4.5, and that’s the last iteration.
As observed earlier, we can solve these linear programs using an extended longest path algorithm. In JavaScript:
// Dual numbers represent linear functions of time.
function Dual(p, v) {
return { position: p, velocity: v };
}
// Adds dual numbers x and y.
function dualPlus(x, y) {
return Dual(x.position + y.position, x.velocity + y.velocity);
}
const epsilon = Math.sqrt(Number.EPSILON);
// Compares dual numbers x and y by their position at a time infinitesimally
// after now.
function dualLessThan(x, y) {
let d = x.position - y.position;
return d < -epsilon || (Math.abs(d) <= epsilon && x.velocity < y.velocity);
}
// Tracks delta, the length of time for which every pair of arguments to
// .dualLessThan() maintains the present relation.
function DeltaTracker() {
return {
delta: Infinity,
dualLessThan: function (x, y) {
let lessThan = dualLessThan(x, y);
if (lessThan) {
[x, y] = [y, x];
}
if (x.velocity < y.velocity) {
this.delta = Math.min(
this.delta,
(x.position - y.position) / (y.velocity - x.velocity)
);
}
return lessThan;
},
};
}
// Converts the adjacency matrix to an adjacency list representation.
function graphFromMatrix(n, matrix) {
let graph = [];
for (let j = 0; j < n; j++) {
graph.push([]);
for (let i = 0; i < j; i++) {
if (matrix[i][j] > 0) {
graph[j].push({ i: i, length: Dual(matrix[i][j], 0) });
}
}
}
return graph;
}
// Essentially the usual algorithm for longest path, but tracks delta.
function longestPathTable(graph, gaps) {
let tracker = DeltaTracker();
let maximum = Dual(0, 0);
let table = [];
for (let j = 0; j < graph.length; j++) {
let argument = null;
if (j > 0) {
maximum = dualPlus(maximum, gaps[j - 1]);
}
for (let edge of graph[j]) {
let x = dualPlus(table[edge.i].maximum, edge.length);
if (tracker.dualLessThan(maximum, x)) {
argument = edge.i;
maximum = x;
}
}
table.push({ argument: argument, maximum: maximum });
}
return [tracker.delta, table];
}
// Essentially the usual algorithm for decoding the longest path from the
// dynamic program table.
function freezeCriticalGaps(table, graph, gaps) {
let j = table.length - 1;
while (j > 0) {
let critical = true;
let argument = table[j].argument;
if (argument !== null) {
j = argument;
} else {
j--;
gaps[j].velocity = 0;
}
}
}
// Changes the time from now to now + delta.
function advanceGaps(gaps, delta) {
for (let i = 0; i < gaps.length; i++) {
gaps[i].position += gaps[i].velocity * delta;
}
}
// Extracts the final result from the dynamic program table.
function positionsFromTable(table) {
let positions = [];
for (let entry of table) {
positions.push(entry.maximum.position);
}
return positions;
}
// Entry point for the layout algorithm.
function layOut(n, matrix) {
let graph = graphFromMatrix(n, matrix);
let gaps = [];
for (let j = 1; j < n; j++) {
gaps.push(Dual(0, 1));
}
while (true) {
let [delta, table] = longestPathTable(graph, gaps);
if (delta == Infinity) {
return positionsFromTable(table);
}
if (table[n - 1].maximum.velocity > 0) {
freezeCriticalGaps(table, graph, gaps);
} else {
advanceGaps(gaps, delta);
}
}
}
// Various test cases.
console.log(
layOut(4, [
[0, 200, 0, 900],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
])
);
console.log(
layOut(5, [
[0, 0, 2, 0, 0],
[0, 0, 0, 5.5, 0],
[0, 0, 0, 0, 8],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
])
);
console.log(
layOut(4, [
[0, 0, 200, 0],
[0, 0, 0, 150],
[0, 0, 0, 0],
[0, 0, 0, 0],
])
);
console.log(
layOut(4, [
[0, 0, 200, 0],
[0, 0, 0, 200],
[0, 0, 0, 0],
[0, 0, 0, 0],
])
);
console.log(
layOut(4, [
[0, 0, 200, 0],
[0, 0, 0, 300],
[0, 0, 0, 0],
[0, 0, 0, 0],
])
);
