---
title: "Working without data copies"
output:
  rmarkdown::html_vignette:
    toc: true
    css: vignette_styles.css
bibliography: references.bib
csl: ieee.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{Working without data copies}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.align = "center",
                      fig.width = 6, fig.height = 4)
options(width = 80)
library(TRexSelector)
old_threads <- trex_set_num_threads(min(2L, trex_get_max_threads()))
```

# The rule

The backend works on the caller's data wherever possible. An R matrix is bound by its memory address,
a memory-mapped file by its mapping, and the preprocessing a method needs is applied in place and undone
afterwards. A design of several gigabytes is therefore held once, by R or by the file, and not a second
time by the backend. This page states the contract that makes it safe.

# Pass double matrices

R stores a matrix of whole numbers as integers. The backend binds double storage, so an integer matrix,
a data frame or a logical matrix is converted first, and the conversion is a copy by R's own rules. The
conversion happens visibly in R before the call, once, and a double matrix passes through untouched.
To keep a large design copy-free, create it as double.

```{r double}
Xi <- matrix(sample(0:2, 60 * 20, replace = TRUE), 60, 20)
storage.mode(Xi)
storage.mode(Xi) <- "double"     # one explicit conversion, then zero-copy from here on
storage.mode(Xi)
```

# Selectors: in place at construction, back at the end

A selector centers and scales `X` in place when it is constructed and writes the original values back
when `select()` returns. In between, `X` holds the working design. `restore()` writes back early for
a selector that will not run, and a selector that is garbage collected without running writes back in
its finalizer.

```{r selector}
set.seed(1)
n <- 80
p <- 40
X <- matrix(rnorm(n * p), n, p)
y <- as.numeric(X[, 1] * 2 + rnorm(n))
X_before <- X + 0
sel <- TRexSelector$new(X, y, tFDR = 0.1, seed = 2, verbose = FALSE, control = trex_control(K = 10))
round(sqrt(colSums(X^2))[1:3], 6)     # the working design: unit columns
sel$select()
max(abs(X - X_before))                 # the original values are back
```

# T-solvers: writers and readers

A T-solver with `normalize = TRUE` or `intercept = TRUE` is a writer. It transforms `X` and `D` in
place, keeps the statistics, and holds a claim on both matrices until `restore()` or garbage collection
writes them back. A second writer on a claimed matrix is refused, because a second transformation would
move the data by rounding.

```{r writers}
D <- matrix(rnorm(n * p), n, p)
w1 <- TLARS_Solver$new(X, D, y)
tryCatch(TLASSO_Solver$new(X, D, y), error = function(e) cat("error:", conditionMessage(e), "\n"))
w1$restore()
```

A solver with `normalize = FALSE` and `intercept = FALSE` is a reader. It runs on the data as it finds
them and never writes, so any number of readers can share one design. A pipeline scales the design once
with a scaler and constructs readers.

```{r readers}
sx <- MLLpNormScaler$new(norm_type = 2, center = TRUE)
invisible(sx$fit_transform_inplace(X))
sd <- MLLpNormScaler$new(norm_type = 2, center = TRUE)
invisible(sd$fit_transform_inplace(D))
y_c <- y - mean(y)
opts <- tsolver_options(normalize = FALSE, intercept = FALSE, precentered = TRUE)
r1 <- TLARS_Solver$new(X, D, y_c, options = opts)
r2 <- TLASSO_Solver$new(X, D, y_c, options = opts)
r1$execute_step(T_stop = 2)
r2$execute_step(T_stop = 2)
invisible(sx$inverse_transform_inplace(X))
invisible(sd$inverse_transform_inplace(D))
max(abs(X - X_before))
```

`y` is never written. The solver centers an internal copy of the response.

# Loading a saved solver

`load()` restores a solver into one constructed on the same data in the same state. Load into a reader
on pre-scaled data, or into a writer constructed with the same options on the same raw data, which then
keeps the duty to write back.

# Mapped files are never read into memory

A memory-mapped matrix is bound through the same interface as an R matrix. The backend touches the
pages it needs and the file keeps the data, so the working set of a run is the backend's own buffers,
never a copy of the design. The R-level accessors (`[` and `as.matrix()`) do copy the requested block
into R memory, by design.

# The low-memory engines

The ridge and elastic-net classes have a `low_memory` engine that avoids every allocation of the size of
`X`.

| Class | `low_memory = FALSE` | `low_memory = TRUE` |
| --- | --- | --- |
| `MLRidgeCV` | each fold gathers and standardizes its training rows and takes a thin SVD, about three to four copies of `X` in flight | the dual (kernel) form, $n \times n$ kernels accumulated from column blocks read in place, memory $O(n^2)$ |
| `MLElasticNet` | coordinate descent on a centered and standardized working copy, one $n \times p$ matrix | coordinate descent on `X` itself, no allocation of the size of `X` |
| `MLElasticNetCV` | one working copy per fit | every fold on `X` itself with a training mask |

The two engines compute the same estimator on the same grid and agree to rounding. The selectors that
tune $\lambda_2$ by cross-validation (GVS, SPCA, Tikhonov) use these engines, and the elastic-net tuner
of the sparse PCA selector keeps its peak memory under 0.3 times the size of `X`.

```{r low-memory}
cv_a <- MLRidgeCV$new()
cv_a$fit(X, y, n_folds = 5, n_lambda = 30)
cv_b <- MLRidgeCV$new(low_memory = TRUE)
cv_b$fit(X, y, n_folds = 5, n_lambda = 30)
all.equal(cv_a$get_cv_errors(), cv_b$get_cv_errors())
```

# Measuring the footprint

`get_memory_bytes()` reports what a T-solver holds on the heap. `X` and `D` are borrowed and not
counted.

```{r bytes}
r1$get_memory_bytes()
```

```{r teardown, include = FALSE}
trex_set_num_threads(old_threads)
```

# References