[Numpy-discussion] Radix sort?

[Jon Wright]

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