[Numpy-discussion] Radix sort?
Jon Wright
wright@esrf...
Tue Jun 19 13:08:36 CDT 2007
Dear numpy experts,
I see from the docs that there seem to be 3 sorting algorithms for array
data (quicksort, mergesort and heapsort). After hearing a rumour about
radix sorts and floats I google'd and now I'm wondering about a radix
sort for numpy (and Numeric) scalars? See:
http://www.stereopsis.com/radix.html
http://en.wikipedia.org/wiki/Radix_sort
The algorithm is apparently a trick where no comparisons are used. A
shockingly bad benchmark of a swig wrapped test function below suggests
it really is quicker than the array.sort() numpy method for uint8. At
least with >256 element uint8 test arrays (numpy1.0.3, mingw32, winXP),
for me, today; of course ymmv.
With large enough arrays of all of the integer numeric types and also
ieee reals, appropriate versions of the radix sort might be able to:
"... kick major ass in the speed department."
[http://www.cs.berkeley.edu/~jrs/61b/lec/36]
Have these sorts already been implemented in scipy? Can someone share
some expertise in this area? There is a problem about the sorting not
being done in-place (eg: better for argsort than sort). I see there is a
lexsort, but it is not 100% obvious to me how to use that to sort
scalars, especially ieee floats. If numpy is a library for handling big
chunks of numbers with knowable binary representations, I guess there
might be some general interest in having this radix sort compiled for
the platforms where it works?
Thanks for any opinions,
Jon
---
==setup.py==
from distutils.core import setup, Extension
e = Extension("_sort_radix_uint8",sources=["radix_wrap.c"])
setup(ext_modules=[e])
==radix.i==
%module sort_radix_uint8
%{
#define SWIG_FILE_WITH_INIT
void sort_radix_uint8(unsigned char * ar, int len){
int hits[256];
int i, j;
for(i=0;i<256;i++) hits[i]=0;
for(i=0;i<len;i++) hits[(int)ar[i]]+=1;
i=0; j=0;
while(i<256){
while(hits[i]>0){
/* shortcut for uint8 */
ar[j++] = (unsigned char) i;
hits[i]--;
}
i++;
}
}
%}
%init %{
import_array();
%}
%include "numpy.i"
%apply (unsigned char* INPLACE_ARRAY1, int DIM1) { (unsigned char * ar,
int len)}
void sort_radix_uint8(unsigned char * ar, int len);
== test.py ==
import numpy,sort_radix_uint8
from time import clock
def clearcache():
t = numpy.random.random(1024*1024)
t=t*t
def test(i):
t = numpy.random.random(i)
a1 = (t*255).astype(numpy.uint8)
a2 = (t*255).astype(numpy.uint8)
clearcache()
tick=clock()
a1.sort()
t1 = clock()-tick
clearcache()
tick = clock()
sort_radix_uint8.sort_radix_uint8(a2)
t2 = clock() - tick
assert a1.all() == a2.all()
clearcache()
tick=clock()
a1.sort() # already sorted...
t3 = clock()-tick
return t1,t2,t3
for j in range(25):
r = test(pow(2,j))
print numpy.argmin(r), j,pow(2,j),r
More information about the Numpy-discussion
mailing list