---
title: "Detailed Documentation"
author: "Konrad Kraemer"
output: html_document
# pdf_document
#classoption:
#  - twocolumn
# rmarkdown::html_vignette

vignette: >
  %\VignetteIndexEntry{Detailed Documentation}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
```

```{css echo=FALSE}
/* Define a margin before h2 element */
h2  {
  margin-top: 6em;
}

/* Define a margin after every first p elements */
p:first-of-type {
  margin-bottom: 3em;
}
```

<style>
body {
text-align: justify}
</style>

## Overview


The *ast2ast* package translates R functions into C++ functions,
returning either an external pointer (XPtr) or an R function.
This package is particularly useful for tasks requiring frequent function evaluations,
such as solving ODE systems or optimization problems.
Using the external pointer generated by C++ can significantly enhance performance,
as shown in the benchmark below.

![Benchmark](benchmark.png)

Supported objects:

*  scalar values
*  vectors
*  matrices
*  arrays
*  custom types (structs) and collections thereof, see `new_type` below

Supported functions:

*  assignment: `=`, `<-`, and `$` (struct field access/assignment)
*  allocation: `vector`, `logical`, `integer`, `numeric`, `:`, `c`, `rep`, `matrix`, and `array`
*  some helpers: `length`, `dim`, `nrow`, `ncol`, `is.na`, `is.nan`, `is.finite`, `is.infinite`, `print`, and `stop`
*  Math functions: `+`, `-`, `*`, `/`, `%%`, `%/%`, `%*%`, `sin`, `asin`, `sinh`, `cos`, `acos`, `cosh`, `tan`, `atan`, `tanh`, `log`, `sqrt`, `abs`, `sign`, `^`, `exp`, `floor`, `ceiling`, `trunc`, and `round`
*  reductions: `sum`, `prod`, `mean`, `max`, `min`, `which.max`, `which.min`, `which`, `all`, `any`
*  vector transforms: `rev`, `sort`, `cumsum`
*  matrix margins: `colSums`, `rowSums`, `colMeans`, `rowMeans`
*  selection: `ifelse`
*  linear algebra: `t`, `chol`, `crossprod`, `tcrossprod`, `diag`, `get_diag`, `solve`, `backsolve`, `forwardsolve`, `rbind`, `cbind`
*  root finding / least squares: `uniroot`, `nnls`
*  functionals (take an `fn`): `map`, `Reduce`, `Filter`, `apply`
*  optimizers / Jacobians of an `fn`: `jacobian`, `lbfgsb`, `pso`
*  casts: `as.numeric`, `as.integer`, `as.logical`
*  sequences: `seq_len`, `seq_along`
*  Subsetting: `[]`, `[[]]`, and `at()` (alias for `[[]]`)
*  control flow: `for`, `while`, `repeat`, `if`, `else if`, and `else`
*  comparison: `==`, `!=`, `>`, `<`, `>=`, and `<=`
*  logical: `&&`, `||`, `&`, `|`, `!`
*  derivative: `seed`, `unseed`, `get_dot`, `deriv`
*  inner functions: `fn`, see below

### Type system in ast2ast

### Overview

R is dynamically typed, while C++ is statically typed. When translating an R function, **ast2ast** must decide static C++ types for:

1. the **arguments** of `f`, and
2. the **variables created inside** `f`.

In **ast2ast**, every type is a combination of:

- a **base type**: `logical`, `integer` (or `int`), `double`
- a **structure**: `scalar`, `vector` / `vec`, `matrix` / `mat`

Typical types are therefore:

- `double` (scalar double)
- `vec(double)` / `vector(double)` (double vector)
- `mat(double)` / `matrix(double)` (double matrix)

A key difference to R: **ast2ast scalars are true scalars**, not length-1 vectors. Scalars cannot be subset using `[]` or `[[ ]]`.

Another difference to R: **negative indices are not supported** (R's `v[-1L]`
drop-element form). Indices must be positive and 1-based.

Another difference to R: **variables are declared and zero-initialised up
front** (Fortran-style: all declarations first, then the body). A variable that
is only assigned on one branch but read unconditionally therefore does *not*
raise an "object not found" error the way R would -- it reads its zero value.
ast2ast does not currently perform use-before-assignment analysis, so make sure
every path that reads a variable also assigns it.

```{r, eval = FALSE}
f <- function(a) {
  if (a > 0) x <- 1.0   # x is only assigned here...
  return(x)             # ...but read here regardless: returns 0.0 when a <= 0
}
```

---

## Setting types for function arguments

### Default behavior

If no `argtypes(...)` block is present, all arguments default to `matrix(double)`, as this is most convinient for numeric code.

---

### Using `argtypes()`

To control argument types, make `argtypes(...)` the **first statement** of `f`'s body. Each entry names one of `f`'s arguments and assigns its type with `type()`:

```{r, eval = FALSE}
f <- function(a, b, c) {
  argtypes(
    a |> type(vec(double)),
    b |> type(mat(double)),
    c |> type(double)
  )
  # ... body ...
}
f_cpp <- ast2ast::translate(f)
```

---

### Borrowing, constness, and references

For arguments, you can additionally control how values are passed:

- `borrow_vec(...)`, `borrow_mat(...)`: borrow memory (no copy)
- `const()`: disallow modification
- `ref()`: pass by reference (only valid when `output = "XPtr"`)

Example:

```{r, eval = FALSE}
f <- function(a, b, c) {
  argtypes(
    a |> type(borrow_vec(double)) |> ref(),            # mutable, passed by reference
    b |> type(borrow_mat(double)) |> ref() |> const(), # read-only matrix reference
    c |> type(double) |> ref()                         # scalar reference (XPtr only)
  )
  # ... body ...
}
```

Notes:

- Borrowed arguments are useful for avoiding allocations and modifying inputs in place.
- `const()` is enforced by the static checker at translation time ("You cannot assign to a constant variable"), before any C++ compilation happens.
- `ref()` is primarily intended for the external pointer interface.

---

### Setting types for variables inside the function

Types inside `f` are often inferred automatically from:

- the first assignment, and/or
- the constructor used (`numeric()`, `matrix()`, `integer()`, etc.)

You can override inference using explicit annotations.

---

### Static types vs. R's dynamic typing

In R, reassigning a variable to a different type mid-function is completely normal:

```{r, eval = FALSE}
a <- 1L
a <- 2.5
```

The generated C++ can't do that. A variable has exactly **one** type for its entire lifetime,
declared once, before any of the function's logic runs. So when a variable's type isn't pinned
down with an explicit `type()` annotation, ast2ast has to *decide* that one type -- and it does so
by looking at **every** assignment to that variable across the whole function, not just the first
one, and picking the type wide enough to hold all of them. That decision is then applied
retroactively to *all* of the variable's assignments -- including ones that come before the line
that actually forced the wider type:

```{r, eval = FALSE}
f <- function() {
  a <- 1L      # looks like an integer assignment...
  a <- 2.5     # ...but a is reassigned a double two lines later
  return(a)
}
```

Because `a` is also assigned a `double` further down, ast2ast decides `a` is `double` for the
*whole* function. The generated C++ declares `a` as `double` up front, and even the first
assignment (`a <- 1L`) compiles to storing `1.0` -- not `1L`. A warning is raised whenever this
happens (`Promoted the type of variable a from ... to ...`), so it's visible rather than silent.

This "pick the wider of two types" logic is the same one used to combine two operands within a
single expression (e.g. `integer + double`), along two independent precedence orders:

- base type: `logical` < `integer` (`int`) < `double`
- structure: `scalar` < `vector` (`vec`) < `matrix` (`mat`) < `array`

The result takes the higher-precedence base type and the higher-precedence structure of the two
operands (or, for a variable, of everything ever assigned to it), e.g.:

- `logical + integer` -> `integer`
- `integer + double` -> `double`
- `double` (scalar) `+ vec(double)` -> `vec(double)`

It only ever widens -- never narrows, and never reconciles genuinely incompatible kinds (e.g. a
custom type and a plain vector; those are hard errors, not promotions).

A few operators promote unconditionally, regardless of what they're given:

- `/` and `^` always produce `double` (R never returns an integer from division or
  exponentiation)
- transcendental functions (`sin`, `sqrt`, `log`, `exp`, ...) always return `double`
- `sum()` keeps `double` as `double`, but promotes `logical`/`integer` to `integer`
  (matches R's own `sum()` behavior)

**This promotion mechanism only applies to inferred variables.** A variable whose type was fixed
explicitly -- via a `type()` annotation, or as a function argument -- is truly immutable: it does
not widen, and assigning it something outside that fixed type is a translation-time error, not a
promotion. If you need a genuinely different type or structure somewhere in the function, create a
new variable with a new name rather than relying on an existing one to change shape.

## Derivatives

The ast2ast package provides built-in support for automatic differentiation (AD) in both forward mode and reverse mode.
Derivative support is enabled when translating a function via the derivative argument:

```{r, eval = FALSE}
fcpp <- ast2ast::translate(f, derivative = "forward")
fcpp <- ast2ast::translate(f, derivative = "reverse")
```

Unlike many high-level AD frameworks, ast2ast intentionally exposes a low-level and explicit interface.
Derivative computations are assembled from a small set of primitive operations (`seed`, `unseed`,
`get_dot`, `deriv`), which keeps the behavior transparent, predictable, and close to the generated
C++ code. A built-in `jacobian(f, x)` wrapper is also available for the common case (see
*Optimizers* below); it is implemented on top of the same primitives.

### Forward mode

In forward mode, derivatives are propagated alongside values. Internally, each scalar carries both its value and its directional derivative (also called its dot value).
The following functions are available:

* seed(x, i): Activates the i-th component of x as the differentiation direction (sets its derivative to 1).
* unseed(x, i): Resets the derivative state of the i-th component.
* get_dot(y): Extracts the directional derivatives of y.

A typical pattern is to compute Jacobians column-by-column by looping over the input variables:

```{r, eval = FALSE}
f <- function(y, x) {
  jac <- matrix(0.0, length(y), length(x))
  for (i in 1L:length(x)) {
    seed(x, i)

    y[[1L]] <- x[[1L]] * x[[2L]]
    y[[2L]] <- x[[1L]] + x[[2L]] * x[[2L]]

    d <- get_dot(y)
    jac[TRUE, i] <- d

    unseed(x, i)
  }
  return(jac)
}

