## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(rzig)
compiler <- find_zig(required = FALSE)
has_zig <- isTRUE(compiler$supported)
# LAPACK and BLAS need the Fortran runtime that R links ($(FLIBS)).
flibs <- suppressWarnings(system2(file.path(R.home("bin"), "R"), c("CMD", "config", "FLIBS"), stdout = TRUE))
flib_dirs <- sub("^-L", "", grep("^-L", strsplit(paste(flibs, collapse = " "), "[[:space:]]+")[[1L]], value = TRUE))
has_flibs <- length(flib_dirs) == 0L || all(dir.exists(flib_dirs))

## ----create-------------------------------------------------------------------
path <- file.path(tempdir(), "numericspkg")
create_rzig_package(path, tests = FALSE)
main_zig <- file.path(path, "src", "rzig", "src", "main.zig")
use_rzig_lapack(path)

## ----rmath, comment = ""------------------------------------------------------
cat(
  "",
  "/// The upper tail of a gamma distribution with shape 2.",
  "/// @param x Quantiles.",
  "/// @return Probabilities.",
  "/// @export",
  "pub fn gamma_upper(ctx: *rzig.Ctx, x: []const f64) rzig.Error![]f64 {",
  "    const out = try ctx.alloc(f64, x.len);",
  "    for (x, out) |value, *slot| {",
  "        slot.* = try rzig.Rmath.cdf(.gamma, value, .{ .shape = 2 }, .{ .lower_tail = false });",
  "    }",
  "    return out;",
  "}",
  sep = "\n", file = main_zig, append = TRUE
)

## ----lapack, comment = ""-----------------------------------------------------
cat(
  "",
  "/// Solve a symmetric positive-definite system with a Cholesky factor.",
  "/// @param a A symmetric positive-definite matrix.",
  "/// @param b A right-hand side.",
  "/// @return The solution.",
  "/// @export",
  "pub fn chol_solve(ctx: *rzig.Ctx, a: rzig.Matrix, b: []const f64) rzig.Error![]f64 {",
  "    const factor = try ctx.matrix(a.nrow, a.ncol);",
  "    @memcpy(factor.data, a.data);",
  "    try rzig.Lapack.potrf(ctx, .upper, factor);",
  "    const x = try ctx.matrix(b.len, 1);",
  "    @memcpy(x.data, b);",
  "    try rzig.Lapack.potrs(ctx, .upper, factor.view(), x);",
  "    return x.data;",
  "}",
  sep = "\n", file = main_zig, append = TRUE
)

## ----optim, comment = ""------------------------------------------------------
cat(
  "",
  "const RObjective = struct {",
  "    ctx: *rzig.Ctx,",
  "    objective: rzig.Function,",
  "",
  "    pub fn value(self: *const RObjective, x: []const f64) rzig.Error!f64 {",
  "        return self.objective.call(self.ctx, f64, .{x});",
  "    }",
  "};",
  "",
  "/// Minimize an R function with BFGS.",
  "/// @param f An R function of a numeric vector.",
  "/// @param start Starting values.",
  "/// @return The minimizer.",
  "/// @export",
  "pub fn minimize_r(ctx: *rzig.Ctx, f: rzig.Function, start: []const f64) rzig.Error![]f64 {",
  "    const x = try ctx.dupe(f64, start);",
  "    _ = try rzig.Optimize.minimize(ctx, .bfgs, &RObjective{ .ctx = ctx, .objective = f }, x, .{});",
  "    return x;",
  "}",
  sep = "\n", file = main_zig, append = TRUE
)

## ----integrate, comment = ""--------------------------------------------------
cat(
  "",
  "const Exponential = struct {",
  "    rate: f64,",
  "",
  "    pub fn value(self: *const Exponential, x: f64) f64 {",
  "        return self.rate * @exp(-self.rate * x);",
  "    }",
  "};",
  "",
  "/// The probability of an exponential variable between two limits.",
  "/// @param rate The rate.",
  "/// @param lower,upper The limits.",
  "/// @return A probability.",
  "/// @export",
  "pub fn exponential_mass(ctx: *rzig.Ctx, rate: f64, lower: f64, upper: f64) rzig.Error!f64 {",
  "    const result = try rzig.Integrate.integrate(ctx, &Exponential{ .rate = rate }, lower, upper, .{});",
  "    return result.value;",
  "}",
  sep = "\n", file = main_zig, append = TRUE
)
document(path)

## ----build, eval = has_zig && has_flibs---------------------------------------
build <- build_zig(path)
library(build$package, character.only = TRUE, lib.loc = build$library)
all.equal(gamma_upper(c(0.5, 2)), pgamma(c(0.5, 2), shape = 2, lower.tail = FALSE))
a <- crossprod(matrix(c(2, 1, 0, 1, 3, 1, 0, 1, 4), 3)) + diag(3)
all.equal(chol_solve(a, c(1, 2, 3)), solve(a, c(1, 2, 3)))
try(chol_solve(matrix(c(1, 2, 2, 1), 2), c(1, 1)))
rosenbrock <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
minimize_r(rosenbrock, c(-1.2, 1))
optim(c(-1.2, 1), rosenbrock, method = "BFGS")$par
exponential_mass(2, 0, Inf)
exponential_mass(2, 0, 1) - pexp(1, 2)
tryCatch(minimize_r(function(x) stop("the objective failed"), c(1, 1)),
  error = conditionMessage)

## ----cleanup, include = FALSE-------------------------------------------------
if (has_zig && has_flibs) unloadNamespace("numericspkg")
unlink(path, recursive = TRUE)