Numpy: znajdź szybko pierwszy indeks wartości

Jak mogę znaleźć indeks pierwszego wystąpienia liczby w tablicy Numpy? Szybkość jest dla mnie ważna. Nie interesują mnie następujące odpowiedzi, ponieważ skanują całą tablicę i nie zatrzymują się, gdy znajdą pierwsze wystąpienie:

itemindex = numpy.where(array==item)[0][0]
nonzero(array == item)[0][0]

Uwaga 1: żadna z odpowiedzi z tego pytania nie wydaje się trafna. Czy istnieje funkcja Numpy zwracająca pierwszy indeks czegoś w tablicy?

Uwaga 2: użycie metody skompilowanej w języku C jest preferowane niż pętla Pythona.

Pojawiła się prośba o funkcję zaplanowaną dla Numpy 2.0.0: https://github.com/numpy/numpy/issues/2269

Chociaż dla ciebie jest o wiele za późno, ale do wykorzystania w przyszłości: użycie numba ( 1 ) jest najłatwiejszym sposobem, dopóki numpy go nie zaimplementuje. Jeśli używasz dystrybucji anaconda python, powinna być już zainstalowana. Kod zostanie skompilowany, więc będzie szybki.

@jit(nopython=True)
def find_first(item, vec):
    """return the index of the first occurence of item in vec"""
    for i in xrange(len(vec)):
        if item == vec[i]:
            return i
    return -1

i wtedy:

>>> a = array([1,7,8,32])
>>> find_first(8,a)
2

Zrobiłem benchmark dla kilku metod:

• argwhere
• nonzero jak w pytaniu
• .tostring() jak w odpowiedzi @Rob Reilink
• pętla Pythona
• Pętla Fortran

Python i Fortran kodu są dostępne. Pominąłem te mało obiecujące, jak konwersja na listę.

Wyniki w skali logarytmicznej. Oś X to pozycja igły (znalezienie, czy jest dalej w dół tablicy, zajmuje więcej czasu); ostatnia wartość to igła, której nie ma w tablicy. Oś Y to czas, aby ją znaleźć.

Tablica miała 1 milion elementów, a testy były wykonywane 100 razy. Wyniki wciąż nieco się wahają, ale trend jakościowy jest jasny: Python i f2py kończą pracę przy pierwszym elemencie, więc skalują się inaczej. Python działa zbyt wolno, jeśli igła nie znajduje się w pierwszym 1%, podczas gdy f2pyjest szybki (ale musisz go skompilować).

Podsumowując, f2py to najszybsze rozwiązanie , zwłaszcza jeśli igła pojawia się dość wcześnie.

Nie jest wbudowany, co jest denerwujące, ale to tylko 2 minuty pracy. Dodaj to do pliku o nazwie search.f90:

subroutine find_first(needle, haystack, haystack_length, index)
    implicit none
    integer, intent(in) :: needle
    integer, intent(in) :: haystack_length
    integer, intent(in), dimension(haystack_length) :: haystack
!f2py intent(inplace) haystack
    integer, intent(out) :: index
    integer :: k
    index = -1
    do k = 1, haystack_length
        if (haystack(k)==needle) then
            index = k - 1
            exit
        endif
    enddo
end

Jeśli szukasz czegoś innego niż integer, po prostu zmień typ. Następnie skompiluj używając:

f2py -c -m search search.f90

po czym możesz zrobić (z Pythona):

import search
print(search.find_first.__doc__)
a = search.find_first(your_int_needle, your_int_array)

Możesz przekonwertować tablicę logiczną na ciąg znaków w języku Python za pomocą, array.tostring()a następnie używając metody find ():

(array==item).tostring().find('\x01')

Wiąże się to jednak z kopiowaniem danych, ponieważ ciągi znaków w Pythonie muszą być niezmienne. Zaletą jest to, że możesz również wyszukać np. Zbocze narastające, znajdując\x00\x01

W przypadku tablic posortowanych np.searchsorteddziała.

Myślę, że napotkałeś problem, w którym inna metoda i pewna wiedza a priori o tablicy naprawdę by pomogła. To coś, w przypadku którego istnieje prawdopodobieństwo X, że znajdziesz odpowiedź w pierwszych Y procentach danych. Rozdzielenie problemu z nadzieją na szczęście, a następnie zrobienie tego w Pythonie z zagnieżdżonym zrozumieniem list lub czymś podobnym.

Pisanie funkcji C wykonującej tę brutalną siłę nie jest zbyt trudne przy użyciu ctypów .

