add simulated binary crossover recombinator (in C)

R/recombinator.sbx.R

@@ -0,0 +1,45 @@
+#' @title
+#' Generator of the Simulated Binary Crossover (SBX) recombinator.
+#'
+#' @description
+#' ...
+#FIXME: add docs
+#'
+#' @param eta [\code{numeric(1)}]\cr
+#'   Parameter eta, i.e., the distance parameters of the crossover distribution.
+#' @param p [\code{numeric(1)}]\cr
+#'   Crossover probability for each gene. Default is \code{1.0}.
+#' @return [\code{ecr_recombinator}]
+#' @family recombinators
+#' @export
+makeSBXRecombinator = function(eta = 5, p = 1.0) {
+  assertNumber(eta, lower = 1, na.ok = FALSE)
+  assertNumber(p, lower = 0, upper = 1, na.ok = FALSE)
+
+  force(eta)
+  force(p)
+
+  recombinator = function(inds, task, control) {
+    par.set = task$par.set
+    n.params = getParamLengths(par.set)
+    lower = getLower(par.set)
+    upper = getUpper(par.set)
+
+    # convert parents to d x 2 matrix for C
+    inds = do.call(cbind, inds)
+
+    # SBX produces two children
+    children = .Call("simulatedBinaryCrossoverC", inds, as.numeric(lower), as.numeric(upper), p, eta)
+
+    return(wrapChildren(children[, 1L], children[, 2L]))
+  }
+
+  makeRecombinator(
+    recombinator = recombinator,
+    name = "Simulated Binary Crossover (SBX) recombinator",
+    description = "Performs simulated binary crossover.",
+    n.parents = 2L,
+    supported = c("float"),
+    params = list(p = p, eta = eta)
+  )
+}

src/helpers.c

@@ -0,0 +1,18 @@
+/* Force value to stay within box constraints.
+ *
+ * @param x [numeric(1)]
+ *   Numeric value.
+ * @param lower [numeric(1)]
+ *   Lower box constraint.
+ * @param upper [numeric(1)]
+ *   Upper box constraint.
+ * @return [numeric(1)]
+ */
+double forceToBounds(const double x, const double lower, const double upper) {
+  if (x < lower) {
+    return(lower);
+  } else if (x > upper) {
+    return(upper);
+  }
+  return(x);
+}

src/helpers.h

@@ -0,0 +1,6 @@
+#ifndef ECR_HELPERS
+#define ECR_HELPERS
+
+double forceToBounds(const double x, const double lower, const double upper);
+
+#endif

src/macros.h

@@ -1,10 +1,11 @@
 #ifndef SEXP_MACROS
 #define SEXP_MACROS
 
-#include <R.h>
-#include <Rinternals.h>
+// basic math helpers
+#define MIN(A, B) ((A < B) ? (A) : (B))
+#define MAX(A, B) ((A > B) ? (A) : (B))
 
-// define some macros
+// R specific macros
 #define EXTRACT_NUMERIC_MATRIX(S_EXP, C_DATA, N_ROW, N_COL) \
     double *C_DATA = REAL(S_EXP); \
     const R_len_t N_ROW = nrows(S_EXP); \
@@ -17,11 +18,12 @@
 #define EXTRACT_INTEGER(S_EXP, I) \
     int I = INTEGER(S_EXP)[0];
 
+#define EXTRACT_REAL(S_EXP, D) \
+    double D = REAL(S_EXP)[0];
+
 #define ALLOC_VECTOR(type, size) (PROTECT(allocVector(type, size)))
 #define ALLOC_REAL_VECTOR(size) (ALLOC_VECTOR(REALSXP, size))
 #define ALLOC_INTEGER_VECTOR(size) (ALLOC_VECTOR(INTSXP, size))
 #define ALLOC_LIST(size) (ALLOC_VECTOR(VECSXP, size))
 
-
-
 #endif

src/simulatedBinaryCrossover.c

@@ -0,0 +1,92 @@
+#include <R.h>
+#include <Rinternals.h>
+
+#include "macros.h"
+#include "helpers.h"
+
+/* Calculate betaq parameter based on beta and eta.
+ *
+ * @param beta [numeric(1)]
+ *   Beta parameter.
+ * @param eta [numeric(1)]
+ *   Eta parameter.
+ * @return betaq
+ */
+static R_INLINE double calculateBetaQ(const double beta, const double eta) {
+  const double rand = unif_rand();
+  const double alpha = 2.0 - pow(beta, -(eta + 1.0));
+  if (rand <= (1.0 / alpha)) {
+    return(pow(rand * alpha, 1.0 / (eta + 1.0)));
+  }
+  return(pow(1.0 / (2.0 - rand * alpha), 1.0 / (eta + 1.0)));
+}
+
+/* Simulated Binary Crossover (SBX).
+ *
+ * Given two individuals a simulated binary crossover is performed.
+ * The operator returns two new individuals.
+ *
+ * @param parents [matrix]
+ *   Parent genomes as a matrix.
+ * @param lower [numeric]
+ *   Vector of lower bounds.
+ * @param upper [numeric]
+ *   Vector of upper bounds.
+ * @param p [numeric(1)]
+ *   Probability in [0,1].
+ * @param eta [numeric(1)]
+ *   Paramter eta.
+ * @return [matrix]
+ */
+SEXP simulatedBinaryCrossoverC(SEXP parents, SEXP lower, SEXP upper, SEXP p, SEXP eta) {
+  SEXP res;
+  double betaq;
+
+  // unpack imcoming R data
+  EXTRACT_NUMERIC_MATRIX(parents, c_parents, d, n);
+  const double *parent1 = c_parents;
+  const double *parent2 = c_parents + d; /* in C we need to compute rows by hand */
+  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);
+
+  /* Allocate result matrix: */
+  PROTECT(res = allocMatrix(REALSXP, d, 2));
+  double *child1 = REAL(res);
+  double *child2 = REAL(res) + d;
+
+  GetRNGstate();
+  for (int i = 0; i < d; ++i) {
+  /* Crossover for dim i is performed with probability p, if
+   * the two parents differ by at least a very small number.
+   */
+    if (unif_rand() <= c_p && fabs(parent1[i] - parent2[i]) > 1.0e-14) {
+      const double y1 = MIN(parent1[i], parent2[i]);
+      const double y2 = MAX(parent1[i], parent2[i]);
+      const double yl = c_lower[i];
+      const double yu = c_upper[i];
+
+      /* Calculate offsprint: */
+      betaq = calculateBetaQ(1.0 + (2.0 * (y1 - yl) / (y2 - y1)), c_eta);
+      const double c1 = forceToBounds(0.5 * ((y1 + y2) - betaq * (y2 - y1)), yl, yu);
+
+      betaq = calculateBetaQ(1.0 + (2.0 * (yu - y2) / (y2 - y1)), c_eta);
+      const double c2 = forceToBounds(0.5 * ((y1 + y2) + betaq * (y2 - y1)), yl, yu);
+
+      if (unif_rand() > 0.5) {
+        child1[i] = c2;
+        child2[i] = c1;
+      } else {
+        child1[i] = c1;
+        child2[i] = c2;
+      }
+    } else {
+      child1[i] = parent1[i];
+      child2[i] = parent2[i];
+    }
+  }
+  PutRNGstate();
+  UNPROTECT(1); /* res */
+  return(res);
+}