matrix_smooth.c rewritten in C++

R/RcppExports.R

@@ -1,6 +1,10 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+matrix_smooth <- function(mat) {
+    .Call(`_oce_matrix_smooth`, mat)
+}
+
 trap <- function(x, y, type) {
     .Call(`_oce_trap`, x, y, type)
 }

R/misc.R

@@ -3906,7 +3906,8 @@ matrixSmooth <- function(m, passes=1)
     storage.mode(m) <- "double"
     if (passes > 0) {
         for (pass in seq.int(1, passes, 1)) {
-            m <- .Call("matrix_smooth", m)
+            message("pass=", pass)
+            m <- .Call(`_oce_matrix_smooth`, m)
         }
     } else {
         warning("matrixSmooth given passes<=0, so returning matrix unmodified")

src/RcppExports.cpp

@@ -5,6 +5,17 @@
 
 using namespace Rcpp;
 
+// matrix_smooth
+NumericMatrix matrix_smooth(NumericMatrix mat);
+RcppExport SEXP _oce_matrix_smooth(SEXP matSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type mat(matSEXP);
+    rcpp_result_gen = Rcpp::wrap(matrix_smooth(mat));
+    return rcpp_result_gen;
+END_RCPP
+}
 // trap
 NumericVector trap(NumericVector x, NumericVector y, NumericVector type);
 RcppExport SEXP _oce_trap(SEXP xSEXP, SEXP ySEXP, SEXP typeSEXP) {

src/matrix_smooth.cpp

@@ -0,0 +1,32 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// This is a utility function called by matrixSmooth,
+// so it does not get a wrapper inserted into the R namespace;
+// that's why there is no roxygen @export here, and no roxygen
+// documentation, either. However, we need to tell Rcpp to export
+// it, so matrixSmooth() can .Call() it.
+//
+// [[Rcpp::export]]
+
+NumericMatrix matrix_smooth(NumericMatrix mat)
+{
+    int ni = mat.nrow(), nj=mat.ncol();
+    NumericMatrix res(ni, nj);
+    // copy edges (FIXME: could use 1D smoother here)
+    for (int j = 0; j < nj; j++) {
+        res(0, j) = mat(0, j);
+        res(ni-1, j) = mat(ni-1, j);
+    }
+    for (int i = 0; i < ni; i++) {
+        res(i, 0) = mat(i, 0);
+        res(i, nj-1) = mat(i, nj-1);
+    }
+    // smooth middle 
+    for (int i = 1; i < ni - 1; i++)
+        for (int j = 1; j < nj - 1; j++)
+            res(i, j) = (2.0*mat(i, j) +
+                    mat(i-1, j) + mat(i+1, j) + mat(i, j-1) + mat(i, j+1)) / 6.0;
+    return(res);
+}
+

src/registerDynamicSymbol.c

@@ -3,8 +3,10 @@
 #include <R_ext/Rdynload.h>
 
 extern SEXP _oce_trap(SEXP, SEXP, SEXP);
+extern SEXP _oce_matrix_smooth(SEXP);
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_oce_matrix_smooth", (DL_FUNC) &_oce_matrix_smooth, 1},
     {"_oce_trap", (DL_FUNC) &_oce_trap, 3},
     {NULL, NULL, 0}
 };

src/trap.cpp

@@ -1,8 +1,12 @@
 #include <Rcpp.h>
 using namespace Rcpp;
 
-// This is a utility function called by \code{\link{integrateTrapezoid}},
-// not intended for direct use.
+// This is a utility function called by integrateTrapeziod(),
+// so it does not get a wrapper inserted into the R namespace;
+// that's why there is no roxygen @export here, and no roxygen
+// documentation, either. However, we need to tell Rcpp to export
+// it, so integreateTrapezoid() can .Call() it.
+//
 // [[Rcpp::export]]
 NumericVector trap(NumericVector x, NumericVector y, NumericVector type)
 {

tests/testthat/test_misc.R

@@ -239,3 +239,21 @@ test_that("get_bit (unused in oce)", {
           expect_equal(c(0, 0, 1, 1, 1, 0, 1, 0), bits)
 })
 
+test_that("matrixSmooth", {
+          data(volcano)
+          v <- matrixSmooth(volcano)
+          ## These values are from the C version of src/matrix_smooth.c as it
+          ## existed on 2018-03-03, and the point of the present check is to
+          ## test the code rewriting in C++ as src/matrix_smooth.cpp, added
+          ## on that date. The first check is to be sure the edge treatmet
+          ## still works, and the others are somewhat arbitrary, built on the
+          ## assumption that if the C++ translation is faulty, the results
+          ## will be crazy-wrong (e.g. transposed matrix).
+          ve <- matrix(c(100, 101, 102, 100, 101.1666667, 102.1666667, 101,
+                         101.8333333, 102.8333333), byrow=FALSE, nrow=3)
+          expect_equal(ve, v[1:3, 1:3])
+          expect_equal(mean(v), 130.1787262)
+          expect_equal(mean(v[,10]), 121.3429119)
+          expect_equal(mean(v[10,]), 127.9808743)
+})
+