Message Boards Message Boards

Get a pair with smallest total distance from one list

Posted 7 years ago

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.

POSTED BY: Yode Japhe
18 Replies
Posted 7 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 BY: Paul Cleary

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 BY: Sander Huisman
Posted 7 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 BY: Yode Japhe

source of the webpage...

POSTED BY: Sander Huisman
Posted 7 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/1077817

and it works! No need to inspect the source. :)

POSTED BY: Alexey Popkov

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 BY: Sander Huisman
Posted 7 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.

POSTED BY: Yode Japhe
Posted 7 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}

Mathematica graphics

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 BY: Yode Japhe

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 BY: Sander Huisman
Posted 7 years ago
POSTED BY: Yode Japhe
Posted 7 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]}]

Mathematica graphics

It is not a perfect solution,because a global minimum pair always is ecpected actually.

POSTED BY: Yode Japhe

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

enter image description here

where the black are the old lines, and the blue ones the new lines.

POSTED BY: Sander Huisman
Posted 7 years ago

The blue ones always is the real minimum pair?

POSTED BY: Yode Japhe

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 BY: Sander Huisman
Posted 7 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 BY: Yode Japhe

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

enter image description here

Total sum distance of distance: 30.8036

We 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 BY: Sander Huisman
Posted 7 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 BY: Yode Japhe

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

enter image description here

One needs to swap all 4 of them together to get the real minimum:

enter image description here

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 BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract