---
title: "Time Complexity, Inferential uncertainty and Spacekime Analytics"
subtitle: "Appendix 03: Kime-Phase Tomography (KPT V.4: 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)
```

# Overview & Conventions

This notebook demonstrates *KPT-GEM* and *KPT-FFT* on two simulations:

 - *Double-slit* physics experiment (latent phase drives interference).

 - *fMRI ON vs. OFF* event-related design (latent phase modulates BOLD fMRI neuroscience study).

We use the circular phase normalization 
$\int_0^{2\pi}\varphi_t(\theta)\,\frac{d\theta}{2\pi}=1$. On a uniform grid of length 
$L_\theta$, this is equivalent to the *row-mean equals 1* constraint. The spectral updates
use the [Wiener-Sobolev filter](https://doi.org/10.1016/S0022-1236(03)00236-2).

$$\widehat{\Phi}(\omega,t)=\frac{\overline{\widehat F(\omega,t)}}{|\widehat F(\omega,t)|^2+\lambda\,\Lambda_p(\omega)}\,\widehat M(\omega,t),\quad \Lambda_p(\omega)=(2\sin(\omega/2))^{2p},$$

with simplex projection in $\theta$ (nonnegativity + unit mass).

*Algorithmic parity.* KPT-FFT is an alternating spectral method (phase-step *and*
surface-step) implemented with `fft()` for the DFT/IDFT. KPT-GEM is the same 
penalized tomography specialized to updating $\varphi_t$ given a stabilized 
surface model $\mathcal S$ (practical when a robust $\mathcal S$ proxy is available
from the data, e.g., ON-OFF construction).

# Utilities (Fourier, projection, metrics, visuals)

We start by defining the core utility functions.

```{r utilities}
# --- In Utilities chunk (right before defining plotting helpers) ---
if (exists("plot_phase_overlay", inherits = TRUE)) {
  rm(plot_phase_overlay)
}

# Dimension-checked overlay (unique name to avoid conflicts)
kpt_plot_phase_overlay <- function(theta, phi_true_row, phi_est_row, phi_est_row2=NULL, t_lab = "", title = "") {
  n_theta <- length(theta)
  stopifnot(length(phi_true_row) == n_theta, length(phi_est_row) == n_theta)

  if (is.null(phi_est_row2)) {
    df <- tibble::tibble(
      theta   = rep(theta, 2L),                            # length = 2*n
      density = c(phi_true_row, phi_est_row),              # length = 2*n
      kind    = rep(c("True", "Estimated"), each = n_theta) # length = 2*n
    )
  } else {
    df <- tibble::tibble(
      theta   = rep(theta, 3L),                            # length = 3*n
      density = c(phi_true_row, phi_est_row, phi_est_row2),              # length = 3*n
      kind    = rep(c("True", "FFT-Estimated", "GEM-Estimated"), each = n_theta) # length = 3*n
    )
  }

  ggplot2::ggplot(df, ggplot2::aes(theta, density, color = kind)) +
    ggplot2::geom_line(linewidth = 0.9) +
    ggplot2::labs(title = title, subtitle = t_lab,
                  x = expression(theta), y = expression(varphi(t,theta))) +
    ggplot2::theme_minimal()
}

# ---- Grids & helpers ----
theta_grid <- function(L) seq(0, 2*pi, length.out = L + 1L)[1L:L]
make_nvec   <- function(J) seq.int(-J, J, by = 1L)

# Work with uniform grid on [0,2π)
wrap_0_2pi <- function(x) (x %% (2*pi))

# Posterior softmax (row-wise, stable)
softmax_rows <- function(Lmat) {
  m <- apply(Lmat, 1L, max)
  W <- exp(Lmat - m)
  W / pmax(rowSums(W), .Machine$double.eps)
}

# Enforce probability simplex per row on uniform theta grid:
#   nonnegativity + ∫ φ dθ/(2π)=1  <=> row mean equals 1
row_project_simplex <- function(P, eps = 1e-12) {
  P <- pmax(P, 0)
  m <- rowMeans(P)
  m[m < eps] <- 1
  sweep(P, 1L, m, `/`)
}

# ---- FFT-friendly mapping between coefficients (n=-J..J) and DFT order (m=0..N-1) ----
coeffs_to_dft_order <- function(x_n, nvec) {
  N <- length(nvec); stopifnot(N == length(x_n))
  x_m <- complex(length = N)
  for (i in seq_along(nvec)) {
    n <- nvec[i]
    m <- (n + N) %% N          # map -J..J -> 0..N-1
    x_m[m + 1L] <- x_n[i]
  }
  x_m
}
dft_order_to_coeffs <- function(x_m, nvec) {
  N <- length(nvec); stopifnot(N == length(x_m))
  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
}

# Evaluate generating function on unit circle via FFT
coeffs_to_Xomega_fft <- function(x_n, nvec) {
  x_m <- coeffs_to_dft_order(x_n, nvec)
  fft(x_m, inverse = FALSE)    # no 1/N scaling
}
# Invert: X(ω) -> coefficients in n=-J..J via IFFT
Xomega_to_coeffs_fft <- function(X_w, nvec) {
  N <- length(nvec); stopifnot(N == length(X_w))
  x_m <- fft(X_w, inverse = TRUE) / N
  dft_order_to_coeffs(x_m, nvec)
}

# Synthesize S(t,θ) from f_n(t) on θ-grid
surface_from_coeffs <- function(f_n_row, nvec, theta) {
  # S(t,θ) = Re(Σ_n f_n e^{i n θ})
  E <- exp(1i * outer(nvec, theta))
  as.vector(Re(drop(f_n_row %*% E)))
}

# Reality constraint on Fourier coeffs: f_-n = Conj(f_n)
enforce_reality_coeffs <- function(f_n, nvec) {
  N <- length(nvec); J <- (N - 1L) / 2L
  # indices: n=0..J and negative mates
  for (j in 1:J) {
    pos <- which(nvec ==  j)
    neg <- which(nvec == -j)
    avg <- 0.5 * (f_n[pos] + Conj(f_n[neg]))
    f_n[pos] <- avg
    f_n[neg] <- Conj(avg)
  }
  # n=0 coefficient should be real
  f_n[nvec==0] <- Re(f_n[nvec==0])
  f_n
}

# ---- Phase metrics per time (KL, JSD, Hellinger, TV) ----
row_normalize_phi <- function(P, eps = 1e-12) {
  P <- pmax(P, eps); m <- rowMeans(P); m[m < eps] <- 1; sweep(P, 1, m, `/`)
}
phase_metrics_df <- function(phi_true, phi_est, theta, time, base2 = TRUE) {
  P <- row_normalize_phi(phi_true); Q <- row_normalize_phi(phi_est)
  K <- nrow(P); L <- ncol(P); stopifnot(K == nrow(Q), L == ncol(Q))
  logb <- if (base2) function(x) log(x, base = 2) else log
  tibble(
    time = time,
    KL_true_est = double(K),
    KL_est_true = double(K),
    JSD_bits    = double(K),
    Hellinger   = double(K),
    TV          = double(K)
  ) %>%
    mutate(
      KL_true_est = rowMeans(P * logb(P / pmax(Q, 1e-12))),
      KL_est_true = rowMeans(Q * logb(Q / pmax(P, 1e-12))),
      {
        M <- 0.5 * (P + Q)
        JSD_bits <- 0.5 * rowMeans(P * logb(P / pmax(M, 1e-12))) +
                    0.5 * rowMeans(Q * logb(Q / pmax(M, 1e-12)))
        Hellinger <- sqrt(0.5 * rowMeans((sqrt(P) - sqrt(Q))^2))
        TV        <- 0.5 * rowMeans(abs(P - Q))
        data.frame(JSD_bits, Hellinger, TV)
      }
    )
}

