¿Cuál es la forma más rápida de asignar nombres de grupos de matrices numpy a índices?

9

Estoy trabajando con 3D pointcloud de Lidar. Los puntos están dados por una matriz numpy que se ve así:

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

Me gustaría mantener mis datos agrupados en cubos de tamaño 50*50*50para que cada cubo conserve algún índice hashable e índices numpy de mi pointsque contiene . Para dividir, asigno cubes = points \\ 50las salidas a:

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

Mi salida deseada se ve así:

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

Mi nube de puntos real contiene hasta unos pocos cientos de millones de puntos 3D. ¿Cuál es la forma más rápida de hacer este tipo de agrupación?

He probado la mayoría de varias soluciones. Aquí hay una comparación del consumo de tiempo asumiendo que el tamaño de los puntos es de alrededor de 20 millones y el tamaño de los cubos distintos es de alrededor de 1 millón:

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

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

Por defecto [elem.tobytes () o tupla -> lista]

#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 + reducción de dimensionalidad [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 posible descargar el cubes.npzarchivo aquí y usar un comando

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

para verificar el tiempo de rendimiento.

mathfux
fuente
¿Siempre tiene el mismo número de índices en cada lista en su resultado?
Mykola Zotko
Sí, siempre es lo mismo: 983234 cubos distintos para todas las soluciones mencionadas anteriormente.
mathfux
1
Es poco probable que una solución Pandas tan simple sea superada por un enfoque simple, ya que se ha dedicado un gran esfuerzo para optimizarla. Un enfoque basado en Cython probablemente podría abordarlo, pero dudo que lo supere.
norok2
1
@mathfux ¿Tiene que tener la salida final como diccionario o estaría bien tener los grupos y sus índices como dos salidas?
Divakar
@ norok2 numpy_indexedsolo se acerca a él también. Supongo que es correcto. Utilizo pandaspara mis procesos de clasificación actualmente.
mathfux

Respuestas:

6

Número constante de índices por grupo

Enfoque n. ° 1

Podemos realizar dimensionality-reductionpara reducir cubesa una matriz 1D. Esto se basa en un mapeo de los datos de cubos dados en una cuadrícula n-dim para calcular los equivalentes de índice lineal, discutidos en detalle here. Luego, en función de la unicidad de esos índices lineales, podemos segregar grupos únicos y sus índices correspondientes. Por lo tanto, siguiendo esas estrategias, tendríamos una solución, así:

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

Alternativa n. ° 1: si los valores enteros en cubesson demasiado grandes, es posible que deseemos hacer dimensionality-reductionque las dimensiones con menor extensión se elijan como ejes primarios. Por lo tanto, para esos casos, podemos modificar el paso de reducción para obtener c1D, así:

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

Enfoque n. ° 2

A continuación, podemos usar la Cython-powered kd-treebúsqueda rápida del vecino más cercano para obtener los índices vecinos más cercanos y, por lo tanto, resolver nuestro caso así:

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]

Caso genérico: número variable de índices por grupo

Extendiremos el método basado en argsort con un poco de división para obtener el resultado deseado, así:

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

Usando versiones 1D de grupos de cubescomo claves

Ampliaremos el método enumerado anteriormente con los grupos de cubesclaves para simplificar el proceso de creación de diccionarios y también hacerlo eficiente con él, de esta manera:

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

A continuación, utilizaremos el numbapaquete para iterar y llegar a la salida final del diccionario hashable. Para ello, habría dos soluciones: una que obtiene las claves y los valores por separado numbay la llamada principal se comprimirá y convertirá a dict, mientras que la otra creará un numba-supportedtipo de dict y, por lo tanto, no se requiere trabajo adicional por parte de la función de llamada principal .

Por lo tanto, tendríamos la primera numbasolución:

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

Y segunda numbasolución como:

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

Tiempos con cubes.npzdatos -

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)

Alternativa n. ° 1: Podemos lograr una mayor velocidad con computadoras numexprgrandes para calcular c1D, así:

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)

Esto sería aplicable en todos los lugares que lo requieran c1D.

Divakar
fuente
Muchas gracias por la respuesta! No esperaba que el uso de cKDTree sea posible aquí. Sin embargo, todavía hay algunos problemas con su # Enfoque1. La longitud de salida es solo 915791. Supongo que esto es algún tipo de conflicto entre dtypes int32yint64
mathfux
@mathfux Supongo number of indices per group would be a constant numberque reuní los comentarios. ¿Sería una suposición segura? Además, ¿estás probando cubes.npzla salida de 915791?
Divakar
Sí. No probé el número de índices por grupo porque el orden de los nombres de los grupos puede ser diferente. Solo pruebo la longitud del diccionario de salida cubes.npzy fue 983234para los otros enfoques que sugerí.
mathfux
1
@mathfux Echa un vistazo a Approach #3 ese caso genérico de número variable de índices.
Divakar
1
@mathfux Sí, esa compensación es necesaria generalmente si el mínimo es inferior a 0. ¡Buena captura de precisión!
Divakar
5

Puede iterar y agregar el índice de cada elemento a la lista correspondiente.

from collections import defaultdict

res = defaultdict(list)

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

El tiempo de ejecución se puede mejorar aún más mediante el uso de tobytes () en lugar de convertir la clave en una tupla.

a B C
fuente
Estoy tratando de hacer una revisión del tiempo de rendimiento en este momento (por 20 millones de puntos). Parece que mi solución es más eficiente en términos de tiempo porque se evita la iteración. Estoy de acuerdo, el consumo de memoria es enorme.
mathfux
Otra propuesta res[tuple(elem)].append(idx)tardó 50 segundos en comparación con su edición, res[elem[0], elem[1], elem[2]].append(idx)que tardó 30 segundos.
mathfux
3

Puedes usar 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

pero no lo hará más rápido que lo que hace Pandas, aunque es el más rápido después de eso (y tal vez la numpy_indexsolución basada), y no viene con la penalización de la memoria. Una colección de lo que se ha propuesto hasta ahora está aquí .

En la máquina de OP que debería acercarse a ~ 12 segundos de tiempo de ejecución.

norok2
fuente
1
Muchas gracias, lo probaré más tarde.
mathfux