Jaki jest najszybszy sposób mapowania nazw grup tablic numpy na indeksy?

9

Pracuję z 3D pointcloud firmy Lidar. Punkty są przyznawane przez tablicę numpy, która wygląda następująco:

points = np.array([[61651921, 416326074, 39805], [61605255, 416360555, 41124], [61664810, 416313743, 39900], [61664837, 416313749, 39910], [61674456, 416316663, 39503], [61651933, 416326074, 39802], [61679969, 416318049, 39500], [61674494, 416316677, 39508], [61651908, 416326079, 39800], [61651908, 416326087, 39802], [61664845, 416313738, 39913], [61674480, 416316668, 39503], [61679996, 416318047, 39510], [61605290, 416360572, 41118], [61605270, 416360565, 41122], [61683939, 416313004, 41052], [61683936, 416313033, 41060], [61679976, 416318044, 39509], [61605279, 416360555, 41109], [61664837, 416313739, 39915], [61674487, 416316666, 39505], [61679961, 416318035, 39503], [61683943, 416313004, 41054], [61683930, 416313042, 41059]])

Chciałbym, aby moje dane pogrupowane w kostkę wielkości 50*50*50tak, że każda kostka zachowuje pewną hashable indeksu i NumPy indeksy z moich pointszawiera . Aby uzyskać podział, przypisuję cubes = points \\ 50które wyjścia do:

cubes = np.array([[1233038, 8326521, 796], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233599, 8326360, 790], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233038, 8326521, 796], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1232105, 8327211, 822], [1232105, 8327211, 822], [1233678, 8326260, 821], [1233678, 8326260, 821], [1233599, 8326360, 790], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1233678, 8326260, 821], [1233678, 8326260, 821]])

Moje pożądane wyniki wyglądają następująco:

{(1232105, 8327211, 822): [1, 13, 14, 18]), 
(1233038, 8326521, 796): [0, 5, 8, 9], 
(1233296, 8326274, 798): [2, 3, 10, 19], 
(1233489, 8326333, 790): [4, 7, 11, 20], 
(1233599, 8326360, 790): [6, 12, 17, 21], 
(1233678, 8326260, 821): [15, 16, 22, 23]}

Moja prawdziwa chmura punktów zawiera do kilkuset milionów punktów 3D. Jaki jest najszybszy sposób na grupowanie tego rodzaju?

Wypróbowałem większość różnych rozwiązań. Oto porównanie obliczania czasu przy założeniu, że wielkość punktów wynosi około 20 milionów, a wielkość odrębnych kostek wynosi około 1 miliona:

Pandas [tuple (elem) -> np.array (dtype = int64)]

import pandas as pd
print(pd.DataFrame(cubes).groupby([0,1,2]).indices)
#takes 9sec

Defauldict [elem.tobytes () or tuple -> list]

#thanks @abc:
result = defaultdict(list)
for idx, elem in enumerate(cubes):
    result[elem.tobytes()].append(idx) # takes 20.5sec
    # result[elem[0], elem[1], elem[2]].append(idx) #takes 27sec
    # result[tuple(elem)].append(idx) # takes 50sec

numpy_indexed [int -> np.array]

# thanks @Eelco Hoogendoorn for his library
values = npi.group_by(cubes).split(np.arange(len(cubes)))
result = dict(enumerate(values))
# takes 9.8sec

Pandy + redukcja wymiarowości [int -> np.array (dtype = int64)]

# thanks @Divakar for showing numexpr library:
import numexpr as ne
def dimensionality_reduction(cubes):
    #cubes = cubes - np.min(cubes, axis=0) #in case some coords are negative 
    cubes = cubes.astype(np.int64)
    s0, s1 = cubes[:,0].max()+1, cubes[:,1].max()+1
    d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
    c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)
    return c1D
cubes = dimensionality_reduction(cubes)
result = pd.DataFrame(cubes).groupby([0]).indices
# takes 2.5 seconds

Można pobrać cubes.npzplik tutaj i użyć polecenia

cubes = np.load('cubes.npz')['array']

sprawdzić czas działania.

