Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
2  
Related question –  Joe Feb 27 '12 at 19:32
    
Also see my question on MathGroup 2 years back: groups.google.com/forum/#!msg/comp.soft-sys.math.mathematica/… –  s0rce Jul 23 '12 at 14:10
add comment

5 Answers

up vote 15 down vote accepted

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 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 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.

share|improve this answer
    
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. –  s0rce Jul 14 '13 at 13:58
    
@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. –  Leonid Shifrin Jul 14 '13 at 14:12
add comment

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]]}]]
share|improve this answer
2  
@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. –  Leonid Shifrin Feb 27 '12 at 10:50
3  
@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. –  kguler Feb 27 '12 at 11:46
    
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 –  Szabolcs Feb 27 '12 at 12:59
add comment

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)]
  ]
share|improve this answer
    
You meant to use this inside Pick, I guess? –  Leonid Shifrin Feb 27 '12 at 14:30
    
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. –  Szabolcs Feb 27 '12 at 14:41
1  
@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. –  Leonid Shifrin Feb 27 '12 at 15:04
1  
@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. –  Szabolcs Feb 27 '12 at 15:10
1  
@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. –  user21 Feb 27 '12 at 15:18
show 6 more comments

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"
]
share|improve this answer
    
This is still bound to be linear time, although perhaps with a small time constant. A good one for unsorted lists though. –  Leonid Shifrin Feb 27 '12 at 10:06
    
@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. –  Szabolcs Feb 27 '12 at 10:12
2  
@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. –  Leonid Shifrin Feb 27 '12 at 10:15
    
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. –  Leonid Shifrin Feb 27 '12 at 11:55
    
@Szabolcs Please see my edit and let me know the timings on your data (in your question) –  rm -rf Feb 27 '12 at 13:22
add comment

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.

share|improve this answer
    
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. –  Szabolcs Feb 27 '12 at 9:50
2  
It is documented here: reference.wolfram.com/mathematica/tutorial/UsingNearest.html –  faleichik Feb 27 '12 at 9:56
1  
I was always curious about how Nearest worked, and I don't know much about data structures, so I finally asked on SO. –  Szabolcs Feb 27 '12 at 10:13
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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