Kod C, który zhakowałem razem (index.c):

long index(long val, long *data, long length){
    long ans, i;
    for(i=0;i<length;i++){
        if (data[i] == val)
            return(i);
    }
    return(-999);
}

i python:

# to compile (mac)
# gcc -shared index.c -o index.dylib
import ctypes
lib = ctypes.CDLL('index.dylib')
lib.index.restype = ctypes.c_long
lib.index.argtypes = (ctypes.c_long, ctypes.POINTER(ctypes.c_long), ctypes.c_long)

import numpy as np
np.random.seed(8675309)
a = np.random.random_integers(0, 100, 10000)
print lib.index(57, a.ctypes.data_as(ctypes.POINTER(ctypes.c_long)), len(a))

i otrzymuję 92.

Zapakuj Pythona w odpowiednią funkcję i gotowe.

Wersja C jest dużo (~ 20x) szybsza dla tego materiału siewnego (ostrzeżenie, że nie jestem dobry z czasem)

import timeit
t = timeit.Timer('np.where(a==57)[0][0]', 'import numpy as np; np.random.seed(1); a = np.random.random_integers(0, 1000000, 10000000)')
t.timeit(100)/100
# 0.09761879920959472
t2 = timeit.Timer('lib.index(57, a.ctypes.data_as(ctypes.POINTER(ctypes.c_long)), len(a))', 'import numpy as np; np.random.seed(1); a = np.random.random_integers(0, 1000000, 10000000); import ctypes; lib = ctypes.CDLL("index.dylib"); lib.index.restype = ctypes.c_long; lib.index.argtypes = (ctypes.c_long, ctypes.POINTER(ctypes.c_long), ctypes.c_long) ')
t2.timeit(100)/100
# 0.005288000106811523

@tal przedstawił już numbafunkcję znajdującą pierwszy indeks, ale działa ona tylko dla tablic 1D. Za pomocą np.ndenumeratemożesz również znaleźć pierwszy indeks w tablicy dowolnie wymiarowej:

from numba import njit
import numpy as np

@njit
def index(array, item):
    for idx, val in np.ndenumerate(array):
        if val == item:
            return idx
    return None

Przykładowy przypadek:

>>> arr = np.arange(9).reshape(3,3)
>>> index(arr, 3)
(1, 0)

Czasy pokazują, że jego wydajność jest podobna do rozwiązania tals :

arr = np.arange(100000)
%timeit index(arr, 5)           # 1000000 loops, best of 3: 1.88 µs per loop
%timeit find_first(5, arr)      # 1000000 loops, best of 3: 1.7 µs per loop

%timeit index(arr, 99999)       # 10000 loops, best of 3: 118 µs per loop
%timeit find_first(99999, arr)  # 10000 loops, best of 3: 96 µs per loop

Jeśli Twoja lista jest posortowana , możesz uzyskać bardzo szybkie przeszukiwanie indeksu za pomocą pakietu „bisect”. To O (log (n)) zamiast O (n).

bisect.bisect(a, x)

znajduje x w tablicy a, zdecydowanie szybciej w posortowanym przypadku niż jakakolwiek procedura C przeglądająca wszystkie pierwsze elementy (dla wystarczająco długich list).

Warto czasem wiedzieć.

O ile wiem tylko np.any i np. Wszystkie na tablicach boolowskich są zwarte.

W twoim przypadku numpy musi dwukrotnie przejść przez całą tablicę, raz, aby utworzyć warunek boolowski, a drugi raz, aby znaleźć indeksy.

Moją rekomendacją w tym przypadku byłoby użycie cython. Myślę, że dostosowanie przykładu do tego przypadku powinno być łatwe, zwłaszcza jeśli nie potrzebujesz dużej elastyczności dla różnych typów i kształtów.

Potrzebowałem tego do mojej pracy, więc nauczyłem się interfejsu C w języku Python i Numpy oraz napisałem własny. http://pastebin.com/GtcXuLyd Dotyczy tylko tablic 1-D, ale działa dla większości typów danych (int, float lub strings), a testy wykazały, że jest około 20 razy szybszy niż oczekiwane podejście w czystym Pythonie. tępy.

Ten problem można skutecznie rozwiązać w czystym numpy, przetwarzając tablicę w kawałkach:

def find_first(x):
    idx, step = 0, 32
    while idx < x.size:
        nz, = x[idx: idx + step].nonzero()
        if len(nz): # found non-zero, return it
            return nz[0] + idx
        # move to the next chunk, increase step
        idx += step
        step = min(9600, step + step // 2)
    return -1

Tablica jest przetwarzana w porcji o rozmiarze step. Im stepdłuższy krok, tym szybsze jest przetwarzanie zerowanej tablicy (najgorszy przypadek). Im mniejszy, tym szybsze przetwarzanie tablicy z wartością niezerową na początku. Sztuczka polega na tym, aby zacząć od małego stepi zwiększyć go wykładniczo. Ponadto nie ma potrzeby podwyższania go powyżej pewnego progu ze względu na ograniczone korzyści.

Porównałem to rozwiązanie z czystym rozwiązaniem ndarary.nonzero i numba z 10 milionami tablic typu float.

import numpy as np
from numba import jit
from timeit import timeit

def find_first(x):
    idx, step = 0, 32
    while idx < x.size:
        nz, = x[idx: idx + step].nonzero()
        if len(nz):
            return nz[0] + idx
        idx += step
        step = min(9600, step + step // 2)
    return -1

@jit(nopython=True)
def find_first_numba(vec):
    """return the index of the first occurence of item in vec"""
    for i in range(len(vec)):
        if vec[i]:
            return i
    return -1


SIZE = 10_000_000
# First only
x = np.empty(SIZE)

find_first_numba(x[:10])

print('---- FIRST ----')
x[:] = 0
x[0] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=1000), 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=1000), 'ms')

print('---- LAST ----')
x[:] = 0
x[-1] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

print('---- NONE ----')
x[:] = 0
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

print('---- ALL ----')
x[:] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

I wyniki na moim komputerze:

---- FIRST ----
ndarray.nonzero 54.733994480002366 ms
find_first 0.0013148509997336078 ms
find_first_numba 0.0002839310000126716 ms
---- LAST ----
ndarray.nonzero 54.56336712999928 ms
find_first 25.38929685000312 ms
find_first_numba 8.022820680002951 ms
---- NONE ----
ndarray.nonzero 24.13432420999925 ms
find_first 25.345200140000088 ms
find_first_numba 8.154927100003988 ms
---- ALL ----
ndarray.nonzero 55.753537260002304 ms
find_first 0.0014760300018679118 ms
find_first_numba 0.0004358099977253005 ms

Czysty ndarray.nonzerojest zdecydowanie luźniejszy. Rozwiązanie numba jest około 5 razy szybsze w najlepszym przypadku. W najgorszym przypadku jest około 3 razy szybsza.

Jeśli szukasz pierwszego niezerowego elementu, możesz użyć następującego hacka:

idx = x.view(bool).argmax() // x.itemsize
idx = idx if x[idx] else -1

Jest to bardzo szybkie rozwiązanie typu „numpy-pure”, ale zawodzi w niektórych przypadkach omówionych poniżej.

Rozwiązanie wykorzystuje fakt, że prawie cała reprezentacja zera dla typów liczbowych składa się z 0bajtów. Dotyczy to również numpy's bool. W ostatnich wersjach numpy, argmax()funkcja używa logiki zwarcia podczas przetwarzania booltypu. Rozmiar boolto 1 bajt.

Więc trzeba:

• utwórz widok tablicy jako bool. Żadna kopia nie jest tworzona
• użyj, argmax()aby znaleźć pierwszy niezerowy bajt za pomocą logiki zwarcia
• przeliczyć przesunięcie tego bajtu do indeksu pierwszego niezerowego elementu przez dzielenie całkowite (operator //) przesunięcia o rozmiar pojedynczego elementu wyrażone w bajtach ( x.itemsize)
• sprawdź, czy x[idx]faktycznie jest niezerowe, aby zidentyfikować przypadek, w którym nie ma wartości niezerowej

Zrobiłem test porównawczy w stosunku do rozwiązania numba i zbudowałem go np.nonzero.

import numpy as np
from numba import jit
from timeit import timeit

def find_first(x):
    idx = x.view(bool).argmax() // x.itemsize
    return idx if x[idx] else -1

@jit(nopython=True)
def find_first_numba(vec):
    """return the index of the first occurence of item in vec"""
    for i in range(len(vec)):
        if vec[i]:
            return i
    return -1


SIZE = 10_000_000
# First only
x = np.empty(SIZE)

find_first_numba(x[:10])

print('---- FIRST ----')
x[:] = 0
x[0] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=1000), 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=1000), 'ms')

print('---- LAST ----')
x[:] = 0
x[-1] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

print('---- NONE ----')
x[:] = 0
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

