10
$\begingroup$

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?

Thanks for your help

$\endgroup$
5
  • 1
    $\begingroup$ How can it be less than O(n^2) in memory or speed? That's the information-theoretic lower bound. $\endgroup$ Jun 4 '12 at 20:48
  • $\begingroup$ Well I am not asking for the impossible ;-). Are you implying the code I wrote is already at the theoretical lower bound? $\endgroup$
    – chris
    Jun 4 '12 at 20:53
  • 2
    $\begingroup$ "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 :-). $\endgroup$
    – whuber
    Jun 5 '12 at 15:22
  • $\begingroup$ @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.) $\endgroup$
    – user52181
    May 2 '19 at 8:31
  • $\begingroup$ @user929304 it seems you just need to replace Norm by something anisotropic ? As in Norm1 := Sqrt[3] Norm[#/{1, 2}] & ?? I have not checked. $\endgroup$
    – chris
    May 2 '19 at 9:01
13
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – chris
    Jun 4 '12 at 20:53
  • $\begingroup$ 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. $\endgroup$
    – Heike
    Jun 4 '12 at 20:57
  • 2
    $\begingroup$ @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. $\endgroup$ Jun 4 '12 at 21:32
10
$\begingroup$

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

(* This answer *)

PeriodicDistanceMatrix[data]; // AbsoluteTiming

(* 0.265632 *)

(* Heike's faster version *)

PeriodicDistanceB[data]; //AbsoluteTiming

(* 2.599099 *)

(* J.M.'s answer *)

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

Mathematica graphics

$\endgroup$
6
  • $\begingroup$ Thanks! will look at it after flying to canada :-) $\endgroup$
    – chris
    May 13 '14 at 20:53
  • 1
    $\begingroup$ @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. $\endgroup$
    – RunnyKine
    May 13 '14 at 21:17
  • $\begingroup$ 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] $\endgroup$
    – Jesse
    Aug 8 '14 at 16:38
  • $\begingroup$ 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}] $\endgroup$
    – Jesse
    Aug 8 '14 at 17:19
  • $\begingroup$ @Jesse, do you want it Compiled like in my answer? $\endgroup$
    – RunnyKine
    Aug 8 '14 at 17:33
9
$\begingroup$

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]
$\endgroup$
2
  • 1
    $\begingroup$ It 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.}} $\endgroup$
    – chris
    Jun 5 '12 at 6:18
  • $\begingroup$ 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. $\endgroup$ Jun 5 '12 at 15:34
5
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – chris
    Jun 5 '12 at 15:15
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Jun 5 '12 at 15:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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