---
title: "Appendix 03: Kime-Phase Tomography (KPT V6)"
subtitle: "Self-contained computational protocol aligned with KPT_revised_paper_body_v1.tex"
author: "SOCR Team"
date: "`r format(Sys.Date(), '%Y-%m-%d')`"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
    code_folding: hide
    df_print: paged
---

This notebook is a complete revision of the previous KPT appendix. It is aligned with the rewritten manuscript `KPT_revised_paper_body_v1.tex` and is organized to mirror the paper's structure: introduction, mathematical foundations, KPT algorithms, statistical inference, applications, and conclusions. The notebook is fully self-contained, uses the **time-specific latent phase model** throughout, adopts the normalization convention based on the Haar probability measure $d\theta/(2\pi)$, and removes the inconsistencies and exploratory sections present in earlier drafts.

The core generative model used everywhere in this protocol is
$$
Y_{j,k} = \mathcal S(t_k, \Theta_j(t_k)) + \varepsilon_{j,k},
\qquad
\Theta_j(t_k) \sim \varphi_{t_k},
\qquad
\varepsilon_{j,k} \sim \mathcal N(0,\sigma^2),
$$
where each phase law $\varphi_{t_k}$ is a density on the circle with respect to $d\theta/(2\pi)$. On a discrete phase grid $\{\theta_\ell\}_{\ell=0}^{L_\theta-1}$, the corresponding normalization is
$$
\frac{1}{L_\theta}\sum_{\ell=0}^{L_\theta-1} \varphi_{t_k}(\theta_\ell) = 1.
$$
Accordingly, the uniform phase law is represented by a vector of ones.

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 9,
  fig.height = 5.5,
  dpi = 120
)

required_pkgs <- c("ggplot2", "dplyr", "tidyr", "purrr", "knitr", "tibble")
missing_pkgs <- required_pkgs[!vapply(required_pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_pkgs) > 0) {
  stop(
    paste0(
      "Please install the following packages before rendering this notebook: ",
      paste(missing_pkgs, collapse = ", ")
    )
  )
}

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(purrr)
  library(knitr)
  library(tibble)
})

set.seed(20260322)
```

# Introduction {#sec-introduction}

This appendix implements the computational procedures described in the revised paper. The emphasis here is practical rather than theorem-proving. In particular, the notebook implements the **projected practical algorithms** used in the numerical experiments:

- **KPT-GEM (practical phase-update variant):** Wiener--Sobolev phase updates with a fixed initialized surface.
- **KPT-FFT (practical alternating variant):** projected phase updates alternating with surface updates.

The revised paper distinguishes these computational routines from the **exact constrained EM reference algorithm**, which is the appropriate object for strict monotonicity claims. That exact reference algorithm is not the default numerical engine in this notebook. This notebook therefore reports empirical convergence diagnostics for the practical projected routines and does **not** interpret them as globally monotone EM algorithms.

# Mathematical Foundations {#sec-methods}

## Model, normalization, and alignment conventions

Throughout the notebook, phase densities are normalized with respect to the Haar probability measure $d\theta/(2\pi)$. On the discrete grid this means **row mean = 1**. All simulation, estimation, alignment, prediction, and evaluation steps below use the same convention.

Because the latent phase origin is a gauge degree of freedom, comparisons between true and estimated phase laws are performed after an explicit anchoring step. In the numerical studies below, we use the first trigonometric moment anchor whenever the first harmonic is non-negligible; if the first moment is nearly zero, the row is left unchanged.

## Simulation and fitting parameters copied from the paper

```{r parameter-table}
simulation_parameters <- tibble::tribble(
  ~Setting, ~Double_slit, ~fMRI,
  "Observation times",
  "K = 50, equally spaced on [0, 1]",
  "K = 150, sampled every 2 seconds on [0, 298]",
  "Phase grid",
  "L_theta = 256",
  "L_theta = 256",
  "Replicates per time point",
  "N = 100",
  "P x R = 10 x 20 = 200",
  "Noise standard deviation",
  "sigma = 0.15",
  "sigma = 0.20",
  "Latent phase law",
  "Two-component von Mises mixture with kappa_phase = 4",
  "von Mises with kappa_on = 6 and kappa_off = 1",
  "Surface generator",
  "Baseline + Gaussian envelope + first/second harmonics with kappa_wave = 12",
  "Baseline drift + HRF-convolved block design + first/second harmonics",
  "KPT-GEM tuning",
  "J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = FALSE, max_iter = 30, tol = 1e-3",
  "J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = FALSE, max_iter = 25, tol = 8e-4",
  "KPT-FFT tuning",
  "J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = TRUE, max_iter = 30, tol = 1e-3",
  "J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = TRUE, max_iter = 25, tol = 8e-4",
  "Evaluation alignment",
  "First-moment anchor; the same rotation is applied to the surface",
  "First-moment anchor; the same rotation is applied to the surface"
)

knitr::kable(
  simulation_parameters,
  caption = "Simulation and fitting parameters used throughout this notebook. These match the revised paper."
)
```

# Kime-Phase Tomography (KPT) {#sec-kpt}

## Core utilities

```{r utilities}
# ---- Grid 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)

freq_index <- function(L) {
  if (L %% 2L == 0L) {
    c(0:(L/2L), -(L/2L - 1L):-1L)
  } else {
    c(0:((L - 1L)/2L), -((L - 1L)/2L):-1L)
  }
}

trapz_vec <- function(x, y) {
  sum(0.5 * (y[-1L] + y[-length(y)]) * diff(x))
}

# ---- Stable row-wise softmax --------------------------------------------------
softmax_rows <- function(log_mat) {
  m <- apply(log_mat, 1L, max)
  w <- exp(log_mat - m)
  w / pmax(rowSums(w), .Machine$double.eps)
}

# ---- Exact Euclidean projection onto {x >= 0, mean(x) = 1} --------------------
euclidean_simplex_project <- function(v, z = length(v)) {
  n <- length(v)
  u <- sort(v, decreasing = TRUE)
  cssv <- cumsum(u) - z
  rho <- max(which(u - cssv / seq_len(n) > 0))
  theta <- cssv[rho] / rho
  pmax(v - theta, 0)
}

row_project_simplex <- function(P) {
  P <- as.matrix(P)
  out <- t(apply(P, 1L, euclidean_simplex_project, z = ncol(P)))
  rownames(out) <- rownames(P)
  colnames(out) <- colnames(P)
  out
}

# ---- Fourier utilities --------------------------------------------------------
discrete_fourier_coeffs <- function(x_theta, theta, nvec) {
  Eneg <- exp(-1i * outer(theta, nvec))
  as.vector(drop(t(x_theta) %*% Eneg) / length(theta))
}

synthesize_from_coeffs <- function(coeffs_n, theta, nvec) {
  Epos <- exp(1i * outer(nvec, theta))
  as.vector(Re(drop(coeffs_n %*% Epos)))
}

enforce_reality_coeffs <- function(coeffs_n, nvec) {
  coeffs_n <- as.complex(coeffs_n)
  coeffs_n[nvec == 0] <- Re(coeffs_n[nvec == 0])
  J <- max(abs(nvec))
  for (j in seq_len(J)) {
    pos <- which(nvec == j)
    neg <- which(nvec == -j)
    if (length(pos) == 1L && length(neg) == 1L) {
      avg <- 0.5 * (coeffs_n[pos] + Conj(coeffs_n[neg]))
      coeffs_n[pos] <- avg
      coeffs_n[neg] <- Conj(avg)
    }
  }
  coeffs_n
}

coeffs_to_dft_order <- function(x_n, nvec) {
  N <- length(nvec)
  out <- complex(length = N)
  for (i in seq_along(nvec)) {
    m <- (nvec[i] + N) %% N
    out[m + 1L] <- x_n[i]
  }
  out
}

dft_order_to_coeffs <- function(x_m, nvec) {
  N <- length(nvec)
  out <- complex(length = N)
  for (i in seq_along(nvec)) {
    m <- (nvec[i] + N) %% N
    out[i] <- x_m[m + 1L]
  }
  out
}

coeffs_to_Xomega_fft <- function(x_n, nvec) {
  fft(coeffs_to_dft_order(x_n, nvec), inverse = FALSE)
}

Xomega_to_coeffs_fft <- function(X_w, nvec) {
  N <- length(nvec)
  dft_order_to_coeffs(fft(X_w, inverse = TRUE) / N, nvec)
}

sobolev_weights <- function(nvec, p = 1L) {
  abs(nvec)^(2 * p)
}

# ---- Periodic interpolation and circular shifts -------------------------------
periodic_interp <- function(theta, y, theta_new) {
  stopifnot(length(theta) == length(y))
  theta_ext <- c(theta, 2*pi)
  y_ext <- c(y, y[1L])
  xout <- theta_new %% (2*pi)
  approx(theta_ext, y_ext, xout = xout, rule = 2L, ties = mean)$y
}

shift_periodic_row <- function(x, delta) {
  L <- length(x)
  n <- freq_index(L)
  X <- fft(x)
  Re(fft(X * exp(-1i * n * delta), inverse = TRUE) / L)
}

shift_periodic_matrix <- function(X, delta_vec) {
  stopifnot(nrow(X) == length(delta_vec))
  out <- matrix(0, nrow = nrow(X), ncol = ncol(X))
  for (k in seq_len(nrow(X))) out[k, ] <- shift_periodic_row(X[k, ], delta_vec[k])
  out
}

# ---- Alignment helpers --------------------------------------------------------
first_trig_moment <- function(phi_mat, theta, harmonic = 1L) {
  apply(phi_mat, 1L, function(row) mean(row * exp(1i * harmonic * theta)))
}

anchor_by_first_moment <- function(phi_mat, theta, harmonic = 1L, tol = 1e-8) {
  m <- first_trig_moment(phi_mat, theta, harmonic = harmonic)
  delta <- ifelse(Mod(m) > tol, Arg(m) / harmonic, 0)
  phi_rot <- shift_periodic_matrix(phi_mat, delta)
  phi_rot <- row_project_simplex(phi_rot)
  list(phi = phi_rot, delta = delta, moment = m)
}

apply_anchor_to_surface <- function(S_mat, delta_vec) {
  shift_periodic_matrix(S_mat, delta_vec)
}

# ---- Surface initialization ---------------------------------------------------
initialize_surface <- function(Y, time, theta, method = c("empirical", "constant")) {
  method <- match.arg(method)
  K <- ncol(Y)
  L <- length(theta)
  if (method == "constant") return(matrix(mean(Y), nrow = K, ncol = L))

  mu_k <- colMeans(Y)
  sd_k <- apply(Y, 2L, stats::sd)

  mu_s <- if (K >= 5L) stats::smooth.spline(time, mu_k, spar = 0.6)$y else mu_k
  sd_s <- if (K >= 5L) stats::smooth.spline(time, sd_k, spar = 0.6)$y else sd_k

  S <- matrix(0, nrow = K, ncol = L)
  for (k in seq_len(K)) {
    S[k, ] <- mu_s[k] + sd_s[k] * cos(theta) + 0.25 * sd_s[k] * cos(2 * theta)
  }
  S
}

