# Find the nearest locations for multiple points

Assume that there are many holes with their locations fixed, and the same number of balls distributed randomly. What is the smallest total distance for the balls fitting into the holes on the precondition that each hole can only fit one ball? For instance, the holes(black dots) are regularly distributed, and the balls(red dots) are randomly distributed. The nearest holes of the individual balls are indicated by arrows.

holes = Tuples[Range[1, 2, 1], 2];
balls = RandomReal[{1, 2}, Dimensions[holes]];
Graphics[{PointSize[Large], Point[holes], Red,PointSize[Medium],Point[balls]}] Thanks for all the helps and answers. The problem is called 'The Euclidean matching problem' or 'Euclidean minimum weight matching problem' 1. I recently found an approximate alogrithm which achieves nearly O(n) time complexity .

 A Near-Linear Constant-Factor Approximation for Euclidean Bipartite Matching

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;
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 -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole, ball -> hole, ball -> hole,
ball -> hole} *)


Visualize the result.

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls],
Line[assignments /. ball[i_] \[DirectedEdge] hole[j_] :> {balls[[i]], holes[[j]]}]}] • By the way, is there anything we can do about all the ugly instances of \[DirectedEdge] up there? – Rahul Oct 23 '14 at 15:53
• -> instead of \[DirectedEdge]? – ybeltukov Oct 23 '14 at 18:29
• By the way, can Sqrt@Total[Outer[Plus,balls,-holes,1]]^2,{3}] increase the speed? – ybeltukov Oct 23 '14 at 19:40
• @ybeltukov: No, if I do edges = Flatten@Table[...]~Join~... it is instantaneous while graph = Graph[edges] is the slow part. – Rahul Oct 23 '14 at 20:19
• Indeed... May be AdjacencyGraph? Because you almost construct the adjacency matrix. Your solution is very good but I believe it is possible to make it even better. – ybeltukov Oct 23 '14 at 20:35

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.

• A small bug detected~ When a point p exists in both balls and holes then the EuclideanDistance between this two ps would be 0, indicating there's no connection between them! It's better to add a small distance offset in the distance matrix: Outer[EuclideanDistance, pointsA, pointsB, 1]+1.*^-10 – Wjx Oct 11 '17 at 8:14
• There is a bug in MMA 11.3 and earlier that for certain ball/hole matrices, the FindMinimumCostFlow never finishes its calculations [CASE:4156292]. – berniethejet Sep 11 '18 at 2:24
• The problem occurs when there are identical balls or holes. The work-around is before generating the costMatrix to do If[Length[balls]>Length[Union[balls]],balls=balls+RandomReal[{-0.001,0.001},Dimensions[balls]]]. And likewise for the holes. – berniethejet Sep 13 '18 at 15:24

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.)

• It seems to fail on big problems (above 20 by 20 seems to be enough): "Could not process the constraints..." and "Lists {min,vals} and FindMinimum[{obj,Flatten[{c1,c2,c3}]},fvars] are not the same shape. >>" – P. Fonseca Jan 5 '15 at 12:34
• @P.Fonseca Some experimentation indicates you might do better using NMinimize in place of FindMinimum. I do not know if the result is guaranteed to be the same and I am looking into this. Also I would use Dispatch[vals] in the replacement because otherwise it is quite slow. – Daniel Lichtblau Jan 5 '15 at 19:55
• @P.Fonseca What version of Mathematica are you using? – Daniel Lichtblau Jan 5 '15 at 19:56
• 10.0.2 win64 sp1. I don't think i copied something wrong (since for less than 20, it worked fine), but if you can't repeat my problem (I first started with 50, and then went down to a rouded working limit), i can give it another try tomorrow, for a better insight. – P. Fonseca Jan 5 '15 at 21:08
• @P.Fonseca I'm looking into some issues that arose in 10.0.2, most notably a slowdown. I am not sure if the failure to converge that I can make happen (with some obscure settings at least) is indicative of anything wrong or just a difficult optimization. At 20x20 we have 400^2 variables, which is not necessarily overwhelming for IP, but certainly not a small problem. – Daniel Lichtblau Jan 6 '15 at 15:31

Minimization of the total Euclidean distance from balls to holes is quite difficult optimization problem. As a start point I propose a simple greedy algorithm:

• Find the nearest hole for every ball.
• Put balls to holes starting from the closest ball-hole pair while the corresponding holes are empty
• Repeat

It is not the best algorithm and may be I will update my post in the future.