print('---- ALL ----')
x[:] = 1
print('ndarray.nonzero', timeit(lambda: x.nonzero()[0][0], number=100)*10, 'ms')
print('find_first', timeit(lambda: find_first(x), number=100)*10, 'ms')
print('find_first_numba', timeit(lambda: find_first_numba(x), number=100)*10, 'ms')

Wynik na moim komputerze to:

---- FIRST ----
ndarray.nonzero 57.63976670001284 ms
find_first 0.0010841979965334758 ms
find_first_numba 0.0002308919938514009 ms
---- LAST ----
ndarray.nonzero 58.96685277999495 ms
find_first 5.923203580023255 ms
find_first_numba 8.762269750004634 ms
---- NONE ----
ndarray.nonzero 25.13398071998381 ms
find_first 5.924289370013867 ms
find_first_numba 8.810063839919167 ms
---- ALL ----
ndarray.nonzero 55.181210660084616 ms
find_first 0.001246920000994578 ms
find_first_numba 0.00028766007744707167 ms

Rozwiązanie jest o 33% szybsze niż numba i jest „czyste-numpy”.

Wady:

• nie działa dla numpy akceptowanych typów, takich jak object
• kończy się niepowodzeniem dla ujemnego zera, które czasami pojawia się w obliczeniach floatlubdouble

Jako wieloletni użytkownik Matlaba od dłuższego czasu szukałem skutecznego rozwiązania tego problemu. Na koniec, zmotywowany dyskusjami i propozycjami w tym wątku , starałem się wymyślić rozwiązanie, które implementuje API podobne do tego, co tu sugerowano , obsługujące na razie tylko tablice 1D.

Używałbyś tego w ten sposób

import numpy as np
import utils_find_1st as utf1st
array = np.arange(100000)
item = 1000
ind = utf1st.find_1st(array, item, utf1st.cmp_larger_eq)

Obsługiwane operatory warunków to: cmp_equal, cmp_not_equal, cmp_larger, cmp_smaller, cmp_larger_eq, cmp_smaller_eq. Dla wydajności rozszerzenie jest napisane w c.

Źródło, testy porównawcze i inne szczegóły znajdziesz tutaj:

https://pypi.python.org/pypi?name=py_find_1st&:action=display

do użytku w naszym zespole (anaconda na linux i macos) Zrobiłem instalator anaconda, który upraszcza instalację, możesz go używać zgodnie z opisem tutaj

https://anaconda.org/roebel/py_find_1st

Tylko uwaga, że ​​jeśli wykonujesz sekwencję wyszukiwań, wzrost wydajności wynikający z zrobienia czegoś sprytnego, takiego jak konwersja na ciąg znaków, może zostać utracony w zewnętrznej pętli, jeśli wymiar wyszukiwania nie jest wystarczająco duży. Zobacz, jak działa iteracja find1 wykorzystująca sztuczkę konwersji ciągów zaproponowaną powyżej i find2, która używa argmax wzdłuż osi wewnętrznej (plus korekta zapewniająca, że ​​niedopasowanie zwraca wartość -1)

import numpy,time
def find1(arr,value):
    return (arr==value).tostring().find('\x01')

def find2(arr,value): #find value over inner most axis, and return array of indices to the match
    b = arr==value
    return b.argmax(axis=-1) - ~(b.any())


for size in [(1,100000000),(10000,10000),(1000000,100),(10000000,10)]:
    print(size)
    values = numpy.random.choice([0,0,0,0,0,0,0,1],size=size)
    v = values>0

    t=time.time()
    numpy.apply_along_axis(find1,-1,v,1)
    print('find1',time.time()-t)

    t=time.time()
    find2(v,1)
    print('find2',time.time()-t)

wyjścia

(1, 100000000)
('find1', 0.25300002098083496)
('find2', 0.2780001163482666)
(10000, 10000)
('find1', 0.46200013160705566)
('find2', 0.27300000190734863)
(1000000, 100)
('find1', 20.98099994659424)
('find2', 0.3040001392364502)
(10000000, 10)
('find1', 206.7590000629425)
('find2', 0.4830000400543213)

To powiedziawszy, znalezienie napisane w C byłoby co najmniej trochę szybsze niż którekolwiek z tych podejść

co powiesz na to

import numpy as np
np.amin(np.where(array==item))

Możesz ukryć swoją tablicę w a listi użyć jej index()metody:

i = list(array).index(item)

O ile mi wiadomo, jest to metoda skompilowana w C.