# ---- Criterion tracking (observed-data log-likelihood) ------------------------
observed_loglik_gaussian <- function(Y, Sgrid, phi_grid, sigma) {
  N <- nrow(Y)
  K <- ncol(Y)
  ll <- 0
  for (k in seq_len(K)) {
    weights <- pmax(phi_grid[k, ], 0)
    weights <- weights / sum(weights)
    dens_mat <- outer(Y[, k], Sgrid[k, ], function(y, mu) stats::dnorm(y, mean = mu, sd = sigma))
    mix_dens <- drop(dens_mat %*% weights)
    ll <- ll + sum(log(pmax(mix_dens, 1e-300)))
  }
  ll
}

# ---- Metrics -----------------------------------------------------------------
normalize_rows_prob <- function(P) {
  P <- pmax(P, 1e-15)
  sweep(P, 1L, rowSums(P), "/")
}

phase_metrics_df <- function(phi_true, phi_est, time, base2 = TRUE) {
  P <- normalize_rows_prob(phi_true)
  Q <- normalize_rows_prob(phi_est)
  log_fun <- if (base2) function(x) log2(pmax(x, 1e-15)) else function(x) log(pmax(x, 1e-15))
  K <- nrow(P)
  out <- tibble(
    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 seq_len(K)) {
    M <- 0.5 * (P[k, ] + Q[k, ])
    out$KL_true_est[k] <- sum(P[k, ] * log_fun(P[k, ] / Q[k, ]))
    out$KL_est_true[k] <- sum(Q[k, ] * log_fun(Q[k, ] / P[k, ]))
    out$JSD_bits[k] <- 0.5 * sum(P[k, ] * log_fun(P[k, ] / M)) +
      0.5 * sum(Q[k, ] * log_fun(Q[k, ] / M))
    out$Hellinger[k] <- sqrt(0.5 * sum((sqrt(P[k, ]) - sqrt(Q[k, ]))^2))
    out$TV[k] <- 0.5 * sum(abs(P[k, ] - Q[k, ]))
  }
  out
}

surface_L2_cone <- function(S_true, S_est, time) {
  stopifnot(all(dim(S_true) == dim(S_est)))
  err_theta <- rowMeans((S_true - S_est)^2)
  true_theta <- rowMeans(S_true^2)
  integrand_err <- time * err_theta
  integrand_true <- time * true_theta
  abs_sq <- trapz_vec(time, integrand_err)
  true_sq <- trapz_vec(time, integrand_true)
  list(
    abs_L2 = sqrt(max(abs_sq, 0)),
    rel_L2 = sqrt(max(abs_sq, 0) / max(true_sq, 1e-15)),
    integrand = tibble(time = time, mse_theta = err_theta)
  )
}

# ---- Predictive summaries -----------------------------------------------------
kpt_predictive_bands <- function(fit, B = 4000L,
                                 probs = c(0.05, 0.25, 0.5, 0.75, 0.95),
                                 seed = 123) {
  set.seed(seed)
  theta <- fit$theta
  Sgrid <- fit$anchored$Sgrid
  phi_grid <- fit$anchored$varphi_theta
  sigma <- fit$parameters$sigma
  time <- fit$time

  out <- vector("list", length(time))
  for (k in seq_along(time)) {
    w <- normalize_rows_prob(matrix(phi_grid[k, ], nrow = 1L))[1L, ]
    idx <- sample.int(length(theta), size = B, replace = TRUE, prob = w)
    mu <- Sgrid[k, idx]
    ysim <- mu + rnorm(B, mean = 0, sd = sigma)
    qs <- stats::quantile(ysim, probs = probs, names = FALSE)
    out[[k]] <- c(time = time[k], mean = mean(ysim), setNames(qs, paste0("q", probs)))
  }
  as.data.frame(do.call(rbind, out))
}
```

## Practical KPT algorithms

```{r kpt-algorithms}
# The practical estimator implemented here mirrors the revised paper:
# - phase densities live on the discrete simplex with row mean = 1
# - KPT-GEM uses projected phase updates with a fixed initialized surface
# - KPT-FFT alternates projected phase updates and surface updates
# - convergence is tracked using real differences, not overwritten objects
# - any optional log-likelihood trace is explicitly labeled as observed log-likelihood

kpt_fit <- 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 = NULL,
                    max_iter = 30, tol = 1e-3,
                    sigma = NULL,
                    init_surface = c("empirical", "constant"),
                    track_loglik = FALSE,
                    verbose = TRUE) {
  algorithm <- match.arg(algorithm)
  init_surface <- match.arg(init_surface)
  if (is.null(update_surface)) update_surface <- identical(algorithm, "FFT")

  I <- nrow(Y)
  K <- ncol(Y)
  if (is.null(time)) time <- seq_len(K)

  theta <- theta_grid(L_theta)
  nvec <- make_nvec(J)
  Nn <- length(nvec)
  omega <- 2*pi * (0:(Nn - 1L)) / Nn
  Lambda_phi <- (2 * sin(omega / 2))^(2 * p_phi)
  Lambda_F <- (2 * sin(omega / 2))^(2 * p_F)

  if (is.null(sigma)) {
    sigma <- 0.5 * stats::mad(as.vector(Y), center = 0) + 1e-6
  }

  phi_grid <- matrix(1, nrow = K, ncol = L_theta)
  Sgrid <- initialize_surface(Y, time, theta, method = init_surface)

  f_n <- matrix(0 + 0i, nrow = K, ncol = Nn)
  for (k in seq_len(K)) {
    f_n[k, ] <- enforce_reality_coeffs(discrete_fourier_coeffs(Sgrid[k, ], theta, nvec), nvec)
  }

  delta_phi <- delta_f <- numeric(max_iter)
  obs_loglik <- rep(NA_real_, max_iter)

  for (it in seq_len(max_iter)) {
    phi_old <- phi_grid
    f_old <- f_n

    # E-step: posterior weights and mixed moments
    m_hat <- matrix(0 + 0i, nrow = K, ncol = Nn)
    for (k in seq_len(K)) {
      prior_k <- normalize_rows_prob(matrix(phi_grid[k, ], nrow = 1L))[1L, ]
      loglik_mat <- outer(Y[, k], Sgrid[k, ], function(y, mu) -(y - mu)^2 / (2 * sigma^2))
      logw <- sweep(loglik_mat, 2L, log(pmax(prior_k, 1e-15)), "+")
      W <- softmax_rows(logw)
      Eneg <- exp(-1i * outer(theta, nvec))
      xi <- W %*% Eneg
      m_hat[k, ] <- colMeans(Y[, k] * xi)
    }

    # Phase update: regularized deconvolution, then exact Euclidean projection
    phi_hat_n <- matrix(0 + 0i, nrow = K, ncol = Nn)
    phi_grid_new <- matrix(0, nrow = K, ncol = L_theta)
    for (k in seq_len(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)
      phi_grid_new[k, ] <- synthesize_from_coeffs(phi_hat_n[k, ], theta, nvec)
    }
    phi_grid_new <- row_project_simplex(phi_grid_new)

    # Recompute phase coefficients from the projected densities
    phi_n_projected <- matrix(0 + 0i, nrow = K, ncol = Nn)
    for (k in seq_len(K)) {
      phi_n_projected[k, ] <- discrete_fourier_coeffs(phi_grid_new[k, ], theta, nvec)
    }

    # Optional surface update
    if (isTRUE(update_surface)) {
      f_n_new <- matrix(0 + 0i, nrow = K, ncol = Nn)
      Sgrid_new <- matrix(0, nrow = K, ncol = L_theta)
      for (k in seq_len(K)) {
        Phi_w <- coeffs_to_Xomega_fft(phi_n_projected[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)
        Sgrid_new[k, ] <- synthesize_from_coeffs(f_n_new[k, ], theta, nvec)
      }
      f_n <- f_n_new
      Sgrid <- Sgrid_new
    }

    phi_grid <- phi_grid_new

    delta_phi[it] <- sqrt(mean((phi_grid - phi_old)^2))
    delta_f[it] <- sqrt(mean(Mod(f_n - f_old)^2))

    if (isTRUE(track_loglik)) {
      obs_loglik[it] <- observed_loglik_gaussian(Y, Sgrid, phi_grid, sigma)
    }

    if (isTRUE(verbose)) {
      if (isTRUE(update_surface)) {
        message(sprintf("[%s] iter %02d: dphi = %.3e, df = %.3e", algorithm, it, delta_phi[it], delta_f[it]))
      } else {
        message(sprintf("[%s] iter %02d: dphi = %.3e", algorithm, it, delta_phi[it]))
      }
    }

    if (max(delta_phi[it], delta_f[it]) < tol) break
  }

  fit <- list(
    varphi_theta = phi_grid,
    Sgrid = Sgrid,
    phi_n = phi_n_projected,
    f_n = f_n,
    theta = theta,
    time = time,
    convergence = list(
      iterations = it,
      delta_phi = delta_phi[seq_len(it)],
      delta_f = delta_f[seq_len(it)],
      observed_loglik = obs_loglik[seq_len(it)]
    ),
    parameters = list(
      algorithm = algorithm,
      J = J,
      L_theta = L_theta,
      lambda_phi = lambda_phi,
      p_phi = p_phi,
      lambda_F = lambda_F,
      p_F = p_F,
      sigma = sigma,
      update_surface = update_surface,
      max_iter = max_iter,
      tol = tol,
      normalization = "row mean = 1"
    )
  )

  anchored <- anchor_by_first_moment(fit$varphi_theta, fit$theta)
  fit$anchored <- list(
    varphi_theta = anchored$phi,
    Sgrid = apply_anchor_to_surface(fit$Sgrid, anchored$delta),
    delta = anchored$delta
  )
  fit
}

fit_kpt_gem <- function(Y, time, ...) {
  kpt_fit(
    Y = Y,
    time = time,
    algorithm = "GEM",
    update_surface = FALSE,
    ...
  )
}

fit_kpt_fft <- function(Y, time, ...) {
  kpt_fit(
    Y = Y,
    time = time,
    algorithm = "FFT",
    update_surface = TRUE,
    ...
  )
}
```

## Plotting helpers

```{r plotting-helpers}
plot_phase_slice <- function(theta, truth, est_list, title = "", subtitle = "") {
  df <- tibble(theta = theta, Truth = truth)
  for (nm in names(est_list)) df[[nm]] <- est_list[[nm]]
  df_long <- tidyr::pivot_longer(df, -theta, names_to = "Series", values_to = "Density")
  ggplot(df_long, aes(theta, Density, color = Series)) +
    geom_line(linewidth = 1) +
    labs(title = title, subtitle = subtitle, x = expression(theta), y = expression(varphi(theta))) +
    theme_minimal(base_size = 12)
}

plot_surface_heatmap <- function(Sgrid, time, theta, title = "") {
  df <- expand.grid(time = time, theta = theta)
  df$value <- as.vector(Sgrid)
  ggplot(df, aes(theta, time, fill = value)) +
    geom_raster(interpolate = FALSE) +
    scale_y_continuous(expand = c(0, 0)) +
    labs(title = title, x = expression(theta), y = "time") +
    theme_minimal(base_size = 12)
}

plot_phase_metric_time <- function(metric_df, title = "") {
  df <- metric_df |>
    tidyr::pivot_longer(-time, names_to = "Metric", values_to = "Value")
  ggplot(df, aes(time, Value, color = Metric)) +
    geom_line(linewidth = 0.9) +
    labs(title = title, x = "time", y = "metric") +
    theme_minimal(base_size = 12)
}

plot_predictive_overlay <- function(Y, time, pred_list, sample_frac = 0.2) {
  n <- nrow(Y)
  K <- ncol(Y)
  obs_df <- expand.grid(rep = seq_len(n), time = time)
  obs_df$y <- as.vector(Y)
  if (sample_frac < 1) {
    keep <- sample.int(nrow(obs_df), size = max(1L, floor(sample_frac * nrow(obs_df))))
    obs_df <- obs_df[keep, ]
  }

  pred_df <- purrr::imap_dfr(pred_list, function(df, nm) {
    tibble(
      time = df$time,
      mean = df$mean,
      q05 = df$q0.05,
      q25 = df$q0.25,
      q75 = df$q0.75,
      q95 = df$q0.95,
      Method = nm
    )
  })

  ggplot() +
    geom_point(data = obs_df, aes(time, y), alpha = 0.18, shape = 1) +
    geom_ribbon(data = pred_df, aes(time, ymin = q05, ymax = q95, fill = Method), alpha = 0.15) +
    geom_ribbon(data = pred_df, aes(time, ymin = q25, ymax = q75, fill = Method), alpha = 0.25) +
    geom_line(data = pred_df, aes(time, mean, color = Method), linewidth = 1) +
    labs(title = "Posterior predictive overlays", x = "time", y = "response") +
    theme_minimal(base_size = 12)
}
```

# Statistical Inference {#sec-inference}

The revised paper recommends bootstrap-based uncertainty quantification whenever the phase projection or weak spectral identifiability may materially affect finite-sample inference. The following helper implements a **row bootstrap**: entire replicate trajectories are resampled with replacement, the estimator is re-fitted, and the anchored phase laws are collected.

This is computationally expensive and is therefore provided as an optional protocol component. It is not run by default in the full experiments below.

```{r bootstrap-helpers}
bootstrap_phase_bands <- function(Y, time, fit_fun, B = 25, seed = 123, probs = c(0.025, 0.5, 0.975), ...) {
  set.seed(seed)
  N <- nrow(Y)
  fits <- vector("list", B)
  for (b in seq_len(B)) {
    idx <- sample.int(N, size = N, replace = TRUE)
    fits[[b]] <- fit_fun(Y[idx, , drop = FALSE], time = time, ...)
  }

  phi_array <- simplify2array(lapply(fits, function(f) f$anchored$varphi_theta))
  # [K, L, B]
  bands <- apply(phi_array, c(1, 2), quantile, probs = probs, na.rm = TRUE)
  list(fits = fits, bands = bands, probs = probs)
}
```

Example usage for local uncertainty studies:

```{r bootstrap-example, eval=FALSE}
# Example (not run by default because it can be slow):
# ds_boot <- bootstrap_phase_bands(
#   Y = ds$Y,
#   time = ds$time,
#   fit_fun = fit_kpt_gem,
#   B = 50,
#   J = 16,
#   L_theta = 256,
#   lambda_phi = 5e-3,
#   p_phi = 1,
#   lambda_F = 5e-3,
#   p_F = 1,
#   sigma = 0.15,
#   max_iter = 30,
#   tol = 1e-3,
#   verbose = FALSE
# )
```

# Applications {#sec-applications}

## Simulation study I: phase-modulated double-slit experiment

### Double-slit simulator

```{r ds-simulator}
dvonmises <- function(theta, mu, kappa) {
  exp(kappa * cos(theta - mu)) / (2*pi * besselI(kappa, 0))
}

simulate_double_slit <- function(N = 100, K = 50, L_theta = 256,
                                 kappa_phase = 4, noise_sd = 0.15,
                                 kappa_wave = 12, seed = 2025) {
  set.seed(seed)
  theta <- theta_grid(L_theta)
  time <- seq(0, 1, length.out = K)

  S_true <- matrix(0, nrow = K, ncol = L_theta)
  phi_true <- matrix(0, nrow = K, ncol = L_theta)

  for (k in seq_len(K)) {
    t <- time[k]
    envelope <- 0.8 * exp(-20 * (t - 0.5)^2)
    baseline <- 0.2 * sin(2*pi * t)
    S_true[k, ] <- baseline +
      envelope * cos(kappa_wave * t + theta) +
      0.24 * envelope * cos(2 * (kappa_wave * t + theta) + pi / 6)

    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)
    p1 <- dvonmises(theta, mu1, kappa_phase)
    p2 <- dvonmises(theta, mu2, kappa_phase)
    dens <- w * p1 + (1 - w) * p2
    phi_true[k, ] <- dens / mean(dens)
  }

  Y <- matrix(0, nrow = N, ncol = K)
  for (k in seq_len(K)) {
    probs <- normalize_rows_prob(matrix(phi_true[k, ], nrow = 1L))[1L, ]
    idx <- sample.int(L_theta, size = N, replace = TRUE, prob = probs)
    Y[, k] <- S_true[k, idx] + rnorm(N, 0, noise_sd)
  }

  truth_anchor <- anchor_by_first_moment(phi_true, theta)
  list(
    Y = Y,
    time = time,
    theta = theta,
    phi_true = phi_true,
    S_true = S_true,
    phi_true_anchored = truth_anchor$phi,
    S_true_anchored = apply_anchor_to_surface(S_true, truth_anchor$delta),
    params = list(
      N = N, K = K, L_theta = L_theta,
      kappa_phase = kappa_phase, noise_sd = noise_sd,
      kappa_wave = kappa_wave, seed = seed
    )
  )
}

ds <- simulate_double_slit()
```

### Fit KPT-GEM and KPT-FFT to the double-slit data

```{r ds-fits}
fit_ds_gem <- fit_kpt_gem(
  Y = 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,
  sigma = 0.15,
  max_iter = 30,
  tol = 1e-3,
  track_loglik = FALSE,
  verbose = TRUE
)

fit_ds_fft <- fit_kpt_fft(
  Y = 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,
  sigma = 0.15,
  max_iter = 30,
  tol = 1e-3,
  track_loglik = FALSE,
  verbose = TRUE
)
```

### Double-slit evaluation after alignment

```{r ds-metrics}
phase_ds_gem <- phase_metrics_df(ds$phi_true_anchored, fit_ds_gem$anchored$varphi_theta, ds$time)
phase_ds_fft <- phase_metrics_df(ds$phi_true_anchored, fit_ds_fft$anchored$varphi_theta, ds$time)

surf_ds_gem <- surface_L2_cone(ds$S_true_anchored, fit_ds_gem$anchored$Sgrid, ds$time)
surf_ds_fft <- surface_L2_cone(ds$S_true_anchored, fit_ds_fft$anchored$Sgrid, ds$time)

summary_ds <- tibble(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(phase_ds_gem$JSD_bits), mean(phase_ds_fft$JSD_bits)),
  Mean_Hellinger = c(mean(phase_ds_gem$Hellinger), mean(phase_ds_fft$Hellinger)),
  Mean_TV = c(mean(phase_ds_gem$TV), mean(phase_ds_fft$TV)),
  Rel_L2_surface = c(surf_ds_gem$rel_L2, surf_ds_fft$rel_L2)
)

knitr::kable(summary_ds, digits = 4, caption = "Double-slit experiment: aligned reconstruction metrics.")
```

### Double-slit representative phase laws and surfaces

```{r ds-figures}
rep_idx_ds <- c(1, 13, 25, 37, 50)
for (k in rep_idx_ds) {
  print(
    plot_phase_slice(
      theta = ds$theta,
      truth = ds$phi_true_anchored[k, ],
      est_list = list(
        `KPT-GEM` = fit_ds_gem$anchored$varphi_theta[k, ],
        `KPT-FFT` = fit_ds_fft$anchored$varphi_theta[k, ]
      ),
      title = "Double-slit phase recovery",
      subtitle = sprintf("time = %.3f", ds$time[k])
    )
  )
}

plot_surface_heatmap(ds$S_true_anchored, ds$time, ds$theta, title = "Double-slit truth: anchored kime surface")
plot_surface_heatmap(fit_ds_gem$anchored$Sgrid, ds$time, ds$theta, title = "Double-slit KPT-GEM: anchored estimated surface")
plot_surface_heatmap(fit_ds_fft$anchored$Sgrid, ds$time, ds$theta, title = "Double-slit KPT-FFT: anchored estimated surface")
```

### Double-slit phase metrics over time

```{r ds-phase-metrics-time}
plot_phase_metric_time(phase_ds_gem, title = "Double-slit KPT-GEM phase metrics versus time")
plot_phase_metric_time(phase_ds_fft, title = "Double-slit KPT-FFT phase metrics versus time")
```

### Posterior predictive time overlays

```{r ds-predictive-overlays}
pred_ds_gem <- kpt_predictive_bands(fit_ds_gem, B = 3000, seed = 11)
pred_ds_fft <- kpt_predictive_bands(fit_ds_fft, B = 3000, seed = 12)
plot_predictive_overlay(ds$Y, ds$time, list(`KPT-GEM` = pred_ds_gem, `KPT-FFT` = pred_ds_fft), sample_frac = 0.25)
```

### Predictive mapping from the kime representation to screen position

The revised paper keeps the double-slit screen-position mapping and now treats it as a **proper normalized likelihood**. The implementation below follows the paper directly.

```{r ds-screen-model}
ds_physics <- function(lambda = 532e-9,
                       d = 20e-6,
                       L = 0.5,
                       x_window = c(-0.02, 0.02),
                       nx = 800,
                       envelope_sigma = 0.005) {
  x <- seq(x_window[1], x_window[2], length.out = nx)
  envelope <- exp(-0.5 * (x / envelope_sigma)^2)
  alpha <- 2*pi * d / (lambda * L)
  list(
    lambda = lambda,
    d = d,
    L = L,
    x = x,
    envelope = envelope,
    alpha = alpha
  )
}

rectify_surface <- function(u) {
  pmax(u - min(u), 0)
}

screen_density_from_kpt <- function(Sbar, phibar, theta, physics, mu = 0, rectify = rectify_surface) {
  x <- physics$x
  alpha <- physics$alpha
  env <- physics$envelope
  phibar <- phibar / mean(phibar)

  intensity <- numeric(length(x))
  for (i in seq_along(x)) {
    theta_query <- (alpha * x[i] + mu - theta) %% (2*pi)
    sval <- periodic_interp(theta, Sbar, theta_query)
    intensity[i] <- env[i] * mean(rectify(sval) * phibar)
  }
  intensity <- pmax(intensity, 0)
  density <- intensity / max(trapz_vec(x, intensity), 1e-15)
  tibble(x = x, intensity = intensity, density = density)
}

simulate_screen_hits <- function(Sbar, phibar, theta, physics, n_hits = 1000, mu = 0, seed = 1) {
  set.seed(seed)
  dens_df <- screen_density_from_kpt(Sbar, phibar, theta, physics, mu = mu)
  probs <- dens_df$density / sum(dens_df$density)
  sample(dens_df$x, size = n_hits, replace = TRUE, prob = probs)
}

posterior_mu_discrete <- function(x_hits, Sbar, phibar, theta, physics, mu_grid = theta) {
  prior <- pmax(phibar, 1e-15)
  prior <- prior / sum(prior)

  logpost <- numeric(length(mu_grid))
  for (m in seq_along(mu_grid)) {
    dens_df <- screen_density_from_kpt(Sbar, phibar, theta, physics, mu = mu_grid[m])
    px <- approx(dens_df$x, dens_df$density, xout = x_hits, rule = 2L)$y
    logpost[m] <- log(prior[m]) + sum(log(pmax(px, 1e-300)))
  }
  logpost <- logpost - max(logpost)
  post <- exp(logpost)
  post <- post / sum(post)
  tibble(mu = mu_grid, posterior = post)
}

plot_screen_density <- function(hit_df, density_list) {
  dens_df <- purrr::imap_dfr(density_list, function(df, nm) dplyr::mutate(df, Method = nm))
  ggplot() +
    geom_histogram(data = hit_df, aes(x, y = after_stat(density)), bins = 60,
                   fill = "grey80", color = "grey50", alpha = 0.5) +
    geom_line(data = dens_df, aes(x, density, color = Method), linewidth = 1) +
    labs(title = "Double-slit screen-position density", x = "screen position x", y = "density") +
    theme_minimal(base_size = 12)
}
```

```{r ds-screen-analysis}
physics_ds <- ds_physics()

Sbar_true <- colMeans(ds$S_true_anchored)
phibar_true <- colMeans(ds$phi_true_anchored)
Sbar_gem  <- colMeans(fit_ds_gem$anchored$Sgrid)
phibar_gem <- colMeans(fit_ds_gem$anchored$varphi_theta)
Sbar_fft  <- colMeans(fit_ds_fft$anchored$Sgrid)
phibar_fft <- colMeans(fit_ds_fft$anchored$varphi_theta)

screen_true <- screen_density_from_kpt(Sbar_true, phibar_true, ds$theta, physics_ds, mu = 0)
screen_gem  <- screen_density_from_kpt(Sbar_gem, phibar_gem, ds$theta, physics_ds, mu = 0)
screen_fft  <- screen_density_from_kpt(Sbar_fft, phibar_fft, ds$theta, physics_ds, mu = 0)

x_hits <- simulate_screen_hits(Sbar_true, phibar_true, ds$theta, physics_ds, n_hits = 1500, mu = 0, seed = 99)
hit_df <- tibble(x = x_hits)

plot_screen_density(hit_df, list(Truth = screen_true, `KPT-GEM` = screen_gem, `KPT-FFT` = screen_fft))
```

```{r ds-mu-posterior}
mu_post_gem <- posterior_mu_discrete(x_hits, Sbar_gem, phibar_gem, ds$theta, physics_ds)
mu_post_fft <- posterior_mu_discrete(x_hits, Sbar_fft, phibar_fft, ds$theta, physics_ds)

mu_df <- bind_rows(
  mutate(mu_post_gem, Method = "KPT-GEM"),
  mutate(mu_post_fft, Method = "KPT-FFT")
)

ggplot(mu_df, aes(mu, posterior, color = Method)) +
  geom_line(linewidth = 1) +
  labs(title = "Posterior distribution for the source offset mu", x = expression(mu), y = "posterior probability") +
  theme_minimal(base_size = 12)
```

## Simulation study II: phase-dependent event-related fMRI

### fMRI simulator matched to the revised paper

```{r fmri-simulator}
canonical_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
}

simulate_fmri_phase_dependent <- 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)
  K <- as.integer(duration / TR)
  time <- seq(0, by = TR, length.out = K)
  theta <- theta_grid(L_theta)

  block_len <- 30 / TR
  design <- rep(c(0, 1), each = block_len)
  design <- rep_len(design, K)

  hrf_vals <- canonical_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
  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)
  }

  baseline <- 0.05 * sin(2*pi * time / max(time))
  mu_theta <- 0.2 * sin(2*pi * time / max(time))

  S_true <- matrix(0, nrow = K, ncol = L_theta)
  phi_true <- matrix(0, nrow = K, ncol = L_theta)

  for (k in seq_len(K)) {
    S_true[k, ] <- baseline[k] + bold[k] * cos(theta) + 0.25 * bold[k] * cos(2 * theta)
    kappa <- if (design[k] == 1) kappa_on else kappa_off
    dens <- dvonmises(theta, mu_theta[k], kappa)
    phi_true[k, ] <- dens / mean(dens)
  }

  I <- P * R
  Y <- matrix(0, nrow = I, ncol = K)
  for (k in seq_len(K)) {
    probs <- normalize_rows_prob(matrix(phi_true[k, ], nrow = 1L))[1L, ]
    idx <- sample.int(L_theta, size = I, replace = TRUE, prob = probs)
    Y[, k] <- S_true[k, idx] + rnorm(I, 0, noise_sd)
  }

  truth_anchor <- anchor_by_first_moment(phi_true, theta)
  list(
    Y = Y,
    time = time,
    theta = theta,
    design = design,
    bold = bold,
    baseline = baseline,
    phi_true = phi_true,
    S_true = S_true,
    phi_true_anchored = truth_anchor$phi,
    S_true_anchored = apply_anchor_to_surface(S_true, truth_anchor$delta),
    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,
      seed = seed
    )
  )
}

fm <- simulate_fmri_phase_dependent()
```

### Fit KPT-GEM and KPT-FFT to the fMRI data

```{r fmri-fits}
fit_fm_gem <- fit_kpt_gem(
  Y = 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,
  sigma = 0.2,
  max_iter = 25,
  tol = 8e-4,
  track_loglik = FALSE,
  verbose = TRUE
)

fit_fm_fft <- fit_kpt_fft(
  Y = 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,
  sigma = 0.2,
  max_iter = 25,
  tol = 8e-4,
  track_loglik = FALSE,
  verbose = TRUE
)
```

### fMRI evaluation after alignment

```{r fmri-metrics}
phase_fm_gem <- phase_metrics_df(fm$phi_true_anchored, fit_fm_gem$anchored$varphi_theta, fm$time)
phase_fm_fft <- phase_metrics_df(fm$phi_true_anchored, fit_fm_fft$anchored$varphi_theta, fm$time)

surf_fm_gem <- surface_L2_cone(fm$S_true_anchored, fit_fm_gem$anchored$Sgrid, fm$time)
surf_fm_fft <- surface_L2_cone(fm$S_true_anchored, fit_fm_fft$anchored$Sgrid, fm$time)

summary_fm <- tibble(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(phase_fm_gem$JSD_bits), mean(phase_fm_fft$JSD_bits)),
  Mean_Hellinger = c(mean(phase_fm_gem$Hellinger), mean(phase_fm_fft$Hellinger)),
  Mean_TV = c(mean(phase_fm_gem$TV), mean(phase_fm_fft$TV)),
  Rel_L2_surface = c(surf_fm_gem$rel_L2, surf_fm_fft$rel_L2)
)

knitr::kable(summary_fm, digits = 4, caption = "Phase-dependent fMRI experiment: aligned reconstruction metrics.")
```

### fMRI raw traces, phase laws, and surfaces

```{r fmri-figures}
raw_df <- tibble(
  time = rep(fm$time, each = min(12, nrow(fm$Y))),
  value = as.vector(t(fm$Y[seq_len(min(12, nrow(fm$Y))), ])),
  run = factor(rep(seq_len(min(12, nrow(fm$Y))), times = length(fm$time)))
)

ggplot(raw_df, aes(time, value, group = run, color = run)) +
  geom_line(alpha = 0.8) +
  labs(title = "Synthetic fMRI raw time series (subset of runs)", x = "time (s)", y = "signal") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

plot_phase_metric_time(phase_fm_gem, title = "fMRI KPT-GEM phase metrics versus time")
plot_phase_metric_time(phase_fm_fft, title = "fMRI KPT-FFT phase metrics versus time")

k_rep <- which.min(abs(fm$time - 28))
plot_phase_slice(
  theta = fm$theta,
  truth = fm$phi_true_anchored[k_rep, ],
  est_list = list(
    `KPT-GEM` = fit_fm_gem$anchored$varphi_theta[k_rep, ],
    `KPT-FFT` = fit_fm_fft$anchored$varphi_theta[k_rep, ]
  ),
  title = "fMRI phase recovery at a representative time point",
  subtitle = sprintf("time = %.1f seconds", fm$time[k_rep])
)

plot_surface_heatmap(fm$S_true_anchored, fm$time, fm$theta, title = "fMRI truth: anchored kime surface")
plot_surface_heatmap(fit_fm_gem$anchored$Sgrid, fm$time, fm$theta, title = "fMRI KPT-GEM: anchored estimated surface")
plot_surface_heatmap(fit_fm_fft$anchored$Sgrid, fm$time, fm$theta, title = "fMRI KPT-FFT: anchored estimated surface")
```

### fMRI posterior predictive overlays

```{r fmri-predictive-overlays}
pred_fm_gem <- kpt_predictive_bands(fit_fm_gem, B = 3000, seed = 21)
pred_fm_fft <- kpt_predictive_bands(fit_fm_fft, B = 3000, seed = 22)
plot_predictive_overlay(fm$Y, fm$time, list(`KPT-GEM` = pred_fm_gem, `KPT-FFT` = pred_fm_fft), sample_frac = 0.10)
```

# Conclusions {#sec-conclusions}

This revised protocol resolves the main issues identified in the previous notebook version.

1. It uses the **time-specific latent phase model** everywhere.
2. It enforces the correct normalization convention based on $d\theta/(2\pi)$ and the discrete row-mean-equals-one rule.
3. It distinguishes the **exact EM reference algorithm** from the **practical projected routines** used numerically.
4. It fixes the earlier convergence-diagnostic bug by computing differences before overwriting state variables.
5. It replaces the earlier truncation-and-renormalization rule with the **exact Euclidean projection** onto the discrete simplex.
6. It stores and uses `sigma` consistently through the predictive functions.
7. It implements the **phase-dependent fMRI simulation** used in the revised paper.
8. It uses one synchronized simulation-parameter table matching the manuscript.
9. It keeps the double-slit screen-position mapping and turns it into a **proper normalized likelihood**, including a posterior for the source offset $\mu$.
10. It removes unrelated experimental debris and external dependencies so the protocol is self-contained and locally executable.

For peer-review submission, this notebook should now serve as a clean computational companion to the revised paper. The next iteration could add larger-scale bootstrap studies, repeated-seed Monte Carlo summaries, or sensitivity analyses for the regularization parameters.

```{r session-info}
sessionInfo()
```
