---
title: "Time Complexity, Inferential uncertainty and Spacekime Analytics"
subtitle: "Appendix 03: Kime-Phase Tomography (KPT V.5: KPT-GEM & KPT-FFT): fMRI and Double-Slit Validation"
author: "SOCR Team"
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output:
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
    pandoc_args: ["+RTS", "-K1024m", "-RTS"]
  word_document:
    toc: true
    toc_depth: '3'
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amssymb}
  - \usepackage{physics}
  - \usepackage{mathrsfs}
  - \usepackage{bm}
---


This [TCIU Appendix](https://tciu.predictive.space/) offers a self-contained implementation of the kime phase tompgraphy algorithms, *KPT-GEM* (expectation maximization) and *KPT-FFT (FFT alternating spectral)* for estimating
the phase distribution and reconstructing the kime-surface representations from observed repeated measurement data from longitudinal processes. The KPT algorihtm is validated using two simulations; *double-slit physics experiment* and *fMRI ON vs. OFF neuroscience experiment*. We report quyantitative performance metrics (summary tables) and provide qualitative results as surface/heatmap overlays contrasting for the *true* vs. *estimated* phase distributions and kime-surfaces.

The key steps of the KPT algorithm include E-step, Wiener-Sobolev filtering, simplex projection, ans reality constraints, subject to normalization, operator conventions.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(Matrix)
  library(plotly)
  library(knitr)
  library(purrr)
})
set.seed(2025)
```

# Introduction

This notebook presents a unified implementation of *Kime-Phase Tomography (KPT)*, a framework for analyzing repeated measurement longitudinal data by decomposing observations into

 1. A deterministic *kime-surface* $\mathcal{S}(t,\theta)$ over complex time $\kappa = te^{i\theta}$, and
 2. A stochastic *phase distribution* $\varphi_t(\theta)$ capturing intrinsic variability.

The forward model for observations is
$$y_{j,k} = \mathcal{S}(t_k, \Theta_j(t_k)) + \varepsilon_{j,k},$$
where $\Theta_j(t_k) \sim \varphi_{t_k}(\cdot)$ is a latent phase and $\varepsilon_{j,k} \sim \mathcal{N}(0,\sigma^2)$ is measurement noise.

# Mathematical Framework

## KPT Inverse Problem

Given repeated measurements $\{y_{j,k}\}_{j=1}^N$ at times $\{t_k\}_{k=1}^K$, KPT solves the inverse problem:

*Recover:* 

 - Phase distributions $\{\varphi_{t_k}\}$
 - Kime-surface $\mathcal{S}(t,\theta)$

*Key insight:* This deconvolution becomes tractable in the Fourier domain where:
$$M(z,t) = F(z,t) \cdot \Phi(z,t), \quad |z|=1$$

where $F(z,t) = \sum_n f_n(t)z^n$ (surface transform), $\Phi(z,t) = \sum_n \hat{\varphi}_t(n)z^n$ (phase characteristic function), and $M(z,t) = \sum_n m_n(t)z^n$ (observable mixed moments).

## KPT Algorithms

Both algorithms implement the same core spectral deconvolution with different update strategies:

### Algorithm 1: KPT-GEM (Generalized EM)

*Inputs:* Data $\{y_{j,k}\}$, regularization $(\lambda_\Phi, p_\Phi)$, $(\lambda_F, p_F)$

*Initialize:* $\varphi^{(0)}_{t_k} \equiv \frac{1}{2\pi}$, $f^{(0)}_n(t_k)$ from data

*Repeat until convergence:*

1. *E-step:* Compute posterior expectations
   $$\xi_{j,k}^{(\ell)}(n) = \mathbb{E}[e^{-in\Theta_j} | y_{j,k}] = \int_0^{2\pi} e^{-in\theta} w_{j,k}^{(\ell)}(\theta) \frac{d\theta}{2\pi}$$
   where $w_{j,k}^{(\ell)}(\theta) \propto \varphi^{(\ell-1)}_{t_k}(\theta) \exp\left(-\frac{(y_{j,k} - \mathcal{S}^{(\ell-1)}(t_k,\theta))^2}{2\sigma^2}\right)$
   
2. *Mixed moments:* $m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)$

3. *Φ-step (phase update):* Apply Wiener filter
   $$\hat{\Phi}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{F}^{(\ell)}(\omega,t_k)}}
   {|\hat{F}^{(\ell)}(\omega,t_k)|^2 + \lambda_\Phi \Lambda_{p_\Phi}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$
   
4. *Optional F-step (surface update):*
   $$\hat{F}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{\Phi}^{(\ell+1)}(\omega,t_k)}}
   {|\hat{\Phi}^{(\ell+1)}(\omega,t_k)|^2 + \lambda_F \Lambda_{p_F}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$

5. *Project:* Ensure $\varphi^{(\ell+1)}_{t_k} \geq 0$ and $\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1$

### Algorithm 2: KPT-FFT (Fully Alternating)

Same as KPT-GEM but *always* alternates between Φ-step and F-step updates at each iteration.
The Sobolev penalty $\Lambda_p(\omega) = (2\sin(\omega/2))^{2p}$ controls smoothness.

<!-- ### Algorithm 1: KPT-GEM (Generalized EM) -->
*Inputs:* Data $\{y_{j,k}\}$, regularization $(\lambda_\Phi, p_\Phi)$, $(\lambda_F, p_F)$

*Initialize:* $\varphi^{(0)}_{t_k} \equiv \frac{1}{2\pi}$, $f^{(0)}_n(t_k)$ from data

*Repeat until convergence:*

1. *E-step:* Compute posterior expectations
   $$\xi_{j,k}^{(\ell)}(n) = \mathbb{E}[e^{-in\Theta_j} | y_{j,k}] = \int_0^{2\pi} e^{-in\theta} w_{j,k}^{(\ell)}(\theta) \frac{d\theta}{2\pi}$$
   where $w_{j,k}^{(\ell)}(\theta) \propto \varphi^{(\ell-1)}_{t_k}(\theta) \exp\left(-\frac{(y_{j,k} - \mathcal{S}^{(\ell-1)}(t_k,\theta))^2}{2\sigma^2}\right)$
   
2. *Mixed moments:* $m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)$

3. $\Phi$*-step (phase update):* Apply Wiener filter
   $$\hat{\Phi}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{F}^{(\ell)}(\omega,t_k)}}
   {|\hat{F}^{(\ell)}(\omega,t_k)|^2 + \lambda_\Phi \Lambda_{p_\Phi}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$
   
4. *F-step (surface update):*
   $$\hat{F}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{\Phi}^{(\ell+1)}(\omega,t_k)}}
   {|\hat{\Phi}^{(\ell+1)}(\omega,t_k)|^2 + \lambda_F \Lambda_{p_F}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$

5. *Project:* Ensure $\varphi^{(\ell+1)}_{t_k} \geq 0$ and $\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1$


# Implementation

## Core Utilities

```{r utilities}
# ---- Grid generation ----
theta_grid <- function(L) seq(0, 2*pi, length.out = L + 1L)[1L:L]
make_nvec <- function(J) seq.int(-J, J, by = 1L)

