# What can be improved in the following custom interpolation?

Consider some test data: a grid of x, y, and values of some function on this grid:

func[x_, y_] =
Exp[-((y - 5)/(10.*Cos[x]^4))] ((Sin[x] Cos[x]^5)/y^2 + Sin[x]^2/(
y^6 + 1) + Sin[x]*Cos[x]^4);
gridx = Table[10^x, {x, -5., Log10[0.2], (Log10[0.2] + 5.)/30}];
gridy = Table[10^
x, {x, Log10[0.05], Log10[150], (Log10[150] - Log10[0.05])/80}];
{xminmax, yminmax} = MinMax[#] & /@ {gridx, gridy};
grid = Tuples[{gridx, gridy}];
distrvals = func @@@ grid;
data = Join[grid, Partition[distrvals, 1], 2];


I would like to "map" it onto random points belonging to the domain within the grid:

nev = 10^6;
{xrand, yrand} = RandomReal[#, nev] & /@ {xminmax, yminmax};
randompts = Join[Partition[xrand, 1], Partition[yrand, 1], 2];


The straightforward but slow way is to use built-in Interpolation:

intBuiltIn[x_, y_] =
Exp[Interpolation[Log[data], InterpolationOrder -> 1][Log[x],
Log[y]]]
vals1 = intBuiltIn @@@ randompts; // AbsoluteTiming


{4.15255,Null}

Here, I interpolated the logarithmized data to account for the structure of the initial grid and the exponential behavior of the initial function. This way is very slow, so I want to speed up it. I implement a bilinear compilable interpolation:

linint =
Compile[{{xlogrand, _Real}, {ylogrand, _Real}, {gridxlog, _Real,
1}, {gridylog, _Real,
1}, {nearestindexx, _Integer}, {nearestindexy, _Integer}, \
{distrlog, _Real, 1}},
Module[{xel, yel, x1, x2, y1, y2, z11, z21, z12,
z22, \[CapitalDelta], leny, indexlargerx, indexlargery},
(*Finding the index of the grids with the elements larged than the \
given random point*)
xel = CompileGetElement[gridxlog, nearestindexx];
indexlargerx = nearestindexx + UnitStep[xlogrand - xel];
yel = CompileGetElement[gridylog, nearestindexy];
indexlargery = nearestindexy + UnitStep[ylogrand - yel];
(*The two points x1, x2 with x1 < x < x2*)
x1 = CompileGetElement[gridxlog, indexlargerx - 1];
x2 = CompileGetElement[gridxlog, indexlargerx];
y1 = CompileGetElement[gridylog, indexlargery - 1];
y2 = CompileGetElement[gridylog, indexlargery];
(*Corresponding values of the function*)
leny = Length[gridylog];
z11 =
CompileGetElement[
distrlog, (indexlargerx - 2)*leny + indexlargery - 1];
z21 =
CompileGetElement[
distrlog, (indexlargerx - 1)*leny + indexlargery - 1];
z12 =
CompileGetElement[
distrlog, (indexlargerx - 2)*leny + indexlargery];
z22 =
CompileGetElement[
distrlog, (indexlargerx - 1)*leny + indexlargery];
\[CapitalDelta] = (x2 - x1) (y2 - y1);
Exp[(*{xlogrand,x1,x2,ylogrand,y1,y2,z11,z21,z12,
z22,*)((x2 - xlogrand)*(y2 - ylogrand))/\[CapitalDelta]*
z11 + ((xlogrand - x1)*(y2 -
ylogrand))/\[CapitalDelta] z21 + ((x2 - xlogrand) (ylogrand -
y1))/\[CapitalDelta]*
z12 + ((xlogrand - x1) (ylogrand - y1))/((x2 - x1) (y2 - y1))
z22(*}*)]
], CompilationTarget -> "C", RuntimeOptions -> "Speed",
RuntimeAttributes -> {Listable}, Parallelization -> True]


and then use Nearest:

(*Logarithmized grids and random points*)
{gridxlog, gridylog, xlogrand, ylogrand, distrlog} =
Log[#] & /@ {gridx, gridy, xrand, yrand,
distrvals}; // AbsoluteTiming
(*Finding the element of the grid which is nearest to the random \
points*)
nearestx = Nearest[gridxlog -> "Index"];
nearesty = Nearest[gridylog -> "Index"];
nearestindexx = nearestx[xlogrand] // Flatten; // AbsoluteTiming
nearestindexy = nearesty[ylogrand] // Flatten; // AbsoluteTiming


Once I create nearestx, nearesty (this has to be done only once), the approach is more than 100 times faster than with using the built-in interpolation:

vals2 = linint[xlogrand, ylogrand, gridxlog, gridylog, nearestindexx,
nearestindexy, distrlog]; // AbsoluteTiming
Abs[(vals2 - vals1)/vals1] // Max


{0.0328157,Null}

1.77041*10^-14

What may be improved in this algorithm? Can it be made faster?

• So your grid is structured: It is product grid. So you can use simple, independent binary searches for the x and y coordinates. This would allow you to circumvent the complicated handling of Nearest's output. Not sure whether that will be faster. Commented Dec 12, 2023 at 1:42
• @HenrikSchumacher : how do you think, if I have the initial data in the form {x,y,z,pdf[x,y,z]} and want to map it to pdf[X,xrand,yrand], where X is some point within the domain of definition of x and xrand, yrand are random points, would it be efficient to use first one of the methods from the question mathematica.stackexchange.com/questions/282637/… to map x,y,z,pdf[x,y,z] $\to$ y,z,pdf[X,y,z] and then the interpolation from this question? Commented Dec 12, 2023 at 22:27
• Phh.... don't know. I actually don't understand what you want to do. Sorry. When in doubt, code both approaches and measure the performance. Commented Dec 12, 2023 at 22:35