---
title: 'Shared-Basis HRF Matching (SBHM): Efficient Voxel-Specific HRF Estimation'
output:
  rmarkdown::html_vignette:
    mathjax: default
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html

vignette: |
  %\VignetteIndexEntry{Shared-Basis HRF Matching (SBHM): Efficient Voxel-Specific HRF Estimation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, include = FALSE}
if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("albersdown", quietly = TRUE)) {
  ggplot2::theme_set(
    albersdown::theme_albers(family = params$family, preset = params$preset)
  )
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  dev = "png",
  dpi = 150,
  warning = FALSE,
  message = FALSE
)

# Make vignette render in pkgdown by using a lightweight fast mode
fast_mode <- identical(Sys.getenv("IN_PKGDOWN"), "true") || identical(Sys.getenv("CI"), "true")

# Always evaluate chunks, but keep computations small in fast mode
eval_chunks <- TRUE

# Global settings used throughout
TR <- 1.0
n_time <- if (fast_mode) 160 else 200
set.seed(123)
```

```{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-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

## The Problem

Real fMRI data shows substantial HRF variability across brain regions. A single
canonical HRF often fits poorly, but estimating fully unconstrained voxel-wise
HRFs is expensive and noisy. SBHM offers a middle path: learn a shared basis
from a library of plausible HRFs, then match each voxel to its best library
member in a low-dimensional space.

The pipeline has four steps:

1. **Build** a shared basis from an HRF library via SVD (`sbhm_build()`)
2. **Prepass** aggregate model fit per voxel (`sbhm_prepass()`)
3. **Match** voxels to library members by cosine similarity (`sbhm_match()`)
4. **Project** trial-wise coefficients to scalar amplitudes (`sbhm_project()`)

In practice, `lss_sbhm()` runs all four steps in a single call.

```{r setup, message=FALSE, warning=FALSE}
library(fmrihrf)
library(fmrilss)

# Derived convenience defaults for examples
n_voxels_default <- if (fast_mode) 20 else 50
n_trials_default <- if (fast_mode) 6 else 10
ranks_default    <- if (fast_mode) 2:6 else 2:10
```

## Step 1: Build the Shared Basis

We start by defining a grid of gamma HRF parameters spanning physiological
variability in peak time and width.

```{r param-grid, eval=eval_chunks}
shapes <- if (fast_mode) seq(6, 10, by = 2) else seq(5, 11, by = 1.5)
rates  <- if (fast_mode) seq(0.8, 1.2, by = 0.2) else seq(0.7, 1.3, by = 0.15)
param_grid <- expand.grid(shape = shapes, rate = rates)
cat("Library size:", nrow(param_grid), "HRFs\n")
```

Each grid row becomes a gamma HRF. We wrap this in a factory function that
`sbhm_build()` can evaluate across the grid.

```{r gamma-fun, eval=eval_chunks}
gamma_fun <- function(shape, rate) {
  f <- function(t) fmrihrf::hrf_gamma(t, shape = shape, rate = rate)
  fmrihrf::as_hrf(f, name = sprintf("gamma(s=%.2f,r=%.2f)", shape, rate), span = 32)
}
```

Now we build the SBHM basis. The SVD decomposes the library into a shared time
basis `B`, singular values `S`, and per-HRF coordinates `A`.

```{r build-sbhm, eval=eval_chunks}
sframe <- sampling_frame(blocklens = n_time, TR = TR)
sbhm <- sbhm_build(
  library_spec = list(fun = gamma_fun, pgrid = param_grid, span = 32,
                      precision = 0.1, method = "conv"),
  r = 6, sframe = sframe,
  baseline = c(0, 0.5), normalize = TRUE, ref = "mean"
)
```

```{r report-dims, eval=eval_chunks}
cat("B (time basis):", dim(sbhm$B), "\n")
cat("A (library coords):", dim(sbhm$A), "\n")
cat("Singular values:", round(sbhm$S, 2), "\n")
```

### Visualizing the Basis

Each basis function captures a principal mode of HRF variation: the main shape,
timing shifts, width differences, and undershoot features.

```{r visualize-basis, fig.width=10, fig.height=6, eval=eval_chunks, fig.alt="Six panels showing SBHM shared basis functions capturing principal modes of HRF variation."}
par(mfrow = c(2, 3), mar = c(3, 3, 2, 1))
for (i in 1:ncol(sbhm$B)) {
  plot(sbhm$tgrid, sbhm$B[, i], type = "l", col = "navy", lwd = 2,
       main = paste0("Basis ", i, " (s=", round(sbhm$S[i], 2), ")"),
       xlab = "Time (s)", ylab = "Amplitude")
  abline(h = 0, col = "gray", lty = 2)
}
```

### Choosing the Rank

Aim for 90--95% variance explained. We can sweep ranks to find the elbow.

```{r choose-rank, fig.width=8, fig.height=5, eval=eval_chunks}
ranks <- ranks_default
ve <- sapply(ranks, function(ri) {
  sum(sbhm_build(library_spec = list(fun = gamma_fun, pgrid = param_grid, span = 32),
                 r = ri, sframe = sframe, normalize = TRUE)$S^2)
})
```

```{r plot-rank, fig.width=8, fig.height=5, eval=eval_chunks, fig.alt="Line plot of variance explained vs rank showing elbow around rank 6."}
plot(ranks, ve / max(ve) * 100, type = "b", pch = 19, col = "navy",
     lwd = 2, xlab = "Rank (r)", ylab = "Variance Explained (%)",
     main = "Choosing SBHM Rank")
abline(h = 95, col = "red", lty = 2)
grid()
```

### Library Coverage

A sample of rank-r reconstructed HRFs shows the range of shapes the library
spans.

```{r library-coverage, fig.width=8, fig.height=5, eval=eval_chunks, fig.alt="Overlaid HRF curves showing the range of shapes spanned by the library reconstruction."}
H_hat <- sbhm$B %*% sbhm$A
sel <- unique(round(seq(1, ncol(H_hat), length.out = min(ncol(H_hat), 12))))
matplot(sbhm$tgrid, H_hat[, sel, drop = FALSE], type = "l", lty = 1,
        col = colorRampPalette(c("#6baed6", "#08519c"))(length(sel)),
        lwd = 1.5, xlab = "Time (s)", ylab = "Amplitude",
        main = "Sample of Library HRFs (rank-r reconstruction)")
abline(h = 0, col = "gray80", lty = 2)
```

## Step 2: Generate Synthetic Data

To demonstrate recovery, we create data where each voxel uses a known HRF from
the library. First, set up the experimental design.

```{r design-setup, eval=eval_chunks}
n_voxels <- n_voxels_default
n_trials <- n_trials_default
safe_end <- max(sbhm$tgrid) - 30
onsets <- seq(20, safe_end, length.out = n_trials)
```

Assign each voxel a random library HRF and build per-trial regressors. The
fresh `set.seed()` here keeps the HRF assignment reproducible independently of
earlier random draws.

```{r assign-hrfs, eval=eval_chunks}
set.seed(456)
true_hrf_idx <- sample(ncol(sbhm$A), n_voxels, replace = TRUE)
design_spec <- list(
  sframe = sframe,
  cond = list(onsets = onsets, duration = 0, span = 30)
)
hrf_B <- sbhm_hrf(sbhm$B, sbhm$tgrid, sbhm$span)
```

```{r build-regressors, eval=eval_chunks}
regressors_by_trial <- lapply(onsets, function(ot) {
  rr_t <- regressor(onsets = ot, hrf = hrf_B, duration = 0, span = 30, summate = FALSE)
  evaluate(rr_t, grid = sbhm$tgrid, precision = 0.1, method = "conv")
})
```

Now generate the signal. Each voxel's response is the sum of per-trial
regressors projected through that voxel's HRF coordinates, scaled by a random
amplitude, plus noise.

```{r generate-signal, eval=eval_chunks}
Y <- matrix(rnorm(n_time * n_voxels, sd = 0.5), n_time, n_voxels)
true_amplitudes <- matrix(rnorm(n_trials * n_voxels, mean = 2, sd = 0.5),
                          n_trials, n_voxels)
for (v in 1:n_voxels) {
  alpha_true <- sbhm$A[, true_hrf_idx[v]]
  for (t in 1:n_trials)
    Y[, v] <- Y[, v] + true_amplitudes[t, v] * (regressors_by_trial[[t]] %*% alpha_true)
}
```

```{r data-summary, eval=eval_chunks}
cat("Y:", dim(Y), "\n")
cat("Unique true HRFs:", length(unique(true_hrf_idx)), "\n")
```

## Step 3: Run the SBHM Pipeline

A single call to `lss_sbhm()` runs the entire pipeline. All control lists
(`prepass`, `match`, `oasis`, `amplitude`) are optional; pass only what you want
to override. See `?lss_sbhm` for full parameter documentation.

```{r run-sbhm, eval=eval_chunks}
res_sbhm <- lss_sbhm(
  Y = Y, sbhm = sbhm, design_spec = design_spec,
  return = "both"
)
```

```{r sbhm-dims, eval=eval_chunks}
cat("Amplitudes:", dim(res_sbhm$amplitude), "\n")
cat("Coefficients:", dim(res_sbhm$coeffs_r), "\n")
cat("Matched indices:", length(res_sbhm$matched_idx), "\n")
```

If you use `fmridesign`, prefer `lss_sbhm_design()` for an even simpler
interface that mirrors `lss_design()`. See `?lss_sbhm_design`.

## Step 4: Evaluate Recovery

### Matching Accuracy

```{r evaluate-recovery, eval=eval_chunks}
accuracy <- mean(res_sbhm$matched_idx == true_hrf_idx)
cat("HRF matching accuracy:", round(100 * accuracy, 1), "%\n")
cat("Confused voxels:", sum(res_sbhm$matched_idx != true_hrf_idx),
    "/", n_voxels, "\n")
```

### Matching Confidence

The `margin` (top-1 minus top-2 cosine score) indicates how unambiguous each
assignment is. Higher is better.

```{r margin-stats, eval=eval_chunks}
cat("Margin -- mean:", round(mean(res_sbhm$margin), 3),
    " median:", round(median(res_sbhm$margin), 3), "\n")
```

```{r visualize-margin, fig.width=8, fig.height=5, eval=eval_chunks, fig.alt="Histogram of matching confidence margins with median line."}
hist(res_sbhm$margin, breaks = 20, col = "skyblue", border = "white",
     main = "Matching Confidence (Margin)",
     xlab = "Margin (top1 - top2 cosine score)")
abline(v = median(res_sbhm$margin), col = "red", lwd = 2, lty = 2)
grid()
```

In this simulation the median margin is small, which means many voxels sit close
to multiple library members. That is exactly the regime where soft blending (see
below) can help.

### Amplitude Recovery

```{r compare-amplitudes, fig.width=8, fig.height=5, eval=eval_chunks, fig.alt="Scatter plot of true vs estimated amplitudes showing recovery quality."}
cor_amp <- cor(as.vector(res_sbhm$amplitude), as.vector(true_amplitudes))
plot(as.vector(true_amplitudes), as.vector(res_sbhm$amplitude),
     pch = 19, col = adjustcolor("navy", alpha.f = 0.3),
     xlab = "True Amplitude", ylab = "Estimated Amplitude",
     main = paste0("Amplitude Recovery (r = ", round(cor_amp, 3), ")"))
abline(0, 1, col = "red", lwd = 2, lty = 2)
grid()
```

```{r recovery-summary, eval=eval_chunks}
recovery_summary <- data.frame(
  HRFMatchingAccuracy = accuracy,
  AmplitudeCorrelation = cor_amp,
  MedianMargin = median(res_sbhm$margin)
)
recovery_summary
```

```{r recovery-summary-check, include = FALSE, eval = eval_chunks}
stopifnot(all(is.finite(as.matrix(recovery_summary))))
```

### Recovered HRF Shapes

For the most confidently matched voxels, the recovered and true HRFs should
overlap closely.

```{r hrf-prep, eval=eval_chunks}
H_hat <- sbhm$B %*% sbhm$A
vox_show <- head(order(-res_sbhm$margin), n = min(6, n_voxels))
```

```{r matched-vs-true, fig.width=10, fig.height=6, eval=eval_chunks, fig.alt="Six panels comparing true and matched HRF shapes for the most confident voxels."}
par(mfrow = c(2, 3), mar = c(3, 3, 2, 1))
for (v in vox_show) {
  rng <- range(c(H_hat[, true_hrf_idx[v]], H_hat[, res_sbhm$matched_idx[v]]))
  plot(sbhm$tgrid, H_hat[, true_hrf_idx[v]], type = "l", col = "#2c7fb8",
       lwd = 2, ylim = rng, main = paste0("Voxel ", v), xlab = "Time (s)", ylab = "HRF")
  lines(sbhm$tgrid, H_hat[, res_sbhm$matched_idx[v]], col = "#d95f02", lwd = 2, lty = 2)
  legend("topright", bty = "n", cex = 0.8, legend = c("True", "Matched"),
         lty = 1:2, lwd = 2, col = c("#2c7fb8", "#d95f02"))
}
```

### Library Manifold

PCA of the library coordinates shows which HRFs were present (true) versus
selected (matched).

```{r pca-library, fig.width=8, fig.height=5, eval=eval_chunks, fig.alt="PCA scatter of library coordinates showing true and matched HRF positions."}
pca <- prcomp(t(sbhm$A), center = TRUE, scale. = TRUE)
pc <- pca$x[, 1:2, drop = FALSE]
plot(pc, pch = 16, col = "gray70", xlab = "PC1", ylab = "PC2",
     main = "Library Coordinate Space (PCA)")
points(pc[unique(true_hrf_idx), , drop = FALSE], pch = 1, col = "#2c7fb8", lwd = 2)
points(pc[unique(res_sbhm$matched_idx), , drop = FALSE], pch = 4, col = "#d95f02", lwd = 2)
legend("topright", bty = "n",
       legend = c("Library", "True", "Matched"),
       pch = c(16, 1, 4), col = c("gray60", "#2c7fb8", "#d95f02"))
```

## Soft Assignment for Ambiguous Voxels

Hard assignment picks the single best HRF per voxel. When the margin is small,
blending the top-K candidates can reduce variance. The built-in approach
requires just two extra arguments.

```{r soft-blend, eval=eval_chunks}
res_soft <- lss_sbhm(
  Y = Y, sbhm = sbhm, design_spec = design_spec,
  match = list(topK = 3, soft_blend = TRUE,
               blend_margin = median(res_sbhm$margin)),
  return = "amplitude"
)
```

```{r soft-compare, eval=eval_chunks}
cor_sv <- cor(as.vector(res_soft$amplitude), as.vector(res_sbhm$amplitude))
cat("Soft vs hard amplitude correlation:", round(cor_sv, 3), "\n")
```

Blending only applies to voxels whose margin falls below `blend_margin`;
confident voxels keep their hard assignment. For manual control over the
blending weights, use `sbhm_prepass()`, `sbhm_match(topK = K)`, and
`sbhm_project()` directly. See `?sbhm_match` for details.

## Shape Estimation Strategies

By default, SBHM estimates each voxel's HRF shape from an aggregate prepass
fit (`alpha_source = "prepass"`). Two alternative strategies are available when
the prepass is unreliable:

- `"trial_projection"` estimates shape from per-trial projection coefficients in
  the shared basis, using the leading singular vector across trials.
- `"oasis_rank1"` extracts shape from the rank-1 approximation of the
  trial-wise OASIS betas. Set `rank1_min` to gate voxels with low rank-1
  variance fraction back to prepass.

```{r alpha-source, eval=FALSE}
res_proj <- lss_sbhm(
  Y, sbhm, design_spec,
  match = list(alpha_source = "trial_projection")
)
```

### Low-Confidence Gating

For voxels where matching is uncertain (low margin or weak signal), you can gate
them to a fallback shape using `min_margin` and `min_beta_norm`. Gated voxels
use the reference coordinate instead of the matched one.

```{r gating-example, eval=FALSE}
res_gated <- lss_sbhm(
  Y, sbhm, design_spec,
  match = list(min_margin = 0.05, min_beta_norm = 1e-3)
)
cat("Voxels gated to fallback:", res_gated$diag$fallback_low_conf_n, "\n")
```

## Amplitude Policy

The final amplitude stage estimates per-trial scalars from the matched HRF. The
`amplitude$method` argument selects the estimator:
`"lss1"` (default, per-trial 2x2 LSS---robust under overlap),
`"global_ls"` (fast ridge-regularized GLM across all trials), or
`"oasis_voxel"` (full OASIS per voxel, heaviest but returns SEs).

The default `"lss1"` provides the best balance of robustness and accuracy for
typical rapid event-related designs.  For slow designs with well-separated
trials, `"global_ls"` can be faster.  Use `cond_gate` to fall back
automatically when the design is too collinear for a given method.

```{r amplitude-policy, eval=FALSE}
out <- lss_sbhm(
  Y, sbhm, design_spec,
  amplitude = list(method = "lss1",
                   ridge_frac = list(x = 0.02, b = 0.02),
                   cond_gate = list(metric = "rho", thr = 0.999, fallback = "global_ls"))
)
```

See `?lss_sbhm` for guidance on choosing by ISI regime (slow, moderate, fast ER).

## Returning Coefficients for Custom Analysis

When you need the full r-dimensional trial-wise coefficients instead of scalar
amplitudes, request `return = "coefficients"`.

```{r coeffs-demo, eval=eval_chunks}
res_c <- lss_sbhm(Y = Y[, 1:5], sbhm = sbhm,
                   design_spec = design_spec, return = "coefficients")
cat("Coefficients:", dim(res_c$coeffs_r), " [r x trials x voxels]\n")
```

## Parameter Quick Reference

| Parameter | Default | Recommended | Notes |
|-----------|---------|-------------|-------|
| `r` (rank) | -- | 6--12 | Aim for 90--95% variance explained |
| `topK` | 3 | 1--5 | Use 3--5 with `soft_blend = TRUE` for ambiguous cases |
| `soft_blend` | TRUE | TRUE | Blend top-K candidates for uncertain voxels |
| `blend_margin` | 0.08 | 0.05--0.15 | Only blend voxels with margin below this |
| `alpha_source` | `"prepass"` | `"prepass"` | Also `"trial_projection"` or `"oasis_rank1"` |
| `prepass$ridge` | NULL | `list(mode = "fractional", lambda = 0.01)` | Stabilizes noisy/collinear designs |
| `match$shrink$tau` | 0 | 0--0.2 | Increase for low SNR |
| `match$whiten` | FALSE | FALSE | Set TRUE with `whiten_power = 0.5` for partial whitening |
| `match$min_margin` | NULL | 0.05--0.1 | Gate low-confidence voxels to fallback shape |
| `prewhiten` | NULL | `list(method = "ar", p = 1L)` | Use for TR < 2s |
| `amplitude$method` | `"lss1"` | varies by ISI | `"global_ls"` for slow ER, `"lss1"` otherwise |

For full details on every parameter, see `?sbhm_build`, `?sbhm_match`, and
`?lss_sbhm`.

## Performance Considerations

SBHM's cost scales as O(T*r*N*V) for the LSS step, where N is the number of
trials. Compared to unconstrained voxel-wise HRF estimation this is typically
10--50x faster, since r << K*N.

For very large datasets (V > 50,000), PCA factorization reduces memory by
fitting the prepass on q "meta-voxels" instead of V voxels.

```{r factorized, eval=FALSE}
pca_Y <- prcomp(Y, center = TRUE, rank. = 100)
res <- lss_sbhm(
  Y = pca_Y$x, sbhm = sbhm, design_spec = design_spec,
  prepass = list(data_fac = list(scores = pca_Y$x, loadings = pca_Y$rotation))
)
```

For ROI-based analyses, simply subset columns of Y before calling `lss_sbhm()`.
A benchmark script is available via
`system.file("benchmarks", "bench_sbhm.R", package = "fmrilss")`.

## Prewhitening

If your TR is short (< 2s) or residuals show autocorrelation, add prewhitening
to the final OASIS step.

```{r prewhiten-example, eval=FALSE}
res_pw <- lss_sbhm(
  Y = Y, sbhm = sbhm, design_spec = design_spec,
  prewhiten = list(method = "ar", p = 1L, pooling = "global", exact_first = "ar1")
)
```

Note: the factorized prepass intentionally skips prewhitening for efficiency.
The final OASIS step still applies it when `prewhiten` is provided.

## References

- Mumford et al. (2012). "Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses." *NeuroImage*.
- Lindquist et al. (2009). "Modeling the hemodynamic response function in fMRI." *NeuroImage*.

## Summary

SBHM provides efficient, interpretable voxel-specific HRF estimation by
learning a shared basis from a plausible library and matching each voxel in a
low-dimensional coefficient space. The `lss_sbhm()` function runs the full
pipeline in a single call, returning trial-wise amplitudes with per-voxel HRF
assignments and confidence scores.

**Next steps:**
see `?sbhm_build` for library construction,
`?lss_sbhm` for the end-to-end pipeline,
and the "Voxel-wise HRF" vignette for unconstrained alternatives.