---
title: "Advanced HRF Modeling and Design"
author: "Bradley R. Buchsbaum"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced HRF Modeling and Design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, 
  fig.height = 5,
  message = FALSE,
  warning = FALSE
)
library(fmrihrf)
library(dplyr)
library(ggplot2)
library(tidyr)
library(Matrix)
if (requireNamespace("viridis", quietly = TRUE)) {
  library(viridis)
} else {
  scale_color_viridis_d <- function(...) ggplot2::scale_color_brewer(palette = "Set1")
}
```

## Introduction

This vignette explores advanced features of `fmrihrf` for systematic HRF modeling, regularization, and experimental design. We'll cover five key functions that extend the basic HRF framework:

- **`hrf_library()`**: Creating systematic collections of HRF variants
- **`reconstruction_matrix()`**: Converting basis coefficients back to HRF shapes
- **`regressor_set()`**: Managing multi-condition experimental designs
- **`regressor_design()`**: Building design matrices for complex experimental blocks

These tools are essential for advanced fMRI modeling where you need flexibility in HRF specification, robust estimation with limited data, or complex experimental designs.

## HRF Libraries: Systematic Parameter Exploration

The `hrf_library()` function creates collections of HRF variants by systematically varying parameters. This is useful for exploring how different HRF assumptions affect your model or for building data-driven HRF basis sets.

### Example 1: Library of Gamma HRFs

Let's create a library of gamma HRFs with different shape and rate parameters:

```{r gamma_library}
# Define parameter grid for gamma HRFs
gamma_params <- expand.grid(
  shape = c(4, 6, 8),
  rate = c(0.8, 1.0, 1.2)
)
print(gamma_params)

# Create a generator function for gamma HRFs
make_gamma_hrf <- function(shape, rate) {
  gen_hrf(hrf_gamma, shape = shape, rate = rate, name = paste0("Gamma_", shape, "_", rate))
}

# Create HRF library
gamma_lib <- hrf_library(make_gamma_hrf, gamma_params)
print(gamma_lib)
nbasis(gamma_lib) # 9 HRFs total (3 x 3 grid)

# Evaluate and visualize
time_points <- seq(0, 20, by = 0.1)
gamma_responses <- gamma_lib(time_points)

# Convert to long format for plotting
gamma_df <- as.data.frame(gamma_responses)
names(gamma_df) <- paste0("Shape", gamma_params$shape, "_Rate", gamma_params$rate)
gamma_df$Time <- time_points

gamma_long <- pivot_longer(gamma_df, -Time, names_to = "Parameters", values_to = "Response")

# Create a more informative plot
ggplot(gamma_long, aes(x = Time, y = Response, color = Parameters)) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d() +
  labs(title = "Library of Gamma HRFs",
       subtitle = "Systematic variation of shape and rate parameters",
       x = "Time (seconds)",
       y = "HRF Response") +
  theme_minimal() +
  theme(legend.position = "right")
```

### Example 2: Library of Lagged SPM HRFs

Here's how to create a library of the SPM canonical HRF with different temporal lags:

```{r spm_lag_library}
# Parameter grid for temporal lags
lag_params <- data.frame(lag = seq(-2, 4, by = 1))
print(lag_params)

# Create library using a helper function that applies lag_hrf
create_lagged_spm <- function(lag) {
  lag_hrf(HRF_SPMG1, lag = lag)
}

spm_lag_lib <- hrf_library(create_lagged_spm, lag_params)
print(spm_lag_lib)

# Evaluate and plot
spm_lag_responses <- spm_lag_lib(time_points)
spm_lag_df <- as.data.frame(spm_lag_responses)
names(spm_lag_df) <- paste0("Lag_", lag_params$lag, "s")
spm_lag_df$Time <- time_points

spm_lag_long <- pivot_longer(spm_lag_df, -Time, names_to = "Lag", values_to = "Response")

ggplot(spm_lag_long, aes(x = Time, y = Response, color = Lag)) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d() +
  labs(title = "Library of Lagged SPM Canonical HRFs",
       subtitle = "Temporal lags from -2 to +4 seconds",
       x = "Time (seconds)",
       y = "HRF Response") +
  theme_minimal()
```


## Reconstruction Matrices: From Coefficients to HRF Shapes

The reconstruction process converts a set of basis coefficients into a continuous HRF shape. Understanding this transformation is key to interpreting estimated HRFs from fMRI analyses.

### How Reconstruction Works

```{r reconstruction_demo}
# Use a small basis for clear visualization
basis_set <- gen_hrf(hrf_bspline, N = 5, degree = 3, span = 30)
eval_times <- seq(0, 30, by = 0.1)

# The reconstruction matrix: each column is a basis function evaluated at time points
recon_matrix <- basis_set(eval_times)
print(paste("Reconstruction matrix dimensions:", nrow(recon_matrix), "time points x", 
            ncol(recon_matrix), "basis functions"))

