Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Suppose we have a point set, the points are labeled 1,2,3,....n. And p[i] is the coordinate of the point. Now I want the label of nearest and next nearest point (maybe several) corresponding every point in this point set.

For example:

p[1] = {0, 0};
p[2] = {0, 1};
p[3] = {1, 0};
p[4] = {0, 5};

and I want a list for nearest neighbour like this:

{{2, 3}, {1}, {1}, {2}}

it means that p[2] and p[3] is nearest to p[1], p[1] is nearest to p[2], etc.

And I need a similar next nearest neighbour list.

I write the following code:

num = 5000;
Do[p[i] = RandomInteger[{1, 100}, 2], {i, 1, num}];

coolist = Table[p[i], {i, 1, num}];

Clear[position];
position[expr_] := 
  With[{positionData = 
     SortBy[#[[1, 1]] -> #[[All, 2]] & /@ 
        GatherBy[
         Extract[expr, #, Verbatim] -> # & /@ 
          Position[expr, _, Depth[expr]], First], 
       Min[Length /@ #[[2]]] &] // Dispatch}, 
   Replace[#, positionData] &];

poscoolist = position[coolist];

Clear[nnsite];
nnsite[k_, coolist_] := Module[{nncoolist},
   nncoolist = Nearest[Complement[coolist, {p[k]}], p[k]];
   Flatten@
    Table[poscoolist[nncoolist[[i]]], {i, 1, Length[nncoolist]}]];

nearestlist = Table[nnsite[i, coolist], {i, 1, num}]; // AbsoluteTiming

Clear[nnlabel];
Thread[Evaluate@Array[nnlabel, num] = nearestlist];

Clear[nnnsite];
nnnsite[k_, coolist_] := Module[{nnncoolist},
   nnncoolist = 
    Nearest[Complement[
      Delete[coolist, Partition[nnlabel[k], 1]], {p[k]}], p[k]];
   Flatten@
    Table[poscoolist[nnncoolist[[i]]], {i, 1, Length[nnncoolist]}]];

nextnearestlist = 
   Table[nnnsite[i, coolist], {i, 1, num}]; // AbsoluteTiming

the function position in the above code is provide by Mr.Wizard (see here). nearestlist and nextnearestlist give the result.

Notice to use my code, supply every coordinate p[i] of points.

For 5000 points, it takes 1 minutes. But for 10000 points, it takes 6 minutes. Quite long!

I feel that my code is so simple, there must be better ways that are more efficient.

share|improve this question
5  
It looks like you're calling Nearest several times. Why not just call it once to get a precomputed NearestFunction, and use that to find the three nearest sites to any given point? The first one will be the site itself, the next two will be the points you want. It will be much more efficient. (I'd post an answer but I don't have access to Mathematica right now.) –  Rahul Narain Dec 10 '13 at 0:52
    
@RahulNarain Oh, thank you. NearestFunction! Let me see it. –  matheorem Dec 10 '13 at 0:56
1  
I don't understand what you're asking. –  Rahul Narain Dec 10 '13 at 1:55
1  
I see, you want all the points at the nearest distance to $x$ in a single list. And then you want another list containing all the points at the second nearest distance. I suppose you could write a function that calls NearestFunction with larger and larger $n$ until it encounters a point at the third nearest distance... –  Rahul Narain Dec 10 '13 at 4:16
1  
Nearest also supports a Nearest[data, x, {n, r}] syntax which returns $n$ nearest points within radius $r$, see: mathematica.stackexchange.com/a/34895/131 –  Yves Klett Dec 10 '13 at 7:11
show 12 more comments

3 Answers

Here is a way to do it. At my first attempt I did not compile the functions and was a bit slower than your reported time. With the compiled function it is faster (at least on a 2Ghz/i7 Macbook Pro) : 5000 takes 23 secs at 10000 is 95 secs.

X := RandomInteger[{0, 100}];
n = 10000; data = Table[{X, X}, {i, 1, n}] ;

tStart = AbsoluteTime[];

"First, compute matrix of square distances (faster than taking the \
square roots) -- compiling is way faster"
dm = Compile[
   {{D, _Real, 2}, n},
   Table[
    (D[[i, 1]] - D[[j, 1]])^2 + (D[[i, 2]] - D[[j, 2]])^2,
    {i, 1, n}, {j, 1, n}]
   ];
distMatrix = dm[data, n]; // AbsoluteTiming

"Now get the first and second closest distances :"
firstDistances = 
   Table[Min[Select[distMatrix[[i]], # != 0 &]], {i, 1, 
     n}]; // AbsoluteTiming
secondDistances = 
   Table[Min[
     Select[distMatrix[[i]], # != 0 && # != 
         firstDistances[[i]] &]], {i, 1, n}]; // AbsoluteTiming

"Finally find the neighbours :"
firstDistanceIndices = 
   Table[Position[distMatrix[[i]], firstDistances[[i]]], {i, 1, 
     n}]; // AbsoluteTiming
secondDistanceIndices = 
   Table[Position[distMatrix[[i]], secondDistances[[i]]], {i, 1, 
     n}]; // AbsoluteTiming

"Done. Time elapsed :"
AbsoluteTime[] - tStart

Clearly the code is $O(n^2)$ and the performance reflects that : doubling $n$ about quadruples the execution time. $n^2/10^6$ is a reasonable approximation to the actual time, at least for $n\leq10000$.

share|improve this answer
    
Thank you for your method! I also posted mine :) –  matheorem Dec 10 '13 at 8:31
add comment

Here's a simple solution using a single precomputed NearestFunction; much faster than $O(n^2)$. I've written it assuming the sites are in a list ps, rather than embedded inside a function p, because I think this way is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]] (* so it returns the index of the site *)

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

Edit: I forgot that you want the neighbours of all the sites collected in a big list. Well, then you just do

nearestSites = nthNearestSites[#, 1] & /@ Range[Length[ps]];
nextNearestSites = nthNearestSites[#, 2] & /@ Range[Length[ps]];

On a hundred thousand sites, these take 2.9 and 6.8 seconds on my machine respectively. On a million, they will probably take a couple of minutes.

P.S.

  1. You could just define nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]], but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the NearestFunction.

  2. I'm not too happy about having to put the N around EuclideanDistance. Unfortunately, NearestFunction doesn't accept something like $\sqrt2$ as the search radius.

share|improve this answer
    
Thank you for your answer! I am testing now, and I also found the EuclideanDistance problem, so I use N and add small offset 0.0000001 to it, I found this small offset is necessary, see my answer. –  matheorem Dec 10 '13 at 8:57
    
I see your answer but it doesn't tell me why the offset is necessary. –  Rahul Narain Dec 10 '13 at 9:03
    
Well, you can try this ps={{0, 0}, {Sqrt[2], 0}, {0, Sqrt[2]}, {Sqrt[2], Sqrt[2]}}. Very odd,.but it's true, offset is necessary. –  matheorem Dec 10 '13 at 9:49
    
Ah... Well, another workaround is that if you set ps = N@{{0, 0}, ...} instead then it works. –  Rahul Narain Dec 10 '13 at 10:45
    
Ok, are you considering to improve your answer with N? Though my method is identical to yours, but I make the label a little complex, so I'll accept yours, if there is no better method –  matheorem Dec 10 '13 at 11:04
show 3 more comments

I figured a way myself using the NearestFunction feature, Thank you for those suggestions :)

num = 5000;
Do[p[i] = RandomInteger[{1, 100}, 2], {i, 1, num}];

coolist = Table[p[i], {i, 1, num}];

Clear[position];
position[expr_] := 
  With[{positionData = 
     SortBy[#[[1, 1]] -> #[[All, 2]] & /@ 
        GatherBy[
         Extract[expr, #, Verbatim] -> # & /@ 
          Position[expr, _, Depth[expr]], First], 
       Min[Length /@ #[[2]]] &] // Dispatch}, 
   Replace[#, positionData] &];

poscoolist = position[coolist];

nf = Nearest[coolist];
(
  Clear[nnsite, nnlabel];
  Do[nnsite[i] = 
    Drop[nf[p[i], {Infinity, 
       N@EuclideanDistance[Last[nf[p[i], 2]], p[i]] + 0.0000001}], 
     1], {i, 1, num}];
  Do[nnlabel[i] = 
    DeleteCases[DeleteDuplicates[Flatten[poscoolist /@ nnsite[i]]], 
     i], {i, 1, num}];

  nnlist = Table[Length@nnlabel[i], {i, 1, num}];

  Clear[nnnsite, nnnlabel];
  Do[nnnsite[i] = 
    Drop[nf[p[i], {Infinity, 
       N@EuclideanDistance[Last[nf[p[i], nnlist[[i]] + 2]], p[i]] + 
        0.0000001}], nnlist[[i]] + 1], {i, 1, num}];
  Do[nnnlabel[i] = 
    DeleteDuplicates[Flatten[poscoolist /@ nnnsite[i]]], {i, 1, num}];

  nearestlist = Table[nnlabel[i], {i, 1, num}];
  nextnearestlist = Table[nnnlabel[i], {i, 1, num}];
  ) // AbsoluteTiming

The performance is 5000 points takes 6 seconds, 10000 points takes 50 seconds.

a small offset such as 0.0000001 is necessary.

You can try {{0, 0}, {Sqrt[2], 0}, {0, Sqrt[2]}, {Sqrt[2], Sqrt[2]}} with and without offset to see the difference.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.