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 Jul 16 '18 at 5:45
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 Jul 15 '18 at 22:20
  • $\begingroup$ Also, DistanceMatrix does not take into account the periodic boundary condition in my code. $\endgroup$ – Fluid Jul 15 '18 at 22:33
  • $\begingroup$ Where in your problem do you state that the problem is periodic? $\endgroup$ – David G. Stork Jul 15 '18 at 22:41
  • 1
    $\begingroup$ Outer[Plus, positions, -positions, 1] faster than Outer[Subtract, positions, positions, 1] :) $\endgroup$ – matrix89 Jul 16 '18 at 3:16
  • 1
    $\begingroup$ @mathe: Interesting... I wonder why it is faster... $\endgroup$ – David G. Stork Jul 16 '18 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 Jul 15 '18 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$ – Henrik Schumacher Jul 16 '18 at 9:05

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.