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 Overview & Conventions

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

Double-slit physics experiment (latent phase drives interference).
fMRI ON vs. OFF event-related design (latent phase modulates BOLD fMRI neuroscience study).

We use the circular phase normalization \(\int_0^{2\pi}\varphi_t(\theta)\,\frac{d\theta}{2\pi}=1\). On a uniform grid of length \(L_\theta\), this is equivalent to the row-mean equals 1 constraint. The spectral updates use the Wiener-Sobolev filter.

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

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

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

2 Utilities (Fourier, projection, metrics, visuals)

We start by defining the core utility functions.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 KPT Algorithms

All methematical, computaitonal, and statistical details are provided in the KPT paper.

3.1 KPT-GEM (penalized tomography with stabilized surface)

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

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

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

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

Repeat:

E-step: compute \(\xi^{(\ell)}_{j,k}(n)\) and \(\widehat m^{(\ell)}_n(t_k)=\tfrac{1}{N}\sum_j y_{j,k}\xi^{(\ell)}_{j,k}(n)\)
\(\Phi\)-step (phase update): For each \(t_k\), FFT in \(n\) to get \(\widehat M^{(\ell)}(\omega,t_k)\) and \(\widehat F^{(\ell)}(\omega,t_k)\); set

\[\widehat\Phi^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat F^{(\ell)}(\omega,t_k)}}{|\widehat F^{(\ell)}(\omega,t_k)|^2+\lambda_\Phi\,\Lambda_{p_\Phi}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k).\]

IFFT in \(\omega\) and project: \(\widehat\varphi^{(\ell+1)}_{t_k}\leftarrow \Pi_{\mathcal P}\big[\mathcal{F}^{-1}\widehat\Phi^{(\ell+1)}(\cdot,t_k)\big]\). Optionally smooth in \(t\).
\(F\)-step (surface update): For each \(t_k\), recompute \(\widehat\Phi^{(\ell+1)}\) and set

\[\widehat F^{(\ell+1)}(\omega,t_k)=\frac{\overline{\widehat \Phi^{(\ell+1)}(\omega,t_k)}}{|\widehat \Phi^{(\ell+1)}(\omega,t_k)|^2+\lambda_F\,\Lambda_{p_F}(\omega)}\,\widehat M^{(\ell)}(\omega,t_k).\]

IFFT in \(\omega\) to obtain \(\widehat f^{(\ell+1)}_n(t_k)\); enforce \(f_{-n}=\overline{f_n}\). Optionally smooth in \(t\).
Synthesis: \(\widehat{\mathcal S}^{(\ell+1)}(t_k,\theta)=\sum_{n=-J}^J \widehat f^{(\ell+1)}_n(t_k)e^{in\theta}\)
Convergence check: stop when \(\max_k\|\widehat\varphi^{(\ell+1)}_{t_k}-\widehat\varphi^{(\ell)}_{t_k}\|_{L^2}<\varepsilon\) and \(\max_{k,n}|\widehat f^{(\ell+1)}_n(t_k)-\widehat f^{(\ell)}_n(t_k)|<\varepsilon\)

Until converged

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

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

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

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

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

    # M-step: spectral deconvolution per time using FFT evaluation on unit circle
    phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      F_w <- coeffs_to_Xomega_fft(f_n[k, ], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(F_w)^2) + lambda * Lambda_p + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      phi_hat_n[k, ] <- Xomega_to_coeffs_fft(Phi_w, nvec)
    }
    # Synthesize φ and project
    phi_new <- matrix(0, nrow = K, ncol = L_theta)
    Epos <- exp(1i * outer(nvec, theta))
    for (k in 1:K) {
      phi_row <- Re(drop(t(Epos) %*% phi_hat_n[k, ]))
      phi_new[k, ] <- phi_row
    }
    phi_new <- row_project_simplex(phi_new)

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

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

3.2 KPT-FFT (alternating phase & surface updates)

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

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

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

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

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

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

