¿Cómo ordenar una matriz de enteros más rápido que quicksort?

La clasificación de una matriz de enteros con la ordenación rápida de Numpy se ha convertido en el cuello de botella de mi algoritmo. Desafortunadamente, numpy no tiene todavía orden de radix . Aunque la ordenación de conteo sería una de una sola línea en números:

np.repeat(np.arange(1+x.max()), np.bincount(x)) 

vea la respuesta aceptada a la pregunta ¿Cómo puedo vectorizar este tipo de conteo de Python para que sea absolutamente lo más rápido posible? Pregunta, los enteros en mi aplicación pueden correr de 0 a 2**32 .

¿Estoy atascado con quicksort?


Esta publicación fue motivada principalmente por la agrupación Numpy usando itertools.groupby pregunta de rendimiento .
También tenga en cuenta que no es simplemente correcto preguntar y responder su propia pregunta, se recomienda explícitamente.

Related of "¿Cómo ordenar una matriz de enteros más rápido que quicksort?"

No, usted no está atascado con quicksort. Podría usar, por ejemplo, integer_sort de Boost.Sort o u4_sort de usort . Al ordenar esta matriz:

 array(randint(0, high=1<<32, size=10**8), uint32) 

Obtengo los siguientes resultados:

  NumPy quicksort: 8.636 s 1.0 (línea de base)
 Boost.Sort integer_sort: 4.327 s 2.0x speedup
 usort u4_sort: 2.065 s 4.2x speedup

No saltaría a conclusiones basadas en este experimento único y usaría a usort ciegas. Me gustaría probar con mis datos reales y medir lo que sucede. Su kilometraje variará dependiendo de sus datos y de su máquina. El integer_sort en Boost.Sort tiene un amplio conjunto de opciones para el ajuste, consulte la documentación .

A continuación describo dos formas de llamar a una función C o C ++ nativa desde Python. A pesar de la larga descripción, es bastante fácil hacerlo.


Boost.Sort

Ponga estas líneas en el archivo spreadsort.cpp:

 #include  #include "boost/sort/spreadsort/spreadsort.hpp" using namespace boost::sort::spreadsort; extern "C" { void spreadsort(std::uint32_t* begin, std::size_t len) { integer_sort(begin, begin + len); } } 

Básicamente crea una instancia del integer_sort plantilla para enteros sin signo de 32 bits; la parte extern "C" garantiza el enlace C al deshabilitar la manipulación de nombres. Suponiendo que está utilizando gcc y que los archivos de boost necesarios para incluir están en el directorio /tmp/boost_1_60_0 , puede comstackrlo:

 g++ -O3 -std=c++11 -march=native -DNDEBUG -shared -fPIC -I/tmp/boost_1_60_0 spreadsort.cpp -o spreadsort.so 

Los indicadores clave son -fPIC para generar un código independiente de la posición y -shared para generar un archivo -shared de objeto compartido . (Lea los documentos de gcc para más detalles.)

Luego, envuelve la spreadsort() C ++ spreadsort() en Python usando ctypes :

 from ctypes import cdll, c_size_t, c_uint32 from numpy import uint32 from numpy.ctypeslib import ndpointer __all__ = ['integer_sort'] # In spreadsort.cpp: void spreadsort(std::uint32_t* begin, std::size_t len) lib = cdll.LoadLibrary('./spreadsort.so') sort = lib.spreadsort sort.restype = None sort.argtypes = [ndpointer(c_uint32, flags='C_CONTIGUOUS'), c_size_t] def integer_sort(arr): assert arr.dtype == uint32, 'Expected uint32, got {}'.format(arr.dtype) sort(arr, arr.size) 

Alternativamente, puedes usar cffi :

 from cffi import FFI from numpy import uint32 __all__ = ['integer_sort'] ffi = FFI() ffi.cdef('void spreadsort(uint32_t* begin, size_t len);') C = ffi.dlopen('./spreadsort.so') def integer_sort(arr): assert arr.dtype == uint32, 'Expected uint32, got {}'.format(arr.dtype) begin = ffi.cast('uint32_t*', arr.ctypes.data) C.spreadsort(begin, arr.size) 

En las cdll.LoadLibrary() y ffi.dlopen() asumí que la ruta al archivo ./spreadsort.so es ./spreadsort.so . Alternativamente, puedes escribir

 lib = cdll.LoadLibrary('spreadsort.so') 

o

 C = ffi.dlopen('spreadsort.so') 

si spreadsort.so la ruta a spreadsort.so a la LD_LIBRARY_PATH entorno LD_LIBRARY_PATH . Véase también bibliotecas compartidas .

Uso. En ambos casos, simplemente llame a la función de envoltura de Python anterior integer_sort() con su matriz numpy de enteros sin signo de 32 bits.


usort

En cuanto a u4_sort , puede comstackrlo de la siguiente manera:

 cc -DBUILDING_u4_sort -I/usr/include -I./ -I../ -I../../ -I../../../ -I../../../../ -std=c99 -fgnu89-inline -O3 -g -fPIC -shared -march=native u4_sort.c -o u4_sort.so 

