Merge pull request #149 from schaepermeier/mpm2-fast-evaluation

Faster evaluation for MPM2

DESCRIPTION

@@ -19,7 +19,13 @@ Authors@R: c(
     family = "Kerschke",
     email = "pascal.kerschke@tu-dresden.de",
     role = "ctb",
-    comment = c(ORCID = "0000-0003-2862-1418")))
+    comment = c(ORCID = "0000-0003-2862-1418")),
+  person(
+    given = "Lennart",
+    family = "Schäpermeier",
+    email = "lennart.schaepermeier@tu-dresden.de",
+    role = "ctb",
+    comment = c(ORCID = "0000-0003-3929-7465")))
 Maintainer: Jakob Bossek <j.bossek@gmail.com>
 URL: https://jakobbossek.github.io/smoof/, https://github.com/jakobbossek/smoof
 BugReports: https://github.com/jakobbossek/smoof/issues

NEWS

@@ -4,6 +4,8 @@ smoof 1.6.0.4
 * Fixed bug in getLoggedValues when logging of x-values was set to FALSE in addLoggingWrapper
 * Fixed bug with instance ID mapping in makeBiObjBBOBFunction
 * Added support for the extended bi-objective BBOB functions (FIDs 56-92) in makeBiObjBBOBFunction
+* Added extra helper functions to extract all problem instance-specific data from MPM2 generator
+* Added R-based evaluation environment for MPM2 generator which greatly accelerates evaluation in multi-objective settings
 
 smoof 1.6.0.3
 =============

NEWS.md

@@ -3,6 +3,8 @@
 ## New features
 
 * Added support for the extended bi-objective BBOB functions (FIDs 56-92) in makeBiObjBBOBFunction
+* Added extra helper functions to extract all problem instance-specific data from MPM2 generator
+* Added R-based evaluation environment for MPM2 generator which greatly accelerates evaluation in multi-objective settings
 
 ## Bugfixes 
 

R/sof.mpm2.R

@@ -12,6 +12,12 @@
 #'   of elliptically shaped peaks.
 #' @param peak.shape [\code{character(1)}]\cr
 #'   Shape of peak(s). Possible values are \dQuote{ellipse} and \dQuote{sphere}.
+#' @param evaluation.env [\code{character(1)}]\cr
+#'   Evaluation environment after the function was created. Possible
+#'   (case-insensitive) values are \dQuote{R} (default) and \dQuote{Python}. The
+#'   original generation of the problem is always done in the original Python
+#'   environment. However, evaluation in R is faster, especially if multiple
+#'   MPM2 functions are used in a multi-objective setting.
 #' @return [\code{smoof_single_objective_function}]
 #' @examples
 #' \dontrun{
@@ -33,7 +39,8 @@
 #' @author \R interface by Jakob Bossek. Original python code provided by the Simon Wessing.
 #'
 #' @export
