Fast application of InterpolatingFunction over transformed coordinates

I have an InterpolatingFunction that takes $$(x, y)$$ coordinates as input and I would like to apply it to a series of data that is rotated by $$\frac{\pi}{4}$$ into $$(\frac{1}{\sqrt{2}}(x+y), \frac{1}{\sqrt{2}}(x-y)))$$ coordinates.

Here's a sample function:

FunctionInterpolation[
Sin[θ]*Cos[ϕ],
{θ, 0, 2 π},
{ϕ, 0, π}
]


Furthermore, if the value is off the interpolation grid I don't want to extrapolate, but would rather provide a default value. In general my use case has very steep walls and so I basically make this infinity.

Now, I've written code to do this quickly using Compile, but I'd like to generalize this and potentially do better. The Compile code is long (and still has a "MainEvaluate" call) so I'll post it as an answer at some later time.

So what's the fastest way to apply that InterpolatingFunction to a transformed set of coordinates (we can assume the transformation can be compiled) with a default value off the standard grid?

Setting up the function and the default value. Notice that I force Mathematica to use machine precision numbers in f.

f = FunctionInterpolation[Sin[θ]*Cos[ϕ], {θ, 0., 2. π}, {ϕ, 0., 1. π}];
default = 0.;


Generating several example points.

n = 100000;
x = RandomReal[{0, 2 Pi}, {n}];
y = RandomReal[{0, Pi}, {n}];
{x, y} = 1/Sqrt[2.] {x + y, x - y};


Now, let's find the points within the domain of f:

inregionQ = UnitStep[x] UnitStep[2 Pi - x] UnitStep[y] UnitStep[Pi - y];
idx = RandomPrivatePositionsOf[inregionQ, 1];


And let's apply f only to these points:

result = ConstantArray[default, Length[x]];
result[[idx]] = f[x[[idx]], y[[idx]]]; // AbsoluteTiming // First


0.741009

It is still not super fast but the reason is that InterpolatingFunctions are rather slow. As we may use a regular tensor grid for f, applying it to many points at once could be way faster with a more appropriate binning and parallelization strategy...

Just for comparison: This is the timing of the evalation of the original function:

result2 = ConstantArray[default, Length[x]];
result2[[idx]] = Sin[x[[idx]]] Cos[y[[idx]]]; // AbsoluteTiming // First


0.001445

It is more than 500 times faster. This tells us (i) that the function to interpolate has to be very complex in order to take advantage of FunctionInterpolation and (ii) that there must be a lot potential to improve the performance of InterpolationFunction.

• Interesting approach. Never seen RandomPrivatePositionsOf but it looks useful, esp. with UnitStep for picking indices. – b3m2a1 Dec 12 '18 at 9:52
• Yeah, Carl Woll used it a couple of times. I really like it. It is much faster than Position for finding integers and it is also faster than Pick[Range[Length[inregionQ]],inregionQ,1]. – Henrik Schumacher Dec 12 '18 at 9:54
• Ah, I see. Yeah. Dunno. It is a black box for me and I only use it for integer lists anyway. I suspect a compiled function behind it, so I do not expect much flexibility from it. – Henrik Schumacher Dec 12 '18 at 10:01
• Per your comparison, I have no analytic form for my original thing and probably no hope ever to have an analytic form. I could try some fits but what fit to use it a mystery as it's a high dimensional system. Hence my hope to deal with interpolation really nicely. – b3m2a1 Dec 12 '18 at 10:34
• Yeah. I did not mean to poke you about that. My hope was more towards a random WRI developer stumbling upon this and feeling an itching... ;) – Henrik Schumacher Dec 12 '18 at 11:04

The following produces an InterpolatingFunction that does the desired extrapolation. It's about the same speed as @Henrik's, but I think it's convenient.

InternalInheritedBlock[{Interpolation},
Unprotect[Interpolation];
Interpolation[args___] /; ! TrueQ[$$in] := Block[{$$in = True},
With[{default = 0.},
Interpolation[args, "ExtrapolationHandler" -> {default &, "WarningMessage" -> False}]
]];
Protect[Interpolation];
ifn = FunctionInterpolation[
Sin[\[Theta]]*Cos[\[Phi]], {\[Theta], 0., 2. \[Pi]}, {\[Phi], 0., N@\[Pi]}]
]


Henrik's example with SeedRandom[1] (Length@idx --> 65673), followed by the ifn above on Henrik's input and the full x,y input:

result[[idx]] = f[x[[idx]], y[[idx]]]; // RepeatedTiming // First
(*  0.801  *)

result[[idx]] = ifn[x[[idx]], y[[idx]]]; // RepeatedTiming // First
(*  0.753  *)

