diff --git a/R/RcppExports.R b/R/RcppExports.R index d0108ac5d6..6f814c9307 100644 --- a/R/RcppExports.R +++ b/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) } diff --git a/R/misc.R b/R/misc.R index 1fc1db6986..5c41223422 100644 --- a/R/misc.R +++ b/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") diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 44542cfd41..86eef4dee5 100644 --- a/src/RcppExports.cpp +++ b/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) { diff --git a/src/matrix_smooth.c b/src/matrix_smooth.c deleted file mode 100644 index eaefa00b75..0000000000 --- a/src/matrix_smooth.c +++ /dev/null @@ -1,46 +0,0 @@ -#include -#include -#include - -SEXP matrix_smooth(SEXP mat) -{ -#define ij(i, j) ((i) + (ni) * (j)) - /* Note: the 2d data are stored in column order */ - SEXP res; - int ni = INTEGER(GET_DIM(mat))[0]; - int nj = INTEGER(GET_DIM(mat))[1]; - int i, j; - double *matp, *resp; - if (!isMatrix(mat)) - error("'mat' must be a matrix"); - if (!isReal(mat)) - error("'mat' must be numeric, not integer"); - matp = REAL(mat); - if (length(mat) != ni * nj) - error("'ni'*'nj' must equal number of elements in 'mat'"); - PROTECT(res = allocMatrix(REALSXP, ni, nj)); - resp = REAL(res); - for (i = 0; i < ni*nj; i++) - resp[i] = 99.99; - // copy edges (FIXME: coiuld use 1D smoother here) - for (j = 0; j < nj; j++) { - resp[ij(0, j)] = matp[ij(0, j)]; - resp[ij(ni-1, j)] = matp[ij(ni-1, j)]; - } - for (i = 0; i < ni; i++) { - resp[ij(i, 0)] = matp[ij(i, 0)]; - resp[ij(i, nj-1)] = matp[ij(i, nj-1)]; - } - // smooth middle - for (i = 1; i < ni - 1; i++) - for (j = 1; j < nj - 1; j++) - resp[ij(i, j)] = (2.0*matp[ij(i, j)] + - matp[ij(i-1, j)] + - matp[ij(i+1, j)] + - matp[ij(i, j-1)] + - matp[ij(i, j+1)]) / 6.0; - UNPROTECT(1); - return(res); -#undef ix -} - diff --git a/src/matrix_smooth.cpp b/src/matrix_smooth.cpp new file mode 100644 index 0000000000..55cf8f6f56 --- /dev/null +++ b/src/matrix_smooth.cpp @@ -0,0 +1,32 @@ +#include +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); +} + diff --git a/src/registerDynamicSymbol.c b/src/registerDynamicSymbol.c index 8c215e5ca9..00a4a6a2a2 100644 --- a/src/registerDynamicSymbol.c +++ b/src/registerDynamicSymbol.c @@ -3,8 +3,10 @@ #include 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} }; diff --git a/src/trap.cpp b/src/trap.cpp index 1214f4910b..c984af234e 100644 --- a/src/trap.cpp +++ b/src/trap.cpp @@ -1,8 +1,12 @@ #include 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) { diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R index 67dee3f530..15359dc641 100644 --- a/tests/testthat/test_misc.R +++ b/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) +}) +