add polynomial mutation operator (in C)

R/mutator.polynomial.R

@@ -0,0 +1,36 @@
+#' @title
+#' Generator of the polynomial mutation operator.
+#'
+#' @description
+#' Performs an polynomial mutation as used in the SMS-EMOA algorithm.
+#'
+#' @param p [\code{numeric(1)}]\cr
+#'   Probability of mutation of each gene.
+#' @param eta [\code{numeric(1)}\cr
+#'   Distance parameter of the mutation distribution.
+#' @return [\code{ecr_mutator}]
+#' @family mutators
+#' @export
+makePolynomialMutation = function(p = 0.2, eta = 10) {
+  assertNumber(p, lower = 0, upper = 1, na.ok = FALSE)
+  assertNumber(eta, lower = 1, finite = TRUE, na.ok = FALSE)
+
+  force(p)
+  force(eta)
+
+  mutator = function(ind, task, control) {
+    par.set = task$par.set
+    lower = getLower(par.set)
+    upper = getUpper(par.set)
+    child = .Call("polynomialMutationC", ind, lower, upper, p, eta)
+    return(child)
+  }
+
+  makeMutator(
+    mutator = mutator,
+    name = "Polynomial mutation operator",
+    description = "Apply polynomial mutation to an individual.",
+    supported = "float",
+    params = list(p = p, eta = eta)
+  )
+}

src/polynomialMutation.c

@@ -0,0 +1,61 @@
+#include <R.h>
+#include <Rinternals.h>
+
+#include "macros.h"
+#include "helpers.h"
+
+/* Polynomial mutation operator.
+ *
+ * Performs polynomial mutation of an individual with probability p for each gene.
+ *
+ * @param x [numeric]
+ *   Individual.
+ * @param lower [numeric]
+ *   Lower box constraints.
+ * @param upper [numeric]
+ *   Upper box constraints.
+ * @param p [numeric(1)]
+ *   Probability for mutation of each gene.
+ * @param eta [numeric(1)]
+ *   Parameter eta.
+ * @return [numeric] Mutated individual.
+ */
+SEXP polynomialMutationC(SEXP x, SEXP lower, SEXP upper, SEXP p, SEXP eta) {
+  SEXP res;
+
+  // unpack imcoming R data
+  EXTRACT_NUMERIC_VECTOR(x, c_x, d);
+  EXTRACT_NUMERIC_VECTOR(lower, c_lower, c_lower_length);
+  EXTRACT_NUMERIC_VECTOR(upper, c_upper, c_upper_length);
+  EXTRACT_REAL(p, c_p);
+  EXTRACT_REAL(eta, c_eta);
+
+  const double mpow = 1.0 / (c_eta + 1.0);
+
+  PROTECT(res = allocVector(REALSXP, d));
+  double *c_res = REAL(res);
+
+  GetRNGstate();
+  double rnd, deltaq;
+  for (int i = 0; i < d; ++i) {
+    // do all the complicated stuff here
+    if (unif_rand() < c_p) {
+      const double delta = c_upper[i] - c_lower[i];
+      rnd = unif_rand();
+      if (rnd <= 0.5) {
+        const double xy = 1.0 - (c_x[i] - c_lower[i]) / delta;
+        deltaq = pow(2.0 * rnd + (1.0 - 2.0 * rnd) * pow(xy, c_eta + 1.0), mpow) - 1.0;
+      } else {
+        const double xy = 1.0 - (c_upper[i] - c_x[i]) / delta;
+        deltaq = 1.0 - pow(2.0 * (1.0 - rnd) + 2.0 * (rnd - 0.5) * pow(xy, c_eta + 1.0), mpow);
+      }
+      c_res[i] = forceToBounds(c_x[i] + deltaq * delta, c_lower[i], c_upper[i]);
+    } else {
+      c_res[i] = c_x[i];
+    }
+  }
+
+  PutRNGstate();
+  UNPROTECT(1); /* res */
+  return(res);
+}