Quickly Find the Index in an Array Closest to Some Value

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 ind开发者_Python百科ex (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)

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

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)

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