params <-
list(family = "red", preset = "homage")

## ----setup, include = FALSE---------------------------------------------------
if (requireNamespace("ggplot2", quietly = TRUE)) ggplot2::theme_set(ggplot2::theme_minimal())
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")
}

## ----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")

## ----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_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")

## ----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_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_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
# Use global=TRUE to get continuous time across blocks
time_points <- samples(sframe, global = TRUE)
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 (block boundary at 240 seconds)
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)