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

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

fr <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
grr <- function(x) c(-400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1]),
                     200 * (x[2] - x[1]^2))
hr <- function(x) matrix(c(1200 * x[1]^2 - 400 * x[2] + 2, -400 * x[1],
                           -400 * x[1], 200), 2, 2)

r <- steihaug(c(-1.2, 1), fr, grr, hr)
r[c("par", "value", "iterations", "message", "counts")]

## -----------------------------------------------------------------------------
hv <- function(x, v) drop(hr(x) %*% v)
m <- steihaug(c(-1.2, 1), fr, grr, hessvec = hv)
identical(m$par, r$par)
m$counts

## -----------------------------------------------------------------------------
tr <- steihaug(c(-1.2, 1), fr, grr, hr, control = list(trace = TRUE))$trace
head(tr)

## -----------------------------------------------------------------------------
fr <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
bobyqa(c(1, 2), fr, lower = c(0, 0), upper = c(4, 4))
newuoa(c(-1.2, 1), fr)

## ----eval = requireNamespace("minqa", quietly = TRUE)-------------------------
extRosen <- function(x) {
  n <- length(x)
  sum(100 * (x[-1] - x[-n]^2)^2 + (1 - x[-n])^2)
}
ours <- bobyqa(rep(0.5, 6), extRosen, lower = -2, upper = 2,
               control = list(npt = 13))
theirs <- minqa::bobyqa(rep(0.5, 6), extRosen, lower = -2, upper = 2,
                        control = list(npt = 13))
identical(unclass(ours), unclass(theirs))

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

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

// Derivative-free objective for NEWUOA.
extern "C" int rosen_f(int n, const double *x, double *f, void *) {
  double t = x[1] - x[0] * x[0];
  *f = 100.0 * t * t + (1.0 - x[0]) * (1.0 - x[0]);
  return 0;
}

// Value/gradient on request, for Steihaug. The Hessian is never requested
// in matrix-free mode.
extern "C" int rosen_vg(int n, const double *x, double *value,
                        double *gradient, double *hessian, void *) {
  double t = x[1] - x[0] * x[0];
  if (value) *value = 100.0 * t * t + (1.0 - x[0]) * (1.0 - x[0]);
  if (gradient) {
    gradient[0] = -400.0 * x[0] * t - 2.0 * (1.0 - x[0]);
    gradient[1] = 200.0 * t;
  }
  return 0;
}

extern "C" int rosen_hv(int n, const double *x, const double *v, double *hv,
                        void *) {
  double h11 = 1200.0 * x[0] * x[0] - 400.0 * x[1] + 2.0, h12 = -400.0 * x[0];
  hv[0] = h11 * v[0] + h12 * v[1];
  hv[1] = h12 * v[0] + 200.0 * v[1];
  return 0;
}

// [[Rcpp::export]]
Rcpp::NumericMatrix fit_many(SEXP ptrTable, int nStarts) {
  _vignette2_iniRcppTrustPtrs(ptrTable);  // normally done once in .onLoad()
  Rcpp::NumericMatrix out(nStarts, 4);
  std::vector<double> res(4 * nStarts);
  // Build with OpenMP (e.g. -fopenmp in the package Makevars) to run
  // this loop in parallel; without it the pragma is ignored.
  #pragma omp parallel for num_threads(2)
  for (int i = 0; i < nStarts; i++) {
    double start[2] = {-1.2 + 0.05 * i, 1.0};

    minqa_options_t mo = minqa_options_default(2, start);
    minqa_result_t mr;
    newuoa_solve_c_ptr(2, start, rosen_f, nullptr, &mo, &mr);

    steihaug_options_t so = steihaug_options_default();
    steihaug_result_t sr;
    steihaug_solve_c_ptr(2, start, rosen_vg, rosen_hv, nullptr, &so, &sr);

    res[4 * i] = mr.par[0];
    res[4 * i + 1] = mr.feval;
    res[4 * i + 2] = sr.argument[0];
    res[4 * i + 3] = sr.iterations;
    minqa_result_free_ptr(&mr);
    steihaug_result_free_ptr(&sr);
  }
  for (int i = 0; i < nStarts; i++)
    for (int j = 0; j < 4; j++) out(i, j) = res[4 * i + j];
  Rcpp::colnames(out) = Rcpp::CharacterVector::create(
      "newuoa_x1", "newuoa_feval", "steihaug_x1", "steihaug_iter");
  return out;
}
'
Rcpp::sourceCpp(code = cpp_code)
head(fit_many(RcppTrust:::.RcppTrustPtr(), 16L))