# A Compiled Function for the distance between a set of point

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}

• 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. – b3m2a1 Jul 16 '18 at 5:45

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

• This is great, but I need the unit vector between the points too. – Fluid Jul 15 '18 at 22:20
• Also, DistanceMatrix does not take into account the periodic boundary condition in my code. – Fluid Jul 15 '18 at 22:33
• Where in your problem do you state that the problem is periodic? – David G. Stork Jul 15 '18 at 22:41
• Outer[Plus, positions, -positions, 1] faster than Outer[Subtract, positions, positions, 1] :) – matrix89 Jul 16 '18 at 3:16
• @mathe: Interesting... I wonder why it is faster... – David G. Stork Jul 16 '18 at 3:21

Have look at CompiledFunctionToolsCompilePrint@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.

• 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}} – Fluid Jul 15 '18 at 22:25
• 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? – Henrik Schumacher Jul 16 '18 at 9:05