# ---- Probability simplex projection ----
row_project_simplex <- function(P, eps = 1e-12) {
  P <- pmax(P, 0)
  m <- rowMeans(P)
  m[m < eps] <- 1
  sweep(P, 1L, m, `/`)
}

# ---- Stable softmax for posterior weights ----
softmax_rows <- function(Lmat) {
  m <- apply(Lmat, 1L, max)
  W <- exp(Lmat - m)
  W / pmax(rowSums(W), .Machine$double.eps)
}

# ---- FFT coefficient mapping (n=-J..J to DFT order) ----
coeffs_to_dft_order <- function(x_n, nvec) {
  N <- length(nvec)
  x_m <- complex(length = N)
  for (i in seq_along(nvec)) {
    n <- nvec[i]
    m <- (n + N) %% N
    x_m[m + 1L] <- x_n[i]
  }
  x_m
}

dft_order_to_coeffs <- function(x_m, nvec) {
  N <- length(nvec)
  x_n <- complex(length = N)
  for (i in seq_along(nvec)) {
    n <- nvec[i]
    m <- (n + N) %% N
    x_n[i] <- x_m[m + 1L]
  }
  x_n
}

# ---- FFT-based evaluation on unit circle ----
coeffs_to_Xomega_fft <- function(x_n, nvec) {
  x_m <- coeffs_to_dft_order(x_n, nvec)
  fft(x_m, inverse = FALSE)
}

Xomega_to_coeffs_fft <- function(X_w, nvec) {
  N <- length(nvec)
  x_m <- fft(X_w, inverse = TRUE) / N
  dft_order_to_coeffs(x_m, nvec)
}

# ---- Surface synthesis from Fourier coefficients ----
surface_from_coeffs <- function(f_n_row, nvec, theta) {
  E <- exp(1i * outer(nvec, theta))
  as.vector(Re(drop(f_n_row %*% E)))
}

# ---- Reality constraint: f_{-n} = conj(f_n) ----
enforce_reality_coeffs <- function(f_n, nvec) {
  N <- length(nvec)
  J <- (N - 1L) / 2L
  
  # n=0 must be real
  f_n[nvec == 0] <- Re(f_n[nvec == 0])
  
  # Conjugate symmetry
  for (j in 1:J) {
    pos_idx <- which(nvec == j)
    neg_idx <- which(nvec == -j)
    if (length(pos_idx) > 0 && length(neg_idx) > 0) {
      avg <- 0.5 * (f_n[pos_idx] + Conj(f_n[neg_idx]))
      f_n[pos_idx] <- avg
      f_n[neg_idx] <- Conj(avg)
    }
  }
  f_n
}

# ---- Phase distribution metrics ----
phase_metrics_df <- function(phi_true, phi_est, theta, time, base2 = TRUE) {
  # Normalize
  P <- pmax(phi_true, 1e-12)
  Q <- pmax(phi_est, 1e-12)
  P <- sweep(P, 1, rowMeans(P), `/`)
  Q <- sweep(Q, 1, rowMeans(Q), `/`)
  
  K <- nrow(P)
  logb <- if (base2) function(x) log2(pmax(x, 1e-12)) else log
  
  metrics <- data.frame(
    time = time,
    KL_true_est = numeric(K),
    KL_est_true = numeric(K),
    JSD_bits = numeric(K),
    Hellinger = numeric(K),
    TV = numeric(K)
  )
  
  for (k in 1:K) {
    # KL divergences
    metrics$KL_true_est[k] <- mean(P[k,] * logb(P[k,]/Q[k,]))
    metrics$KL_est_true[k] <- mean(Q[k,] * logb(Q[k,]/P[k,]))
    
    # Jensen-Shannon divergence
    M <- 0.5 * (P[k,] + Q[k,])
    metrics$JSD_bits[k] <- 0.5 * mean(P[k,] * logb(P[k,]/M)) + 
                           0.5 * mean(Q[k,] * logb(Q[k,]/M))
    
    # Hellinger distance
    metrics$Hellinger[k] <- sqrt(0.5 * mean((sqrt(P[k,]) - sqrt(Q[k,]))^2))
    
    # Total variation
    metrics$TV[k] <- 0.5 * mean(abs(P[k,] - Q[k,]))
  }
  
  metrics
}

# ---- Surface L2 error (cone-weighted) ----
surface_L2_cone <- function(S_true, S_est, time, theta) {
  K <- nrow(S_true)
  L <- ncol(S_true)
  dtheta <- 2*pi / L
  
  # Time weights
  w_t <- time
  if (length(w_t) == 0 || K == 1) w_t <- rep(1, K)
  
  # Compute L2 norm
  diff2 <- (S_true - S_est)^2
  per_t_L2_sq <- rowSums(diff2) * dtheta * w_t
  total_L2 <- sqrt(sum(per_t_L2_sq))
  
  # Relative error
  true_norm_sq <- rowSums(S_true^2) * dtheta * w_t
  rel_L2 <- total_L2 / max(sqrt(sum(true_norm_sq)), 1e-12)
  
  list(total_L2 = total_L2, rel_L2 = rel_L2)
}

