Fix Issue #407 scalarize parego infill (#408)

* try to fix this mess

* finalize fix

* rename convertOptPathToDf

* remove test file

* documentation

* fix

* allow less vectors in infill crits

* get design from task data

* implement paego plot (WIP)

* Make exampleRunMultiobj plots work again

R/OptState_getter.R

@@ -34,6 +34,11 @@ getOptStateTasks = function(opt.state) {
   tasks
 }
 
+getOptStateDesigns = function(opt.state) {
+  tasks = getOptStateTasks(opt.state)
+  lapply(tasks, getTaskData)
+}
+
 getOptStateTimeModel = function(opt.state) {
   opt.path = getOptStateOptPath(opt.state)
   time.model = opt.state$time.model

R/SMBO.R

@@ -69,7 +69,7 @@ updateSMBO = function(opt.state, x, y) {
   assertNumeric(y[[1]], len = control$n.objectives)
 
 
-  infill.values = control$infill.crit$fun(points = x, models = getOptStateModels(opt.state)[[1]], control = control, design = convertOptPathToDf(opt.path, control), attributes = TRUE, iter = getOptStateLoop(opt.state))
+  infill.values = control$infill.crit$fun(points = x, models = getOptStateModels(opt.state)[[1]], control = control, designs = getOptStateDesigns(opt.state), attributes = TRUE, iter = getOptStateLoop(opt.state))
 
   prop = makeProposal(
     control = control,

R/convertOptPathToDf.R → R/convertToDesign.R

@@ -1,13 +1,17 @@
 # Takes x.vals and y.vals from the opt.path, possibly imputes missing x.vals and
 # returns a data.frame.
 #
+# @param opt.state [OptState]
 # @param opt.path [\code{\link[ParamHelpers]{optPath}}]\cr
 #   Optimization path.
+#   Might be altered for multipoint proposals etc.
 # @param control [\code{\link{MBOControl}}]\cr
 #   MBO control object.
+#   Might be altered for multipoint proposals etc.
 # @return [\code{data.frame}]
-convertOptPathToDf = function(opt.path, control) {
+#   With x and y columns
+
+convertToDesign = function(opt.path, control) {
   df = as.data.frame(opt.path, include.rest = FALSE)
-  df = convertDataFrameCols(df, ints.as.num = TRUE, logicals.as.factor = TRUE)
-  return(df)
-}
+  convertDataFrameCols(df, ints.as.num = TRUE, logicals.as.factor = TRUE)
+}

R/infillOptCMAES.R

@@ -5,15 +5,15 @@
 # @param infill.crit [\code{function}]\cr
 #   Infill criterion function.
 # @param models [\code{\link{WrappedModel}}]\cr
-#   Model fitted on design.
+#   Model fitted on designs.
 # @param control [\code{\link{MBOControl}}]\cr
 #   Control object.
 # @param par.set [\code{\link[ParamHelpers]{ParamSet}}]\cr
 #   Parameter set.
 # @param opt.path [\code{\link[ParamHelpers{OptPath}}]\cr
 #   Optimization path / archive.
-# @param design [\code{data.frame}]\cr
-#   Design of already visited points.
+# @param designs [list of \code{data.frame}]\cr
+#   List with usually one element containing the Design of already visited points.
 # @param iter [\integer(1)]
 #   Current iteration
 # @param ... [\code{ANY}]\cr
@@ -24,15 +24,15 @@
 
 # the first start is always at the best point of the current opt.path.
 # works only for numerics and integers, latter are simply rounded.
-infillOptCMAES = function(infill.crit, models, control, par.set, opt.path, design, iter, ...) {
+infillOptCMAES = function(infill.crit, models, control, par.set, opt.path, designs, iter, ...) {
   rep.pids = getParamIds(par.set, repeated = TRUE, with.nr = TRUE)
   cmaes.control = control$infill.opt.cmaes.control
 
   fn = smoof::makeSingleObjectiveFunction(
     fn = function(x) {
       newdata = as.data.frame(t(x))
       colnames(newdata) = rep.pids
-      infill.crit(newdata, models, control, par.set, design, iter, ...)
+      infill.crit(newdata, models, control, par.set, designs, iter, ...)
     },
     par.set = par.set,
     vectorized = TRUE

R/infillOptEA.R

@@ -2,7 +2,7 @@
 # kind of mimics our multi-objective approach, so we can
 # compare more honestly.
 # See infillOptCMAES.R for interface explanation.
-infillOptEA = function(infill.crit, models, control, par.set, opt.path, design, iter, ...) {
+infillOptEA = function(infill.crit, models, control, par.set, opt.path, designs, iter, ...) {
   requirePackages("emoa", why = "infillOptEA")
 
   # get constants and init shit
@@ -23,7 +23,7 @@ infillOptEA = function(infill.crit, models, control, par.set, opt.path, design,
 
     # random inital population:
     X = generateDesign(mu, par.set, fun = randomLHS)
-    y = infill.crit(X, models, control, par.set, design, iter, ...)
+    y = infill.crit(X, models, control, par.set, designs, iter, ...)
 
     for (i in seq_len(control$infill.opt.ea.maxit)) {
       # Create new individual (mu + lambda)
@@ -44,7 +44,7 @@ infillOptEA = function(infill.crit, models, control, par.set, opt.path, design,
         # Add new individual:
         X[nrow(X) + 1, ] = child1
         child2 = setColNames(as.data.frame(as.list(child1)), repids)
-        y[length(y) + 1] = infill.crit(child2, models, control, par.set, design, iter, ...)
+        y[length(y) + 1] = infill.crit(child2, models, control, par.set, designs, iter, ...)
       }
 
       # get elements we want to remove from current population

R/infillOptFocus.R

@@ -6,7 +6,7 @@
 #FIXME should we shrink if a local value is NA (dependent param)
 #
 # See infillOptCMAES.R for interface explanation.
-infillOptFocus = function(infill.crit, models, control, par.set, opt.path, ...) {
+infillOptFocus = function(infill.crit, models, control, par.set, opt.path, designs, iter, ...) {
   global.y = Inf
 
   # restart the whole crap some times
@@ -22,7 +22,7 @@ infillOptFocus = function(infill.crit, models, control, par.set, opt.path, ...)
       # convert to param encoding our model was trained on and can use
       newdesign = convertDataFrameCols(newdesign, ints.as.num = TRUE, logicals.as.factor = TRUE)
 
-      y = infill.crit(newdesign, models, control, ps.local, ...)
+      y = infill.crit(newdesign, models, control, ps.local, designs, iter, ...)
 
 
       # get current best value

R/infillOptNSGA2.R

@@ -2,7 +2,7 @@
 # NSGA2 is used to create set of candidate soluation
 # Greedy cypervolume contribution selection is used to select prop.point from
 # from candidate points
-infillOptMultiObjNSGA2 = function(infill.crit, models, control, par.set, opt.path, design, iter, ...) {
+infillOptMultiObjNSGA2 = function(infill.crit, models, control, par.set, opt.path, designs, iter, ...) {
 
   # build target function for vectorized nsga2 and run it
   rep.pids = getParamIds(par.set, repeated = TRUE, with.nr = TRUE)
@@ -13,8 +13,9 @@ infillOptMultiObjNSGA2 = function(infill.crit, models, control, par.set, opt.pat
     asMatrixRows(lapply(seq_len(m), function(i) {
       # we need to make sure mininimize in control is a scalar, so we can multiply it in infill crits...
       control$minimize = control$minimize[i]
+      control$y.name = control$y.name[i]
       infill.crit(points = newdata, models = models[i], control = control,
-        par.set = par.set, design = design, iter = iter, ...)
+        par.set = par.set, designs = designs[i], iter = iter, ...)
     }))
   }
   res = mco::nsga2(fun.tmp, idim = getParamNr(par.set, devectorize = TRUE), odim = control$n.objectives,
@@ -34,17 +35,18 @@ infillOptMultiObjNSGA2 = function(infill.crit, models, control, par.set, opt.pat
   candidate.vals = asMatrixCols(lapply(seq_len(m), function(i) {
     # we need to make sure mininimize in control is a scalar, so we can multiply it in infill crits...
     control$minimize = control$minimize[i]
-
+    control$y.name = control$y.name[i]
     newdata = as.data.frame(res$par)
     colnames(newdata) = rep.pids
     hv.contr.crit(points = newdata, models = models[i], control = control,
-      par.set = par.set, design = design, iter = iter, ...)
+      par.set = par.set, designs = designs[i], iter = iter, ...)
   }))
 
   prop.points = matrix(nrow = 0, ncol = ncol(candidate.points))
   prop.vals = matrix(nrow = 0, ncol = ncol(candidate.vals))
   colnames(prop.vals) = control$y.name
-  ys = design[, control$y.name]
+  ys = Map(function(i, y.name) designs[[i]][,y.name], i = seq_along(control$y.name), y.name = control$y.name)
+  ys = do.call(cbind, ys)
 
   ref.point = getMultiObjRefPoint(ys, control)
   prop.hv.contrs = numeric(control$propose.points)
@@ -62,6 +64,6 @@ infillOptMultiObjNSGA2 = function(infill.crit, models, control, par.set, opt.pat
 
   # FIXME: cleanup - i'm reall unsure how to set the names of prop.points technically
   prop.points = as.data.frame(prop.points)
-  colnames(prop.points) = names(design[, colnames(design) %nin% control$y.name])
+  colnames(prop.points) = getParamIds(par.set, repeated = TRUE, with.nr = TRUE)
   list(prop.points = prop.points, prop.hv.contrs = prop.hv.contrs)
 }

R/infill_crits.R

@@ -54,7 +54,7 @@ NULL
 #' @rdname infillcrits
 makeMBOInfillCritMeanResponse = function() {
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       ifelse(control$minimize, 1, -1) * predict(models[[1L]], newdata = points)$data$response
     },
     name = "Mean response",
@@ -67,7 +67,7 @@ makeMBOInfillCritMeanResponse = function() {
 #' @rdname infillcrits
 makeMBOInfillCritStandardError = function() {
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
        -predict(models[[1L]], newdata = points)$data$se
     },
     name = "Standard error",
@@ -83,10 +83,13 @@ makeMBOInfillCritEI = function(se.threshold = 1e-6) {
   assertNumber(se.threshold, lower = 1e-20)
   force(se.threshold)
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       model = models[[1L]]
-      maximize.mult = ifelse(control$minimize, 1, -1)
+      design = designs[[1]]
+      maximize.mult = if (control$minimize) 1 else -1
+      assertString(control$y.name)
       y = maximize.mult * design[, control$y.name]
+      assertNumeric(y, any.missing = FALSE)
       p = predict(model, newdata = points)$data
       p.mu = maximize.mult * p$response
       p.se = p$se
@@ -117,9 +120,9 @@ makeMBOInfillCritCB = function(cb.lambda = NULL) {
   assertNumber(cb.lambda, lower = 0, null.ok = TRUE)
   force(cb.lambda)
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       model = models[[1L]]
-      maximize.mult = ifelse(control$minimize, 1, -1)
+      maximize.mult = if (control$minimize) 1 else -1
       p = predict(model, newdata = points)$data
       #FIXME: removed cb.inflate.se for now (see issue #309)
       # if (cb.inflate.se) {
@@ -159,9 +162,10 @@ makeMBOInfillCritAEI = function(aei.use.nugget = FALSE, se.threshold = 1e-6) {
   force(se.threshold)
 
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       model = models[[1L]]
-      maximize.mult = ifelse(control$minimize, 1, -1)
+      design = designs[[1L]]
+      maximize.mult = if (control$minimize) 1 else -1
       p = predict(model, newdata = points)$data
       p.mu = maximize.mult * p$response
       p.se = p$se
@@ -205,9 +209,10 @@ makeMBOInfillCritEQI = function(eqi.beta = 0.75, se.threshold = 1e-6) {
   force(se.threshold)
 
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       model = models[[1L]]
-      maximize.mult = ifelse(control$minimize, 1, -1)
+      design = designs[[1L]]
+      maximize.mult = if (control$minimize) 1 else -1
       # compute q.min
       design_x = design[, (colnames(design) %nin% control$y.name), drop = FALSE]
       p.current.model = predict(object = model, newdata = design_x)$data
@@ -258,10 +263,11 @@ makeMBOInfillCritDIB = function(cb.lambda = 1, sms.eps = NULL) {
   if (!is.null(sms.eps))
     assertNumber(sms.eps, lower = 0, finite = TRUE)
   makeMBOInfillCrit(
-    fun = function(points, models, control, par.set, design, iter, progress, attributes = FALSE) {
+    fun = function(points, models, control, par.set, designs, iter, progress, attributes = FALSE) {
       # get ys and cb-value-matrix for new points, minimize version
       maximize.mult = ifelse(control$minimize, 1, -1)
-      ys = as.matrix(design[, control$y.name]) %*% diag(maximize.mult)
+      ys = Map(function(i, y.name) designs[[i]][, y.name], i = seq_along(control$y.name), y.name = control$y.name)
+      ys = do.call(cbind, ys) %*% diag(maximize.mult)
 
       ps = lapply(models, predict, newdata = points)
       means = extractSubList(ps, c("data", "response"), simplify = "cols")

R/makeMBOInfillCrit.R

@@ -52,7 +52,7 @@ makeMBOInfillCrit = function(fun, name, id,
   requires.se = FALSE) {
   assertFunction(
     fun,
-    args = c("points", "models", "control", "par.set", "design", "iter", "progress", "attributes"),
+    args = c("points", "models", "control", "par.set", "designs", "iter", "progress", "attributes"),
     ordered = TRUE)
 
   assertString(name)

R/makeTaskSingleObj.R

@@ -10,7 +10,7 @@
 #   MBO control object.
 # @return [\code{\link[mlr]{SupervisedTask}}]
 makeTaskSingleObj = function(opt.path, control) {
-  data = convertOptPathToDf(opt.path, control)
+  data = convertToDesign(opt.path, control)
   data$dob = data$eol = NULL
 
   # user selected to (log)-transform the y-column

R/makeTasks.R

@@ -8,7 +8,7 @@ makeTasks = function(opt.state) {
     tasks = list(makeTaskSingleObj(opt.path, control))
   } else {
     if (control$multiobj.method == "parego")
-      tasks = makeTasksParEGO(opt.path, control, all.possible.weights = getOptProblemAllPossibleWeights(opt.problem))
+      tasks = makeTasksParEGO(opt.state)
     else
       tasks = makeTasksMultiObj(opt.path, control)
   }

R/makeTasksMultiObj.R

@@ -10,7 +10,7 @@
 #   MBO control object.
 # @return [\code{list(\link[mlr]{SupervisedTask}})]
 makeTasksMultiObj = function(opt.path, control) {
-  data = convertOptPathToDf(opt.path, control)
+  data = convertToDesign(opt.path, control)
 
   # FIXME: trafo.y.fun should be a list of length y.name
   # user selected to (log)-transform the y-column

R/makeTasksParEGO.R

@@ -12,10 +12,28 @@
 # @return [\code{list}] List with elements
 #   tasks: list of tasks
 #   weights: matrix of used weight vectors
-makeTasksParEGO = function(opt.path, control, all.possible.weights) {
+makeTasksParEGO = function(opt.state) {
+
+  res = generateParEgoDfData(opt.state)
+  # create the scalarized regression tasks
+  tasks = apply(res$lambdas, 1, function(lambda) {
+    data = generateParEgoDf(res, lambda)
+    makeRegrTask(target = "y.scalar", data = data) 
+  })
+  attr(tasks, "weight.mat") = res$lambdas
+  return(tasks)
+}
+
+generateParEgoDfData = function(opt.state) {
+
+  opt.path = getOptStateOptPath(opt.state)
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  all.possible.weights = getOptProblemAllPossibleWeights(opt.problem)
+
   n.points = control$propose.points
   # get data + normalize the targets to [0, 1] + drop them from data
-  data = convertOptPathToDf(opt.path, control)
+  data = convertToDesign(opt.path, control)
   data = dropNamed(data, control$y.name)
   y = getOptPathY(opt.path)
   if (control$multiobj.parego.normalize == "standard") {
@@ -58,15 +76,18 @@ makeTasksParEGO = function(opt.path, control, all.possible.weights) {
   }
   lambdas = rbind(margin.points, lambdas)
 
-  # create the scalarized regression tasks
-  tasks = vector(mode = "list", length = n.points)
-  for (loop in seq_len(n.points)) {
-    # make sure to minimize, then create scalarized response
-    lambda = lambdas[loop,, drop = TRUE] * ifelse(control$minimize, 1, -1)
-    y2 = y %*% diag(lambda)
-    data$y.scalar = apply(y2, 1, max) + control$multiobj.parego.rho * rowSums(y2)
-    tasks[[loop]] = makeRegrTask(target = "y.scalar", data = data)
-  }
-  attr(tasks, "weight.mat") = lambdas
-  return(tasks)
+  list(
+    lambdas = lambdas,
+    minimize = control$minimize,
+    multiobj.parego.rho = control$multiobj.parego.rho,
+    y = y,
+    data = data
+  )
 }
+
+generateParEgoDf = function(res, lambda) {
+  lambda = lambda * ifelse(res$minimize, 1, -1)
+  y2 = res$y %*% diag(lambda)
+  res$data$y.scalar = apply(y2, 1, max) + res$multiobj.parego.rho * rowSums(y2)
+  return(res$data)
+}

R/multiobj_helpers.R

@@ -50,10 +50,11 @@ getMultiObjRefPoint = function (ys, control, minimize = control$minimize) {
 
 # evaluate an infill crit on multiple models (one per objective in multi-objective)
 # returns matrix of crit vals, rows = points, cols = crits
-evalCritFunForMultiObjModels = function(infill.crit.fun, points, models, control, par.set, design, iter) {
+evalCritFunForMultiObjModels = function(infill.crit.fun, points, models, control, par.set, designs, iter) {
   asMatrixCols(lapply(seq_along(models), function(i) {
       # we need to make sure mininimize in control is a scalar, so we can multiply it in infill crits...
       control$minimize = control$minimize[i]
-      infill.crit.fun(points, models[i], control, par.set, design, iter)
+      control$y.name = control$y.name[i]
+      infill.crit.fun(points, models[i], control, par.set, designs[i], iter)
   }))
 }