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
sort for numpy (and Numeric) scalars? See:

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
setup(ext_modules=[e])

%{
#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 ==
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()
t2 = clock() - tick

assert a1.all() == a2.all()

clearcache()
tick=clock()
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

```