Nelder–Mead
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 880
library(shiny)
## ---------------------------------------------------------------- objectives
mk <- function(label, fv, xlim, ylim, start, xstar) {
force(fv)
list(label = label, fv = fv, f = function(p) fv(p[1], p[2]),
xlim = xlim, ylim = ylim, start = start,
xstar = matrix(xstar, ncol = 2))
}
FNS <- list(
quad = mk(
"(a) Well-conditioned quadratic",
function(x, y) 0.5 * (2 * (x - 1)^2 + 1.6 * (x - 1) * (y - 1) +
2 * (y - 1)^2),
c(-3, 4), c(-3, 4), c(-2, 3), c(1, 1)
),
banana = mk(
"(b) Rosenbrock banana (narrow curved ridge)",
function(x, y) (1 - x)^2 + 100 * (y - x^2)^2,
c(-2, 2), c(-1, 3), c(-1.2, 1), c(1, 1)
),
mix = mk(
"(c) Highly correlated, two modes",
function(x, y) {
rho <- 0.9; den <- 1 - rho^2
q <- function(u, v) (u^2 - 2 * rho * u * v + v^2) / den
-log(exp(-0.5 * q(x - 1.2, y - 1.2)) +
0.6 * exp(-0.5 * q(x + 1.5, y + 1.5) / 1.5) + 1e-12)
},
c(-4, 4), c(-4, 4), c(-3, -0.5), c(1.2, 1.2)
),
himmel = mk(
"(d) Himmelblau (four global minima)",
function(x, y) (x^2 + y - 11)^2 + (x + y^2 - 7)^2,
c(-5, 5), c(-5, 5), c(-4, 4),
rbind(c(3, 2), c(-2.805118, 3.131312),
c(-3.779310, -3.283186), c(3.584428, -1.848126))
),
beale = mk(
"(e) Beale (flat plateau, sharp valley)",
function(x, y) (1.5 - x + x * y)^2 + (2.25 - x + x * y^2)^2 +
(2.625 - x + x * y^3)^2,
c(-4.5, 4.5), c(-4.5, 4.5), c(-1, 1), c(3, 0.5)
),
rastrigin = mk(
"(f) Rastrigin (many local minima)",
function(x, y) 20 + x^2 - 10 * cos(2 * pi * x) +
y^2 - 10 * cos(2 * pi * y),
c(-5.12, 5.12), c(-5.12, 5.12), c(-3.1, 4.2), c(0, 0)
)
)
## ------------------------------------------------------------- Nelder-Mead
## coefficients: reflection a, expansion g, contraction r, shrink s
run_nm <- function(FN, x0, h, a = 1, g = 2, r = 0.5, s = 0.5,
maxit = 80, tol = 1e-7) {
f <- FN$f
S <- rbind(x0, x0 + c(h, 0), x0 + c(0, h))
fS <- apply(S, 1, f)
steps <- list()
note <- "maximum number of iterations reached"
for (k in seq_len(maxit)) {
o <- order(fS); S <- S[o, , drop = FALSE]; fS <- fS[o]
dia <- max(c(sqrt(sum((S[1, ] - S[2, ])^2)),
sqrt(sum((S[1, ] - S[3, ])^2)),
sqrt(sum((S[2, ] - S[3, ])^2))))
if (dia < tol || diff(range(fS)) < 1e-12) {
steps[[k]] <- list(S = S, fS = fS, dia = dia, op = "converged",
cen = colMeans(S[1:2, , drop = FALSE]),
xr = NULL, xe = NULL, xc = NULL, Sn = S)
note <- "converged (simplex collapsed)"
break
}
cen <- colMeans(S[1:2, , drop = FALSE]) # centroid of the best two
xr <- cen + a * (cen - S[3, ]); fr <- f(xr)
xe <- NULL; xc <- NULL
Sn <- S; fn <- fS
if (fr < fS[1]) { # better than the best
xe <- cen + g * (xr - cen); fe <- f(xe)
if (fe < fr) { Sn[3, ] <- xe; fn[3] <- fe; op <- "expansion" }
else { Sn[3, ] <- xr; fn[3] <- fr; op <- "reflection" }
} else if (fr < fS[2]) { # middling: accept
Sn[3, ] <- xr; fn[3] <- fr; op <- "reflection"
} else {
if (fr < fS[3]) { # outside contraction
xc <- cen + r * (xr - cen); fc <- f(xc)
if (fc <= fr) { Sn[3, ] <- xc; fn[3] <- fc; op <- "outside contraction" }
else op <- "shrink"
} else { # inside contraction
xc <- cen + r * (S[3, ] - cen); fc <- f(xc)
if (fc < fS[3]) { Sn[3, ] <- xc; fn[3] <- fc; op <- "inside contraction" }
else op <- "shrink"
}
if (op == "shrink") {
Sn[2, ] <- S[1, ] + s * (S[2, ] - S[1, ])
Sn[3, ] <- S[1, ] + s * (S[3, ] - S[1, ])
fn[2] <- f(Sn[2, ]); fn[3] <- f(Sn[3, ])
}
}
steps[[k]] <- list(S = S, fS = fS, dia = dia, op = op, cen = cen,
xr = xr, xe = xe, xc = xc, Sn = Sn)
S <- Sn; fS <- fn
}
list(steps = steps, n = length(steps), note = note)
}
## ---------------------------- widen the data box to the device aspect ratio
expand_box <- function(xlim, ylim, pin) {
if (length(pin) != 2 || any(!is.finite(pin)) || any(pin <= 0))
return(list(xlim = xlim, ylim = ylim))
rr <- pin[1] / pin[2]
wd <- diff(xlim); ht <- diff(ylim)
if (wd / ht < rr) {
w <- ht * rr; m <- mean(xlim); xlim <- m + c(-0.5, 0.5) * w
} else {
hh <- wd / rr; m <- mean(ylim); ylim <- m + c(-0.5, 0.5) * hh
}
list(xlim = xlim, ylim = ylim)
}
COL <- c(refl = "#E8820C", exp = "#7B3FBF", con = "#1F6FB4",
new = "#2E8B57", best = "#2E8B57", worst = "#C81E1E")
## --------------------------------------------------------------------- UI
ui <- fluidPage(
titlePanel("Shinylive App for the Nelder–Mead Algorithm"),
fluidRow(
column(
3,
selectInput("fn", "Objective function",
setNames(names(FNS), sapply(FNS, `[[`, "label"))),
sliderInput("k", "Iteration k", min = 0, max = 1, value = 0, step = 1),
sliderInput("h", "Initial simplex size (% of range)",
min = 2, max = 40, value = 12, step = 1),
checkboxInput("showpath", "Keep past simplices", TRUE),
htmlOutput("info")
),
column(
2,
div(
style = "margin-top:26px;",
actionButton("play", "Play", class = "btn-primary btn-lg",
style = "width:100%; margin-bottom:10px;"),
actionButton("step", "Next", class = "btn-lg",
style = "width:100%; margin-bottom:10px;"),
actionButton("back", "Prev", class = "btn-lg",
style = "width:100%;"),
hr(),
checkboxInput("autoplay", "Autoplay on change", FALSE),
actionButton("reset", "Reset initial",
style = "width:100%; margin-bottom:8px;"),
helpText("Click inside the contour plot to choose a new initial value.")
)
),
column(
7,
plotOutput("contour", click = "click", height = "540px"),
plotOutput("conv", height = "190px")
)
)
)
## ----------------------------------------------------------------- server
server <- function(input, output, session) {
FN <- reactive(FNS[[input$fn]])
start <- reactiveVal(FNS[[1]]$start)
playing <- reactiveVal(FALSE)
observeEvent(input$fn, start(FNS[[input$fn]]$start))
observeEvent(input$reset, start(FN()$start))
observeEvent(input$click, {
p <- c(input$click$x, input$click$y)
if (all(is.finite(p))) start(p)
})
path <- reactive({
req(start())
FNc <- FN()
run_nm(FNc, start(), h = input$h / 100 * diff(FNc$xlim))
})
observeEvent(path(), {
updateSliderInput(session, "k", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
observeEvent(input$play, {
if (isTRUE(playing())) {
playing(FALSE)
} else {
if (as.integer(input$k) >= path()$n - 1)
updateSliderInput(session, "k", value = 0)
playing(TRUE)
}
})
observeEvent(playing(), {
updateActionButton(session, "play",
label = if (isTRUE(playing())) "Pause" else "Play")
})
observe({
if (!isTRUE(playing())) return()
kmax <- isolate(path()$n) - 1
k <- isolate(as.integer(input$k))
if (k >= kmax) { playing(FALSE); return() }
invalidateLater(900, session)
updateSliderInput(session, "k", value = k + 1)
})
observeEvent(input$step, {
playing(FALSE)
updateSliderInput(session, "k",
value = min(as.integer(input$k) + 1, max(0, path()$n - 1)))
})
observeEvent(input$back, {
playing(FALSE)
updateSliderInput(session, "k", value = max(as.integer(input$k) - 1, 0))
})
output$contour <- renderPlot({
FNc <- FN(); P <- path()
k <- min(as.integer(input$k), P$n - 1) + 1L
st <- P$steps[[k]]
par(mar = c(4, 4, 3, 1))
plot.new()
bx <- expand_box(FNc$xlim, FNc$ylim, par("pin"))
plot.window(xlim = bx$xlim, ylim = bx$ylim, xaxs = "i", yaxs = "i")
xs <- seq(bx$xlim[1], bx$xlim[2], length.out = 180)
ys <- seq(bx$ylim[1], bx$ylim[2], length.out = 180)
z <- outer(xs, ys, FNc$fv)
zf <- z[is.finite(z)]
lv <- unique(quantile(zf, probs = seq(0, 1, length.out = 30)^1.7))
zr <- matrix(rank(z, na.last = "keep", ties.method = "average"),
nrow = length(xs))
image(xs, ys, zr, add = TRUE,
col = colorRampPalette(c("#f7fbff", "#9ec4e3"))(128))
contour(xs, ys, z, levels = lv, add = TRUE,
col = "grey45", drawlabels = FALSE)
axis(1); axis(2); box()
title(main = paste0("Nelder-Mead | ", FNc$label),
xlab = expression(x[1]), ylab = expression(x[2]), cex.main = 1.0)
mtext(paste0("iteration ", k - 1, ": ", st$op), side = 3, line = 0.1,
cex = 1.1, font = 2)
points(FNc$xstar[, 1], FNc$xstar[, 2], pch = 8, col = "red",
cex = 1.3, lwd = 2)
## past simplices
if (isTRUE(input$showpath) && k > 1)
for (j in seq_len(k - 1))
polygon(P$steps[[j]]$S[, 1], P$steps[[j]]$S[, 2],
border = "grey55", lty = 3)
## current simplex
S <- st$S
polygon(S[, 1], S[, 2], border = "black", lwd = 2,
col = adjustcolor("white", alpha.f = 0.35))
points(st$cen[1], st$cen[2], pch = 3, col = "black", lwd = 2, cex = 1.1)
points(S[, 1], S[, 2], pch = 21, cex = 1.6, lwd = 2, bg = "white",
col = c(COL["best"], "grey30", COL["worst"]))
text(S[, 1], S[, 2], labels = c("B", "G", "W"), pos = 3, offset = 0.6,
font = 2, col = c(COL["best"], "grey30", COL["worst"]))
## candidate points generated this iteration
if (!is.null(st$xr)) {
segments(S[3, 1], S[3, 2], st$xr[1], st$xr[2],
col = COL["refl"], lwd = 2, lty = 2)
points(st$xr[1], st$xr[2], pch = 19, col = COL["refl"], cex = 1.5)
}
if (!is.null(st$xe)) {
segments(st$xr[1], st$xr[2], st$xe[1], st$xe[2],
col = COL["exp"], lwd = 2, lty = 2)
points(st$xe[1], st$xe[2], pch = 19, col = COL["exp"], cex = 1.5)
}
if (!is.null(st$xc))
points(st$xc[1], st$xc[2], pch = 19, col = COL["con"], cex = 1.5)
## the simplex that results
if (st$op != "converged") {
polygon(st$Sn[, 1], st$Sn[, 2], border = COL["new"], lwd = 3, lty = 2)
if (st$op == "shrink")
arrows(S[2:3, 1], S[2:3, 2], st$Sn[2:3, 1], st$Sn[2:3, 2],
col = COL["con"], lwd = 2, length = 0.10)
}
legend("topleft", cex = 0.85, bg = "#ffffffcc", box.col = NA,
pch = c(21, 21, 3, 19, 19, 19, NA),
lty = c(NA, NA, NA, NA, NA, NA, 2),
lwd = c(2, 2, 2, NA, NA, NA, 3),
col = c(COL["best"], COL["worst"], "black", COL["refl"],
COL["exp"], COL["con"], COL["new"]),
legend = c("best vertex B", "worst vertex W",
"centroid of B and G", "reflection", "expansion",
"contraction / shrink", "next simplex"))
})
output$conv <- renderPlot({
P <- path(); k <- min(as.integer(input$k), P$n - 1) + 1L
fb <- sapply(P$steps, function(s) s$fS[1])
fw <- sapply(P$steps, function(s) s$fS[3])
par(mar = c(4, 4.5, 1.5, 1))
matplot(seq_len(P$n) - 1, cbind(fb, fw), type = "b", pch = 20, lty = 1,
col = c(COL["best"], COL["worst"]),
xlab = "iteration k", ylab = "f at simplex vertices")
abline(v = k - 1, col = "grey60", lty = 2)
legend("topright", bty = "n", cex = 0.9, lty = 1, pch = 20,
col = c(COL["best"], COL["worst"]), legend = c("best", "worst"))
})
output$info <- renderUI({
P <- path(); k <- min(as.integer(input$k), P$n - 1) + 1L
st <- P$steps[[k]]
fmt <- function(v, d = 4) formatC(v, format = "g", digits = d)
HTML(paste0(
"<hr><b>iteration</b> ", k - 1, " of ", P$n - 1, "<br>",
"<b>operation</b>: ", st$op, "<br><br>",
"<b>B</b> = (", fmt(st$S[1, 1]), ", ", fmt(st$S[1, 2]), "), f = ",
fmt(st$fS[1], 6), "<br>",
"<b>G</b> = (", fmt(st$S[2, 1]), ", ", fmt(st$S[2, 2]), "), f = ",
fmt(st$fS[2], 6), "<br>",
"<b>W</b> = (", fmt(st$S[3, 1]), ", ", fmt(st$S[3, 2]), "), f = ",
fmt(st$fS[3], 6), "<br><br>",
"<b>simplex diameter</b> = ", fmt(st$dia), "<br>",
"<b>f range</b> = ", fmt(diff(range(st$fS))),
"<br><br><i>", P$note, "</i>"
))
})
}
shinyApp(ui, server)About the app
Shows how the Nelder-Mead method finds a minimum without derivatives, by repeatedly reflecting, expanding, contracting, and shrinking a simplex. Reflection, expansion, contraction and shrink steps are colour-coded on the contour plot as the simplex moves.
The app above animates every reflection, expansion, contraction, and shrink on six test objectives; click inside the panel to place a new starting simplex.
NoteR source for this app
library(shiny)
## ---------------------------------------------------------------- objectives
mk <- function(label, fv, xlim, ylim, start, xstar) {
force(fv)
list(label = label, fv = fv, f = function(p) fv(p[1], p[2]),
xlim = xlim, ylim = ylim, start = start,
xstar = matrix(xstar, ncol = 2))
}
FNS <- list(
quad = mk(
"(a) Well-conditioned quadratic",
function(x, y) 0.5 * (2 * (x - 1)^2 + 1.6 * (x - 1) * (y - 1) +
2 * (y - 1)^2),
c(-3, 4), c(-3, 4), c(-2, 3), c(1, 1)
),
banana = mk(
"(b) Rosenbrock banana (narrow curved ridge)",
function(x, y) (1 - x)^2 + 100 * (y - x^2)^2,
c(-2, 2), c(-1, 3), c(-1.2, 1), c(1, 1)
),
mix = mk(
"(c) Highly correlated, two modes",
function(x, y) {
rho <- 0.9; den <- 1 - rho^2
q <- function(u, v) (u^2 - 2 * rho * u * v + v^2) / den
-log(exp(-0.5 * q(x - 1.2, y - 1.2)) +
0.6 * exp(-0.5 * q(x + 1.5, y + 1.5) / 1.5) + 1e-12)
},
c(-4, 4), c(-4, 4), c(-3, -0.5), c(1.2, 1.2)
),
himmel = mk(
"(d) Himmelblau (four global minima)",
function(x, y) (x^2 + y - 11)^2 + (x + y^2 - 7)^2,
c(-5, 5), c(-5, 5), c(-4, 4),
rbind(c(3, 2), c(-2.805118, 3.131312),
c(-3.779310, -3.283186), c(3.584428, -1.848126))
),
beale = mk(
"(e) Beale (flat plateau, sharp valley)",
function(x, y) (1.5 - x + x * y)^2 + (2.25 - x + x * y^2)^2 +
(2.625 - x + x * y^3)^2,
c(-4.5, 4.5), c(-4.5, 4.5), c(-1, 1), c(3, 0.5)
),
rastrigin = mk(
"(f) Rastrigin (many local minima)",
function(x, y) 20 + x^2 - 10 * cos(2 * pi * x) +
y^2 - 10 * cos(2 * pi * y),
c(-5.12, 5.12), c(-5.12, 5.12), c(-3.1, 4.2), c(0, 0)
)
)
## ------------------------------------------------------------- Nelder-Mead
## coefficients: reflection a, expansion g, contraction r, shrink s
run_nm <- function(FN, x0, h, a = 1, g = 2, r = 0.5, s = 0.5,
maxit = 80, tol = 1e-7) {
f <- FN$f
S <- rbind(x0, x0 + c(h, 0), x0 + c(0, h))
fS <- apply(S, 1, f)
steps <- list()
note <- "maximum number of iterations reached"
for (k in seq_len(maxit)) {
o <- order(fS); S <- S[o, , drop = FALSE]; fS <- fS[o]
dia <- max(c(sqrt(sum((S[1, ] - S[2, ])^2)),
sqrt(sum((S[1, ] - S[3, ])^2)),
sqrt(sum((S[2, ] - S[3, ])^2))))
if (dia < tol || diff(range(fS)) < 1e-12) {
steps[[k]] <- list(S = S, fS = fS, dia = dia, op = "converged",
cen = colMeans(S[1:2, , drop = FALSE]),
xr = NULL, xe = NULL, xc = NULL, Sn = S)
note <- "converged (simplex collapsed)"
break
}
cen <- colMeans(S[1:2, , drop = FALSE]) # centroid of the best two
xr <- cen + a * (cen - S[3, ]); fr <- f(xr)
xe <- NULL; xc <- NULL
Sn <- S; fn <- fS
if (fr < fS[1]) { # better than the best
xe <- cen + g * (xr - cen); fe <- f(xe)
if (fe < fr) { Sn[3, ] <- xe; fn[3] <- fe; op <- "expansion" }
else { Sn[3, ] <- xr; fn[3] <- fr; op <- "reflection" }
} else if (fr < fS[2]) { # middling: accept
Sn[3, ] <- xr; fn[3] <- fr; op <- "reflection"
} else {
if (fr < fS[3]) { # outside contraction
xc <- cen + r * (xr - cen); fc <- f(xc)
if (fc <= fr) { Sn[3, ] <- xc; fn[3] <- fc; op <- "outside contraction" }
else op <- "shrink"
} else { # inside contraction
xc <- cen + r * (S[3, ] - cen); fc <- f(xc)
if (fc < fS[3]) { Sn[3, ] <- xc; fn[3] <- fc; op <- "inside contraction" }
else op <- "shrink"
}
if (op == "shrink") {
Sn[2, ] <- S[1, ] + s * (S[2, ] - S[1, ])
Sn[3, ] <- S[1, ] + s * (S[3, ] - S[1, ])
fn[2] <- f(Sn[2, ]); fn[3] <- f(Sn[3, ])
}
}
steps[[k]] <- list(S = S, fS = fS, dia = dia, op = op, cen = cen,
xr = xr, xe = xe, xc = xc, Sn = Sn)
S <- Sn; fS <- fn
}
list(steps = steps, n = length(steps), note = note)
}
## ---------------------------- widen the data box to the device aspect ratio
expand_box <- function(xlim, ylim, pin) {
if (length(pin) != 2 || any(!is.finite(pin)) || any(pin <= 0))
return(list(xlim = xlim, ylim = ylim))
rr <- pin[1] / pin[2]
wd <- diff(xlim); ht <- diff(ylim)
if (wd / ht < rr) {
w <- ht * rr; m <- mean(xlim); xlim <- m + c(-0.5, 0.5) * w
} else {
hh <- wd / rr; m <- mean(ylim); ylim <- m + c(-0.5, 0.5) * hh
}
list(xlim = xlim, ylim = ylim)
}
COL <- c(refl = "#E8820C", exp = "#7B3FBF", con = "#1F6FB4",
new = "#2E8B57", best = "#2E8B57", worst = "#C81E1E")
## --------------------------------------------------------------------- UI
ui <- fluidPage(
titlePanel("Shinylive App for the Nelder–Mead Algorithm"),
fluidRow(
column(
3,
selectInput("fn", "Objective function",
setNames(names(FNS), sapply(FNS, `[[`, "label"))),
sliderInput("k", "Iteration k", min = 0, max = 1, value = 0, step = 1),
sliderInput("h", "Initial simplex size (% of range)",
min = 2, max = 40, value = 12, step = 1),
checkboxInput("showpath", "Keep past simplices", TRUE),
htmlOutput("info")
),
column(
2,
div(
style = "margin-top:26px;",
actionButton("play", "Play", class = "btn-primary btn-lg",
style = "width:100%; margin-bottom:10px;"),
actionButton("step", "Next", class = "btn-lg",
style = "width:100%; margin-bottom:10px;"),
actionButton("back", "Prev", class = "btn-lg",
style = "width:100%;"),
hr(),
checkboxInput("autoplay", "Autoplay on change", FALSE),
actionButton("reset", "Reset initial",
style = "width:100%; margin-bottom:8px;"),
helpText("Click inside the contour plot to choose a new initial value.")
)
),
column(
7,
plotOutput("contour", click = "click", height = "540px"),
plotOutput("conv", height = "190px")
)
)
)
## ----------------------------------------------------------------- server
server <- function(input, output, session) {
FN <- reactive(FNS[[input$fn]])
start <- reactiveVal(FNS[[1]]$start)
playing <- reactiveVal(FALSE)
observeEvent(input$fn, start(FNS[[input$fn]]$start))
observeEvent(input$reset, start(FN()$start))
observeEvent(input$click, {
p <- c(input$click$x, input$click$y)
if (all(is.finite(p))) start(p)
})
path <- reactive({
req(start())
FNc <- FN()
run_nm(FNc, start(), h = input$h / 100 * diff(FNc$xlim))
})
observeEvent(path(), {
updateSliderInput(session, "k", max = max(1, path()$n - 1), value = 0)
playing(isTRUE(input$autoplay) && path()$n > 1)
})
observeEvent(input$play, {
if (isTRUE(playing())) {
playing(FALSE)
} else {
if (as.integer(input$k) >= path()$n - 1)
updateSliderInput(session, "k", value = 0)
playing(TRUE)
}
})
observeEvent(playing(), {
updateActionButton(session, "play",
label = if (isTRUE(playing())) "Pause" else "Play")
})
observe({
if (!isTRUE(playing())) return()
kmax <- isolate(path()$n) - 1
k <- isolate(as.integer(input$k))
if (k >= kmax) { playing(FALSE); return() }
invalidateLater(900, session)
updateSliderInput(session, "k", value = k + 1)
})
observeEvent(input$step, {
playing(FALSE)
updateSliderInput(session, "k",
value = min(as.integer(input$k) + 1, max(0, path()$n - 1)))
})
observeEvent(input$back, {
playing(FALSE)
updateSliderInput(session, "k", value = max(as.integer(input$k) - 1, 0))
})
output$contour <- renderPlot({
FNc <- FN(); P <- path()
k <- min(as.integer(input$k), P$n - 1) + 1L
st <- P$steps[[k]]
par(mar = c(4, 4, 3, 1))
plot.new()
bx <- expand_box(FNc$xlim, FNc$ylim, par("pin"))
plot.window(xlim = bx$xlim, ylim = bx$ylim, xaxs = "i", yaxs = "i")
xs <- seq(bx$xlim[1], bx$xlim[2], length.out = 180)
ys <- seq(bx$ylim[1], bx$ylim[2], length.out = 180)
z <- outer(xs, ys, FNc$fv)
zf <- z[is.finite(z)]
lv <- unique(quantile(zf, probs = seq(0, 1, length.out = 30)^1.7))
zr <- matrix(rank(z, na.last = "keep", ties.method = "average"),
nrow = length(xs))
image(xs, ys, zr, add = TRUE,
col = colorRampPalette(c("#f7fbff", "#9ec4e3"))(128))
contour(xs, ys, z, levels = lv, add = TRUE,
col = "grey45", drawlabels = FALSE)
axis(1); axis(2); box()
title(main = paste0("Nelder-Mead | ", FNc$label),
xlab = expression(x[1]), ylab = expression(x[2]), cex.main = 1.0)
mtext(paste0("iteration ", k - 1, ": ", st$op), side = 3, line = 0.1,
cex = 1.1, font = 2)
points(FNc$xstar[, 1], FNc$xstar[, 2], pch = 8, col = "red",
cex = 1.3, lwd = 2)
## past simplices
if (isTRUE(input$showpath) && k > 1)
for (j in seq_len(k - 1))
polygon(P$steps[[j]]$S[, 1], P$steps[[j]]$S[, 2],
border = "grey55", lty = 3)
## current simplex
S <- st$S
polygon(S[, 1], S[, 2], border = "black", lwd = 2,
col = adjustcolor("white", alpha.f = 0.35))
points(st$cen[1], st$cen[2], pch = 3, col = "black", lwd = 2, cex = 1.1)
points(S[, 1], S[, 2], pch = 21, cex = 1.6, lwd = 2, bg = "white",
col = c(COL["best"], "grey30", COL["worst"]))
text(S[, 1], S[, 2], labels = c("B", "G", "W"), pos = 3, offset = 0.6,
font = 2, col = c(COL["best"], "grey30", COL["worst"]))
## candidate points generated this iteration
if (!is.null(st$xr)) {
segments(S[3, 1], S[3, 2], st$xr[1], st$xr[2],
col = COL["refl"], lwd = 2, lty = 2)
points(st$xr[1], st$xr[2], pch = 19, col = COL["refl"], cex = 1.5)
}
if (!is.null(st$xe)) {
segments(st$xr[1], st$xr[2], st$xe[1], st$xe[2],
col = COL["exp"], lwd = 2, lty = 2)
points(st$xe[1], st$xe[2], pch = 19, col = COL["exp"], cex = 1.5)
}
if (!is.null(st$xc))
points(st$xc[1], st$xc[2], pch = 19, col = COL["con"], cex = 1.5)
## the simplex that results
if (st$op != "converged") {
polygon(st$Sn[, 1], st$Sn[, 2], border = COL["new"], lwd = 3, lty = 2)
if (st$op == "shrink")
arrows(S[2:3, 1], S[2:3, 2], st$Sn[2:3, 1], st$Sn[2:3, 2],
col = COL["con"], lwd = 2, length = 0.10)
}
legend("topleft", cex = 0.85, bg = "#ffffffcc", box.col = NA,
pch = c(21, 21, 3, 19, 19, 19, NA),
lty = c(NA, NA, NA, NA, NA, NA, 2),
lwd = c(2, 2, 2, NA, NA, NA, 3),
col = c(COL["best"], COL["worst"], "black", COL["refl"],
COL["exp"], COL["con"], COL["new"]),
legend = c("best vertex B", "worst vertex W",
"centroid of B and G", "reflection", "expansion",
"contraction / shrink", "next simplex"))
})
output$conv <- renderPlot({
P <- path(); k <- min(as.integer(input$k), P$n - 1) + 1L
fb <- sapply(P$steps, function(s) s$fS[1])
fw <- sapply(P$steps, function(s) s$fS[3])
par(mar = c(4, 4.5, 1.5, 1))
matplot(seq_len(P$n) - 1, cbind(fb, fw), type = "b", pch = 20, lty = 1,
col = c(COL["best"], COL["worst"]),
xlab = "iteration k", ylab = "f at simplex vertices")
abline(v = k - 1, col = "grey60", lty = 2)
legend("topright", bty = "n", cex = 0.9, lty = 1, pch = 20,
col = c(COL["best"], COL["worst"]), legend = c("best", "worst"))
})
output$info <- renderUI({
P <- path(); k <- min(as.integer(input$k), P$n - 1) + 1L
st <- P$steps[[k]]
fmt <- function(v, d = 4) formatC(v, format = "g", digits = d)
HTML(paste0(
"<hr><b>iteration</b> ", k - 1, " of ", P$n - 1, "<br>",
"<b>operation</b>: ", st$op, "<br><br>",
"<b>B</b> = (", fmt(st$S[1, 1]), ", ", fmt(st$S[1, 2]), "), f = ",
fmt(st$fS[1], 6), "<br>",
"<b>G</b> = (", fmt(st$S[2, 1]), ", ", fmt(st$S[2, 2]), "), f = ",
fmt(st$fS[2], 6), "<br>",
"<b>W</b> = (", fmt(st$S[3, 1]), ", ", fmt(st$S[3, 2]), "), f = ",
fmt(st$fS[3], 6), "<br><br>",
"<b>simplex diameter</b> = ", fmt(st$dia), "<br>",
"<b>f range</b> = ", fmt(diff(range(st$fS))),
"<br><br><i>", P$note, "</i>"
))
})
}
shinyApp(ui, server)This app accompanies Optimization for Maximum Likelihood Estimation in the book.