Wydajność numpy.searchsorted jest słaba w przypadku strukturowanych tablic

Question

_{Przepraszam z góry, jeśli używam niewłaściwych warunków, nie krępuj się.}

Mam posortowaną tablicę z dtype'<f16, |S30'. Kiedy używam searchsorted na pierwszym polu, działa bardzo wolno (około 0,4 sekundy dla 3 milionów elementów). To znacznie dłużej niż bisect, aby zrobić to samo na zwykłej liście krotek w Pythonie.

%timeit a['f0'].searchsorted(400.) 
1 loops, best of 3: 398 ms per loop 

Jednakże, jeśli mogę skopiować część pływaka na inny, odrębny tablicy, wyszukiwanie jest szybsze niż bisect:

b = a['f0'].copy() 

%timeit b.searchsorted(400.) 
1000000 loops, best of 3: 945 ns per loop 

moje pytania są następujące:

1. Czy robię coś źle czy jest to regresja w NumPy?
2. Czy istnieje sposób obejścia tego bez powielania danych?

Answer 1

Pamiętam, że widziałem to jakiś czas temu. Jeśli dobrze pamiętam, myślę, że searchsorted tworzy tymczasową kopię danych, gdy dane nie są ciągłe. Jeśli będę miał czas później, zajrzę do kodu, aby potwierdzić, że to się dzieje (a może ktoś, kto zna ten kod, może to potwierdzić).

W międzyczasie, jeśli nie chcesz zrestrukturyzować kodu, aby uniknąć korzystania z tablicy strukturalnej, najlepiej jest użyć numeru bisect_left(a['f0'], 400.). Na moim komputerze jest on 8 razy wolniejszy niż wyszukiwane na ciągłej tablicy, ale jest 1000 razy szybszy niż wyszukiwanie w niezespolonej tablicy.

In [5]: a = np.arange((6e6)).view([('f0', float), ('f1', float)]) 

In [6]: timeit a['f0'].searchsorted(400.) 
10 loops, best of 3: 51.1 ms per loop 

In [7]: timeit a['f0'].copy() 
10 loops, best of 3: 51 ms per loop 

In [8]: timeit bisect_left(a['f0'], 400.) 
10000 loops, best of 3: 52.8 us per loop 

In [9]: f0 = a['f0'].copy() 

In [10]: timeit f0.searchsorted(400.) 
100000 loops, best of 3: 7.85 us per loop 

Answer 2

Oto kod ilustrujący rozmiar problemu (stan na 11 maja 2015 r.) I sposób jego "naprawienia".

import numpy as np 
import bisect 
import timeit 
from random import randint 

dtype = np.dtype([ ('pos','<I'),('sig','<H') ])    # my data is unsigned 32bit, and unsigned 16bit 
data1 = np.fromfile('./all2/840d.0a9b45e8c5344abf6ac761017e93b5bb.2.1bp.binary', dtype) 

dtype2 = np.dtype([('pos',np.uint32),('sig',np.uint32)]) # convert data to both unsigned 32bit 
data2 = data1.astype(dtype2) 

data3 = data2.view(('uint32', len(data2.dtype.names))) # convert above to a normal array (not structured array) 

print data1.dtype.descr # [('pos', '<u4'), ('sig', '<u2')] 
print data2.dtype.descr # [('pos', '<u4'), ('sig', '<u4')] 
print data3.dtype.descr # [('', '<u4')] 

print data1.nbytes # 50344494 
print data2.nbytes # 67125992 
print data3.nbytes # 67125992 

print data1['pos'].max() # 2099257376 
print data2['pos'].max() # 2099257376 
print data3[:,0].max() # 2099257376 

def b1(): return bisect.bisect_left(data1['pos'],   randint(100000000,200000000)) 
def b2(): return bisect.bisect_left(data2['pos'],   randint(100000000,200000000)) 
def b3(): return bisect.bisect_left(data3[:,0],    randint(100000000,200000000)) 
def ss1(): return np.searchsorted(data1['pos'],    randint(100000000,200000000)) 
def ss2(): return np.searchsorted(data2['pos'],    randint(100000000,200000000)) 
def ss3(): return np.searchsorted(data3[:,0],    randint(100000000,200000000)) 

def ricob1(): return bisect.bisect_left(data1['pos'], np.uint32(randint(100000000,200000000))) 
def ricob2(): return bisect.bisect_left(data2['pos'], np.uint32(randint(100000000,200000000))) 
def ricob3(): return bisect.bisect_left(data3[:,0], np.uint32(randint(100000000,200000000))) 
def ricoss1(): return np.searchsorted(data1['pos'], np.uint32(randint(100000000,200000000))) 
def ricoss2(): return np.searchsorted(data2['pos'], np.uint32(randint(100000000,200000000))) 
def ricoss3(): return np.searchsorted(data3[:,0],  np.uint32(randint(100000000,200000000))) 

print timeit.timeit(b1,number=300) # 0.0085117816925 
print timeit.timeit(b2,number=300) # 0.00826191902161 
print timeit.timeit(b3,number=300) # 0.00828003883362 
print timeit.timeit(ss1,number=300) # 6.57477498055 
print timeit.timeit(ss2,number=300) # 5.95308589935 
print timeit.timeit(ss3,number=300) # 5.92483091354 

print timeit.timeit(ricob1,number=300) # 0.00120902061462 
print timeit.timeit(ricob2,number=300) # 0.00120401382446 
print timeit.timeit(ricob3,number=300) # 0.00120711326599 
print timeit.timeit(ricoss1,number=300) # 4.39265394211 
print timeit.timeit(ricoss2,number=300) # 0.00116586685181 
print timeit.timeit(ricoss3,number=300) # 0.00108909606934 

Aktualizacja! Dzięki komentarzom Rico wygląda na to, że ustawiasz typ dla numeru, który chcesz przeszukać/bisect jest naprawdę importowany! Jednak na tablicy strukturalnej z 32-bitowym i 16-bitowym numerem int, nadal jest wolna (choć nie jest tak blisko jak poprzednio).