## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5)

## ----setup--------------------------------------------------------------------
library(paretoinfer)

## ----problem------------------------------------------------------------------
means <- rbind(A1 = c(0.15, 0.70), A2 = c(0.40, 0.40), A3 = c(0.70, 0.15),
               A4 = c(0.42, 0.43), A5 = c(0.60, 0.60), A6 = c(0.30, 0.80))
simulate <- function(alternative, n, h = 0.1) {
  mu <- means[alternative, ]
  data.frame(y1 = runif(n, mu[1] - h, mu[1] + h), y2 = runif(n, mu[2] - h, mu[2] + h))
}
design <- pareto_design(alternatives = rownames(means), objectives = c("y1", "y2"),
                        bounds = c(0, 1), alpha = 0.05, epsilon = 0.05, max_evaluations = 3000)

## ----adaptive-----------------------------------------------------------------
set.seed(1)
adaptive <- pareto_run(design, simulate, batch = 5, initial = 5)
glance(adaptive)[, c("n_evaluations", "n_certified_optimal", "n_dominated", "identified")]

## ----uniform------------------------------------------------------------------
set.seed(1)
uniform <- initialize_pareto(design)
while (is.na(uniform$stopping_reason)) {
  for (k in design$alternatives) {
    ev <- simulate(k, 5); ev$alternative <- k
    uniform <- update_objectives(uniform, ev)
    if (!is.na(uniform$stopping_reason)) break
  }
}
glance(uniform)[, c("n_evaluations", "n_certified_optimal", "n_dominated", "identified")]

## ----per-alternative----------------------------------------------------------
data.frame(alternative = design$alternatives,
           adaptive = unname(adaptive$n), uniform = unname(uniform$n))

## ----guarantee----------------------------------------------------------------
frontier <- rbind(A1 = c(0.30, 0.60), A2 = c(0.50, 0.40), A3 = c(0.45, 0.45))
sim_frontier <- function(alternative, n) {
  mu <- frontier[alternative, ]
  data.frame(y1 = runif(n, mu[1] - 0.2, mu[1] + 0.2), y2 = runif(n, mu[2] - 0.2, mu[2] + 0.2))
}
design0 <- pareto_design(alternatives = rownames(frontier), objectives = c("y1", "y2"),
                         bounds = c(0, 1), alpha = 0.05, epsilon = 0, max_evaluations = 300)
set.seed(2)
discards <- replicate(12, {
  state <- suppressWarnings(pareto_run(design0, sim_frontier, batch = 10, initial = 10))
  length(plausible_pareto_set(state)$dominated)
})
table(discards)

## ----epsilon------------------------------------------------------------------
set.seed(3)
by_epsilon <- sapply(c(0.02, 0.05, 0.10), function(eps) {
  d <- pareto_design(alternatives = rownames(means), objectives = c("y1", "y2"),
                     bounds = c(0, 1), alpha = 0.05, epsilon = eps, max_evaluations = 4000)
  s <- pareto_run(d, simulate, batch = 5, initial = 5)
  c(evaluations = s$n_evaluations, dominated = length(plausible_pareto_set(s)$dominated))
})
colnames(by_epsilon) <- paste0("eps = ", c(0.02, 0.05, 0.10))
by_epsilon

## ----alpha--------------------------------------------------------------------
set.seed(4)
by_alpha <- sapply(c(0.10, 0.05, 0.01), function(a) {
  d <- pareto_design(alternatives = rownames(means), objectives = c("y1", "y2"),
                     bounds = c(0, 1), alpha = a, epsilon = 0.05, max_evaluations = 4000)
  pareto_run(d, simulate, batch = 5, initial = 5)$n_evaluations
})
data.frame(alpha = c(0.10, 0.05, 0.01), evaluations = by_alpha)