7

This is similar to this post, but actually more complicated. I have two lists in that question, and I want to get a pair. But in here, I want to get a pair from one list (even points). Given these 20 points:

SeedRandom[1]
pts = RandomInteger[20, {20, 2}]

{{5,0},{7,0},{2,3},{0,0},{16,14},{3,8},{19,5},{18,16},{12,0},{19,4},{7,3},{0,4},{20,3},{5,12},{19,8},{11,2},{3,10},{4,2},{17,11},{15,6}}

I can use FindIndependentEdgeSet to get one pair like this:

g = CompleteGraph[20];
List @@@ FindIndependentEdgeSet[VertexReplace[g, Thread[VertexList[g] -> pts]]]

{{{5,0},{7,0}},{{2,3},{3,8}},{{0,0},{16,14}},{{19,5},{4,2}},{{18,16},{3,10}},{{12,0},{7,3}},{{19,4},{11,2}},{{0,4},{20,3}},{{5,12},{19,8}},{{17,11},{15,6}}}

The total distance is 

Total[EuclideanDistance @@@ pair] // N

113.859

But I'm sure it is not the smallest pair. Actually, I think I need a FindIndependentEdgeSet of edge weight version, but it seems the FindIndependentEdgeSet regards weighted graph as unweighted directly. Can anyone give me advice? Of course, I'm happy to know other methods that can do this which aren't based on Graph Theory.

This post related this question.

yode
  • 26,686
  • 4
  • 62
  • 167
  • Are you interested only in a perfect solution, i.e. the minimal weight pairing, or would an optimization approach that produces a "good" pairing be useful too? – Mr.Wizard Apr 09 '17 at 07:36
  • @Mr.Wizard I just expected the optimal paring and hope to speed up it. Thanks. :) – yode Apr 09 '17 at 09:21
  • 2
    @Mr.Wizard A brute force look shows the min distance with SeedRandom[1] to be 30.8036. This takes ~100s to run for the 20 point case and scales badly (Binomial[2 n, n])... – Quantum_Oli Apr 09 '17 at 11:43

5 Answers5

5

There is a solution based on IGLargestIndependentVertexSets,which from Szabolcs's package IGraphM,when the point is less 15 or much less.I mean,the effeciency it very poor.I will produce just 10 points to test it.

Built a weight graph with EuclideanDistance 

SeedRandom[1]
pts = RandomInteger[20, {10, 2}];
g = CompleteGraph[10]

Mathematica graphics

Find a smallest independent vertex sets