Ejecute este comando en el directorio donde se encuentra el archivo u4_sort.c . (Probablemente haya una forma menos intrincada, pero no lo logré. Acabo de mirar el archivo deps.mk en el directorio de usuarios para encontrar las marcas de comstackción necesarias e incluir las rutas).

A continuación, puede ajustar la función C de la siguiente manera:

 from cffi import FFI from numpy import uint32 __all__ = ['integer_sort'] ffi = FFI() ffi.cdef('void u4_sort(unsigned* a, const long sz);') C = ffi.dlopen('u4_sort.so') def integer_sort(arr): assert arr.dtype == uint32, 'Expected uint32, got {}'.format(arr.dtype) begin = ffi.cast('unsigned*', arr.ctypes.data) C.u4_sort(begin, arr.size) 

En el código anterior, asumí que la ruta a u4_sort.so se ha añadido a la LD_LIBRARY_PATH entorno LD_LIBRARY_PATH .

Uso. Como antes con Boost.Sort, simplemente llama a la función de envoltura de Python anterior integer_sort() con su matriz numpy de enteros sin signo de 32 bits.

Una base de clasificación de radix 256 (1 byte) puede generar una matriz de conteos utilizada para determinar el tamaño del depósito en 1 paso de los datos, y luego toma 4 pases para hacer la clasificación. En mi sistema, Intel 2600K 3.4ghz, usando la comstackción de lanzamiento de Visual Studio para un progtwig C ++, toma alrededor de 1.15 segundos ordenar los 10 ^ 8 psuedo enteros sin signo de 32 bits sin signo.

Mirando el código u4_sort mencionado en la respuesta de Ali, los progtwigs son similares, pero u4_sort usa tamaños de campo de {10,11,11}, toma 3 pases para ordenar los datos y 1 paso para volver a copiar, mientras que este ejemplo usa tamaños de campo de { 8,8,8,8}, tomando 4 pases para ordenar los datos. El proceso probablemente tenga un ancho de banda de memoria limitado debido a las escrituras de acceso aleatorio, por lo que las optimizaciones en u4_sort, como las macros para el cambio, un ciclo con cambio fijo por campo, no ayudan mucho. Mis resultados son mejores probablemente debido a las diferencias del sistema y / o comstackdor. (Tenga en cuenta que u8_sort es para enteros de 64 bits).

Código de ejemplo:

 // a is input array, b is working array void RadixSort(uint32_t * a, uint32_t *b, size_t count) { size_t mIndex[4][256] = {0}; // count / index matrix size_t i,j,m,n; uint32_t u; for(i = 0; i < count; i++){ // generate histograms u = a[i]; for(j = 0; j < 4; j++){ mIndex[j][(size_t)(u & 0xff)]++; u >>= 8; } } for(j = 0; j < 4; j++){ // convert to indices m = 0; for(i = 0; i < 256; i++){ n = mIndex[j][i]; mIndex[j][i] = m; m += n; } } for(j = 0; j < 4; j++){ // radix sort for(i = 0; i < count; i++){ // sort by current lsb u = a[i]; m = (size_t)(u>>(j<<3))&0xff; b[mIndex[j][m]++] = u; } std::swap(a, b); // swap ptrs } } 

Un radix-sort con python/numba(0.23) , de acuerdo con el progtwig @rcgldr C, con multiproceso en un procesador de 2 núcleos.

Primero, el orden de radix en numba, con dos matrices globales para la eficiencia.

 from threading import Thread from pylab import * from numba import jit n=uint32(10**8) m=n//2 if 'x1' not in locals() : x1=array(randint(0,1<<16,2*n),uint16); #to avoid regeneration x2=x1.copy() x=frombuffer(x2,uint32) # randint doesn't work with 32 bits here :( y=empty_like(x) nbits=8 buffsize=1<>dec)& mask]+=1 k+=1 j=t=0 for j in range(buffsize): b=u[j] u[j]=t t+=b v=u.copy() k=0 while k>dec)&mask y[u[j]]=x[k] u[j]+=1 k+=1 x,y=y,x dec+=nbits 

Luego la paralelización, posible con la opción nogil .

 def para(nthreads=2): threads=[Thread(target=radix, args=(x[i*n//nthreads(i+1)*n//nthreads], y[i*n//nthreads:(i+1)*n//nthreads])) for i in range(nthreads)] for t in threads: t.start() for t in threads: t.join() @jit def fuse(x,y): kl=0 kr=n//2 k=0 while k 

puntos de referencia:

 In [24]: %timeit x2=x1.copy();x=frombuffer(x2,uint32) # time offset 1 loops, best of 3: 431 ms per loop In [25]: %timeit x2=x1.copy();x=frombuffer(x2,uint32);x.sort() 1 loops, best of 3: 37.8 s per loop In [26]: %timeit x2=x1.copy();x=frombuffer(x2,uint32);para(1) 1 loops, best of 3: 5.7 s per loop In [27]: %timeit x2=x1.copy();x=frombuffer(x2,uint32);sort() 1 loops, best of 3: 4.02 s per loop 

Así que una solución de python pura con una ganancia de 10x (37s-> 3.5s) en mi máquina pobre de 1GHz. Se puede mejorar con más núcleos y multifusión.