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

Improve this question
$\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$
    – Daniel Lichtblau
    Commented Jun 4, 2012 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
    Commented Jun 4, 2012 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
    Commented Jun 5, 2012 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
    Commented May 2, 2019 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
    Commented May 2, 2019 at 9:01

4 Answers 4

Reset to default
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.

Improve this answer
$\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
    Commented Jun 4, 2012 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
    Commented Jun 4, 2012 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$
    – Leonid Shifrin
    Commented Jun 4, 2012 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

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks! will look at it after flying to canada :-) $\endgroup$
    – chris
    Commented May 13, 2014 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
    Commented May 13, 2014 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
    Commented Aug 8, 2014 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
    Commented Aug 8, 2014 at 17:19
  • $\begingroup$ @Jesse, do you want it Compiled like in my answer? $\endgroup$
    – RunnyKine
    Commented Aug 8, 2014 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]
Improve this answer
$\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
    Commented Jun 5, 2012 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$
    – Daniel Lichtblau
    Commented Jun 5, 2012 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.

Improve this answer
$\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
    Commented Jun 5, 2012 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 missing motivation ♦
    Commented Jun 5, 2012 at 15:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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