fcpp_forward <- ast2ast::translate(f, derivative = "forward")
```

Forward mode is most efficient when the number of inputs is small relative to the number of outputs.

### Reverse mode

In reverse mode, derivatives are accumulated by propagating sensitivities backward from
the outputs to the inputs. This is particularly efficient when the number of outputs is small relative to the number of inputs.

Reverse mode provides the function:

* deriv(y, x): Computes the Jacobian of y with respect to x.

Example:

```{r, eval = FALSE}
f <- function(y, x) {
  y[[1L]] <- x[[1L]] * x[[2L]]
  y[[2L]] <- x[[1L]] + x[[2L]] * x[[2L]]
  jac <- deriv(y, x)
  return(jac)
}

fcpp_reverse <- ast2ast::translate(f, derivative = "reverse")
```
The call to deriv() must appear explicitly in your function body. No automatic differentiation is performed unless requested.

### Design philosophy

Derivative computation in ast2ast is explicit by design.
The full control flow—loops, seeding, unseeding, derivative extraction, and accumulation—is written directly in R and translated into C++.

This approach:
* avoids hidden performance costs,
* makes derivative logic easy to inspect and debug,
* gives full control over memory and evaluation order,
* maps naturally to high-performance C++ code.

Rather than hiding differentiation behind abstractions, ast2ast treats derivatives as first-class values that can be manipulated like any other object.

## Inner functions

Functions can be defined inside `f` using `fn()`. This is required whenever a function
needs to be passed as a value, e.g. to `uniroot()`. An inner function is declared with
three positional parts: `argtypes(...)` (types of its arguments, same form as the outer
`argtypes(...)` block), `return(...)` (its return type), and a `{ }` block (the R code).
A single-statement body does not need the `{ }`, but wrapping it is fine too.

```{r, eval = FALSE}
f <- function(a) {
  argtypes(a |> type(int))
  factorial <- fn(
    argtypes(a |> type(int) |> const()),
    return(int),
    {
      if (a == 1L) return(a) else return(a * factorial(a - 1L))
    }
  )
  return(factorial(a))
}
fcpp <- ast2ast::translate(f)
```

An inner function's non-`const` parameters only bind to bare variables passed at the
call site, not to arbitrary expressions (e.g. `x + 1`, `x[[1L]]`) -- a non-`const`
parameter is a mutable reference to the caller's argument, and an expression has no
addressable storage to reference. Declare the parameter `const()` if you need to pass
an expression:

```{r, eval = FALSE}
sq <- fn(
  argtypes(x |> type(double) |> const()),  # const -- accepts expressions
  return(double),
  return(x * x)
)
# sq(a + b) is fine; without const() on x, only sq(a) (a bare variable) would be.
```

Inner functions can call each other (including mutual recursion) and can be passed to
functions expecting a function argument, e.g. `uniroot`:

```{r, eval = FALSE}
f <- function(interval) {
  argtypes(interval |> type(vec(double)))
  g <- fn(
    argtypes(x |> type(double)),
    return(double),
    {
      return(x^2 - 4)
    }
  )
  res <- uniroot(g, interval, 1e-10, 1000)
  return(res$root)
}
fcpp <- ast2ast::translate(f)
```

`uniroot(f, interval, tol, maxiter)` returns a struct with fields `root`, `f_root`,
`iter`, and `estim_prec` (accessed via `$`, see custom types below). `f` must take a
single `double` and return a `double`; an optional fifth argument is passed through to
`f` as extra data (`f` then takes two arguments).

`nnls(A, b)` solves the non-negative least squares problem and returns the solution
vector directly.

## Functionals

These take an inner function (`fn`) as their first argument.

* `map(f, x, ...)` -- apply `f` element-wise over the given vectors; scalars
  broadcast. The result shape follows `f`'s return type: scalar produces a
  vector, vector a matrix, matrix/array an array (with `n` as the last axis),
  and a `new_type` a `collection`.
* `Reduce(f, x)` -- left fold seeded with `x[[1]]`; `f` is `fn(acc, elem) -> acc`
  and must return the same type as its accumulator argument.
* `Filter(f, x)` -- keep the elements of vector `x` for which `f(elem)` is `TRUE`
  (`f` returns a logical scalar).
* `apply(f, MARGIN, x)` -- function first, unlike base R's
  `apply(X, MARGIN, FUN)`. `x` is a matrix, `MARGIN` is `1` (rows) or `2`
  (columns), and `f` maps a vector to a scalar or a vector.

```{r, eval = FALSE}
f <- function(x) {
  argtypes(
    x |> type(vec(double))
  )
  sq <- fn(
    argtypes(
      a |> type(double) |> const()
    ),
    return(double),
    return(a * a)
  )
  return(map(sq, x))
}
fcpp <- ast2ast::translate(f)
fcpp(1:5)
```

## Optimizers

* `jacobian(f, x[, data])` -- the `m`-by-`n` Jacobian of `f` (double vector to
  double vector) evaluated at `x`. Requires `derivative = "forward"` or
  `"reverse"`.
* `lbfgsb(f, x, lower, upper, maxit, factr, pgtol, lmm[, data])` --
  bound-constrained L-BFGS-B via R's own C routine. `f` maps a double vector to a
  scalar double. Under `derivative = "forward"` / `"reverse"` the gradient is
  exact (computed with `jacobian`); otherwise it is a central-difference
  approximation. Returns a struct with `$par`, `$value`, `$convergence`,
  `$counts`.
* `pso(f, lower, upper, ngen, npop, error_threshold, global[, data])` --
  derivative-free particle-swarm optimisation. `f` maps a double vector to a
  scalar double; returns the best parameter vector found.

For `lbfgsb` and `pso` the optional trailing `data` argument (any non-function,
non-character value) is passed to `f` unchanged as a second argument, so `f` then
takes two arguments.

```{r, eval = FALSE}
rosen <- function(p) {
  argtypes(
    p |> type(vec(double))
  )
  loss <- fn(
    argtypes(
      x |> type(vec(double)) |> const()
    ),
    return(double),
    {
      a <- 1.0 - x[[1L]]
      b <- x[[2L]] - x[[1L]] * x[[1L]]
      return(a * a + 100.0 * b * b)
    }
  )
  lo <- c(-5.0, -5.0)
  up <- c(5.0, 5.0)
  res <- lbfgsb(loss, p, lo, up, 100L, 1e7, 1e-8, 5L)
  return(res$par)
}
fcpp <- ast2ast::translate(rosen, derivative = "reverse")
fcpp(c(-1.2, 1.0))
```

## Custom types (`new_type`)

Besides scalars, vectors, matrices and arrays, **ast2ast** supports user-defined struct
types. They are declared in a `types_f` helper function passed to `translate()`:

```{r, eval = FALSE}
types_f <- function() {
  new_type(Point, slots(x |> type(double), y |> type(double)))
}