-makeMPM2Function = function(n.peaks, dimensions, topology, seed, rotated = TRUE, peak.shape = "ellipse") {
+makeMPM2Function = function(n.peaks, dimensions, topology, seed, rotated = TRUE,
+                            peak.shape = "ellipse", evaluation.env = "R") {
   # if (isWindows()) {
   #  stopf("No support for the multiple peaks model 2 generator at the moment.")
   # }
@@ -64,9 +71,11 @@ makeMPM2Function = function(n.peaks, dimensions, topology, seed, rotated = TRUE,
   # import reticulate namespace
   BBmisc::requirePackages("_reticulate", why = "smoof::makeMultiplePeaksModel2Function")
 
-  # initialize 3 functions from mpm2.py as NULL such that they have a visible binding when checking the pkg
+  # initialize helper functions from mpm2.py as NULL such that they have a visible binding when checking the pkg
   evaluateProblem = getGlobalOptimaParams = getLocalOptimaParams = NULL
-
+  getCovarianceMatrices = getAllPeaks = getAllHeights = getAllShapes = NULL
+  getAllRadii = NULL
+# 
   # load funnel generator to global environment
   eval(reticulate::py_run_file(system.file("mpm2.py", package = "smoof")), envir = .GlobalEnv)
   eval(reticulate::source_python(system.file("mpm2.py", package = "smoof"), envir = .GlobalEnv, convert = TRUE), envir = .GlobalEnv)
@@ -78,19 +87,77 @@ makeMPM2Function = function(n.peaks, dimensions, topology, seed, rotated = TRUE,
   global.opt.params = eval(getGlobalOptimaParams(n.peaks, dimensions, topology, seed, rotated, peak.shape), envir = .GlobalEnv)
   global.opt.params = matrix(global.opt.params[[1L]], nrow = 1L)
 
+  if (tolower(evaluation.env) == "python") {
+    evalFn = function(x) {
+      evaluateProblem(x, n.peaks, dimensions, topology, seed, rotated, peak.shape)
+    }
+  } else if (tolower(evaluation.env) == "r") {
+    # helper functions
+    getPeakMetadata = function(n.peaks, dimensions, topology, seed, rotated, peak.shape) {
+
+      xopt = getAllPeaks(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+      colnames(xopt) = paste0("x", 1:dimensions)
+
+      cov.mats = getCovarianceMatrices(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+      cov = lapply(cov.mats, function(x) matrix(x[[1]], nrow = dimensions, ncol = dimensions))
+
+      height = getAllHeights(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+
+      shape = getAllShapes(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+
+      radius = getAllRadii(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+
+      peak_fns = lapply(1:nrow(xopt), function(i) {
+        createPeakFunction(cov[[i]], xopt[i,], height[i], shape[i], radius[i])
+      })
+
+      fn = function(x) {
+        min(sapply(peak_fns, function(f) f(x)))
+      }
+
+      list(
+        xopt = xopt,
+        cov = cov,
+        height = height,
+        shape = shape,
+        radius = radius,
+        peak_fns = peak_fns,
+        fn = fn
+      )
+    }
+
+    createPeakFunction = function(cov, xopt, height, shape, radius) {
+      function(x) {
+        md = sqrt(t(x - xopt) %*% cov %*% (x - xopt))
+        g = height / (1 + md**shape / radius)
+        return(1 - g)
+      }
+    }
+
+    peakData = getPeakMetadata(n.peaks, dimensions, topology, seed, rotated, peak.shape)
+    evalFn = function(x) {
+      if (is.matrix(x)) {
+        apply(x, 2, peakData$fn)
+      } else {
+        peakData$fn(x)
+      }
+    }
+  } else {
+    warning(paste0("Unknown evaluation.env \"", evaluation.env, "\""))
+  }
+
   smoof.fn = makeSingleObjectiveFunction(
     name = sprintf("Funnel_%i_%i_%i_%s_%s%s", n.peaks, dimensions, seed, topology, peak.shape, ifelse(rotated, "_rotated", "")),
     description = sprintf("Funnel-like function\n(n.peaks: %i, dimension: %i, topology: %s, seed: %i, rotated: %s, shape: %s)",
-      n.peaks, dimensions, topology, seed, rotated, peak.shape),
-    fn = function(x) {
-      evaluateProblem(x, n.peaks, dimensions, topology, seed, rotated, peak.shape)
-    },
+                          n.peaks, dimensions, topology, seed, rotated, peak.shape),
+    fn = evalFn,
     par.set = par.set,
     vectorized = TRUE,
     tags = c("non-separable", "scalable", "continuous", "multimodal"),
     local.opt.params = local.opt.params,
     global.opt.params = global.opt.params
   )
+
   return(smoof.fn)
 }
 

inst/mpm2.py

@@ -280,6 +280,19 @@ def getCovMatrices(self):
             X = X.reshape(len(X) * len(X[0]))
             res.append(np.matrix(X).tolist())
         return res
+
+    def getAllPeaks(self):
+        return np.vstack(self.peaks)
+
+    def getAllHeights(self):
+        return [peak.height for peak in self.peaks]
+
+    def getAllShapes(self):
+        return [peak.shape for peak in self.peaks]
+
+    def getAllRadii(self):
+        return [peak.radius for peak in self.peaks]
+
 
 
 
@@ -335,6 +348,26 @@ def getCovarianceMatrices(npeaks, dimension, topology, randomSeed, rotated, peak
     initProblem(npeaks, dimension, topology, randomSeed, rotated, peakShape)
     return currentProblem.getCovMatrices()
 
+def getAllPeaks(npeaks, dimension, topology, randomSeed, rotated, peakShape):
+    global currentProblem
+    initProblem(npeaks, dimension, topology, randomSeed, rotated, peakShape)
+    return currentProblem.getAllPeaks()
+
+def getAllHeights(npeaks, dimension, topology, randomSeed, rotated, peakShape):
+    global currentProblem
+    initProblem(npeaks, dimension, topology, randomSeed, rotated, peakShape)
+    return currentProblem.getAllHeights()
+
+def getAllShapes(npeaks, dimension, topology, randomSeed, rotated, peakShape):
+    global currentProblem
+    initProblem(npeaks, dimension, topology, randomSeed, rotated, peakShape)
+    return currentProblem.getAllShapes()
+
+def getAllRadii(npeaks, dimension, topology, randomSeed, rotated, peakShape):
+    global currentProblem
+    initProblem(npeaks, dimension, topology, randomSeed, rotated, peakShape)
+    return currentProblem.getAllRadii()
+
 if __name__ == "__main__":
     # nothing to do here
     pass

tests/testthat/test_soofuns.R

@@ -74,10 +74,42 @@ test_that("Multiple peaks model 2 (MPM2) functions work", {
         for (topology in c("funnel", "random")) {
           for (rotated in c(TRUE, FALSE)) {
             for (peak.shape in c("ellipse", "sphere")) {
-              fn = makeMPM2Function(n.peaks = n.peaks, dimension = dimension, topology = topology, seed = 123, rotated = rotated, peak.shape = peak.shape)
-              expect_is(fn, "smoof_single_objective_function")
-              y = fn(rep(0.1, dimension))
-              expect_true(is.numeric(y))
+              fnp = makeMPM2Function(n.peaks = n.peaks, dimensions = dimension,
+                                     topology = topology, seed = 123, rotated = rotated,
+                                     peak.shape = peak.shape, evaluation.env = "Python")
+              fnr = makeMPM2Function(n.peaks = n.peaks, dimensions = dimension,
+                                     topology = topology, seed = 123, rotated = rotated,
+                                     peak.shape = peak.shape, evaluation.env = "R")
+
+              # confirm that both evaluation environments can be created
+              expect_is(fnp, "smoof_single_objective_function")
+              yp = fnp(rep(0.1, dimension))
+              expect_true(is.numeric(yp))
+
+              expect_is(fnr, "smoof_single_objective_function")
+              yr = fnr(rep(0.1, dimension))
+              expect_true(is.numeric(yr))
+
+              # confirm that results are identical between evaluation environments
+              expect_identical(yp, yr)
+
+              # confirm vectorization works as expected
+              expect_true(isVectorized(fnr))
+              expect_true(isVectorized(fnp))
+
+              par1 = rep(0.1, dimension)
+              par2 = rep(0.2, dimension)
+
+              res.seq = c(fnr(par1), fnr(par2))
+              res.vec = fnr(cbind(par1, par2))
+              expect_true(all(res.seq == res.vec),
+                          info = sprintf("Sequential and vectorized input not equal for %s (R)", getID(fnr)))
+
+              res.seq = c(fnr(par1), fnr(par2))
+              res.vec = fnr(cbind(par1, par2))
+              expect_true(all(res.seq == res.vec),
+                          info = sprintf("Sequential and vectorized input not equal for %s (Python)", getID(fnr)))
+
             }
           }
         }