This TCIU Appendix 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.

1 Introduction

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

A deterministic kime-surface \(\mathcal{S}(t,\theta)\) over complex time \(\kappa = te^{i\theta}\), and
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.

2 Mathematical Framework

2.1 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).

2.2 KPT Algorithms

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

2.2.1 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:

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)\)
Mixed moments: \(m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)\)
Φ-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)\]
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)\]
Project: Ensure \(\varphi^{(\ell+1)}_{t_k} \geq 0\) and \(\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1\)

2.2.2 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.

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:

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)\)
Mixed moments: \(m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)\)
\(\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)\]
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)\]
Project: Ensure \(\varphi^{(\ell+1)}_{t_k} \geq 0\) and \(\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1\)

3 Implementation

3.1 Core 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")
}

3.2 Unified KPT Implementation

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

4 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.

4.1 Setup and 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)
}

4.2 Double-Slit Simulation

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)

4.3 Run KPT with Proper Alignment

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

## Alignment verification:

cat("Mean rotation angle (FFT):", mean(delta_fft), "rad\n")

## Mean rotation angle (FFT): -0.1223663 rad

cat("Mean rotation angle (GEM):", mean(delta_gem), "rad\n")

## Mean rotation angle (GEM): -0.1884956 rad

cat("Phase normalization check (should be ~1):\n")

## Phase normalization check (should be ~1):

cat("  FFT aligned:", mean(rowMeans(phi_fft_aligned)), "\n")

##   FFT aligned: 1

cat("  GEM aligned:", mean(rowMeans(phi_gem_aligned)), "\n")

##   GEM aligned: 1

4.4 Quantitative Evaluation (Post-Alignment)

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

Double-Slit: Performance Metrics (After Alignment)
Algorithm	Mean_JSD_bits	Mean_Hellinger	Mean_TV	Surface_RelL2
KPT-GEM	0.5839	0.7095	0.7382	1.3126
KPT-FFT	0.5129	0.6541	0.6877	1.3486

# Compare before/after alignment
cat("\n=== Alignment Impact ===\n")

## 
## === Alignment Impact ===

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

## FFT Mean JSD before alignment: 0.5074208 bits

cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n")

## FFT Mean JSD after alignment: 0.5129105 bits

cat("Improvement:", 
    mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n")

## Improvement: -0.005489682 bits

4.5 Visualization (Aligned Results)

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

4.6 Kime-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

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

## 
## === FINAL RESULTS (Double-Slit with Alignment) ===

print(summary_ds)

##   Algorithm Mean_JSD_bits Mean_Hellinger   Mean_TV Surface_RelL2
## 1   KPT-GEM     0.5839381      0.7094548 0.7382164      1.312560
## 2   KPT-FFT     0.5129105      0.6540909 0.6877467      1.348592

cat("\nInterpretation:\n")

## 
## Interpretation:

cat("- JSD near 0: Perfect recovery\n")

## - JSD near 0: Perfect recovery

cat("- Hellinger in [0,1]: Lower is better\n")

## - Hellinger in [0,1]: Lower is better

cat("- Total Variation in [0,1]: Lower is better\n")

## - Total Variation in [0,1]: Lower is better

cat("- Relative L2 < 1: Good surface recovery\n")

## - Relative L2 < 1: Good surface recovery

4.8 Key Observations

The alignment step is critical for proper evaluation because:

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.
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.
Impact: Without alignment, the metrics can be artificially inflated. The aligned results show the true recovery performance.
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.

5 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).

5.1 fMRI Simulation with Phase-Dependent BOLD Signal

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

## fMRI simulation complete:

cat("  Subjects × Runs:", fm$params$P, "×", fm$params$R, "\n")

##   Subjects × Runs: 10 × 20

cat("  Time points:", ncol(fm$Y), "\n")

##   Time points: 150

cat("  Phase resolution:", length(fm$theta), "\n")

##   Phase resolution: 256

5.2 Run KPT Analysis on fMRI Data

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

## 
## === Alignment Results ===

cat("Mean rotation FFT:", mean(delta_fft), "rad\n")

## Mean rotation FFT: 0.1679142 rad

cat("Mean rotation GEM:", mean(delta_gem), "rad\n")

## Mean rotation GEM: -0.02094395 rad

cat("Phase normalization check:\n")

## Phase normalization check:

cat("  FFT:", mean(rowMeans(phi_fft_aligned)), "\n")

##   FFT: 1

cat("  GEM:", mean(rowMeans(phi_gem_aligned)), "\n")

##   GEM: 1

5.3 Quantitative Evaluation

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