# ---- Visualization helpers ----
plot_phase_overlay <- function(theta, phi_true_row, phi_est_row1, 
                              phi_est_row2 = NULL, t_lab = "", title = "") {
  if (is.null(phi_est_row2)) {
    df <- data.frame(
      theta = rep(theta, 2),
      density = c(phi_true_row, phi_est_row1),
      Method = rep(c("True", "Estimated"), each = length(theta))
    )
  } else {
    df <- data.frame(
      theta = rep(theta, 3),
      density = c(phi_true_row, phi_est_row1, phi_est_row2),
      Method = rep(c("True", "KPT-GEM", "KPT-FFT"), each = length(theta))
    )
  }
  
  ggplot(df, aes(theta, density, color = Method)) +
    geom_line(linewidth = 1.2) +
    labs(title = title, subtitle = t_lab,
         x = expression(theta), y = expression(varphi(theta))) +
    theme_minimal() +
    theme(legend.position = "bottom")
}
```

## Unified KPT Implementation

```{r kpt-unified}
# ---- Initialize surface from data ----
initialize_surface <- function(Y, time, theta, method = "empirical") {
  K <- ncol(Y)
  L <- length(theta)
  
  if (method == "empirical") {
    # Use mean and variance structure
    mu_k <- colMeans(Y)
    A_k <- apply(Y, 2, sd)
    
    # Smooth estimates
    if (K > 4) {
      mu_s <- smooth.spline(time, mu_k, spar = 0.6)$y
      A_s <- smooth.spline(time, A_k, spar = 0.6)$y
    } else {
      mu_s <- mu_k
      A_s <- A_k
    }
    
    # Build surface: baseline + amplitude * harmonics
    S <- matrix(0, K, L)
    for (k in 1:K) {
      S[k,] <- mu_s[k] + A_s[k] * cos(theta) + 0.25 * A_s[k] * cos(2*theta)
    }
  } else {
    # Simple constant initialization
    S <- matrix(mean(Y), K, L)
  }
  
  S
}

# ---- Main unified KPT function ----
kpt_unified <- function(Y, time = NULL,
                       J = 12, L_theta = 256,
                       lambda_phi = 1e-2, p_phi = 1,
                       lambda_F = 1e-2, p_F = 1,
                       algorithm = c("GEM", "FFT"),
                       update_surface = TRUE,
                       max_iter = 30, tol = 1e-3,
                       sigma = NULL, verbose = TRUE) {
  
  algorithm <- match.arg(algorithm)
  I <- nrow(Y)  # Number of subjects/runs
  K <- ncol(Y)  # Number of time points
  
  # Default time grid
  if (is.null(time)) time <- seq_len(K)
  
  # Setup grids
  theta <- theta_grid(L_theta)
  nvec <- make_nvec(J)
  Nn <- length(nvec)
  omega <- 2*pi * (0:(Nn-1)) / Nn
  
  # Estimate noise if not provided
  if (is.null(sigma)) {
    sigma <- 0.5 * mad(as.vector(Y), center = 0) + 1e-6
  }
  
  # Initialize phase distribution (uniform)
  phi_grid <- matrix(1.0, nrow = K, ncol = L_theta)
  
  # Initialize surface
  Sgrid <- initialize_surface(Y, time, theta)
  
  # Compute initial Fourier coefficients
  f_n <- matrix(0+0i, nrow = K, ncol = Nn)
  for (k in 1:K) {
    Eneg <- exp(-1i * outer(theta, nvec))
    f_n[k,] <- as.vector((t(Sgrid[k,]) %*% Eneg) / L_theta)
    f_n[k,] <- enforce_reality_coeffs(f_n[k,], nvec)
  }
  
  # Sobolev penalties
  Lambda_phi <- (2 * sin(omega/2))^(2*p_phi)
  Lambda_F <- (2 * sin(omega/2))^(2*p_F)
  
  # Storage for convergence
  obj <- numeric(max_iter)
  
  # Main iteration loop
  for (it in 1:max_iter) {
    
    # ==== E-step: Compute posterior expectations ====
    m_hat <- matrix(0+0i, nrow = K, ncol = Nn)
    
    for (k in 1:K) {
      # Prior at time k
      prior_k <- phi_grid[k,]
      
      # Likelihood: p(y|theta)
      Srow <- matrix(Sgrid[k,], nrow = I, ncol = L_theta, byrow = TRUE)
      LogL <- -(Y[,k] - Srow)^2 / (2 * sigma^2)
      
      # Posterior weights (properly normalized)
      LogW <- sweep(LogL, 2L, log(pmax(prior_k, 1e-12)), `+`)
      W <- softmax_rows(LogW)  # [I x L_theta]
      
      # Compute expectations E[e^{-in*Theta} | y]
      Eneg <- exp(-1i * outer(theta, nvec))
      xi <- W %*% Eneg  # [I x Nn]
      
      # Mixed moments (corrected)
      m_hat[k,] <- colSums(Y[,k] * xi) / I
    }
    
    # ==== Φ-step: Update phase distribution ====
    phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn)
    
    for (k in 1:K) {
      # FFT to evaluate on unit circle
      F_w <- coeffs_to_Xomega_fft(f_n[k,], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k,], nvec)
      
      # Wiener filter
      denom <- Mod(F_w)^2 + lambda_phi * Lambda_phi + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      
      # IFFT back to coefficients
      phi_hat_n[k,] <- Xomega_to_coeffs_fft(Phi_w, nvec)
    }
    
    # Synthesize and project to simplex
    phi_new <- matrix(0, nrow = K, ncol = L_theta)
    Epos <- exp(1i * outer(nvec, theta))
    for (k in 1:K) {
      phi_row <- Re(drop(t(Epos) %*% phi_hat_n[k,]))
      phi_new[k,] <- phi_row
    }
    phi_new <- row_project_simplex(phi_new)
    
    # ==== F-step: Update surface (if requested) ====
    if (update_surface && (algorithm == "FFT" || it > 5)) {
      f_n_new <- matrix(0+0i, nrow = K, ncol = Nn)
      
      for (k in 1:K) {
        # Use updated phase
        Phi_w <- coeffs_to_Xomega_fft(phi_hat_n[k,], nvec)
        M_w <- coeffs_to_Xomega_fft(m_hat[k,], nvec)
        
        # Wiener filter for surface
        denom <- Mod(Phi_w)^2 + lambda_F * Lambda_F + 1e-10
        F_w <- Conj(Phi_w) * M_w / denom
        
        # Reality constraint
        f_n_new[k,] <- enforce_reality_coeffs(Xomega_to_coeffs_fft(F_w, nvec), nvec)
      }
      
      # Synthesize new surface
      S_new <- matrix(0, nrow = K, ncol = L_theta)
      for (k in 1:K) {
        S_new[k,] <- surface_from_coeffs(f_n_new[k,], nvec, theta)
      }
      
      # Update
      f_n <- f_n_new
      Sgrid <- S_new
    }
    
    # ==== Check convergence ====
    diff_phi <- sqrt(mean((phi_new - phi_grid)^2))
    diff_f <- if (update_surface && (algorithm == "FFT" || it > 5)) {
      sqrt(mean(Mod(f_n_new - f_n)^2))
    } else {
      0
    }
    
    if (verbose) {
      if (algorithm == "GEM") {
        message(sprintf("[%s] iter %02d: Δφ=%.3e", algorithm, it, diff_phi))
      } else {
        message(sprintf("[%s] iter %02d: Δφ=%.3e, Δf=%.3e", algorithm, it, diff_phi, diff_f))
      }
    }
    
    # Update state
    phi_grid <- phi_new
    obj[it] <- mean(Mod(m_hat)^2)
    
    # Convergence check
    if (algorithm == "GEM" && diff_phi < tol) break
    if (algorithm == "FFT" && max(diff_phi, diff_f) < tol) break
  }
  
  # Return results
  list(
    varphi_theta = phi_grid,
    Sgrid = Sgrid,
    f_n = f_n,
    theta = theta,
    time = time,
    convergence = list(
      iterations = it,
      objective = obj[1:it]
    ),
    parameters = list(
      algorithm = algorithm,
      lambda_phi = lambda_phi, p_phi = p_phi,
      lambda_F = lambda_F, p_F = p_F,
      sigma = sigma, J = J
    )
  )
}
```


# Double-Slit Experiment (with Phase Alignment)

Since the phase distributions are only identified up to a *global rotation*, we need
to align them prior to comparison, visualizaiton, contextual interpretation, or 
spacekime analytic prediction.

## Setup and Utilities

```{r setup-utilities}
# Core utilities for phase alignment
anchor_phase_by_m1 <- function(varphi_mat, theta) {
  L <- length(theta)
  dtheta <- 2*pi / L
  
  # Compute first trigonometric moment per row
  m1 <- apply(varphi_mat, 1L, function(row) {
    sum(exp(1i * theta) * row) * dtheta / (2*pi)
  })
  
  # Extract rotation angles
  delta <- Arg(m1)
  
  # Rotate each row by -delta using FFT
  rotate_row <- function(row, delta_k) {
    F <- fft(row)
    n <- c(0:(L/2), -(L/2-1):-1)
    row_rot <- Re(fft(F * exp(-1i * n * delta_k), inverse = TRUE) / L)
    pmax(row_rot, 0)  # Ensure non-negative
  }
  
  # Apply rotation
  out <- mapply(function(idx, d) rotate_row(varphi_mat[idx,], d),
                idx = seq_len(nrow(varphi_mat)), 
                d = delta, 
                SIMPLIFY = FALSE)
  out <- do.call(rbind, out)
  
  # Renormalize
  out / rowMeans(out)
}

