Find the nearest locations for multiple points
Note: Please use Quantum_Oli's answer instead, which is a much faster implementation.
This is an instance of the assignment problem, which is a special case of the minimum-cost flow problem, which can be solved directly in Mathematica.
n = {5, 5};
SeedRandom[1234];
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;
Construct the bipartite graph between balls and holes with edge costs equal to the distances between them, and add two dummy "source" and "target" vertices. Strangely, this is the most time-consuming part.
graph = Graph[
Flatten@Table[
Property[ball[i] \[DirectedEdge] hole[j],
EdgeCost -> EuclideanDistance[balls[[i]], holes[[j]]]],
{i, Length@balls}, {j, Length@holes}]
~Join~
Table[Property[source \[DirectedEdge] ball[i], EdgeCost -> 0], {i, Length@balls}]
~Join~
Table[Property[hole[j] \[DirectedEdge] target, EdgeCost -> 0], {j, Length@holes}]];
Solve the minimum-cost flow problem.
assignments =
Cases[FindMinimumCostFlow[graph, source, target, "EdgeList"],
ball[_] \[DirectedEdge] hole[_]]
(*{ball[1] -> hole[18], ball[2] -> hole[15], ball[3] -> hole[1],
ball[4] -> hole[8], ball[5] -> hole[2], ball[6] -> hole[25],
ball[7] -> hole[16], ball[8] -> hole[11], ball[9] -> hole[10],
ball[10] -> hole[22], ball[11] -> hole[23], ball[12] -> hole[5],
ball[13] -> hole[6], ball[14] -> hole[24], ball[15] -> hole[12],
ball[16] -> hole[4], ball[17] -> hole[19], ball[18] -> hole[9],
ball[19] -> hole[21], ball[20] -> hole[13], ball[21] -> hole[3],
ball[22] -> hole[14], ball[23] -> hole[17], ball[24] -> hole[20],
ball[25] -> hole[7]} *)
Visualize the result.
Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls],
Line[assignments /. ball[i_] \[DirectedEdge] hole[j_] :> {balls[[i]], holes[[j]]}]}]
A faster version of Rahul's Answer
This question and the answers herein have just helped me solve a very similar problem using a solution based on @Rahul's answer which I find very elegant.
However, as is discussed briefly in the comments on that answer the code given (though very easy to understand) starts to run slow for larger numbers of balls and holes. I needed to perform solve tens of assignment problems with upwards of 50 balls and holes in each. Rahul's code was taking around 12 seconds to construct the Graph
for the 25 ball case, the solution below runs in 0.0012 seconds - I think faster than any of the others - I post it here for people looking for a fast solution in the future!
The principle is as @ybeltukov suggests that rather than generate the graph we simply write down the weighted adjacency matrix for the situation and pass it to FindMinimumCostFlow
as a cost matrix:
SourceTargetCostMatrix[pointsA_, pointsB_] :=
Module[{lA = Length[pointsA], lB = Length[pointsB]},
ArrayFlatten@{
{0, ConstantArray[1, {1, lA}], ConstantArray[0, {1, lB}], 0},
{ConstantArray[0, {lA, 1}], ConstantArray[0, {lA, lA}],
Outer[EuclideanDistance, pointsA, pointsB, 1],
ConstantArray[0, {lA, 1}]},
{ConstantArray[0, {lB, 1}], ConstantArray[0, {lB, lA}],
ConstantArray[0, {lB, lB}], ConstantArray[1, {lB, 1}]},
{0, ConstantArray[0, {1, lA}], ConstantArray[0, {1, lB}], 0}
}
]
costMatrix = SourceTargetCostMatrix[balls, holes];
assignments = Cases[
FindMinimumCostFlow[costMatrix, 1, Length[costMatrix], "EdgeList"],
x_ \[DirectedEdge] y_ /; x != 1 && y != Length[costMatrix]
];
Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls],
Line[assignments /. i_ \[DirectedEdge] j_ :> {balls[[i - 1]], holes[[j - Length[balls] - 1]]}]
}]
It can solve a 25x25 grid in just over a second, roughly a quarter of the time is for SourceTargetCostMatrix
, the remaining is FindMinimumCostFlow
:
Different Cost Functions
My problem required more priority placed on assigning those balls nearest to holes to that corresponding hole at the cost of having a few balls a very long way from a hole. I therefore used the Log
of the EuclideanDistance
which worked very well. To do so obviously simply replace EuclideanDistance
in SourceTargetCostMatrix
by whatever cost function you would like.
I will crib shamelessly from example and code for illustrating by @ybeltukov.
The example:
n = {5, 5};
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;
We can solve this as a linear programming problem. It looks like an integer linear program, but these are known to be solvable as relaxations thereof, that is, solutions to the relaxed LP will be integer valued (provided the solution is unique).
We set up the problem as below, to use FindMinimum
. That will invoke linear programming. I find it easier to formulate in terms of variables rather than explicit matrix and vector constraints.
len = Length[holes];
vars = Array[x, {len, len}];
fvars = Flatten[vars];
c1 = Thread[Total[vars] == 1];
c2 = Thread[Total[vars, {2}] == 1];
c3 = Map[0 <= # <= 1 &, fvars];
dists = Table[
vars[[j, k]]*EuclideanDistance[balls[[j]], holes[[k]]], {j,
len}, {k, len}];
obj = Total[dists, 2];
Now we solve it. The option setting is for speed. When the problem size is a bit larger than this, it will use interior point anyway, but for this size the automatic mode makes a slower choice.
SetOptions[LinearProgramming, Method -> "InteriorPoint"];
{min, vals} = FindMinimum[{obj, Flatten[{c1, c2, c3}]}, fvars];
res = Position[Round[vars /. vals], 1, 2]
(* {{1, 10}, {2, 8}, {3, 7}, {4, 23}, {5, 20}, {6, 2}, {7, 6}, {8,
11}, {9, 18}, {10, 24}, {11, 1}, {12, 19}, {13, 25}, {14, 14}, {15,
17}, {16, 22}, {17, 4}, {18, 15}, {19, 5}, {20, 12}, {21, 21}, {22,
13}, {23, 16}, {24, 3}, {25, 9}} *)
The picture:
Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ res]}]
(Disclosure: Had this looked incorrect, I would have thrown transposes into the formulation of the objective until I got it right.)