# Let's visualize the basis functions themselves first
basis_df <- as.data.frame(recon_matrix)
names(basis_df) <- paste0("B", 1:5)
basis_df$Time <- eval_times

basis_long <- pivot_longer(basis_df, -Time, names_to = "Basis", values_to = "Value")

ggplot(basis_long, aes(x = Time, y = Value, color = Basis)) +
  geom_line(linewidth = 1.2) +
  scale_color_viridis_d(option = "turbo") +
  labs(title = "B-spline Basis Functions",
       subtitle = "Each basis function covers a different time window",
       x = "Time (seconds)",
       y = "Basis Function Value") +
  theme_minimal()

# Now demonstrate reconstruction with different coefficient patterns
coefficient_sets <- list(
  "Early Peak" = c(0.2, 1.0, 0.3, 0.0, 0.0),
  "Canonical" = c(0.0, 0.3, 1.0, 0.4, -0.1),
  "Late Peak" = c(0.0, 0.0, 0.3, 1.0, 0.2),
  "Double Peak" = c(0.0, 0.8, 0.2, 0.9, 0.0)
)

# Reconstruct HRFs for each coefficient set
reconstruction_df <- data.frame()
for (name in names(coefficient_sets)) {
  coefs <- coefficient_sets[[name]]
  hrf_values <- as.vector(recon_matrix %*% coefs)
  
  df <- data.frame(
    Time = eval_times,
    HRF = hrf_values,
    Pattern = name
  )
  reconstruction_df <- rbind(reconstruction_df, df)
}

ggplot(reconstruction_df, aes(x = Time, y = HRF, color = Pattern)) +
  geom_line(linewidth = 1.5) +
  scale_color_manual(values = c("Early Peak" = "#E69F00", 
                               "Canonical" = "#009E73",
                               "Late Peak" = "#0072B2",
                               "Double Peak" = "#D55E00")) +
  labs(title = "Different HRF Shapes from Same Basis Set",
       subtitle = "Varying coefficients produces diverse HRF patterns",
       x = "Time (seconds)",
       y = "HRF Response") +
  theme_minimal() +
  theme(legend.position = "bottom")
```

### Interactive Visualization: Building an HRF Step by Step

```{r reconstruction_interactive}
# Let's build up a canonical HRF step by step
canonical_coefs <- c(0.0, 0.3, 1.0, 0.4, -0.1)

# Create data for cumulative reconstruction
cumulative_df <- data.frame()
for (i in 1:5) {
  # Zero out coefficients after position i
  temp_coefs <- canonical_coefs
  if (i < 5) temp_coefs[(i+1):5] <- 0
  
  # Calculate cumulative HRF
  cumulative_hrf <- as.vector(recon_matrix %*% temp_coefs)
  
  # Store individual contribution
  individual_coefs <- rep(0, 5)
  individual_coefs[i] <- canonical_coefs[i]
  individual_contribution <- as.vector(recon_matrix %*% individual_coefs)
  
  df <- data.frame(
    Time = rep(eval_times, 2),
    Value = c(cumulative_hrf, individual_contribution),
    Type = rep(c("Cumulative", "Individual"), each = length(eval_times)),
    Step = i,
    Basis = paste0("Adding B", i, " (coef=", round(canonical_coefs[i], 2), ")")
  )
  cumulative_df <- rbind(cumulative_df, df)
}

# Create faceted plot showing the build-up
ggplot(cumulative_df, aes(x = Time, y = Value, color = Type)) +
  geom_line(linewidth = 1.2) +
  facet_wrap(~Basis, ncol = 5) +
  scale_color_manual(values = c("Cumulative" = "black", "Individual" = "red")) +
  labs(title = "Building an HRF: Sequential Addition of Weighted Basis Functions",
       subtitle = "Red: individual contribution, Black: cumulative sum",
       x = "Time (seconds)",
       y = "Value") +
  theme_minimal() +
  theme(legend.position = "bottom",
        strip.text = element_text(size = 9))

# Show coefficient importance
coef_importance <- data.frame(
  Basis = paste0("B", 1:5),
  Coefficient = canonical_coefs,
  `Absolute Value` = abs(canonical_coefs)
)

