---
title: Large data
output:
  rmarkdown::html_vignette:
    css: albers.css
    includes:
      in_header:
      - albers-header.inc
      - albers-header.html
    toc: yes
    toc_depth: 2.0
resource_files:
- albers.css
- albers-fonts.css
- albers.js
- albers-header.inc
- albers-header.html
- fonts
params:
  family: red
  preset: interaction

vignette: |
  %\VignetteIndexEntry{Large data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("systemfonts", quietly = TRUE)) albersdown::albers_register_fonts()
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
source("_common.R")
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

```{r load}
library(neuroim2)
```

Dense is the default because it is the simplest thing that works, not because it
is the right thing at scale. An R array holds `double`s whatever the file said,
so a 56 MB `int16` run is 236 MB the moment you read it, and a session of runs
will not fit. The two out-of-core routes below are not a consolation prize for
running out of memory: for the scattered, voxel-at-a-time access that
searchlight, ROI and connectivity work actually perform, they are also
*faster* than loading everything first.

On an ordinary run (64 × 64 × 36 × 200), pulling 5,000 scattered voxel time
courses measured **44 ms** through `read_vec(path, mask = brain)` and **144 ms**
through `mode = "mmap"` including the open, against **4.8 s** for reading the
whole image and then indexing it. So:

- **Pass a `mask`** whenever the analysis is confined to one — which is nearly
  always. It is the cheapest option on both axes.
- **Use `mode = "mmap"`** when the access order is scattered and unknown ahead of
  time.
- **Stay dense** when you genuinely need every voxel of a small image in memory
  at once.

What follows measures each backend so the choice is yours rather than a
recommendation you have to take on faith. The backends differ enormously in
*which* access pattern they make cheap, and picking the wrong one is worse than
any of the defaults.

Everything below measures the same dataset, written to an uncompressed NIfTI so
that every backend can read it:

```{r data}
mask <- demo_mask()
brain <- as.mask(mask)
bold <- demo_bold(n_time = 60)

path <- tempfile(fileext = ".nii")
write_vec(bold, path)

dim(bold)
round(file.size(path) / 1e6, 1)
```

The file is uncompressed. Every backend can read `.nii.gz` too, but two of them
pay for it: memory mapping silently decompresses to a temporary uncompressed
copy, and file-backed access constructs happily and then fails at the first
read. Write uncompressed if you intend to map.

## What each backend costs

```{r backends}
dense <- read_vec(path)
sparse <- read_vec(path, mask = brain)
mapped <- read_vec(path, mode = "mmap")
filebacked <- read_vec(path, mode = "filebacked")

mb <- function(x) round(as.numeric(object.size(x)) / 1e6, 1)

data.frame(
  backend = c("DenseNeuroVec", "SparseNeuroVec", "MappedNeuroVec", "FileBackedNeuroVec"),
  memory_MB = c(mb(dense), mb(sparse), mb(mapped), mb(filebacked))
)
```

They all hold the same data:

```{r agreement}
ref <- as.numeric(series(dense, 32, 32, 12))

vapply(
  list(sparse = sparse, mapped = mapped, filebacked = filebacked),
  function(x) isTRUE(all.equal(as.numeric(series(x, 32, 32, 12)), ref)),
  logical(1)
)
```

## The number that actually decides it

Memory is the easy axis. The one that catches people is access *pattern*. Here
is the same work — twenty whole volumes, then twenty voxel time courses — on
each backend:

```{r timing}
per_call <- function(f, n) unname(system.time(for (i in seq_len(n)) f(i))["elapsed"] / n)

bench <- function(x, n = 20) {
  c(
    volumes = per_call(function(i) invisible(x[[i]]), n),
    series = per_call(function(i) series(x, 32, 32, 12), n)
  )
}

timings <- rbind(
  dense = bench(dense),
  sparse = bench(sparse),
  mapped = bench(mapped),
  filebacked = bench(filebacked, n = 5)
)

signif(timings, 2)
```

```{r timing-ratio}
round(timings["filebacked", "series"] / timings[c("dense", "sparse", "mapped"), "series"])
```

Seconds per call, and the ratio is computed rather than quoted, so it stays
honest on your machine. Fetching one voxel's time course from a
`FileBackedNeuroVec` costs three to four orders of magnitude more than it costs
anywhere else, because each query materialises the whole 4D
array as doubles before discarding all but one row — the disk read is a small
part of it. **Use it only when you genuinely work volume at a time.**
`MappedNeuroVec` is cheap on both patterns and is the better default.

## The ladder

**Dense** is what you get with no arguments. Predictable arithmetic, indexing
and plotting, and nothing to think about — at the cost of holding every voxel as
a `double`. Right for a single small image, wrong for a session.

**Sparse** stores only the voxels a mask admits. Here the brain is 29% of the
grid, and the object is about a third the size:

```{r sparse}
c(in_mask = sum(brain), total = prod(dim(brain)))
c(dense_MB = mb(dense), sparse_MB = mb(sparse))
```

Dimensions are unchanged — sparsity is storage, not shape. Reach for it when a
mask defines your analysis anyway.

**Memory-mapped** leaves the data on disk and lets the operating system page it
in. Pointed at a `.nii.gz` it decompresses to a temporary uncompressed copy
first, so you pay the disk either way — better to store uncompressed
deliberately than to have a full copy appear in `tempdir()`.

```{r mmap}
class(mapped)[1]
dim(mapped)
```

`as_mmap()` converts an in-memory object, writing an uncompressed file if it has
to:

```{r as-mmap}
mmap_path <- tempfile(fileext = ".nii")
converted <- as_mmap(dense, file = mmap_path, overwrite = TRUE)

class(converted)[1]
isTRUE(all.equal(as.numeric(series(converted, 32, 32, 12)), ref))
```

Give it an explicit `file` when you can. Left to itself it writes into
`tempdir()` and you have a second full copy of the data with no obvious owner.

**BigNeuroVec** backs a masked series with an on-disk `bigstatsr` matrix, which
is the option for data larger than RAM that still needs voxelwise access:

```{r bigvec}
big <- read_vec(path, mode = "bigvec", mask = brain)

class(big)[1]
round(mb(big), 2)
```

### What object.size() does not see

The two on-disk backends report near-zero above, and that is `object.size()`
telling you the *R object* is a handle. It is not telling you the data is free:
mapped pages are charged to the process as you touch them, and `BigNeuroVec`
keeps a backing file — about 14 MB here — beside the original. Measure resident
memory and disk when the difference matters.

## Dropping to parcels

Every backend so far keeps one time course per voxel. The largest reduction
available comes from deciding you do not need them: a `ClusteredNeuroVec` stores
one series per *parcel*.

Parcels have to be spatially contiguous for any of this to mean anything.
Assigning labels at random would put all 100 parcels everywhere at once and
collapse their centroids onto the same point, so cluster them by position:

```{r clustered}
set.seed(1)
coords_in <- index_to_grid(mask, which(as.vector(brain)))
km <- kmeans(coords_in, centers = 100, iter.max = 30)

parcels <- ClusteredNeuroVol(brain, km$cluster)
cvec <- ClusteredNeuroVec(bold, parcels)

num_clusters(cvec)
dim(cvec)
round(apply(centroids(cvec), 2, sd), 1)
```

Those centroid standard deviations are voxels: the parcels are spread across the
brain, which is what makes centroid distances meaningful below.

It still behaves like a 4D image — `dim()` is unchanged and indexing works — but
only 100 time courses exist behind it:

```{r clustered-size}
dim(as.matrix(cvec, by = "cluster"))

c(dense_MB = mb(dense), clustered_MB = mb(cvec))
```

```{r clustered-access}
length(series(cvec, 32, 32, 12))
dim(cvec[, , , 1])

c(voxel_series = sum(brain), stored_series = num_clusters(cvec),
  fold = round(sum(brain) / num_clusters(cvec)))
```

Every voxel in a parcel returns that parcel's series: roughly a 300-fold
reduction in stored series here, and no within-parcel detail at all. For
atlas-based connectivity that is exactly right; for anything needing voxelwise
variation it is not.

One sharp edge: voxels outside the mask have no parcel, and `series()` returns
`NA` for them rather than zero or an error. Restrict to the mask before doing
arithmetic over the whole grid.

Cluster-level searchlights come with it, using centroid distances rather than
voxel neighbourhoods:

```{r cluster-searchlight}
windows <- cluster_searchlight_series(cvec, k = 3)

length(windows)
dim(values(windows[[1]]))
```

```{r cleanup, include = FALSE}
unlink(c(path, mmap_path, big@data$backingfile))
```

## Choosing

| Situation | Backend | Why |
|:--|:--|:--|
| It fits in memory | `DenseNeuroVec` | simplest, fastest, no caveats |
| A mask defines the analysis | `SparseNeuroVec` | stores only admitted voxels |
| Large uncompressed file, mixed access | `MappedNeuroVec` | cheap on both access patterns |
| Large file, strictly volume at a time | `FileBackedNeuroVec` | never for voxel series |
| Larger than RAM, needs voxel access | `BigNeuroVec` | on-disk matrix |
| Parcel-level analysis is enough | `ClusteredNeuroVec` | orders of magnitude fewer series |

If you are unsure: read with `read_vec()`, and change only when a measurement
tells you to. `object.size()` and `system.time()` will get you most of the way,
as long as you remember that the first does not see mapped pages or backing
files.

## Where to go next

- `vignette("volumes-and-vectors")` — the containers these all implement
- `vignette("regions-and-searchlights")` — parcels and searchlights in depth
- `?read_vec`, `?as_mmap`, `?ClusteredNeuroVec`, `?BigNeuroVec`