Given an ordered set S, find the points in S corresponding to another list

Question

Given a sorted list of numbers $S$, I want to create a function that accepts a list of numbers $L$ and for each number $l \in L$ it returns the index of the largest number $s \in S$ such that $s<l$. For simplicity, assume that all elements of $L$ are greater than the minimum value of $S$. Here is an example of a sorted list $S$:

SeedRandom[13];
S = Sort @ RandomReal[10, 5]

{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}

And, here are a couple examples of the argument to the function:

SeedRandom[10];
list1 = RandomReal[10, 5]
list2 = RandomReal[10, 3]

{6.67917, 8.33874, 4.61316, 4.83263, 9.52033}

{6.0669, 1.22425, 6.13959}

Then, I want to create a function f:

f = findIndices[S];

such that

f[list1]
f[list2]

return:

(*
{2, 4, 2, 2, 5}
{2, 1, 2}
*)

One possibility is to use:

findIndices[s_] := Interpolation[
    Thread[{s, Range@Length@s-1}],
    InterpolationOrder->0,
    "ExtrapolationHandler" -> {Evaluate[Length[s]]&, "WarningMessage" -> False}
]

But, this approach is quite slow when dealing with large arguments:

f = findIndices[S];

tst = RandomReal[{.5, 10}, 10^6];
f[tst]; //AbsoluteTiming

{1.15025, Null}

I'm interested in arguments on the order of 10^6 elements, and ordered sets $S$ on the order 10^4 elements. Is there a faster method?

Answer 1

Small $S$

Simple and relatively fast:

