Get the position of points on a line segment with constraints-CodePudding

The
I have implemented the following which satisfies constraints #1~#4. Can some please help with constraint #5? It can be partially met (like b in the example) or fully met (like c in the example). We can use standard deviation to measure the quality of the solution.

// Know minimum distance is stored in a matrix // e.g. // m = [ // [0, 100, 0, 900], // [0, 0, 0, 0], // [0, 0, 0, 0], // [0, 0, 0, 0], // ] let m = Array<Array<number>>(); let min_distance = (i: number, j: number) => { return m[i][j]; } let final_distance = (i: number, j: number): number => { if (j === i) { return 0; } else if(j - i === 1) { let temp = []; for(let k = 0; k <= i; k ) { temp.push(min_distance(k, j) - final_distance(k, i)); } return Math.max(...temp) } else { let sum = 0; for(let l = i; l < j; l ) { sum = final_distance(l, l 1); } return sum; } } let final_pos = (i: number): number => { return final_distance(0, i); }

CodePudding user response：

Ooh, neat problem.

We can express 1–4 as a linear program. (This looks like a step backward, but I promise I’ll pay it off in the end with code.) 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.

I bring up linear programming for two reasons. First, if you can find a convex optimization library for your preferred programming environment, then you can use it to solve this linear program, add a new constraint on the primary objective like x4 - x1 <= 900, and then solve the new program as a quadratic program that minimizes the sum of the squares of the gaps, which is equivalent to minimizing their standard deviation.

Second, I want to talk about complementary slackness. This means here that an arc (i.e., a minimum distance constraint) is on some longest path if and only if, in all optimal solutions, the distance between its endpoints is the minimum.

Let me propose a different interpretation of criterion 5: that we maximize the minimum gap (breaking ties by the second minimum, breaking ties by the third, etc.; the tie breaking is really important since otherwise a single small gap can mess up even spacing elsewhere). I would argue that it’s better to use a maximally outlier-sensitive objective for a problem like this. Then we can use complementary slackness to design a primal-dual algorithm:

Use the given minimum distances to form a directed graph. Add “gap” arcs from one point to the next, starting at length zero.

Compute a longest path in the current graph.

If there are gap arcs that have not appeared in some longest path yet, lengthen each of them by some amount ε > 0 and go to 2.

In (somewhat hurried) Python:

import collections EPSILON = 1 # A graph maps each vertex to a list of pairs (outbound neighbor, length). def make_graph(min_distances): graph = collections.defaultdict(list) for u, v, l in min_distances: graph[u].append((v, l)) return graph # Returns a longest path from first to last. def longest_path(graph, first, last): # Dictionary from each vertex u reachable from vertex 1 to the length of the # longest path from u. table = {} stack = [first] while stack: u = stack.pop() if u not in table: # Mark the computation of u as in progress. None will cause an error # if used below (giving us cycle detection). table[u] = None # Put u back on the stack so that we revisit it in post order. stack.append(u) stack.extend(v for (v, l) in graph[u]) elif table[u] is None: # All of the inputs are ready. table[u] = max((l table[v] for (v, l) in graph[u]), default=0) # Read off the longest path. path = [first] u = first while u != last: _, u = max((l table[v], v) for (v, l) in graph[u]) path.append(u) return path, table def layout(min_distances, n): # If gap_lengths[u] exists, then the value is the length of the gap after u. gap_lengths = {} # Minimum length of gaps whose length is to be determined. min_gap = 0 # Loop until all gap lengths are known. while not all(u in gap_lengths for u in range(1, n)): graph = make_graph(min_distances) for u in range(1, n): graph[u].append((u 1, gap_lengths.get(u, min_gap))) path, table = longest_path(graph, 1, n) for i in range(1, len(path)): u = path[i - 1] v = path[i] if u == v - 1: if u not in gap_lengths and table[u] - table[v] == min_gap: gap_lengths[u] = min_gap min_gap = EPSILON return [gap_lengths[u] for u in range(1, n)] if __name__ == "__main__": print(layout([(1, 3, 200), (1, 4, 900)], 4)) print(layout([(1, 3, 200), (2, 4, 300)], 4))

This algorithm is a little boring on the given example; we lengthen each of the gap arcs to 301 and then stop. It’s also slow if we make ε small but imprecise if we make it large. We can fix the problem by varying ε (which I’ll rename δ since it may be large now), however, at the cost of making the code more complicated.

Each iteration of the algorithm implemented below is O(n m) where m is the number of constraints, with a reasonable constant. The number of iterations should be O(n) in practice, but I think it can be larger on carefully constructed instances. Hopefully this is fast enough for your purposes.

