Krdytkyu

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 = Random`Private`PositionsOf[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.

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, 

          Compile`GetElement[crds, j++], 

          #

          ] &,

        coords

        ];

    j = 1;

    (* bounds checking for test *)

    With[

     {

      test =

       With[

        {

         tests =

          MapThread[

           #[[1]] <=

             If[NumericQ[#2], 

              j++; #2, 

              Compile`GetElement[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

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.

Internal`InheritedBlock[{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  *)