iden = IGLargestIndependentVertexSets[LineGraph[g]];
pair = List @@@ First[MinimalBy[EdgeList[g][[##]] & /@ iden, 
    Total[EuclideanDistance @@@ N[#]] &]]

{{{5,0},{0,0}},{{7,0},{12,0}},{{2,3},{3,8}},{{16,14},{18,16}},{{19,5},{19,4}}}

N[Total[EuclideanDistance @@@ pair]]

18.9274

I hope this is a good start..

Community
  • 1
yode
  • 26,686
  • 4
  • 62
  • 167
4

I don't know if this is of any value at all but I was curious to see how an optimization approach would work, and this is what I came up with:

SeedRandom[1]
pts = RandomInteger[20, {20, 2}];

tdist = N @ Total[EuclideanDistance @@@ #] &;

fn[pt_List, x_ /; VectorQ[x, NumberQ]] := tdist @ Partition[pt[[Ordering @ x]], 2]

vars = Unique @ Table["p", {Length @ pts}];

{min, sol} = 
  NMinimize[fn[pts, vars], vars, Method -> "SimulatedAnnealing", 
   MaxIterations -> 1000, PrecisionGoal -> 1];

pts[[Ordering[vars /. sol]]] ~Partition~ 2

tdist @ %

{{{12, 0}, {11, 2}}, {{0, 0}, {0, 4}}, {{2, 3}, {4, 2}}, {{3, 10}, {5, 12}}, 
{{18, 16}, {16, 14}}, {{17, 11}, {19, 8}}, {{3, 8}, {7, 3}}, {{15, 6}, 
 {19, 5}}, {{5, 0}, {7, 0}}, {{19, 4}, {20, 3}}}

31.675

The same code starting with pts = RandomInteger[20, {10, 2}]; gives:

{{{19, 5}, {19, 4}}, {{2, 3}, {3, 8}}, {{12, 0}, {7, 0}},
 {{0, 0}, {5, 0}}, {{18, 16}, {16, 14}}}

18.9274
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
4

FindShortestTour appears to be a sub-optimal but fast method.

SeedRandom[1]
pts = RandomInteger[20, {20, 2}];

tdist = N@Total[EuclideanDistance @@@ #] &;

o = FindShortestTour[pts][[2]];

Partition[# @ pts[[o]], 2] & /@ {Most, Rest};

pairs = MinimalBy[%, tdist, 1][[1]]

tdist @ pairs

ListLinePlot @ pairs

{{{5, 0}, {4, 2}}, {{2, 3}, {0, 0}}, {{0, 4}, {3, 8}}, {{3, 10}, {5, 12}},
 {{16, 14}, {18, 16}}, {{17, 11}, {19, 8}}, {{19, 5}, {20, 3}},
 {{19, 4}, {15, 6}}, {{12, 0}, {11, 2}}, {{7, 3}, {7, 0}}}

32.0483

enter image description here

The same code with pts = RandomInteger[20, {10, 2}]; gives:

{{{7, 0}, {12, 0}}, {{19, 4}, {19, 5}}, {{18, 16}, {16, 14}},
 {{3, 8}, {2, 3}}, {{0, 0}, {5, 0}}}

18.9274
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • @yode I did not look at your self-answer closely enough! Why does your output only have ten points in it? Shouldn't it have all 20? I think I misunderstand your problem. :-( – Mr.Wizard Apr 09 '17 at 08:57
  • Sorry,my solution is very slow,so I use RandomInteger[20, {10, 2}] to instead of 20. – yode Apr 09 '17 at 09:00
  • @yode I see. :-) I guess you are interested in only the optimal paring and not one that is "close" -- is that correct? If so I should delete my answers as neither method is optimal. – Mr.Wizard Apr 09 '17 at 09:04
  • Yes,I just expected the optimal paring and hope to speed up it. :) – yode Apr 09 '17 at 09:06
  • @yode I believe your self-answer is also not an optimal pairing, so I shall not delete my answers at this time. – Mr.Wizard Apr 09 '17 at 09:10
  • See my update please. :) – yode Apr 09 '17 at 09:19
  • Actually I am very happy to see the FindShortestTour is used here,but I'm just little suspect that when we get the shortest loop,then we can get the shorest pair from it?But I don't very sure indeed. – yode Apr 09 '17 at 09:27
  • 1
    @yode No, I do not think this will always be the shortest pairing, which is why I wrote *sub-optimal*. I suspect that to find the true shortest pairing every time will require exhaustive testing and be very slow, as this is true for the Traveling Salesman Problem to the best of my knowledge. – Mr.Wizard Apr 09 '17 at 09:36
  • As your hint,actually I have tried FindHamiltonianPath,it is often give a optimal pairing,but not always. – yode Apr 13 '17 at 14:12
3

Might not be terribly efficient but this can be done with integer linear programming. It is quite similar to a method shown here but unfortunately the relaxed problem in this case need not deliver an integer solution.

Here is the example.

SeedRandom[1]
len = 20;
pts = RandomInteger[20, {len, 2}];

For each 1<=i<j<=20 we'll create a variable to denote whether or not to create an edge between the ith and jth points. So these are intended as 0-1 variables. We impose constraints that every point be hit by exactly one edge. The objective function will be the sum of all distances between points linked by an edge.

ovars = Array[x, {len, len}];
vars = MapIndexed[#1[[#2[[1]] + 1 ;; -1]] &, ovars];
fvars = Flatten[vars];
zeroRule = Thread[Complement[Flatten[ovars], fvars] -> 0];

dists = Table[
   vars[[j, k]]*EuclideanDistance[pts[[j]], pts[[j + k]]], {j, 
    len - 1}, {k, 1, Length[vars[[j]]]}];
obj = Total[dists, 2];
c1 = Thread[
   Map[Total, vars] + (Map[Total, Transpose[ovars]] /. zeroRule) == 1];
c2 = Map[0 <= # <= 1 &, fvars];

We solve it with FindMinimum.

{min, vals} = 
  FindMinimum[{obj, Flatten[{c1, c2, Element[fvars, Integers]}]}, 
   fvars];

Now deduce which edges get used.

ores = Position[Round[vars /. vals], 1, 2];
res = Map[{#[[1]], #[[2]] + #[[1]]} &, ores]

Again cribbing from someone else's graphics, we show the result.

Graphics[{PointSize[Large], Point[pts], Red, PointSize[Medium], 
  Arrow[{pts[[#]], pts[[#2]]} & @@@ res]}]

enter image description here

This method handles 150 points in a half minute or so. One bottleneck is the computation of the objective function.

Daniel Lichtblau
  • 58,970
  • 2
  • 101
  • 199
  • I have seen your answer many times,give me some time to check the FindMinimum always give a global minimal pair.I will response you after some days. – yode May 04 '17 at 00:18
3

This gives the same (optimal) result as Daniel Lichtblau's answer for the given example. Integer linear programming is also used, but in this answer this is done by FindShortestTour.

max = 20;
nn = 20;
SeedRandom[1]
pts = RandomInteger[max, {nn, 2}]
boolPts =
  Join[
   Transpose@Join[{-Range[nn]}, Transpose@ pts], 
   Table[{i, 0, 0}, {i,  nn}]];
bound = 2 max^2;
tour = FindShortestTour[
   boolPts,
    DistanceFunction -> (If[#[[1]] +  #2[[1]] == 0,
       -bound,
       If[
        #[[1]] #2[[1]] < 0 
        ,
        bound
        ,
        EuclideanDistance[#[[2 ;;]], #2[[2 ;;]]]
        ]] &), Method -> "IntegerLinearProgramming"];

resultLength = nn*bound + N@First@tour

30.8036

pairs = pts[[#]] & /@ Partition[tour[[2]], 2][[1 ;; -1 ;; 2]]
resultLength == Total[EuclideanDistance @@@ pairs]

{{{5, 0}, {7, 0}}, {{5, 12}, {3, 10}}, {{0, 0}, {2, 3}}, {{11, 
   2}, {12, 0}}, {{19, 4}, {20, 3}}, {{17, 11}, {19, 8}}, {{15, 
   6}, {19, 5}}, {{16, 14}, {18, 16}}, {{7, 3}, {4, 2}}, {{3, 8}, {0, 
   4}}}
 True
Jacob Akkerboom
  • 12,215
  • 45
  • 79
  • Some question should be listed.Firstly,the boolPts can be Catenate[MapIndexed[{PadRight[#2, 3], Prepend[#1, -First[#2]]} &, pts]]?Secondly,why we set the bound is 2 max^2?.Thirdly,Do you sure the FindShortestTour always will give Global minimal pairs and have you see this?fourthly,could you explain your DistanceFunction? – yode May 04 '17 at 00:39
  • The result is not same with this comment.And how to get the pair as your tour? – yode May 04 '17 at 00:55
  • @yode 1) yes, you can use that Catenate code if you prefer. 2) bound is set heuristically to 2 max^2. It can be seen as a reward/penalty that has is so big that if you can get it as a reward, you will always choose to do so and similarly you always avoid it as a penalty. If you do not include a reward, my procedure may find triplets instead of pairs. It is also used as a penalty to "remove edges" between vertices. 3) I think it will always find the optimum, I think the reduction to a FindshortestTour (traveling salesman) problem is correct and... – Jacob Akkerboom May 04 '17 at 07:02
  • FindShortestTour with the option Method -> "IntegerLinearProgramming" will give the optimal solution for the problem it is considering. 4) See my explanation under "2)". Re: "The result is not same...", well I get the same decimals right (30.8036)? Re: "how to get the pair(s)..", it is stored in the variable pairs. – Jacob Akkerboom May 04 '17 at 07:05
  • @yode it was probably a bit silly of me to not use a Graph. I felt it might be more efficient to not have to calculate all distances between all points, but perhaps the time this takes is negligible anyway. If I make a graph the intended connections between vertices may be more apparent. – Jacob Akkerboom May 04 '17 at 07:10
  • 1)If I use the Catenate version,I will get a wrong pairs. 2)I don't very understand your DistanceFunction still.I think it define two rules.First rule,it will connected that {n,*,*} and {-n,0,0}.Second rule,it will prefer to connect those same winding points.But I don't know why it will get a optimal pair.. – yode May 05 '17 at 00:48
  • @yode if Catenate gives the wrong pairs, then probably I should use Select instead of Partition and Part to get the pairs from the tour, I will test that in a minute. – Jacob Akkerboom May 05 '17 at 06:50
  • About why it gives an optimal pair, maybe it helps to look at the variable tour. You see that we always have two integers <= 20 then two integers >20. The integers <= 20 correspond to the original points and they come in pairs in tour in that they are separated by integers >20 (imaginary points). To see why we get pairs, note that we have 20 edges with a reward of bound, so the algorithm wants to jump exactly 20 times between points {n,*,*} and {-n,0,0}... – Jacob Akkerboom May 05 '17 at 06:50
  • What ends up happening is that from an imaginary point {j,*,*}, if possible, we jump to its corresponding real point {-j,*,*}. After that we jump to another real point {-k,*,*} (the only connected imaginary point is the point we came from, that is already visited so not allowed) and then we jump back to an imaginary point... – Jacob Akkerboom May 05 '17 at 06:51
  • If we were not to immediately jump to the imaginary point {k,*,*}, then then edge {{-k,*,*}, {k,*,*}} can never be used, so that we miss a reward of bound, which is not optimal. This is why we get a tour with this structure of pairs of real points and pairs of imaginary points. To see that it is optimal, note the following. All configurations of pairs correspond to a correct tour, for example {{1,5},{3,2}} corresponds to {1,5 25, 23, 3, 2, 22, 21, 1}. Furthermore, the cost of any such tour is 20*bound + costOfPairs. As 20*bounds is constant, we are optimizing costOfPairs – Jacob Akkerboom May 05 '17 at 06:51
  • This can help?And Options[GraphComputation`LegacyFindShortestTour] also accept a custom DistanceFunction – yode May 09 '17 at 15:09