octave-maintainers
[Top][All Lists]

## Re: benchmarks - sort

 From: David Bateman Subject: Re: benchmarks - sort Date: Wed, 21 Jan 2004 23:09:52 +0100 User-agent: Mutt/1.4.1i

```Hi Thomas,

Well the fact is I've come a bit further in testing some sorting
algorithms.  I've tried several different cases. Firstly, what Alois
looking at radix sorting. Then Paul Kienzle strongly suggested the
Python mergesort algorithm as having great abilities especially for
partially ordered lists. I've taken the Python code and implemented it
as a class and use it to sort either directly on the doubles or
recasting as uint64 and use the same trick I used for the radix
sort.

One problem I have with the Python code, is that when calling
"[b,bi] = sort(a)", I have to create a new class containing the element and
the index like

template <class T>
class
vec_index
{
public:
vec_index (void) : vec (0), indx (-1) { }
vec_index (T v, int i) : vec (v), indx (i) { }
vec_index<T> &operator = (const vec_index<T> &x);
T vec;
int indx;
};

template vec_index<double>;

bool
operator < (const vec_index<double> &a, const vec_index<double> &b)
{
return (xisnan(b.vec) || (a.vec < b.vec));
}

This is due to the fact that the calls in the class to sort the
elements are of the form "tmp = *a; *a = *b; *b = tmp;". This imposes
quite a penalty on the mergesort. I can really see another way to do
this without writing specific versions of the code for all
instantiations of the mergesort class for with and without indexing
(Ugly)... Any suggestions are welcome.

In any case there are 5 test case, with and without indexing

1a) Matlab R12 "b = sort(a)" - baseline for unindexed octave tests
1b) Matlab R12 "[b,bi] = sort(a)" - baseline for indexed octave tests
2a) Octave "b = sort(a)" - orginal octave code
2b) Octave "[b,bi] = sort(a)" - orginal octave code with indexing
3b) Octave "[b,bi] = radix(a)" - Radix sort test code with indexing
4a) Octave "b = merge_double(a)" - Mergesort on doubles
4b) Octave "[b,bi] = merge_double(a)" - Mergesort on doubles with indexing
5a) Octave "b = merge_ieee7543(a)" - Mergesort on IEEE754 case as uint64
5b) Octave "[b,bi] = merge_double(a)" - Mergesort on IEEE754 case as uint64
with indexing

I've benchmarked the matlab code seperately (in "matlab.log") and used
this to benchmark the other four versions of the code, with and without
returning the index array. All tests were run on the same machine under
the same conditions. I also trying to do a bit of averaging in the tests
to remove variation in the runtimes. I attach the code I used in the tarball
in this mail.

The names of the tests I used are the same as those in the Python code and
are

*sort     Randomly distributed lists
\sort     Descending list
/sort     Ascending list
3sort     Ascending list with 3 values randomly interchanged
+sort     Ascending list with the last 10 values replace with random data
=sort     All values equal

There are two tables in "octave.log" for each test. The first table is the
runtime of each iteration, while the second table is the relative time wrt
Matlab. So as to not bias the results against Octave, with its problems
with "for" loops, I've assumed that

a = randn(n,rep)
t = cputime;
b = sort(a);
time (cputime - t) / rep;

gives the runtime of "b = sort(a)", where "a = randn(n,1)".

The basic result is that the radix sort compares well against Matlab
for randomly distributed values (1.10 times the speed of Matlab), but
as this sort takes no advantage of partial ordering, it breaks down in
these cases. The mergesort with double values is faster than Matlab
for partailly ordered lists, but about 2.5 to 3 times slower than
Matlab for random lists. Interestingly, although merge_double is
significantly slower with indexing, it maintains the same relationship
with Matlab.

The clear winner for me is the mergesort with integer comparison of
IEEE754 values. This is roughly two times faster than Matlab for
partially ordered lists, while being only slightly slower (1.29
relative to Matlab) than the radix sort (1.10 relative to Matlab).
This algorithm is roughly 6 times faster than the current Octave
sort code. The performance degrades slightly with indexing (roughly
2.0 relative to Matlab).

Note that my test code, checked the correctness of the sorting
algorithms, that it sorted NaN, Inf and -Inf correctly and that
it was stable.

There are several questions to answer before this might be rewritten
for inclusion either in octave or octave-forge.

* Does octave support any platform with 8-byte integers? If so the
mergesort class could be used with doubles in this case. If not I
can simplify the code. In any case having mergesort as a class
means that other parts of octave might use this class (sort char
matrices for instance).

* Is there any other way to win back the speed loss I get for sorting
with indexing with the mergesort class? Maybe someone could look at
merge.cc and oct-sort.cc and give their thoughts?

* Is there any chance of this going into octave, or should I be targeting
octave-forge?

* Should I bother implementing complex sorting (done it for radix sort
but not in the attached tar-ball)?

Cheers
David

Daprès THOMAS Paul Richard <address@hidden> (le 21/01/2004):
> For what it is worth ( 0.02 euros?), in the course of educating myself about
> the C++ standard library, I did a test of the sorting capability of
> multimaps.  It turns out to be a bit off the performance discussion, in that
> the timing is EXACTLY the same as that of octave's sort and about a factor
> of five slower than Matlab R13 on the same machine (not very surprising,
> since the standard library uses exactly the same algorithm as octave!).
> However, it is interesting just how concise a routine, exploiting the
> standard library can become.  Please find the routine below.
>
> Paul Thomas
>
>  //Test of auto-sorting properties of standard library multi-map
> //In input vector x is inserted element by element into a
> //multi-map vmap, together with the input index of the element.
> //Since multi-map insert sorts according to the value of the first
> //element, reading the map pairs back, using the multi-map iterator,
> //yields the sorted values and the sorted index.
> //
> //Paul Thomas                                   17/01/04
>
> #include <octave/oct.h>
> #include <octave/variables.h>
> #include <map>
> using namespace std;
>
> DEFUN_DLD (mysort, args, ,
>    "mysort. Call using \
>    [a,b]=mysort(x)")
> {
>    typedef multimap<double,double> ddmap;       //<value,index>
>    ColumnVector vin(args(0).vector_value());    //turn input into
> ColumnVector
>    int vinlen=vin.length();                     //number of input elements
>    ddmap vmap;                                  //multimap container for x
> and its index
>    ddmap::iterator pos;                         //iterator for the multimap
>    double didx=1;                               //double index
>    for (int idx=0;idx<vinlen;idx++)             //let multimap insert sort
> by value
>    {
>       vmap.insert(make_pair(vin(idx),didx++));  //insert value,index
>    }
>    int idx=0;                                   //output counter
>    ColumnVector idxout(vinlen);                 //sorted index
>    for (pos=vmap.begin();pos!=vmap.end();pos++,idx++)
>    {
>       vin(idx)=pos->first;                      //sorted value
>       idxout(idx)=pos->second;                  //sorted index
>    }
>    octave_value_list retval(vin);               //sorted value to output
>    retval.append(octave_value(idxout));         //sorted index to output
>    return retval;                               //return
> }
>
>
> ---
> Outgoing mail is certified Virus Free.
> Checked by AVG anti-virus system (http://www.grisoft.com).
> Version: 6.0.564 / Virus Database: 356 - Release Date: 19/01/04
>

--
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: