goldstein_price <- function(x) {
x1 <- x[1]
x2 <- x[2]
u1 <- (x1 + x2 + 1)^2
u2 <- 19 - 14 * x1 + 3 * x1^2 - 14 * x2 + 6 * x1 * x2 + 3 * x2^2
u3 <- (2 * x1 - 3 * x2)^2
u4 <- 18 - 32 * x1 + 12 * x1^2 + 48 * x2 - 36 * x1 * x2 +
27 * x2^2
(1 + u1 * u2) * (30 + u3 * u4)
}
rosenbrock <- function(x) {
100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
}
six_hump_camelback <- function(x) {
x1 <- x[1]
x2 <- x[2]
(4 - 2.1 * x1^2 + x1^4 / 3) * x1^2 + x1 * x2 +
(-4 + 4 * x2^2) * x2^2
}
rastrigin_duan <- function(x) {
x[1]^2 + x[2]^2 - cos(18 * x[1]) - cos(18 * x[2])
}
griewank_duan <- function(x) {
divisor <- if (length(x) == 2L) 200 else 4000
sum(x^2 / divisor) - prod(cos(x / sqrt(seq_along(x)))) + 1
}
shekel <- function(x) {
a <- rbind(
c(4, 1, 8, 6, 3, 2, 5, 8, 6, 7),
c(4, 1, 8, 6, 7, 9, 5, 1, 2, 3.6),
c(4, 1, 8, 6, 3, 2, 3, 8, 6, 7),
c(4, 1, 8, 6, 7, 9, 3, 1, 2, 3.6)
)
c_values <- c(0.1, 0.2, 0.2, 0.4, 0.4, 0.6, 0.3, 0.7, 0.5, 0.5)
-sum(vapply(seq_along(c_values), function(i) {
1 / (sum((x - a[, i])^2) + c_values[i])
}, numeric(1)))
}
hartman <- function(x) {
a6 <- rbind(
c(10, 0.05, 3, 17),
c(3, 10, 3.5, 8),
c(17, 17, 1.7, 0.05),
c(3.5, 0.1, 10, 10),
c(1.7, 8, 17, 0.1),
c(8, 14, 8, 14)
)
c6 <- c(1, 1.2, 3, 3.2)
p6 <- rbind(
c(0.1312, 0.2329, 0.2348, 0.4047),
c(0.1696, 0.4135, 0.1451, 0.8828),
c(0.5569, 0.8307, 0.3522, 0.8732),
c(0.0124, 0.3736, 0.2883, 0.5743),
c(0.8283, 0.1004, 0.3047, 0.1091),
c(0.5886, 0.9991, 0.6650, 0.0381)
)
-sum(vapply(seq_along(c6), function(i) {
c6[i] * exp(-sum(a6[, i] * (x - p6[, i])^2))
}, numeric(1)))
}
duan_control <- list(
seed = 1L,
complexes = 5L,
max_evaluations = 10000L,
kstop = 5L,
pcento = 0,
parameter_epsilon = 1e-6
)
duan_problems <- list(
goldstein_price = list(
fn = goldstein_price,
lower = c(-2, -2),
upper = c(2, 2),
initial = c(1, 1),
optimum = 3
),
rosenbrock = list(
fn = rosenbrock,
lower = c(-5, -5),
upper = c(5, 5),
initial = c(-3, 4),
optimum = 0
),
six_hump_camelback = list(
fn = six_hump_camelback,
lower = c(-5, -2),
upper = c(5, 8),
initial = c(0, 0),
optimum = -1.031628453489877
),
rastrigin_duan = list(
fn = rastrigin_duan,
lower = c(-1, -1),
upper = c(1, 1),
initial = c(0.5, 0.5),
optimum = -2
),
griewank_duan = list(
fn = griewank_duan,
lower = rep(-600, 10),
upper = rep(600, 10),
initial = rep(100, 10),
optimum = 0
),
shekel = list(
fn = shekel,
lower = rep(0, 4),
upper = rep(10, 4),
initial = rep(1, 4),
optimum = -10.5364098252
),
hartman = list(
fn = hartman,
lower = rep(0, 6),
upper = rep(1, 6),
initial = rep(0.5, 6),
optimum = -3.322368011415515
)
)
sceua_duan_value <- function(problem) {
set.seed(duan_control$seed)
n <- length(problem$lower)
sceua::sceua(
fn = problem$fn,
lower = problem$lower,
upper = problem$upper,
initial = problem$initial,
max_evaluations = duan_control$max_evaluations,
kstop = duan_control$kstop,
pcento = duan_control$pcento,
complexes = duan_control$complexes,
points_per_complex = 2L * n + 1L,
simplex_size = n + 1L,
evolution_steps = 2L * n + 1L,
min_complexes = duan_control$complexes,
parameter_epsilon = duan_control$parameter_epsilon
)$value
}
rtop_duan_value <- function(problem) {
set.seed(duan_control$seed)
n <- length(problem$lower)
silent_rtop(
problem$fn,
pars = problem$initial,
lower = problem$lower,
upper = problem$upper,
maxn = duan_control$max_evaluations,
kstop = duan_control$kstop,
pcento = duan_control$pcento,
ngs = duan_control$complexes,
npg = 2L * n + 1L,
nps = n + 1L,
nspl = 2L * n + 1L,
mings = duan_control$complexes,
peps = duan_control$parameter_epsilon,
iniflg = 1
)
}
soilhyp_duan_value <- function(problem) {
set.seed(duan_control$seed)
n <- length(problem$lower)
silent_soilhyp(
FUN = problem$fn,
par = problem$initial,
lower = problem$lower,
upper = problem$upper,
control = list(
maxeval = duan_control$max_evaluations,
ncomplex = duan_control$complexes,
cce.iter = 2L * n + 1L,
reltol = duan_control$pcento,
tolsteps = duan_control$kstop,
initsample = "latin"
)
)
}
safe_duan_value <- function(package, problem) {
value <- tryCatch(
switch(
package,
sceua = sceua_duan_value(problem),
rtop = rtop_duan_value(problem),
soilhyp = soilhyp_duan_value(problem)
),
error = function(e) NA_real_
)
as.numeric(value)
}
# Number of timed runs per (problem, package). bench::mark runs each
# expression exactly this many times.
duan_iterations <- 5L
# Dependency-free cell visual: a fixed-width unicode bar for a scalar value
# scaled to a per-problem maximum. A positive sub-maximum value always gets
# at least one filled cell so the fastest run still shows a sliver.
bar_text <- function(value, max_value, width = 10L,
fill = "\u2588", empty = "\u2591") {
if (is.na(value) || is.na(max_value) || max_value == 0) return("")
n <- round(width * value / max_value)
if (value > 0 && value < max_value) n <- max(n, 1L)
n <- min(max(n, 0L), width)
paste0(strrep(fill, n), strrep(empty, width - n))
}
fmt_time <- function(seconds) {
if (is.na(seconds)) return("\u2014")
if (seconds >= 1) return(sprintf("%.2f s", seconds))
if (seconds >= 1e-3) return(sprintf("%.2f ms", seconds * 1e3))
sprintf("%.0f \u00b5s", seconds * 1e6)
}
# Compute the best value for every (problem, package) pair. `safe_duan_value`
# returns NA when a package errors, so a failed run still has a row.
packages <- c("sceua", "rtop", "soilhyp")
duan_values <- do.call(rbind, lapply(names(duan_problems), function(name) {
problem <- duan_problems[[name]]
values <- vapply(packages, safe_duan_value, numeric(1), problem = problem)
data.frame(
problem = name,
package = packages,
value = values,
abs_error = abs(values - problem$optimum),
stringsAsFactors = FALSE
)
}))
# Time every pair. `check = FALSE` disables bench's result-equality check; the
# tables below compare values to the optimum explicitly. A pair that errors is
# recorded with `NA` median time rather than aborting the vignette.
benchmark_duan_row <- function(name, package, problem) {
exprs <- list(
sceua = quote(sceua_duan_value(problem)),
rtop = quote(rtop_duan_value(problem)),
soilhyp = quote(soilhyp_duan_value(problem))
)
set.seed(duan_control$seed)
res <- tryCatch(
bench::mark(
eval(exprs[[package]]),
iterations = duan_iterations,
relative = FALSE,
check = FALSE
),
error = function(e) NULL
)
median_time <- if (is.null(res)) NA_real_ else as.numeric(res$median[[1]])
data.frame(
problem = name,
package = package,
median_time = median_time,
stringsAsFactors = FALSE
)
}
duan_timing <- do.call(rbind, lapply(names(duan_problems), function(name) {
problem <- duan_problems[[name]]
do.call(rbind, lapply(packages, function(pkg) {
benchmark_duan_row(name, pkg, problem)
}))
}))
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
duan_summary <- merge(
duan_values,
duan_timing,
by = c("problem", "package"),
sort = FALSE
)
duan_summary <- duan_summary[order(
match(duan_summary$problem, names(duan_problems)),
match(duan_summary$package, packages)
), ]
duan_summary$optimum <- vapply(
duan_summary$problem,
function(n) duan_problems[[n]]$optimum,
numeric(1)
)
fmt_value <- function(x) ifelse(is.na(x), "\u2014", sprintf("%.6g", x))
fmt_err <- function(x) ifelse(is.na(x), "\u2014", sprintf("%.3e", x))
# Per-problem table. The Timing bar is scaled to the slowest run within the
# problem, so a full bar marks the slowest package and shorter bars mark
# faster ones. A package that errors shows em dashes.
duan_problem_table <- function(name) {
rows <- duan_summary[duan_summary$problem == name, , drop = FALSE]
max_time <- if (all(is.na(rows$median_time))) {
NA_real_
} else {
max(rows$median_time, na.rm = TRUE)
}
timing_bar <- vapply(
rows$median_time, bar_text, character(1), max_value = max_time
)
data.frame(
Package = rows$package,
`Best value` = fmt_value(rows$value),
Optimum = sprintf("%.6g", rows$optimum[1]),
MAE = fmt_err(rows$abs_error),
`Median time` = vapply(rows$median_time, fmt_time, character(1)),
Timing = ifelse(timing_bar == "", "\u2014", timing_bar),
check.names = FALSE,
stringsAsFactors = FALSE
)
}
duan_problem_caption <- function(name) {
p <- duan_problems[[name]]
sprintf(
"%s: bounds [%s], optimum %g, %d parameters",
name,
paste(sprintf("[%g, %g]", p$lower, p$upper), collapse = ", "),
p$optimum,
length(p$lower)
)
}