# ---- Surface L2 (cone-weighted) ----
surface_L2_cone <- function(S_true, S_est, time, theta) {
  stopifnot(all(dim(S_true) == dim(S_est)))
  K <- nrow(S_true); L <- ncol(S_true)
  dtheta <- 2*pi / L
  dt <- numeric(K)
  if (K == 1) dt[] <- 1 else {
    dt[1] <- (time[2]-time[1])/2; dt[K] <- (time[K]-time[K-1])/2
    if (K > 2) dt[2:(K-1)] <- (time[3:K]-time[1:(K-2)])/2
  }
  w_t <- time
  diff2 <- (S_true - S_est)^2
  per_t_L2_sq <- rowSums(diff2) * dtheta * w_t
  total_L2 <- sqrt(sum(per_t_L2_sq * dt))
  true_norm_sq <- rowSums(S_true^2) * dtheta * w_t
  rel_L2 <- total_L2 / max(sqrt(sum(true_norm_sq * dt)), .Machine$double.eps)
  list(total_L2 = total_L2, rel_L2 = rel_L2,
       per_time_L2 = sqrt(rowSums(diff2) * dtheta * w_t))
}

# ---- Visualization helpers ----
plot_phase_overlay <- function(theta, phi_true_row, phi_est_row, t_lab = "", title = "") {
  df <- tibble(
    theta = rep(theta, 2),
    density = c(phi_true_row, phi_est_row),
    kind = rep(c("True","Estimated"), each = length(theta))
  )
  ggplot(df, aes(theta, density, color = kind)) +
    geom_line(linewidth = 0.9) +
    labs(title = title, subtitle = t_lab, x = expression(theta), y = expression(varphi(t,theta))) +
    theme_minimal()
}
# plot_surface_3d <- function(time, theta, S, title = "Surface", opacity=0.9, colorscale='Viridis') {
#   plot_ly(x = time, y = theta, z = S, type = "surface", opacity=opacity, 
#           colorscale=colorscale, name=title) %>%
#     layout(title = title,
#            scene = list(xaxis = list(title = "t"),
#                         yaxis = list(title = "θ"), # expression(theta)),
#                         zaxis = list(title = "Intensity")))
# }
plot_surface_3d <- function(time, theta, S, title = "Surface", opacity=0.9, 
                           colorscale='Viridis', colorbar_title = "Intensity",
                           showscale = TRUE) {
  plot_ly(x = time, y = theta, z = S, type = "surface", opacity=opacity, 
          colorscale=colorscale, 
          colorbar = list(title = colorbar_title),
          showscale = showscale) %>%
    layout(title = title,
           scene = list(xaxis = list(title = "t"),
                        yaxis = list(title = "θ"),
                        zaxis = list(title = colorbar_title)))
}
```

We need to avoid potential phase-periodicity inconsistencies or phase-shift ambiguity (anchoring)
between phase domains $[-\pi,\pi)$ and $[0,2\pi)$. As KPT recovers 
$M=F\Phi$, if  $\Phi(\cdot,t)$ is rotated by $\alpha(t)$ and $F(\cdot,t)$ is counter-rotated
by $-\alpha(t)$, then the product doesn’t change. Without an anchoring rule,
$\hat{\varphi}^t(\theta)$ can drift by a constant rotation $\theta \to \theta + \alpha$
making the performance-validation plots shifted (not necessarily by $\pi$, just by some 
$\alpha$. This is expected and harmless unless we compare the estimated phases to
to the fixed (gold-standard) ground truth phases.

# KPT Algorithms

All methematical, computaitonal, and statistical details are provided in the 
[KPT paper](https://www.overleaf.com/project/68a72229f8e907137629d964).

## KPT-GEM (penalized tomography with stabilized surface)

This KPT algorithmic variant (*KPT-GEM*, generalized expectation maximization estimation) 
updates $\widehat{\varphi}_t$ via Wiener-Sobolev filtering given a stable proxy 
$\mathcal S(t,\theta)$ (from the data), then projects to the simplex.

### Kime-Phase Tomography generalized EM algorithm (KPT-GEM) Algorithm

**Inputs:** $\{y_{j,k}\}$, harmonics $n\in\{-J,\ldots,J\}$, tolerance $\varepsilon$, $(\lambda_\Phi,p_\Phi)$, $(\lambda_F,p_F)$

**Initialize:** $\widehat\varphi^{(0)}_{t_k}\equiv \tfrac{1}{2\pi}$; initialize $\widehat f^{(0)}_n(t_k)$ (e.g., from cross‑sectional averages); set iteration $\ell=0$

**Repeat:**

1.  **E-step:** compute $\xi^{(\ell)}_{j,k}(n)$ and $\widehat m^{(\ell)}_n(t_k)=\tfrac{1}{N}\sum_j y_{j,k}\xi^{(\ell)}_{j,k}(n)$

2.  **$\Phi$-step (phase update):** For each $t_k$, FFT in $n$ to get $\widehat M^{(\ell)}(\omega,t_k)$ and $\widehat F^{(\ell)}(\omega,t_k)$; set
    
    $$\widehat\Phi^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat F^{(\ell)}(\omega,t_k)}}{|\widehat F^{(\ell)}(\omega,t_k)|^2+\lambda_\Phi\,\Lambda_{p_\Phi}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k).$$
    
    IFFT in $\omega$ and project: $\widehat\varphi^{(\ell+1)}_{t_k}\leftarrow \Pi_{\mathcal P}\big[\mathcal{F}^{-1}\widehat\Phi^{(\ell+1)}(\cdot,t_k)\big]$. Optionally smooth in $t$.

3.  **$F$-step (surface update):** For each $t_k$, recompute $\widehat\Phi^{(\ell+1)}$ and set
    
    $$\widehat F^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat \Phi^{(\ell+1)}(\omega,t_k)}}{|\widehat \Phi^{(\ell+1)}(\omega,t_k)|^2+\lambda_F\,\Lambda_{p_F}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k).$$
    
    IFFT in $\omega$ to obtain $\widehat f^{(\ell+1)}_n(t_k)$; enforce $f_{-n}=\overline{f_n}$. Optionally smooth in $t$.

4.  **Synthesis:** $\widehat{\mathcal S}^{(\ell+1)}(t_k,\theta)=\sum_{n=-J}^J \widehat f^{(\ell+1)}_n(t_k)e^{in\theta}$

5.  **Convergence check:** stop when $\max_k\|\widehat\varphi^{(\ell+1)}_{t_k}-\widehat\varphi^{(\ell)}_{t_k}\|_{L^2}<\varepsilon$ and $\max_{k,n}|\widehat f^{(\ell+1)}_n(t_k)-\widehat f^{(\ell)}_n(t_k)|<\varepsilon$

**Until** converged

<!-- \begin{algorithm}[Penalized Alternating Minimization (KPT-GEM)] -->
<!-- \begin{algorithmic}[1] -->
<!-- \caption{Kime-Phase Tomography generalized EM algorithm (KPT-GEM)}\label{alg:KPT_EM} -->
<!-- \STATE \textbf{Inputs:} $\{y_{j,k}\}$, harmonics $n\in\{-J,\ldots,J\}$, tolerance $\varepsilon$, $(\lambda_\Phi,p_\Phi)$, $(\lambda_F,p_F)$ -->
<!-- \STATE \textbf{Initialize:} $\widehat\varphi^{(0)}_{t_k}\equiv \tfrac{1}{2\pi}$; initialize $\widehat f^{(0)}_n(t_k)$ (e.g., from cross‑sectional averages); set iteration $\ell=0$ -->
<!-- \REPEAT -->
<!-- \STATE \textbf{E-step:} compute $\xi^{(\ell)}_{j,k}(n)$ and $\widehat m^{(\ell)}_n(t_k)=\tfrac{1}{N}\sum_j y_{j,k}\xi^{(\ell)}_{j,k}(n)$ -->
<!-- \STATE \textbf{\(\Phi\)-step (phase update):} For each $t_k$, FFT in $n$ to get $\widehat M^{(\ell)}(\omega,t_k)$ and $\widehat F^{(\ell)}(\omega,t_k)$; set -->
<!-- \[ -->
<!-- \widehat\Phi^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat F^{(\ell)}(\omega,t_k)}}{|\widehat F^{(\ell)}(\omega,t_k)|^2+\lambda_\Phi\,\Lambda_{p_\Phi}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k). -->
<!-- \] -->
<!-- IFFT in $\omega$ and project: $\widehat\varphi^{(\ell+1)}_{t_k}\leftarrow \Pi_{\mathcal P}\big[\mathcal{F}^{-1}\widehat\Phi^{(\ell+1)}(\cdot,t_k)\big]$. Optionally smooth in $t$. -->
<!-- \STATE \textbf{\(F\)-step (surface update):} For each $t_k$, recompute $\widehat\Phi^{(\ell+1)}$ and set -->
<!-- \[ -->
<!-- \widehat F^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat \Phi^{(\ell+1)}(\omega,t_k)}}{|\widehat \Phi^{(\ell+1)}(\omega,t_k)|^2+\lambda_F\,\Lambda_{p_F}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k). -->
<!-- \] -->
<!-- IFFT in $\omega$ to obtain $\widehat f^{(\ell+1)}_n(t_k)$; enforce $f_{-n}=\overline{f_n}$. Optionally smooth in $t$. -->
<!-- \STATE \textbf{Synthesis:} $\widehat{\mathcal S}^{(\ell+1)}(t_k,\theta)=\sum_{n=-J}^J \widehat f^{(\ell+1)}_n(t_k)e^{in\theta}$ -->
<!-- \STATE \textbf{Convergence check:} stop when $\max_k\|\widehat\varphi^{(\ell+1)}_{t_k}-\widehat\varphi^{(\ell)}_{t_k}\|_{L^2}<\varepsilon$ and $\max_{k,n}|\widehat f^{(\ell+1)}_n(t_k)-\widehat f^{(\ell)}_n(t_k)|<\varepsilon$ -->
<!-- \UNTIL{converged} -->
<!-- \end{algorithmic} -->
<!-- \end{algorithm} -->

```{r kpt-gem}
build_S_from_data <- function(Y, time, theta, alpha2 = 0.2, smooth_spar = 0.6, offset_frac = 0.1) {
  # Robust ON-OFF style proxy: baseline + amplitude * (cos θ + α2 cos 2θ)
  K <- ncol(Y); L <- length(theta)
  mu_k <- apply(Y, 2L, mean)
  A_k  <- apply(Y, 2L, sd)
  mu_s <- stats::smooth.spline(time, mu_k, spar = smooth_spar)$y
  A_s  <- stats::smooth.spline(time, A_k,  spar = smooth_spar)$y
  off  <- offset_frac * median(abs(A_s) + 1e-8)
  cos1 <- cos(outer(rep(1, K), theta))
  cos2 <- cos(2 * outer(rep(1, K), theta))
  S <- mu_s %o% rep(1, L) + (A_s + off) %o% rep(1, L) * cos1 + (alpha2 * A_s) %o% rep(1, L) * cos2
  S
}

