# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1265624212 -3600
# Node ID 3a8c13b71612471c5db8f03e98f9ae98f6094c84
# Parent  272179888089ab8ce3c36139c6d474c67aae1e02
implement special-case optimization for sort of index vectors

diff -r 272179888089 -r 3a8c13b71612 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Mon Feb 08 01:36:18 2010 -0500
+++ b/liboctave/ChangeLog	Mon Feb 08 11:16:52 2010 +0100
@@ -1,3 +1,18 @@
+2010-02-08  Jaroslav Hajek  <highegg@gmail.com>
+
+	* idx-vector.h (idx_vector::idx_base_rep::sort_idx): New pure virtual
+	function.
+	(idx_vector::idx_colon_rep::sort_idx,
+	idx_vector::idx_range_rep::sort_idx,
+	idx_vector::idx_scalar_rep::sort_idx,
+	idx_vector::idx_vector_rep::sort_idx,
+	idx_vector::idx_mask_rep::sort_idx): New override decls.
+	idx_vector::sort (Array<octave_idx_type>&): New method.
+	* idx-vector.cc (idx_vector::idx_range_rep::sort_idx,
+	idx_vector::idx_vector_rep::sort_idx): New methods.
+	(idx_vector::idx_vector_rep::sort_uniq_clone): Rewrite
+	to use bucket sort under plausible circumstances.
+
 2010-02-08  Jaroslav Hajek  <highegg@gmail.com>
 
 	* Array.cc (Array<T>::permute): Fix result dimensions when inv=true.
diff -r 272179888089 -r 3a8c13b71612 liboctave/idx-vector.cc
--- a/liboctave/idx-vector.cc	Mon Feb 08 01:36:18 2010 -0500
+++ b/liboctave/idx-vector.cc	Mon Feb 08 11:16:52 2010 +0100
@@ -28,6 +28,7 @@
 #endif
 
 #include <cstdlib>
+#include <memory>
 
 #include <iostream>
 
@@ -85,6 +86,16 @@
     return i;
 }
 
+idx_vector::idx_base_rep *
+idx_vector::idx_colon_rep::sort_idx (Array<octave_idx_type>&)
+{
+  (*current_liboctave_error_handler)
+    ("internal error: idx_colon_rep::sort_idx");
+
+  count++;
+  return this;
+}
+
 std::ostream& 
 idx_vector::idx_colon_rep::print (std::ostream& os) const
 {
@@ -162,6 +173,26 @@
     }
 }
 
+idx_vector::idx_base_rep *
+idx_vector::idx_range_rep::sort_idx (Array<octave_idx_type>& idx)
+{
+  if (step < 0 && len > 0)
+    {
+      idx.clear (1, len);
+      for (octave_idx_type i = 0; i < len; i++)
+        idx.xelem (i) = len - 1 - i;
+      return new idx_range_rep (start + (len - 1)*step, len, -step, DIRECT);
+    }
+  else
+    {
+      idx.clear (1, len);
+      for (octave_idx_type i = 0; i < len; i++)
+        idx.xelem (i) = i;
+      count++;
+      return this;
+    }
+}
+
 std::ostream& 
 idx_vector::idx_range_rep::print (std::ostream& os) const
 {
@@ -238,6 +269,15 @@
   return data;
 }
 
