2
$\begingroup$

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];

Having this tabulated data, I want to make a fast interpolation for the random points belonging to the domain defined by the x,y 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 = 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 = 
    Compile`GetElement[
     distrlog, (indexlargerx - 2)*leny + indexlargery - 1];
   z21 = 
    Compile`GetElement[
     distrlog, (indexlargerx - 1)*leny + indexlargery - 1];
   z12 = 
    Compile`GetElement[
     distrlog, (indexlargerx - 2)*leny + indexlargery];
   z22 = 
    Compile`GetElement[
     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?

Improve this question
$\endgroup$
4
  • 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$
    – Henrik Schumacher
    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$
    – John Taylor
    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$
    – Henrik Schumacher
    Commented Dec 12, 2023 at 22:35
  • $\begingroup$ @HenrikSchumacher : my goal is to quickly sample events with hypothetical new physics particles with some mass mass having the tabulated angle-energy distribution {mass, polar angle, energy, distribution}, that describes the probability of their production at any experimental facility. In the very first step, I generate the angles and energies using the importance sampling method. For this, I first uniformly generate angles and energies within the range interesting to me, and then evaluate the selection weights using the custom interpolation. $\endgroup$
    – John Taylor
    Commented Jul 6 at 22:20

1 Answer 1

Reset to default
1
$\begingroup$

There are two identical (by timings) codes. The first one is

linint = 
 Compile[{{rand\[Theta]log, _Real}, {randElog, _Real}, \
{grid\[Theta]log, _Real, 1}, {gridElog, _Real, 
    1}, {index\[Theta], _Integer}, {indexE, _Integer}, {distrlog, \
_Real, 1}}, 
  Module[{\[Theta]el, Eel, x1, x2, y1, y2, z11, z21, z12, 
    z22, \[CapitalDelta], indexlarger\[Theta], indexlargerE, lenE},
   (*Finding the index of the grids with the elements larged than the \
given random point*)
   \[Theta]el = Compile`GetElement[grid\[Theta]log, index\[Theta]];
   indexlarger\[Theta] = 
    index\[Theta] + UnitStep[rand\[Theta]log - \[Theta]el];
   Eel = Compile`GetElement[gridElog, indexE];
   indexlargerE = indexE + UnitStep[randElog - Eel];
   (*The two points x1, x2 with x1 < x < x2*)
   x1 = Compile`GetElement[grid\[Theta]log, indexlarger\[Theta] - 1];
   x2 = Compile`GetElement[grid\[Theta]log, indexlarger\[Theta]];
   y1 = Compile`GetElement[gridElog, indexlargerE - 1];
   y2 = Compile`GetElement[gridElog, indexlargerE];
   (*Corresponding values of the function*)
   lenE = Length[gridElog];
   z11 = 
    Compile`GetElement[
     distrlog, (indexlarger\[Theta] - 2)*lenE + indexlargerE - 1];
   z21 = 
    Compile`GetElement[
     distrlog, (indexlarger\[Theta] - 1)*lenE + indexlargerE - 1];
   z12 = 
    Compile`GetElement[
     distrlog, (indexlarger\[Theta] - 2)*lenE + indexlargerE];
   z22 = 
    Compile`GetElement[
     distrlog, (indexlarger\[Theta] - 1)*lenE + indexlargerE];
   \[CapitalDelta] = (x2 - x1) (y2 - y1);
   Exp[(*{rand\[Theta]log,x1,x2,randElog,y1,y2,z11,z21,z12,
    z22,*)((x2 - rand\[Theta]log)*(y2 - randElog))/\[CapitalDelta]*
      z11 + ((rand\[Theta]log - x1)*(y2 - 
         randElog))/\[CapitalDelta] z21 + ((x2 - 
         rand\[Theta]log) (randElog - y1))/\[CapitalDelta]*
      z12 + ((rand\[Theta]log - x1) (randElog - y1))/((x2 - x1) (y2 - 
         y1)) z22](*}*)
   ], CompilationTarget -> "C", RuntimeOptions -> "Speed", 
  RuntimeAttributes -> {Listable}, Parallelization -> True]
<< CompiledFunctionTools`
CompilePrint@linint;

which uses the indices of the nearest neighbors of the grids computed using Nearest (which is not formally compilable), index\[Theta],indexE, as the external input:

nearestx = Nearest[gridxlog -> "Index"];
nearesty = Nearest[gridylog -> "Index"];
index\[Theta] = nearestx[xlogrand] // Flatten; // AbsoluteTiming
indexE = nearesty[ylogrand] // Flatten; // AbsoluteTiming

with the definitions gridxlog,gridylog,xlogrand,ylogrand provided in the question.

The second one is

nearestIndexSearch = 
  Compile[{{arr, _Real, 1}, {val, _Real}}, 
   Module[{i = 1}, While[i <= Length[arr] && val >= arr[[i]], i++];
    i - 1], CompilationTarget -> "C", RuntimeOptions -> "Speed"];

(*Compiled Bilinear Interpolation with integrated Nearest \
functionality*)
linintWithNearest = 
  With[{nf = nearestIndexSearch}, 
   Compile[{{xlogrand, _Real}, {ylogrand, _Real}, {gridxlog, _Real, 
      1}, {gridylog, _Real, 1}, {distrlog, _Real, 1}}, 
    Module[{xel, yel, x1, x2, y1, y2, z11, z21, z12, 
      z22, \[CapitalDelta], leny, indexlargerx, indexlargery, 
      nearestindexx, 
      nearestindexy},(*Finding the index of the grids with the \
elements larger than the given random point*)
     nearestindexx = nf[gridxlog, xlogrand];
     nearestindexy = nf[gridylog, ylogrand];
     (*Adjust indices to ensure the correct bounding values*)
     indexlargerx = nearestindexx + 1;
     indexlargery = nearestindexy + 1;
     (*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 = 
      Compile`GetElement[
       distrlog, (indexlargerx - 2)*leny + indexlargery - 1];
     z21 = 
      Compile`GetElement[
       distrlog, (indexlargerx - 1)*leny + indexlargery - 1];
     z12 = 
      Compile`GetElement[
       distrlog, (indexlargerx - 2)*leny + indexlargery];
     z22 = 
      Compile`GetElement[
       distrlog, (indexlargerx - 1)*leny + indexlargery];
     \[CapitalDelta] = (x2 - x1) (y2 - y1);
     Exp[(*{rand\[Theta]log,x1,x2,randElog,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, 
    CompilationOptions -> {"InlineCompiledFunctions" -> True}]];

which replaces Nearest with the simple linear search of the nearest neighbor that may be embedded into the compiled interpolator. For the test example provided in the question, it is as fast as the first one; this is if assuming that the indices nearestx,nearesty will need to be recomputed each time when calling linint (this is my case).

Unfortunately, I did not manage to make a faster interpolator inside Mathematica, which would be desirable.

Regarding external libraries (as it was suggested in the comments to the question), I am afraid that this cannot be done in an OS-independent way, which is crucial for my demands (see an attempt for the Windows case here).

Improve this answer
$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.