fMRI ON/OFF: Performance Metrics (After Alignment)
Algorithm	Mean_JSD_bits	Mean_Hellinger	Mean_TV	Surface_RelL2
KPT-GEM	0.4445	0.5942	0.6491	1.5013
KPT-FFT	0.3794	0.5534	0.5588	1.4574

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

## 
## === Alignment Impact ===

cat("FFT Mean JSD before alignment:", mean(pm_fft_unaligned$JSD_bits), "bits\n")

## FFT Mean JSD before alignment: 0.3847786 bits

cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n")

## FFT Mean JSD after alignment: 0.3794252 bits

cat("Improvement:", mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n")

## Improvement: 0.005353388 bits

5.4 ON vs OFF Block 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")

fMRI: Performance by Block Type
Algorithm	Block	Mean_JSD	Mean_TV
KPT-GEM	ON	0.5422	0.7290
KPT-FFT	ON	0.5147	0.6812
KPT-GEM	OFF	0.3467	0.5692
KPT-FFT	OFF	0.2442	0.4363

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

## 
## === Selected Time Points ===

cat("Peak ON at t =", fm$time[idx_on_peak], "s\n")

## Peak ON at t = 42 s

cat("Baseline OFF at t =", fm$time[idx_off_baseline], "s\n")

## Baseline OFF at t = 72 s

5.5 Visualization

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

5.6 Phase Distribution 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")

5.7 Kime-Surface Comparison

# 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

5.8 Time-Resolved 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")

## 
## === Surface Correlation ===

cat("GEM: mean =", mean(cors_gem), ", sd =", sd(cors_gem), "\n")

## GEM: mean = 0.7364706 , sd = 0.6553402

cat("FFT: mean =", mean(cors_fft), ", sd =", sd(cors_fft), "\n")

## FFT: mean = 0.6433351 , sd = 0.41099

5.9 Final Summary

# Consolidated results
cat("\n=== FINAL RESULTS (fMRI with Alignment) ===\n")

## 
## === FINAL RESULTS (fMRI with Alignment) ===

print(summary_fmri)

##   Algorithm Mean_JSD_bits Mean_Hellinger   Mean_TV Surface_RelL2
## 1   KPT-GEM     0.4444587      0.5942159 0.6490989      1.501273
## 2   KPT-FFT     0.3794252      0.5533982 0.5587536      1.457355

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

fMRI: Performance by Block Type
Algorithm	Mean_JSD_bits	Mean_Hellinger	Mean_TV	Surface_RelL2
KPT-GEM	0.4445	0.5942	0.6491	1.5013
KPT-FFT	0.3794	0.5534	0.5588	1.4574

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

## 
## === Performance by Block Type ===

print(summary_by_block)

##   Algorithm Block  Mean_JSD   Mean_TV
## 1   KPT-GEM    ON 0.5422401 0.7290193
## 2   KPT-FFT    ON 0.5146552 0.6812139
## 3   KPT-GEM   OFF 0.3466774 0.5691786
## 4   KPT-FFT   OFF 0.2441952 0.4362933

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

fMRI: Performance by Block Type
Algorithm	Block	Mean_JSD	Mean_TV
KPT-GEM	ON	0.5422	0.7290
KPT-FFT	ON	0.5147	0.6812
KPT-GEM	OFF	0.3467	0.5692
KPT-FFT	OFF	0.2442	0.4363

cat("\n=== Key Findings ===\n")

## 
## === Key Findings ===

cat("1. Phase concentration higher during ON blocks (kappa =", fm$params$kappa_on, ")\n")

## 1. Phase concentration higher during ON blocks (kappa = 6 )

cat("2. Phase concentration lower during OFF blocks (kappa =", fm$params$kappa_off, ")\n")

## 2. Phase concentration lower during OFF blocks (kappa = 1 )

cat("3. Both algorithms recover this structure after alignment\n")

## 3. Both algorithms recover this structure after alignment

cat("4. FFT (with surface update) performs slightly better overall\n")

## 4. FFT (with surface update) performs slightly better overall

# Performance interpretation
cat("\n=== Metric Interpretation ===\n")

## 
## === Metric Interpretation ===

cat("- JSD = 0: Perfect recovery (identical distributions)\n")

## - JSD = 0: Perfect recovery (identical distributions)

cat("- JSD < 0.1 bits: Excellent recovery\n")

## - JSD < 0.1 bits: Excellent recovery

cat("- JSD < 0.5 bits: Good recovery\n")

## - JSD < 0.5 bits: Good recovery

cat("- Relative L2 < 1: Better than baseline\n")

## - Relative L2 < 1: Better than baseline

Time Complexity, Inferential uncertainty and Spacekime Analytics

Appendix 03: Kime-Phase Tomography (KPT V.5: KPT-GEM & KPT-FFT): fMRI and Double-Slit Validation

SOCR Team

10/13/2025