## Distances between points in periodic cube

10

2

Question

How can one implement more efficiently/elegantly/memory savvily the following function which returns a matrix of all Euclidian distances between points in 3D within a cube of width size, while accounting for the periodicity of the cube?

I have tried

Clear[PeriodicDistance];
PeriodicDistance[x_, size_: 1] := Module[{xm, ym, zm},
xm = {Outer[EuclideanDistance, x[[;; , 1]], x[[;; , 1]]],
size - Outer[EuclideanDistance, x[[;; , 1]], x[[;; , 1]]]};
xm = Map[Min[Abs[#]] &, Transpose[xm, {3, 1, 2}], {2}];
ym = {Outer[EuclideanDistance, x[[;; , 2]], x[[;; , 2]]],
size - Outer[EuclideanDistance, x[[;; , 2]], x[[;; , 2]]]};
ym = Map[Min[Abs[#]] &, Transpose[ym, {3, 1, 2}], {2}];
zm = {Outer[EuclideanDistance, x[[;; , 3]], x[[;; , 3]]],
size - Outer[EuclideanDistance, x[[;; , 3]], x[[;; , 3]]]};
zm = Map[Min[Abs[#]] &, Transpose[zm, {3, 1, 2}], {2}];
Sqrt[xm^2 + ym^2 + zm^2]
]


As a side question, is it possible to carry out such computation within mathematica using e.g. octree so that it could be scaled up to millions on points?

How can it be less than O(n^2) in memory or speed? That's the information-theoretic lower bound. – Daniel Lichtblau – 2012-06-04T20:48:28.800

Well I am not asking for the impossible ;-). Are you implying the code I wrote is already at the theoretical lower bound? – chris – 2012-06-04T20:53:02.887

2"Millions of points" would produce a matrix with tens of trillions of (nonzero) entries. I don't see scaling up to that in Mathematica any time soon :-). – whuber – 2012-06-05T15:22:03.127

@chris Hi Chris, great question and very instructive answers. Quick question about a slight extension, how could we e.g. extend Heike's method to non-cubic boxes? (say a cubic box where for instance one dimension is elongated.) – None – 2019-05-02T08:31:02.503

@user929304 it seems you just need to replace Norm by something anisotropic ? As in Norm1 := Sqrt[3] Norm[#/{1, 2}] & ?? I have not checked. – chris – 2019-05-02T09:01:54.123

13

I'm assuming here that x is a list of points between which you want to calculate the distance. If so, then I think your code can be condensed to something like

PeriodicDistance[x_, size_: 1] := Outer[Norm@Mod[#2 - #1, size, -size/2] &, x, x, 1]


Edit

A faster version of the code above is something like

PeriodicDistance3[x_, size_: 1] :=
Map[Norm, Mod[Outer[Subtract, x, x, 1], size, -size/2], {-2}]


On my system this is about as efficient as the code in the original question, whereas the previous version is a factor 2 to 3 slower.

I sort of knew someone would come with something shorter! Impressive. I didn't realize Mod could take two arguments. The beauty of your solution is that it works in any dimensions as well. – chris – 2012-06-04T20:53:50.063

It's shorter, but it seems to be slower as well. For 200 points it's three times as slow as the original code. I'll see if I can speed things up. – Heike – 2012-06-04T20:57:54.233

2

@Heike The main source of inefficiency is that Outer is unable to utilize Listability of built-in functions. In a very similar situation here, I used a different method utilizing that. Unfortunately, Norm also suffers from that drawback, which is why I did not use it in that answer.

– Leonid Shifrin – 2012-06-04T21:32:45.607

10

Very late to the party, but I'll show a method that's faster than anything posted so far and will be hard to beat.

First let's define our PeriodicDistance:

PeriodicDistance = Compile[{{x, _Real, 1}, {y, _Real, 2}, {size, _Real}},
Sqrt@Total[((x - y) - size*Round[(x - y)/size])^2],
CompilationTarget :> "C", RuntimeOptions -> "Speed"];


Now our distance matrix:

PeriodicDistanceMatrix[dat_?MatrixQ, size_Real:1.0] :=
With[{tr = Transpose[dat]}, Map[PeriodicDistance[dat[[#]], tr, size] &, Range@Length@dat]]


Some Timings:

data = RandomReal[1, {2000, 3}];

PeriodicDistanceMatrix[data]; // AbsoluteTiming

(* 0.265632 *)

(* Heike's faster version *)

PeriodicDistanceB[data]; //AbsoluteTiming

(* 2.599099 *)

PeriodicDistance[data]; // AbsoluteTiming

(* 6.246279 *)

(* O.P.'s version *)

PeriodicDistance2[data]; // AbsoluteTiming

(* 10.026386 *)


And yes, all the outputs are identical.

Edit

Needs["GeneralUtilities"] (* for v10 *)

BenchmarkPlot[{RunnyKine, Heike, JM}, RandomReal[1, {#, 3}] &, 2^Range[6, 14]]


Thanks! will look at it after flying to canada :-) – chris – 2014-05-13T20:53:06.377

1@chris, You're Welcome and fly safe :). I didn't think anyone was going to look at this since it's been almost 2 years since you posed this question. – RunnyKine – 2014-05-13T21:17:43.063

I know this is slightly unrelated, because the OP asked for distances on a periodic cube, but how could you implement your method for a periodic rectangle? As in, one side goes from [0,Pi] and the other [0, 2Pi] – Jesse – 2014-08-08T16:38:07.317

Is there a cleaner way to do it than this?

PeriodicDistance[v1_, v2_, bound_] := Norm[bound[[#]]/2 - Abs[bound[[#]]/2 - Abs[v1[[#]] - v2[[#]]]] & /@ {1, 2}] – Jesse – 2014-08-08T17:19:01.427

@Jesse, do you want it Compiled like in my answer? – RunnyKine – 2014-08-08T17:33:27.460

@RunnyKine I was more looking for a sexier way to access the arguments than using [[#]] and mapping to {1,2}. – Jesse – 2014-08-11T23:06:14.193

9

You want something like

PeriodicDistance[pts_, size_: 1] :=
Outer[EuclideanDistance, size*FractionalPart[pts/size],
size*FractionalPart[pts/size], 1]


But what you did is a bit different in terms of how you wrap things periodically. If it is really as you intended, replace size everywhere by 2*size.

PeriodicDistanceB[pts_, size_: 1] :=
Outer[EuclideanDistance, 2*size*FractionalPart[pts/(2*size)],
2*size*FractionalPart[pts/(2*size)], 1]


1It seems (?) neither give the answer I want: e.g. PeriodicDistanceB[{{0.1, 0, 0}, {0.9, 0, 0}}] should return {{0., 0.2}, {0.2, 0.}} not {{0., 0.8}, {0.8, 0.}} – chris – 2012-06-05T06:18:42.160

Ah. I see the issue. Roughly, you are using the closest distance possible, allowing for taking points to be in neighboring cubes rather than the same one if that happens to make them closer. At this point I'd say use one of the other proposed methods, assuming they do what you want. – Daniel Lichtblau – 2012-06-05T15:34:45.420

5

I've had to do this operation once when I was doing research on Voronoi diagrams. Heike's method is nice. Here is one possible alternative, making use of a fold-over technique to account for periodicity:

cPerDist = Compile[{{v1, _Real, 1}, {v2, _Real, 1}, {size, _Real}},
Norm[size/2 - Abs[size/2 - Abs[v1 - v2]]]];

PeriodicDistance[x_?MatrixQ, size_:1] := Outer[cPerDist[#1, #2, size] &, x, x, 1]


(You might want to look at a plot of the function $\frac12-\left|\frac12-\left|x\right|\right|$ to get an idea of how the fold-over works.)

Here, we generate and use a compiled function, since this will seem to be used many times on inexact numbers. If you need exact expressions for the lengths, it should be straightforward to write an uncompiled version of cPerDist.

Try it out:

PeriodicDistance[{{0.1, 0, 0}, {0.9, 0, 0}}]
{{0., 0.2}, {0.2, 0.}}

pts = BlockRandom[SeedRandom[42, Method -> "Legacy"]; RandomReal[1, {12, 3}]];

PeriodicDistance[pts]
{{0.,0.369061,0.557349,0.273739,0.580937,0.607056,0.630916,0.414241,0.149751,0.214448,0.577028,0.467806},
{0.369061,0.,0.431869,0.636656,0.611335,0.628004,0.641477,0.584125,0.465563,0.34617,0.454007,0.574909},
{0.557349,0.431869,0.,0.653288,0.526924,0.323731,0.230079,0.450194,0.538027,0.589348,0.581119,0.669271},
{0.273739,0.636656,0.653288,0.,0.521814,0.433284,0.539068,0.278055,0.201587,0.392315,0.607604,0.52001},
{0.580937,0.611335,0.526924,0.521814,0.,0.333371,0.527569,0.529642,0.484959,0.691898,0.237093,0.474143},
{0.607056,0.628004,0.323731,0.433284,0.333371,0.,0.308742,0.343472,0.590752,0.69726,0.468426,0.596956},
{0.630916,0.641477,0.230079,0.539068,0.527569,0.308742,0.,0.486608,0.564332,0.543762,0.589842,0.62605},
{0.414241,0.584125,0.450194,0.278055,0.529642,0.343472,0.486608,0.,0.396313,0.377811,0.538572,0.380586},
{0.149751,0.465563,0.538027,0.201587,0.484959,0.590752,0.564332,0.396313,0.,0.307534,0.602186,0.56417},
{0.214448,0.34617,0.589348,0.392315,0.691898,0.69726,0.543762,0.377811,0.307534,0.,0.600932,0.352568},
{0.577028,0.454007,0.581119,0.607604,0.237093,0.468426,0.589842,0.538572,0.602186,0.600932,0.,0.276919},
{0.467806,0.574909,0.669271,0.52001,0.474143,0.596956,0.62605,0.380586,0.56417,0.352568,0.276919,0.}}


Some limited tests I did indicate that the method is at least as fast as Heike's.

Thanks. For 2000 3D points I get PeriodicDistance[pts]; // Timing {9.88636,Null} whereas Heike's code yields PeriodicDistance3[pts]; // Timing {5.40553,Null}` Did you do 3D Voronoi Triangulation in mathematica? – chris – 2012-06-05T15:15:12.097

Hmm, they took about the same length of time (give or take a hundredths of a second) on my box. Oh well. I was looking at both two and three dimensions, so I wrote my routines to work for arbitrary dimensions. – J. M.'s ennui – 2012-06-05T15:22:13.440