Faster groupby (#1096)

* Tweaked sort.max_radix_bits constant to speed up sorting
* Groupby operations are now significantly faster (especially function first) because columns are no longer materialized before performing the reduction.

Speed improvements:
* sort(): 221.671 ms -> 85.602 ms
* first(): 76.152 ms -> 0.084 ms
* sum(): 179.940 ms -> 69.260 ms

Author: st-pasha

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 - during grouping, group keys are now added automatically to the select list.
 - implement `sum()` reducer.
 - `==` operator now works for string columns too.
+- Improved performance of groupby operations.
 
 #### Fixed
 - fread will no longer emit an error if there is an NA string in the header.

c/datatable.cc

@@ -31,10 +31,9 @@ DataTable::DataTable(Column** cols)
   nrows = cols[0]->nrows;
 
   for (Column* col = cols[++ncols]; cols[ncols] != nullptr; ++ncols) {
-    // TODO: restore, once Column also uses RowIndex
-    // if (rowindex != col->rowindex()) {
-    //   throw ValueError() << "Mismatched RowIndex in Column " << ncols;
-    // }
+    if (rowindex != col->rowindex()) {
+      throw ValueError() << "Mismatched RowIndex in Column " << ncols;
+    }
     if (nrows != col->nrows) {
       throw ValueError() << "Mismatched length in Column " << ncols << ": "
                          << "found " << col->nrows << ", expected " << nrows;

c/expr/binaryop.cc

@@ -475,6 +475,8 @@ static mapperfn resolve0(SType lhs_type, SType rhs_type, int opcode, void** para
 
 Column* binaryop(int opcode, Column* lhs, Column* rhs)
 {
+  lhs->reify();
+  rhs->reify();
   int64_t lhs_nrows = lhs->nrows;
   int64_t rhs_nrows = rhs->nrows;
   SType lhs_type = lhs->stype();

c/expr/py_expr.cc

@@ -58,7 +58,7 @@ PyObject* expr_column(PyObject*, PyObject* args)
     PyErr_Format(PyExc_ValueError, "Invalid column index %lld", index);
   }
   Column* col = dt->columns[index]->shallowcopy(ri);
-  col->reify();
+  // col->reify();
   return pycolumn::from_column(col, nullptr, 0);
 }
 

c/expr/reduceop.cc

@@ -40,13 +40,14 @@ static Column* reduce_first(Column* arg, const Groupby& groupby) {
     return Column::new_data_column(arg->stype(), 0);
   }
   size_t ngrps = groupby.ngroups();
-  arr32_t indices(ngrps);
-  // TODO: avoid copy (by allowing RowIndex to be created from a MemoryRange)
-  std::memcpy(indices.data(), groupby.offsets_r(), ngrps * sizeof(int32_t));
-  RowIndex ri = RowIndex::from_array32(std::move(indices), true);
-  Column* res = arg->shallowcopy(ri);
-  res->reify();
-  return res;
+  // groupby.offsets array has length `ngrps + 1` and contains offsets of the
+  // beginning of each group. We will take this array and reinterpret it as a
+  // RowIndex (taking only the first `ngrps` elements). Applying this rowindex
+  // to the column will produce the vector of first elements in that column.
+  arr32_t indices(ngrps, groupby.offsets_r());
+  RowIndex ri = RowIndex::from_array32(std::move(indices), true)
+                .uplift(arg->rowindex());
+  return arg->shallowcopy(ri);
 }
 
 
@@ -65,11 +66,12 @@ static void sum_skipna(const int32_t* groups, int32_t grp, void** params) {
   OT sum = 0;
   int32_t row0 = groups[grp];
   int32_t row1 = groups[grp + 1];
-  for (int32_t i = row0; i < row1; ++i) {
-    IT x = inputs[i];
-    if (ISNA<IT>(x)) continue;
-    sum += static_cast<OT>(x);
-  }
+  col0->rowindex().strided_loop(row0, row1, 1,
+    [&](int64_t i) {
+      IT x = inputs[i];
+      if (!ISNA<IT>(x))
+        sum += static_cast<OT>(x);
+    });
   outputs[grp] = sum;
 }
 
@@ -90,15 +92,16 @@ static void mean_skipna(const int32_t* groups, int32_t grp, void** params) {
   OT delta = 0;
   int32_t row0 = groups[grp];
   int32_t row1 = groups[grp + 1];
-  for (int32_t i = row0; i < row1; ++i) {
-    IT x = inputs[i];
-    if (ISNA<IT>(x)) continue;
-    OT y = static_cast<OT>(x) - delta;
-    OT t = sum + y;
-    delta = (t - sum) - y;
-    sum = t;
-    cnt++;
-  }
+  col0->rowindex().strided_loop(row0, row1, 1,
+    [&](int64_t i) {
+      IT x = inputs[i];
+      if (ISNA<IT>(x)) return;
+      OT y = static_cast<OT>(x) - delta;
+      OT t = sum + y;
+      delta = (t - sum) - y;
+      sum = t;
+      cnt++;
+    });
   outputs[grp] = cnt == 0? GETNA<OT>() : sum / cnt;
 }
 
@@ -120,15 +123,16 @@ static void stdev_skipna(const int32_t* groups, int32_t grp, void** params) {
   int64_t cnt = 0;
   int32_t row0 = groups[grp];
   int32_t row1 = groups[grp + 1];
-  for (int32_t i = row0; i < row1; ++i) {
-    IT x = inputs[i];
-    if (ISNA<IT>(x)) continue;
-    cnt++;
-    OT t1 = x - mean;
-    mean += t1 / cnt;
-    OT t2 = x - mean;
-    m2 += t1 * t2;
-  }
+  col0->rowindex().strided_loop(row0, row1, 1,
+    [&](int64_t i) {
+      IT x = inputs[i];
+      if (ISNA<IT>(x)) return;
+      cnt++;
+      OT t1 = x - mean;
+      mean += t1 / cnt;
+      OT t2 = x - mean;
+      m2 += t1 * t2;
+    });
   outputs[grp] = cnt <= 1? GETNA<OT>() : std::sqrt(m2 / (cnt - 1));
 }
 
@@ -147,12 +151,13 @@ static void min_skipna(const int32_t* groups, int32_t grp, void** params) {
   T res = infinity<T>();
   int32_t row0 = groups[grp];
   int32_t row1 = groups[grp + 1];
-  for (int32_t i = row0; i < row1; ++i) {
-    T x = inputs[i];
-    if (!ISNA<T>(x) && x < res) {
-      res = x;
-    }
-  }
+  col0->rowindex().strided_loop(row0, row1, 1,
+    [&](int64_t i) {
+      T x = inputs[i];
+      if (!ISNA<T>(x) && x < res) {
+        res = x;
+      }
+    });
   outputs[grp] = res;
 }
 
@@ -171,12 +176,13 @@ static void max_skipna(const int32_t* groups, int32_t grp, void** params) {
   T res = -infinity<T>();
   int32_t row0 = groups[grp];
   int32_t row1 = groups[grp + 1];
-  for (int32_t i = row0; i < row1; ++i) {
-    T x = inputs[i];
-    if (!ISNA<T>(x) && x > res) {
-      res = x;
-    }
-  }
+  col0->rowindex().strided_loop(row0, row1, 1,
+    [&](int64_t i) {
+      T x = inputs[i];
+      if (!ISNA<T>(x) && x > res) {
+        res = x;
+      }
+    });
   outputs[grp] = res;
 }
 

c/expr/unaryop.cc

@@ -135,6 +135,7 @@ static mapperfn resolve0(SType stype, int opcode) {
 Column* unaryop(int opcode, Column* arg)
 {
   if (opcode == OpCode::Plus) return arg->shallowcopy();
+  arg->reify();
 
   SType arg_type = arg->stype();
   SType res_type = arg_type;

c/py_columnset.cc

@@ -140,6 +140,7 @@ PyObject* columns_from_columns(PyObject*, PyObject* args)
       return nullptr;
     }
     reinterpret_cast<pycolumn::obj*>(elem)->ref = nullptr;
+    columns[i]->reify();
   }
   columns[ncols] = nullptr;
 

c/utils/array.h

@@ -48,10 +48,14 @@ template <typename T> class array
   private:
     T* x;
     size_t n;
+    bool owned;
+    int64_t : 56;
 
   public:
-    array(size_t len = 0) : x(nullptr), n(0) { resize(len); }
-    ~array() { std::free(x); }
+    array(size_t len = 0) : x(nullptr), n(0), owned(true) { resize(len); }
+    array(size_t len, const T* ptr)
+      : x(const_cast<T*>(ptr)), n(len), owned(false) {}
+    ~array() { if (owned) std::free(x); }
     // copy-constructor and assignment are forbidden
     array(const array<T>&) = delete;
     array<T>& operator=(const array<T>&) = delete;
@@ -64,12 +68,14 @@ template <typename T> class array
       using std::swap;
       swap(first.x, second.x);
       swap(first.n, second.n);
+      swap(first.owned, second.owned);
     }
 
     template <typename S> array<S> cast() {
       array<S> res;
       res.n = n * sizeof(T) / sizeof(S);
       res.x = reinterpret_cast<S*>(x);
+      res.owned = owned;
       x = nullptr;
       return res;
     }
@@ -79,7 +85,7 @@ template <typename T> class array
       size_t size = sizeof(T) * n;
       x = nullptr;
       n = 0;
-      return MemoryRange(size, ptr, /* own = */ true);
+      return MemoryRange(size, ptr, /* own = */ owned);
     }
 
     // Standard operators
@@ -95,6 +101,9 @@ template <typename T> class array
 
     void resize(size_t newn) {
       if (newn == n) return;
+      if (!owned) {
+        throw MemoryError() << "Cannot resize array: not owned";
+      }
       if (newn == 0) {
         std::free(x);
         x = nullptr;

datatable/frame.py

@@ -981,7 +981,7 @@ def data_viewer(row0, row1, col0, col1):
     "sort.max_chunk_length", xtype=int, default=1024, core=True)
 
 options.register_option(
-    "sort.max_radix_bits", xtype=int, default=8, core=True)
+    "sort.max_radix_bits", xtype=int, default=12, core=True)
 
 options.register_option(
     "sort.over_radix_bits", xtype=int, default=8, core=True)

tests/fread/test_fread_api.py

@@ -592,7 +592,7 @@ def test_fread_skip_blank_lines_true():
     assert d0.topython() == [[1, 3], [2, 4]]
 
 
-@pytest.mark.xfail()
+@pytest.mark.skip("Issue #838")
 def test_fread_skip_blank_lines_false():
     inp = "A,B\n1,2\n  \n\n3,4\n"
     with pytest.warns(DatatableWarning) as ws: