Noisy Optimization: OCBA and Incumbent

NAMESPACE

@@ -58,6 +58,7 @@ export(renderExampleRunPlot)
 export(setMBOControlInfill)
 export(setMBOControlMultiObj)
 export(setMBOControlMultiPoint)
+export(setMBOControlNoisy)
 export(setMBOControlTermination)
 export(trafoLog)
 export(trafoSqrt)

R/OptState_getter.R

@@ -144,3 +144,8 @@ getOptStateValidTerminationStates = function() {
   c("term.iter", "term.time", "term.exectime", "term.yval", "term.feval", "term.custom")
 }
 
+getOptStateIntensification = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  return(!is.null(control$noisy.method) && control$noisy.method %in% c("incumbent", "ocba"))
+}

R/evalTargetFun.R

@@ -25,13 +25,30 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
   # short names and so on
   nevals = length(xs)
   ny = control$n.objectives
+
+  # trafo X points
+  xs.trafo = lapply(xs, trafoValue, par = par.set)
+
+  # handle noisy instances
+  if (isTRUE(control$noisy.instances > 1L)) {
+    nevals = nevals * control$noisy.instances
+    xs = rep(xs, each = control$noisy.instances)
+    extras = rep(extras, each = control$noisy.instances)
+    if (!control$noisy.self.replicating) {
+      xs.trafo = rep(xs.trafo, each = control$noisy.instances)
+      if (!is.na(control$noisy.instance.param)) {
+        inst.param = lapply(seq_len(control$noisy.instances), function(x) setNames(list(x), control$noisy.instance.param))
+        xs.trafo = Map(c, xs.trafo, inst.param)
+      }
+    }
+  }
+
+
   num.format = control$output.num.format
   num.format.string = paste("%s = ", num.format, sep = "")
   dobs = ensureVector(asInteger(getOptStateLoop(opt.state)), n = nevals, cl = "integer")
   imputeY = control$impute.y.fun
 
-  # trafo X points
-  xs.trafo = lapply(xs, trafoValue, par = par.set)
 
   # function to measure of fun call
     wrapFun = function(x) {
@@ -43,6 +60,9 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
         user.extras = attr(y, "extras")
         y = setAttribute(y, "extras", NULL)
       }
+      if (!is.null(control$noisy.instance.param) && !is.na(control$noisy.instance.param) && !control$noisy.self.replicating) {
+        user.extras = c(user.extras, x[control$noisy.instance.param])
+      }
       st = proc.time() - st
       list(y = y, time = st[3], user.extras = user.extras)
     }
@@ -56,6 +76,21 @@ evalTargetFun.OptState = function(opt.state, xs, extras) {
   res = parallelMap(wrapFun, xs.trafo, level = "mlrMBO.feval",
     impute.error = if (is.null(imputeY)) NULL else identity)
 
+  # handle noisy instances of self.replicating functions
+  if (isTRUE(control$noisy.instances > 1L) && control$noisy.self.replicating) {
+    xs.trafo = rep(xs.trafo, each = control$noisy.instances)
+    res = lapply(res, function(r) {
+      if (is.error(r)) {
+        rep(list(r), control$noisy.instances)
+      } else {
+        lapply(seq_along(r$y), function(i) {
+          list(y = r$y[i], time = r$time / length(r$y), user.extras = c(r$user.extras, setNames(list(i), control$noisy.instance.param)))
+        })
+      }
+    })
+    res = unlist(res, recursive = FALSE)
+  }
+
   # loop evals and to some post-processing
   for (i in seq_len(nevals)) {
     r = res[[i]]; x = xs[[i]]; x.trafo = xs.trafo[[i]]; dob = dobs[i]

R/intensifyOptState.R

@@ -0,0 +1,170 @@
+intensifyOptState = function(opt.state) {
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+
+  switch(control$noisy.method, 
+      "incumbent" = intensifyIncumbent(opt.state),
+      "ocba" = intensifyOCBA(opt.state)
+    )
+}
+
+intensifyIncumbent = function(opt.state) {
+
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+
+  # get a summary of the design  
+  ds = getOptPathSummary(opt.state, par.names)
+  nds = nrow(ds)
+
+  # incumbent: current best point w. r. t. mean over all function evaluations
+  # the newest point cannot be the incumbent, it is always a challenger
+  # DOES NOT WORK FOR MULTIPOINT PROPOSAL YET
+  inc = which.min(ds[- nds, ]$y) 
+  # incumbent is replicated once in each iteration
+  replicatePoint(opt.state, x = ds[inc, ..par.names], type = "incumbent", reps = 1L)
+
+  # determine a set of challengers
+  if (control$noisy.incumbent.nchallengers == 0L) {
+    cls = c(nds)
+  } else {
+    # determine set of p points to be challenged against incumbent
+    # incumbent is excluded (cannot be challenged against itself)
+    # and new point is always set as a challenger
+    # points are drawn randomly without replacement with probability prop. to their function value
+    cls = setdiff(seq_len(nds), c(inc, nds))
+    p = min(control$noisy.incumbent.nchallengers, nds - 2)
+    probs = exp(- ds[cls, ]$y) / sum(exp(- ds[cls, ]$y))
+    cls = sample(cls, size = p, prob = probs, replace = FALSE)
+    cls = c(cls, nds)
+  }
+
+  # start the race
+  for (cl in cls) {
+
+    r = 1L
+    replicatePoint(opt.state, x = ds[cl, ..par.names], type = paste("challenger"), reps = r)
+    ds = getOptPathSummary(opt.state, par.names)
+
+    # proceed as long as challenger has less runs than incumbent and is better than incumbent
+    while((ds[cl, "runs"] < ds[inc, "runs"]) && (ds[cl, "y"] < ds[inc, "y"])) {
+      r = 2L * r
+      replicatePoint(opt.state, x = ds[cl, ..par.names], type = paste("challenger"), reps = r)
+      ds = getOptPathSummary(opt.state, par.names)
+    }
+
+  }
+  return(opt.state)
+}
+
+replicatePoint = function(opt.state, x, type, reps = 1L) {
+
+  # replicate rows according to the number of desired replicates
+  xs = seq_len(nrow(x))
+  xrep = x[rep(xs, reps), ]
+
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+
+  prop = makeProposal(control, xrep, prop.type = rep(type, nrow(xrep)))
+  evalProposedPoints.OptState(opt.state, prop)
+
+  return(opt.state)
+}
+
+getOptPathSummary = function(opt.state, par.names) {
+    op = as.data.table(getOptStateOptPath(opt.state))
+    ds = op[, .(y = mean(y), ysd = sd(y), runs = .N), by = par.names]
+    return(ds)
+}
+
+
+intensifyOCBA = function(opt.state) {
+
+  # some intialization
+  opt.problem = getOptStateOptProblem(opt.state)
+  control = getOptProblemControl(opt.problem)
+  par.set = getOptProblemParSet(opt.problem)
+
+  op = as.data.table(getOptStateOptPath(opt.state))
+  par.names = colnames(op)[1:(which(colnames(op) == "y") - 1)] #FIXME: This sucks
+
+  # minimum number of replicates at each point 
+  minrep = max(control$noisy.ocba.initial, 2L)
+
+  # calculate summary of the dsign 
+  ds = getOptPathSummary(opt.state, par.names)
+  nds = nrow(ds)
+
+  # make sure that initially, each point is evaluated at least minrep times 
+  xinit = rep(seq_len(nds), pmax(minrep - ds$runs, 0))
+  opt.state = replicatePoint(opt.state, x = ds[xinit, ..par.names], type = paste("initeval"))
+
+  ds = getOptPathSummary(opt.state, par.names)
+  add = distributeOCBA(ds, budget = control$noisy.ocba.budget)
+  reps = rep(seq_len(nds), add)
+
+  replicatePoint(opt.state, x = ds[reps, ..par.names], type = paste("OCBA"))
+
+  return(opt.state)
+}
+
+
+distributeOCBA = function(ds, budget) {
+
+  nds = nrow(ds)
+
+  # TODO: until now only minimization possible 
+  tbudget = budget + sum(ds$runs)
+
+  # search for the best and second-best dsign
+  b = order(ds$y)[1]
+  s = order(ds$y)[2]
+
+  # vector of ratios
+  ratio = rep(0, nds)
+  ratio[s] = 1
+
+  # calculate ratios
+  tmp = (ds[b, ]$y - ds[s, ]$y) / (ds[b, ]$y - ds[- c(s, b), ]$y)
+  ratio[- c(s, b)] = tmp^2 * ds[- c(s, b), ]$ysd^2 / ds[s, ]$ysd^2
+  ratio[b] = ds[b, ]$ysd * sqrt(sum(ratio^2 / ds$ysd^2))
+
+  # additional replications
+  add = rep(0, nds)
+
+  # do not disable any dsign
+  disabled = rep(FALSE, nds)
+
+  more_alloc = TRUE
+
+  while (more_alloc) {
+
+    add[!disabled] = roundPreserveSum(tbudget / sum(ratio[!disabled]) * ratio[!disabled])
+
+    # disable designs that have been run too much
+    disabled = disabled | (ds$runs > add)
+    more_alloc = any(ds$runs > add)
+
+    # set additional replications s.t. already run replications are set
+    add[disabled] = ds[disabled, ]$runs
+
+    # decrease total budget correspondingly
+    tbudget = budget + sum(ds$runs) - sum(add[disabled])
+  }
+
+  add = add - ds$runs
+
+  return(add)
+}
+
+roundPreserveSum = function(x, digits = 0) {
+  up = 10 ^ digits
+  x = x * up
+  y = floor(x)
+  indices = tail(order(x-y), round(sum(x)) - sum(y))
+  y[indices] = y[indices] + 1
+  y / up
+}

R/makeMBOControl.R

@@ -187,11 +187,21 @@ print.MBOControl = function(x, ...) {
   } else {
     if (x$propose.points == 1) {
       print(x$infill.crit)
+      catf("")
       catf("Infill optimizer              : %s", x$infill.opt)
       catf("Infill optimizer restarts     : %i", x$infill.opt.restarts)
     } else {
       catf("Multi-point method            : %s", x$multipoint.method)
     }
     catf("Final point by                : %s", x$final.method)
   }
+
+  # if (!is.null(x$noisy.method))
+  #   if (x$noisy.method == "fixed")
+  #     catf("Replication strategy          : fixed (noisy instances = %i)", x$noisy.instances)
+  #   if (x$noisy.method == "ocba")
+  #     catf("Replication strategy          : ocba (initial budget = %i, replication budget = %i)", 
+  #       x$noisy.ocba.budget, x$noisy.ocba.initial)
+  #   if(x$noisy.method == "incumbent")
+  #     catf("Replication strategy          : incumbent (number of challengers = (1 + %i)", x$noisy.incumbent.nchallengers)
 }

R/makeTaskSingleObj.R

@@ -21,5 +21,13 @@ makeTaskSingleObj = function(opt.path, control) {
     data[[y.name]] = trafo.y.fun(data[[y.name]])
   }
 
+  agg = control$noisy.instance.aggregation
+
+  if (!is.null(agg)) {
+    par.names = colnames(data)[1:(which(colnames(data) == "y") - 1)] #FIXME: This sucks
+    data = setDT(data)[, .(y = agg(y)), by = par.names]
+    data = data.frame(data)
+  }
+
   makeRegrTask(target = control$y.name, data = data)
 }

