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

## -----------------------------------------------------------------------------
library(RcppTrust)

# Rosenbrock's function, exactly the example from ?trust::trust
objfun <- function(x) {
  f <- expression(100 * (x2 - x1^2)^2 + (1 - x1)^2)
  g1 <- D(f, "x1"); g2 <- D(f, "x2")
  h11 <- D(g1, "x1"); h12 <- D(g1, "x2"); h22 <- D(g2, "x2")
  x1 <- x[1]; x2 <- x[2]
  list(
    value = eval(f), gradient = c(eval(g1), eval(g2)),
    hessian = rbind(c(eval(h11), eval(h12)), c(eval(h12), eval(h22)))
  )
}

out <- trust(objfun, c(3, 1), 1, 5)
out[c("value", "argument", "converged", "iterations")]

## ----eval = requireNamespace("trust", quietly = TRUE) && requireNamespace("microbenchmark", quietly = TRUE)----
library(microbenchmark)

# same Rosenbrock objfun as above: 2 parameters, ~20 iterations
mb_small <- microbenchmark(
  trust = trust::trust(objfun, c(3, 1), 1, 5),
  RcppTrust = RcppTrust::trust(objfun, c(3, 1), 1, 5),
  times = 50
)
print(mb_small)

## ----eval = requireNamespace("trust", quietly = TRUE) && requireNamespace("microbenchmark", quietly = TRUE)----
d <- 30
mu <- seq_len(d)
barrier_objfun <- function(x) {
  normxsq <- sum(x^2)
  omnormxsq <- 1 - normxsq
  if (normxsq >= 1) return(list(value = Inf))
  f <- sum(x * mu) - log(omnormxsq)
  g <- mu + 2 * x / omnormxsq
  B <- 4 * outer(x, x) / omnormxsq^2 + 2 * diag(d) / omnormxsq
  list(value = f, gradient = g, hessian = B)
}

r1 <- trust::trust(barrier_objfun, rep(0, d), 1, 100)
r2 <- RcppTrust::trust(barrier_objfun, rep(0, d), 1, 100)
# same solution, different number of steps to get there (see above)
c(trust_iterations = r1$iterations, RcppTrust_iterations = r2$iterations)
max(abs(r1$argument - r2$argument))

mb_large <- microbenchmark(
  trust = trust::trust(barrier_objfun, rep(0, d), 1, 100),
  RcppTrust = RcppTrust::trust(barrier_objfun, rep(0, d), 1, 100),
  times = 30
)
print(mb_large)

## -----------------------------------------------------------------------------
cpp_code <- '
// [[Rcpp::depends(RcppTrust)]]
#include <Rcpp.h>

extern "C" {
#define iniRcppTrustPtrs _vignette_iniRcppTrustPtrs
#include <RcppTrust.h>
iniRcppTrust
}

// A thread-safe C objective function: Rosenbrock again, this time
// filling value/gradient/hessian directly instead of returning a list.
extern "C" int rosenbrock_c(int n, const double *par, double *value,
                             double *gradient, double *hessian, void *ud) {
  double x1 = par[0], x2 = par[1];
  double t = x2 - x1 * x1;
  *value = 100.0 * t * t + (1.0 - x1) * (1.0 - x1);
  gradient[0] = -400.0 * x1 * t - 2.0 * (1.0 - x1);
  gradient[1] = 200.0 * t;
  hessian[0] = -400.0 * x2 + 1200.0 * x1 * x1 + 2.0;
  hessian[1] = hessian[2] = -400.0 * x1;
  hessian[3] = 200.0;
  return 0;
}

// [[Rcpp::export]]
Rcpp::List fit_rosenbrock(SEXP ptrTable, Rcpp::NumericVector parinit) {
  // normally done once, in .onLoad() -- see above
  _vignette_iniRcppTrustPtrs(ptrTable);

  trust_options_t opts = trust_options_default(1.0, 5.0);
  trust_result_t res;
  trust_solve_c_ptr(parinit.size(), parinit.begin(), rosenbrock_c, nullptr, &opts, &res);

  Rcpp::List out = Rcpp::List::create(
    Rcpp::_["argument"] = Rcpp::NumericVector(res.argument, res.argument + res.n),
    Rcpp::_["value"] = res.value,
    Rcpp::_["converged"] = res.converged != 0,
    Rcpp::_["iterations"] = res.iterations);
  trust_result_free_ptr(&res);
  return out;
}
'
Rcpp::sourceCpp(code = cpp_code)

fit_rosenbrock(RcppTrust:::.RcppTrustPtr(), c(3, 1))