# Get a pair with smallest total distance from one list

Posted 5 years ago
10136 Views
|
18 Replies
|
12 Total Likes
|
 Cross post in Stack Exchange in here. I want to get a pair with smallest total distance from pts in following 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 to get a minimum cost perfect matching, but it seems the FindIndependentEdgeSet regards weighted graph as unweighted directly.It often make me depressed. Of course, I'm happy to know other methods that can do this which aren't based on Graph Theory.And I found a ready-made algorithm for min cost perfect matching here.If anybody can implement it in Mathematica,that will be an excited thing.Because as I know,that is a very efficient algorithm.I hope the answer can meet a additional demand:it can use similar method to find a maximal cost perfect matching,which is promising.Help to imporove Mathematica,Help me out,Please.
18 Replies
Sort By:
Posted 5 years ago
 My first try would just be swapping pairs of connections until we have an optimum: ClearAll[TrySwitch, TotalDistance, TrySwitch, FindPairs] TotalDistance[pts_List, conns_List] := Total[Norm[N[pts[[#1]] - pts[[#2]]]] & @@@ conns] TrySwitch[pts_List, conns_List, indices_List] := Module[{orig, newconns, subsets, ids}, orig = Sort[Sort /@ conns[[indices]]]; ids = Partition[#, 2] & /@ Permutations[Flatten[orig]]; ids = Union[Sort /@ Map[Sort, ids, {2}]]; ids = TakeSmallestBy[ids, TotalDistance[pts, #] &, 1][[1]]; newconns = Delete[conns, List /@ indices]; newconns~Join~ids ] TrySwitchAll[pts_List, conns_List] := Module[{subsets, newerror, j = 0, olderror, newconns = conns}, subsets = Subsets[Range[Length[conns]], {2}]; newerror = 1; olderror = -1; While[olderror =!= newerror \[And] j < 1000 (* safeguard *), j++; Do[ newconns = TrySwitch[pts, newconns, s]; , {s, subsets} ]; olderror = newerror; newerror = TotalDistance[pts, newconns]; ]; newconns ] FindPairs[pts_List] := Module[{candidates, out, tmp, left}, If[EvenQ[Length[pts]] \[And] Length[pts] >= 2, TrySwitchAll[pts, Partition[Range[Length[pts]], 2]] , Print["Even number of points >=2 needed!"] ] ] Which works with your example: pts = {{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}}; out = FindPairs[pts] TotalDistance[pts, out] Graphics[{Red, Point[pts], Black, Line[pts[[#]] & /@ out]}] Total sum distance of distance: 30.8036We can make it slightly smarter by starting with a more intelligent initial choice of pairs. Here I look for the nearest 7 neighbours of each point forming 7*n candidate pairs, we sort these candidates by distance, and then greedily pick them; shortest distance first. Only picking each point once... ClearAll[FindPairs] FindPairs[pts_List] := Module[{candidates, out, tmp, left}, If[EvenQ[Length[pts]] \[And] Length[pts] >= 2, candidates = Transpose[{Range[Length[pts]], Nearest[pts -> "Index", pts, Min[Length[pts], 8]]}]; candidates = With[{x = #1, y = Rest[#2]}, {x, #} & /@ y] & @@@ candidates; candidates = Join @@ candidates; candidates = {#, TotalDistance[pts, #]} & /@ candidates; candidates = SortBy[candidates, Last][[All, 1]]; out = Reap[While[Length[candidates] > 0, tmp = candidates[[1]]; candidates = Rest[candidates]; candidates = DeleteCases[candidates, {___, tmp[[1]], ___} | {___, tmp[[2]], ___}]; Sow[tmp]; ]]; out = out[[2, 1]]; left = Complement[Range[Length[pts]], Flatten[out]]; out = out~Join~Partition[left, 2]; out = TrySwitchAll[pts, out]; out , Print["Even number of points >=2 needed!"] ] ] Which gives the same answer. Note that any good solution, that there are no lines crossing; for if there were crossing lines we could swap the points to get a smaller total. See also https://www.youtube.com/watch?v=GX3URhx6i2E around the 28:30... for a similar problem.
Posted 5 years ago
 Brilliant try.Thanks very much.But if I give some test pts for test,I found it not always give a right answer?Such as SeedRandom[1] pts = RandomReal[100, {8, 2}]  {{81.7389,11.142},{78.9526,18.7803},{24.1361,6.57388},{54.2247,23.1155},{39.6006,70.0474},{21.1826,74.8657},{42.2851,24.7495},{97.7172,82.5163}} If we use a brute force method to find a definitely right answer: First[pairs = MinimalBy[Partition[#, 2] & /@ Permutations[pts], Total[EuclideanDistance @@@ N[#]] &]] First[Total[EuclideanDistance @@@ N[#]] & /@ pairs]  {{{81.7389,11.142},{78.9526,18.7803}},{{24.1361,6.57388},{42.2851,24.7495}},{{54.2247,23.1155},{97.7172,82.5163}},{{39.6006,70.0474},{21.1826,74.8657}}}126.475(This is total distance) But If we use your FindPairs pts[[#]] & /@ FindPairs[pts] Total[EuclideanDistance @@@ N[pts[[#]] & /@ FindPairs[pts]]]  {{{24.1361,6.57388},{54.2247,23.1155}},{{81.7389,11.142},{78.9526,18.7803}},{{39.6006,70.0474},{21.1826,74.8657}},{{42.2851,24.7495},{97.7172,82.5163}}}141.565(This is total distance)
Posted 5 years ago
 Correct, that is partly because the end-condition was not proper. Here I fixed that: ClearAll[CanonicalPairs,TrySwitch,TotalDistance,TrySwitch,TryAllSubsets,FindPairs] CanonicalPairs[conns_List]:=Sort[Sort/@conns] TotalDistance[pts_List,conns_List]:=Total[Norm[N[pts[[#1]]-pts[[#2]]]]&@@@conns] TrySwitch[pts_List,conns_List,indices_List]:=Module[{orig,newconns,subsets,ids}, orig=CanonicalPairs[conns[[indices]]]; ids=Partition[#,2]&/@Permutations[Flatten[orig]]; ids=DeleteDuplicates[CanonicalPairs/@ids]; ids=TakeSmallestBy[ids,TotalDistance[pts,#]&,1][[1]]; newconns=Delete[conns,List/@indices]; newconns~Join~ids ] TryAllSubsets[pts_List,conns_List,subsets_List]:=Module[{inconns,tmp}, inconns=CanonicalPairs[conns]; Do[ tmp=CanonicalPairs[TrySwitch[pts,inconns,s]]; If[tmp=!=inconns,Return[tmp,Module]]; , {s,subsets} ]; Return[inconns] ] TrySwitchAll[pts_List,conns_List]:=Module[{subsets}, subsets=Subsets[Range[Length[conns]],{2}]; NestWhile[TryAllSubsets[pts,#,subsets]&,conns,(UnsameQ@@(TotalDistance[pts,#]&/@{##}))&,3,1000] ] FindPairs[pts_List]:=Module[{candidates,out,tmp,left},If[EvenQ[Length[pts]]\[And]Length[pts]>=2,TrySwitchAll[pts,Partition[Range[Length[pts]],2]],Print["Even number of points >=2 needed!"]]] This is much neater code, but works in a similar way. It tries switching pairs (that is the '2' inside TrySwitchAll) until it stabilises. But sometimes, one needs to swap a triplet (3), or quadruplet (4) to get to the minimum (rare, but can happen). Basically I'm tickling the partial solution, but what one sometimes needs is a huge kick to get it out of its local minimum.See e.g. this example: pts = {{57.526218026559405, 15.570452206207037}, {90.54332967253814, 52.25381104673323}, {32.476740815657735, 7.28878961340456}, {91.39816425481379, 37.81652745685858}, {19.24680608590677, 95.04441305558967}, {83.85131355308184, 58.07174181385142}, {91.42216182515247, 25.28663683844907}, {71.81082383731945, 10.795462957661243}}; out = FindPairs[pts] TotalDistance[pts, out] Graphics[{Red, Point[pts], Black, Line[pts[[#]] & /@ out]}] One needs to swap all 4 of them together to get the real minimum:O btw, I never claimed that this algorithm was perfect, but it should give reasonable (local) minima quite quickly. I was aware that there can be local minima... I'm not sure if you can get the real minimum without a very convoluted algorithm that takes lots and lots of time...
Posted 5 years ago
 Yep,I had realized it is very hard before arriving here. :) So I'm hope to implement that ready-made algorithm for min cost perfect matching in Mathemtica.But it is a pitty that I cannot read any C++ code.And as my try,we cannot find a pair even in 70 points if we use FindPairs,but that algorithm can make it in 1000 points within one second as my friend's test.
Posted 5 years ago
 I have a method based on FindHamiltonianPath.It has a not very bad efficiency,which can find a pair(it is not always a global minimum pair as my test) from 100 points within 1s SeedRandom[1] pts = RandomReal[100, {100, 2}]; g = CompleteGraph[Length[pts]]; verGraph = VertexReplace[g, Thread[VertexList[g] -> pts]]; weighGraph = Graph[verGraph, EdgeWeight -> EuclideanDistance @@@ EdgeList[verGraph]]; pair = Partition[FindHamiltonianPath[weighGraph], 2]; Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}] It is not a perfect solution,because a global minimum pair always is ecpected actually.
Posted 5 years ago
 I tried 100, 200 and 300 points, they all work fine. Good luck translating that huge code, the code is not hard, but without comments you have no idea what each subroutine accomplishes...Using FindHamiltonianPath/FindShortestTour will not give the minimum, because you only consider half of the edges of the tour. But you CAN use it as start for my FindPairs function: FindPairs[pts_List] := Module[{candidates, out, tmp, left, start}, If[EvenQ[Length[pts]] \[And] Length[pts] >= 2, start = Most[FindShortestTour[pts][[2]]]; TrySwitchAll[pts, Partition[start, 2]] , Print["Even number of points >=2 needed!"] ] ] As for the above example, there are clearly improvements that can be made to this solution: SeedRandom[1] pts = RandomReal[100, {100, 2}]; g = CompleteGraph[Length[pts]]; verGraph = VertexReplace[g, Thread[VertexList[g] -> pts]]; weighGraph = Graph[verGraph, EdgeWeight -> EuclideanDistance @@@ EdgeList[verGraph]]; pair = Partition[FindHamiltonianPath[weighGraph], 2]; pts = Join @@ pair; conns = Partition[Range[Length[pts]], 2]; orig = Graphics[{Line[pts[[#]]] & /@ conns, Red, PointSize[.015], Point[pts]}]; newconn = TrySwitchAll[pts, conns]; new = Graphics[{Blue, Line[pts[[#]]] & /@ newconn}]; Show[{orig, new}] You method, then improved using mine...where the black are the old lines, and the blue ones the new lines.
Posted 5 years ago
 The blue ones always is the real minimum pair?
Posted 5 years ago
 A better minimum as compared to the black one. I'm not sure if it is THE minimum. You can see that some pairs 'switched', especially clear in the top right with the 'triangle'.
Posted 5 years ago
 I have changed the title,which is more clear for what I want directly.With the help of the help of Mr. Sander Huisman, I can even find the global minimum pair(I cannot sure it is the global minimum actually) from 200 points AbsoluteTiming[newconn = TrySwitchAll[pts, conns];] new = Graphics[{Blue, Line[pts[[#]]] & /@ newconn, PointSize[.01], Red, Point[pts]}]  {25.5138, Null} Well,I have to say there is long way to go.My original target is find the global minimum pair from about 3000+ points. So a better solution is ecpected still.
Posted 5 years ago
 What is this actually used for? I'm aware of connecting two sets of points, but i have not heard of this problem. How critical is it to get the global minimum versus sufficiently nice (practical) local minimum. I have a feeling this problem is very much akin to the travelling salesman problem, which for 50 or so points also already get tricky to ensure to get the global optimum. The one in mathematica switches to a 'good enough' algorithm for large number of points. The reference you cited in the opening post; do they really ensure it is a global minimum? the paper is hard to read (for me), as I'm not really into the field...
Posted 5 years ago
 Oh,a minimum cost perfect matching can serve many case,just we cannot realize it sometimes.I can think out right now another two examples I have encoutered. [Share a method to patch the discontinuous line in a image](https://mathematica.stackexchange.com/questions/109870/share-a-method-to-patch-the-discontinuous-line-in-a-image) [Reconstruct text pages cut by shredder](https://mathematica.stackexchange.com/a/140407/21532) Actually I can consider it as a minimum cost perfect matching question,and in such case,a local minimum, but not a global minimum, will lead to a disastrous result.And the link,I cited, work for this indeed as I know.
Posted 5 years ago
 Would FindShortestTour be a better choice, (if I have read this correctly) as an example AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; pair = Partition[pts[[Last[FindShortestTour[pts]]]], 2]; Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; Total[Table[ Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]] or by rotating the points by 1 to get the other possible set of pairs AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; pair = Partition[pts[[Last[FindShortestTour[pts]]]], 2]; pair = Partition[RotateLeft[Flatten[pair, 1]], 2]; Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; Total[Table[ Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]] in comparison to this method AbsoluteTiming[SeedRandom[1]; pts = RandomReal[100, {100, 2}]; g = CompleteGraph[Length[pts]]; verGraph = VertexReplace[g, Thread[VertexList[g] -> pts]]; weighGraph = Graph[verGraph, EdgeWeight -> EuclideanDistance @@@ EdgeList[verGraph]]; pair = Partition[FindHamiltonianPath[weighGraph], 2]; Print[Graphics[{Line[pair], Red, PointSize[.015], Point[pts]}]]; Total[Table[ Sqrt[(pair[[q, 1, 1]] - pair[[q, 2, 1]])^2 + (pair[[q, 1, 2]] - pair[[q, 2, 2]])^2], {q, 1, Length[pair]}]]] The shortest tour option is also much faster, 3000 points is under a second on my machine.
Posted 5 years ago
 This does not give the minimum pairs, see my reply above: http://community.wolfram.com/groups/-/m/t/1075278 Mainly you're minimizing the total length rather then the even (or odd) pieces. But it can be a good starting point for other algorithms to improve on it. But certainly does not guarantee minimal solution.
Posted 5 years ago
 Could I know how do you get that url(http://community.wolfram.com/groups/-/m/t/1075278) about a reply under a post?
Posted 5 years ago
 source of the webpage...
Posted 5 years ago
 source of the webpage... It looks like we can get the message ID by copying the "Reply" hyperlink under the post. For example, for your post by copying I get javascript:_19_addQuickReply('reply',%20'1077817');setQuickReplyTo(1077817);  where 1077817 obviously is the message ID. Then for generating the link we should just prepend http://community.wolfram.com/groups/-/m/t/:http://community.wolfram.com/groups/-/m/t/1077817and it works! No need to inspect the source. :)
Posted 5 years ago
 That's kinda what I meant, if you click 'inspect element' on most parts of the post, you can see the ID all over the place...
Posted 5 years ago
 Thanks for your endeavor for my this question.And I think a good message deserve to bring to you.I receive a good solution here based on FindShortestTour.