# Listing over multiple inputs (tensor rank > 1) in a compiled function

I'm trying to optimize the solution given by Michael E2 to my previous problem (Efficiently determining if 3D points are within a surface composed of polygons) which I figured would be relatively straight forward to do with Compile but I'm having trouble listing over two inputs to a compiled function.

Starting with Michael's code:

data = Flatten[
Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]},
{x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];

plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None];

polygonCoords =
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}];

g = Graphics3D[{Opacity[.09], EdgeForm[Opacity[.3]],
Polygon[#,
VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400, Axes -> True];

side[{P_, Q_, R_, ___}, X_] := Det@Differences[{X, P, Q, R}];
insideQ[polyhedron_, point_] := And @@ Positive[side[#, point] & /@ polyhedron];

points = RandomReal[{-1.5, 1.5}, {1000, 3}];
insideQ[polygonCoords, #] & /@ points; // AbsoluteTiming


Too slow: 0.5 sec for 1000 points and I have millions.

I inlined the formula for the det and replaced And + Positive with Min + UnitStep:

sideC = Compile[{{poly, _Real, 2}, {X, _Real, 1}},
Block[{p1 = poly[[1, 1]], p2 = poly[[1, 2]], p3 = poly[[1, 3]],
q1 = poly[[2, 1]], q2 = poly[[2, 2]], q3 = poly[[2, 3]],
r1 = poly[[3, 1]], r2 = poly[[3, 2]], r3 = poly[[3, 3]],
x1 = X[[1]], x2 = X[[2]], x3 = X[[3]]},
UnitStep[-p3 q2 r1 + p2 q3 r1 + p3 q1 r2 - p1 q3 r2 - p2 q1 r3 +
p1 q2 r3 + p3 q2 x1 - p2 q3 x1 - p3 r2 x1 + q3 r2 x1 +
p2 r3 x1 - q2 r3 x1 - p3 q1 x2 + p1 q3 x2 + p3 r1 x2 -
q3 r1 x2 - p1 r3 x2 + q1 r3 x2 + p2 q1 x3 - p1 q2 x3 -
p2 r1 x3 + q2 r1 x3 + p1 r2 x3 - q1 r2 x3]
], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];

insideQC[polyhedron_, point_] := Min@sideC[polygonCoords[[All, ;; 3]], point]


This lists over all the polygons. Lets test it,

insideQC[polygonCoords, #] & /@ points; // AbsoluteTiming


Now this is over 10x faster, 0.04 seconds. Compiling this again, because of the CompiledFunctionCall provides only a tiny speed improvment.

I think it should be faster to list over both the list of polygon coordinates and the points. I suspect this question is related: Threading a compiled function over multiple arguments of different lengths, specifically ruebenko's last answer which works fine on MMA 9 on Win7 x64, It should be possible?

Why doesn't this work:

sideC[polygonCoords, points]


Shouldn't this list over both inputs? How can I make this work?

I tried to force inline-ing of the compiled function definition but it doesn't appear to be supported?

insideQC =
Compile[{{polyhedron, _Real, 3}, {points, _Real, 2}},
Min@sideC[polyhedron, #] & /@ points, CompilationTarget -> "C",
RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True},
RuntimeAttributes -> {Listable}, Parallelization -> True];


Compile::cfinll: "The CompiledFunction could not be inlined because its use requires threading with the Listable runtime attribute."

-
Have you perhaps found a solution since you've asked this question? – shrx Jun 2 '15 at 10:02