From 741a2732c7a55abd75ba7939374241ff2980b026 Mon Sep 17 00:00:00 2001
From: Jakob Bossek <info@jakobbossek.de>
Date: Sat, 23 Jan 2016 00:21:27 +0100
Subject: [PATCH] add polynomial mutation operator (in C)

---
 R/mutator.polynomial.R        | 36 +++++++++++++++++++++
 man/makePolynomialMutation.Rd | 30 +++++++++++++++++
 src/polynomialMutation.c      | 61 +++++++++++++++++++++++++++++++++++
 3 files changed, 127 insertions(+)
 create mode 100644 R/mutator.polynomial.R
 create mode 100644 man/makePolynomialMutation.Rd
 create mode 100644 src/polynomialMutation.c

diff --git a/R/mutator.polynomial.R b/R/mutator.polynomial.R
new file mode 100644
index 0000000..cd477fc
--- /dev/null
+++ b/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)
+  )
+}
diff --git a/man/makePolynomialMutation.Rd b/man/makePolynomialMutation.Rd
new file mode 100644
index 0000000..b132fd8
--- /dev/null
+++ b/man/makePolynomialMutation.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mutator.polynomial.R
+\name{makePolynomialMutation}
+\alias{makePolynomialMutation}
+\title{Generator of the polynomial mutation operator.}
+\usage{
+makePolynomialMutation(p = 0.2, eta = 10)
+}
+\arguments{
+\item{p}{[\code{numeric(1)}]\cr
+Probability of mutation of each gene.}
+
+\item{eta}{[\code{numeric(1)}\cr
+Distance parameter of the mutation distribution.}
+}
+\value{
+[\code{ecr_mutator}]
+}
+\description{
+Performs an polynomial mutation as used in the SMS-EMOA algorithm.
+}
+\seealso{
+Other mutators: \code{\link{makeBitFlipMutator}},
+  \code{\link{makeGaussMutator}},
+  \code{\link{makeInsertionMutator}},
+  \code{\link{makeScrambleMutator}},
+  \code{\link{makeSwapMutator}},
+  \code{\link{makeUniformMutator}}
+}
+
diff --git a/src/polynomialMutation.c b/src/polynomialMutation.c
new file mode 100644
index 0000000..a4d9c29
--- /dev/null
+++ b/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);
+}