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 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
