Подтвердить что ты не робот

Как отсортировать массив целых чисел быстрее, чем quicksort?

Сортировка массива целых чисел с помощью quicksort numpy стала узкое место моего алгоритма. К сожалению, у numpy нет еще не отсортировано. Хотя подсчет сортировки будет однострочным в numpy:

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

см. принятый ответ Как я могу векторизовать этот счетчик python, так что это абсолютно так же быстро, как это может быть? вопрос, целые числа в моем приложении может работать от 0 до 2**32.

Я застрял в быстрой сортировке?


<суб > Этот пост был в первую очередь мотивирован Группировка Numpy с использованием производительности itertools.groupby вопрос.
Также отметим, что не просто нормально спросить и ответить на ваш собственный вопрос, это явно поощряется. Суб >

4b9b3361

Ответ 1

Нет, вы не застряли в быстрой сортировке. Вы могли бы использовать, например, integer_sort от Boost.Sort или u4_sort из usort. При сортировке этого массива:

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

Получаю следующие результаты:

NumPy quicksort:         8.636 s  1.0  (baseline)
Boost.Sort integer_sort: 4.327 s  2.0x speedup
usort u4_sort:           2.065 s  4.2x speedup

Я бы не стал делать выводы на основе этого единственного эксперимента и использовать usort слепо. Я бы проверил свои фактические данные и измерил, что произойдет. Ваш пробег будет изменяться в зависимости от ваших данных и вашего устройства. integer_sort в Boost.Sort имеет богатый набор настроек для настройки, см. .

Ниже я описываю два способа вызова собственной C или С++ функции из Python. Несмотря на длинное описание, это довольно легко сделать.


Boost.Sort

Поместите эти строки в файл spreadsort.cpp:

#include <cinttypes>
#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);
    }
}

Он в основном создает шаблоны integer_sort для 32-битного целые числа без знака; часть extern "C" обеспечивает связь C путем отключения имя mangling. Предполагая, что вы используете gcc и что необходимые файлы включения boost находятся в каталоге /tmp/boost_1_60_0, вы можете скомпилировать его:

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

Ключевыми флагами являются -fPIC для генерации код, не зависящий от позиции и -shared для генерации общий объект .so файл. (Читайте документы gcc для получения дополнительной информации.)

Затем вы завершаете функцию spreadsort() С++ в Python, используя 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)

В качестве альтернативы вы можете использовать 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)

При вызовах cdll.LoadLibrary() и ffi.dlopen() я предположил, что путь к файлу spreadsort.so ./spreadsort.so. С другой стороны, вы можете написать

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

или

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

если вы добавите путь к spreadsort.so в среду LD_LIBRARY_PATH переменная. См. Также Общие библиотеки.

Использование. В обоих случаях вы просто вызываете вышеуказанную функцию-оболочку Python integer_sort() с вашим массивом numpy из 32-битных целых чисел без знака.


usort

Что касается u4_sort, вы можете скомпилировать его следующим образом:

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

Выполните эту команду в каталоге, где находится файл u4_sort.c. (Наверное, есть менее хакерский способ, но я не понял этого. просто просмотрел файл deps.mk в каталоге usort, чтобы узнать необходимые флаги компилятора и включить пути.)

Затем вы можете обернуть функцию C следующим образом:

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)

В приведенном выше коде я предположил, что путь к u4_sort.so был добавляется к переменной среды LD_LIBRARY_PATH.

Использование. Как и раньше, с помощью Boost.Sort вы просто вызываете вышеуказанную функцию-оболочку Python integer_sort() с помощью массива numpy из 32-битных целых чисел без знака.

Ответ 2

База 256 сортировки оснований (1 байт) может генерировать матрицу подсчетов, используемых для определения размера ведра в 1 проход данных, а затем занимает 4 прохода, чтобы выполнить сортировку. В моей системе Intel 2600K 3.4ghz, используя сборку Visual Studio для С++-программы, требуется около 1,15 секунды для сортировки беспорядочных 32-битных целых чисел в 10 ^ 8 psuedo.

Глядя на код u4_sort, упомянутый в ответе Али, программы похожи, но u4_sort использует размеры полей {10,11,11}, принимая 3 прохода для сортировки данных и 1 проход для копирования назад, в то время как в этом примере используется поле размеры {8,8,8,8}, взяв 4 прохода для сортировки данных. Вероятно, этот процесс ограничивает пропускную способность памяти из-за записи произвольного доступа, поэтому оптимизации в u4_sort, такие как макросы для сдвига, один цикл с фиксированным сдвигом на поле, не очень помогают. Мои результаты лучше, вероятно, из-за различий в системе и/или компиляторах. (Примечание: u8_sort для 64-битных целых чисел).

Пример кода:

//  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
    }
}

Ответ 3

Радиус-сортировка с python/numba(0.23), в соответствии с программой @rcgldr C, с многопоточным процессором на 2 ядра.

сначала сортировка radix в numba, с двумя глобальными массивами для эффективности.

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<<nbits
mask=buffsize-1

@jit(nopython=True,nogil=True)
def radix(x,y):
    xs=x.size
    dec=0
    while dec < 32 :
        u=np.zeros(buffsize,uint32)
        k=0
        while k<xs:
            u[(x[k]>>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<xs:
            j=(x[k]>>dec)&mask
            y[u[j]]=x[k]
            u[j]+=1
            k+=1
        x,y=y,x
        dec+=nbits

то параллелизация возможна с опцией 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<n:
        if y[kl]<x[kr] :
            x[k]=y[kl]
            kl+=1
            if kl==m : break
        else :
            x[k]=x[kr]
            kr+=1
        k+=1

def sort():
    para(2)
    y[:m]=x[:m]
    fuse(x,y)

:

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      

Итак, чистое решение python с выигрышем 10x (37s- > 3,5s) на моей бедной машине с частотой 1 ГГц. Может быть увеличено с помощью большего количества ядер и многоточечной обработки.