36
$\begingroup$

Suppose we have a sorted list of values. Let's use list = Sort@RandomReal[1, 1000000]; for this example.

I need a fast function window[list, {xmin, xmax}] which will return all list elements $x$ for which $x_\textrm{min} \le x \le x_\textrm{max}$.

How can this be implemented in Mathematica? I am looking both for fast and for elegant solutions. The direct solution is implementing binary search, but there are several ways to do this, and perhaps Mathematica already has something built in that I am not aware of.

Here's the most naïve implementation:

window[list_, {xmin_, xmax_}] := 
 list[[LengthWhile[list, # < xmin &] + 1 ;; LengthWhile[list, # <= xmax &]]]

Summary:

Here are the timings I get for the different solutions for some random data of a million machine reals which also contains duplicates:

  • My original naive solution: 3.85 s

  • Leonid, using binary search: 0.01 s (close to the measurable limit, $\log n$ complexity)

  • R.M., using Clip: 0.59 s (linear time, no sorting required)

  • faleichik, using Nearest: 1.29 s (strangely, this also runs in linear time, by measurement)

  • kguler, using Map (which autocompiles) and Pick: 0.30 s (also linear time, the fastest simple linear time solution so far, does not require sorting either)

For sorted data, the fastest solution is Leonid's, which uses binary search and has logarithmic complexity.

For unsorted data, the fastest (and also one of the simplest) is kguler's. A not so obvious trick was using Boole on the condition to allow it to be automatically compiled.

Improve this question
$\endgroup$
2

6 Answers 6

Reset to default
27
$\begingroup$

As you said, essentially you need binary search, since you have a sorted list and binary search has a logarithmic complexity. However, since

  • the limiting numbers may not be present in the list
  • some numbers may be present more than once

we'd need modified binary search. Here is a possible implementation:

(* maximum number smaller than or equal to the limit *)
bsearchMin[list_List, elem_] :=
  Module[{n0 = 1, n1 = Length[list], m},
    While[n0 <= n1,
     m = Floor[(n0 + n1)/2];
     If[list[[m]] == elem, 
         While[list[[m]] == elem, m++]; 
         Return[m - 1]];
     If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]
    ];
    If[list[[m]] < elem, m, m - 1] 
  ];

and

(* minimum number larger than or equal to the limit *)
bsearchMax[list_List, elem_] :=
  Module[{n0 = 1, n1 = Length[list], m},
    While[n0 <= n1,
      m = Floor[(n0 + n1)/2];
      If[list[[m]] == elem, 
         While[list[[m]] == elem, m--]; 
         Return[m + 1]];
      If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]
    ];
    If[list[[m]] > elem, m, m + 1] 
  ];

With the help of these:

window[list_, {xmin_, xmax_}] :=
  With[{minpos = bsearchMax[list, xmin], maxpos =  bsearchMin[list, xmax]},
    Take[list, {minpos, maxpos}] /; ! MemberQ[{minpos, maxpos}, -1]
  ];
window[__] := {};

For example:

lst = {1, 4, 4, 4, 6, 7, 7, 11, 11, 11, 11, 13, 15, 18, 19, 22, 23, 25, 27, 30}

window[lst, {4, 11}]

(* ==> {4, 4, 4, 6, 7, 7, 11, 11, 11, 11} *)

You can Compile functions bsearchMin and bsearchMax, if you expect lots of duplicate elements (this will speed an inner While loop). Compiling them per se won't improve the speed much (unless you call these very often), since the complexity is logarithmic in any case.

This is certainly generally more efficient than Nearest for sorted lists (perhaps unless you have lots of repeated elements, but then you can compile), because Nearest is a general algorithm which can not take into account the sorted nature of the list. I tried on 10^7 elements list, and it took something 0.0003 seconds for that.

Compiled version

Compiled versions speed up bsearchMin and bsearchMax, but seem not to improve the performance of window[]. See discussion in comments section. 

bsearchMax = Compile[{{list, _Complex, 1}, {elem, _Real}},
  Block[{n0 = 1, n1 = Length[list], m = 0},
    While[n0 <= n1,
      m = Floor[(n0 + n1)/2];
      If[list[[m]] == elem,
        While[m >= n0 && list[[m]] == elem, m--]; Return[m + 1]  ];
      If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]];
    If[list[[m]] > elem, m, m + 1]
  ]
  ,
  RuntimeAttributes -> {Listable},
  CompilationTarget -> "C"
]

