## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(arcopt)

## ----mixture-saddle-----------------------------------------------------------
set.seed(42)
y <- c(rnorm(100, -2, 1), rnorm(100, 2, 1))

mixture_nll <- function(theta) {
  -sum(log(0.5 * dnorm(y, theta[1], 1) + 0.5 * dnorm(y, theta[2], 1)))
}
mixture_gr <- function(theta) {
  p1 <- 0.5 * dnorm(y, theta[1], 1)
  p2 <- 0.5 * dnorm(y, theta[2], 1)
  w1 <- p1 / (p1 + p2)
  w2 <- p2 / (p1 + p2)
  -c(sum(w1 * (y - theta[1])), sum(w2 * (y - theta[2])))
}
mixture_hess <- function(theta, h = 1e-5) {
  d <- 2L
  h_mat <- matrix(0, d, d)
  for (i in seq_len(d)) {
    e_i <- numeric(d)
    e_i[i] <- h
    h_mat[, i] <- (mixture_gr(theta + e_i) - mixture_gr(theta - e_i)) /
      (2 * h)
  }
  0.5 * (h_mat + t(h_mat))
}

res <- arcopt(c(0.01, 0.01), mixture_nll, mixture_gr, mixture_hess,
              control = list(trace = 0))
res$par               # near (-2, 2) modulo label permutation
res$diagnostics$solver_mode_final
min(eigen(res$hessian, only.values = TRUE)$values)  # > 0: local min

## ----tr-dormant---------------------------------------------------------------
rosen_fn <- function(x) (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
rosen_gr <- function(x) {
  c(-2 * (1 - x[1]) - 400 * x[1] * (x[2] - x[1]^2),
    200 * (x[2] - x[1]^2))
}
rosen_hess <- function(x) {
  matrix(c(1200 * x[1]^2 - 400 * x[2] + 2, -400 * x[1],
           -400 * x[1],                     200), 2, 2)
}

res <- arcopt(c(-1.2, 1), rosen_fn, rosen_gr, rosen_hess,
              control = list(trace = 0))
res$diagnostics$solver_mode_final
res$diagnostics$ridge_switches

## ----eval = FALSE-------------------------------------------------------------
# arcopt(x0, fn, gr, hess,
#        control = list(tr_fallback_enabled = FALSE))

## ----polish-------------------------------------------------------------------
fn   <- function(x) sum(0.5 * x^2 + x^4 / 12)
gr   <- function(x) x + x^3 / 3
hess <- function(x) diag(1 + x^2)
x0   <- rep(10, 5)

res_default <- arcopt(x0, fn, gr, hess,
                      control = list(maxit = 50, trace = 0))
res_polish  <- arcopt(x0, fn, gr, hess,
                      control = list(maxit = 50, trace = 0,
                                     qn_polish_enabled = TRUE))

c(default_hess = res_default$evaluations$hess,
  polish_hess  = res_polish$evaluations$hess)

res_polish$diagnostics$solver_mode_final     # "qn_polish"
res_polish$diagnostics$qn_polish_switches    # 1