5
$\begingroup$

Consider some list

list=RandomReal[{0,1},10^6];

I would like to find the closest element to some value of x that is larger than it, and also extract its position. I may do this straightforward:

xval=RandomReal[{0,1}];
el=Select[list,#>xval&]//Min;
pos=Position[list,el]

but it is slow. Is it possible to do it somehow else? Using, e.g., Nearest (which is fast).

$\endgroup$
1
  • 1
    $\begingroup$ Classical application for a binary search. Try ResourceFunction["BinarySearch"]. $\endgroup$ Dec 10, 2023 at 1:06

3 Answers 3

7
$\begingroup$

Nearest will either give us the value we're looking for or the value to the left in a sorted list.

So we can simply find the nearest element, if it's bigger than the input we're done, otherwise we take the element to the right in the sorted list. One caveat is to handle the case where the input is larger than the largest element of list.

SeedRandom[1];
list = Sort[RandomReal[{0, 1}, 10^6]];
xvals = RandomReal[{0, 1}, 10^6];

nf = Nearest[list -> "Index"];

inds = Flatten[nf[xvals]];

indsup = Clip[inds + UnitStep[xvals - list[[inds]]], {1, Length[list]}];

nearsup = list[[indsup]];
$\endgroup$
5
$\begingroup$

If unsorted, then a straightforward procedural routine might be good:

mingt = Compile[{{a, _Real, 1}, x},
   Module[{y = (1. + $MachineEpsilon) $MaxMachineNumber}, (* Inf` *)
    Do[
     If[x < z < y, y = z],
     {z, a}];
    If[y > $MaxMachineNumber, x, y]
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "CompareWithTolerance" -> False,
   CompilationOptions -> {"InlineExternalDefinitions" -> True}
   ];

SeedRandom[2];
data = RandomReal[1, 10^6];

mingt[data, 0.5] // AbsoluteTiming

(* {0.003916, 0.500002} *)

It's faster than sorting:

Nearest[Sort@data, 0.5] // AbsoluteTiming

(* {0.083824, {0.500002}} *)

It returns x, if there are no greater numbers.

$\endgroup$
2
  • 3
    $\begingroup$ Erm. You compare the time of one lookup of you method with the time to sort and create a search tree with Nearest. The latter is however a one-time cost; just create a NearestFunction nf = Nearest[Sort@data->"Index"]. It will return the position (so that you can potentially increment it). Then calls to nf are very quick, in particular if you query many positions at once. $\endgroup$ Dec 10, 2023 at 0:57
  • $\begingroup$ That's quite right. I answered without thinking that through. Thanks. I suppose if x is fixed and data varies is less likely. Though it could happen, depending on what one is searching for. $\endgroup$
    – Goofy
    Dec 10, 2023 at 1:20
1
$\begingroup$

For the present problem, a compiled version of a binary search should work well. At least it would free you from post processing the output of Nearest, which might be tedious.

cBinarySearchAfter = Compile[{{sortedlist, _Real, 1}, {y, _Real}},
  Block[{a, fa, b, fb, c, fc},
   a = 1;
   b = Length[sortedlist];
   
   fa = Compile`GetElement[sortedlist, a];
   fb = Compile`GetElement[sortedlist, b];
   
   If[y < fa, Return[1];];
   
   If[y >= fb, Return[b + 1];];
   
   (*Always fa<=y<fb;*)
   While[b > a + 1,
    c = a + Quotient[b - a, 2];
    fc = Compile`GetElement[sortedlist, c];
    
    If[fc <= y,
     a = c; fa = fc;
     ,
     b = c; fb = fc;
     ];
    ];
   
   Return[b]
   ],
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  RuntimeOptions -> "Speed"
  ];

Here a usage examle.

n = 1000000;
m = 1000000;
list = RandomReal[{0, 1}, n];
sortedlist = Join[{0.}, Sort[list], {1.}]; // AbsoluteTiming // First

xlist = RandomReal[{0, 1}, m];

pos = cBinarySearchAfter[sortedlist, xlist]; // AbsoluteTiming // First

And @@ (Thread[sortedlist[[pos - 1]] <= xlist])
And @@ (Thread[xlist < sortedlist[[pos]]])

0.076439

0.098603

True

True

It is only half as fast as the Nearest approach without postprocessing:

nf = Nearest[Sort[list] -> "Index"]; // AbsoluteTiming // First
Flatten[nf[xlist]] + 1; // AbsoluteTiming // First

0.101369

0.057467

I think once written in C/C++ and properly parallelized, it would be on par with Nearest or sightly faster.

$\endgroup$
8
  • $\begingroup$ Thanks! How do you think, is it efficient to use the 1D binary search to construct the linear interpolation of some nD data, considering some unstructured grid to which it must be mapped? $\endgroup$ Dec 10, 2023 at 8:43
  • 1
    $\begingroup$ That depends. Is the grid fully unstructured? Like a triangle or tetrahedral mesh? Then you actually want to find the (or a) closest d-dimensional cell that contains the point. If R is you MeshRegion that contains the grid and if xlist is the list of query points, then you get all top-dimensional cells closest to this point by Region`Mesh`MeshNearestCellIndex[R, xlist][[All, 2]]. $\endgroup$ Dec 10, 2023 at 12:41
  • $\begingroup$ The grid is a set of random points within the region covered by the interpolation function. $\endgroup$ Dec 10, 2023 at 12:42
  • 1
    $\begingroup$ I assume you have one function value (scalar or vector) per vertex of the grid to interpolate. Is that right? So if you have a point x, then you first look up in which top-dimensional cell you are. Then you request the coordinates of the vertices of the cell and use them to express x in terms of barycentric coordinates with respect to the vertices. Next to request the function values at the vertices and use the barycentric coordinates to interpolate them linearly. $\endgroup$ Dec 10, 2023 at 14:49
  • 1
    $\begingroup$ Also the barycentric coordinates tell you whether x is really within the cell. (If x lies outside the domain there is still be a closest cell, but one of the barycentric coordinates of x w.r.t. that cell will be negative. $\endgroup$ Dec 10, 2023 at 14:50

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.