kpt_gem <- function(Y, time, J = 12, L_theta = 256, lambda = 1e-2, p = 1,
                    max_iter = 30, tol = 1e-3, sigma = NULL, verbose = TRUE) {
  I <- nrow(Y); K <- ncol(Y)
  theta <- theta_grid(L_theta); nvec <- make_nvec(J); Nn <- length(nvec)
  if (is.null(sigma)) sigma <- 0.5 * stats::mad(as.numeric(Y), center = 0) + 1e-6
  Sgrid <- build_S_from_data(Y, seq_len(K), theta)     # proxy surface on time index (unit steps)

  # Precompute f_n(t_k) by IFFT across θ for each time (DFT by matrix for clarity)
  f_n <- matrix(0+0i, nrow = K, ncol = Nn)
  for (k in 1:K) {
    Eneg <- exp(-1i * outer(theta, nvec))      # [Lθ x (2J+1)]
    f_n[k, ] <- as.vector((t(Sgrid[k, ]) %*% Eneg) / L_theta)
    f_n[k, ] <- enforce_reality_coeffs(f_n[k, ], nvec)
  }

  phi_grid <- matrix(1.0, nrow = K, ncol = L_theta)
  Lambda_p <- (2 * sin(pi * (0:(Nn-1))/Nn))^(2*p)  # Λ_p(ω) with ω=2πk/Nn
  obj <- numeric(max_iter)

  for (it in 1:max_iter) {
    # E-step: posterior weights & mixed moments
    m_hat <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      prior_k <- phi_grid[k, ]
      Srow <- matrix(Sgrid[k, ], nrow = I, ncol = L_theta, byrow = TRUE)
      LogL <- - ( (Y[, k] - Srow)^2 ) / (2 * sigma^2)
      LogW <- sweep(LogL, 2L, log(pmax(prior_k, 1e-12)), `+`)
      W <- softmax_rows(LogW)                              # [I x Lθ]
      Eneg <- exp(-1i * outer(theta, nvec))                # [Lθ x (2J+1)]
      xi <- W %*% Eneg                                     # [I x (2J+1)]
      m_hat[k, ] <- colMeans(xi * Y[, k])
    }

    # M-step: spectral deconvolution per time using FFT evaluation on unit circle
    phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      F_w <- coeffs_to_Xomega_fft(f_n[k, ], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(F_w)^2) + lambda * Lambda_p + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      phi_hat_n[k, ] <- Xomega_to_coeffs_fft(Phi_w, nvec)
    }
    # Synthesize φ and project
    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)

    diff_norm <- sqrt(mean((phi_new - phi_grid)^2))
    phi_grid <- phi_new
    obj[it] <- mean(Mod(m_hat)^2)
    if (verbose) message(sprintf("[GEM] iter %02d: Δ=%.3e", it, diff_norm))
    if (diff_norm < tol) { obj <- obj[1:it]; break }
  }

  list(varphi_theta = phi_grid, phi_hat_n = phi_hat_n,
       f_n = f_n, Sgrid = Sgrid, theta = theta, nvec = nvec,
       omega = 2*pi*(0:(Nn-1))/Nn, time = seq_len(K), lambda = lambda, p = p,
       sigma = sigma, objective = obj)
}
```

## KPT-FFT (alternating phase & surface updates)

This is the fully alternating spectral algorithm consistent with the $\Phi$-step / $F$-step updates

$$\Phi \leftarrow \frac{\overline F}{|F|^2+\lambda_\Phi\Lambda_{p_\Phi}}\,M,\qquad
F \leftarrow \frac{\overline \Phi}{|\Phi|^2+\lambda_F\Lambda_{p_F}}\,M,$$

with IFFT to coefficients and synthesis after each step (plus reality constraint for $f_n$).

## Algorithm: Fast Kime‑Phase Tomography via FFT (KPT-FFT)

**Input:** $\{y_{j,k}\}$, $J$, $(\lambda_\Phi,p_\Phi)$, $(\lambda_F,p_F)$, tolerance $\epsilon$

**Init:** $\widehat\varphi^{(0)}_{t_k}\equiv\tfrac{1}{2\pi}$; initialize $\widehat f^{(0)}_n(t_k)$

**Repeat:**

1.  **E-step (vectorized):** compute $\xi^{(\ell)}_{j,k}(n)$ and $\widehat m^{(\ell)}_n(t_k)$
2.  **FFT:** $\widehat M^{(\ell)}(\omega,t_k)\leftarrow \mathrm{FFT}_n\{\widehat m^{(\ell)}_n(t_k)\}$
3.  **Phase update:** 
    $$\widehat\Phi^{(\ell+1)}(\omega,t_k)\leftarrow \dfrac{\overline{\widehat F^{(\ell)}(\omega,t_k)}}{|\widehat F^{(\ell)}(\omega,t_k)|^2+\lambda_\Phi\Lambda_{p_\Phi}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k)$$
4.  **IFFT & projection:** $\widehat\varphi^{(\ell+1)}_{t_k}\leftarrow \Pi_{\mathcal P}[\mathrm{IFFT}_\omega\{\widehat\Phi^{(\ell+1)}(\omega,t_k)\}]$
5.  **Surface update:** 
    $$\widehat F^{(\ell+1)}(\omega,t_k)\leftarrow \dfrac{\overline{\widehat \Phi^{(\ell+1)}(\omega,t_k)}}{|\widehat \Phi^{(\ell+1)}(\omega,t_k)|^2+\lambda_F\Lambda_{p_F}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k)$$
6.  **IFFT & reality:** $\widehat f^{(\ell+1)}_n(t_k)\leftarrow \mathrm{IFFT}_\omega\left\{\widehat F^{(\ell+1)}(\omega,t_k)\right\}$ and enforce $f_{-n}=\overline{f_n}$

**Until** $\|\widehat\varphi^{(\ell+1)}-\widehat\varphi^{(\ell)}\|_{L^2}<\epsilon$ and $\max_{k,n}|\widehat f^{(\ell+1)}_n(t_k)-\widehat f^{(\ell)}_n(t_k)|<\epsilon$

**Return** $\widehat\varphi_t$ and $\widehat{\mathcal S}(t,\theta)$.

```{r kpt-fft}
kpt_fft <- function(Y, time,
                    J = 12, L_theta = 256,
                    lambda_phi = 1e-2, p_phi = 1,
                    lambda_F   = 1e-2, p_F   = 1,
                    max_iter = 30, tol = 1e-3, sigma = NULL, verbose = TRUE) {
  I <- nrow(Y); K <- ncol(Y)
  theta <- theta_grid(L_theta); nvec <- make_nvec(J); Nn <- length(nvec)
  omega <- 2*pi*(0:(Nn-1))/Nn
  if (is.null(sigma)) sigma <- 0.5 * stats::mad(as.numeric(Y), center = 0) + 1e-6

  # Initialize: uniform φ, simple S proxy
  phi_grid <- matrix(1.0, nrow = K, ncol = L_theta)
  Sgrid <- build_S_from_data(Y, seq_len(K), theta)
  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)
  }

  Lambda_phi <- (2 * sin(0.5 * omega))^(2*p_phi)
  Lambda_F   <- (2 * sin(0.5 * omega))^(2*p_F)

  obj <- numeric(max_iter)

  for (it in 1:max_iter) {
    # ---- E-step (posterior stats & mixed moments) ----
    m_hat <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      prior_k <- phi_grid[k, ]
      Srow <- matrix(Sgrid[k, ], nrow = I, ncol = L_theta, byrow = TRUE)
      LogL <- - ( (Y[, k] - Srow)^2 ) / (2 * sigma^2)
      LogW <- sweep(LogL, 2L, log(pmax(prior_k, 1e-12)), `+`)
      W <- softmax_rows(LogW)
      Eneg <- exp(-1i * outer(theta, nvec))
      xi <- W %*% Eneg
      m_hat[k, ] <- colMeans(xi * Y[, k])
    }

    # ---- Φ-step (phase update in spectral domain via FFT) ----
    phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      F_w <- coeffs_to_Xomega_fft(f_n[k, ], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(F_w)^2) + lambda_phi * Lambda_phi + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      phi_hat_n[k, ] <- Xomega_to_coeffs_fft(Phi_w, nvec)
    }
    # Synthesize φ(θ) & project
    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 (surface update in spectral domain via FFT) ----
    f_n_new <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      Phi_w <- coeffs_to_Xomega_fft(phi_hat_n[k, ], nvec)
      M_w   <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(Phi_w)^2) + lambda_F * Lambda_F + 1e-10
      F_w   <- Conj(Phi_w) * M_w / denom
      f_n_new[k, ] <- enforce_reality_coeffs(Xomega_to_coeffs_fft(F_w, nvec), nvec)
    }
    # Synthesize S(t,θ)
    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)

    # ---- Convergence & update ----
    diff_phi <- sqrt(mean((phi_new - phi_grid)^2))
    diff_f   <- sqrt(mean(Mod(f_n_new - f_n)^2))
    if (verbose) message(sprintf("[FFT] iter %02d: Δφ=%.3e, Δf=%.3e", it, diff_phi, diff_f))
    phi_grid <- phi_new; f_n <- f_n_new; Sgrid <- S_new
    obj[it] <- mean(Mod(m_hat)^2)
    if (max(diff_phi, diff_f) < tol) { obj <- obj[1:it]; break }
  }

  list(varphi_theta = phi_grid, phi_hat_n = phi_hat_n,
       f_n = f_n, Sgrid = Sgrid, theta = theta, nvec = nvec, omega = omega,
       time = seq_len(K),
       lambda_phi = lambda_phi, p_phi = p_phi, lambda_F = lambda_F, p_F = p_F,
       sigma = sigma, objective = obj)
}
```

# Simulation 1: Double-Slit Physics Experiment

We simulate an interference-like kime-surface

$$\mathcal S(t,\theta)=B(t)+A_1(t)\cos(\kappa t+\theta)+A_2(t)\cos(2(\kappa t+\theta)+\phi_2),$$
e.g., $\mu(t) + A(t)\cos(\theta) + 0.25A(t)\cos(2\theta)$,
with a time-varying *mixture of von Mises* phase laws (two lobes drifting over time), then generate repeated measurements
$$\underbrace{y_{j,k}}_{observed}=\overbrace{\mathcal 
S(t_k,\underbrace{\Theta_j(t_k)}_{intrinsic\\ phase\ domain\\ stochasticity})}^{kimesurface}
+\underbrace{\varepsilon_{j,k}}_{extrinsic\ range\\ stochasticity}.$$

```{r sim-double-slit}
# von Mises pdf on grid (no external deps)
dvonmises <- function(theta, mu, kappa) {
  (1/(2*pi* besselI(kappa, 0))) * exp(kappa * cos(theta - mu))
}

