Was ist der schnellste Weg, um Gruppennamen von Numpy-Arrays Indizes zuzuordnen?

9

Ich arbeite mit 3D Pointcloud von Lidar. Die Punkte werden durch ein numpy-Array angegeben, das folgendermaßen aussieht:

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]])

Ich möchte meine Daten in Würfel mit einer Größe gruppieren, 50*50*50damit jeder Würfel einige Hash-Indizes und Numpy-Indizes meiner pointsenthaltenen Würfel beibehält . Um eine Aufteilung zu erhalten, ordne ich folgende cubes = points \\ 50Ausgänge zu:

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]])

Meine gewünschte Ausgabe sieht folgendermaßen aus:

{(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]}

Meine echte Punktwolke enthält bis zu einige hundert Millionen 3D-Punkte. Was ist der schnellste Weg, um diese Art der Gruppierung durchzuführen?

Ich habe die meisten verschiedenen Lösungen ausprobiert. Hier ist ein Vergleich der Zeitaufnahme unter der Annahme, dass die Größe der Punkte etwa 20 Millionen und die Größe der einzelnen Würfel etwa 1 Million beträgt:

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

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

Defauldict [elem.tobytes () oder Tupel -> Liste]

#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

Pandas + Dimensionsreduktion [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

Es ist möglich , zum Download - cubes.npzDatei hier und einen Befehl

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

um die Leistungszeit zu überprüfen.

mathfux
quelle
Haben Sie in Ihrem Ergebnis immer die gleiche Anzahl von Indizes in jeder Liste?
Mykola Zotko
Ja, es ist immer dasselbe: 983234 verschiedene Würfel für alle oben genannten Lösungen.
Mathfux
1
Es ist unwahrscheinlich, dass eine so einfache Pandas-Lösung durch einen einfachen Ansatz übertroffen wird, da große Anstrengungen unternommen wurden, um sie zu optimieren. Ein Cython-basierter Ansatz könnte sich ihm wahrscheinlich nähern, aber ich bezweifle, dass er ihn übertreffen würde.
Norok2
1
@mathfux Müssen Sie die endgültige Ausgabe als Wörterbuch haben oder wäre es in Ordnung, die Gruppen und ihre Indizes als zwei Ausgaben zu haben?
Divakar
@ norok2 numpy_indexednähert sich auch nur. Ich denke es ist richtig. Ich verwende pandasderzeit für meine Klassifizierungsprozesse.
Mathfux

Antworten:

6

Konstante Anzahl von Indizes pro Gruppe

Ansatz Nr. 1

Wir können auftreten dimensionality-reduction , um cubesauf ein 1D-Array zu reduzieren . Dies basiert auf einer Abbildung der gegebenen Würfeldaten auf ein n-dim-Gitter, um die im Detail diskutierten linearen Indexäquivalente zu berechnen here. Basierend auf der Eindeutigkeit dieser linearen Indizes können wir dann eindeutige Gruppen und ihre entsprechenden Indizes trennen. Wenn wir diesen Strategien folgen, hätten wir eine Lösung wie diese -

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))

Alternative Nr. 1: Wenn die ganzzahligen Werte in cubeszu groß sind, möchten wir dies möglicherweise so tun dimensionality-reduction, dass die Dimensionen mit kürzerer Ausdehnung als primäre Achsen ausgewählt werden. Daher können wir für diese Fälle den Reduktionsschritt modifizieren, um Folgendes zu erhalten c1D:

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

Ansatz Nr. 2

Als nächstes können wir verwenden Cython-powered kd-tree schnelle Suche nach dem nächsten Nachbarn verwenden, um die nächsten Nachbarindizes zu erhalten und damit unseren Fall wie folgt zu lösen:

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]

Allgemeiner Fall: Variable Anzahl von Indizes pro Gruppe

Wir werden die argsort-basierte Methode mit einigen Aufteilungen erweitern, um die gewünschte Ausgabe zu erhalten.

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))

