| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
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
01/03/2026
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:## Mean rotation angle (FFT): -0.1223663 rad## Mean rotation angle (GEM): -0.1884956 rad## Phase normalization check (should be ~1):## FFT aligned: 1## GEM aligned: 14.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)")| 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 |
##
## === 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## FFT Mean JSD after alignment: 0.5129105 bits## Improvement: -0.005489682 bits4.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)")