rotate_surface_by_delta <- function(S_mat, theta, delta_vec) {
  stopifnot(nrow(S_mat) == length(delta_vec))
  L <- length(theta)
  n <- c(0:(L/2), -(L/2-1):-1)
  
  rotate_row <- function(row, delta_k) {
    F <- fft(row)
    Re(fft(F * exp(-1i * n * delta_k), inverse = TRUE) / L)
  }
  
  out <- mapply(function(idx, d) rotate_row(S_mat[idx,], d),
                idx = seq_len(nrow(S_mat)), 
                d = delta_vec, 
                SIMPLIFY = FALSE)
  do.call(rbind, out)
}
```

## Double-Slit Simulation

```{r simulate-ds}
simulate_double_slit_aligned <- function(N = 100, K = 50, L_theta = 256,
                                        kappa_phase = 4, noise_sd = 0.15,
                                        kappa_wave = 12, seed = 2025) {
  set.seed(seed)
  
  # Grids
  theta <- theta_grid(L_theta)
  time <- seq(0, 1, length.out = K)
  
  # Build true kime-surface: quantum interference pattern
  S_true <- matrix(0, K, L_theta)
  for (k in 1:K) {
    t <- time[k]
    # Gaussian envelope for wave packet
    envelope <- 0.8 * exp(-20*(t - 0.5)^2)
    # Base drift
    baseline <- 0.2 * sin(2*pi*t)
    # Interference pattern with two harmonics
    S_true[k,] <- baseline + 
                  envelope * cos(kappa_wave*t + theta) + 
                  0.3 * envelope * cos(2*(kappa_wave*t + theta) + pi/6)
  }
  
  # Time-varying phase distribution (mixture of von Mises)
  # This represents the quantum state distribution
  phi_true <- matrix(0, K, L_theta)
  for (k in 1:K) {
    t <- time[k]
    # Two modes representing different quantum paths
    mu1 <- (2*pi * (0.2 + 0.6*t)) %% (2*pi)
    mu2 <- (2*pi * (0.9 - 0.6*t)) %% (2*pi)
    w <- 0.5 + 0.4*sin(2*pi*t)  # Time-varying mixture weight
    
    # von Mises mixture
    p1 <- exp(kappa_phase * cos(theta - mu1)) / (2*pi * besselI(kappa_phase, 0))
    p2 <- exp(kappa_phase * cos(theta - mu2)) / (2*pi * besselI(kappa_phase, 0))
    p <- w * p1 + (1-w) * p2
    phi_true[k,] <- p / mean(p)  # Normalize to mean=1
  }
  
  # Generate observations: Y = S(t, Theta) + noise
  Y <- matrix(0, N, K)
  for (k in 1:K) {
    # Sample phases from the true distribution
    probs <- pmax(phi_true[k,], 1e-12)
    idx <- sample.int(L_theta, N, replace = TRUE, prob = probs)
    
    # Generate measurements
    Y[,k] <- S_true[k, idx] + rnorm(N, 0, noise_sd)
  }
  
  list(
    Y = Y,
    time = time,
    theta = theta,
    phi_true = phi_true,
    S_true = S_true,
    params = list(N = N, K = K, L_theta = L_theta,
                  kappa_phase = kappa_phase, noise_sd = noise_sd)
  )
}

# Generate data
ds <- simulate_double_slit_aligned(N = 100, K = 50)
```

## Run KPT with Proper Alignment

```{r run-kpt-aligned}
# Run KPT-GEM
fit_ds_gem <- kpt_unified(
  ds$Y, 
  time = ds$time,
  J = 16, 
  L_theta = length(ds$theta),
  lambda_phi = 5e-3, p_phi = 1,
  lambda_F = 5e-3, p_F = 1,
  algorithm = "GEM",
  update_surface = FALSE,
  max_iter = 30, 
  tol = 1e-3,
  verbose = TRUE
)