mathfux
źródło
Czy zawsze masz taką samą liczbę indeksów na każdej liście w wyniku?
Mykoła Zotko,
Tak, zawsze jest tak samo: 983234 odrębne kostki dla wszystkich wyżej wymienionych rozwiązań.
mathfux,
1
Jest mało prawdopodobne, aby takie proste rozwiązanie Pandas zostało pokonane przez proste podejście, ponieważ wiele wysiłku włożono w jego optymalizację. Podejście oparte na Cython prawdopodobnie mogłoby do niego podejść, ale wątpię, czy przyniosłoby to lepsze rezultaty.
norok2
1
@mathfux Czy musisz mieć końcowy wynik jako słownik, czy byłoby dobrze, gdyby grupy i ich indeksy były dwoma wyjściami?
Divakar
@ norok2 też numpy_indexeddo tego podchodzi. Myślę, że to prawda. pandasObecnie używam do swoich procesów klasyfikacji.
mathfux,

Odpowiedzi:

6

Stała liczba wskaźników na grupę

Podejście nr 1

Możemy wykonać dimensionality-reductionredukcję cubesdo tablicy 1D. Jest to oparte na odwzorowaniu danych danych kostek na siatkę n-dim w celu obliczenia ekwiwalentów indeksu liniowego, omówionych szczegółowo here. Następnie, w oparciu o wyjątkowość tych liniowych wskaźników, możemy segregować unikalne grupy i odpowiadające im wskaźniki. Dlatego zgodnie z tymi strategiami mielibyśmy jedno rozwiązanie, takie jak -

N = 4 # number of indices per group
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
indices = sidx.reshape(-1,N)
unq_groups = cubes[indices[:,0]]

# If you need in a zipped dictionary format
out = dict(zip(map(tuple,unq_groups), indices))

Alternatywa # 1: Jeśli wartości całkowite w cubessą zbyt duże, możemy chcieć zrobić dimensionality-reductiontak, aby wymiary o mniejszym zasięgu były wybierane jako osie podstawowe. Dlatego w tych przypadkach możemy zmodyfikować krok redukcji, aby uzyskać c1D-

s1,s2 = cubes[:,:2].max(0)+1
s = np.r_[s2,1,s1*s2]
c1D = cubes.dot(s)

Podejście nr 2

Następnie możemy użyć Cython-powered kd-treeszybkiego wyszukiwania najbliższego sąsiada, aby uzyskać indeksy najbliższego sąsiedztwa, a tym samym rozwiązać nasz przypadek w ten sposób -

from scipy.spatial import cKDTree

idx = cKDTree(cubes).query(cubes, k=N)[1] # N = 4 as discussed earlier
I = idx[:,0].argsort().reshape(-1,N)[:,0]
unq_groups,indices = cubes[I],idx[I]

Przypadek ogólny: Zmienna liczba wskaźników na grupę

Rozszerzymy metodę opartą na argsort z pewnym podziałem, aby uzyskać pożądaną wydajność, tak jak -

c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)

sidx = c1D.argsort()
c1Ds = c1D[sidx]
split_idx = np.flatnonzero(np.r_[True,c1Ds[:-1]!=c1Ds[1:],True])
grps = cubes[sidx[split_idx[:-1]]]

indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
# If needed as dict o/p
out = dict(zip(map(tuple,grps), indices))

Korzystanie z wersji 1D grup cubesjako kluczy

Rozszerzymy wcześniej wymienioną metodę o grupy cubeskluczy as, aby uprościć proces tworzenia słownika, a także zwiększyć jego efektywność, tak jak to -

def numpy1(cubes):
    c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)        
    sidx = c1D.argsort()
    c1Ds = c1D[sidx]
    mask = np.r_[True,c1Ds[:-1]!=c1Ds[1:],True]
    split_idx = np.flatnonzero(mask)
    indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
    out = dict(zip(c1Ds[mask[:-1]],indices))
    return out

Następnie wykorzystamy numbapakiet do iteracji i przejdziemy do ostatecznego wyjścia słownika haszującego. Idąc za tym, będą dwa rozwiązania - jedno, które pobiera klucze i wartości oddzielnie, numbaa główne wywołanie zostanie skompresowane i zamienione na dyktowanie, a drugie stworzynumba-supported typ dykta, a zatem główna funkcja wywoływania nie wymaga dodatkowej pracy .

Zatem mielibyśmy pierwsze numbarozwiązanie:

from numba import  njit

@njit
def _numba1(sidx, c1D):
    out = []
    n = len(sidx)
    start = 0
    grpID = []
    for i in range(1,n):
        if c1D[sidx[i]]!=c1D[sidx[i-1]]:
            out.append(sidx[start:i])
            grpID.append(c1D[sidx[start]])
            start = i
    out.append(sidx[start:])
    grpID.append(c1D[sidx[start]])
    return grpID,out

def numba1(cubes):
    c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
    sidx = c1D.argsort()
    out = dict(zip(*_numba1(sidx, c1D)))
    return out

I drugie numbarozwiązanie jako:

from numba import types
from numba.typed import Dict

int_array = types.int64[:]

@njit
def _numba2(sidx, c1D):
    n = len(sidx)
    start = 0
    outt = Dict.empty(
        key_type=types.int64,
        value_type=int_array,
    )
    for i in range(1,n):
        if c1D[sidx[i]]!=c1D[sidx[i-1]]:
            outt[c1D[sidx[start]]] = sidx[start:i]
            start = i
    outt[c1D[sidx[start]]] = sidx[start:]
    return outt

def numba2(cubes):
    c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)    
    sidx = c1D.argsort()
    out = _numba2(sidx, c1D)
    return out

Czasy z cubes.npzdanymi -

In [4]: cubes = np.load('cubes.npz')['array']

In [5]: %timeit numpy1(cubes)
   ...: %timeit numba1(cubes)
   ...: %timeit numba2(cubes)
2.38 s ± 14.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2.13 s ± 25.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.8 s ± 5.95 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Alternatywa # 1: Możemy osiągnąć dalsze przyspieszenie przy numexprobliczaniu dużych tablic c1D, tak jak -

import numexpr as ne

s0,s1 = cubes[:,0].max()+1,cubes[:,1].max()+1
d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)

Miałoby to zastosowanie we wszystkich miejscach, które tego wymagają c1D.

Divakar
źródło
Wielkie dzięki za odpowiedź! Nie spodziewałem się, że użycie cKDTree jest możliwe tutaj. Nadal jednak występują problemy z Twoim # podejściem1. Długość wyjścia wynosi tylko 915791. Myślę, że to jakiś konflikt między dtypes int32iint64
mathfux,
@mathfux Zakładam, number of indices per group would be a constant numberże zebrałem komentarze. Czy byłoby to bezpieczne założenie? Czy testujesz cubes.npztakże dane wyjściowe 915791?
Divakar
Tak. Nie testowałem liczby indeksów na grupę, ponieważ kolejność nazw grup może być inna. Testuję długość słownika wyjściowego zcubes.npz tylko i dotyczy to 983234innych podejść, które zasugerowałem.
mathfux,
1
@mathfux Sprawdź Approach #3 ogólny przypadek zmiennej liczby indeksów.
Divakar
1
@mathfux Tak, to przesunięcie jest ogólnie potrzebne, jeśli minimum jest mniejsze niż 0. Dobry chwyt na precyzji!
Divakar
5

Możesz po prostu iterować i dodać indeks każdego elementu do odpowiedniej listy.

from collections import defaultdict

res = defaultdict(list)

for idx, elem in enumerate(cubes):
    #res[tuple(elem)].append(idx)
    res[elem.tobytes()].append(idx)

Środowisko wykonawcze można dodatkowo ulepszyć za pomocą tobytes () zamiast konwertowania klucza na krotkę.

ABC
źródło
W tej chwili próbuję dokonać przeglądu czasu wykonania (za 20 milionów punktów). Wydaje się, że moje rozwiązanie jest bardziej wydajne pod względem czasu, ponieważ unika się iteracji. Zgadzam się, zużycie pamięci jest ogromne.
mathfux,
kolejna propozycja res[tuple(elem)].append(idx)zajęła 50 sekund w porównaniu do jej edycji, res[elem[0], elem[1], elem[2]].append(idx)która zajęła 30 sekund.
mathfux,
3

Możesz użyć Cython:

%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True

import math
import cython as cy

cimport numpy as cnp


cpdef groupby_index_dict_cy(cnp.int32_t[:, :] arr):
    cdef cy.size_t size = len(arr)
    result = {}
    for i in range(size):
        key = arr[i, 0], arr[i, 1], arr[i, 2]
        if key in result:
            result[key].append(i)
        else:
            result[key] = [i]
    return result

ale nie sprawi, że będziesz szybszy niż to, co robi Panda, chociaż po tym jest najszybszy (i być może numpy_index oparte na rozwiązaniu) i nie wiąże się z karą pamięci. Zbiór tego, co dotychczas zaproponowano, znajduje się tutaj .

W maszynie OP, która powinna zbliżyć się do ~ 12 sekund czasu wykonania.

norok2
źródło
1
Wielkie dzięki, przetestuję to później.
mathfux