---
title: "Numerical routines of R in Zig code"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Numerical routines of R in Zig code}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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))
```

Zig code in an R package can use the numerical routines that R itself uses:
the distributions of R's maths library, the BLAS and LAPACK that R links, the
optimizers of `optim()`, the integration of `integrate()`, and R functions
passed as arguments. Each routine runs with R's errors turned into errors of
the Zig call, so the Zig code returns them and R reports them after the Zig
frames have returned. The package below uses each of them; the build and the
calls run when a supported Zig compiler is found (this copy
`r if (has_zig) "ran" else "skipped"` them).

```{r 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)
```

`use_rzig_lapack()` asks `configure` to link LAPACK, BLAS and the Fortran
runtime, as *Writing R Extensions* prescribes.

## Distributions

`rzig.Rmath` evaluates the density, the distribution function and the
quantile function of R's distributions with R's parameter names and defaults:

```{r 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
)
```

## Linear algebra

`rzig.Blas` and `rzig.Lapack` call R's BLAS and LAPACK on column-major
matrices and check the shapes before each call. A Cholesky solve:

```{r 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
)
```

A matrix that is not positive definite is an R error with LAPACK's reason.

## Optimization and R functions as arguments

`rzig.Optimize.minimize` runs the optimizers of `optim()` on a problem with a
`value` method and, optionally, a `gradient` method; without one, it uses the
central differences of `optim()`. An `rzig.Function` argument receives an R
function, and its `call` method evaluates it, so the objective can be written
in R:

```{r 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
)
```

An R error in the objective stops the optimizer and reaches R unchanged, with
its class, after the Zig frames have run their `defer` blocks.

## Integration

`rzig.Integrate.integrate` integrates a Zig function as `integrate()` does,
also over infinite ranges:

```{r 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 and run

```{r 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)
```

```{r cleanup, include = FALSE}
if (has_zig && has_flibs) unloadNamespace("numericspkg")
unlink(path, recursive = TRUE)
```