+idx_vector::idx_base_rep *
+idx_vector::idx_scalar_rep::sort_idx (Array<octave_idx_type>& idx)
+{
+  idx.clear (1, 1);
+  idx.fill (0);
+  count++;
+  return this;
+}
+
 std::ostream& idx_vector::idx_scalar_rep::print (std::ostream& os) const
 {
   return os << data;
@@ -385,16 +425,134 @@
 idx_vector::idx_base_rep *
 idx_vector::idx_vector_rep::sort_uniq_clone (bool uniq)
 {
-  octave_idx_type *new_data = new octave_idx_type[len];
-  std::copy (data, data + len, new_data);
-  std::sort (new_data, new_data + len);
-  octave_idx_type new_len;
-  if (uniq)
-    new_len = std::unique (new_data, new_data + len) - new_data;
-  else 
-    new_len = len;
+  if (len == 0)
+    {
+      count++;
+      return this;
+    }
+
+  // This is wrapped in auto_ptr so that we don't leak on out-of-memory.
+  std::auto_ptr<idx_vector_rep> new_rep ( 
+    new idx_vector_rep (0, len, ext, orig_dims, DIRECT));
+
+  if (ext > len*xlog2 (1.0 + len))
+    {
+      // Use standard sort via octave_sort.
+      octave_idx_type *new_data = new octave_idx_type[len];
+      new_rep->data = new_data;
+
+      std::copy (data, data + len, new_data);
+      octave_sort<octave_idx_type> lsort;
+      lsort.set_compare (ASCENDING);
+      lsort.sort (new_data, len);
+
+      if (uniq)
+        {
+          octave_idx_type new_len = std::unique (new_data, new_data + len) - new_data;
+          new_rep->len = new_len;
+          if (new_rep->orig_dims.length () == 2 && new_rep->orig_dims(0) == 1)
+            new_rep->orig_dims = dim_vector (1, new_len);
+          else
+            new_rep->orig_dims = dim_vector (new_len, 1);
+        }
+    }
+  else if (uniq)
+    {
+      // Use two-pass bucket sort (only a mask array needed).
+      OCTAVE_LOCAL_BUFFER_INIT (bool, has, ext, false);
+      for (octave_idx_type i = 0; i < len; i++)
+        has[data[i]] = true;
+
+      octave_idx_type new_len = 0;
+      for (octave_idx_type i = 0; i < ext; i++)
+        new_len += has[i];
+
+      new_rep->len = new_len;
+      if (new_rep->orig_dims.length () == 2 && new_rep->orig_dims(0) == 1)
+        new_rep->orig_dims = dim_vector (1, new_len);
+      else
+        new_rep->orig_dims = dim_vector (new_len, 1);
+
+      octave_idx_type *new_data = new octave_idx_type[new_len];
+      new_rep->data = new_data;
+
+      for (octave_idx_type i = 0, j = 0; i < ext; i++)
+        if (has[i])
+          new_data[j++] = i;
+    }
+  else
+    {
+      // Use two-pass bucket sort.
+      OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, cnt, ext, 0);
+      for (octave_idx_type i = 0; i < len; i++)
+        cnt[data[i]]++;
 
-  return new idx_vector_rep (new_data, new_len, ext, orig_dims, DIRECT);
+      octave_idx_type *new_data = new octave_idx_type[len];
+      new_rep->data = new_data;
+
+      for (octave_idx_type i = 0, j = 0; i < ext; i++)
+        {
+          for (octave_idx_type k = 0; k < cnt[i]; k++)
+            new_data[j++] = i;
+        }
+    }
+
+  return new_rep.release ();
+}
+
+idx_vector::idx_base_rep *
+idx_vector::idx_vector_rep::sort_idx (Array<octave_idx_type>& idx)
+{
+  // This is wrapped in auto_ptr so that we don't leak on out-of-memory.
+  std::auto_ptr<idx_vector_rep> new_rep (
+    new idx_vector_rep (0, len, ext, orig_dims, DIRECT));
+
+  if (ext > len*xlog2 (1.0 + len))
+    {
+      // Use standard sort via octave_sort.
+      idx.clear (orig_dims);
+      octave_idx_type *idx_data = idx.fortran_vec ();
+      for (octave_idx_type i = 0; i < len; i++)
+        idx_data[i] = i;
+
+      octave_idx_type *new_data = new octave_idx_type [len];
+      new_rep->data = new_data;
+      std::copy (data, data + len, new_data);
+
+      octave_sort<octave_idx_type> lsort;
+      lsort.set_compare (ASCENDING);
+      lsort.sort (new_data, idx_data, len);
+    }
+  else
+    {
+      // Use two-pass bucket sort.
+      OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, cnt, ext, 0);
+
+      for (octave_idx_type i = 0; i < len; i++)
+        cnt[data[i]]++;
+
+      idx.clear (orig_dims);
+      octave_idx_type *idx_data = idx.fortran_vec ();
+
+      octave_idx_type *new_data = new octave_idx_type [len];
+      new_rep->data = new_data;
+
+      for (octave_idx_type i = 0, k = 0; i < ext; i++)
+        {
+          octave_idx_type j = cnt[i];
+          cnt[i] = k;
+          k += j;
+        }
+
+      for (octave_idx_type i = 0; i < len; i++)
+        {
+          octave_idx_type j = data[i], k = cnt[j]++;
+          new_data[k] = j;
+          idx_data[k] = i;
+        }
+    }
+
+  return new_rep.release ();
 }
 
 std::ostream& 
@@ -518,6 +676,17 @@
     }
 }
 
+idx_vector::idx_base_rep *
+idx_vector::idx_mask_rep::sort_idx (Array<octave_idx_type>& idx)
+{
+  idx.clear (len, 1);
+  for (octave_idx_type i = 0; i < len; i++)
+    idx.xelem (i) = i;
+
+  count++;
+  return this;
+}
+
 const idx_vector idx_vector::colon (new idx_vector::idx_colon_rep ());
 
 idx_vector::idx_vector (const Array<bool>& bnda)
diff -r 272179888089 -r 3a8c13b71612 liboctave/idx-vector.h
--- a/liboctave/idx-vector.h	Mon Feb 08 01:36:18 2010 -0500
+++ b/liboctave/idx-vector.h	Mon Feb 08 11:16:52 2010 +0100
@@ -89,6 +89,8 @@
 
     // Sorts, maybe uniqifies, and returns a clone object pointer.
     virtual idx_base_rep *sort_uniq_clone (bool uniq = false) = 0;
