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.
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 theNearestFunction
.I'm not too happy about having to put the
N
aroundEuclideanDistance
. Unfortunately,NearestFunction
doesn't accept something like $\sqrt2$ as the search radius.