Repeat:

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

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

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

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

  # Initialize: uniform φ, simple S proxy
  phi_grid <- matrix(1.0, nrow = K, ncol = L_theta)
  Sgrid <- build_S_from_data(Y, seq_len(K), theta)
  f_n <- matrix(0+0i, nrow = K, ncol = Nn)
  for (k in 1:K) {
    Eneg <- exp(-1i * outer(theta, nvec))
    f_n[k, ] <- as.vector((t(Sgrid[k, ]) %*% Eneg) / L_theta)
    f_n[k, ] <- enforce_reality_coeffs(f_n[k, ], nvec)
  }

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

  obj <- numeric(max_iter)

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

    # ---- Φ-step (phase update in spectral domain via FFT) ----
    phi_hat_n <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      F_w <- coeffs_to_Xomega_fft(f_n[k, ], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(F_w)^2) + lambda_phi * Lambda_phi + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      phi_hat_n[k, ] <- Xomega_to_coeffs_fft(Phi_w, nvec)
    }
    # Synthesize φ(θ) & project
    phi_new <- matrix(0, nrow = K, ncol = L_theta)
    Epos <- exp(1i * outer(nvec, theta))
    for (k in 1:K) {
      phi_row <- Re(drop(t(Epos) %*% phi_hat_n[k, ]))
      phi_new[k, ] <- phi_row
    }
    phi_new <- row_project_simplex(phi_new)

    # ---- F-step (surface update in spectral domain via FFT) ----
    f_n_new <- matrix(0+0i, nrow = K, ncol = Nn)
    for (k in 1:K) {
      Phi_w <- coeffs_to_Xomega_fft(phi_hat_n[k, ], nvec)
      M_w   <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- (Mod(Phi_w)^2) + lambda_F * Lambda_F + 1e-10
      F_w   <- Conj(Phi_w) * M_w / denom
      f_n_new[k, ] <- enforce_reality_coeffs(Xomega_to_coeffs_fft(F_w, nvec), nvec)
    }
    # Synthesize S(t,θ)
    S_new <- matrix(0, nrow = K, ncol = L_theta)
    for (k in 1:K) S_new[k, ] <- surface_from_coeffs(f_n_new[k, ], nvec, theta)

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

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

4 Simulation 1: Double-Slit Physics Experiment

We simulate an interference-like kime-surface

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

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

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

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

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

  # Time-varying von Mises mixture over theta (two drifting lobes)
  mu1 <- (2*pi*(0.2 + 0.6*t)) %% (2*pi)
  mu2 <- (2*pi*(0.9 - 0.6*t)) %% (2*pi)
  w   <- 0.5 + 0.4*sin(2*pi*t)             # mixture weight
  phi_true <- matrix(0, nrow = K, ncol = L_theta)
  for (k in 1:K) {
    p <- w[k]*dvonmises(theta, mu1[k], kappa_phase) +
         (1-w[k])*dvonmises(theta, mu2[k], kappa_phase)
    phi_true[k, ] <- p / mean(p)           # mean == 1 (∫dθ/(2π)=1)
  }

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

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

## List of 5
##  $ Y       : num [1:256, 1:256] -0.24023 0.04252 0.11409 -0.09725 0.00726 ...
##  $ time    : num [1:256] 0 0.00392 0.00784 0.01176 0.01569 ...
##  $ theta   : num [1:256] 0 0.0245 0.0491 0.0736 0.0982 ...
##  $ phi_true: num [1:256, 1:256] 1.001 0.932 0.865 0.803 0.744 ...
##  $ S_true  : num [1:256, 1:256] 0.00714 0.01253 0.01791 0.02327 0.0286 ...

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

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

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

4.2 Double-Slit Results: Quantitative & Qualitative Evaluation

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

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

summary_ds <- tibble(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(pm_gem$JSD_bits), mean(pm_fft$JSD_bits)),
  Mean_Hellinger = c(mean(pm_gem$Hellinger), mean(pm_fft$Hellinger)),
  Mean_TV       = c(mean(pm_gem$TV),       mean(pm_fft$TV)),
  Rel_L2_Surface= c(m_gem$rel_L2,          m_fft$rel_L2)
)
kable(summary_ds, digits = 4, caption = "Double-Slit: Summary metrics: True Kimesurface vs. KPT-GEM and KPT-FFT")

Double-Slit: Summary metrics: True Kimesurface vs. KPT-GEM and KPT-FFT
Algorithm	Mean_JSD_bits	Mean_Hellinger	Mean_TV	Rel_L2_Surface
KPT-GEM	0.6072	0.7266	0.7532	1.2514
KPT-FFT	0.6101	0.7302	0.7524	1.5167

# Additional debugging
cat("Length check:\n")

## Length check:

cat("ds$theta:", length(ds$theta), "\n")

## ds$theta: 256

cat("fit_ds_fft$theta:", length(fit_ds_fft$theta), "\n")

## fit_ds_fft$theta: 256

cat("ds$phi_true row:", length(ds$phi_true[1,]), "\n")

## ds$phi_true row: 256

cat("fit_ds_fft$varphi_theta row:", length(fit_ds_fft$varphi_theta[1,]), "\n")

## fit_ds_fft$varphi_theta row: 256

# Check if they're the same
all.equal(length(ds$theta), length(fit_ds_fft$theta))

## [1] TRUE

all.equal(length(ds$phi_true[1,]), length(fit_ds_fft$varphi_theta[1,]))

## [1] TRUE

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

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

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

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

Time Complexity, Inferential uncertainty and Spacekime Analytics

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

SOCR Team

10/06/2025