ecr on fire: non-dominated sorting in C

R/doNondominatedSorting.R

@@ -23,6 +23,14 @@
 #'   }
 #' @export
 doNondominatedSorting = function(x) {
+  assertMatrix(x, min.rows = 2L, min.cols = 2L, any.missing = FALSE, all.missing = FALSE, mode = "numeric")
+  ranks = .Call("doNondominatedSortingC", x)
+  #FIXME: return ranks also
+  return(list(ranks = ranks, dom.counter = rep(1L, ncol(x))))
+}
+
+# Old pure R implementation
+doNondominatedSortingR = function(x) {
   # initialize domination front wrapper
   fronts = list()
   fronts[[1L]] = list()

src/dominates.c

@@ -85,3 +85,87 @@ SEXP dominatedC(SEXP r_points) {
   UNPROTECT(1);
   return(r_res);
 }
+
+/*
+ * Nondominated sorting algorithm as used in NSGA EMOA.
+ *
+ * @description
+ *  Returns a list of two vectors:
+ *  - front (integer): Ranks/front of the corresponding points.
+ *  - domcount (integer): Number of points the corresponding points dominates.
+ *
+ * @param r_points [matrix]
+ *   Numeric (n x m) matrix.
+ * @return [list]
+ */
+SEXP doNondominatedSortingC(SEXP r_points) {
+  // unwrap R structure into C objects
+  EXTRACT_NUMERIC_MATRIX(r_points, c_points, dim, n_points);
+
+  // rank, i.e., number of front
+  SEXP r_rank = ALLOC_INTEGER_VECTOR(n_points);
+  int* c_rank = INTEGER(r_rank);
+
+  for (int i = 0; i < n_points; ++i) {
+    c_rank[i] = 0;
+  }
+
+  // count how many points are already correctly sorted
+  int n_sorted = 0;
+
+  // count number of how any individuals i is dominated
+  SEXP r_domcounter = ALLOC_INTEGER_VECTOR(n_points);
+  int* c_domcounter = INTEGER(r_domcounter);
+
+  for (int i = 0; i < n_points; ++i) {
+    c_domcounter[i] = 0;
+  }
+
+  for (unsigned int i = 0; i < n_points; ++i) {
+    for (unsigned int j = (i + 1); j < n_points; ++j) {
+      int isDominated = dominance(c_points, i, j, dim);
+      if (isDominated > 0) {
+        //FIXME: save which point(s) are dominated by which
+        ++c_domcounter[j];
+      } else if (isDominated < 0) {
+        ++c_domcounter[i];
+      }
+    }
+
+    // assign rank 1 to all points that are nondominated
+    if (c_domcounter[i] == 0) {
+      c_rank[i] = 1;
+      ++n_sorted;
+    }
+  }
+
+  unsigned int cur_rank = 1;
+  while (n_sorted < n_points) {
+    //Rprintf("Sorted: %i of %i\n", n_sorted, n_points);
+    // iterate over all points and search for the ones with the currently
+    // "active" rank
+    for (int i = 0; i < n_points; ++i) {
+      if (c_rank[i] != cur_rank) {
+        continue;
+      }
+      // otherwise check if current point dominates another one
+      for (int j = 0; j < n_points; ++j) {
+        //FIXME: replace with non-redundant dominance check
+        if (dominance(c_points, i, j, dim) > 0) { // i dominates j
+          --c_domcounter[j];
+          if (c_domcounter[j] == 0) { /* now on first front */
+            c_rank[j] = (cur_rank + 1);
+            ++n_sorted;
+          }
+        }
+      }
+    }
+    ++cur_rank;
+  }
+
+  //FIXME: make copy of c_domcounter
+  UNPROTECT(1);
+  //Free(n_sorted);
+  //Free(c_domcounter);
+  return(r_rank);
+}