simulate_double_slit <- function(I = 256, K = 256, L_theta = 256,
                                 kappa_phase = 6, noise_sd = 0.20,
                                 kappa_wave = 12, seed = 42) {
  set.seed(seed)
  theta <- theta_grid(L_theta)
  t <- seq(0, 1, length.out = K)

  # Surface components
  B  <- 0.2 * sin(2*pi*t)                  # small baseline
  A1 <- 0.8 * exp(-20*(t-0.5)^2)           # Gaussian envelope
  A2 <- 0.3 * exp(-20*(t-0.5)^2)
  phi2 <- pi/6

  S_true <- sapply(t, function(tt)
    B[which.min(abs(t-tt))] +
      A1[which.min(abs(t-tt))] * cos(kappa_wave*tt + theta) +
      A2[which.min(abs(t-tt))] * cos(2*(kappa_wave*tt + theta) + phi2)
  ) %>% t()

  # Time-varying von Mises mixture over theta (two drifting lobes)
  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)             # mixture weight
  phi_true <- matrix(0, nrow = K, ncol = L_theta)
  for (k in 1:K) {
    p <- w[k]*dvonmises(theta, mu1[k], kappa_phase) +
         (1-w[k])*dvonmises(theta, mu2[k], kappa_phase)
    phi_true[k, ] <- p / mean(p)           # mean == 1 (∫dθ/(2π)=1)
  }

  # Draw Θ_jk by discrete sampling on theta grid according to phi_true[k,]
  sample_theta <- function(p_row) {
    sample(theta, size = I, replace = TRUE, prob = pmax(p_row, 1e-12))
  }
  Y <- matrix(NA_real_, nrow = I, ncol = K)
  for (k in 1:K) {
    thetas <- sample_theta(phi_true[k, ])
    mu_k <- sapply(thetas, function(th) S_true[k, which.min(abs(theta - th))])
    Y[, k] <- mu_k + rnorm(I, 0, noise_sd)
  }
  list(Y = Y, time = t, theta = theta, phi_true = phi_true, S_true = S_true)
}