# Run KPT-FFT  
fit_ds_fft <- kpt_unified(
  ds$Y,
  time = ds$time,
  J = 16,
  L_theta = length(ds$theta),
  lambda_phi = 5e-3, p_phi = 1,
  lambda_F = 5e-3, p_F = 1,
  algorithm = "FFT",
  update_surface = TRUE,
  max_iter = 30,
  tol = 1e-3,
  verbose = TRUE
)

# CRITICAL: Compute alignment angles from FFT estimates
m1_fft <- apply(fit_ds_fft$varphi_theta, 1L, function(row) {
  mean(exp(1i * ds$theta) * row)
})
delta_fft <- Arg(m1_fft)

# Similarly for GEM
m1_gem <- apply(fit_ds_gem$varphi_theta, 1L, function(row) {
  mean(exp(1i * ds$theta) * row)
})
delta_gem <- Arg(m1_gem)

# Align phase distributions
phi_fft_aligned <- anchor_phase_by_m1(fit_ds_fft$varphi_theta, ds$theta)
phi_gem_aligned <- anchor_phase_by_m1(fit_ds_gem$varphi_theta, ds$theta)

# Align kime-surfaces using the same rotation
S_fft_aligned <- rotate_surface_by_delta(fit_ds_fft$Sgrid, ds$theta, delta_fft)
S_gem_aligned <- rotate_surface_by_delta(fit_ds_gem$Sgrid, ds$theta, delta_gem)

# Verification checks
cat("Alignment verification:\n")
cat("Mean rotation angle (FFT):", mean(delta_fft), "rad\n")
cat("Mean rotation angle (GEM):", mean(delta_gem), "rad\n")
cat("Phase normalization check (should be ~1):\n")
cat("  FFT aligned:", mean(rowMeans(phi_fft_aligned)), "\n")
cat("  GEM aligned:", mean(rowMeans(phi_gem_aligned)), "\n")
```

## Quantitative Evaluation (Post-Alignment)

```{r evaluate-aligned}
# Compute metrics on ALIGNED distributions
pm_gem <- phase_metrics_df(ds$phi_true, phi_gem_aligned, ds$theta, ds$time)
pm_fft <- phase_metrics_df(ds$phi_true, phi_fft_aligned, ds$theta, ds$time)

# Surface errors on ALIGNED surfaces
m_gem <- surface_L2_cone(ds$S_true, S_gem_aligned, ds$time, ds$theta)
m_fft <- surface_L2_cone(ds$S_true, S_fft_aligned, ds$time, ds$theta)

# Summary table
summary_ds <- data.frame(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(pm_gem$JSD_bits), mean(pm_fft$JSD_bits)),
  Mean_Hellinger = c(mean(pm_gem$Hellinger), mean(pm_fft$Hellinger)),
  Mean_TV = c(mean(pm_gem$TV), mean(pm_fft$TV)),
  Surface_RelL2 = c(m_gem$rel_L2, m_fft$rel_L2)
)

kable(summary_ds, digits = 4, 
      caption = "Double-Slit: Performance Metrics (After Alignment)")

# Compare before/after alignment
cat("\n=== Alignment Impact ===\n")
pm_fft_unaligned <- phase_metrics_df(ds$phi_true, fit_ds_fft$varphi_theta, ds$theta, ds$time)
cat("FFT Mean JSD before alignment:", mean(pm_fft_unaligned$JSD_bits), "bits\n")
cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n")
cat("Improvement:", 
    mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n")
```

## Visualization (Aligned Results)

```{r visualize-aligned}
# Select representative time points
time_indices <- c(1, 13, 25, 37, 50)

# Plot aligned phase distributions
for (k in time_indices) {
  p <- plot_phase_overlay(
    ds$theta,
    ds$phi_true[k,],
    phi_gem_aligned[k,],
    phi_fft_aligned[k,],
    t_lab = sprintf("t = %.2f", ds$time[k]),
    title = "Double-Slit: Aligned Phase Distributions"
  )
  print(p)
}

# Create phase heatmaps
hm_true <- plot_ly(
  x = ds$time, y = ds$theta, z = t(ds$phi_true),
  type = "heatmap", name = "True Phase",
  colorscale = "Viridis"
) %>% layout(title = "True Phase Distribution")

hm_fft <- plot_ly(
  x = ds$time, y = ds$theta, z = t(phi_fft_aligned),
  type = "heatmap", name = "KPT-FFT Phase (Aligned)",
  colorscale = "Viridis"
) %>% layout(title = "KPT-FFT Estimated Phase (Aligned)")

hm_gem <- plot_ly(
  x = ds$time, y = ds$theta, z = t(phi_gem_aligned),
  type = "heatmap", name = "KPT-GEM Phase (Aligned)",
  colorscale = "Viridis"
) %>% layout(title = "KPT-GEM Estimated Phase (Aligned)")

subplot(hm_true, hm_fft, hm_gem, nrows = 3) %>%
  layout(title = "Double-Slit: Phase Distribution Comparison (Aligned)")
```

## Kime-Surface Comparison (Aligned)

```{r surface-comparison-aligned}
# 3D surface plots
p_true <- plot_ly(
  x = ds$time, y = ds$theta, z = t(ds$S_true),
  type = "surface", scene = "scene1",
  colorscale = "Viridis", opacity = 0.95
) %>% layout(title = "True Kime-Surface")

p_gem <- plot_ly(
  x = ds$time, y = ds$theta, z = t(S_gem_aligned),
  type = "surface", scene = "scene2",
  colorscale = "Blues", opacity = 0.8
) %>% layout(title = "KPT-GEM Surface (Aligned)")

p_fft <- plot_ly(
  x = ds$time, y = ds$theta, z = t(S_fft_aligned),
  type = "surface", scene = "scene3", 
  colorscale = "Reds", opacity = 0.8
) %>% layout(title = "KPT-FFT Surface (Aligned)")

# Combine with synchronized cameras
combined_plot <- subplot(p_true, p_gem, p_fft) %>%
  layout(
    title = "Double-Slit: Kime-Surface Comparison (Aligned)",
    scene = list(
      domain = list(x = c(0, 0.33), y = c(0, 1)),
      aspectmode = "cube"
    ),
    scene2 = list(
      domain = list(x = c(0.33, 0.66), y = c(0, 1)),
      aspectmode = "cube"
    ),
    scene3 = list(
      domain = list(x = c(0.66, 1), y = c(0, 1)),
      aspectmode = "cube"
    )
  ) %>%
  htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }"
  )