f[list_] := Total[UnitStep[list - #] & /@ S]

is very fast for small $S$:

tst = RandomReal[{.5, 10}, 10^6];
f[tst]; // AbsoluteTiming

{0.0305546, Null}

However, the running time is $O(\operatorname{length}(S) \cdot \operatorname{length}(list))$ - for larger $S$, this will get slow, fast.

---

Large $S$

For large $S$, we have to beat GeometricFunctions'BinarySearch, which is a kernel function. Time to get creative! BinarySearch isn't actually that fast when it's called a lot of times for a short list:

GeometricFunctions`BinarySearch[S, #] & /@ tst; // AbsoluteTiming

{0.615308, Null}

It would be much faster if we could process one long list instead of many short ones. Now, imagine we could somehow use Ordering[Join[tst,S]] to find out which index goes where - that would be about 4 times faster!

Let's look at Ordering[Join[list1, S]]:

{6, 7, 3, 4, 1, 8, 9, 2, 10, 5}

The indices up to Length[list1] point to locations in list1, higher values point to locations in S. So we can go through this list one by one, and count the number of indices into S and assign the current count to the indices in list:

accumulateAndGather = Compile[{{order, _Integer, 1}, {n, _Integer}},
   Module[{result = ConstantArray[0, n], current = 0},
    Do[
     If[o <= n, result[[o]] = current, current++],
     {o, order}];
    result]];

And the lookup would then look like this:

f[list_] := Module[{order},
  order = Ordering[Join[list, S]];
  accumulateAndGather[order, Length[list]]]

This is slightly slower for small S:

tst = RandomReal[{.5, 10}, 10^6];
f[tst]; // AbsoluteTiming

{0.181935, Null}

But much faster for larger S:

S = Sort@RandomReal[10, 10^5];
tst = RandomReal[{.5, 10}, 10^6];
f[tst]; // AbsoluteTiming

{0.188625, Null}

Almost no change! Most of the time is spent in Ordering, and sorting a list with 1,100,000 elements doesn't take much longer than sorting one with 1,000,000 elements.

For comparison, with this S, binary search takes much longer:

GeometricFunctions`BinarySearch[S, #] & /@ tst; // AbsoluteTiming

{64.0738, Null}

The running time for Ordering should be $O(n\log(n))$ where $n:=\operatorname{length}(S) + \operatorname{length}(list)$. Running time for binary search is $O(\log(\operatorname{length}(S)) \operatorname{length}(list))$. If $S$ and $list$ are about the same size, this simplifies to $O(n\log(n))$, too - but with larger constants. Only if $S$ is much bigger than $list$, binary search will be faster.

Answer 2

I wrote down a listable implementation of binary search. For each element $x \in L$, HenrikfindIndices[S,L] returns the indices of the largest $s \in S$ with $s<x$ or $0$ if no such index exists.

HenrikfindIndices = Compile[{{S, _Real, 1}, {x, _Real}},
   Block[{a, b, c, aval, bval, cval},
    a = 1;
    b = Length[S];
    aval = S[[a]];
    bval = S[[b]];
    If[x <= aval,
     0,
     While[b > a + 1, 
      c = a + Quotient[b - a, 2];
      cval = S[[c]];
      If[x < cval, 
       b = c; bval = cval, 
       a = c; aval = cval];
      ];
     If[x <= bval, a, b]
     ]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

For S and tst of same size, it is comparable to nikie's solution on my machine with respect to speed. But for S much larger that tst, this seems to be much faster.

Here is nikie's variant in compiled form.

nikiefindIndices = Compile[{{S, _Real, 1}, {x, _Real, 1}},
   Module[{result, current = 0, n = Length[x]},
    result = Table[0, n];
    Do[
     If[o <= n, 
      result[[o]] = current, 
      current++
     ], 
     {o,Ordering[Join[x, S]]}
    ];
    result
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Test data

S = Sort@RandomReal[{0., 10.}, 10^8];
tst = RandomReal[{-1.0, 11.}, 10^6];
aa = nikiefindIndices[S, tst]; // AbsoluteTiming // First
bb = HenrikfindIndices[S, tst]; // AbsoluteTiming // First

(* 2.44242 *)
(* 0.352638 *)

For the orders of magnitude requested by the OP:

S = Sort@RandomReal[{0., 10.}, 10^4];
tst = RandomReal[{-1.0, 11.}, 10^6];
aa = nikiefindIndices[S, tst]; // AbsoluteTiming // First
bb = HenrikfindIndices[S, tst]; // AbsoluteTiming // First

(* 0.141677 *)
(* 0.090989 *)

Edit: Putting the check for bounds in front saves time. Moreover, the right boundary value is not needed any more. So, we can remove some register juggling. However, this leads only to minor speed-up (about 5-10 %). Here is the reworked code

HenrikfindIndices2 = Compile[{{S, _Real, 1}, {x, _Real}},
  Block[{a, b, c, aval, cval},
   a = 1;
   b = Length[S];
   If[Compile`GetElement[S, b] < x,
    b,
    aval = Compile`GetElement[S, a];
    If[x <= aval,
     0,
     While[b > a + 1,
      c = a + Quotient[b - a, 2];
      cval = Compile`GetElement[S, c];
      If[x <= cval,
       b = c,
       a = c; aval = cval
       ];
      ];
     a
     ]
    ]
   ],
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  RuntimeOptions -> "Speed"
  ]

Answer 3

You could split interval between minimal and maximal element of $S$ to $n$ equal bins, and create list of positions of largest elements of $S$ that are smaller than beginnings of subsequent bins, this operation is $O(n + n_S)$ in time and $O(n)$ in memory and need to be performed only one time for given $S$.

Then for each number $l$ you could test to which bin it falls into, in $O(1)$ time (by shifting and calculating quotient), then you would need to search for appropriate position considering elements of single bin not of whole $S$.

If elements of $S$ are approximately uniformly distributed then taking $n$ equal to number of elements of $S$ will give reasonable distribution, otherwise you may want to sacrifice some memory to get better "bin resolution".

---

Below is a simple implementation that uses linear search starting from position acquired from appropriate bin. More sophisticated version might store also positions of smallest elements of $S$ larger than endings of subsequent bins and, depending on number of elements in bin, could choose between linear and binary search.

Compiled function returning list of positions of largest elements of given ordered list s, that are smaller than beginnings of subsequent bins of given size d:

lessThanBinStartPos = Compile[{{s, _Real, 1}, {d, _Real}},
  Module[{j = 2},
    Table[
      While[Compile`GetElement[s, j] < x, ++j];
      j - 1
      ,
      {x, First@s, Last@s - .5 d, d}
    ]
  ],
  CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

Function returning position of largest element of s that is smaller than given l. It uses list of positions sBinPos that can be calculated using above lessThanBinStartPos.

maxLessPosInternal = Compile[{{s, _Real, 1}, {sBinPos, _Integer, 1}, {d, _Real}, {l, _Real}},
  Module[{min, max, result},
    min = Compile`GetElement[s, 1];
    If[l <= min, Return@0];
    max = Compile`GetElement[s, Length@s];
    If[max < l, Return@Length@s];
    If[max == l, Return[Length@s - 1]];

    result = Length@s;
    Do[
      If[Compile`GetElement[s, i] >= l,
        result = i - 1;
        Break[]
      ],
      {i, Compile`GetElement[sBinPos, Quotient[l - min, d] + 1], Length@s}
    ];
    result
  ],
  CompilationTarget -> "C", RuntimeOptions -> "Speed",
  RuntimeAttributes -> {Listable}, Parallelization -> True
];

Helper function for special case of one-element s that have same "listability behavior" as maxLessPosInternal.

maxLessPosInternalSingle = Compile[{{x, _Real}, {l, _Real}},
  If[x < l, 1, 0],
  CompilationTarget -> "C", RuntimeOptions -> "Speed",
  RuntimeAttributes -> {Listable}, Parallelization -> True
];

Final function that takes ordered list s and, optionally, number of bins, and returns function calculating position of largest element of s that is smaller than given argument.

maxLessPos // ClearAll
maxLessPos@s_ := maxLessPos[s, Length@s]
maxLessPos[{}, _] = Compile[{{l, _Real}}, 0,
  CompilationTarget -> "C", RuntimeOptions -> "Speed",
  RuntimeAttributes -> {Listable}, Parallelization -> True
];
maxLessPos[{x_}, _] := maxLessPosInternalSingle[x, #]&;
maxLessPos[s_List, n_Integer?Positive] :=
  With[{d = (Last@s - First@s) / n},
    With[{sBinPos = lessThanBinStartPos[s, d]},
      maxLessPosInternal[s, sBinPos, d, #]&
    ]
  ]

Basic tests:

SeedRandom@13;
S = Sort@RandomReal[10, 5];
SeedRandom@10;
list1 = RandomReal[10, 5];
list2 = RandomReal[10, 3];

f = maxLessPos@S
f@list1
f@list2
(* maxLessPosInternal[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {1, 1, 1, 2, 2}, 1.65542, #1] & *)
(* {2, 4, 2, 2, 5} *)
(* {2, 1, 2} *)

Basic benchmarks:

SeedRandom@0
S = Sort@RandomReal[{0., 10.}, 10^4];
tst = RandomReal[{-1.0, 11.}, 10^6];

(resNikie = nikiefindIndices[S, tst]) // MaxMemoryUsed // RepeatedTiming
(resHenrik = HenrikfindIndices[S, tst]) // MaxMemoryUsed // RepeatedTiming
(resLS = LSIndices[S, tst]) // MaxMemoryUsed // RepeatedTiming
(resWiz = fn2[S, tst]) // MaxMemoryUsed // RepeatedTiming
(resJkuczm = maxLessPos[S]@tst) // MaxMemoryUsed // RepeatedTiming
resNikie === resHenrik === resLS === resWiz === resJkuczm
(* {0.19, 40160656} *)
(* {0.21,  8001568} *)
(* {0.35, 32278888} *)
(* {0.21, 40402880} *)
(* {0.060, 8081904} *)
(* True *)

On my computer maxLessPos was fastest, nikiefindIndices was next - over three times slower.

Answer 4

I'll have to spend time digesting the other answers here. In the mean time, here is an approach using Nearest, which is only competitive with the other answers in M11.1+.

The basic idea is that Nearest[S -> "Index"] will produce a NearestFunction object. When this NearestFunction object is applied to a number, it will return the index of the nearest point in $S$, and this nearest point will be either the nearest point less than the input, or the nearest point greater than the input. If it is the nearest point less than the input, than we are done. If it is the nearest point greater than the input, than we just subtract 1 from the index. The only possible issue is when the NearestFunction returns more than one point when multiple points are equally near, but in this case we can just use the first point given. If $S$ has no duplicates, than one point will be to the left, and one point will be to the right.

Now, what makes this approach really fast is that NearestFunction objects can accept lists, and is designed to work very quickly for lists (it is designed to be vectorized). Determining whether the number is to the left or right and subtracting 1 if it is to the right can also be vectorized. Finally, the process of creating the NearestFunction need only be done once for any particular ordered set $S$.

LeftNeighbor[s_] := LeftNeighborFunction[s, Nearest[s->"Index"]]
LeftNeighbor[s_, list_] := LeftNeighbor[s][list]

LeftNeighborFunction[s_, nf_][list_] := With[{n = nf[list][[All,1]]},
    n - UnitStep[s[[n]] - list]
]

MakeBoxes[i : LeftNeighborFunction[s_, nf_], StandardForm] ^:= Module[
    {
    len = Length[s],
    g = FirstCase[ToBoxes[nf], _GraphicsBox, GraphicsBox[Point[{0,0}]], Infinity]
    },

    BoxForm`ArrangeSummaryBox[
        LeftNeighborFunction,
        i,
        RawBoxes@g,
        {
        BoxForm`MakeSummaryItem[{"Data points: ", Length[s]}, StandardForm],
        BoxForm`MakeSummaryItem[{"Range: ", MinMax[s]}, StandardForm]
        },
        {},
        StandardForm,
        "Interpretable"->True
    ]
]

Let's see LeftNeighbor in action. Here is some data:

SeedRandom[1]
S = Sort @ RandomReal[10, 10^5];
tst = RandomReal[{-1, 11}, 10^6];

Here is how to use LeftNeighbor:

f = LeftNeighbor[S]; //AbsoluteTiming
r1 = f[tst]; //AbsoluteTiming

{0.003054, Null}

{0.113312, Null}

(alternatively, one could use LeftNeighbor[S, tst], but in my use case, $S$ doesn't change, so computing the LeftNeighborFunction only once is better). Finally, let's compare LeftNeighbor to some of the other answers:

r2 = HenrikfindIndices[S, tst]; //AbsoluteTiming
r3 = LSIndices[S, tst]; //AbsoluteTiming
jfunc = maxLessPos[S]; //AbsoluteTiming
r4 = jfunc @ tst; //AbsoluteTiming

r1 === r2 === r3 === r4

{0.139673, Null}

{0.24236, Null}

{0.001829, Null}

{0.043159, Null}

True

The answer by @jkuczm is the clear winner, with the Nearest approach second. I very much like the approach of preprocessing the ordered set $S$ to create a function, and then using the function on test sets.

Answer 5

I would like to contribute another version of vectorized binary search, which is an adaptation of my answer from this old thread. I post this not so much due to a speed advantage (which isn't that significant, and isn't even always there), but to show how one can implement vectorized control flow inside Compile - which in my opinion is an interesting technique.

Code

The main function looks as follows:

ClearAll[vectorizedBinarySearch];
vectorizedBinarySearch = 
  Compile[{{list, _Real, 1}, {elems, _Real, 1}}, 
    Module[{n1 = {1}, n0 = {1}, m = {1}, un1 = {0}, prevm = {0}, 
      len = Length[list], otherlen = Length @ elems},
      n0 = m = Table[1, {otherlen}];
      n1 = Table[len, {otherlen}];
      While[ True,
        prevm = m;
        m = Floor[(n0 + n1)/2] ;
        m = m + BitXor[Unitize[m], 1]; 
        un1 = UnitStep[list[[m]] - elems];
        If[m == prevm, Break[]];
        n1 = n1  + (m - n1 - 1) * un1;
        n0 = (n0 - m - 1) * un1 + m + 1;
      ];
      m * (1 - un1)
    ], 
    CompilationTarget -> "C", RuntimeOptions -> "Speed" , 
    RuntimeAttributes -> {Listable}, Parallelization -> True
  ];

and then the function that computes the indices looks like

ClearAll[LSIndices];
LSIndices[s_List, tst_List, shunkSize_: 5000] :=
  Join @@ vectorizedBinarySearch[
    s, 
    Partition[tst, shunkSize, shunkSize, {1, 1}, {}]
  ];

Explanation

Let us ignore for a moment the Listable attribute of vectorizedBinarySearch (which is not essential for its operation, but simply allows one to use several kernels), and focus on the body of the function (While loop). At any given loop iteration, variables n0, n1 and m store low and high binary search thresholds, and candidate element positions, for all elements elems from the second list at once.

What happens then is that every iteration of the vectorized While loop performs an iteration of a binary search for all elements elems at once. The If statement in a classical binary search also has been vectorized, using mask

un1 = UnitStep[list[[m]] - elems]

The line

m = m + BitXor[Unitize[m], 1]

is needed since in this approach, we want to avoid zero positions in intermediate iterations. We recover the zero positions in the last line

m * (1 - un1)

Also in this approach, I had to maintain the value of m from the previous iteration of the loop, to know when to break out of it.

The top-level function LSIndices simply splits the second list in sublists of size chunkSize (which is an adjustable parameter), and then joins the result back to a single list. This is done to make use of the Listable attribute of vectorizedBinarySearch and at the same time avoid massive allocations / deallocations on the heap, which would happen if we feed the entire large second list of elements to vectorizedBinarySearch at once. Both of these measures improve the performance.

Benchmarks

The benchmarks I performed based on the code / data from the answer of Henrik Schumacher, show that on my system (Mac OS X 10.10.5, Intel Core i7), for smaller data, LSIndices performs roughly twice slower than the other two, while for larger data, it has a slight edge (5 - 10 %) over HenrikfindIndices.

Notes

As I said at the top of the post, the main reason I posted this is to show-case vectorized control flow inside Compile, as an interesting (IMO) alternative to more traditional techniques involving Compile.

Answer 6

For your consideration, based on (75582) and (30648).

On your specific tst I find this seven times faster than findIndices[S].

fn2[S_, L_] :=
  Module[{o, u, a},
    o = Ordering @ Join[L, S];
    u = UnitStep[o - Length@L - 1];
    a = Accumulate[u];
    a[[o]] = a;
    Take[a, Length@L]
  ]

fn2[S, list1]
fn2[S, list2]

{2, 4, 2, 2, 5}

{2, 1, 2}

Variation

A variation that pre-sorts $L$ is useful for short $L$ but comes at a slight cost for long $L$.

fn3[S_, L_] :=
  Module[{oL, o, u},
    oL = Ordering[L];
    o = Ordering @ Join[S, L[[oL]] ];
    u = BitXor[1, UnitStep[o - Length@S - 1]];
    Pick[Accumulate @ u, u, 0][[oL]]
  ]

Performance

nikie already posted an Ordering method; here is a simplistic benchmark including it.

nikie[S_, list_] := Module[{order}, order = Ordering[Join[list, S]];
  accumulateAndGather[order, Length[list]]]

big = Sort @ RandomReal[1, 1*^5];

Needs["GeneralUtilities`"]

BenchmarkPlot[
  {fn2[big, #] &, fn3[big, #] &, nikie[big, #] &},
  RandomReal[1, #] &,
  10^Range[7],
  Joined -> True
]

Answer 7

Eric Towers was having trouble setting up a recursion in the Mathematica language. Here is my take on what I believe he wanted to do, with a little help from Henrik.

find[{}] := 0;

find[{s_}] := Boole[# >= s];

find[s_List] := 
  2^⌊Log2 @ N @ Length @ s⌋ /. p_ :>
    If @@ {# >= s[[p]], p + find[ s[[p + 1 ;;]] ], find[ s[[;; p - 1]] ]}

f = Evaluate[find[S]] &

If[#1 >= 7.95001, 4 + Boole[#1 >= 8.6823], 
  If[#1 >= 4.56535, 2 + Boole[#1 >= 7.04274], Boole[#1 >= 0.405196]]] &

f /@ list1

{2, 4, 2, 2, 5}

Answer 8

---

Added in Edit: Let's decide to do exactly what the OP says: produce a functor that takes a list, processes it, and returns a Listable function taking floating point numbers to indices. What we really want to do is construct and return something like

If[# >= 7.95001, 4 + If[# >= 8.6823,1,0], 
  If[# >= 4.56535, 2 + If[# >= 7.04274,1,0], If[# >= 0.405196,1,0]]] &

when handed the example list of the Question. This should be easy enough by composing partially specialized If[-,-,-] functions recursively down to the leaves of the binary search expression tree.

---

You might think something like this would work

Clear[findIndices];
findIndices[S_] := Module[{
   binSearch, f, x},

  binSearch[{}] :> 0;
  binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0];
  binSearch[s_List /; Length[s] > 1] :> Module[{
     po2, lower, upper},
    po2 = 2^Floor[Log[2, Length[s]]];
    lower = Take[s, po2 - 1];
    upper = Take[s, {po2 + 1, -1}];
    If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
     binSearch[Evaluate[lower]]]
   ];

  f = Function[x,
    Evaluate[
     FixedPoint[
      Evaluate,
      binSearch[S]
     ]
    ]
   ];
  SetAttributes[f, Listable];
  f
 ]

But the utterly inscrutable evaluation rules ensure that you can never actually evaluate a recursive function except when doing so is a mistake. (I would like to say this is humour, but it is not.)

---

Added in Edit: (We do eventually overcome Mathematica's unintelligible evaluation by circumventing it entirely. This was present in the original answer and appears below this long trek through M'ma's failure to be a functional programming tool.) What do I mean by "except when doing so is a mistake"? Forget to put in your base case for a recursion:

Clear[f]
f[n_] = f[n - 1]
  (*  $IterationLimit::itlim: Iteration limit of 4096 exceeded. >>  *)
  (*  Hold[f[-1 + (-4095 + n)]]  *)

Mathematica is quite ready to evaluate that recursive expansion. Now let's compare that behaviour with the semantics of the above code. We have a bunch of delayed rules for binSearch expressions (and the relevant base cases are present). We repeatedly Evaluate binSearch[S] and the resulting expressions in the vain hope that this will result in recursive expansion of binSearch. Let's replace f = ... with

Print[binSearch[S]];
Print[Evaluate[binSearch[S]]];

Then follow the modified findIndices with

findIndices[S]
  (*  binSearch$59538[{0.405196,4.56535,7.04274,7.95001,8.6823}]  *)
  (*  binSearch$59538[{0.405196,4.56535,7.04274,7.95001,8.6823}]  *)

so those delayed rules don't do anything and Evaluate doesn't seem to evaluate anything. Let's try assignments (Some of the following variants appear, possibly spliced with other variants, in the first edit to this post.)

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x}, binSearch[{}] = 0;
  binSearch[s_List /; Length[s] == 1] = If[x >= s[[1]], 1, 0];
  binSearch[s_List /; Length[s] > 1] = 
   Module[{po2, lower, upper}, po2 = 2^Floor[Log[2, Length[s]]];
    lower = Take[s, po2 - 1];
    upper = Take[s, {po2 + 1, -1}];
    If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
     binSearch[Evaluate[lower]]]];

  Print[binSearch[S]];
  Print[Evaluate[binSearch[S]]];
 ]
findIndices[S]
  (*  Part::partd: Part specification s$[[1]] is longer than depth of object. >>  *)
  (*  ...  *)

Unsurprisingly, the immediate evaluation of the RHSs of the Sets fails. SetDelayed?

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x}, binSearch[{}] := 0;
  binSearch[s_List /; Length[s] == 1] := If[x >= s[[1]], 1, 0];
  binSearch[s_List /; Length[s] > 1] := 
   Module[{po2, lower, upper}, po2 = 2^Floor[Log[2, Length[s]]];
    lower = Take[s, po2 - 1];
    upper = Take[s, {po2 + 1, -1}];
    If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
     binSearch[Evaluate[lower]]]];

  Print[binSearch[S]];
  Print[Evaluate[binSearch[S]]];
 ]
findIndices[S]
  (*  If[x$61633 >= 7.95001, po2$61634 + binSearch$61633[Evaluate[upper$61634]], binSearch$61633[Evaluate[lower$61634]]]  *)
  (*  If[x$61633 >= 7.95001, po2$61635 + binSearch$61633[Evaluate[upper$61635]], binSearch$61633[Evaluate[lower$61635]]]  *)

Well, that's somewhat better. binSearch is expanded exactly once, but the variable binding in the inner Module is ignored so that the result references the Globally unresolvable names po2, binSearch, upper, and lower.

Well hmm... Relying on function expansion is a non-starter. Maybe we can get this by rewriting the root expression binSearch[S], ReplaceRepeateding until the result stops changing...

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :> 
      Module[{po2, lower, upper}, po2 = 2^Floor[Log[2, Length[s]]];
       lower = Take[s, po2 - 1];
       upper = Take[s, {po2 + 1, -1}];
       If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
        binSearch[Evaluate[lower]]]]
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$61671 >= 7.95001, po2$61672 + binSearch$61671[Evaluate[upper$61672]], binSearch$61671[Evaluate[lower$61672]]]  *)
  (*  If[x$61671 >= 7.95001, po2$61672 + binSearch$61671[Evaluate[upper$61672]], binSearch$61671[Evaluate[lower$61672]]]  *)

... so apparently RuleDelayed and Module have the option to not bother doing any of that variable binding they're documented to do. Let's force that binding to be a little more prompt with With ... and we have to nest Withs since the second and third local depend on the value of the first local ...

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :>
      With[{po2 = 2^Floor[Log[2, Length[s]]]},
       With[{lower = Take[s, po2 - 1],
         upper = Take[s, {po2 + 1, -1}]},
        If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
         binSearch[Evaluate[lower]]]
        ]]
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$62332 >= 7.95001, 4 + binSearch$62332[Evaluate[{8.6823}]], binSearch$62332[Evaluate[{0.405196, 4.56535, 7.04274}]]]  *)
  (*  If[x$62332 >= 7.95001, 4 + binSearch$62332[Evaluate[{8.6823}]], binSearch$62332[Evaluate[{0.405196, 4.56535, 7.04274}]]]  *)

... slightly better. At least With can do Module's job for it and we get evaluated copies of upper and lower, but not evaluated enough. Perhaps strip the Evaluates?

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :>
      With[{po2 = 2^Floor[Log[2, Length[s]]]},
       With[{lower = Take[s, po2 - 1],
         upper = Take[s, {po2 + 1, -1}]},
        If[x >= s[[po2]], po2 + binSearch[upper], binSearch[lower]]
        ]]
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$62345 >= 7.95001, 4 + If[x$62345 >= {8.6823}[[1]], 1, 0],
       With[{po2$ = 2^Floor[Log[2, Length[{0.405196, 4.56535, 7.04274}]]]},
        With[{lower$ = Take[{0.405196, 4.56535, 7.04274}, po2$ - 1],
          upper$ = Take[{0.405196, 4.56535, 7.04274}, {po2$ + 1, -1}]},
         If[x$62345 >= {0.405196, 4.56535, 7.04274}[[po2$]],
          po2$ + binSearch$62345[upper$], binSearch$62345[lower$]]]]]  *)      (*  If[x$62345 >= 7.95001, 4 + If[x$62345 >= {8.6823}[[1]], 1, 0],
       With[{po2$ = 2^Floor[Log[2, Length[{0.405196, 4.56535, 7.04274}]]]},
        With[{lower$ = Take[{0.405196, 4.56535, 7.04274}, po2$ - 1],
          upper$ = Take[{0.405196, 4.56535, 7.04274}, {po2$ + 1, -1}]},
         If[x$62345 >= {0.405196, 4.56535, 7.04274}[[po2$]],
          po2$ + binSearch$62345[upper$], binSearch$62345[lower$]]]]]  *)

so With isn't actually binding any values. Maybe we can help it out by not using locals...

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :>
      With[{po2 = 2^Floor[Log[2, Length[s]]]},
       If[x >= s[[po2]], po2 + binSearch[Take[s, {po2 + 1, -1}]], 
        binSearch[Take[s, po2 - 1]]]
      ]
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$62362 >= 7.95001, 4 + binSearch$62362[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {4 + 1, -1}]], binSearch$62362[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]  *)
  (*  If[x$62362 >= 7.95001, 4 + binSearch$62362[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {4 + 1, -1}]], binSearch$62362[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]  *)

... so now With has changed its mind and actually binds po2 to a value. Mysterious. However, there's no change in expanding the once nested calls to binSearch; they're still unexpanded. Why? Because the head of Take[...] is not List. So why has po2 been evaluated in place, but 4-1, Take[...], et al, not? Beats me. Let's try telling them to be evaluated.

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :>
      With[{po2 = 2^Floor[Log[2, Length[s]]]},
       If[x >= s[[po2]], 
        po2 + binSearch[Evaluate[Take[s, {po2 + 1, -1}]]], 
        binSearch[Evaluate[Take[s, po2 - 1]]]]
       ]
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$62371 >= 7.95001, 4 + binSearch$62371[Evaluate[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {4+1, -1}]]], binSearch$62371[Evaluate[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]]  *)
  (*  If[x$62371 >= 7.95001, 4 + binSearch$62371[Evaluate[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {4+1, -1}]]], binSearch$62371[Evaluate[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]]  *)

That certainly made the result longer. Didn't evaluate anything, though. Maybe we should catch that non-List head and force it to evaluate.

Clear[findIndices];
findIndices[S_] := Module[
  {binSearch, f, x},
  f = binSearch[S] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :>
      With[{po2 = 2^Floor[Log[2, Length[s]]]},
       If[x >= s[[po2]], po2 + binSearch[Take[s, {po2 + 1, -1}]], 
        binSearch[Take[s, po2 - 1]]]
       ],
     binSearch[other_ /; Head[other] =!= List] :> 
      (Print["Hi."]; binSearch[Evaluate[other]])
     };

  Print[f];
  Print[Evaluate[f]];
  ]
findIndices[S]
  (*  If[x$543 >= 7.95001, 4 + binSearch$543[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823},{4+1, -1}]], binSearch$543[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]  *)
  (*  If[x$543 >= 7.95001, 4 + binSearch$543[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823},{4+1, -1}]], binSearch$543[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, 4-1]]]  *)