#include <algorithm> #include <iostream> #include <limits> #include <optional> #include <utility> #include <vector> namespace { //// // Problem representation. // Represents the index of a point or (later) a vertex. using Index = std::size_t; // Constrains points left and right to be at least length apart. struct Constraint { Index left; Index right; double length; }; // Represents a problem instance. struct Instance { Index n; std::vector<Constraint> constraints; }; //// // Represents a dual number. We use dual numbers when we need to track both the // value and the derivative of a quantity. struct DualNumber { double value; double derivative; }; DualNumber operator (DualNumber x, DualNumber y) { return {x.value y.value, x.derivative y.derivative}; } DualNumber operator-(DualNumber x, DualNumber y) { return {x.value - y.value, x.derivative - y.derivative}; } std::ostream &operator<<(std::ostream &o, DualNumber x) { return o << x.value << " " << x.derivative << "*epsilon"; } // Defined so that std::max favors the longest path whose length is growing // fastest. bool operator<(DualNumber x, DualNumber y) { return x.value < y.value || (x.value == y.value && x.derivative < y.derivative); } //// // Directed graph representation. As noted in my answer, this problem is closely // related to the longest path problem. struct Arc { Index left; // Index right is given by the index of the vertex in which this arc appears. DualNumber length; }; struct Vertex { // Arcs for which this vertex is the right endpoint. std::vector<Arc> arcs; DualNumber greatest_path_length = {0, 0}; std::optional<Index> next_arc_on_longest_path = std::nullopt; }; using Graph = std::vector<Vertex>; Graph GraphFromInstance(const Instance &problem) { Graph graph(problem.n); // Insert the gap arcs. for (Index u = 1; u < graph.size(); u ) { graph.at(u).arcs.push_back({u - 1, {0, 1}}); } // Insert the constraint arcs. for (Constraint c : problem.constraints) { graph.at(c.right).arcs.push_back({c.left, {c.length, 0}}); } // Compute initial longest paths. for (Index u = 0; u < graph.size(); u ) { Vertex &v = graph.at(u); for (Index i = 0; i < v.arcs.size(); i ) { Arc arc = v.arcs.at(i); DualNumber length = graph.at(arc.left).greatest_path_length arc.length; if (v.greatest_path_length < length) { v.next_arc_on_longest_path = i; v.greatest_path_length = length; } } } return graph; } void RecomputePathLengths(Graph &graph) { for (Index u = 1; u < graph.size(); u ) { Vertex &v = graph.at(u); Arc arc = v.arcs.at(v.next_arc_on_longest_path.value()); v.greatest_path_length = graph.at(arc.left).greatest_path_length arc.length; } } void StopGrowthOnLongestPath(Graph &graph) { for (Index u = graph.size() - 1; u != 0;) { Vertex &v = graph.at(u); Arc &arc = v.arcs.at(v.next_arc_on_longest_path.value()); arc.length.derivative = 0; u = arc.left; } RecomputePathLengths(graph); } std::vector<double> LayOut(const Instance &problem) { Graph graph = GraphFromInstance(problem); for (;;) { StopGrowthOnLongestPath(graph); // Compute delta, the time to the next gap becoming tight. double delta = std::numeric_limits<double>::infinity(); std::optional<std::pair<Index, Index>> became_tight = std::nullopt; for (Index u = 1; u < graph.size(); u ) { const Vertex &v = graph.at(u); for (Index i = 0; i < v.arcs.size(); i ) { Arc arc = v.arcs.at(i); DualNumber greatest_via_arc = graph.at(arc.left).greatest_path_length arc.length; DualNumber slack = v.greatest_path_length - greatest_via_arc; if (slack.derivative < 0) { double x = slack.value / -slack.derivative; if (x < delta) { became_tight = {u, i}; delta = x; } } } } if (!became_tight.has_value()) { break; } // Grow gaps. for (Index u = 1; u < graph.size(); u ) { for (Arc &arc : graph.at(u).arcs) { arc.length.value = arc.length.derivative * delta; } } // Record the gap that became tight. auto [u, i] = became_tight.value(); graph.at(u).next_arc_on_longest_path = i; } // Extract the final gaps. std::vector<double> gaps(graph.size() - 1); for (Index u = 1; u < graph.size(); u ) { gaps.at(u - 1) = graph.at(u).arcs.front().length.value; } return gaps; } std::ostream &operator<<(std::ostream &o, const std::vector<double> &v) { o << "["; bool comma = false; for (const double &x : v) { if (comma) { o << ", "; } comma = true; o << x; } return o << "]"; } } // namespace int main() { std::cout << LayOut({4, {{0, 2, 200}, {0, 3, 900}}}) << "\n"; std::cout << LayOut({4, {{0, 2, 200}, {1, 3, 150}}}) << "\n"; std::cout << LayOut({4, {{0, 2, 200}, {1, 3, 200}}}) << "\n"; std::cout << LayOut({4, {{0, 2, 200}, {1, 3, 300}}}) << "\n"; }