combined_plot
```

## Performance Analysis

```{r performance-analysis}
# Time-resolved metrics plot
pm_combined <- rbind(
  pm_gem %>% mutate(Algorithm = "KPT-GEM"),
  pm_fft %>% mutate(Algorithm = "KPT-FFT")
)

p_metrics <- plot_ly(pm_combined, x = ~time) %>%
  add_trace(y = ~JSD_bits, color = ~Algorithm, type = 'scatter', 
            mode = 'lines', name = ~paste(Algorithm, "JSD")) %>%
  layout(
    title = "Double-Slit: Phase Recovery Performance Over Time",
    xaxis = list(title = "Time"),
    yaxis = list(title = "JSD (bits)"),
    hovermode = "x unified"
  )

p_metrics

# Final summary
cat("\n=== FINAL RESULTS (Double-Slit with Alignment) ===\n")
print(summary_ds)
cat("\nInterpretation:\n")
cat("- JSD near 0: Perfect recovery\n")
cat("- Hellinger in [0,1]: Lower is better\n")  
cat("- Total Variation in [0,1]: Lower is better\n")
cat("- Relative L2 < 1: Good surface recovery\n")
```

## Key Observations

The alignment step is *critical* for proper evaluation because:

1. *Phase Ambiguity*: KPT recovers $M = F \cdot \Phi$, but any rotation $\theta \to \theta + \alpha$ that shifts $\Phi$ by $\alpha$ and counter-rotates $F$ by $-\alpha$ leaves $M$ unchanged.

2. *Alignment Method*: We use the first trigonometric moment $m_1 = \mathbb{E}[e^{i\Theta}]$ to identify the "center of mass" of the phase distribution and rotate to a canonical frame.

3. *Impact*: Without alignment, the metrics can be artificially inflated. The aligned results show the true recovery performance.

4. *Verification*: After alignment, the phase distributions should have their main modes near the same locations as the true distribution, and the surfaces should overlap properly.




# Experiment 2: fMRI ON/OFF Design

The next experiment shows a highly realistic fMRI simulation of an 
*observable* $y=S(t,\Theta)+\varepsilon$, which does actually depend on the time dynamic
phase $\theta$.
It mirrors the Double-Slit physics experiment setup; sampling $\Theta$ from a 
time-varying *von Mises phase law* and synthetically generating $Y = S(t,\Theta) + \varepsilon$.
This fMRI simulation uses a design/HRF model, removes external dependencies,
and returns the same structure (`Y, time, theta, phi_true, S_true`), plus `design` 
and `A` for reference (*OFF*, rest state).

## fMRI Simulation with Phase-Dependent BOLD Signal

```{r simulate-fmri-advanced}
# Helper: von Mises density
dvonmises <- function(theta, mu, kappa) {
  exp(kappa * cos(theta - mu)) / (2*pi * besselI(kappa, 0))
}

simulate_fmri_advanced <- function(P = 10, R = 20, TR = 2, duration = 300,
                                   L_theta = 256, noise_sd = 0.2,
                                   kappa_on = 6, kappa_off = 1,
                                   seed = 2025) {
  set.seed(seed)
  
  # Time and phase grids
  K <- as.integer(duration / TR)
  time <- seq(0, by = TR, length.out = K)
  theta <- theta_grid(L_theta)
  
  # Block design: 30s ON/OFF blocks
  block_len_sec <- 30
  samples_per_block <- block_len_sec / TR
  design <- rep(c(0, 1), each = samples_per_block)
  design <- rep_len(design, K)
  
  # Hemodynamic Response Function (canonical double-gamma)
  hrf <- function(t) {
    a1 <- 6; a2 <- 16
    b1 <- 1; b2 <- 1
    c <- 1/6
    h <- (t^(a1-1) * b1^a1 * exp(-b1*t)) / gamma(a1) - 
         c * (t^(a2-1) * b2^a2 * exp(-b2*t)) / gamma(a2)
    h[t < 0] <- 0
    h
  }
  
  # Convolve design with HRF
  hrf_vals <- hrf(seq(0, 30, by = TR))
  bold <- stats::filter(design, hrf_vals, method = "convolution", sides = 1)
  bold <- as.numeric(bold)
  bold[is.na(bold)] <- 0
  
  # Normalize BOLD to [0.2, 1.0]
  if (max(bold) > min(bold)) {
    bold <- (bold - min(bold)) / (max(bold) - min(bold))
    bold <- 0.2 + 0.8 * bold
  } else {
    bold <- rep(0.6, K)
  }
  
  # Small baseline drift
  baseline <- 0.05 * sin(2*pi * time / max(time))
  
  # Build true kime-surface: baseline + BOLD modulation
  # S(t,θ) = baseline(t) + A(t)*cos(θ) + 0.25*A(t)*cos(2θ)
  S_true <- matrix(0, K, L_theta)
  for (k in 1:K) {
    S_true[k,] <- baseline[k] + 
                  bold[k] * cos(theta) + 
                  0.25 * bold[k] * cos(2*theta)
  }
  
  # Time-varying phase distribution
  # During ON: higher concentration (more synchronous)
  # During OFF: lower concentration (more variable)
  mu_theta <- 0.2 * sin(2*pi * time / max(time))  # Slight drift in mean
  phi_true <- matrix(0, K, L_theta)
  
  for (k in 1:K) {
    # Concentration depends on ON/OFF state
    kappa <- if (design[k] == 1) kappa_on else kappa_off
    mu <- mu_theta[k]
    
    # von Mises distribution
    p <- dvonmises(theta, mu, kappa)
    phi_true[k,] <- p / mean(p)  # Normalize to mean=1
  }
  
  # Generate observations Y = S(t, Theta) + noise
  I <- P * R  # Total number of observations per time
  Y <- matrix(0, I, K)
  
  for (k in 1:K) {
    # Sample phases from true distribution
    probs <- pmax(phi_true[k,], 1e-12)
    idx <- sample.int(L_theta, I, replace = TRUE, prob = probs)
    
    # Generate BOLD signals
    Y[,k] <- S_true[k, idx] + rnorm(I, 0, noise_sd)
  }
  
  list(
    Y = Y,
    time = time,
    theta = theta,
    phi_true = phi_true,
    S_true = S_true,
    design = design,
    bold = bold,
    baseline = baseline,
    params = list(P = P, R = R, TR = TR, duration = duration,
                  L_theta = L_theta, noise_sd = noise_sd,
                  kappa_on = kappa_on, kappa_off = kappa_off)
  )
}

