Quickly Find the Index in an Array Closest to Some Value

Question

I have an array of values, t, that is always in increasing order (but not always uniformly spaced). I have another single value, x. I need to find the index in t such that t[index] is closest to x. The function must return zero for x < t.min() and the max index (or -1) for x > t.max().

I've written two functions to do this. The first one, f1, is MUCH quicker in this simple timing test. But I like how the second one is just one line. This calculation will be done on a large array, potentially many times per second.

Can anyone come up with some other function with comparable timing to the first but with cleaner looking code? How about something quicker then the first (speed is most important)?

Thanks!

Code:

import numpy as np
import timeit

t = np.arange(10,100000)         # Not always uniform, but in increasing order
x = np.random.uniform(10,100000) # Some value to find within t

def f1(t, x):
   ind = np.searchsorted(t, x)   # Get index to preserve order
   ind = min(len(t)-1, ind)      # In case x > max(t)
   ind = max(1, ind)             # In case x < min(t)
   if x < (t[ind-1] + t[ind]) / 2.0:   # Closer to the smaller number
      ind = ind-1
   return ind

def f2(t, x):
   return np.abs(t-x).argmin()

print t,           '\n', x,           '\n'
print f1(t, x),    '\n', f2(t, x),    '\n'
print t[f1(t, x)], '\n', t[f2(t, x)], '\n'

runs = 1000
time = timeit.Timer('f1(t, x)', 'from __main__ import f1, t, x')
print round(time.timeit(runs), 6)

time = timeit.Timer('f2(t, x)', 'from __main__ import f2, t, x')
print round(time.timeit(runs), 6)

Since your array is sorted, try a binary search. See the answer to this question: http://stackoverflow.com/questions/212358/binary-search-in-python — payne, May 19 '11 at 23:04
I am just leaving work, but want to look at this later. I think you may be able to improve on your first function by short circuiting once you've tested for x < min(t) and x > max(t), but I haven't had a chance to test it yet. — g.d.d.c, May 19 '11 at 23:05

score 7 · Accepted Answer · answered May 19 '11 at 23:57

7

This seems much quicker (for me, Python 3.2-win32, numpy 1.6.0):

from bisect import bisect_left
def f3(t, x):
    i = bisect_left(t, x)
    if t[i] - x > 0.5:
        i-=1
    return i

Output:

[   10    11    12 ..., 99997 99998 99999]
37854.22200356027
37844
37844
37844
37854
37854
37854
f1 0.332725
f2 1.387974
f3 0.085864

answered May 19 '11 at 23:57

Kabie

10,489
1
38
45

Ya, much quicker. However, it does need to be adjusted slightly to account for the special cases where x < t.min and x > t.max – Scott B May 20 '11 at 17:42
It also won't work correctly if the gap is larger: `f3([100, 200], 199)` says `0` although `200` is closer to `199` than `100` is. See [my answer here](https://stackoverflow.com/a/56335408/240443). – Amadan May 28 '19 at 06:24

score 1 · Answer 2 · edited Jul 20 '16 at 09:37

1

Use searchsorted:

t = np.arange(10,100000)         # Not always uniform, but in increasing order
x = np.random.uniform(10,100000)

print t.searchsorted(x)

Edit:

Ah yes, I see that's what you do in f1. Maybe f3 below is easier to read than f1.

def f3(t, x):
    ind = t.searchsorted(x)
    if ind == len(t):
        return ind - 1 # x > max(t)
    elif ind == 0:
        return 0
    before = ind-1
    if x-t[before] < t[ind]-x:
        ind -= 1
    return ind

edited Jul 20 '16 at 09:37

Steve

8,469
1
26
37

answered May 19 '11 at 23:15

Maxime Rischard

11
2

That just gives the index of where the number would be inserted, not necessarily the closest index. That's why his f1 function requires extra steps. – Dan May 19 '11 at 23:18
Ah yes I didn't see the searchsorted in f1. Sorry. – Maxime Rischard May 19 '11 at 23:29

score 1 · Answer 3 · edited May 23 '17 at 12:30

np.searchsorted is binary search (split the array in half each time). So you have to implement it in a way it return the last value smaller than x instead of returning zero.

Look at this algorithm (from here):

def binary_search(a, x):
    lo=0
    hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        midval = a[mid]
        if midval < x:
            lo = mid+1
        elif midval > x: 
            hi = mid
        else:
            return mid
    return lo-1 if lo > 0 else 0

just replaced the last line (was return -1). Also changed the arguments.

As the loops are written in Python, it may be slower than the first one... (Not benchmarked)

Quickly Find the Index in an Array Closest to Some Value

3 Answers3

Linked