bsearchMin = Compile[{{list, _Complex, 1}, {elem, _Real}},
  Block[{n0=1,n1=Length[list],m = 0},
    While[n0<=n1,
      m=Floor[(n0+n1)/2];
      If[list[[m]]==elem,
        While[m<=n1 && list[[m]]==elem,m++]; Return[m-1]  ];
      If[list[[m]]<elem, n0=m+1, n1=m-1]];
    If[list[[m]]<elem,m,m-1]
  ]
  ,
  RuntimeAttributes -> {Listable},
  CompilationTarget -> "C"
]
Improve this answer
$\endgroup$
19
  • $\begingroup$ I wanted to say that I think this should become a built-in part of Mathematica. I use it all the time, definitely one of my favorite answers on the site. $\endgroup$
    – s0rce
    Commented Jul 14, 2013 at 13:58
  • $\begingroup$ @s0rce Then you should know that it likely has a bug. I seem to remember that I had to fix it in this answer, but I did not propagate the changes back. Please see the edit history of that answer to see which changes I made there. $\endgroup$
    – Leonid Shifrin
    Commented Jul 14, 2013 at 14:12
  • $\begingroup$ A minor quibble: I'd use Quotient[n0 + n1, 2] instead of Floor[(n0 + n1)/2]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 15, 2015 at 14:37
  • $\begingroup$ @karolis Your bsearchMin contains mistake, m isn't initialized. And what sample have you used for testing the compiled function? With lst = RandomReal[1, 10^8]; test = window[lst, {0.3, 0.6}]; // AbsoluteTiming the compiled version becomes slower. $\endgroup$
    – xzczd ♦
    Commented Dec 6, 2015 at 5:50
  • 1
    $\begingroup$ @Karolis It was me who fixed the typo. Thanks for providing the compiled versions. $\endgroup$
    – Leonid Shifrin
    Commented Dec 7, 2015 at 8:14
18
$\begingroup$

Using Pick with Boole selector