Our rule for non-List heads was never called (because it would have printed "Hi." if it had). So the problem is not that we have a non-List head. Although we have a non-List head. I wonder which of the implicit and poorly documented evaluation rules is causing the Takes to not be Takes except when they are.

Notice: At no time did our semantics change: Recursively evaluate to the leaves. We're told Mathematica is talented at repeated string rewriting and recursive function expansion. But apparently not.

Maybe if we could have set binSearch to have attribute EvaluateAll (if that were a thing) the Kernel's randomly applied evaluation and variable binding would have actually done what we said we wanted.

And here's where I quit pretending the language is capable of functional programming. The functional programming fails because we have no control over when those Takes evaluate. The only way to proceed is to be precisely non-functional: hide the functions wrapping the Takes from Mathematica in strings, then re-interpret the strings as expressions after the recursion completes. This is the stupidest possible thing to have to do. The cognitive load in guessing which apply of several nested "this evaluates except ..."s is unworkable.

Even if someone manages to make some variant of this work, I'm not interested. The semantics were straightforward. The evaluation was borked by the Kernel at every step.

---

It would seem the only way to overcome Mathematica's psychic ability to engage in recursion only when you don't want it is to circumvent the evaluation rules by constructing string expressions, so you can actually control evaluation.

Clear[binSearch];
binSearch["{}"] = "0";
binSearch[s_String /; Length[ToExpression[s]] == 1] := 
  "If[#>=" <> ToString[ToExpression[s][[1]]] <> ",1,0]";