result[[idx]] = ifn[x, y]; // RepeatedTiming // First
(*  0.814  *)


Here's my version of this. I use it internally in a discretized function object to provide a nice compiled functional form. One other thing this does is allow you to set some coordinates to a constant value to take slices through the function:

iInterpCompile[interp_, coords_, preProcessor_, def_, min_] :=

Module[{potInterp = interp, dom = interp[[1]], j = 1},
Block[{crds, pt, gps},
Module[
{
getPot,
varBlock,
varBlock2,
preCoords
},
getPot[a_] := Apply[potInterp, a];
preCoords =
Quiet@
preProcessor@
Map[
If[# === Automatic,
CompileGetElement[crds, j++],
#
] &,
coords
];
j = 1;
(* bounds checking for test *)
With[
{
test =
With[
{
tests =
#[[1]] <=
If[NumericQ[#2],
j++; #2,
CompileGetElement[pt, j++]
] <= #[[2]] &,
{dom, preCoords}
]
},
If[MemberQ[tests, False],
False,
And @@ DeleteCases[tests, True]
]
],
crdSpec = preCoords,
means = Mean /@ dom,
fn =
Compile[
{{crds, _Real, 1}},
Evaluate@preCoords,
RuntimeAttributes -> {Listable}
]
},
Compile @@

Hold[(* True Compile body but with Hold for code injection *)
{{gps, _Real, 2}},
Module[
{
pts = fn@gps,
ongrid = Table[0, {i, Length@gps}],
intres,
midpt = means,
interpvals,
interpcounter,
interppts,
i = 1,
barrier = def,
minimum = min
},
interppts = Table[midpt, {i, Length@pts}];
(* Find points in domain *)
interpcounter = 0;
Do[
If[test,
ongrid[[i]] = 1;
interppts[[++interpcounter]] = pt
];
i++,
{pt, pts}
];
(* Apply InterpolatingFunction only once (one MainEval call) *)

If[interpcounter > 0,
intres = interppts[[;; interpcounter]];
interpvals = getPot[Transpose@intres] - minimum;
(* Pick points that are in domain *)

interpcounter = 0;
Map[
If[# == 1, interpvals[[++interpcounter]], barrier] &,
ongrid
],
Table[barrier, {i, Length@gps}]
]
],
{{getPot[_], _Real, 1}},
RuntimeOptions -> {
"EvaluateSymbolically" -> False,
"WarningMessages" -> False
},
CompilationOptions -> {"InlineExternalDefinitions" -> True},
RuntimeAttributes -> {Listable}
]
]
]
]
]

Options[InterpCompile] =
{
"CoordinateTransformation" -> Identity,
"Minimum" -> 0.,
"Default" -> 0.,
"Mode" -> "Vector"
};
InterpCompile[
interp_,
coordSpec : {(Automatic | _?NumericQ) ..} | Automatic : Automatic,
ops : OptionsPattern[]
] :=
If[OptionValue["Mode"] === "Vector",
iInterpCompile[
interp,
Replace[
coordSpec,
Automatic :> ConstantArray[Automatic, Length@interp[[1]]]
],
OptionValue["CoordinateTransformation"],
OptionValue["Default"],
OptionValue["Minimum"]
]
]


Here's a sample usage. First set up the compiled function:

ga =
CoordinateBoundsArray[
{{0., 2. π}, {0., 1. π}},
{π/24., π/24.}
];
fn =
Interpolation@
Flatten[
Join[
ga,
Map[
{Sin[#[[1]]]*Cos[#[[2]]]} &,
ga,
{2}
],
3
],
1
];
cf =
InterpCompile[fn,
"CoordinateTransformation" ->
(1/Sqrt[2.]*
{
#[[2]] + #[[1]], #[[1]] - #[[2]]
}
&
)];


Then test the timing:

cf@testGrid; // RepeatedTiming

{0.47, Null}


And plot the result:

cf@testGrid // ListPlot3D[#, PlotRange -> All] &


It's only barely faster than Henrik's method and multiple orders of magnitude more complex, but on the other hand it's pretty much fully general, which is convenient.

We can also use it to take slices through the interpolation:

With[
{
cf2 =
InterpCompile[fn,
{Automatic, #},
"CoordinateTransformation" ->
(1/Sqrt[2.]*
{
#[[2]] + #[[1]], #[[1]] - #[[2]]
}
&
)]
},
Transpose@{#, cf2@Transpose@{#}} &@Range[-8, 8, .001] //
ListLinePlot[#, ImageSize -> 250] &
] & /@ Range[-3, 3] // Partition[#, 2] & // GraphicsGrid
`