---
title: Reading fmristore HDF5 files
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: interaction
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Reading fmristore HDF5 files}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("albersdown", quietly = TRUE)) {
  if (requireNamespace("systemfonts", quietly = TRUE)) try(albersdown::albers_register_fonts(), silent = TRUE)
  if (requireNamespace("ggplot2", quietly = TRUE)) try(ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset)), silent = TRUE)
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  eval = FALSE
)
```

```{r albers-classes, echo=FALSE, results='asis', eval=TRUE}
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
))
```

## When do you need this?

Many fMRI analysis pipelines store first-level results in HDF5 files using the
fmristore format. These files can contain voxelwise statistical maps (labeled
volumes), parcel-level timeseries, or latent-space decompositions. If you want
to run group-level analyses on data already stored in fmristore, fmrigds reads
it directly---auto-detecting the layout and converting to the GDS structure.

All examples below use `eval = FALSE` because they require HDF5 files that are
not shipped with the package. To run them locally, you need `hdf5r`:

```r
install.packages("hdf5r")
```

## Two paths for reading data

fmrigds supports two kinds of HDF5 input. The path you use depends on how the
file was created:

**Path B (recommended)**: Native `/gds` format
:   Files written by `write_out(..., format = "h5")` contain a `/gds` group
    with full metadata, provenance, and structure. These round-trip perfectly.

**Path A (legacy)**: Original fmristore layouts
:   Files created by the fmristore package use one of three layout types:
    LabeledVolumeSet, Parcellated, or LatentNeuroVec. fmrigds auto-detects
    the layout and converts to GDS.

If you control file creation, use Path B. If you're reading existing fmristore
files, Path A works automatically with `gds(path)`.

## Native `/gds` layout (Path B)

The simplest case---files with a `/gds` group are read directly with full
fidelity:

```r
library(fmrigds)

plan <- gds("group_store.h5")
g <- compute(plan)

# Everything is preserved as written
names(assays(g))
space(g)
metadata(g)$provenance
```

Create files in this format with `write_out()`:

```r
gds(source) |>
  reduce(method = "fixed") |>
  write_out("group_results.h5", format = "h5") |>
  compute()
```

## Labeled volumes (LabeledVolumeSet)

The most common fmristore layout for first-level outputs. It stores
statistical maps in voxel space with one dataset per contrast.

Structure inside the HDF5 file:

```
/mask         3D binary mask (which voxels have data)
/data/
  cope1       Packed vector of values for masked voxels
  cope2       ...
/labels       List of contrast names
```

```r
plan <- gds("subject01_lvs.h5")
g <- compute(plan)

space(g)               # voxel space
names(assays(g))       # "beta", "var", etc.
contrasts(g)           # e.g., c("cope1", "cope2")
dim(assay(g, "beta"))  # [n_voxels x n_subjects x n_contrasts]
```

Contrast names containing special characters are sanitized
(e.g., `"cope-1"` becomes `"cope_1"`).

## Parcellated data (H5ParcellatedScan)

Parcellated layouts store timeseries data for brain parcels (ROIs). Each file
contains:

```
/cluster_map    3D integer array (voxels -> parcel IDs)
/scans/
  001           Matrix: parcels x timepoints
  002           ...
```

The challenge: fmristore stores *timeseries*, but GDS needs *contrasts*. You
must tell fmrigds how to convert timepoints into contrasts using a
**temporal policy**.

### Temporal policies

**`"as_is"` (default)** --- treat each timepoint as a separate contrast:

```r
plan <- gds("parcellated.h5")
g <- compute(plan)
contrasts(g)  # c("t_001", "t_002", "t_003", ...)
```

Use this when you want to analyse the full timecourse.

**`"mean"`** --- average all timepoints into one contrast:

```r
plan <- gds("parcellated.h5", temporal_policy = "mean")
g <- compute(plan)
contrasts(g)  # "mean"
```

Only `beta` is averaged; `var`/`se` become `NA` (no uncertainty is
fabricated).

**`"design"`** --- apply a contrast matrix to extract specific effects:

```r
# Extract two contrasts from 3 timepoints
C <- cbind(early = c(1, 0, 0), late = c(0, 1, 0))

plan <- gds(
  "parcellated.h5",
  temporal_policy = "design",
  contrast_matrix = C,
  contrast_names = c("early", "late")
)

g <- compute(plan)
contrasts(g)           # c("early", "late")
dim(assay(g, "beta"))  # [n_parcels x n_subjects x 2]
```

Variance propagation: if the source file includes `/var`, fmrigds applies
`Var(C'y) = (C^2)' Var(y)` to propagate uncertainty through the linear
combination.

### Space type

After reading, the GDS has parcel-type space:

```r
g <- compute(gds("parcellated.h5"))
space(g)  # parcels space with labels
```

## Latent bases (LatentNeuroVec)

Latent layouts store dimensionality-reduced data where brain activity is
represented by a small number of spatial components (e.g., from PCA or ICA).

Structure:

```
/basis          Matrix: components x voxels
/scans/
  001           Embeddings: how much each component is expressed over time
  002           ...
```

```r
plan <- gds("latent_small.h5")
g <- compute(plan)

space(g)               # basis space
dim(assay(g, "beta"))  # [k_components x n_subjects x n_contrasts]
```

Each "sample" is a spatial component (a brain network). If the original data
had no uncertainty estimates, `var`/`se` will be `NA`.

## Reading subsets efficiently

For large datasets, read specific slices using the `block` parameter instead
of loading everything into memory:

```r
# Read only first 10 samples and first contrast
raw <- compute(
  plan,
  assays = "beta",
  block = list(sample = 1:10, contrast = 1)
)

# Read specific subjects
raw <- compute(
  plan,
  assays = c("beta", "var"),
  block = list(subject = c("sub-01", "sub-03"))
)
```

## Provenance tracking

Every `compute()` operation records what was done:

```r
plan <- gds("data.h5") |>
  reduce(method = "fixed")

g <- compute(plan)
metadata(g)$provenance$log
```

For native `/gds` files (Path B), provenance is preserved across
write--read cycles.

## Quick reference

| Scenario | What to do |
|---|---|
| Creating new files | `write_out(..., format = "h5")` (Path B) |
| Existing fmristore files | `gds(path)` --- layout auto-detected |
| Parcellated timeseries | Set `temporal_policy` (`"as_is"`, `"mean"`, or `"design"`) |
| Large files | Use `block = list(...)` to read subsets |
| Reproducibility | Provenance tracked automatically in `metadata()$provenance` |

## Next steps

- `vignette("fmrigds")` --- getting started with the core pipeline
- `vignette("spatial-operations")` --- masking, alignment, and space
  transformations
- `vignette("as-plan-and-spatial-fdr")` --- post-hoc corrections and spatial
  FDR
- `?gds` for the full list of adapter arguments