ds <- simulate_double_slit()
str(ds, 1)
```

## Testing KPT-GEM and KPT-FFT on Double-Slit Experiment

```{r run-ds}
# KPT-GEM (surface stabilized from data)
fit_ds_gem <- kpt_gem(ds$Y, time = ds$time, J = 16, L_theta = length(ds$theta),
                      lambda = 2e-2, p = 1, max_iter = 30, tol = 5e-4, verbose = TRUE)

# KPT-FFT (alternating)
fit_ds_fft <- kpt_fft(ds$Y, time = ds$time, J = 16, L_theta = length(ds$theta),
                      lambda_phi = 2e-2, p_phi = 1,
                      lambda_F   = 2e-2, p_F   = 1,
                      max_iter = 30, tol = 5e-4, verbose = TRUE)
```

## Double-Slit Results: Quantitative & Qualitative Evaluation

```{r eval-ds}
# Phase metrics per time
pm_gem <- phase_metrics_df(ds$phi_true, fit_ds_gem$varphi_theta, ds$theta, ds$time)
pm_fft <- phase_metrics_df(ds$phi_true, fit_ds_fft$varphi_theta, ds$theta, ds$time)

# Surface L2 under cone measure
S_gem <- fit_ds_gem$Sgrid
S_fft <- fit_ds_fft$Sgrid
m_gem <- surface_L2_cone(ds$S_true, S_gem, ds$time, ds$theta)
m_fft <- surface_L2_cone(ds$S_true, S_fft, ds$time, ds$theta)

summary_ds <- tibble(
  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)),
  Rel_L2_Surface= c(m_gem$rel_L2,          m_fft$rel_L2)
)
kable(summary_ds, digits = 4, caption = "Double-Slit: Summary metrics: True Kimesurface vs. KPT-GEM and KPT-FFT")

# Additional debugging
cat("Length check:\n")
cat("ds$theta:", length(ds$theta), "\n")
cat("fit_ds_fft$theta:", length(fit_ds_fft$theta), "\n")
cat("ds$phi_true row:", length(ds$phi_true[1,]), "\n")
cat("fit_ds_fft$varphi_theta row:", length(fit_ds_fft$varphi_theta[1,]), "\n")

# Check if they're the same
all.equal(length(ds$theta), length(fit_ds_fft$theta))
all.equal(length(ds$phi_true[1,]), length(fit_ds_fft$varphi_theta[1,]))

# Overlays at representative times
pick_t <- unique(round(seq(1, length(ds$time), length.out = 4)))
# plots <- lapply(pick_t, function(kk) {
#   kpt_plot_phase_overlay(
#     ds$theta,
#     ds$phi_true[kk, ],
#     fit_ds_fft$varphi_theta[kk, ],
#     t_lab = paste("t =", round(ds$time[kk], 3)),
#     title = "Double-Slit: Phase (True vs. KPT-FFT)"
#   )
# })
plots <- lapply(pick_t, function(kk) {
    kpt_plot_phase_overlay(
        ds$theta,
        ds$phi_true[kk, ],
        fit_ds_fft$varphi_theta[kk, ],
        fit_ds_gem$varphi_theta[kk, ],
        t_lab = paste("t =", round(ds$time[kk], 3)),
        title = "Double-Slit: Phases (True vs. KPT-Estimated)"
    )
})
plots[[1]]; plots[[2]]; plots[[3]]; plots[[4]]