+    // Sorts, and returns a sorting permutation (aka Array::sort).
+    virtual idx_base_rep *sort_idx (Array<octave_idx_type>&) = 0;
 
     // Checks whether the index is colon or a range equivalent to colon.
     virtual bool is_colon_equiv (octave_idx_type) const
@@ -135,6 +137,8 @@
     idx_base_rep *sort_uniq_clone (bool = false) 
       { count++; return this; }
 
+    idx_base_rep *sort_idx (Array<octave_idx_type>&);
+
     bool is_colon_equiv (octave_idx_type) const
       { return true; }
 
@@ -183,6 +187,8 @@
 
     idx_base_rep *sort_uniq_clone (bool uniq = false);
 
+    idx_base_rep *sort_idx (Array<octave_idx_type>&);
+
     bool is_colon_equiv (octave_idx_type n) const
       { return start == 0 && step == 1 && len == n; }
 
@@ -240,6 +246,8 @@
     idx_base_rep *sort_uniq_clone (bool = false)
       { count++; return this; }
 
+    idx_base_rep *sort_idx (Array<octave_idx_type>&);
+
     bool is_colon_equiv (octave_idx_type n) const
       { return n == 1 && data == 0; }
 
@@ -305,6 +313,8 @@
 
     idx_base_rep *sort_uniq_clone (bool uniq = false);
 
+    idx_base_rep *sort_idx (Array<octave_idx_type>&);
+
     dim_vector orig_dimensions (void) const
       { return orig_dims; }
 
@@ -369,6 +379,8 @@
     idx_base_rep *sort_uniq_clone (bool = false) 
       { count++; return this; }
 
+    idx_base_rep *sort_idx (Array<octave_idx_type>&);
+
     dim_vector orig_dimensions (void) const
       { return orig_dims; }
 
@@ -560,6 +572,9 @@
   idx_vector sorted (bool uniq = false) const
     { return idx_vector (rep->sort_uniq_clone (uniq)); }
 
+  idx_vector sorted (Array<octave_idx_type>& sidx) const
+    { return idx_vector (rep->sort_idx (sidx)); }
+
   dim_vector orig_dimensions (void) const { return rep->orig_dimensions (); }
 
   octave_idx_type orig_rows (void) const
diff -r 272179888089 -r 3a8c13b71612 src/ChangeLog
--- a/src/ChangeLog	Mon Feb 08 01:36:18 2010 -0500
+++ b/src/ChangeLog	Mon Feb 08 11:16:52 2010 +0100
@@ -1,3 +1,8 @@
+2010-02-08  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-re-mat.cc (octave_matrix::sort): Special-case sorting a known
+	index vector.
+
 2010-02-08  John W. Eaton  <jwe@octave.org>
 
 	* ov-class.cc (Fclass): If more than 1 argument, check that
diff -r 272179888089 -r 3a8c13b71612 src/ov-re-mat.cc
--- a/src/ov-re-mat.cc	Mon Feb 08 01:36:18 2010 -0500
+++ b/src/ov-re-mat.cc	Mon Feb 08 11:16:52 2010 +0100
@@ -271,6 +271,34 @@
   return retval;
 }
 
+octave_value 
+octave_matrix::sort (octave_idx_type dim, sortmode mode) const
+{
+  // If this matrix is known to be a valid index, and we're doing a vector sort,
+  // sort via idx_vector which may be more efficient occassionally.
+  if (idx_cache && mode == ASCENDING && ndims () == 2
+      && ((dim == 0 && matrix.columns () == 1) || (dim == 1 && matrix.rows () == 1)))
+    {
+      return index_vector ().sorted ();
+    }
+  else
+    return octave_base_matrix<NDArray>::sort (dim, mode);
+}
+
+octave_value 
+octave_matrix::sort (Array<octave_idx_type> &sidx, octave_idx_type dim,
+                     sortmode mode) const
+{
+  // If this matrix is known to be a valid index, and we're doing a vector sort,
+  // sort via idx_vector which may be more efficient occassionally.
+  if (idx_cache && mode == ASCENDING && ndims () == 2
+      && ((dim == 0 && matrix.columns () == 1) || (dim == 1 && matrix.rows () == 1)))
+    {
+      return index_vector ().sorted (sidx);
+    }
+  else
+    return octave_base_matrix<NDArray>::sort (sidx, dim, mode);
+}
 octave_value
 octave_matrix::convert_to_str_internal (bool, bool, char type) const
 {
diff -r 272179888089 -r 3a8c13b71612 src/ov-re-mat.h
--- a/src/ov-re-mat.h	Mon Feb 08 01:36:18 2010 -0500
+++ b/src/ov-re-mat.h	Mon Feb 08 11:16:52 2010 +0100
@@ -179,6 +179,10 @@
 
   octave_value diag (octave_idx_type k = 0) const;
 
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = ASCENDING) const;
+
   // Use matrix_ref here to clear index cache.
   void increment (void) { matrix_ref () += 1.0; }