Move which_within_distance to C++ (#102)

* First attempt

* Don't coerce matrix to vector

* Use writable namespace

* Try this I guess

* Use cpp11 classes throughout

* Use STL instead

New commit instead of revert to also fix whitespace

* Suggestions from code review

Author: mikemahoney218

DESCRIPTION

@@ -57,3 +57,6 @@ Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.0
+LinkingTo: 
+    cpp11
+SystemRequirements: C++11

NAMESPACE

@@ -78,3 +78,4 @@ importFrom(stats,cutree)
 importFrom(stats,dist)
 importFrom(stats,hclust)
 importFrom(stats,kmeans)
+useDynLib(spatialsample, .registration = TRUE)

R/0_imports.R

@@ -2,6 +2,7 @@
 "_PACKAGE"
 
 ## usethis namespace: start
+#' @useDynLib spatialsample, .registration = TRUE
 ## usethis namespace: end
 NULL
 

R/buffer.R

@@ -44,6 +44,9 @@ buffered_complement <- function(ind, buff_ind, n) {
 
 row_ids_within_dist <- function(distmat, indices, dist) {
   if (units::set_units(dist, NULL) > 0) {
+    # c++ won't implicitly cast, so do it in R
+    mode(distmat) <- "numeric"
+    dist <- as.numeric(dist)
     purrr::map(
       # indices is the output of split_unnamed
       indices,
@@ -55,22 +58,10 @@ row_ids_within_dist <- function(distmat, indices, dist) {
       # not sorted in
       #
       # So here we append the new indices to the old and de-duplicate them
-      ~ unique(c(.x, which_within_dist(distmat, .x, dist)))
+      ~ unique(c(.x, which_within_dist(distmat, as.numeric(.x), dist)))
     )
   } else {
     # initialize to integer(0) in case buffer is <= 0:
     lapply(seq_along(indices), function(x) integer(0))
   }
 }
-
-# Return row IDs for which elements of `data` are within `dist` of `data[idx, ]`
-# Note that data[idx, ] are within any positive distance of themselves
-# and as such are returned by this function
-which_within_dist <- function(distmat, idx, dist) {
-  unlist(
-    purrr::map(
-      idx,
-      ~ which(distmat[.x, ] <= dist)
-    )
-  )
-}

R/cpp11.R

@@ -0,0 +1,5 @@
+# Generated by cpp11: do not edit by hand
+
+which_within_dist <- function(distmat, idx, dist) {
+  .Call(`_spatialsample_which_within_dist`, distmat, idx, dist)
+}

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll

src/code.cpp

@@ -0,0 +1,34 @@
+#include <cpp11.hpp>
+#include <vector>
+using namespace cpp11;
+
+[[cpp11::register]]
+cpp11::writable::integers which_within_dist(doubles_matrix<> distmat, doubles idx, double dist) {
+
+  int n_idx = idx.size();
+  int n_matrix = distmat.ncol();
+  std::vector<bool> comparisons(n_matrix);
+  int cur_row;
+
+  for (int i = 0; i < n_idx; i++) {
+    cur_row = idx[i] - 1;
+    for (int j = 0; j < n_matrix; j++) {
+      if (distmat(cur_row, j) <= dist) {
+        comparisons[j] = true;
+      }
+    }
+  }
+
+  auto n_pos = std::count(comparisons.begin(), comparisons.end(), true);
+  std::vector<int> out(n_pos);
+  int cur_idx = 0;
+  for (int i = 0; i < n_matrix; i++) {
+    if (comparisons[i]) {
+      out[cur_idx] = i + 1;
+      ++cur_idx;
+    }
+  }
+
+  return out;
+
+}

src/cpp11.cpp

@@ -0,0 +1,27 @@
+// Generated by cpp11: do not edit by hand
+// clang-format off
+
+
+#include "cpp11/declarations.hpp"
+#include <R_ext/Visibility.h>
+
+// code.cpp
+cpp11::writable::integers which_within_dist(doubles_matrix<> distmat, doubles idx, double dist);
+extern "C" SEXP _spatialsample_which_within_dist(SEXP distmat, SEXP idx, SEXP dist) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(which_within_dist(cpp11::as_cpp<cpp11::decay_t<doubles_matrix<>>>(distmat), cpp11::as_cpp<cpp11::decay_t<doubles>>(idx), cpp11::as_cpp<cpp11::decay_t<double>>(dist)));
+  END_CPP11
+}
+
+extern "C" {
+static const R_CallMethodDef CallEntries[] = {
+    {"_spatialsample_which_within_dist", (DL_FUNC) &_spatialsample_which_within_dist, 3},
+    {NULL, NULL, 0}
+};
+}
+
+extern "C" attribute_visible void R_init_spatialsample(DllInfo* dll){
+  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+  R_useDynamicSymbols(dll, FALSE);
+  R_forceSymbols(dll, TRUE);
+}