CodePudding user response：

First, we should find the smallest possible length of the segment. This can be done with the greedy algorithm: place the first point at 0 (without loss of generality) and add the following points one by one with a coordinate as small as possible without violating distances to all previous points. We would fix the last point at its position and would not move it any further. I would suppose that correctness is obvious.

Now the tricky part with optimizing the standard deviation of neighboring distances. We would introduce n-1 variables - the distances, I would name them Di. Trying to minimize their standard deviation is the same as minimizing variance.

$\bigsum\limits_{i=1}^n(D_i-E(D))^2=\bigsum\limits_{i=1}^nD_i^2-2E(D)\bigsum\limits_{i=1}^nD_i nE(D)^2=\bigsum\limits_{i=1}^nD_i^2-nE(D)^2$

And because E(D) is constant as we fixed the first and last points, minimizing standard deviation is the same as minimizing the sum of differences.

Constraints on distances are all linear as whenever you have m_d(Pi,Pj) you would add constraint (sum Dk for k=i to j) >= m_d(Pi,Pj). To guarantee that sum of distances is what we need it to be we also would add contracts sum Dk is both less or equal and greater or equal to the minimum distance calculated in step one. Therefore this problem is a quadratic optimization problem, that can be solved with libraries.

Note that the coordinates for the points then would be easily restored from distances as prefix sums.

Even though for 30 points probably any algorithm of quadratic programming would be fast enough if you want to invest time into further optimization if you intend to increase the number of points, as this problem can actually be solved exactly in polynomial time.

The first and most important optimization. Because the target function is the sum of squares, the matrix of target quadratic form is positive definite (as it is the identity matrix), and therefore there exists a polynomial algorithm, as have been shown in this paper.

The last optimization I would like to mention might not be relevant, but it is very good on most of the possible m_d matrices. To reduce the quadratic programming problem size, we can find optimal positions for some points before quadratic programming.

First, we would make reductions in matrix to ensure that for any i<j<k it is true that m_d(Pi,Pj) m_d(Pj,Pk)>=m_d(Pi,Pk) with updating the right side. Note that with this computation m_d(P1,Pn) is also the minimum possible length of the segment. After this, we can say that for any i, if m_d(P1,Pi) m_d(Pi,Pn)=m_d(P1,Pn) is true there is only one possible position for the point Pi, and, this allows us to remove one of the variables as now distance from P(i-1) to Pi and from Pi to P(i 1) can be replaced with one variable - distance from P(i-1) to P(i 1): the first distance is just Pi position - P(i-1) position with can be derived by induction, and the P(i 1) position is P(i-1) position this new distance.

Edit

As noted in @david-eisenstat there exists an interpretation of the 5th point as maximizing minimum gap, but there is a simpler and asymptotically faster solution, namely, answer binary search.

We know that the minimum gap is at least 0 and at most the segment length / n. We can maximize it with binary search, with a checker ensuring that every gap is at least our value. If findMinLength is a function that solves 1-4, the following is the solution that optimizes this metric (pseudocode):

function findOptimal(Matrix md, int n) { minLength = findMinLength(md, n) //The answer to minimize segment length minGapL = 0 //Lower bound in maximized minimum gap length minGapR = ceiling(minLength / n) //Upper bound in maximized minimum gap while minGapL < minGapR { minGapTest = (minGapL minGapR) / 2 mdCopy = m_d.copy() //This loop ensures that every gap is at least minimum gap tested for (int i = 1; i < n; i) { mdCopy[i - 1][i] = max(mdCopy[i - 1][i], miGapTest) } int testLength = findMinLength(mdCopy, n) if (testLength == minLength) { //if it is possible to get tested minimum gap update lower bound minGapL = minGapTest } else { //otherwise update lower bound minGapR = minGapTest - 1; } } return (minLength, minGapL); }

With this function returning both length and maximized minGap one can reconstruct the answer with the matrix modification as in the loop and then solve just problem for 1-4. This solution is polynomial (the 1-4 problem can be solved in linear time as shown above, and binary search is logarithmic of the answer, and with the answer being exponential of the input size it is also linear of the input size).