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.

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?

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

share|improve this question
    
How can it be less than O(n^2) in memory or speed? That's the information-theoretic lower bound. –  Daniel Lichtblau Jun 4 '12 at 20:48
    
Well I am not asking for the impossible ;-). Are you implying the code I wrote is already at the theoretical lower bound? –  chris Jun 4 '12 at 20:53
    
"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 Jun 5 '12 at 15:22

4 Answers 4

up vote 11 down vote accepted

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.

share|improve this answer
    
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 Jun 4 '12 at 20:53
    
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 Jun 4 '12 at 20:57
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 Jun 4 '12 at 21:32

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]
share|improve this answer
1  
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.}} –  chris Jun 5 '12 at 6:18
    
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 Jun 5 '12 at 15:34

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.

share|improve this answer
    
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 Jun 5 '12 at 15:15
    
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. Jun 5 '12 at 15:22

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.

share|improve this answer
    
Thanks! will look at it after flying to canada :-) –  chris May 13 at 20:53
    
@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 May 13 at 21:17

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.