### Plot heatmap of tru and KPT recoverd phases
# ds$phi_true, fit_ds_gem$varphi_theta, ds$theta, ds$time)
# pm_fft <- phase_metrics_df(ds$phi_true, fit_ds_fft$varphi_theta, ds$theta, ds$time

hm_phase_true <- plot_ly(x=ds$time, y=ds$theta, z=ds$phi_true, type="heatmap", name="DS: True Phase",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "True Phase"))
hm_phase_gem  <- plot_ly(x=ds$time, y=ds$theta, z=fit_ds_gem$varphi_theta, type="heatmap", name="DS: KPT-GEM Phase", 
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-GEM Phase")) 
hm_phase_fft  <- plot_ly(x=ds$time, y=ds$theta, z=fit_ds_fft$varphi_theta, type="heatmap",  name="DS: KPT-FFT Phase",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-FFT Phase"))
# Display one interactively (others available in object)
# p_fft

subplot(hm_phase_true, hm_phase_gem, hm_phase_fft, nrows = 3) %>% 
  layout(title="Double-Slit Experiment: Heatmaps of True, KPT-GEM & KPT-FFT Phases")

### FIG. 3. (Double-Slit Experiment) Phase metrics (KL, symmetric JSD in bits, Hellinger, TV) 
# summarizing the performance of the KPT-GEM to recover ϕt at each time. JSD’s bit scale is an 
# information-theoretic distance that is easy to compare across runs
pm_gem_plot <- plot_ly(pm_gem, x = ~time) %>%
  add_trace(y = ~KL_true_est, type = 'scatter', mode = 'lines', 
            name = 'KL(true||est)', line = list(color = '#1f77b4', width = 2)) %>%
  add_trace(y = ~KL_est_true, type = 'scatter', mode = 'lines', 
            name = 'KL(est||true)', line = list(color = '#ff7f0e', width = 2)) %>%
  add_trace(y = ~JSD_bits, type = 'scatter', mode = 'lines', 
            name = 'JSD (bits)', line = list(color = '#2ca02c', width = 2)) %>%
  add_trace(y = ~Hellinger, type = 'scatter', mode = 'lines', 
            name = 'Hellinger', line = list(color = '#d62728', width = 2)) %>%
  add_trace(y = ~TV, type = 'scatter', mode = 'lines', 
            name = 'Total Variation', line = list(color = '#9467bd', width = 2)) %>%
  layout(title = list(text = "KPT-GEM: Distribution Distance Metrics Over Time", 
                      x = 0.5, xanchor = 'center'),
         xaxis = list(title = "Time", gridcolor = '#e0e0e0'),
         yaxis = list(title = "Distance Metric Value", gridcolor = '#e0e0e0'),
         plot_bgcolor = 'white',
         paper_bgcolor = 'white',
         legend = list(x = 0.02, y = 0.98, 
                       bgcolor = 'rgba(255,255,255,0.9)',
                       bordercolor = '#CCCCCC',
                       borderwidth = 1))

pm_gem_plot
```

```{r ds-surfaces-3d}
# Create plots
p_true <- plot_surface_3d(ds$time, ds$theta, ds$S_true, "True Surface",
                          opacity=0.95, colorscale='Viridis', colorbar_title = "True Surface")
p_gem  <- plot_surface_3d(ds$time, ds$theta, S_gem, "KPT-GEM Surface", 
                          opacity=0.4, colorscale='Blues', colorbar_title = "KPT-GEM Surface")
p_fft  <- plot_surface_3d(ds$time, ds$theta, S_fft, "KPT-FFT Surface",
                          opacity=0.1, colorscale='Reds', colorbar_title = "KPT-FFT Surface")

subplot(p_true, p_gem, p_fft) %>% 
  layout(title = "Double-Slit: True Kimesurface vs. KPT-GEM vs. KPT-FFT")

# p_true <- plot_surface_3d(ds$time, ds$theta, ds$S_true, "DS: True Surface",
#                           opacity=0.9, colorscale='Viridis')
# p_gem  <- plot_surface_3d(ds$time, ds$theta, S_gem,     "DS: KPT-GEM Surface", 
#                           opacity=0.7, colorscale='red') 
# p_fft  <- plot_surface_3d(ds$time, ds$theta, S_fft,     "DS: KPT-FFT Surface",
#                           opacity=0.5, colorscale='blue')
# # Display one interactively (others available in object)
# # p_fft
# 
# subplot(p_true, p_gem, p_fft) %>% layout(title="Double-Slit: True Kimesurface vs. KPT-GEM vs. KPT-FFT")

#### Heatmaps
hm_p_true <- plot_ly(x=ds$time, y=ds$theta, z=ds$S_true, type="heatmap", name="DS: True Surface",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "True Surface"))
hm_p_gem  <- plot_ly(x=ds$time, y=ds$theta, z=S_gem, type="heatmap", name="DS: KPT-GEM Surface", 
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-GEM Surface")) 
hm_p_fft  <- plot_ly(x=ds$time, y=ds$theta, z=S_fft, type="heatmap",  name="DS: KPT-FFT Surface",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-FFT Surface"))
# Display one interactively (others available in object)
# p_fft

subplot(hm_p_true, hm_p_gem, hm_p_fft, nrows = 3) %>% 
  layout(title="Double-Slit Experiment: Heatmaps of True Kimesurface, KPT-GEM & KPT-FFT")

## 3 3D scenes with kime-surfaces
# Create individual plots with distinct scene names
# Create the three surfaces using your original approach
p_true <- plot_ly(x = ds$time, y = ds$theta, z = ds$S_true, type = "surface", 
                  scene = 'scene1', name = "True Surface",
                  colorscale = 'Viridis', opacity = 0.95) %>% 
  layout(title = "DS: True Surface")

p_gem <- plot_ly(x = ds$time, y = ds$theta, z = S_gem, type = "surface", 
                 scene = 'scene2', name = "KPT-GEM Surface",
                 colorscale = 'Blues', opacity = 0.4) %>% 
  layout(title = "DS: KPT-GEM Surface")

p_fft <- plot_ly(x = ds$time, y = ds$theta, z = S_fft, type = "surface",
                 scene = 'scene3', name = "KPT-FFT Surface",
                 colorscale = 'Reds', opacity = 0.1) %>% 
  layout(title = "DS: KPT-FFT Surface")

# Scene titles/annotations
annotations <- list(
  list(showarrow = FALSE, text = 'True Surface', 
       xref = 'paper', yref = 'paper', x = 0.15, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14)),
  list(showarrow = FALSE, text = 'KPT-GEM Surface', 
       xref = 'paper', yref = 'paper', x = 0.5, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14)),
  list(showarrow = FALSE, text = 'KPT-FFT Surface', 
       xref = 'paper', yref = 'paper', x = 0.85, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14))
)

# Combine the plots with explicit domain assignments
combined_plot <- subplot(p_true, p_gem, p_fft) %>%
  layout(
    title = "DS: True Kimesurface vs. KPT-GEM vs. KPT-FFT",
    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))
    ),
    annotations = annotations,
    margin = list(l = 0, r = 0, b = 0, t = 50)
  )