binSearch[s_String] := Module[{
   sexpr, po2, lower, upper},
  sexpr = ToExpression[s];
  po2 = 2^Floor[Log[2, Length[sexpr]]];
  lower = Take[sexpr, po2 - 1];
  upper = Take[sexpr, {po2 + 1, -1}];
  "If[#>=" <> ToString[sexpr[[po2]]] <> "," <> ToString[po2] <> 
   " + " <> binSearch[ToString[upper]] <> "," <> 
   binSearch[ToString[lower]] <> "]"
 ]

Clear[findIndices];
findIndices[S_] := Module[{
   f, x},
  f = Function[Evaluate[ToExpression[binSearch[ToString[S]]]]];
  (* SetAttributes[f, Listable];  *) (*  Don't bother.  Does nothing.  *)
  f
 ]

Note the comment about SetAttributes. Nothing you can set in here will cause the f in the Global namespace to be Listable, not even explicitly making it so in the Global namespace.

Clear[f]
f = findIndices[S];
SetAttributes[f, Listable]
f
 (*  If[#1 >= 7.95001, 4 + If[#1 >= 8.6823, 1, 0], 
       If[#1 >= 4.56535, 2 + If[#1 >= 7.04274, 1, 0], 
         If[#1 >= 0.405196, 1, 0]]] &  *)

f /@ {1.1, 5.1, 9.1}
 (*  {1, 2, 5)  *)

What do I mean setting Listable does nothing?

f[{1.1, 5.1, 9.1}]
 (*  If[{1.1, 5.1, 9.1} >= 7.95001, 
      4 + If[{1.1, 5.1, 9.1} >= 8.6823, 1, 0], If[{1.1, 5.1, 9.1} >= 4.56535, 
        2 + If[{1.1, 5.1, 9.1} >= 7.04274, 1, 0], If[{1.1, 5.1, 9.1} >= 0.405196, 1, 0]]]  *)

(sigh) What I wouldn't give for a Mathematica that wasn't an impediment to functional programming...

Oh, right. Timing.

tst = RandomReal[{0.5, 10}, 10^6];
f /@ tst; // AbsoluteTiming
 (*  {0.172239, Null}  *)

I imagine it would be a little faster Compiled, but I'm done fighting with this language for a few days.

---

Edit: 20170910T0429Z

Compilation makes it about 25-times slower.

Clear[cf];
cf = Compile[{{x, _Real}}, f[x]]

cf /@ tst; // AbsoluteTiming
 (*  {4.37982, Null}  *)

As for f, setting cf listable does nothing, so I don't waste space on it.

It's worth pointing out that all of the Compiled versions here risk erroneous output due to precision loss. Test values exceedingly close to but greater than a separating value can compare equal to the separating value when coerced to _Real, so will be reported in the bin one less than their actual bin.

---

---

Since I get the impression Mr. Wizard is unable to generate his own demonastrations of inscrutable evaluation fails...

Clear[findIndices];
findIndices[S_] := Module[{binSearch, f, x},
  binSearch[{}] := 0;
  binSearch[s_List /; Length[s] == 1] := If[x >= s[[1]], 1, 0];
  binSearch[s_List /; Length[s] > 1] := 
   Module[{po2, lower, upper}, po2 = 2^Floor[Log[2, Length[s]]];
    lower = Take[s, po2 - 1];
    upper = Take[s, {po2 + 1, -1}];
    If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
     binSearch[Evaluate[lower]]]];
  Print[binSearch[S]];
  f = Function[x, Evaluate[FixedPoint[Evaluate, binSearch[S]]]];
  SetAttributes[f, Listable];
  f
  ]
f = findIndices[S]
(*  If[x$1452 >= 7.95001, po2$1453 + binSearch$1452[Evaluate[upper$1453]],
     binSearch$1452[Evaluate[lower$1453]]]  *)
(*  Function[x$, 
     If[x$1452 >= 7.95001, po2$1454 + binSearch$1452[Evaluate[upper$1454]], 
      binSearch$1452[Evaluate[lower$1454]]]]  *)


Clear[findIndices];
findIndices[S_] := Module[{binSearch, f, x},
  f = Function[x, binSearch[S]] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :> 
      Module[{po2, lower, upper}, po2 = 2^Floor[Log[2, Length[s]]];
       lower = Take[s, po2 - 1];
       upper = Take[s, {po2 + 1, -1}];
       If[x >= s[[po2]], po2 + binSearch[Evaluate[upper]], 
        binSearch[Evaluate[lower]]]]
     };
  SetAttributes[f, Listable];
  f
  ]
f = findIndices[S]
(*  Function[x$, Module[{po2$, lower$, upper$}, 
     po2$ = 2^Floor[Log[2, Length[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}]]]; 
     lower$ = Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, po2$ - 1]; 
     upper$ = Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, {po2$ + 1, -1}]; 
     If[x$1792 >= {0.405196, 4.56535, 7.04274, 7.95001, 8.6823}[[po2$]], 
      po2$ + binSearch$1792[Evaluate[upper$]], binSearch$1792[Evaluate[lower$]]]]]  *)


