c code for is_dominated

.gitignore

@@ -30,3 +30,5 @@ cran-comments.md
 .vscode/
 docs/
 .vscode
+/src/*.so
+/src/*.o

NAMESPACE

@@ -27,3 +27,4 @@ import(data.table)
 import(mlr3misc)
 import(paradox)
 importFrom(R6,R6Class)
+useDynLib(bbotk,c_is_dominated)

R/helper.R

@@ -9,18 +9,38 @@ terminated_error = function(optim_instance) {
     c("terminated_error", "error", "condition"))
 }
 
-# r implementation of emoa do_is_dominated (emoa/src/dominance.c)
-# for 1000x2 matrix it is actually faster then going into c
-#
 #' @title Calculate which points are dominated
 #' @description
 #' A slow implementation that calculates which points are not dominated,
 #' i.e. points that belong to the Pareto front.
 #'
 #' @param ymat (`matrix()`) \cr
 #'   A numeric matrix. Each column (!) contains one point.
+#' @useDynLib bbotk c_is_dominated
 #' @export
 is_dominated = function(ymat) {
+  assert_matrix(ymat, mode = "double")
+  .Call(c_is_dominated, ymat, PACKAGE = "bbotk")
+}
+
+if (FALSE) {
+  ymat = matrix(c(1, 2, 2, 1, 2, 2), nrow = 2)
+  print(ymat)
+  is_dominated(ymat)
+
+  ymat = matrix(runif(2 * 100), nrow = 2)
+  is_dominated(ymat)
+  ymat
+  bm(is_dominated(ymat), is_dominated_r(ymat))
+
+
+  ymat = matrix(c(0.1, 0.3, 0.2, 0.4, 0.3, 0.1), nrow = 2)
+  ymat
+  is_dominated(ymat)
+
+}
+
+is_dominated_r = function(ymat) {
   #FIXME Replace with C implementation
   # implementation based emoa/src/dominance.c
 

src/init.c

@@ -0,0 +1,17 @@
+#include <R.h>
+#include <Rinternals.h>
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+/* .Call calls */
+extern SEXP c_is_dominated(SEXP);
+
+static const R_CallMethodDef CallEntries[] = {
+    {"c_is_dominated", (DL_FUNC) &c_is_dominated, 1},
+    {NULL, NULL, 0}
+};
+
+void R_init_bbotk(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}

src/is_dominated.c

@@ -0,0 +1,67 @@
+#include <R.h>
+#include <Rinternals.h>
+
+/*
+ * Adapted from https://github.com/olafmersmann/emoa/blob/master/src/dominance.c
+ * written by Olaf Mersmann (OME) <[email protected]>.
+ */
+
+static R_INLINE int dominates(const double * x, const double * y, const R_len_t d) {
+    Rboolean x_flag = 0, y_flag = 0;
+
+    for (R_len_t k = 0; k < d; k++) {
+        if (x[k] < y[k]) {
+            y_flag = 1; // y cannot dominate x
+        } else if (y[k] < x[k]) {
+            x_flag = 1; // x cannot dominate y
+        }
+
+        // Note that we could break as soon as both x_flag and y_flag are true.
+        // However, we assume that the number of dimensions d is usually small
+        // so it is probably faster to just walk over all components
+    }
+
+    return y_flag - x_flag;
+}
+
+SEXP c_is_dominated(SEXP p_) {
+    // accessors for input matrix p_
+    const R_len_t n = ncols(p_);
+    const R_len_t d = nrows(p_);
+    const double * p = REAL(p_);
+
+    // accessors for output vector res_
+    SEXP res_ = PROTECT(allocVector(LGLSXP, n));
+    int * res = LOGICAL(res_);
+    for (R_len_t i = 0; i < n; i++) res[i] = FALSE;
+
+    // iterate over all columns of input
+    for (R_len_t i = 0; i < n; i++) {
+        if (res[i]) {
+            // current column was marked as dominated in a
+            // previous iteration; skip
+            continue;
+        }
+
+        // find a non-dominated column to compare with
+        for (R_len_t j = i + 1; j < n; j++) {
+            if (res[j]) {
+                continue;
+            }
+
+            // compare vector p[, i] with vector p[, j]
+            int dom = dominates(p + (i * d), p + (j * d), d);
+
+            if (dom > 0) {
+                // i dominates j
+                res[j] = TRUE;
+            } else if (dom < 0) {
+                // j dominates i
+                res[i] = TRUE;
+            }
+        }
+    }
+
+    UNPROTECT(1);
+    return res_;
+}