# Add the same camera synchronization JavaScript from your working example
combined_plot <- combined_plot %>% 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
```



# Simulation 2: fMRI ON vs. OFF

We generate an ON/OFF event design, produce an ON-OFF signal per $subject \times run$, and impose a 
phase process that is *more concentrated* during $ON$ (stimulus) blocks (latent synchrony).

```{r sim-fmri}
simulate_fmri_onoff <- function(P = 10, R = 20, TR = 2, duration = 512, # duration = 512 = 2 * 256
                                L_theta = 256, noise_sd = 0.20,
                                kappa_on = 5, kappa_off = 0.5, seed = 77) {
  set.seed(seed)
  K <- as.integer(duration / TR); t <- seq(0, by = TR, length.out = K)
  # Block design: 30s ON/OFF
  block_len <- 30; on_len <- block_len/(2*TR); design <- rep(c(0,1), each = on_len)
  design <- rep_len(design, K)
  # HRF-like shape (simple)
  hrf <- function(x) (x^5 * exp(-x)) / gamma(6) - 0.35 * (x^7 * exp(-x)) / gamma(8)
  h <- hrf(seq(0, 12, length.out = on_len))
  A <- stats::filter(rep_len(design, K), filter = h, method = "convolution", sides = 1)
  A <- as.numeric(A); A[is.na(A)] <- 0; A <- scales::rescale(A, to = c(0.2, 1.0))
  mu <- 0.05 * sin(2*pi*t/max(t))

  # Observations per subject×run: Y = (ON - OFF) proxy style
  I <- P * R
  ON  <- array(0, dim = c(P, R, K))
  OFF <- array(0, dim = c(P, R, K))
  for (p in 1:P) for (r in 1:R) {
    OFF[p, r, ] <- mu + rnorm(K, 0, noise_sd)
    ON[p, r, ]  <- mu + A + rnorm(K, 0, noise_sd)
  }
  Y <- sapply(1:K, function(k) as.numeric(ON[, , k] - OFF[, , k]))  # [I x K]

  # Construct S_true(t,θ) ~ baseline + A(t)*cos(θ) + 0.25 A(t) cos(2θ)
  # vs. S(t,\theta)=B(t)+A_1(t)\cos(\kappa t+\theta)+A_2(t)\cos(2(\kappa t+\theta)+\phi_2)

  theta <- theta_grid(L_theta)
  S_true <- sapply(1:K, function(k) mu[k] + A[k]*cos(theta) + 0.25*A[k]*cos(2*theta)) %>% t()

  # Phase law more concentrated in ON than OFF
  phi_true <- matrix(0, nrow = K, ncol = L_theta)
  for (k in 1:K) {
    kap <- if (design[k] == 1) kappa_on else kappa_off
    p <- (1/(2*pi*besselI(kap, 0))) * exp(kap * cos(theta - 0))  # μ=0
    phi_true[k, ] <- p / mean(p)
  }
  list(Y = Y, time = t, theta = theta, phi_true = phi_true, S_true = S_true)
}

fm <- simulate_fmri_onoff()
str(fm, 1)
```

## Run KPT-GEM and KPT-FFT on fMRI data

```{r run-fmri}
fit_fm_gem <- kpt_gem(fm$Y, time = fm$time, J = 12, L_theta = length(fm$theta),
                      lambda = 1e-2, p = 1, max_iter = 25, tol = 8e-4, verbose = TRUE)

fit_fm_fft <- kpt_fft(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,
                      max_iter = 25, tol = 8e-4, verbose = TRUE)
```

## fMRI: Quantitative & Qualitative Evaluation

```{r eval-fmri}
pm_gem <- phase_metrics_df(fm$phi_true, fit_fm_gem$varphi_theta, fm$theta, fm$time)
pm_fft <- phase_metrics_df(fm$phi_true, fit_fm_fft$varphi_theta, fm$theta, fm$time)

S_gem <- fit_fm_gem$Sgrid
S_fft <- fit_fm_fft$Sgrid
m_gem <- surface_L2_cone(fm$S_true, S_gem, fm$time, fm$theta)
m_fft <- surface_L2_cone(fm$S_true, S_fft, fm$time, fm$theta)

summary_fmri <- tibble(
  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)),
  Rel_L2_Surface= c(m_gem$rel_L2,          m_fft$rel_L2)
)
kable(summary_fmri, digits = 4, caption = "fMRI ON/OFF: Summary metrics: True Kimesurface vs. KPT-GEM vs. KPT-FFT")

# Phase overlays (sample ON and OFF times)
idx_on  <- which(diff(sign(diff(fm$time)))==0)[1]; idx_on <- ifelse(is.na(idx_on), 30, idx_on)
idx_off <- max(1, idx_on - 10)
# p_on  <- plot_phase_overlay(fm$theta, fm$phi_true[idx_on, ],  fit_fm_fft$varphi_theta[idx_on, ],
#                             t_lab = sprintf("ON t=%.0fs", fm$time[idx_on]),
#                             title = "fMRI: Phase (True vs. KPT-FFT)")
# p_off <- plot_phase_overlay(fm$theta, fm$phi_true[idx_off, ], fit_fm_fft$varphi_theta[idx_off, ],
#                             t_lab = sprintf("OFF t=%.0fs", fm$time[idx_off]),
#                             title = "fMRI: Phase (True vs. KPT-FFT)")
# p_on

p_on  <- kpt_plot_phase_overlay(fm$theta, fm$phi_true[idx_on, ], fit_fm_fft$varphi_theta[idx_on, ],
                      fit_ds_gem$varphi_theta[idx_on, ], t_lab = sprintf("ON t=%.0fs", fm$time[idx_on]),
                      title = "fMRI: Phase (True vs. KPT-FFT)")
p_off <- kpt_plot_phase_overlay(fm$theta, fm$phi_true[idx_off, ], fit_fm_fft$varphi_theta[idx_off, ],
                      fit_ds_gem$varphi_theta[idx_on, ], t_lab = sprintf("OFF t=%.0fs", fm$time[idx_off]),
                      title = "fMRI: Phase (True vs. KPT-FFT)")
p_on
p_off

# # Dimension-checked overlay (unique name to avoid conflicts)
# kpt_plot_phase_overlay <- function(theta, phi_true_row, phi_est_row, t_lab = "", title = "") {
#   n_theta <- length(theta)
#   stopifnot(length(phi_true_row) == n_theta, length(phi_est_row) == n_theta)
# 
#   df <- tibble::tibble(
#     theta   = rep(theta, 2L),                            # length = 2*n
#     density = c(phi_true_row, phi_est_row),              # length = 2*n
#     kind    = rep(c("True", "Estimated"), each = n_theta) # length = 2*n
#   )
# 
#   ggplot2::ggplot(df, ggplot2::aes(theta, density, color = kind)) +
#     ggplot2::geom_line(linewidth = 0.9) +
#     ggplot2::labs(title = title, subtitle = t_lab,
#                   x = expression(theta), y = expression(varphi(t,theta))) +
#     ggplot2::theme_minimal()
# }

