Using Indexed variables in Compile?

Question

This follows a: useful post about Compile from about 10 years ago

In that post there is a comment:

You can not really use indexed variables in Compile, although it may appear that you can.

For example,

cf=Compile[{{p, _Real, 3}},Indexed[p, {1, 2, 3}] Indexed[p, {3, 2, 1}]]

will compile and work just fine.

With[{rr = RandomReal[{0, 1}, {3, 3, 3}]},
 {rr[[1, 2, 3]] rr[[3, 2, 1]], cf[rr]}
 ]

I seem to remember, but cannot find, a utility in Needs["CompiledFunctionTools`"] that assists in using Indexed variables in Compile.

Does anyone know if such utility exists; if so what is it, and does it make extraction from a tensor more efficient?

Needs["CompiledFunctionTools`"]
?CompiledFunctionTools`*

doesn't show anything like what I am looking for.

Answer 1

I'll sum up what I've learned.

The approach of using Indexed variables like so:

cf=Compile[{{p, _Real, 3}},Indexed[p, {1, 2, 3}] Indexed[p, {3, 2, 1}]]

is fine and also efficient.

Here is a usage scenario:

Suppose there is a very large set of line segments and you wish to find all intersecting pairs. It needs to be fast.

Suppose you want to call the function like so:

segmentsIntersect[{{{x1B, y1B}, {x1E, y1E}}, {{x2B, y2B}, {x2E, y2E}}}]

where {{{x1B, y1B}, {x1E, y1E}}, {{x2B, y2B}, {x2E, y2E}}} is a list of two segments, each segment being a list of its endpoints and the function returns True if they intersect.

Compute the intersection condition by using a parametric form of the line segments:

I use [FormalP] here because I find it useful in the body of Compile or Function.

s1B = {x1B, y1B} = {Indexed[\[FormalP], {1, 1, 1}], 
    Indexed[\[FormalP], {1, 1, 2}]};
s1E = {x1E, y1E} = {Indexed[\[FormalP], {1, 2, 1}], 
    Indexed[\[FormalP], {1, 2, 2}]};
s2B = {x2B, y2B} = {Indexed[\[FormalP], {2, 1, 1}], 
    Indexed[\[FormalP], {2, 1, 2}]};
s2E = {x2E, y2E} = {Indexed[\[FormalP], {2, 2, 1}], 
    Indexed[\[FormalP], {2, 2, 2}]};

seg1 = s1B + v (s1E - s1B);
seg2 = s2B + u (s2E - s2B);

The condition for segment intersection is that the parameters u and v are both in [0,1].

sol = {{uSol, vSol}} = FullSimplify[
   SolveValues[Thread[seg1 ==  seg2] , {u, v}]

intersectionCondition = And @@ ((0 <= # <= 1) & /@ sol[[1]])

Then the symbolic condition can be used in Compile:

segmentsIntersect = Compile[{{\[FormalP], _Real, 3}}, Evaluate[intersectionCondition]]

Evaluate is needed here because Compile has the attribute HoldAll

For clarity, I've neglected failure modes where a determinant-divisor is zero and the segments are collinear.

This can also be improved by filtering with the condition that the bounding boxes of each segment intersect.

Compile is fairly efficient at recognizing terms that appear twice and writing code that computes them only once. However, in this case, Compile misses the determinant-divisor. The compiled code can be improved by by forcing the computation of the determinant once.

For a large set of segment pairs, it is useful to make the compiled function listable:

segmentsIntersect = 
 Compile[{{\[FormalP], _Real, 3}}, Evaluate[intersectionCondition],
  RuntimeAttributes -> {Listable}]

Then,

segmentPairs = RandomReal[{0, 1}, {10, 2, 2, 2}];
segmentsIntersect[segmentPairs]