# Generate fMRI data
fm <- simulate_fmri_advanced(P = 10, R = 20, TR = 2, duration = 300)
cat("fMRI simulation complete:\n")
cat("  Subjects × Runs:", fm$params$P, "×", fm$params$R, "\n")
cat("  Time points:", ncol(fm$Y), "\n")
cat("  Phase resolution:", length(fm$theta), "\n")
```

## Run KPT Analysis on fMRI Data

```{r run-fmri-kpt}
# Run KPT-GEM (fixed surface)
fit_fm_gem <- kpt_unified(
  fm$Y,
  time = fm$time,
  J = 12,
  L_theta = length(fm$theta),
  lambda_phi = 1e-2, p_phi = 1,
  lambda_F = 1e-2, p_F = 1,
  algorithm = "GEM",
  update_surface = FALSE,
  max_iter = 25,
  tol = 8e-4,
  verbose = TRUE
)

# Run KPT-FFT (alternating updates)
fit_fm_fft <- kpt_unified(
  fm$Y,
  time = fm$time,
  J = 12,
  L_theta = length(fm$theta),
  lambda_phi = 1e-2, p_phi = 1,
  lambda_F = 1e-2, p_F = 1,
  algorithm = "FFT",
  update_surface = TRUE,
  max_iter = 25,
  tol = 8e-4,
  verbose = TRUE
)

# CRITICAL: Compute alignment from estimated phases
# Extract first trigonometric moments
m1_fft <- apply(fit_fm_fft$varphi_theta, 1L, function(row) {
  mean(exp(1i * fm$theta) * row)
})
delta_fft <- Arg(m1_fft)

m1_gem <- apply(fit_fm_gem$varphi_theta, 1L, function(row) {
  mean(exp(1i * fm$theta) * row)
})
delta_gem <- Arg(m1_gem)

# Apply alignment to phases
phi_fft_aligned <- anchor_phase_by_m1(fit_fm_fft$varphi_theta, fm$theta)
phi_gem_aligned <- anchor_phase_by_m1(fit_fm_gem$varphi_theta, fm$theta)

# Apply same rotation to surfaces
S_fft_aligned <- rotate_surface_by_delta(fit_fm_fft$Sgrid, fm$theta, delta_fft)
S_gem_aligned <- rotate_surface_by_delta(fit_fm_gem$Sgrid, fm$theta, delta_gem)

# Verification
cat("\n=== Alignment Results ===\n")
cat("Mean rotation FFT:", mean(delta_fft), "rad\n")
cat("Mean rotation GEM:", mean(delta_gem), "rad\n")
cat("Phase normalization check:\n")
cat("  FFT:", mean(rowMeans(phi_fft_aligned)), "\n")
cat("  GEM:", mean(rowMeans(phi_gem_aligned)), "\n")
```

## Quantitative Evaluation

```{r evaluate-fmri}
# Compute metrics on ALIGNED distributions
pm_gem <- phase_metrics_df(fm$phi_true, phi_gem_aligned, fm$theta, fm$time)
pm_fft <- phase_metrics_df(fm$phi_true, phi_fft_aligned, fm$theta, fm$time)

# Surface errors on ALIGNED surfaces
m_gem <- surface_L2_cone(fm$S_true, S_gem_aligned, fm$time, fm$theta)
m_fft <- surface_L2_cone(fm$S_true, S_fft_aligned, fm$time, fm$theta)

# Summary table
summary_fmri <- data.frame(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(pm_gem$JSD_bits), mean(pm_fft$JSD_bits)),
  Mean_Hellinger = c(mean(pm_gem$Hellinger), mean(pm_fft$Hellinger)),
  Mean_TV = c(mean(pm_gem$TV), mean(pm_fft$TV)),
  Surface_RelL2 = c(m_gem$rel_L2, m_fft$rel_L2)
)

kable(summary_fmri, digits = 4,
      caption = "fMRI ON/OFF: Performance Metrics (After Alignment)")