# Additional debugging
cat("fMRI: Length check:\n")
cat("fm$theta:", length(fm$theta), "\n")
cat("fit_fm_fft$theta:", length(fit_fm_fft$theta), "\n")
cat("fm$phi_true row:", length(fm$phi_true[1,]), "\n")
cat("fit_fm_fft$varphi_theta row:", length(fit_fm_fft$varphi_theta[1,]), "\n")

# Check if they're the same
all.equal(length(fm$theta), length(fit_fm_fft$theta))
all.equal(length(fm$phi_true[1,]), length(fit_fm_fft$varphi_theta[1,]))

# Overlays at representative times
# pick_t <- unique(round(seq(1, length(fm$time), length.out = 4)))
# plots <- lapply(pick_t, function(kk) {
#   kpt_plot_phase_overlay(
#     fm$theta,
#     fm$phi_true[kk, ],
#     fit_fm_fft$varphi_theta[kk, ],
#     fit_fm_fft$varphi_theta[kk, ],
#     t_lab = paste("t =", round(fm$time[kk], 3)),
#     title = "fMRI: Phase (True vs. KPT-FFT)"
#   )
# })
# plots[[1]]; plots[[2]]; plots[[3]]; plots[[4]]
pick_t <- unique(round(seq(1, length(fm$time), length.out = 10)))
plots <- lapply(pick_t, function(kk) {
    kpt_plot_phase_overlay(
        ds$theta,
        ds$phi_true[kk, ],
        fit_ds_fft$varphi_theta[kk, ],
        fit_ds_gem$varphi_theta[kk, ],
        t_lab = paste("t =", round(ds$time[kk], 3)),
        title = "fMRI: Phases (True vs. KPT-Estimated)"
    )
})
plots[[1]]; plots[[2]]; plots[[3]]; plots[[4]]; plots[[10]]
```

```{r fmri-surfaces-3d}
p_true <- plot_surface_3d(fm$time, fm$theta, fm$S_true, "FMRI: True Surface",
                          opacity=0.95, colorscale='Viridis', colorbar_title = "True Surface")
p_gem  <- plot_surface_3d(fm$time, fm$theta, S_gem, "FMRI: KPT-GEM Surface", 
                          opacity=0.4, colorscale='Blues', colorbar_title = "KPT-GEM Surface")
p_fft  <- plot_surface_3d(fm$time, fm$theta, S_fft, "FMRI: KPT-FFT Surface",
                          opacity=0.1, colorscale='Reds', colorbar_title = "KPT-FFT Surface")
subplot(p_true, p_gem, p_fft) %>% 
  layout(title = "fMRI: True Kimesurface vs. KPT-GEM vs. KPT-FFT")

#### Heatmaps
hm_p_true <- plot_ly(x=fm$time, y=fm$theta, z=fm$S_true, type="heatmap", name="FMRI: True Surface",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "True Surface"))
hm_p_gem  <- plot_ly(x=fm$time, y=fm$theta, z=S_gem, type="heatmap", name="FMRI: KPT-GEM Surface", 
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-GEM Surface")) 
hm_p_fft  <- plot_ly(x=fm$time, y=fm$theta, z=S_fft, type="heatmap",  name="FMRI: KPT-FFT Surface",
                          opacity=0.9, colorscale='Viridis', colorbar = list(title = "KPT-FFT Surface"))
# Display one interactively (others available in object) # p_fft
subplot(hm_p_true, hm_p_gem, hm_p_fft, nrows = 3) %>% 
  layout(title="fMRI Experiment: Heatmaps of True Kimesurface, KPT-GEM & KPT-FFT")

## 3 3D scenes with kime-surfaces
# Create individual plots with distinct scene names
# Create the three surfaces using your original approach
p_true <- plot_ly(x = fm$time, y = fm$theta, z = fm$S_true, type = "surface", 
                  scene = 'scene1', name = "True Surface",
                  colorscale = 'Viridis', opacity = 0.95) %>% 
  layout(title = "FMRI: True Surface")

p_gem <- plot_ly(x = fm$time, y = fm$theta, z = S_gem, type = "surface", 
                 scene = 'scene2', name = "KPT-GEM Surface",
                 colorscale = 'Blues', opacity = 0.4) %>% 
  layout(title = "FMRI: KPT-GEM Surface")

p_fft <- plot_ly(x = fm$time, y = fm$theta, z = S_fft, type = "surface",
                 scene = 'scene3', name = "KPT-FFT Surface",
                 colorscale = 'Reds', opacity = 0.1) %>% 
  layout(title = "FMRI: KPT-FFT Surface")

# Scene titles/annotations
annotations <- list(
  list(showarrow = FALSE, text = 'True Surface', 
       xref = 'paper', yref = 'paper', x = 0.15, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14)),
  list(showarrow = FALSE, text = 'KPT-GEM Surface', 
       xref = 'paper', yref = 'paper', x = 0.5, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14)),
  list(showarrow = FALSE, text = 'KPT-FFT Surface', 
       xref = 'paper', yref = 'paper', x = 0.85, y = 0.95,
       xanchor = 'center', yanchor = 'bottom', font = list(size = 14))
)

# Combine the plots with explicit domain assignments
combined_plot <- subplot(p_true, p_gem, p_fft) %>%
  layout(
    title = "fMRI: True Kimesurface vs. KPT-GEM vs. KPT-FFT",
    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))
    ),
    annotations = annotations,
    margin = list(l = 0, r = 0, b = 0, t = 50)
  )

# Add the same camera synchronization JavaScript from your working example
combined_plot <- combined_plot %>% 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
```

# Consolidated Performance Tables

```{r tables}
bind_rows(
  summary_ds  %>% mutate(Experiment = "Double-Slit"),
  summary_fmri%>% mutate(Experiment = "fMRI ON/OFF")
) %>%
  relocate(Experiment) %>%
  arrange(Experiment, Algorithm) %>%
  kable(digits = 4, caption = "Consolidated metrics across experiments")
```

# Notes on Reproducibility & Alignment with Paper

 - *Normalization & geometry.* Phase integrals use $d\theta/(2\pi)$; row-mean-1 enforcement ensures densities integrate to $1$ on a uniform grid. All kime-surface comparisons use the *cone-weighted* $t\,dt\,d\theta$ $L^2$ norm.
 - *E-step statistics.* We use the unified *posterior Fourier expectations* $\xi_{i,k}(n)$ (row-wise softmax), eliminating notation drift and matching the paper’s algorithms.
 - *Spectral steps.* Both KPT-GEM and KPT-FFT algorithms implement the *Wiener–Sobolev* filter; KPT-FFT alternates $\Phi$- and $F$-updates, while KPT-GEM uses a stabilized surface proxy for robustness. The *reality constraint* $f_{-n}=\overline{f_n}$ is enforced in both algorithms.
 -*Anchoring.* A global phase-shift ambiguity is inherent; visualizations here keep the canonical grid (anchoring options can be added if needed).


<!--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-->