f <- function(p) {
  argtypes(p |> type(Point))
  p$x <- p$x + 1
  return(p)
}

fcpp <- ast2ast::translate(f, types_f = types_f)

p <- structure(list(x = 1, y = 2), class = "Point")
fcpp(p)
```

Notes:

- On the R side, a value of a custom type is a named `list` with a matching `class`
  attribute, e.g. `structure(list(x = 1, y = 2), class = "Point")`.
- Slots can themselves be custom types, allowing nested structs.
- `collection(TypeName)` is a vector of a custom type, e.g.
  `slots(points |> type(collection(Point)))` or, as a local variable,
  `coll |> type(collection(Point))`. Allocate one with
  `vector(mode = "Point", n)`.
- Fields are read and written with `$`, including chained access such as
  `s$circles[[1L]]$center$x`.

## Interpolation

To interpolate values, the 'cmr' function can be used. The function needs three arguments.

* the first argument is the point of the independent variable (x) for which the dependent variable should be calculated (y). This has to be a scalar (a length-one vector triggers a warning and only uses its first element).
* the second argument is a vector defining the points of the independent variable (x). This has to be a vector of at least length four.
* the third argument is a vector defining the points of the dependent variable (y). This has to be a vector of at least length four.

```{r, eval = FALSE}
f <- function() {
  dep <- c(0, 1, 0.5, 2.5, 3.5, 4.5, 4)
  indep <- 1:7
  evalpoints <- c(
    0.5, 1, 1.5, 2, 2.5,
    3, 3.5, 4, 4.5, 5,
    5.5, 6, 6.5
  )
  for (i in evalpoints) {
    print(cmr(i, indep, dep))
  }
}
```