window[list_, {xmin_, xmax_}] := 
 Pick[list, Boole[xmin <= # <= xmax] & /@ list, 1]

With

 list = Sort@RandomReal[1, 1000000];
 {min, max} = Sort@RandomReal[1, 2];

Timings:

 Table[ClearSystemCache[]; 
 Timing[window[list, {min, max}];], {50}] // Mean
 (* ==> {0.0674, Null} *)

on a laptop with Vista 64bit OS, Intel Core2 Duo T9600 2.80GHz, 8G memory.

UPDATE: Using Pickwith alternative selector arrays:

UnitStep

  windowUnitStep1[list_, {xmin_, xmax_}] := 
  Pick[list, UnitStep[(list-xmin)(xmax-list)], 1]

or

  windowUnitStep2[list_, {xmin_, xmax_}] := 
  Pick[list, UnitStep[list-xmin]UnitStep[xmax-list], 1]

both are twice as fast as Boole.

UnitStep Compiled (Ruebenko's compiled function win)

  windowUnitStep3[list_, {xmin_, xmax_}] := 
  Pick[list, win[list,xmin,xmax], 1]

is twice as fast as uncompiled UnitStep.

Using GatherBy with Boole:

  windowGatherBy[list_, {xmin_, xmax_}] := Last@GatherBy[list, Boole[xmin <= # <= xmax] &]

gives similar timings to window.

Using SparseArray with Part or Take:

The following alternatives attempt to take into account the fact that input data is sorted, thus the first and the last non-zero positions in SparseArray[UnitStep[(list-min)(max-list)]] give the first and the last positions of the portion of input list that satisfy the bounds.

 windowSparseArray1[list_, xmin_, xmax_] := 
 With[{fromTo = SparseArray[UnitStep[(list - xmin) (xmax - list)]][
  "NonzeroPositions"][[{1, -1}]]}, 
  list[[fromTo[[1, 1]] ;; fromTo[[2, 1]]]]]

or

 windowSparseArray2[list_, xmin_, xmax_] := 
 With[{fromTo = SparseArray[UnitStep[(list - xmin) (xmax - list)]][
  "NonzeroPositions"][[{1, -1}]]}, 
  Take[list, {fromTo[[1, 1]], fromTo[[2, 1]]}]]

both give rougly 50 percent speed improvement over window above. Using Ruebenko's compiled UnitStep to construct the array again doubles the speed:

 windowSparseArray3[list_, xmin_, xmax_] := 
 With[{fromTo = SparseArray[win[list,xmin,xmax]][
  "NonzeroPositions"][[{1, -1}]]}, 
  Take[list, {fromTo[[1, 1]], fromTo[[2, 1]]}]]
Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ @Szabolcs By looking at the timings, which are rather impressive here for a linear-time top-level solution, I would guess that Map auto-compiles the test. Coupled with Pick being optimized on packed arrays, this represents a viable linear-time alternative, to my mind. $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 10:50
  • 3
    $\begingroup$ @Szabolcs, and @Leonid, this was intended as a baseline. It is much faster than alternatives using LengthWhile (1.52662), Position (1.4015), Clip(0.1819) or Nearest (0.42962) on the same data set. Of course, a method that explicitly uses binary search is uncomparably better: Leonid's binary search method gives over 40x better results than plain Pick. $\endgroup$
    – kglr
    Commented Feb 27, 2012 at 11:46
  • $\begingroup$ It is really not so obvious why this is so fast. If you regularly use this pattern then could you include some additional explanations? Select is more natural here, but it's also much slower (due to unpacking?) Related: mathematica.stackexchange.com/a/11/12 and mathematica.stackexchange.com/q/1803/12 $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 12:59
17
$\begingroup$

I think Nearest[] is the most effective way. You don't even need to sort your data. 

a = RandomReal[1, 100];
nf = Nearest@a;
xmin = 0.01; xmax = 0.6;
x0 = (xmin + xmax)/2; dx = (xmax - xmin)/2;
nf[x0, {\[Infinity], dx}] // Sort

{0.0117819, 0.013102, 0.0177269, 0.0356801, 0.040019, 0.0504563, \
0.0627056, 0.0749593, 0.0758206, 0.106541, 0.107941, 0.112281, \
0.117172, 0.132445, 0.143151, 0.157252, 0.166585, 0.179652, 0.217876, \
0.241301, 0.242821, 0.254276, 0.258477, 0.267544, 0.268951, 0.280489, \
0.290386, 0.305346, 0.315458, 0.318908, 0.337006, 0.338169, 0.339338, \
0.362153, 0.366946, 0.371712, 0.386563, 0.396061, 0.416329, 0.426874, \
0.430932, 0.439427, 0.460844, 0.473224, 0.475559, 0.476573, 0.479037, \
0.480472, 0.503684, 0.513969, 0.521916, 0.535221, 0.541562, 0.54198, \
0.554534, 0.558954, 0.563491, 0.565873, 0.582683, 0.58919, 0.592807, \
0.593541}

For array of 100 000 numbers it took 0.062 seconds on my machine. For million -- 0.688.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ I did not know about this syntax of a NearestFunction (i.e. an argument of the form {\[Infinity], dx}). It is not documented on the doc page for Nearest. Your answer shows that it was definitely worth asking this question. I do wonder how Nearest works though. It works on two-dimensional data as well, so I doubt it does a plain sort internally. $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 9:50
  • 2
    $\begingroup$ It is documented here: reference.wolfram.com/mathematica/tutorial/UsingNearest.html $\endgroup$
    – faleichik
    Commented Feb 27, 2012 at 9:56
  • 1
    $\begingroup$ I was always curious about how Nearest worked, and I don't know much about data structures, so I finally asked on SO. $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 10:13
14
$\begingroup$

Here are a few approaches:

1: Using Clip

This should be definitely faster than the naïve implementation and is a good un-compiled option for unsorted lists

 window[list_, {xmin_, xmax_}] :=  Clip[list, {xmin, xmax}, {{}, {}}] // Flatten

However, as Leonid notes, this also unpacks the array (causing a slight drop in speed) because the last argument is not numerical, although this can be handled by clipping differently.

2: Using Pick and IntervalMemberQ

This is a straightforward mathematical implementation of the problem, and is again faster than the naïve approach.

window[list_, {xmin_, xmax_}] := Pick[list, 
    IntervalMemberQ[Interval[{xmin, xmax}], list], True]

This will also unpack the array.

3: Forward-backward search (Compiled)

Since you have a sorted list, the following first searches forward till it hits the first element >=xmin and then searches backward till it hits the first element <= xmax and returns everything in between. Compiling to C and parallellizing it makes it very fast (300x faster than naïve, 30x faster than Clip and 170x faster than IntervalMemberQ on my machine).

window = Compile[{{list, _Real, 1}, {xmin, _Real}, {xmax, _Real}},
    Module[{i, j},
        i = 1; While[list[[i]] < xmin, i++];
        j = 1; While[list[[-j]] > xmax, j++];
        list[[i ;; -j]]
    ],
    CompilationTarget -> "C", Parallelization -> True, 
    "RuntimeOptions" -> "Speed"
]
Improve this answer
$\endgroup$
5
  • $\begingroup$ This is still bound to be linear time, although perhaps with a small time constant. A good one for unsorted lists though. $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 10:06
  • $\begingroup$ @Leonid In Mathematica sometimes the fastest solution for input data of practical sizes is not the same as the best complexity solution, so this is pretty useful. In this case I was interested in using the same list over and over, so a pre-processing (like sorting) is affordable. This is why I asked about sorted lists. $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 10:12
  • 3
    $\begingroup$ @Szabolcs And this is why I used binary search in my solution. You can't beat a log with linear for large lists, even in Mathematica where some time constants (e.g. packed array vs. unpacked arrays) are very different and create rather unnatural performance handicaps. So, in this particular case: the complexity is important. $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 10:15
  • $\begingroup$ The less obvious problem here is that Clip unpacks, because the last argument is not numerical. This seriously degrades the performance. This can be fixed by Clip-ping differently and using other methods of eliminating clipped elements, which won't unpack. $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 11:55
  • $\begingroup$ @Szabolcs Please see my edit and let me know the timings on your data (in your question) $\endgroup$
    – rm -rf ♦
    Commented Feb 27, 2012 at 13:22
12
$\begingroup$

Here is my entry. It's O(n), but quite fast, so if you ever have unsorted data, this is a choice:

win = Compile[{{inVec, _Real, 1}, {min, _Real, 0}, {max, _Real, 0}},
  UnitStep[(inVec - min)*(-inVec + max)]
  ]
Improve this answer
$\endgroup$
11
  • $\begingroup$ You meant to use this inside Pick, I guess? $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 14:30
  • $\begingroup$ One interesting thing about kguler's solution compared to this is that it will autos-specialize to the type of input. I only managed to get explicitly compiled solutions (like this one) working if I specified the type of the input (_Real or _Integer). When we let Map auto-compile, it seems to choose the correct one automatically. Of course in this particular situation using a function that works with reals will also work with integers (in practice). But if it is compiled specifically for integers, it will be faster on an integer vector. $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 14:41
  • 1
    $\begingroup$ @Szabolcs I gave an explicit example of how the memoization can be done, in this answer, section "Making JIT-compiled functions...". All that remains is to write a dispatcher which determines the types and forms the right sets of type declarations. This is not difficult at all. $\endgroup$
    – Leonid Shifrin
    Commented Feb 27, 2012 at 15:04
  • 1
    $\begingroup$ @LeonidShifrin We should really start a site blog (please see the chat room). Would you like to volunteer for a couple of posts on these topics? I think it's better suited for a blog than answers. $\endgroup$
    – Szabolcs
    Commented Feb 27, 2012 at 15:10
  • 1
    $\begingroup$ @LeonidShifrin, the default for auto compilation has to be byte code: 1) M- does not ship with a C compiler. 2) compilation to C takes for ever, compared to compiling for byte code. $\endgroup$
    – user21
    Commented Feb 27, 2012 at 15:18
10
$\begingroup$

Assuming that you are interested in multiple different windows for the same list, then the following approach will be much faster than the other answers. Basically, compute a NearestFunction of the data once, and then use that NearestFunction for each window of interest. Here is a function that does this:

WindowFunction[list_] := With[{s = Sort@list},
    WindowFunction[Nearest[s->"Index"], s]
]

WindowFunction[nf_, list_][min_, max_] := Module[{r,s},
    {r, s} = nf[{min, max}][[All,1]];
    Take[list, {r + Boole[list[[r]] < min], s - Boole[list[[s]] > max]}]
]

Here is a comparison with the accepted answer. Sample data;

list = Sort @ RandomReal[1, 10^6];

Compute the WindowFunction (this step is a bit slow, but only needs to be done once):

wf = WindowFunction[list]; //AbsoluteTiming

{0.044266, Null}

Compare:

r1 = wf[.49, .51]; //RepeatedTiming
r2 = window[list, {.49, .51}]; //RepeatedTiming

r1 === r2

{0.000043, Null}

{0.00018, Null}

True

About 4 times faster. One could also add a summary box format for WindowFunction if desired.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Your answers always make me wonder where is the limit to optimizations in Mathematica. I thought I understood its computation model / optimizations fairly well to squeeze most of the speed out, but your answers regularly put me to shame. Hats off, and of course a +1. $\endgroup$
    – Leonid Shifrin
    Commented Jan 7, 2018 at 17:40
  • 1
    $\begingroup$ Also, I will use an opportunity to say this once again, and I am sure that the community as a whole shares that opinion. Your contributions make this already great resource just that much better, and you also set a great example of following up old questions and providing high-quality answers to them - something that is alas not commonly done. It would be great if more of us would follow your path in this regard too. $\endgroup$
    – Leonid Shifrin
    Commented Jan 7, 2018 at 17:45
  • $\begingroup$ I really appreciate your contribution which helped me so much.but there's a fly in the ointment: the method you provided seems only work for the situation that there is no repeating elements in the list. obviously it is easy to fix this. $\endgroup$
    – Ren Witcher
    Commented Mar 24, 2019 at 18:31

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.