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], 
vals1 = intBuiltIn @@@ randompts; // AbsoluteTiming


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 = Compile`GetElement[gridxlog, nearestindexx];
   indexlargerx = nearestindexx + UnitStep[xlogrand - xel];
   yel = Compile`GetElement[gridylog, nearestindexy];
   indexlargery = nearestindexy + UnitStep[ylogrand - yel];
   (*The two points x1, x2 with x1 < x < x2*)
   x1 = Compile`GetElement[gridxlog, indexlargerx - 1];
   x2 = Compile`GetElement[gridxlog, indexlargerx];
   y1 = Compile`GetElement[gridylog, indexlargery - 1];
   y2 = Compile`GetElement[gridylog, indexlargery];
   (*Corresponding values of the function*)
   leny = Length[gridylog];
   z11 = 
     distrlog, (indexlargerx - 2)*leny + indexlargery - 1];
   z21 = 
     distrlog, (indexlargerx - 1)*leny + indexlargery - 1];
   z12 = 
     distrlog, (indexlargerx - 2)*leny + indexlargery];
   z22 = 
     distrlog, (indexlargerx - 1)*leny + indexlargery];
   \[CapitalDelta] = (x2 - x1) (y2 - y1);
    z22,*)((x2 - xlogrand)*(y2 - ylogrand))/\[CapitalDelta]*
      z11 + ((xlogrand - x1)*(y2 - 
         ylogrand))/\[CapitalDelta] z21 + ((x2 - xlogrand) (ylogrand -
      z12 + ((xlogrand - x1) (ylogrand - y1))/((x2 - x1) (y2 - y1))
   ], 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 \
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



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

  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 12, 2023 at 1:42
  • $\begingroup$ @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? $\endgroup$ Commented Dec 12, 2023 at 22:27
  • $\begingroup$ 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. $\endgroup$ Commented Dec 12, 2023 at 22:35


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.