ggplot(coef_importance, aes(x = Basis, y = Coefficient, fill = Coefficient > 0)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  scale_fill_manual(values = c("FALSE" = "#D55E00", "TRUE" = "#009E73"),
                    labels = c("Negative", "Positive")) +
  labs(title = "Coefficient Values for Canonical HRF",
       subtitle = "B3 dominates the shape, B5 provides the undershoot",
       x = "Basis Function",
       y = "Coefficient Value",
       fill = "Sign") +
  theme_minimal()
```

## Regressor Sets: Multi-Condition Experimental Designs

The `regressor_set()` function simplifies creating regressors for multi-condition experiments where each condition shares the same HRF but has different event timings.

```{r regressor_set_demo}
# Simulate a 3-condition experiment
set.seed(123)
n_events_per_condition <- 8
total_duration <- 240  # 4 minutes

# Generate random onsets for each condition
condition_A_onsets <- sort(runif(n_events_per_condition, 0, total_duration))
condition_B_onsets <- sort(runif(n_events_per_condition, 0, total_duration))
condition_C_onsets <- sort(runif(n_events_per_condition, 0, total_duration))

# Combine all onsets and create factor
all_onsets <- c(condition_A_onsets, condition_B_onsets, condition_C_onsets)
conditions <- factor(rep(c("TaskA", "TaskB", "TaskC"), each = n_events_per_condition))

# Create regressor set
reg_set <- regressor_set(onsets = all_onsets, fac = conditions, hrf = HRF_SPMG1)
print(reg_set)

# Evaluate at scan times (TR = 2s)
TR <- 2
scan_times <- seq(0, total_duration, by = TR)
design_matrix <- evaluate(reg_set, scan_times)

print(dim(design_matrix)) # Time points x 3 conditions

# Visualize the design matrix
design_df <- as.data.frame(design_matrix)
names(design_df) <- c("TaskA", "TaskB", "TaskC")
design_df$Time <- scan_times

design_long <- pivot_longer(design_df, -Time, names_to = "Condition", values_to = "Response")

ggplot(design_long, aes(x = Time, y = Response, color = Condition)) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d() +
  labs(title = "Multi-Condition fMRI Design Matrix",
       subtitle = "Three experimental conditions with shared HRF",
       x = "Time (seconds)",
       y = "Predicted BOLD Response",
       color = "Condition") +
  theme_minimal()

# Add event markers
onset_df <- data.frame(
  Time = all_onsets,
  Condition = conditions,
  Marker = 1
)

ggplot(design_long, aes(x = Time, y = Response, color = Condition)) +
  geom_line(linewidth = 1) +
  geom_point(data = onset_df, aes(x = Time, y = -0.1, color = Condition), 
             size = 2, alpha = 0.7) +
  scale_color_viridis_d() +
  labs(title = "Design Matrix with Event Onsets",
       subtitle = "Points show stimulus onset times",
       x = "Time (seconds)",
       y = "Predicted BOLD Response",
       color = "Condition") +
  theme_minimal()
```

## Regressor Design: Complex Block Designs

For more complex experimental designs with multiple blocks or runs, `regressor_design()` provides a higher-level interface that handles block-relative timing and creates design matrices directly.

```{r regressor_design_demo}
# Create a sampling frame for 2 blocks of 120 seconds each
sframe <- sampling_frame(
  blocklens = c(120, 120),  # Two 4-minute blocks (120 scans each at TR = 2s)
  TR = 2                    # 2-second TR
)
print(sframe)

# Generate block-relative event onsets
# Block 1: Faces at 10, 50, 90; Houses at 30, 70 seconds
# Block 2: Faces at 15, 55, 95; Houses at 35, 75 seconds  
block_onsets <- c(10, 30, 50, 70, 90, 15, 35, 55, 75, 95)
block_ids <- c(rep(1, 5), rep(2, 5))
event_conditions <- factor(c("Faces", "Houses", "Faces", "Houses", "Faces", 
                            "Faces", "Houses", "Faces", "Houses", "Faces"))

# Create design matrix using regressor_design
design_mat <- regressor_design(
  onsets = block_onsets,
  fac = event_conditions,
  block = block_ids,
  sframe = sframe,
  hrf = HRF_SPMG1
)

print(dim(design_mat)) # Total time points across both blocks x 2 conditions

# Convert to data frame for plotting
time_points <- samples(sframe)
design_plot_df <- as.data.frame(design_mat)
names(design_plot_df) <- c("Faces", "Houses")
design_plot_df$Time <- time_points
design_plot_df$Block <- rep(1:2, each = 120) # 120 scans per block

design_plot_long <- pivot_longer(design_plot_df, c("Faces", "Houses"), 
                                names_to = "Condition", values_to = "Response")

# Plot with block separation
ggplot(design_plot_long, aes(x = Time, y = Response, color = Condition)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = 240, linetype = "dashed", alpha = 0.7) +
  scale_color_viridis_d() +
  labs(title = "Multi-Block Experimental Design",
       subtitle = "Two blocks with different event schedules (dashed line = block boundary)",
       x = "Time (seconds)",
       y = "Predicted BOLD Response",
       color = "Condition") +
  theme_minimal()

# Show global vs block-relative timing
timing_df <- data.frame(
  Block = block_ids,
  Block_Relative_Onset = block_onsets,
  Global_Onset = global_onsets(sframe, block_onsets, block_ids),
  Condition = event_conditions
)

print(timing_df)
```