Clear[findIndices];
findIndices[S_] := Module[{binSearch, f, x},
  f = Function[x, binSearch[S]] //. {
     binSearch[{}] :> 0,
     binSearch[s_List /; Length[s] == 1] :> If[x >= s[[1]], 1, 0],
     binSearch[s_List /; Length[s] > 1] :> Module[{po2},
       po2 = 2^Floor[Log[2, Length[s]]];
       If[x >= s[[po2]], 
        po2 + binSearch[Evaluate[Take[s, {po2 + 1, -1}]]], 
        binSearch[Evaluate[Take[s, po2 - 1]]]]]
     };
  SetAttributes[f, Listable];
  f
  ]
f = findIndices[S]
(*  Function[x$, 
     Module[{po2$}, po2$ = 2^Floor[
Log[2, Length[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}]]]; 
      If[x$2307 >= {0.405196, 4.56535, 7.04274, 7.95001, 8.6823}[[po2$]], 
       po2$ + binSearch$2307[Evaluate[Take[{0.405196, 4.56535, 
             7.04274, 7.95001, 8.6823}, {po2$ + 1, -1}]]], binSearch$2307[
        Evaluate[Take[{0.405196, 4.56535, 7.04274, 7.95001, 8.6823}, po2$ - 1]]]]]]  *)


Clear[findIndices];
findIndices[S_] := Module[{binSearch, f, x}, binSearch[{}] := 0;
  binSearch[s_List /; Length[s] == 1] := If[# >= s[[1]], 1, 0];
  binSearch[s_List /; Length[s] > 1] :=
   With[{po2 = 2^Floor[Log[2, Length[s]]]},
    With[{lower = Take[s, po2 - 1],
      upper = Take[s, {po2 + 1, -1}]},
     If[# >= s[[po2]], po2 + binSearch[upper], binSearch[lower]]]];
  binSearch[other_] := binSearch[Evaluate[other]];
  Print[binSearch[S]];
  f = Function[Evaluate[FixedPoint[Evaluate, binSearch[S]]]];
  SetAttributes[f, Listable];
  f]
f = findIndices[S]
(*  If[#1 >= 7.95001, 4 + binSearch$2384[{8.6823}], 
         binSearch$2384[{0.405196,4.56535, 7.04274}]]  *)
(*  If[#1 >= 7.95001, 4 + binSearch$2384[{8.6823}], 
         binSearch$2384[{0.405196, 4.56535, 7.04274}]]] &  *)