R/mboTemplate.R

@@ -37,6 +37,10 @@ mboTemplate.OptState = function(obj) {
   repeat {
     prop = proposePoints(opt.state)
     evalProposedPoints.OptState(opt.state, prop)
+    intensify = getOptStateIntensification(opt.state)
+    if (intensify) {
+      intensifyOptState(opt.state)
+    }
     finalizeMboLoop(opt.state)
     terminate = getOptStateTermination(opt.state)
     if (terminate$term) {

R/setMBOControlNoisy.R

@@ -0,0 +1,57 @@
+#' @title Set options for handling noisy functions.
+#' @description
+#' Extends an MBO control object with options for handling noisy functions. 
+#' @template arg_control
+#' @param method [\code{character(1)}]\cr
+#'   Which of the replication strategies should be used? Possible values are:
+#'   \dQuote{fixed}: Every point is evaluated \code{instances} times. \cr
+#'   \dQuote{incumbent}: Use an incumbent strategy as intensification strategy. 
+#'    The size of the set of additional challengers (apart from the incumbent) can be specified in \code{incumbent.nchallengers}. \cr
+#'   \dQuote{ocba}: Distribution replication budget according to OCBA. 
+#'    The replication budget per iteration is specified in \code{ocba.budget},
+#'    the initial number of iterations per parameter is specified in \code{ocba.initial}. \cr 
+#' @param instances [\code{integer(1)}]\cr
+#'   How many instances of one parameter will be calculated?
+#' @param instance.param [\code{character(1)}]\cr
+#'   What is the name of the function param that defines the instance?
+#' @param self.replicating [\code{logical(1)}]\cr
+#'   TRUE if the function returns a vector of noisy results for one input. Then \code{instances} specifies the length of the result we expect.
+#' @param incumbent.nchallengers [\code{integer(1)}]\cr
+#'   The size of the set of additional challengers (apart from the incumbent), defaults to \code{0}.
+#' @param ocba.budget [\code{integer(1)}]\cr
+#'   The budget that is allocated in each iteration per Optimal Computing Budget Allocation (OCBA) rule, defaults to 10.
+#' @param ocba.initial [\code{integer(1)}]\cr
+#'   The number of initial replications at each new point, defaults to \code{3}. 
+#'   This needs to be larger than 1, since OCBA requires an initial variance estimate at each point.
+#' @param instance.aggregation [\code{function}]\cr
+#'   Should data be aggregated per instance? If yes, a function (e. g. mean) needs to be specified.  
+#' @return [\code{\link{MBOControl}}].
+#' @family MBOControl
+#' @export
+setMBOControlNoisy = function(control,
+  method = NULL,
+  instances = NULL,
+  instance.param = NULL,
+  instance.aggregation = NULL,
+  self.replicating = NULL,
+  incumbent.nchallengers = NULL, 
+  ocba.budget = NULL,
+  ocba.initial = NULL) {
+
+  assertClass(control, "MBOControl")
+  control$noisy.method = coalesce(method, control$noisy.method, "fixed")
+  assertChoice(control$noisy.method, choices = c("fixed", "incumbent", "ocba"))
+  control$noisy.instances = assertInt(instances, lower = 1L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.instances %??% 1L
+  control$noisy.self.replicating = assertFlag(self.replicating, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.self.replicating %??% FALSE
+  control$noisy.instance.param = assertString(instance.param, null.ok = TRUE, na.ok = TRUE) %??% control$noisy.instance.param %??% ifelse(control$noisy.self.replicating, "noisy.repl", NA_character_)
+  control$noisy.instance.aggregation = assertClass(instance.aggregation, "function", null.ok = TRUE) %??% control$noisy.instance.aggregation 
+  control$noisy.ocba.budget = assertInt(ocba.budget, lower = 1L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.ocba.budget %??% 10L
+  control$noisy.ocba.initial = assertInt(ocba.initial, lower = 2L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.ocba.initial %??% 3L
+  control$noisy.incumbent.nchallengers = assertInt(incumbent.nchallengers, lower = 0L, null.ok = TRUE, na.ok = FALSE) %??% control$noisy.incumbent.nchallengers %??% 0L
+
+  if (control$noisy.self.replicating && control$noisy.instance.param != "noisy.repl") {
+    stop("You can not change the instance.param for self replicating functions.")
+  }
+
+  return(control)
+}