n = {5, 5};
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium], Point[balls]}] res = Flatten[#, 2] &@ Last@Reap@
Module[{h = holes, b = balls, nh, nb, nf, num, put, fill,
nondup},
nb = nh = Range@Length@h;
While[h != {},
nf = Nearest[h -> Automatic];
num = nf /@ b // Flatten;
put = Ordering@Total[(h[[num]] - b)^2, {2}];
fill = num[[put]];
nondup = Floor@BinarySearch[Range@Length@fill, 1/2,
1 - Boole@DuplicateFreeQ@fill[[;; #]] &];
put = put[[;; nondup]];
fill = fill[[;; nondup]];
Sow@Transpose@{nb[[put]], nh[[fill]]};
b = Delete[b, Transpose@{put}];
nb = Delete[nb, Transpose@{put}];
h = Delete[h, Transpose@{fill}];
nh = Delete[nh, Transpose@{fill}];
]
]
(* {{8, 24}, {13, 5}, {19, 16}, {23, 10}, {1, 23}, {3, 17}, {5,
14}, {16, 21}, {17, 3}, {9, 8}, {7, 4}, {18, 11}, {4, 22}, {21,
15}, {11, 20}, {2, 6}, {14, 2}, {10, 7}, {15, 9}, {25, 1}, {12,
25}, {6, 18}, {22, 19}, {20, 13}, {24, 12}} *)

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ res]}] • @ ybeltukov, your answer is helpful and enlightening.The first step is ‘Closest pair of points problem’. But a simple counterexample is for 2 hole-ball pairs with the given coordinates: holes = {{0, 0}, {0, 2}}; balls = {{0, 1.01}, {2, 2}}; I think the farest balls should be considered first. – novice Oct 23 '14 at 8:32
• @novice Anyway, greedy algorithm is bad. It can't get the global minimum. I found that Hungarian algorithm is exactly what you want. It is a polynomial algorithm. Unfortunately I have no time to implement it now. This Wolfram demonstration can also be helpful. – ybeltukov Oct 23 '14 at 11:33
• The Hungarian algorithm solves what is known as the assignment problem. Apparently one can solve such problems directly in Mathematica since version 9. – Rahul Oct 23 '14 at 11:55
• @RahulNarain I knew that it should correspond to a graph problem! Very nice :) – ybeltukov Oct 23 '14 at 18:31
• Because of triangle inequality, no two edges of the optimal matching can cross each other. /*A Survey on Algorithms for Euclidean Matching*/ – novice Jan 6 '15 at 0:53

Ok this may end up being a greedy algorithm, though it was based on an old Dijkstra algorithm that I modified, but since balls and holes are not connected in any way is just greedy.

Since I don't have V10 installed here, I don't have the same BinarySearch[] functions (Must be different from Combinatorica), I couldn't compare against @ybeltukov. But I bet my is slower due to the indexing of distances, as a N^2 matrix, and finding the nearest using Min on it. However it might be educative to see another one.

Using his sample data,

n = {5, 5};
SeedRandom;
holes = N@Tuples@Range@n;
balls = RandomReal[{0, # + 1}, Times @@ n] & /@ n // Transpose;


If you find another way to index the distance, it may be faster.

dist = Table[
Norm[holes[[i]] - balls[[j]]], {i, 1, Length[holes]}, {j, 1,
Length[balls]}];
filledholes = ConstantArray[False, Length[holes]];
balldroped = ConstantArray[False, Length[balls]];
parm = Max[dist];

a = Reap[
While[And @@ filledholes != True,
posmin = First@Position[dist, Min[dist]];
If[filledholes[[posmin[]]] != True &&
balldroped[[posmin[]]] != True,
filledholes[[posmin[]]] = True;
balldroped[[posmin[]]] = True;
Sow[{posmin[], posmin[]}];
dist[[First@posmin, Last@posmin]] = parm + 1.0;
,
dist[[First@posmin, Last@posmin]] = parm + 1.0;
Continue[];
];
];]

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ (Last@Last@a)]}] The total distance being

Plus @@ (Norm[balls[[#]] - holes[[#2]]] & @@@ (Last@Last@a))


22.9346

With Timing it is

0.039872

In case you want maximize the distance, then easy to replace Min for Max, vis versa, and the adjustment to the index. In ybeltukov answer you would weighted the results with 1/nf.

dist = Table[
Norm[holes[[i]] - balls[[j]]], {i, 1, Length[holes]}, {j, 1,
Length[balls]}];
filledholes = ConstantArray[False, Length[holes]];
balldroped = ConstantArray[False, Length[balls]];
parm = Min[dist];
b = Reap[
While[And @@ filledholes != True,
posmax = First@Position[dist, Max[dist]];
If[filledholes[[posmax[]]] != True &&
balldroped[[posmax[]]] != True,
filledholes[[posmax[]]] = True;
balldroped[[posmax[]]] = True;
Sow[{posmax[], posmax[]}];
dist[[First@posmax, Last@posmax]] = 0.99*parm;
,
dist[[First@posmax, Last@posmax]] = 0.99*parm;
Continue[];
];
];]

Graphics[{PointSize[Large], Point[holes], Red, PointSize[Medium],
Point[balls], Arrow[{balls[[#]], holes[[#2]]} & @@@ (Last@Last@b)]}] The total distance being

Plus @@ (Norm[balls[[#]] - holes[[#2]]] & @@@ (Last@Last@b))


101.006

This may help:

locations=RandomInteger[{-10,10},{4,2}];
holes=RandomInteger[{-10,10},{4,2}];
ListPlot[{locations, holes}, PlotLegends -> {"locations", "holes"}] distanceFunc = Nearest[holes, DistanceFunction -> EuclideanDistance];
N@Total[EuclideanDistance[#, First@distanceFunc[#]] & /@ locations]

• thanks for your answer, what if each hole can only take one ball? – novice Oct 22 '14 at 7:49
• @novice ahh, I miss that – molekyla777 Oct 22 '14 at 8:06