--- title: "Time Complexity, Inferential uncertainty and Spacekime Analytics" subtitle: "Appendix 03: Kime-Phase Tomography (KPT V.5: KPT-GEM & KPT-FFT): fMRI and Double-Slit Validation" author: "SOCR Team" date: "`r format(Sys.time(),'%m/%d/%Y')`" output: html_document: theme: spacelab highlight: tango includes: before_body: TCIU_header.html toc: true number_sections: true toc_depth: 3 toc_float: collapsed: false smooth_scroll: true code_folding: hide pandoc_args: ["+RTS", "-K1024m", "-RTS"] word_document: toc: true toc_depth: '3' header-includes: - \usepackage{amsmath} - \usepackage{amssymb} - \usepackage{physics} - \usepackage{mathrsfs} - \usepackage{bm} --- This [TCIU Appendix](https://tciu.predictive.space/) offers a self-contained implementation of the kime phase tompgraphy algorithms, *KPT-GEM* (expectation maximization) and *KPT-FFT (FFT alternating spectral)* for estimating the phase distribution and reconstructing the kime-surface representations from observed repeated measurement data from longitudinal processes. The KPT algorihtm is validated using two simulations; *double-slit physics experiment* and *fMRI ON vs. OFF neuroscience experiment*. We report quyantitative performance metrics (summary tables) and provide qualitative results as surface/heatmap overlays contrasting for the *true* vs. *estimated* phase distributions and kime-surfaces. The key steps of the KPT algorithm include E-step, Wiener-Sobolev filtering, simplex projection, ans reality constraints, subject to normalization, operator conventions. ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) suppressPackageStartupMessages({ library(ggplot2) library(dplyr) library(tidyr) library(Matrix) library(plotly) library(knitr) library(purrr) }) set.seed(2025) ``` # Introduction This notebook presents a unified implementation of *Kime-Phase Tomography (KPT)*, a framework for analyzing repeated measurement longitudinal data by decomposing observations into 1. A deterministic *kime-surface* $\mathcal{S}(t,\theta)$ over complex time $\kappa = te^{i\theta}$, and 2. A stochastic *phase distribution* $\varphi_t(\theta)$ capturing intrinsic variability. The forward model for observations is $$y_{j,k} = \mathcal{S}(t_k, \Theta_j(t_k)) + \varepsilon_{j,k},$$ where $\Theta_j(t_k) \sim \varphi_{t_k}(\cdot)$ is a latent phase and $\varepsilon_{j,k} \sim \mathcal{N}(0,\sigma^2)$ is measurement noise. # Mathematical Framework ## KPT Inverse Problem Given repeated measurements $\{y_{j,k}\}_{j=1}^N$ at times $\{t_k\}_{k=1}^K$, KPT solves the inverse problem: *Recover:* - Phase distributions $\{\varphi_{t_k}\}$ - Kime-surface $\mathcal{S}(t,\theta)$ *Key insight:* This deconvolution becomes tractable in the Fourier domain where: $$M(z,t) = F(z,t) \cdot \Phi(z,t), \quad |z|=1$$ where $F(z,t) = \sum_n f_n(t)z^n$ (surface transform), $\Phi(z,t) = \sum_n \hat{\varphi}_t(n)z^n$ (phase characteristic function), and $M(z,t) = \sum_n m_n(t)z^n$ (observable mixed moments). ## KPT Algorithms Both algorithms implement the same core spectral deconvolution with different update strategies: ### Algorithm 1: KPT-GEM (Generalized EM) *Inputs:* Data $\{y_{j,k}\}$, regularization $(\lambda_\Phi, p_\Phi)$, $(\lambda_F, p_F)$ *Initialize:* $\varphi^{(0)}_{t_k} \equiv \frac{1}{2\pi}$, $f^{(0)}_n(t_k)$ from data *Repeat until convergence:* 1. *E-step:* Compute posterior expectations $$\xi_{j,k}^{(\ell)}(n) = \mathbb{E}[e^{-in\Theta_j} | y_{j,k}] = \int_0^{2\pi} e^{-in\theta} w_{j,k}^{(\ell)}(\theta) \frac{d\theta}{2\pi}$$ where $w_{j,k}^{(\ell)}(\theta) \propto \varphi^{(\ell-1)}_{t_k}(\theta) \exp\left(-\frac{(y_{j,k} - \mathcal{S}^{(\ell-1)}(t_k,\theta))^2}{2\sigma^2}\right)$ 2. *Mixed moments:* $m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)$ 3. *Φ-step (phase update):* Apply Wiener filter $$\hat{\Phi}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{F}^{(\ell)}(\omega,t_k)}} {|\hat{F}^{(\ell)}(\omega,t_k)|^2 + \lambda_\Phi \Lambda_{p_\Phi}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$ 4. *Optional F-step (surface update):* $$\hat{F}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{\Phi}^{(\ell+1)}(\omega,t_k)}} {|\hat{\Phi}^{(\ell+1)}(\omega,t_k)|^2 + \lambda_F \Lambda_{p_F}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$ 5. *Project:* Ensure $\varphi^{(\ell+1)}_{t_k} \geq 0$ and $\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1$ ### Algorithm 2: KPT-FFT (Fully Alternating) Same as KPT-GEM but *always* alternates between Φ-step and F-step updates at each iteration. The Sobolev penalty $\Lambda_p(\omega) = (2\sin(\omega/2))^{2p}$ controls smoothness. *Inputs:* Data $\{y_{j,k}\}$, regularization $(\lambda_\Phi, p_\Phi)$, $(\lambda_F, p_F)$ *Initialize:* $\varphi^{(0)}_{t_k} \equiv \frac{1}{2\pi}$, $f^{(0)}_n(t_k)$ from data *Repeat until convergence:* 1. *E-step:* Compute posterior expectations $$\xi_{j,k}^{(\ell)}(n) = \mathbb{E}[e^{-in\Theta_j} | y_{j,k}] = \int_0^{2\pi} e^{-in\theta} w_{j,k}^{(\ell)}(\theta) \frac{d\theta}{2\pi}$$ where $w_{j,k}^{(\ell)}(\theta) \propto \varphi^{(\ell-1)}_{t_k}(\theta) \exp\left(-\frac{(y_{j,k} - \mathcal{S}^{(\ell-1)}(t_k,\theta))^2}{2\sigma^2}\right)$ 2. *Mixed moments:* $m_n^{(\ell)}(t_k) = \frac{1}{N}\sum_{j=1}^N y_{j,k} \xi_{j,k}^{(\ell)}(n)$ 3. $\Phi$*-step (phase update):* Apply Wiener filter $$\hat{\Phi}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{F}^{(\ell)}(\omega,t_k)}} {|\hat{F}^{(\ell)}(\omega,t_k)|^2 + \lambda_\Phi \Lambda_{p_\Phi}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$ 4. *F-step (surface update):* $$\hat{F}^{(\ell+1)}(\omega,t_k) = \frac{\overline{\hat{\Phi}^{(\ell+1)}(\omega,t_k)}} {|\hat{\Phi}^{(\ell+1)}(\omega,t_k)|^2 + \lambda_F \Lambda_{p_F}(\omega)} \hat{M}^{(\ell)}(\omega,t_k)$$ 5. *Project:* Ensure $\varphi^{(\ell+1)}_{t_k} \geq 0$ and $\int \varphi^{(\ell+1)}_{t_k} d\theta/(2\pi) = 1$ # Implementation ## Core Utilities ```{r utilities} # ---- Grid generation ---- theta_grid <- function(L) seq(0, 2*pi, length.out = L + 1L)[1L:L] make_nvec <- function(J) seq.int(-J, J, by = 1L) # ---- Probability simplex projection ---- row_project_simplex <- function(P, eps = 1e-12) { P <- pmax(P, 0) m <- rowMeans(P) m[m < eps] <- 1 sweep(P, 1L, m, `/`) } # ---- Stable softmax for posterior weights ---- softmax_rows <- function(Lmat) { m <- apply(Lmat, 1L, max) W <- exp(Lmat - m) W / pmax(rowSums(W), .Machine$double.eps) } # ---- FFT coefficient mapping (n=-J..J to DFT order) ---- coeffs_to_dft_order <- function(x_n, nvec) { N <- length(nvec) x_m <- complex(length = N) for (i in seq_along(nvec)) { n <- nvec[i] m <- (n + N) %% N x_m[m + 1L] <- x_n[i] } x_m } dft_order_to_coeffs <- function(x_m, nvec) { N <- length(nvec) x_n <- complex(length = N) for (i in seq_along(nvec)) { n <- nvec[i] m <- (n + N) %% N x_n[i] <- x_m[m + 1L] } x_n } # ---- FFT-based evaluation on unit circle ---- coeffs_to_Xomega_fft <- function(x_n, nvec) { x_m <- coeffs_to_dft_order(x_n, nvec) fft(x_m, inverse = FALSE) } Xomega_to_coeffs_fft <- function(X_w, nvec) { N <- length(nvec) x_m <- fft(X_w, inverse = TRUE) / N dft_order_to_coeffs(x_m, nvec) } # ---- Surface synthesis from Fourier coefficients ---- surface_from_coeffs <- function(f_n_row, nvec, theta) { E <- exp(1i * outer(nvec, theta)) as.vector(Re(drop(f_n_row %*% E))) } # ---- Reality constraint: f_{-n} = conj(f_n) ---- enforce_reality_coeffs <- function(f_n, nvec) { N <- length(nvec) J <- (N - 1L) / 2L # n=0 must be real f_n[nvec == 0] <- Re(f_n[nvec == 0]) # Conjugate symmetry for (j in 1:J) { pos_idx <- which(nvec == j) neg_idx <- which(nvec == -j) if (length(pos_idx) > 0 && length(neg_idx) > 0) { avg <- 0.5 * (f_n[pos_idx] + Conj(f_n[neg_idx])) f_n[pos_idx] <- avg f_n[neg_idx] <- Conj(avg) } } f_n } # ---- Phase distribution metrics ---- phase_metrics_df <- function(phi_true, phi_est, theta, time, base2 = TRUE) { # Normalize P <- pmax(phi_true, 1e-12) Q <- pmax(phi_est, 1e-12) P <- sweep(P, 1, rowMeans(P), `/`) Q <- sweep(Q, 1, rowMeans(Q), `/`) K <- nrow(P) logb <- if (base2) function(x) log2(pmax(x, 1e-12)) else log metrics <- data.frame( time = time, KL_true_est = numeric(K), KL_est_true = numeric(K), JSD_bits = numeric(K), Hellinger = numeric(K), TV = numeric(K) ) for (k in 1:K) { # KL divergences metrics$KL_true_est[k] <- mean(P[k,] * logb(P[k,]/Q[k,])) metrics$KL_est_true[k] <- mean(Q[k,] * logb(Q[k,]/P[k,])) # Jensen-Shannon divergence M <- 0.5 * (P[k,] + Q[k,]) metrics$JSD_bits[k] <- 0.5 * mean(P[k,] * logb(P[k,]/M)) + 0.5 * mean(Q[k,] * logb(Q[k,]/M)) # Hellinger distance metrics$Hellinger[k] <- sqrt(0.5 * mean((sqrt(P[k,]) - sqrt(Q[k,]))^2)) # Total variation metrics$TV[k] <- 0.5 * mean(abs(P[k,] - Q[k,])) } metrics } # ---- Surface L2 error (cone-weighted) ---- surface_L2_cone <- function(S_true, S_est, time, theta) { K <- nrow(S_true) L <- ncol(S_true) dtheta <- 2*pi / L # Time weights w_t <- time if (length(w_t) == 0 || K == 1) w_t <- rep(1, K) # Compute L2 norm diff2 <- (S_true - S_est)^2 per_t_L2_sq <- rowSums(diff2) * dtheta * w_t total_L2 <- sqrt(sum(per_t_L2_sq)) # Relative error true_norm_sq <- rowSums(S_true^2) * dtheta * w_t rel_L2 <- total_L2 / max(sqrt(sum(true_norm_sq)), 1e-12) list(total_L2 = total_L2, rel_L2 = rel_L2) } # ---- Visualization helpers ---- plot_phase_overlay <- function(theta, phi_true_row, phi_est_row1, phi_est_row2 = NULL, t_lab = "", title = "") { if (is.null(phi_est_row2)) { df <- data.frame( theta = rep(theta, 2), density = c(phi_true_row, phi_est_row1), Method = rep(c("True", "Estimated"), each = length(theta)) ) } else { df <- data.frame( theta = rep(theta, 3), density = c(phi_true_row, phi_est_row1, phi_est_row2), Method = rep(c("True", "KPT-GEM", "KPT-FFT"), each = length(theta)) ) } ggplot(df, aes(theta, density, color = Method)) + geom_line(linewidth = 1.2) + labs(title = title, subtitle = t_lab, x = expression(theta), y = expression(varphi(theta))) + theme_minimal() + theme(legend.position = "bottom") } ``` ## Unified KPT Implementation ```{r kpt-unified} # ---- Initialize surface from data ---- initialize_surface <- function(Y, time, theta, method = "empirical") { K <- ncol(Y) L <- length(theta) if (method == "empirical") { # Use mean and variance structure mu_k <- colMeans(Y) A_k <- apply(Y, 2, sd) # Smooth estimates if (K > 4) { mu_s <- smooth.spline(time, mu_k, spar = 0.6)$y A_s <- smooth.spline(time, A_k, spar = 0.6)$y } else { mu_s <- mu_k A_s <- A_k } # Build surface: baseline + amplitude * harmonics S <- matrix(0, K, L) for (k in 1:K) { S[k,] <- mu_s[k] + A_s[k] * cos(theta) + 0.25 * A_s[k] * cos(2*theta) } } else { # Simple constant initialization S <- matrix(mean(Y), K, L) } S } # ---- Main unified KPT function ---- kpt_unified <- function(Y, time = NULL, J = 12, L_theta = 256, lambda_phi = 1e-2, p_phi = 1, lambda_F = 1e-2, p_F = 1, algorithm = c("GEM", "FFT"), update_surface = TRUE, max_iter = 30, tol = 1e-3, sigma = NULL, verbose = TRUE) { algorithm <- match.arg(algorithm) I <- nrow(Y) # Number of subjects/runs K <- ncol(Y) # Number of time points # Default time grid if (is.null(time)) time <- seq_len(K) # Setup grids theta <- theta_grid(L_theta) nvec <- make_nvec(J) Nn <- length(nvec) omega <- 2*pi * (0:(Nn-1)) / Nn # Estimate noise if not provided if (is.null(sigma)) { sigma <- 0.5 * mad(as.vector(Y), center = 0) + 1e-6 } # Initialize phase distribution (uniform) phi_grid <- matrix(1.0, nrow = K, ncol = L_theta) # Initialize surface Sgrid <- initialize_surface(Y, time, theta) # Compute initial Fourier coefficients f_n <- matrix(0+0i, nrow = K, ncol = Nn) for (k in 1:K) { Eneg <- exp(-1i * outer(theta, nvec)) f_n[k,] <- as.vector((t(Sgrid[k,]) %*% Eneg) / L_theta) f_n[k,] <- enforce_reality_coeffs(f_n[k,], nvec) } # Sobolev penalties Lambda_phi <- (2 * sin(omega/2))^(2*p_phi) Lambda_F <- (2 * sin(omega/2))^(2*p_F) # Storage for convergence obj <- numeric(max_iter) # Main iteration loop for (it in 1:max_iter) { # ==== E-step: Compute posterior expectations ==== m_hat <- matrix(0+0i, nrow = K, ncol = Nn) for (k in 1:K) { # Prior at time k prior_k <- phi_grid[k,] # Likelihood: p(y|theta) Srow <- matrix(Sgrid[k,], nrow = I, ncol = L_theta, byrow = TRUE) LogL <- -(Y[,k] - Srow)^2 / (2 * sigma^2) # Posterior weights (properly normalized) LogW <- sweep(LogL, 2L, log(pmax(prior_k, 1e-12)), `+`) W <- softmax_rows(LogW) # [I x L_theta] # Compute expectations E[e^{-in*Theta} | y] Eneg <- exp(-1i * outer(theta, nvec)) xi <- W %*% Eneg # [I x Nn] # Mixed moments (corrected) m_hat[k,] <- colSums(Y[,k] * xi) / I } # ==== Φ-step: Update phase distribution ==== phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn) for (k in 1:K) { # FFT to evaluate on unit circle F_w <- coeffs_to_Xomega_fft(f_n[k,], nvec) M_w <- coeffs_to_Xomega_fft(m_hat[k,], nvec) # Wiener filter denom <- Mod(F_w)^2 + lambda_phi * Lambda_phi + 1e-10 Phi_w <- Conj(F_w) * M_w / denom # IFFT back to coefficients phi_hat_n[k,] <- Xomega_to_coeffs_fft(Phi_w, nvec) } # Synthesize and project to simplex phi_new <- matrix(0, nrow = K, ncol = L_theta) Epos <- exp(1i * outer(nvec, theta)) for (k in 1:K) { phi_row <- Re(drop(t(Epos) %*% phi_hat_n[k,])) phi_new[k,] <- phi_row } phi_new <- row_project_simplex(phi_new) # ==== F-step: Update surface (if requested) ==== if (update_surface && (algorithm == "FFT" || it > 5)) { f_n_new <- matrix(0+0i, nrow = K, ncol = Nn) for (k in 1:K) { # Use updated phase Phi_w <- coeffs_to_Xomega_fft(phi_hat_n[k,], nvec) M_w <- coeffs_to_Xomega_fft(m_hat[k,], nvec) # Wiener filter for surface denom <- Mod(Phi_w)^2 + lambda_F * Lambda_F + 1e-10 F_w <- Conj(Phi_w) * M_w / denom # Reality constraint f_n_new[k,] <- enforce_reality_coeffs(Xomega_to_coeffs_fft(F_w, nvec), nvec) } # Synthesize new surface S_new <- matrix(0, nrow = K, ncol = L_theta) for (k in 1:K) { S_new[k,] <- surface_from_coeffs(f_n_new[k,], nvec, theta) } # Update f_n <- f_n_new Sgrid <- S_new } # ==== Check convergence ==== diff_phi <- sqrt(mean((phi_new - phi_grid)^2)) diff_f <- if (update_surface && (algorithm == "FFT" || it > 5)) { sqrt(mean(Mod(f_n_new - f_n)^2)) } else { 0 } if (verbose) { if (algorithm == "GEM") { message(sprintf("[%s] iter %02d: Δφ=%.3e", algorithm, it, diff_phi)) } else { message(sprintf("[%s] iter %02d: Δφ=%.3e, Δf=%.3e", algorithm, it, diff_phi, diff_f)) } } # Update state phi_grid <- phi_new obj[it] <- mean(Mod(m_hat)^2) # Convergence check if (algorithm == "GEM" && diff_phi < tol) break if (algorithm == "FFT" && max(diff_phi, diff_f) < tol) break } # Return results list( varphi_theta = phi_grid, Sgrid = Sgrid, f_n = f_n, theta = theta, time = time, convergence = list( iterations = it, objective = obj[1:it] ), parameters = list( algorithm = algorithm, lambda_phi = lambda_phi, p_phi = p_phi, lambda_F = lambda_F, p_F = p_F, sigma = sigma, J = J ) ) } ``` # Double-Slit Experiment (with Phase Alignment) Since the phase distributions are only identified up to a *global rotation*, we need to align them prior to comparison, visualizaiton, contextual interpretation, or spacekime analytic prediction. ## Setup and Utilities ```{r setup-utilities} # Core utilities for phase alignment anchor_phase_by_m1 <- function(varphi_mat, theta) { L <- length(theta) dtheta <- 2*pi / L # Compute first trigonometric moment per row m1 <- apply(varphi_mat, 1L, function(row) { sum(exp(1i * theta) * row) * dtheta / (2*pi) }) # Extract rotation angles delta <- Arg(m1) # Rotate each row by -delta using FFT rotate_row <- function(row, delta_k) { F <- fft(row) n <- c(0:(L/2), -(L/2-1):-1) row_rot <- Re(fft(F * exp(-1i * n * delta_k), inverse = TRUE) / L) pmax(row_rot, 0) # Ensure non-negative } # Apply rotation out <- mapply(function(idx, d) rotate_row(varphi_mat[idx,], d), idx = seq_len(nrow(varphi_mat)), d = delta, SIMPLIFY = FALSE) out <- do.call(rbind, out) # Renormalize out / rowMeans(out) } rotate_surface_by_delta <- function(S_mat, theta, delta_vec) { stopifnot(nrow(S_mat) == length(delta_vec)) L <- length(theta) n <- c(0:(L/2), -(L/2-1):-1) rotate_row <- function(row, delta_k) { F <- fft(row) Re(fft(F * exp(-1i * n * delta_k), inverse = TRUE) / L) } out <- mapply(function(idx, d) rotate_row(S_mat[idx,], d), idx = seq_len(nrow(S_mat)), d = delta_vec, SIMPLIFY = FALSE) do.call(rbind, out) } ``` ## Double-Slit Simulation ```{r simulate-ds} simulate_double_slit_aligned <- function(N = 100, K = 50, L_theta = 256, kappa_phase = 4, noise_sd = 0.15, kappa_wave = 12, seed = 2025) { set.seed(seed) # Grids theta <- theta_grid(L_theta) time <- seq(0, 1, length.out = K) # Build true kime-surface: quantum interference pattern S_true <- matrix(0, K, L_theta) for (k in 1:K) { t <- time[k] # Gaussian envelope for wave packet envelope <- 0.8 * exp(-20*(t - 0.5)^2) # Base drift baseline <- 0.2 * sin(2*pi*t) # Interference pattern with two harmonics S_true[k,] <- baseline + envelope * cos(kappa_wave*t + theta) + 0.3 * envelope * cos(2*(kappa_wave*t + theta) + pi/6) } # Time-varying phase distribution (mixture of von Mises) # This represents the quantum state distribution phi_true <- matrix(0, K, L_theta) for (k in 1:K) { t <- time[k] # Two modes representing different quantum paths mu1 <- (2*pi * (0.2 + 0.6*t)) %% (2*pi) mu2 <- (2*pi * (0.9 - 0.6*t)) %% (2*pi) w <- 0.5 + 0.4*sin(2*pi*t) # Time-varying mixture weight # von Mises mixture p1 <- exp(kappa_phase * cos(theta - mu1)) / (2*pi * besselI(kappa_phase, 0)) p2 <- exp(kappa_phase * cos(theta - mu2)) / (2*pi * besselI(kappa_phase, 0)) p <- w * p1 + (1-w) * p2 phi_true[k,] <- p / mean(p) # Normalize to mean=1 } # Generate observations: Y = S(t, Theta) + noise Y <- matrix(0, N, K) for (k in 1:K) { # Sample phases from the true distribution probs <- pmax(phi_true[k,], 1e-12) idx <- sample.int(L_theta, N, replace = TRUE, prob = probs) # Generate measurements Y[,k] <- S_true[k, idx] + rnorm(N, 0, noise_sd) } list( Y = Y, time = time, theta = theta, phi_true = phi_true, S_true = S_true, params = list(N = N, K = K, L_theta = L_theta, kappa_phase = kappa_phase, noise_sd = noise_sd) ) } # Generate data ds <- simulate_double_slit_aligned(N = 100, K = 50) ``` ## Run KPT with Proper Alignment ```{r run-kpt-aligned} # Run KPT-GEM fit_ds_gem <- kpt_unified( ds$Y, time = ds$time, J = 16, L_theta = length(ds$theta), lambda_phi = 5e-3, p_phi = 1, lambda_F = 5e-3, p_F = 1, algorithm = "GEM", update_surface = FALSE, max_iter = 30, tol = 1e-3, verbose = TRUE ) # Run KPT-FFT fit_ds_fft <- kpt_unified( ds$Y, time = ds$time, J = 16, L_theta = length(ds$theta), lambda_phi = 5e-3, p_phi = 1, lambda_F = 5e-3, p_F = 1, algorithm = "FFT", update_surface = TRUE, max_iter = 30, tol = 1e-3, verbose = TRUE ) # CRITICAL: Compute alignment angles from FFT estimates m1_fft <- apply(fit_ds_fft$varphi_theta, 1L, function(row) { mean(exp(1i * ds$theta) * row) }) delta_fft <- Arg(m1_fft) # Similarly for GEM m1_gem <- apply(fit_ds_gem$varphi_theta, 1L, function(row) { mean(exp(1i * ds$theta) * row) }) delta_gem <- Arg(m1_gem) # Align phase distributions phi_fft_aligned <- anchor_phase_by_m1(fit_ds_fft$varphi_theta, ds$theta) phi_gem_aligned <- anchor_phase_by_m1(fit_ds_gem$varphi_theta, ds$theta) # Align kime-surfaces using the same rotation S_fft_aligned <- rotate_surface_by_delta(fit_ds_fft$Sgrid, ds$theta, delta_fft) S_gem_aligned <- rotate_surface_by_delta(fit_ds_gem$Sgrid, ds$theta, delta_gem) # Verification checks cat("Alignment verification:\n") cat("Mean rotation angle (FFT):", mean(delta_fft), "rad\n") cat("Mean rotation angle (GEM):", mean(delta_gem), "rad\n") cat("Phase normalization check (should be ~1):\n") cat(" FFT aligned:", mean(rowMeans(phi_fft_aligned)), "\n") cat(" GEM aligned:", mean(rowMeans(phi_gem_aligned)), "\n") ``` ## Quantitative Evaluation (Post-Alignment) ```{r evaluate-aligned} # Compute metrics on ALIGNED distributions pm_gem <- phase_metrics_df(ds$phi_true, phi_gem_aligned, ds$theta, ds$time) pm_fft <- phase_metrics_df(ds$phi_true, phi_fft_aligned, ds$theta, ds$time) # Surface errors on ALIGNED surfaces m_gem <- surface_L2_cone(ds$S_true, S_gem_aligned, ds$time, ds$theta) m_fft <- surface_L2_cone(ds$S_true, S_fft_aligned, ds$time, ds$theta) # Summary table summary_ds <- data.frame( Algorithm = c("KPT-GEM", "KPT-FFT"), Mean_JSD_bits = c(mean(pm_gem$JSD_bits), mean(pm_fft$JSD_bits)), Mean_Hellinger = c(mean(pm_gem$Hellinger), mean(pm_fft$Hellinger)), Mean_TV = c(mean(pm_gem$TV), mean(pm_fft$TV)), Surface_RelL2 = c(m_gem$rel_L2, m_fft$rel_L2) ) kable(summary_ds, digits = 4, caption = "Double-Slit: Performance Metrics (After Alignment)") # Compare before/after alignment cat("\n=== Alignment Impact ===\n") pm_fft_unaligned <- phase_metrics_df(ds$phi_true, fit_ds_fft$varphi_theta, ds$theta, ds$time) cat("FFT Mean JSD before alignment:", mean(pm_fft_unaligned$JSD_bits), "bits\n") cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n") cat("Improvement:", mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n") ``` ## Visualization (Aligned Results) ```{r visualize-aligned} # Select representative time points time_indices <- c(1, 13, 25, 37, 50) # Plot aligned phase distributions for (k in time_indices) { p <- plot_phase_overlay( ds$theta, ds$phi_true[k,], phi_gem_aligned[k,], phi_fft_aligned[k,], t_lab = sprintf("t = %.2f", ds$time[k]), title = "Double-Slit: Aligned Phase Distributions" ) print(p) } # Create phase heatmaps hm_true <- plot_ly( x = ds$time, y = ds$theta, z = t(ds$phi_true), type = "heatmap", name = "True Phase", colorscale = "Viridis" ) %>% layout(title = "True Phase Distribution") hm_fft <- plot_ly( x = ds$time, y = ds$theta, z = t(phi_fft_aligned), type = "heatmap", name = "KPT-FFT Phase (Aligned)", colorscale = "Viridis" ) %>% layout(title = "KPT-FFT Estimated Phase (Aligned)") hm_gem <- plot_ly( x = ds$time, y = ds$theta, z = t(phi_gem_aligned), type = "heatmap", name = "KPT-GEM Phase (Aligned)", colorscale = "Viridis" ) %>% layout(title = "KPT-GEM Estimated Phase (Aligned)") subplot(hm_true, hm_fft, hm_gem, nrows = 3) %>% layout(title = "Double-Slit: Phase Distribution Comparison (Aligned)") ``` ## Kime-Surface Comparison (Aligned) ```{r surface-comparison-aligned} # 3D surface plots p_true <- plot_ly( x = ds$time, y = ds$theta, z = t(ds$S_true), type = "surface", scene = "scene1", colorscale = "Viridis", opacity = 0.95 ) %>% layout(title = "True Kime-Surface") p_gem <- plot_ly( x = ds$time, y = ds$theta, z = t(S_gem_aligned), type = "surface", scene = "scene2", colorscale = "Blues", opacity = 0.8 ) %>% layout(title = "KPT-GEM Surface (Aligned)") p_fft <- plot_ly( x = ds$time, y = ds$theta, z = t(S_fft_aligned), type = "surface", scene = "scene3", colorscale = "Reds", opacity = 0.8 ) %>% layout(title = "KPT-FFT Surface (Aligned)") # Combine with synchronized cameras combined_plot <- subplot(p_true, p_gem, p_fft) %>% layout( title = "Double-Slit: Kime-Surface Comparison (Aligned)", scene = list( domain = list(x = c(0, 0.33), y = c(0, 1)), aspectmode = "cube" ), scene2 = list( domain = list(x = c(0.33, 0.66), y = c(0, 1)), aspectmode = "cube" ), scene3 = list( domain = list(x = c(0.66, 1), y = c(0, 1)), aspectmode = "cube" ) ) %>% htmlwidgets::onRender( "function(x, el) { x.on('plotly_relayout', function(d) { const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key)); if (camera.length) { const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key)); const new_layout = {}; scenes.forEach(key => { new_layout[key] = {...x.layout[key], camera: {...d[camera]}}; }); Plotly.relayout(x, new_layout); } }); }" ) combined_plot ``` ## Performance Analysis ```{r performance-analysis} # Time-resolved metrics plot pm_combined <- rbind( pm_gem %>% mutate(Algorithm = "KPT-GEM"), pm_fft %>% mutate(Algorithm = "KPT-FFT") ) p_metrics <- plot_ly(pm_combined, x = ~time) %>% add_trace(y = ~JSD_bits, color = ~Algorithm, type = 'scatter', mode = 'lines', name = ~paste(Algorithm, "JSD")) %>% layout( title = "Double-Slit: Phase Recovery Performance Over Time", xaxis = list(title = "Time"), yaxis = list(title = "JSD (bits)"), hovermode = "x unified" ) p_metrics # Final summary cat("\n=== FINAL RESULTS (Double-Slit with Alignment) ===\n") print(summary_ds) cat("\nInterpretation:\n") cat("- JSD near 0: Perfect recovery\n") cat("- Hellinger in [0,1]: Lower is better\n") cat("- Total Variation in [0,1]: Lower is better\n") cat("- Relative L2 < 1: Good surface recovery\n") ``` ## Key Observations The alignment step is *critical* for proper evaluation because: 1. *Phase Ambiguity*: KPT recovers $M = F \cdot \Phi$, but any rotation $\theta \to \theta + \alpha$ that shifts $\Phi$ by $\alpha$ and counter-rotates $F$ by $-\alpha$ leaves $M$ unchanged. 2. *Alignment Method*: We use the first trigonometric moment $m_1 = \mathbb{E}[e^{i\Theta}]$ to identify the "center of mass" of the phase distribution and rotate to a canonical frame. 3. *Impact*: Without alignment, the metrics can be artificially inflated. The aligned results show the true recovery performance. 4. *Verification*: After alignment, the phase distributions should have their main modes near the same locations as the true distribution, and the surfaces should overlap properly. # Experiment 2: fMRI ON/OFF Design The next experiment shows a highly realistic fMRI simulation of an *observable* $y=S(t,\Theta)+\varepsilon$, which does actually depend on the time dynamic phase $\theta$. It mirrors the Double-Slit physics experiment setup; sampling $\Theta$ from a time-varying *von Mises phase law* and synthetically generating $Y = S(t,\Theta) + \varepsilon$. This fMRI simulation uses a design/HRF model, removes external dependencies, and returns the same structure (`Y, time, theta, phi_true, S_true`), plus `design` and `A` for reference (*OFF*, rest state). ## fMRI Simulation with Phase-Dependent BOLD Signal ```{r simulate-fmri-advanced} # Helper: von Mises density dvonmises <- function(theta, mu, kappa) { exp(kappa * cos(theta - mu)) / (2*pi * besselI(kappa, 0)) } simulate_fmri_advanced <- function(P = 10, R = 20, TR = 2, duration = 300, L_theta = 256, noise_sd = 0.2, kappa_on = 6, kappa_off = 1, seed = 2025) { set.seed(seed) # Time and phase grids K <- as.integer(duration / TR) time <- seq(0, by = TR, length.out = K) theta <- theta_grid(L_theta) # Block design: 30s ON/OFF blocks block_len_sec <- 30 samples_per_block <- block_len_sec / TR design <- rep(c(0, 1), each = samples_per_block) design <- rep_len(design, K) # Hemodynamic Response Function (canonical double-gamma) hrf <- function(t) { a1 <- 6; a2 <- 16 b1 <- 1; b2 <- 1 c <- 1/6 h <- (t^(a1-1) * b1^a1 * exp(-b1*t)) / gamma(a1) - c * (t^(a2-1) * b2^a2 * exp(-b2*t)) / gamma(a2) h[t < 0] <- 0 h } # Convolve design with HRF hrf_vals <- hrf(seq(0, 30, by = TR)) bold <- stats::filter(design, hrf_vals, method = "convolution", sides = 1) bold <- as.numeric(bold) bold[is.na(bold)] <- 0 # Normalize BOLD to [0.2, 1.0] if (max(bold) > min(bold)) { bold <- (bold - min(bold)) / (max(bold) - min(bold)) bold <- 0.2 + 0.8 * bold } else { bold <- rep(0.6, K) } # Small baseline drift baseline <- 0.05 * sin(2*pi * time / max(time)) # Build true kime-surface: baseline + BOLD modulation # S(t,θ) = baseline(t) + A(t)*cos(θ) + 0.25*A(t)*cos(2θ) S_true <- matrix(0, K, L_theta) for (k in 1:K) { S_true[k,] <- baseline[k] + bold[k] * cos(theta) + 0.25 * bold[k] * cos(2*theta) } # Time-varying phase distribution # During ON: higher concentration (more synchronous) # During OFF: lower concentration (more variable) mu_theta <- 0.2 * sin(2*pi * time / max(time)) # Slight drift in mean phi_true <- matrix(0, K, L_theta) for (k in 1:K) { # Concentration depends on ON/OFF state kappa <- if (design[k] == 1) kappa_on else kappa_off mu <- mu_theta[k] # von Mises distribution p <- dvonmises(theta, mu, kappa) phi_true[k,] <- p / mean(p) # Normalize to mean=1 } # Generate observations Y = S(t, Theta) + noise I <- P * R # Total number of observations per time Y <- matrix(0, I, K) for (k in 1:K) { # Sample phases from true distribution probs <- pmax(phi_true[k,], 1e-12) idx <- sample.int(L_theta, I, replace = TRUE, prob = probs) # Generate BOLD signals Y[,k] <- S_true[k, idx] + rnorm(I, 0, noise_sd) } list( Y = Y, time = time, theta = theta, phi_true = phi_true, S_true = S_true, design = design, bold = bold, baseline = baseline, params = list(P = P, R = R, TR = TR, duration = duration, L_theta = L_theta, noise_sd = noise_sd, kappa_on = kappa_on, kappa_off = kappa_off) ) } # Generate fMRI data fm <- simulate_fmri_advanced(P = 10, R = 20, TR = 2, duration = 300) cat("fMRI simulation complete:\n") cat(" Subjects × Runs:", fm$params$P, "×", fm$params$R, "\n") cat(" Time points:", ncol(fm$Y), "\n") cat(" Phase resolution:", length(fm$theta), "\n") ``` ## Run KPT Analysis on fMRI Data ```{r run-fmri-kpt} # Run KPT-GEM (fixed surface) fit_fm_gem <- kpt_unified( fm$Y, time = fm$time, J = 12, L_theta = length(fm$theta), lambda_phi = 1e-2, p_phi = 1, lambda_F = 1e-2, p_F = 1, algorithm = "GEM", update_surface = FALSE, max_iter = 25, tol = 8e-4, verbose = TRUE ) # Run KPT-FFT (alternating updates) fit_fm_fft <- kpt_unified( fm$Y, time = fm$time, J = 12, L_theta = length(fm$theta), lambda_phi = 1e-2, p_phi = 1, lambda_F = 1e-2, p_F = 1, algorithm = "FFT", update_surface = TRUE, max_iter = 25, tol = 8e-4, verbose = TRUE ) # CRITICAL: Compute alignment from estimated phases # Extract first trigonometric moments m1_fft <- apply(fit_fm_fft$varphi_theta, 1L, function(row) { mean(exp(1i * fm$theta) * row) }) delta_fft <- Arg(m1_fft) m1_gem <- apply(fit_fm_gem$varphi_theta, 1L, function(row) { mean(exp(1i * fm$theta) * row) }) delta_gem <- Arg(m1_gem) # Apply alignment to phases phi_fft_aligned <- anchor_phase_by_m1(fit_fm_fft$varphi_theta, fm$theta) phi_gem_aligned <- anchor_phase_by_m1(fit_fm_gem$varphi_theta, fm$theta) # Apply same rotation to surfaces S_fft_aligned <- rotate_surface_by_delta(fit_fm_fft$Sgrid, fm$theta, delta_fft) S_gem_aligned <- rotate_surface_by_delta(fit_fm_gem$Sgrid, fm$theta, delta_gem) # Verification cat("\n=== Alignment Results ===\n") cat("Mean rotation FFT:", mean(delta_fft), "rad\n") cat("Mean rotation GEM:", mean(delta_gem), "rad\n") cat("Phase normalization check:\n") cat(" FFT:", mean(rowMeans(phi_fft_aligned)), "\n") cat(" GEM:", mean(rowMeans(phi_gem_aligned)), "\n") ``` ## Quantitative Evaluation ```{r evaluate-fmri} # Compute metrics on ALIGNED distributions pm_gem <- phase_metrics_df(fm$phi_true, phi_gem_aligned, fm$theta, fm$time) pm_fft <- phase_metrics_df(fm$phi_true, phi_fft_aligned, fm$theta, fm$time) # Surface errors on ALIGNED surfaces m_gem <- surface_L2_cone(fm$S_true, S_gem_aligned, fm$time, fm$theta) m_fft <- surface_L2_cone(fm$S_true, S_fft_aligned, fm$time, fm$theta) # Summary table summary_fmri <- data.frame( Algorithm = c("KPT-GEM", "KPT-FFT"), Mean_JSD_bits = c(mean(pm_gem$JSD_bits), mean(pm_fft$JSD_bits)), Mean_Hellinger = c(mean(pm_gem$Hellinger), mean(pm_fft$Hellinger)), Mean_TV = c(mean(pm_gem$TV), mean(pm_fft$TV)), Surface_RelL2 = c(m_gem$rel_L2, m_fft$rel_L2) ) kable(summary_fmri, digits = 4, caption = "fMRI ON/OFF: Performance Metrics (After Alignment)") # Compare before/after alignment pm_fft_unaligned <- phase_metrics_df(fm$phi_true, fit_fm_fft$varphi_theta, fm$theta, fm$time) cat("\n=== Alignment Impact ===\n") cat("FFT Mean JSD before alignment:", mean(pm_fft_unaligned$JSD_bits), "bits\n") cat("FFT Mean JSD after alignment:", mean(pm_fft$JSD_bits), "bits\n") cat("Improvement:", mean(pm_fft_unaligned$JSD_bits) - mean(pm_fft$JSD_bits), "bits\n") ``` ## ON vs OFF Block Analysis ```{r on-off-analysis} # Separate ON and OFF blocks idx_on <- which(fm$design == 1) idx_off <- which(fm$design == 0) # Metrics for ON blocks pm_gem_on <- pm_gem[idx_on,] pm_fft_on <- pm_fft[idx_on,] # Metrics for OFF blocks pm_gem_off <- pm_gem[idx_off,] pm_fft_off <- pm_fft[idx_off,] # Summary by condition summary_by_block <- data.frame( Algorithm = rep(c("KPT-GEM", "KPT-FFT"), 2), Block = rep(c("ON", "OFF"), each = 2), Mean_JSD = c(mean(pm_gem_on$JSD_bits), mean(pm_fft_on$JSD_bits), mean(pm_gem_off$JSD_bits), mean(pm_fft_off$JSD_bits)), Mean_TV = c(mean(pm_gem_on$TV), mean(pm_fft_on$TV), mean(pm_gem_off$TV), mean(pm_fft_off$TV)) ) kable(summary_by_block, digits = 4, caption = "fMRI: Performance by Block Type") # Representative time points idx_on_peak <- idx_on[which.max(fm$bold[idx_on])] idx_off_baseline <- idx_off[which.min(fm$bold[idx_off])] cat("\n=== Selected Time Points ===\n") cat("Peak ON at t =", fm$time[idx_on_peak], "s\n") cat("Baseline OFF at t =", fm$time[idx_off_baseline], "s\n") ``` ## Visualization ```{r visualize-fmri-phases} # Plot phases at ON and OFF peaks p_on <- plot_phase_overlay( fm$theta, fm$phi_true[idx_on_peak,], phi_gem_aligned[idx_on_peak,], phi_fft_aligned[idx_on_peak,], t_lab = sprintf("ON Block (t = %.0fs)", fm$time[idx_on_peak]), title = "fMRI: Phase Distribution During Stimulus" ) p_off <- plot_phase_overlay( fm$theta, fm$phi_true[idx_off_baseline,], phi_gem_aligned[idx_off_baseline,], phi_fft_aligned[idx_off_baseline,], t_lab = sprintf("OFF Block (t = %.0fs)", fm$time[idx_off_baseline]), title = "fMRI: Phase Distribution During Rest" ) print(p_on) print(p_off) # Additional time points time_indices <- seq(1, length(fm$time), by = 30) # Every 60s for (k in time_indices[1:3]) { block_type <- if(fm$design[k] == 1) "ON" else "OFF" p <- plot_phase_overlay( fm$theta, fm$phi_true[k,], phi_gem_aligned[k,], phi_fft_aligned[k,], t_lab = sprintf("%s Block (t = %.0fs)", block_type, fm$time[k]), title = "fMRI: Phase Distribution" ) print(p) } ``` ## Phase Distribution Heatmaps ```{r fmri-phase-heatmaps} # Create heatmaps with design overlay hm_true <- plot_ly( x = fm$time, y = fm$theta, z = t(fm$phi_true), type = "heatmap", colorscale = "Viridis", colorbar = list(title = "Density") ) %>% layout(title = "True Phase Distribution") hm_gem <- plot_ly( x = fm$time, y = fm$theta, z = t(phi_gem_aligned), type = "heatmap", colorscale = "Viridis", colorbar = list(title = "Density") ) %>% layout(title = "KPT-GEM Phase (Aligned)") hm_fft <- plot_ly( x = fm$time, y = fm$theta, z = t(phi_fft_aligned), type = "heatmap", colorscale = "Viridis", colorbar = list(title = "Density") ) %>% layout(title = "KPT-FFT Phase (Aligned)") subplot(hm_true, hm_gem, hm_fft, nrows = 3) %>% layout(title = "fMRI: Phase Distributions Over Time") ``` ## Kime-Surface Comparison ```{r fmri-surfaces} # 3D surface plots p_true <- plot_ly( x = fm$time, y = fm$theta, z = t(fm$S_true), type = "surface", scene = "scene1", colorscale = "Viridis", opacity = 0.95 ) %>% layout(title = "True BOLD Surface") p_gem <- plot_ly( x = fm$time, y = fm$theta, z = t(S_gem_aligned), type = "surface", scene = "scene2", colorscale = "Blues", opacity = 0.8 ) %>% layout(title = "KPT-GEM Surface") p_fft <- plot_ly( x = fm$time, y = fm$theta, z = t(S_fft_aligned), type = "surface", scene = "scene3", colorscale = "Reds", opacity = 0.8 ) %>% layout(title = "KPT-FFT Surface") # Combined plot with synchronized cameras combined_surface <- subplot(p_true, p_gem, p_fft) %>% layout( title = "fMRI: Kime-Surface Comparison (Aligned)", scene = list( domain = list(x = c(0, 0.33), y = c(0, 1)), aspectmode = "cube", camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5)) ), scene2 = list( domain = list(x = c(0.33, 0.66), y = c(0, 1)), aspectmode = "cube", camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5)) ), scene3 = list( domain = list(x = c(0.66, 1), y = c(0, 1)), aspectmode = "cube", camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5)) ) ) %>% htmlwidgets::onRender( "function(x, el) { x.on('plotly_relayout', function(d) { const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key)); if (camera.length) { const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key)); const new_layout = {}; scenes.forEach(key => { new_layout[key] = {...x.layout[key], camera: {...d[camera]}}; }); Plotly.relayout(x, new_layout); } }); }" ) combined_surface ``` ## Time-Resolved Performance ```{r fmri-time-performance} # Create design indicator for visualization design_df <- data.frame( time = fm$time, design = fm$design, bold = fm$bold ) # Plot metrics with design overlay time_diff <- mean(diff(design_df$time)) p_metrics <- plot_ly() %>% add_segments( data = design_df[design_df$design == 1,], x = ~time, xend = ~time + time_diff, y = 0, yend = 1, line = list(width = 15, color = 'rgba(200,200,200,0.2)'), showlegend = FALSE, hoverinfo = 'skip' ) %>% # Add metric traces add_trace( data = pm_gem, x = ~time, y = ~JSD_bits, type = 'scatter', mode = 'lines', name = 'KPT-GEM', line = list(color = 'blue', width = 2) ) %>% add_trace( data = pm_fft, x = ~time, y = ~JSD_bits, type = 'scatter', mode = 'lines', name = 'KPT-FFT', line = list(color = 'red', width = 2) ) %>% layout( title = "fMRI: Phase Recovery Performance (ON blocks shaded)", xaxis = list(title = "Time (s)"), yaxis = list(title = "JSD (bits)", range = c(0, max(c(pm_gem$JSD_bits, pm_fft$JSD_bits)))), hovermode = "x unified" ) p_metrics # Surface correlation over time cors_gem <- sapply(1:length(fm$time), function(k) { cor(fm$S_true[k,], S_gem_aligned[k,]) }) cors_fft <- sapply(1:length(fm$time), function(k) { cor(fm$S_true[k,], S_fft_aligned[k,]) }) cat("\n=== Surface Correlation ===\n") cat("GEM: mean =", mean(cors_gem), ", sd =", sd(cors_gem), "\n") cat("FFT: mean =", mean(cors_fft), ", sd =", sd(cors_fft), "\n") ``` ## Final Summary ```{r fmri-final-summary} # Consolidated results cat("\n=== FINAL RESULTS (fMRI with Alignment) ===\n") print(summary_fmri) kable(summary_fmri, digits = 4, caption = "fMRI: Performance by Block Type") cat("\n=== Performance by Block Type ===\n") print(summary_by_block) kable(summary_by_block, digits = 4, caption = "fMRI: Performance by Block Type") cat("\n=== Key Findings ===\n") cat("1. Phase concentration higher during ON blocks (kappa =", fm$params$kappa_on, ")\n") cat("2. Phase concentration lower during OFF blocks (kappa =", fm$params$kappa_off, ")\n") cat("3. Both algorithms recover this structure after alignment\n") cat("4. FFT (with surface update) performs slightly better overall\n") # Performance interpretation cat("\n=== Metric Interpretation ===\n") cat("- JSD = 0: Perfect recovery (identical distributions)\n") cat("- JSD < 0.1 bits: Excellent recovery\n") cat("- JSD < 0.5 bits: Good recovery\n") cat("- Relative L2 < 1: Better than baseline\n") ```