3
$\begingroup$

I need to increase the speed of my code. I have a set of $n$ points in a square box of length $boxLength$.

n = 60;
boxLength = 20;
positions = Table[{Random[], Random[]}*boxLength, {i, 1, n}];

The module

RelativePosition[i_, j_] := 
 Module[{ri, rj, li, lj  , dr, Rij, sizeRij, untiRij},
   ri = positions[[i]];
   rj = positions [[j]];
   li = lmoment [[i]] ;
   lj = lmoment [[j]] ;

   dr = ri - rj;
   dr = dr - Round[dr/boxLength]*boxLength;
   Rij = dr + li - lj;  

   sizeRij = 1.*Norm[Rij];
   untiRij = Rij/sizeRij;

   Return [{sizeRij, untiRij }]
   ]

gives me the relative distance between two points. I used the following Compiled Function instead

RelativePosition2 = Compile[{{i, _Integer}, {j, _Integer}},

  Module[{ii = i, jj = j, ri, rj, li, lj  , dr, Rij, sizeRij, untiRij},

   ri = positions[[ii]];
   rj = positions [[jj]];
   li = lmoment [[ii]] ;
   lj = lmoment [[jj]] ;

   dr = ri - rj;
   dr = dr - Round[dr/boxLength]*boxLength;
   Rij = dr + li - lj;  
   sizeRij = 1.*Norm[Rij];
   untiRij = Rij/sizeRij;

   Return [{sizeRij, untiRij }]
   ]
  ]

RelativePosition[1, 2]

It gives me a result but also tells me

CompiledFunction::cfex: Could not complete external evaluation at instruction 31; proceeding with uncompiled evaluation.

What does the message mean? I cannot figure out what is wrong. Does it mean the Mathematica is not compiling the function?

Moreover, the compiled function is slower. The code

Do[RelativePosition[1, 2], {i, 10000}]; // AbsoluteTiming

gives

{0.27862, Null}

white

Do[RelativePosition2[1, 2], {i, 10000}]; // AbsoluteTiming

gives

{0.367405, Null}
$\endgroup$
1
  • 1
    $\begingroup$ Norm is a super general function so it does not compile. Use CompilePrint (read about it elsewhere here) to get a sense for what is preventing your function from truly compiling. $\endgroup$
    – b3m2a1
    Commented Jul 16, 2018 at 5:45

2 Answers 2

4
$\begingroup$
Parallelize[DistanceMatrix[positions]]

For your example, Timing gives 0.002269 seconds, faster than your code.

If you want a list of the difference vectors:

Outer[Subtract, positions, positions, 1]

To get the normalized (unit-length) vector between each pair of points:

Map[Normalize, Outer[Subtract, positions, positions, 1], {2}]
$\endgroup$
7
  • $\begingroup$ This is great, but I need the unit vector between the points too. $\endgroup$
    – Fluid
    Commented Jul 15, 2018 at 22:20
  • $\begingroup$ Also, DistanceMatrix does not take into account the periodic boundary condition in my code. $\endgroup$
    – Fluid
    Commented Jul 15, 2018 at 22:33
  • $\begingroup$ Where in your problem do you state that the problem is periodic? $\endgroup$ Commented Jul 15, 2018 at 22:41
  • 1
    $\begingroup$ Outer[Plus, positions, -positions, 1] faster than Outer[Subtract, positions, positions, 1] :) $\endgroup$
    – matrix42
    Commented Jul 16, 2018 at 3:16
  • 1
    $\begingroup$ @mathe: Interesting... I wonder why it is faster... $\endgroup$ Commented Jul 16, 2018 at 3:21
2
$\begingroup$

Have look at CompiledFunctionTools`CompilePrint@RelativePosition. Do you see the calls to MainEvaluate? They are due to the global variables positions, lmoment, and boxLength which are unknown to the compiled function. As a quick guess, I would suggest the following:

RelativePosition = Compile[{
   {positions, _Real, 2}, {lmoment, _Real, 2}, {boxLength, _Real},
   {i, _Integer}, {j, _Integer}
   },
  Module[{ii = i, jj = j, ri, rj, li, lj, dr, Rij, sizeRij, untiRij},
   ri = positions[[ii]];
   rj = positions[[jj]];
   li = lmoment[[ii]];
   lj = lmoment[[jj]];
   dr = ri - rj;
   dr = dr - Round[dr/boxLength]*boxLength;
   Rij = dr + li - lj;
   sizeRij = 1.*Norm[Rij];
   untiRij = Rij/sizeRij;
   Return[{sizeRij, untiRij}]]
  ]

However, the return value has a mixed type (it is not an array). That may cause trouble when you run it.

$\endgroup$
2
  • $\begingroup$ Thanks, I used your function and ran the code RelativePosition[positions, lmoment, boxLength, 1, 2] and it gave me: CompiledFunction::cfta: Argument {0.3,0.3,0.3,0.3,0.3} at position 2 should be a rank 2 tensor of machine-size real numbers. {1.22977, {0.820146, -0.572154}} $\endgroup$
    – Fluid
    Commented Jul 15, 2018 at 22:25
  • 1
    $\begingroup$ Since you did not provide all relevant information, that's what had to happen: I had to make some assumptions about the nature of positions, lmoment, and boxLength. Apparently, some of them were wrong. But who is to blame for that? $\endgroup$ Commented Jul 16, 2018 at 9:05

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.