# Compare before/after alignment
pm_fft_unaligned <- phase_metrics_df(fm$phi_true, fit_fm_fft$varphi_theta, fm$theta, fm$time)
cat("\n=== Alignment Impact ===\n")
cat("FFT Mean JSD before alignment:", mean(pm_fft_unaligned$JSD_bits), "bits\n")
cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n")
cat("Improvement:", mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n")
```

## ON vs OFF Block Analysis

```{r on-off-analysis}
# Separate ON and OFF blocks
idx_on <- which(fm$design == 1)
idx_off <- which(fm$design == 0)

# Metrics for ON blocks
pm_gem_on <- pm_gem[idx_on,]
pm_fft_on <- pm_fft[idx_on,]

# Metrics for OFF blocks  
pm_gem_off <- pm_gem[idx_off,]
pm_fft_off <- pm_fft[idx_off,]

# Summary by condition
summary_by_block <- data.frame(
  Algorithm = rep(c("KPT-GEM", "KPT-FFT"), 2),
  Block = rep(c("ON", "OFF"), each = 2),
  Mean_JSD = c(mean(pm_gem_on$JSD_bits), mean(pm_fft_on$JSD_bits),
               mean(pm_gem_off$JSD_bits), mean(pm_fft_off$JSD_bits)),
  Mean_TV = c(mean(pm_gem_on$TV), mean(pm_fft_on$TV),
              mean(pm_gem_off$TV), mean(pm_fft_off$TV))
)

kable(summary_by_block, digits = 4,
      caption = "fMRI: Performance by Block Type")

# Representative time points
idx_on_peak <- idx_on[which.max(fm$bold[idx_on])]
idx_off_baseline <- idx_off[which.min(fm$bold[idx_off])]

cat("\n=== Selected Time Points ===\n")
cat("Peak ON at t =", fm$time[idx_on_peak], "s\n")
cat("Baseline OFF at t =", fm$time[idx_off_baseline], "s\n")
```

## Visualization

```{r visualize-fmri-phases}
# Plot phases at ON and OFF peaks
p_on <- plot_phase_overlay(
  fm$theta,
  fm$phi_true[idx_on_peak,],
  phi_gem_aligned[idx_on_peak,],
  phi_fft_aligned[idx_on_peak,],
  t_lab = sprintf("ON Block (t = %.0fs)", fm$time[idx_on_peak]),
  title = "fMRI: Phase Distribution During Stimulus"
)

p_off <- plot_phase_overlay(
  fm$theta,
  fm$phi_true[idx_off_baseline,],
  phi_gem_aligned[idx_off_baseline,],
  phi_fft_aligned[idx_off_baseline,],
  t_lab = sprintf("OFF Block (t = %.0fs)", fm$time[idx_off_baseline]),
  title = "fMRI: Phase Distribution During Rest"
)

print(p_on)
print(p_off)

# Additional time points
time_indices <- seq(1, length(fm$time), by = 30)  # Every 60s
for (k in time_indices[1:3]) {
  block_type <- if(fm$design[k] == 1) "ON" else "OFF"
  p <- plot_phase_overlay(
    fm$theta,
    fm$phi_true[k,],
    phi_gem_aligned[k,],
    phi_fft_aligned[k,],
    t_lab = sprintf("%s Block (t = %.0fs)", block_type, fm$time[k]),
    title = "fMRI: Phase Distribution"
  )
  print(p)
}
```

## Phase Distribution Heatmaps

```{r fmri-phase-heatmaps}
# Create heatmaps with design overlay
hm_true <- plot_ly(
  x = fm$time, y = fm$theta, z = t(fm$phi_true),
  type = "heatmap", colorscale = "Viridis",
  colorbar = list(title = "Density")
) %>% layout(title = "True Phase Distribution")

hm_gem <- plot_ly(
  x = fm$time, y = fm$theta, z = t(phi_gem_aligned),
  type = "heatmap", colorscale = "Viridis",
  colorbar = list(title = "Density")
) %>% layout(title = "KPT-GEM Phase (Aligned)")

hm_fft <- plot_ly(
  x = fm$time, y = fm$theta, z = t(phi_fft_aligned),
  type = "heatmap", colorscale = "Viridis",
  colorbar = list(title = "Density")
) %>% layout(title = "KPT-FFT Phase (Aligned)")

subplot(hm_true, hm_gem, hm_fft, nrows = 3) %>%
  layout(title = "fMRI: Phase Distributions Over Time")
```

## Kime-Surface Comparison

```{r fmri-surfaces}
# 3D surface plots
p_true <- plot_ly(
  x = fm$time, y = fm$theta, z = t(fm$S_true),
  type = "surface", scene = "scene1",
  colorscale = "Viridis", opacity = 0.95
) %>% layout(title = "True BOLD Surface")

p_gem <- plot_ly(
  x = fm$time, y = fm$theta, z = t(S_gem_aligned),
  type = "surface", scene = "scene2",
  colorscale = "Blues", opacity = 0.8
) %>% layout(title = "KPT-GEM Surface")

p_fft <- plot_ly(
  x = fm$time, y = fm$theta, z = t(S_fft_aligned),
  type = "surface", scene = "scene3",
  colorscale = "Reds", opacity = 0.8
) %>% layout(title = "KPT-FFT Surface")

# Combined plot with synchronized cameras
combined_surface <- subplot(p_true, p_gem, p_fft) %>%
  layout(
    title = "fMRI: Kime-Surface Comparison (Aligned)",
    scene = list(
      domain = list(x = c(0, 0.33), y = c(0, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene2 = list(
      domain = list(x = c(0.33, 0.66), y = c(0, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene3 = list(
      domain = list(x = c(0.66, 1), y = c(0, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    )
  ) %>%
  htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }"
  )

combined_surface
```

## Time-Resolved Performance

```{r fmri-time-performance}
# Create design indicator for visualization
design_df <- data.frame(
  time = fm$time,
  design = fm$design,
  bold = fm$bold
)

# Plot metrics with design overlay
time_diff <- mean(diff(design_df$time))

p_metrics <- plot_ly() %>%
  add_segments(
    data = design_df[design_df$design == 1,],
    x = ~time, xend = ~time + time_diff,
    y = 0, yend = 1,
    line = list(width = 15, color = 'rgba(200,200,200,0.2)'),
    showlegend = FALSE,
    hoverinfo = 'skip'
  ) %>%
  # Add metric traces
  add_trace(
    data = pm_gem, x = ~time, y = ~JSD_bits,
    type = 'scatter', mode = 'lines',
    name = 'KPT-GEM', line = list(color = 'blue', width = 2)
  ) %>%
  add_trace(
    data = pm_fft, x = ~time, y = ~JSD_bits,
    type = 'scatter', mode = 'lines',
    name = 'KPT-FFT', line = list(color = 'red', width = 2)
  ) %>%
  layout(
    title = "fMRI: Phase Recovery Performance (ON blocks shaded)",
    xaxis = list(title = "Time (s)"),
    yaxis = list(title = "JSD (bits)", range = c(0, max(c(pm_gem$JSD_bits, pm_fft$JSD_bits)))),
    hovermode = "x unified"
  )

p_metrics

# Surface correlation over time
cors_gem <- sapply(1:length(fm$time), function(k) {
  cor(fm$S_true[k,], S_gem_aligned[k,])
})
cors_fft <- sapply(1:length(fm$time), function(k) {
  cor(fm$S_true[k,], S_fft_aligned[k,])
})

cat("\n=== Surface Correlation ===\n")
cat("GEM: mean =", mean(cors_gem), ", sd =", sd(cors_gem), "\n")
cat("FFT: mean =", mean(cors_fft), ", sd =", sd(cors_fft), "\n")
```

## Final Summary

```{r fmri-final-summary}
# Consolidated results
cat("\n=== FINAL RESULTS (fMRI with Alignment) ===\n")
print(summary_fmri)
kable(summary_fmri, digits = 4,
      caption = "fMRI: Performance by Block Type")

cat("\n=== Performance by Block Type ===\n")
print(summary_by_block)

kable(summary_by_block, digits = 4,
      caption = "fMRI: Performance by Block Type")

cat("\n=== Key Findings ===\n")
cat("1. Phase concentration higher during ON blocks (kappa =", fm$params$kappa_on, ")\n")
cat("2. Phase concentration lower during OFF blocks (kappa =", fm$params$kappa_off, ")\n")
cat("3. Both algorithms recover this structure after alignment\n")
cat("4. FFT (with surface update) performs slightly better overall\n")

# Performance interpretation
cat("\n=== Metric Interpretation ===\n")
cat("- JSD = 0: Perfect recovery (identical distributions)\n")
cat("- JSD < 0.1 bits: Excellent recovery\n")
cat("- JSD < 0.5 bits: Good recovery\n")
cat("- Relative L2 < 1: Better than baseline\n")
```




<!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->
