Find the nearest locations for multiple points

Question

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 [2].

1 http://dl.acm.org/citation.cfm?id=1882725&CFID=469610786&CFTOKEN=72872074

[2] A Near-Linear Constant-Factor Approximation for Euclidean Bipartite Matching

Answer 1

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]]}]}]

Answer 2

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.

Answer 3

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

Answer 4

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]}]

Answer 5

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[1234];
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[[1]]]] != True && 
               balldroped[[posmin[[2]]]] != True,
                filledholes[[posmin[[1]]]] = True;
                balldroped[[posmin[[2]]]] = True;
                Sow[{posmin[[2]], posmin[[1]]}];
                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[[1]]]] != True && 
               balldroped[[posmax[[2]]]] != True,
                filledholes[[posmax[[1]]]] = True;
                balldroped[[posmax[[2]]]] = True;
                Sow[{posmax[[2]], posmax[[1]]}];
                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

Answer 6

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]