Columns support for matrix class

NEWS.md

@@ -1,5 +1,6 @@
 # cpp11 (development version)
 
+* Allow `cpp11::matrix` to be accessed either row-wise or column-wise (the default) depending on the user's choice (@alyst, #229)
 * Fixed `+` and `+=` operators of `r_vector::[const_]iterator` to conform the *iterators* concept:
   `+=` updates the iterator, and `+` returns the updated copy, while keeping the original unchanged (@alyst, #231)
 * Remove undefined behavior when constructing global `cpp11::sexp`s (#224)

cpp11test/R/cpp11.R

@@ -92,6 +92,10 @@ row_sums <- function(x) {
   .Call(`_cpp11test_row_sums`, x)
 }
 
+col_sums <- function(x) {
+  .Call(`_cpp11test_col_sums`, x)
+}
+
 protect_one_ <- function(x, n) {
   invisible(.Call(`_cpp11test_protect_one_`, x, n))
 }

cpp11test/src/cpp11.cpp

@@ -149,7 +149,7 @@ extern "C" SEXP _cpp11test_gibbs_cpp(SEXP N, SEXP thin) {
   END_CPP11
 }
 // matrix.cpp
-cpp11::doubles_matrix gibbs_cpp2(int N, int thin);
+cpp11::doubles_matrix<> gibbs_cpp2(int N, int thin);
 extern "C" SEXP _cpp11test_gibbs_cpp2(SEXP N, SEXP thin) {
   BEGIN_CPP11
     return cpp11::as_sexp(gibbs_cpp2(cpp11::as_cpp<cpp11::decay_t<int>>(N), cpp11::as_cpp<cpp11::decay_t<int>>(thin)));
@@ -170,10 +170,17 @@ extern "C" SEXP _cpp11test_gibbs_rcpp2(SEXP N, SEXP thin) {
   END_CPP11
 }
 // matrix.cpp
-cpp11::doubles row_sums(cpp11::doubles_matrix x);
+cpp11::doubles row_sums(cpp11::doubles_matrix<cpp11::by_row> x);
 extern "C" SEXP _cpp11test_row_sums(SEXP x) {
   BEGIN_CPP11
-    return cpp11::as_sexp(row_sums(cpp11::as_cpp<cpp11::decay_t<cpp11::doubles_matrix>>(x)));
+    return cpp11::as_sexp(row_sums(cpp11::as_cpp<cpp11::decay_t<cpp11::doubles_matrix<cpp11::by_row>>>(x)));
+  END_CPP11
+}
+// matrix.cpp
+cpp11::doubles col_sums(cpp11::doubles_matrix<cpp11::by_column> x);
+extern "C" SEXP _cpp11test_col_sums(SEXP x) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(col_sums(cpp11::as_cpp<cpp11::decay_t<cpp11::doubles_matrix<cpp11::by_column>>>(x)));
   END_CPP11
 }
 // protect.cpp
@@ -389,6 +396,7 @@ extern "C" {
 extern SEXP run_testthat_tests(SEXP);
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_cpp11test_col_sums",                 (DL_FUNC) &_cpp11test_col_sums,                 1},
     {"_cpp11test_cpp11_add_vec_for_",       (DL_FUNC) &_cpp11test_cpp11_add_vec_for_,       2},
     {"_cpp11test_cpp11_insert_",            (DL_FUNC) &_cpp11test_cpp11_insert_,            1},
     {"_cpp11test_cpp11_release_",           (DL_FUNC) &_cpp11test_cpp11_release_,           1},

cpp11test/src/matrix.cpp

@@ -4,7 +4,7 @@
 using namespace cpp11;
 
 [[cpp11::register]] SEXP gibbs_cpp(int N, int thin) {
-  cpp11::writable::doubles_matrix mat(N, 2);
+  cpp11::writable::doubles_matrix<> mat(N, 2);
   double x = 0, y = 0;
   for (int i = 0; i < N; i++) {
     for (int j = 0; j < thin; j++) {
@@ -18,8 +18,8 @@ using namespace cpp11;
   return mat;
 }
 
-[[cpp11::register]] cpp11::doubles_matrix gibbs_cpp2(int N, int thin) {
-  cpp11::writable::doubles_matrix mat(N, 2);
+[[cpp11::register]] cpp11::doubles_matrix<> gibbs_cpp2(int N, int thin) {
+  cpp11::writable::doubles_matrix<> mat(N, 2);
   double x = 0, y = 0;
   for (int i = 0; i < N; i++) {
     for (int j = 0; j < thin; j++) {
@@ -67,11 +67,11 @@ using namespace Rcpp;
   return (mat);
 }
 
-[[cpp11::register]] cpp11::doubles row_sums(cpp11::doubles_matrix x) {
+[[cpp11::register]] cpp11::doubles row_sums(cpp11::doubles_matrix<cpp11::by_row> x) {
   cpp11::writable::doubles sums(x.nrow());
 
   int i = 0;
-  for (auto& row : x) {
+  for (auto row : x) {
     sums[i] = 0.;
     for (auto&& val : row) {
       if (cpp11::is_na(val)) {
@@ -85,3 +85,22 @@ using namespace Rcpp;
 
   return sums;
 }
+
+[[cpp11::register]] cpp11::doubles col_sums(cpp11::doubles_matrix<cpp11::by_column> x) {
+  cpp11::writable::doubles sums(x.ncol());
+
+  int i = 0;
+  for (auto col : x) {
+    sums[i] = 0.;
+    for (auto&& val : col) {
+      if (cpp11::is_na(val)) {
+        sums[i] = NA_REAL;
+        break;
+      }
+      sums[i] += val;
+    }
+    ++i;
+  }
+
+  return sums;
+}

cpp11test/src/test-matrix.cpp

@@ -7,7 +7,7 @@
 
 context("matrix-C++") {
   test_that("matrix dim attributes are correct for writable matrices") {
-    cpp11::writable::doubles_matrix x(5, 2);
+    cpp11::writable::doubles_matrix<cpp11::by_row> x(5, 2);
 
     cpp11::integers dim(SEXP(x.attr("dim")));
 
@@ -16,42 +16,107 @@ context("matrix-C++") {
 
     expect_true(x.nrow() == 5);
     expect_true(x.ncol() == 2);
+    expect_true(x.nslices() == 5);
+    expect_true(x.slice_size() == 2);
+    expect_true(x.slice_stride() == 5);
+    expect_true(x.slice_offset(0) == 0);
+    expect_true(x.slice_offset(1) == 1);
+    expect_true(x[1].size() == 2);
+    expect_true(x[1].stride() == 5);
   }
   test_that("matrix dim attributes are correct for read only matrices") {
     auto getExportedValue = cpp11::package("base")["getExportedValue"];
 
-    cpp11::doubles_matrix x(SEXP(getExportedValue("datasets", "volcano")));
+    test_that("matrix<by_row> attributes are correct") {
+      cpp11::doubles_matrix<cpp11::by_row> x(getExportedValue("datasets", "volcano"));
 
-    expect_true(x.size() == 5307);
-    expect_true(x.nrow() == 87);
-    expect_true(x.ncol() == 61);
+      expect_true(x.size() == 5307);
+      expect_true(x.nrow() == 87);
+      expect_true(x.ncol() == 61);
+      expect_true(x.nslices() == 87);
+      expect_true(x.slice_size() == 61);
+      expect_true(x.slice_stride() == 87);
+      expect_true(x.slice_offset(0) == 0);
+      expect_true(x.slice_offset(1) == 1);
+      expect_true(x[1].size() == 61);
+      expect_true(x[1].stride() == 87);
+    }
+    test_that("matrix<by_column> attributes are correct") {
+      cpp11::doubles_matrix<cpp11::by_column> x(getExportedValue("datasets", "volcano"));
+
+      expect_true(x.size() == 5307);
+      expect_true(x.nrow() == 87);
+      expect_true(x.ncol() == 61);
+      expect_true(x.nslices() == 61);
+      expect_true(x.slice_size() == 87);
+      expect_true(x.slice_stride() == 1);
+      expect_true(x.slice_offset(0) == 0);
+      expect_true(x.slice_offset(1) == 87);
+      expect_true(x[1].size() == 87);
+      expect_true(x[1].stride() == 1);
+    }
   }
 
   test_that("row based subsetting works") {
     auto getExportedValue = cpp11::package("base")["getExportedValue"];
 
-    cpp11::doubles_matrix x(SEXP(getExportedValue("datasets", "volcano")));
+    cpp11::doubles_matrix<cpp11::by_row> x(getExportedValue("datasets", "volcano"));
+    expect_true(x.nslices() == 87);
+    expect_true(x.slice_size() == 61);
 
     auto r = x[0];
     expect_true(r[0] == 100);
     expect_true(r[60] == 103);
+
+    auto r2 = x[2];
+    expect_true(r2[0] == 102);
+    expect_true(r2[60] == 104);
+  }
+
+  test_that("column based subsetting works") {
+    auto getExportedValue = cpp11::package("base")["getExportedValue"];
+
+    cpp11::doubles_matrix<cpp11::by_column> x(getExportedValue("datasets", "volcano"));
+    expect_true(x.nslices() == 61);
+    expect_true(x.slice_size() == 87);
+
+    auto c = x[0];
+    expect_true(c[0] == 100);
+    expect_true(c[86] == 97);
+
+    auto c2 = x[5];
+    expect_true(c2[0] == 101);
+    expect_true(c2[86] == 99);
   }
 
   test_that("index based subsetting works") {
     auto getExportedValue = cpp11::package("base")["getExportedValue"];
 
-    cpp11::doubles_matrix x(SEXP(getExportedValue("datasets", "volcano")));
+    cpp11::doubles_matrix<cpp11::by_row> xr(getExportedValue("datasets", "volcano"));
+    expect_true(xr(0, 0) == 100);
+    expect_true(xr(0, 60) == 103);
+    expect_true(xr(10, 13) == 121);
 
-    expect_true(x(0, 0) == 100);
-    expect_true(x(0, 60) == 103);
+    cpp11::doubles_matrix<cpp11::by_column> xc(getExportedValue("datasets", "volcano"));
+    expect_true(xc(0, 0) == 100);
+    expect_true(xc(0, 60) == 103);
+    expect_true(xc(10, 13) == 121);
   }
 
   test_that("copy constructor works") {
     auto getExportedValue = cpp11::package("base")["getExportedValue"];
-    cpp11::doubles_matrix x(SEXP(getExportedValue("datasets", "volcano")));
-    cpp11::doubles_matrix y(x);
+    cpp11::doubles_matrix<cpp11::by_row> x(getExportedValue("datasets", "volcano"));
+
+    cpp11::doubles_matrix<cpp11::by_row> yr(x);
+    expect_true(x.nrow() == yr.nrow());
+    expect_true(x.ncol() == yr.ncol());
+    expect_true(yr.nslices() == yr.nrow());
+    expect_true(SEXP(x) == SEXP(yr));
 
-    expect_true(x.nrow() == y.nrow());
-    expect_true(SEXP(x) == SEXP(y));
+    cpp11::doubles_matrix<cpp11::by_column> yc(x);
+    expect_true(x.nrow() == yc.nrow());
+    expect_true(x.ncol() == yc.ncol());
+    expect_true(yc.nslices() == yc.ncol());
+    expect_true(SEXP(x) == SEXP(yc));
   }
 }

cpp11test/tests/testthat/test-matrix.R

@@ -3,7 +3,23 @@ test_that("row_sums gives same result as rowSums", {
   expect_equal(row_sums(x), rowSums(x))
 
   # With missing values
-  y <- cbind(x1 = 3, x2 = c(4:1, 2:5))
-  y[3, ] <- NA; x[4, 2] <- NA
+  x[4, 2] <- NA
   expect_equal(row_sums(x), rowSums(x))
+
+  y <- cbind(x1 = 3, x2 = c(4:1, 2:5))
+  y[3, ] <- NA;
+  expect_equal(row_sums(y), rowSums(y))
+})
+
+test_that("col_sums gives same result as colSums", {
+  x <- cbind(3, c(4:1, 2:5))
+  expect_equal(col_sums(x), colSums(x))
+
+  # With missing values
+  x[4, 2] <- NA
+  expect_equal(col_sums(x), colSums(x))
+
+  y <- cbind(3, c(4:1, 2:5))
+  y[3, ] <- NA;
+  expect_equal(col_sums(y), colSums(y))
 })