Verwenden von 1D-Versionen von Gruppen cubesals Schlüssel

Wir werden die zuvor aufgeführte Methode um die Gruppen von cubesals Schlüssel erweitern, um den Prozess der Wörterbucherstellung zu vereinfachen und sie auch so effizient zu gestalten.

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

Als nächstes werden wir nutzen numba Paket verwenden, um zu iterieren und zur endgültigen Ausgabe des Hashable-Wörterbuchs zu gelangen. Dazu gibt es zwei Lösungen: Eine, bei der die Schlüssel und Werte separat verwendet werden, numbaund der Hauptaufruf werden komprimiert und in Dikt konvertiert, während die andere einen numba-supportedDiktat-Typ erstellt und somit keine zusätzliche Arbeit für die Hauptaufruffunktion erforderlich ist .

Wir hätten also die erste numbaLösung:

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

Und zweite numbaLösung als:

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

Timings mit cubes.npzDaten -

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)

Alternative Nr. 1: Wir können eine weitere Beschleunigung erzielen, indem numexprgroße Arrays c1Dwie folgt berechnet werden -

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)

Dies gilt für alle erforderlichen Stellen c1D.

Divakar
quelle
Vielen Dank für die Antwort! Ich habe nicht erwartet, dass die Verwendung von cKDTree hier möglich ist. Es gibt jedoch immer noch einige Probleme mit Ihrem # Approach1. Die Länge der Ausgabe beträgt nur 915791. Ich denke, dies ist eine Art Konflikt zwischen dtypes int32undint64
mathfux
@mathfux Ich gehe davon aus, number of indices per group would be a constant numberdass ich die Kommentare gesammelt habe. Wäre das eine sichere Annahme? Testen Sie auch cubes.npzdie Ausgabe von 915791?
Divakar
Ja, ich will. Ich habe die Anzahl der Indizes pro Gruppe nicht getestet, da die Reihenfolge der Gruppennamen unterschiedlich sein kann. Ich teste die Länge des Wörterbuchs der Ausgabe cubes.npznur von und es war 983234für die anderen Ansätze, die ich vorgeschlagen habe.
Mathfux
1
@mathfux Suchen Sie Approach #3 nach dem generischen Fall einer variablen Anzahl von Indizes.
Divakar
1
@mathfux Ja, dieser Ausgleich ist im Allgemeinen erforderlich, wenn das Minimum weniger als 0 beträgt. Guter Fang für die Präzision!
Divakar
5

Sie können einfach iterieren und den Index jedes Elements zur entsprechenden Liste hinzufügen.

from collections import defaultdict

res = defaultdict(list)

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

Die Laufzeit kann weiter verbessert werden, indem tobytes () verwendet wird, anstatt den Schlüssel in ein Tupel zu konvertieren.

ABC
quelle
Ich versuche gerade, die Leistungszeit zu überprüfen (für 20 Millionen Punkte). Es scheint, dass meine Lösung zeitlich effizienter ist, da Iterationen vermieden werden. Ich stimme zu, der Speicherverbrauch ist enorm.
Mathfux
Ein anderer Vorschlag res[tuple(elem)].append(idx)dauerte 50 Sekunden gegenüber seiner Ausgabe, res[elem[0], elem[1], elem[2]].append(idx)die 30 Sekunden dauerte.
Mathfux
3

Sie könnten Cython verwenden:

%%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

Aber es wird Sie nicht schneller machen als das, was Pandas tut, obwohl es danach das schnellste ist (und vielleicht die numpy_indexbasierte Lösung) und nicht mit der Speicherstrafe verbunden ist. Eine Sammlung der bisher vorgeschlagenen Vorschläge finden Sie hier .

Auf dem OP-Computer sollte die Ausführungszeit ungefähr 12 Sekunden betragen.

norok2
quelle
1
Vielen Dank, ich werde es später testen.
Mathfux