--- title: "SOCR TCIU: Chapter 3: The Kime-Phase Problem" 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 word_document: toc: true toc_depth: '3' subtitle: "Kime-phase Estimation and Kime-Operators in Spacekime Representation" editor_options: markdown: wrap: 72 --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # The effects of Kime-Magnitudes and Kime-Phases Jointly, the *amplitude spectrum* (magnitudes) and the *phase spectrum* (phases) uniquely describe the spacetime representation of a signal. However, the importance of each of these two spectra is not equivalent. In general, the effect of the phase spectrum is more important compared to the corresponding effects of the amplitude spectrum. In other words, the magnitudes are less susceptible to noise or the accuracy of their estimations. The effects of magnitude perturbations are less critical relative to proportional changes in the phase spectrum. For instance, particularly in terms of spacetime locations where the signal is zero, the signal can be reconstructed (by the IFT) relatively accurately using incorrect magnitudes solely by using the correct phases [REF](https://doi.org/10.1533/9780857099457.1.75). For a real valued signal $f$, suppose the amplitude of its Fourier transform, $FT(f)=\hat{f}$, is $A(\omega) > 0, \forall \omega$, then: $$f(x)=IFT(\hat{f})=Re\left (\frac{1}{2\pi}\int_{R} \underbrace{A(\omega)e^{i\phi(\omega)}}_{\hat{f}(\omega)}\ e^{i\omega x}d\omega \right)= Re\left (\frac{1}{2\pi}\int_{R}A(\omega)e^{i(\phi(\omega)+\omega x)}d\omega\right) = \frac{1}{2\pi}\int_{R} {A(\omega) \cos(\phi(\omega)+\omega x)}d\omega.$$ Thus, the zeros of $f(x)$ occur for $\omega x+ \phi(\omega)=\pm k\frac{\pi}{2}$, $k= 1,2,3,.$. A solely amplitude driven reconstruction $\left ( f_A(x)=IFT(\hat{f})=\frac{1}{2\pi}\int_{R}\underbrace{A(\omega)}_{no\ phase}\ e^{i\omega x}d\omega \right)$ would yield worse results than a solely-phase based reconstruction $\left ( f_{\phi}(x)=IFT(\hat{f})=\frac{1}{2\pi} \int_{R}\underbrace{e^{i\phi(\omega)}}_{no\ amplitude}\ e^{i\omega x}d\omega\right )$. The latter would have a different total energy from the original signal, however, it would include some signal recognizable features as the zeroth-level curves of the original $f$ and the phase-only reconstruction $f_{\phi}$ signals will be preserved. This suggests that the *Fourier phase* of a signal is more informative than the *Fourier amplitude*, i.e., the magnitudes are robust to errors or perturbations. In X-ray crystallography, crystal structures are bombarded by particles/waves, which are diffracted by the crystal to yield the observed diffraction spots or patterns. Each diffraction spot corresponds to a point in the reciprocal lattice and represents a particle wave with some specific amplitude and a relative phase. Probabilistically, as the particles (e.g., gamma-rays or photons) are reflected from the crystal, their scatter directions are proportional to the square of the wave amplitude, i.e., the square of the wave Fourier magnitude. X-rays capture these amplitudes as counts of particle directions, but miss all information about the relative phases of different diffraction patterns. Spacekime analytics are analogous to X-ray crystallography, DNA helix modeling, and other applications, where only the Fourier magnitudes (time), i.e., power spectrum, is only observed, but not the phases (kime-directions), which need to be estimated to correctly reconstruct the intrinsic 3D object structure [REF](https://www.sciencedirect.com/science/article/pii/B9781904275268500051), in our case, the correct spacekime analytical inference. Clearly, signal reconstruction based solely on either the amplitudes or the phases is an ill-posed problem, i.e., there will be many *alternative solutions*. In practice, such *signal* or *inference* reconstructions are always application-specific, rely on some a priori knowledge on the process (or objective function), or depend an information-theoretic criteria to derive conditional solutions. Frequently, such solutions are obtained via least squares, maximum entropy criteria, maximum a posterior distributions, Bayesian estimations, or simply by approximating the unknown amplitudes or phases using prior observations, similar processes, or theoretical models. # Solving the Missing Kime-Phase Problem There are many alternative solutions to the problem of estimating the unobserved kime-phases. All solutions depend on the quality of the data (e.g., noise), the signal energy (e.g., strength of association between covariates and outcomes), and the general experimental design. There can be rather large errors in the phase reconstructions, which will in turn affect the final spacekime analytic results. Most phase-problem solutions are based on the idea that having some *prior knowledge* about the characteristics of the experimental design (case-study phenomenon) and the desired inference (spacekime analytics). For instance, if we artificially *load the energy* of the case-study (e.g., by lowering the noise, increasing the SNR, or increasing the strength of the relation between explanatory and outcome variables), the phases computed from the this stronger-signal dataset will be more accurate representations than the original phase estimates. Examples of phase-problem solutions include *energy modification* and *fitting and refinement* methods. In [TCIU Section 6 (Circular Kime-Phase and Electron Orbit Densities) we cover the basics of the kime-phase representation](https://www.socr.umich.edu/TCIU/HTMLs/Chapter6_Kime_Phases_Circular.html). ![](https://socr.umich.edu/docs/uploads/2025/SOCR_KimePhase_Animation_2025.gif) Three kime-phase distributions at a fixed spatial location. ## Energy Modification Strategies In general, *energy modification* techniques rely on prior knowledge, testable hypotheses, or intuition to modify the dataset by strengthening the *expected* relation we are trying to uncover using spacekime analytics. ### Kime-phase noise distribution flattening In many practical applications, part of the dataset (including both cases and features) include valuable information, whereas the rest of the data may include irrelevant, noisy, or disruptive information. Clearly, we can't explicitly untangle these two components, however, we do expect that the irrelevant data portion would yield uninformative/unimportant kime-phases, which may be used to estimate the kime-phase noise-level and noise-distribution. Intuitively, if we modify the dataset to flatten the irrelevant kime-phases, the estimates of the corresponding true-signal kime-phases may be more accurate or more representative. We can think of this process as using kime-phase information from some known strong features to improve the kime-phase information of other particular features. Kime-phase noise distribution flattening requires that the kime-phases be good enough to detect the boundaries between the strong-features and the rest. ### Multi-sample Kime-Phase Averaging It's natural to assume that multiple instances of the same process would yield similar analytics and inference results. For a large dataset, we can use ensemble methods (e.g., SuperLearner, and CBDA) to iteratively generate independent samples, which would be expected to lead to analogous kime-phase estimated and analytical results. Thus, we expect that when salient features are extracted by spacekime analytics based on independent samples, their kime-phase estimates should be highly associated (e.g., correlated), albeit perhaps not identical. However, weak features would exhibit exactly the opposite effect - their kime-phases may be highly variable (noisy). By averaging the kime-phases, noisy-areas in the dataset may cancel out, whereas, patches of strong-signal may preserve the kime-phase details, which would lead to increased kime forecasting accuracy and reproducibility of the kime analytics. ### Histogram equalization As common experimental designs and similar datasets exhibit analogous characteristics, the corresponding spacekime analytics are also expected to be synergistic. Spacekime inference that does not yield results in some controlled or expected range, may be indicative of incorrect kime-phase estimation. We can use histogram equalization methods to improve the kime-phase estimates. This may be accomplished by altering the distribution of kime-phases to either match the phase distribution of other similar experimental designs or generate more expected spacekime analytical results. ### Fitting and refinement Related to *energy modification* strategies, the *fitting and refinement* technique capitalizes on the fact that strong energy datasets tend to have a smaller set of salient features. So, if we construct case-studies with some strong features, the corresponding kime-phases will be more accurate, and the resulting inference/analytics will be more powerful and highly reproducible. Various classification, regression, supervised and unsupervised methods, and other model-based techniques allow us to both fit a model (estimate coefficients and structure) as well as apply the model for outcome predictions and forecasting. Such models permit control over the characteristics of individual features and multivariate inter-relations, which can be exploited to gather valuable kime-phase information. Starting with a reasonable guess (kime-phase prior), the *fitting and refinement* technique can be applied iteratively to (1) reconstructing the data into spacetime using the kime-phase estimates, (2) fit or estimate the spacekime analytical model, (3) compare the analytical results and inference to expected outcomes, and (4) refine the kime-phase estimator aiming to gain better outcomes (#3). Indeed, other *energy modification* strategies (e.g., averaging or flattening) can be applied before a new iteration to build a new model is initiated (#1 and #2). ## Data Source Type ```{r} library(EBImage) square_arr <- matrix(nrow=256, ncol=256) circle_arr <- matrix(nrow=256, ncol=256) for (i in 1:256) { for (j in 1:256) { if ( abs(i-128) < 30 && abs(j-128) < 30) square_arr[i,j]=1 # sqrt((i-128)^2+(j-128)^2)/30 else square_arr[i,j]=0 if ( sqrt((i-128)^2 + (j-128)^2)<30) circle_arr[i,j]=1 # 1-sqrt((i-128)^2+(j-128)^2)/30 else circle_arr[i,j]=0 } } # FFT SHIFT #' This function is useful for visualizing the Fourier transform with the zero-frequency #' component in the middle of the spectrum. #' #' @param img_ff A Fourier transform of a 1D signal, 2D image, or 3D volume. #' @param dim Number of dimensions (-1, 1, 2, 3). #' @return A properly shifted FT of the array. #' fftshift <- function(img_ff, dim = -1) { rows <- dim(img_ff)[1] cols <- dim(img_ff)[2] # planes <- dim(img_ff)[3] swap_up_down <- function(img_ff) { rows_half <- ceiling(rows/2) return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)])) } swap_left_right <- function(img_ff) { cols_half <- ceiling(cols/2) return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half])) } #swap_side2side <- function(img_ff) { # planes_half <- ceiling(planes/2) # return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half])) #} if (dim == -1) { img_ff <- swap_up_down(img_ff) return(swap_left_right(img_ff)) } else if (dim == 1) { return(swap_up_down(img_ff)) } else if (dim == 2) { return(swap_left_right(img_ff)) } else if (dim == 3) { # Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays # install.packages("abind") library(abind) planes <- dim(img_ff)[3] rows_half <- ceiling(rows/2) cols_half <- ceiling(cols/2) planes_half <- ceiling(planes/2) img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)], img_ff[(1:rows_half), (1:cols), (1:planes)], along=1) img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)], img_ff[1:rows, 1:cols_half, (1:planes)], along=2) img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half], along=3) return(img_ff) } else { stop("Invalid dimension parameter") } } ``` ## Figure 3.8 ```{r} # image(square_arr); image(circle_arr) # display(circle_arr, method = "raster") display(square_arr, method = "raster") ``` ```{r} X1 = fft(square_arr) X1_mag <- sqrt(Re(X1)^2+Im(X1)^2) X1_phase <- atan2(Im(X1), Re(X1)) # FT of Circle # No shift applied here (perhaps should be consistent or just show the difference?) X2 = fft(circle_arr) # display(Re(X2), method = "raster") X2_mag <- sqrt(Re(X2)^2+Im(X2)^2) # display(X2_mag, method = "raster") # magnitude only X2_phase <- atan2(Im(X2), Re(X2)) # display(X2_phase, method = "raster") # phase only display(fftshift(X1_mag), method = "raster") ``` ```{r} # Take 2: IFT Magnitude= Square and Phase = Square + IID noise (N(0,3/2)) set.seed(1234) IID_noise <- matrix(rnorm(prod(dim(X1_phase)), mean=0, sd=1.5), nrow=dim(X1_phase)[1]) # dim(IID_noise) # 256 256 Real = X1_mag * cos(X1_phase + IID_noise) Imaginary = X1_mag * sin(X1_phase + IID_noise) ift_X1mag_X1phase_Noise = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1)) plot(density(IID_noise), xlim=c(-8,8), col="blue", lwd=2) lines(density(X1_phase), col="red", lwd=2) ``` ```{r} mixed_density <- density(X1_phase + IID_noise) mixed_density_mod2pi <- density((X1_phase + IID_noise)%%(2*pi) -pi) display(ift_X1mag_X1phase_Noise, method = "raster") ``` ```{r} plot(mixed_density, xlim=c(-8,8), ylim=c(0,0.17), col="blue", lwd=2, # should density be modulo %%(2*pi)? main="Phase Distributions: Raw Square, Square+IID N(m=0,s=3/2), Mixed Phases mod 2*Pi", cex.main=0.8, xlab = "Phase", ylab = "Density") lines(density(X1_phase), col="red", lwd=4) lines(mixed_density_mod2pi, col="green", lwd=2) text(x=3.2, y=-0.005, expression(pi)) text(x=-3.2, y=-0.005, expression(-pi)) legend("center", legend=c("(Raw) Square Phases", "Square + N(0,1.5)", expression(paste("Square + N(0,1.5) mod 2*", pi))), col=c("red","blue", "green"), lty=1, lwd=c(4,2,2), cex=1.0, y.intersp=1.0, x.intersp=1.0, title = "Phases", bty = "n") ``` ```{r} # randomly sample 10 indices to pairwise plot ind1 <- sample(dim(X1_phase)[1], 10); ind2 <- sample(dim(X1_phase)[2], 10) corr <- cor.test((X1_phase)[ind1,ind2], (X1_phase + IID_noise)[ind1,ind2], method = "pearson", conf.level = 0.99) plot((X1_phase)[ind1,ind2], (X1_phase + IID_noise)[ind1,ind2], main= sprintf("Square Image: Raw vs. IID N(m=0,s=3/2) Noise-corrupted Phases: Corr=%s CI=(%s,%s)", round(cor(as.vector(X1_phase), as.vector(X1_phase + IID_noise)), digits=3), round(corr$conf.int[1], digits=2), round(corr$conf.int[2], digits=2)), cex.main=0.8, xlab = "Raw", ylab = "Phase + N(0,1.5)") abline( lm(as.vector((X1_phase + IID_noise)[ind1,ind2]) ~ as.vector((X1_phase)[ind1,ind2])), col="red", lwd=2) ``` ```{r} # Take 2: Same level of noise on Amplitudes: # IFT Magnitude= Square + 50*IID noise (N(0,3/2)) and Phase = Square set.seed(1234) IID_noise <- matrix(rnorm(prod(dim(X1_mag)), mean=0, sd=3/2), nrow=dim(X1_mag)[1]) # dim(IID_noise) # 256 256 corr <- cor.test((X1_mag+50*IID_noise)[ind1,ind2], (X1_mag)[ind1,ind2], method = "pearson", conf.level = 0.99) plot(density(X1_mag+50*IID_noise), xlim=c(0,40), ylim=c(0,1.5), col="blue", lwd=2) lines(density(X1_mag), col="red", lwd=2) ``` ```{r, fig.height=5, fig.width=8} plot(X1_mag[ind1,ind2], abs(X1_mag+50*IID_noise)[ind1,ind2], xlim=c(0,10), ylim=c(0,200), main= sprintf("Square Image: Mag + 50*IID N(m=0,s=3/2) Noise-corrupted vs. Raw Amplitudes: Corr=%s CI=(%s,%s)", round(cor(as.vector(X1_mag), as.vector(X1_mag + 50*IID_noise)), digits=3), round(corr$conf.int[1], digits=2), round(corr$conf.int[2], digits=2)), cex.main=0.8, xlab = "Raw", ylab = "Amplitude + 50*N(0,1.5)") abline( lm(as.vector(abs(X1_mag+50*IID_noise)[ind1,ind2]) ~ as.vector((X1_mag)[ind1,ind2])), col="red", lwd=2) ``` ```{r} Real = (X1_mag+ 50*IID_noise) * cos(X1_phase) Imaginary = (X1_mag+ 50*IID_noise) * sin(X1_phase) ift_X1mag_X1phase_Noise = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1)) display(ift_X1mag_X1phase_Noise, method = "raster") ``` ```{r} # Magnitude distributions mixed_density <- density(abs(X1_mag + 50*IID_noise)) plot(mixed_density, xlim=c(0,200), col="blue", lwd=2, main="Magnitude Distributions: Raw Square, Square+50*N(m=0,s=3/2)", xlab = "Magnitude", ylab = "Density") lines(density(X1_mag), col="red", lwd=4) legend("topright", legend=c("(Raw) Square Magnitudes", "Square + 50*N(0,1.5)"), col=c("red","blue"), lty=1, lwd=c(4,2), cex=1.0, y.intersp=1.0, x.intersp=1.0, title = "Magnitudes", bty = "n") ``` ```{r} # Take 3: Linear Transform of the Phases: SquarePhase ~ CirclePhase # Take 2: IFT Magnitude= Square and Phase = LM(CirclePhase) lm_Squ_Cir <- lm(as.vector(X1_phase) ~ as.vector(X2_phase)) plot(as.vector(X2_phase), as.vector(X1_phase), col="blue", lwd=2, xlab = "Circle", ylab = "Square", main= sprintf("Linear Phase Transformation (SquarePhase ~ CirclePhase), Corr(Cir, Squ)=%s", round(cor(as.vector(X1_phase) , as.vector(X2_phase)), digits=3))) abline(lm_Squ_Cir, col="red", lwd=2) Real = X1_mag * cos(lm_Squ_Cir$coefficients[1] + lm_Squ_Cir$coefficients[2]*X2_phase) Imaginary = X1_mag * sin(lm_Squ_Cir$coefficients[1] + lm_Squ_Cir$coefficients[2]*X2_phase) ift_X1mag_X2phase_LM = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1)) display(ift_X1mag_X2phase_LM, method = "raster") ``` # QM Phase Estimation by Probing Measurements in Different Bases The fundamental rationale for observing only the *magnitude* of the quantum wavefunction, while phases remain unobservable directly, stems from the structure of quantum mechanics. Physical observables are always *real*, tied to Hermitian operators, and measurement outcomes depend on probabilities derived from the *squared magnitude of the wavefunction*, $|\psi(x)|^2$. Global or relative phases do not affect these probabilities, unless they manifest as relative phase differences in superposition states. *Global phases* cancel out in expectation values, while *relative phases* influence interference patterns, which are still inferred indirectly through probabilities. Thus, *direct phase measurement* is impossible. Instead, phases are estimated by *observing interference effects* across multiple measurement bases. ## Pauli Operators and Bases in Quantum Mechanics The relationship between *Pauli operators* and *bases* is fundamental to quantum mechanics. *Pauli operators* are matrix representations of *Pauli matrices* plus the *identity matrix* that form a basis for $2\times 2$ Hermitian matrices $$\sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\ , \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}\ , \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\ , \quad I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$ The *eigenspectrum* of each Pauli operator includes the *eigenvalues* $\lambda = \pm 1$ with their base-specific corresponding *eigenvectors.* For the computational basis, $\sigma_z$, - the *eigenvector* $|0\rangle = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ corresponds with *eigenvalue* $\lambda = +1$, - and the *eigenvector* $|1\rangle = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$ corresponds with *eigenvalue* $\lambda = -1$. Similarly, for the $X$-basis, $\sigma_x$ - $|+\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix}$ with $\lambda = +1$, - $|-\rangle = \frac{1}{\sqrt{2}}(|0\rangle - |1\rangle) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ -1 \end{pmatrix}$ with $\lambda = -1$. And for the $Y$-basis, $\sigma_y$ - $|+i\rangle = \frac{1}{\sqrt{2}}(|0\rangle + i|1\rangle) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ i \end{pmatrix}$ with $\lambda = +1$, - $|-i\rangle = \frac{1}{\sqrt{2}}(|0\rangle - i|1\rangle) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ -i \end{pmatrix}$ with $\lambda = -1$. The Pauli matrices satisfy important commutation relations: $$[\sigma_k, \sigma_l] = 2\mathrm{i}\epsilon_{klm}\sigma_m,$$ where $[\sigma_k, \sigma_l] = \sigma_k\sigma_l - \sigma_l\sigma_k$ is the *commutator*, $\mathrm{i}$ is the *imaginary unit*, the indices $k,l,m$ run over $\{1,2,3\}$ corresponding to $\{x\equiv 1,y\equiv 2,z\equiv 3\}$, we use the Einstein summation convention with summation over repeated index $m$, and $\epsilon_{klm}$ is the Levi-Civita symbol, $$\epsilon_{klm} = \begin{cases} +1 & \text{for even permutations of } (1,2,3) \\ -1 & \text{for odd permutations of } (1,2,3) \\ 0 & \text{if any indices are repeated} \end{cases} .$$ Explicitly, this gives $$[\sigma_x, \sigma_y] = [\sigma_1, \sigma_2] = 2\mathrm{i}\sigma_z,\quad i.e., \quad [\sigma_1, \sigma_2] = 2\mathrm{i}\sigma_3$$ $$[\sigma_y, \sigma_z] = [\sigma_2, \sigma_3] = 2\mathrm{i}\sigma_x,\quad i.e., \quad [\sigma_2, \sigma_3] = 2\mathrm{i}\sigma_1,$$ $$[\sigma_z, \sigma_x] = [\sigma_3, \sigma_1] = 2\mathrm{i}\sigma_y\quad i.e., \quad [\sigma_3, \sigma_1] = 2\mathrm{i}\sigma_2.$$ In 2D, these relations can be verified directly using the matrix representations $$\sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_2 = \begin{pmatrix} 0 & -\mathrm{i} \\ \mathrm{i} & 0 \end{pmatrix}, \quad \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}.$$ In QM, this cyclic structure reflects the rotational symmetry of the Pauli matrices and is fundamental to their role in describing spin-$1/2$ systems and qubit operations. They also satisfy *anticommutation relations* $\{\sigma_i, \sigma_j\} = 2\delta_{ij}I$, where $\delta_{ij}$ is the Kronecker delta. For each basis we also have the *completeness relations*, in the $\sigma_z$ basis, $|0\rangle\langle 0| + |1\rangle\langle 1| = I$, in the $\sigma_x$ basis, $|+\rangle\langle +| + |-\rangle\langle -| = I$, and in the $\sigma_y$ basis $|+i\rangle\langle +i| + |-i\rangle\langle -i| = I$. Any single-qubit unitary operation can be expressed as $$U = e^{i\alpha}e^{i(\theta_x\sigma_x + \theta_y\sigma_y + \theta_z\sigma_z)},$$ where $\alpha, \theta_x, \theta_y, \theta_z$ are real parameters. The non-commutativity leads to *uncertainty relations* $\Delta A \Delta B \geq \frac{1}{2}|\langle [A,B] \rangle|$, where the lower bound holds in a state-dependent manner, e.g., consider an example state $\|\psi\rangle = \|+\rangle \|\psi\rangle =\|+\rangle$ where $\langle \sigma_z \rangle \not= 0$. For Pauli operators, this means $\Delta \sigma_x \Delta \sigma_y \geq |\langle \sigma_z \rangle|$, with cyclic permutations. A measurement in any Pauli basis corresponds to projectors; for the $\sigma_z$ basis, $P_0 = |0\rangle\langle 0|$, $P_1 = |1\rangle\langle 1|$, for the $\sigma_x$ basis, $P_+ = |+\rangle\langle +|$, $P_- = |-\rangle\langle -|$, and for the $\sigma_y$ basis, $P_{+i} = |+i\rangle\langle +i|$, $P_{-i} = |-i\rangle\langle -i|$. The Bloch sphere representation indicates that any single-qubit state can be written as $|\psi\rangle = \cos(\theta/2)|0\rangle + e^{i\phi}\sin(\theta/2)|1\rangle$, where $\theta$ and $\phi$ are spherical coordinates on the Bloch sphere, with radius $r=1$. Finally, the *expectation values of Pauli operators* reflect the relation between the classical Cartesian coordinates and spherical coordinates $$x=\langle \sigma_x \rangle = \sin\theta\cos\phi, \quad y=\langle \sigma_y \rangle = \sin\theta\sin\phi, \quad z=\langle \sigma_z \rangle = \cos\theta.$$ ## Example: Qubit Phase Estimation Consider a qubit in the state $|\psi\rangle = \frac{1}{\sqrt{2}}\left(|0\rangle + e^{i\phi}|1\rangle\right),$ where $\phi$ is the *relative phase.* We will show how to probe the phase $\phi$ by collecting process measurement in 3 *different bases.* 1. *Measurement in the Computational Basis*, $Z$-basis, $|0\rangle, |1\rangle$: The probabilities $P(0) = |\langle 0|\psi\rangle|^2 = \frac{1}{2}$ and $P(1) = |\langle 1|\psi\rangle|^2=\frac{1}{2}$ are independent of $\phi$, and hence, no phase information is obtained purely by measurements in the $Z$-basis. 2. *Measurement in the* $X$-Basis, $|+\rangle, |-\rangle$: The states $|+\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ and $|-\rangle = \frac{1}{\sqrt{2}}(|0\rangle - |1\rangle)$ transform the phase into measurable probabilities, which reveals $\cos\phi$ $$P(+) = \frac{1 + \cos\phi}{2}, \quad P(-) = \frac{1 - \cos\phi}{2}.$$ 3. *Measurement in the* $Y$-Basis, $|y_+\rangle, |y_-\rangle$: Measuring the states $|y_+\rangle = \frac{1}{\sqrt{2}}(|0\rangle + i|1\rangle)$ and $|y_-\rangle = \frac{1}{\sqrt{2}}(|0\rangle - i|1\rangle)$ reveals $\sin\phi$ $$P(y_+) = \frac{1 + \sin\phi}{2}, \quad P(y_-) = \frac{1 - \sin\phi}{2}.$$ By combining results from the $X$- and $Y$-bases, both $\cos\phi$ and $\sin\phi$ are determined, allowing reconstruction of $\phi$ modulo $2\pi$. In essence, the phases are *encoded in interference patterns* that depend on the measurement basis. By choosing *multiple non-commuting (incompatible) bases* (e.g., $X$, $Y$, and $Z$), the relative phase information is *rotated* into observable probabilities. This underpins techniques like *quantum state tomography* and *interferometry*, where phase estimation relies on coherence across complementary observables. ### Adding a Commuting Basis to the $Y$-Basis Adds No New Phase Information When two observables *commute*, their eigenbases are identical (or *compatible*). Measuring in a commuting basis is equivalent to measuring in the original basis. Phase estimation requires measurements in complementary bases. To resolve $\phi$, we need measurements in *non-commuting bases* (e.g., $X$- and $Y$-bases). Complementary bases “probe” the phase through interference in distinct directions. For a single qubit, all observables commuting with $Y$ are trivial linear combinations of $Y$ and $I$, so no new eigenbases exist. In larger systems, degeneracy allows commuting observables with distinct eigenbases, but even there, relative phases are resolved via non-commuting measurements. Let's demonstrate that adding measurements in a new basis that commutes with the $Y$-basis adds no new phase information, still using the same *qubit example* and utilizing the Pauli-$Y$ operator. For a single qubit, since $Y$ has no degeneracy in its spectrum, any observable commuting with $Y$ must share its eigenbasis. Let’s define a trivial commuting operator, e.g., $\mathcal{O} = Y + \alpha I$, where $\alpha \in \mathbb{R}$ and $I$ is the identity operator. The new operator $\mathcal{O}$ shares the same *eigenstates* $|y_+\rangle$ and $|y_-\rangle$ as $Y$. Hence, measuring in the $\mathcal{O}$-basis is equivalent to measuring in the $Y$-basis itself. Again, let the qubit state be $|\psi\rangle = \frac{1}{\sqrt{2}}\left(|0\rangle + e^{i\phi}|1\rangle\right).$ A measurement in the $Y$-*Basis* $\left (|y_+\rangle, |y_-\rangle\right )$ gives probabilities that reveal $\sin\phi$, $$P(y_+) = \frac{1 + \sin\phi}{2}, \quad P(y_-) = \frac{1 - \sin\phi}{2}.$$ On the other hand, a measurement in the *new commuting basis*, $\mathcal{O}$ $\left (|y_+\rangle, |y_-\rangle\right )$, reflects observations in the eigenbasis of $\mathcal{O} = Y + \alpha I$. Since $\mathcal{O}$ shares its eigenstates with $Y$, the probabilities are identical $$P(y_+^\mathcal{O}) \equiv P(y_+)= \frac{1 + \sin\phi}{2}, \quad P(y_-^\mathcal{O}) \equiv P(y_-)= \frac{1 - \sin\phi}{2},$$ which indicates that *no new information* about $\phi$ is obtained by jointly measuring in the commuting basis. That is, measurements in $\mathcal{O}$ merely replicate the $Y$-basis result. Adding a commuting basis (e.g., $\mathcal{O} = Y + \alpha I$) provides *no independent data* about $\phi$. Only by measuring in *non-commuting bases* (e.g., $X$- and $Y$-bases) can we reconstruct the phase. This principle is foundational to quantum state tomography and phase estimation protocols like quantum Fourier transform. # From QM Wavefunction Phase Estimation to Kime-phase Estimation in Spacekime Representation Recall the earlier spacekime formulation of the complex-time representation (kime, $\kappa = t e^{i\theta}$) of repeated measurements corresponding to random draws from a time t-dependent kime-phase distribution ($\theta \sim \Phi(t)$). As the kime-phase is unobservable, we will explore strategies for introducing kime-measurement schemes, which are inspired by the quantum mechanical approach for recovering (or estimating) of the wavefunction phase, by taking repeated measurements in different non-commutative bases. Specifically, we'll extend "commutation" of variables (operators) in QM to distribution-action in kime-representation. This will provide a way to probe the enigmatic kime-phase by evaluating the action of different (time-dependent) kime-phase distributions $\Phi(t)$ on the class of kime-test functions, where the phase-distributions resemble the different bases in QM. ## Approach 1 *Approach 1* is a conceptual roadmap for designing a *kime-measurement* protocol that *mimics* quantum mechanical strategies for recovering wavefunction phase via repeated measurements in non-commuting bases. Quantum mechanics uses multiple, incompatible measurement bases to glean partial information about the complex amplitudes (including phases) of a quantum state. This approach attempts an *analogous* scheme in the kime framework by defining multiple families of kime-test functions - *“kime operators”* - whose “actions” on a hidden phase distribution $\Phi(t)$ can reveal partial information about $\Phi$. Effectively, we probe the “kime-phase” distribution $\Phi(t)$ by evaluating *several distinct* functionals (i.e., the actions of $\Phi$ on different test functions). If these families of test functions are in some sense “*incompatible*” (akin to *non-commuting bases*), then multiple measurements reveal more details (information) about the phase distribution. This is analogous to how QM measurements in multiple bases proble the phase of the quantum wavefunction. *In standard QM*, we measure a state $\lvert\psi\rangle$ in different bases, e.g., the $\{\lvert x\rangle\}$ basis for position, the $\{\lvert p\rangle\}$ basis for momentum. These bases are represented by *non-commuting observables* $\hat{X}$ and $\hat{P}$. *Each measurement* yields a probability distribution $\lvert\langle x\mid \psi\rangle\rvert^2$ or $\lvert\langle p\mid \psi\rangle\rvert^2$. Together, these distributions provide partial constraints on the full complex wavefunction (including the relative phases). *The non-commutation* ensures we *cannot* measure both bases simultaneously with arbitrary precision; hence we need repeated trials (fresh copies of the state $\lvert\psi\rangle$) to measure across different bases.\ Thus, *phase retrieval* in QM is intimately linked to measuring multiple incompatible observables. Each “basis” picks up a different interference pattern, revealing complementary aspects of the wavefunction’s phase structure. In the **Kime Picture**, we similarly probe the phase distribution $\Phi(t)$ by posing that the *kime-magnitude*, i.e., the usual time, $t = |\kappa|$, and the *kime-phase* $\theta$ is drawn from an unknown distribution $\Phi(\theta; t)$. We can treat $\theta$ as a “random quantity” unobservable in a direct sense - similar to how in standard QM, the wavefunction’s global phase is unobservable, but relative phases manifest in interference patterns. Our goal is to devise a scheme akin to “probing $\theta$” by making repeated measurements that “test” the distribution $\Phi$ with different *functionals* or *operators*, each of which partially encodes phase information. Following the classical distribution theory formalism, we'll consider the phase distirbution $\Phi(\theta; t)$ as a generalized function or measure. We only “see” $\Phi$ indirectly via its *action* on test functions $f(\theta)$, e.g., $\langle \Phi,\, f\rangle \;=\;\int f(\theta)\,\Phi(\theta;t)\,d\theta.$ A single *test function* $f$ can be viewed as a single “measurement setting,” from which we glean $\text{(Measured value)} \;=\; \int f(\theta)\,\Phi(\theta;t)\,d\theta.$ That is analogous to measuring an observable $\hat{A}$ in quantum mechanics and reading out $\langle \hat{A}\rangle_\psi$. To capture *kime-phase* information, we need multiple families of test functions that do *not* commute in some sense, i.e. cannot be simultaneously diagonalized or that “probe” $\Phi$ in *incompatible* ways. For instance, we can define two bases: - *Basis 1*: A set of *wavelet*-like test functions $\{f_\alpha(\theta)\}$.\ - *Basis 2*: A set of *harmonic* or Fourier-like test functions $\{g_\beta(\theta)\}$. The notion of “*non-commutation*” can be mirrored by the fact that the expansions in $f_\alpha$ vs. $g_\beta$ do *not* trivially coincide or diagonalize $\Phi$. In practice, we can impose a cross-condition such that $$\int f_\alpha(\theta)\,g_\beta(\theta)\,d\theta \;\neq\; 0 \quad\text{(overlaps in expansions)}.$$ Then measuring the “action” of $\Phi$ on $\{f_\alpha\}$ or on $\{g_\beta\}$ yields *complementary* (but not redundant) constraints on $\Phi$. To mimic *quantum tomography*, assume we have many repeated data sets (or repeated time-series) $\{f_i(t)\}$. In each data set, we effectively “test” $\Phi$ with one set of functions $\{f_\alpha\}$. In another set of repeated measurements, we test $\Phi$ with $\{g_\beta\}$. Those different test families might be *incompatible* in the sense that measuring one precludes measuring the other on the *same* sample - like a quantum system being collapsed by a certain measurement. But with multiple independent samples (fresh draws from $\theta$), we can systematically measure *both* families across runs. Pragmatic implementations of an operational kime-phase “measurement” scheme for real experiments requires further considerations. With *actual data*, the “kime phase” $\theta$ is not literally an angle that can be dialed up on an apparatus. Rather, each repeated measurement or repeated trial corresponds to a random draw from $\theta\sim\Phi(\theta;t)$. To “test” $\Phi(\theta)$ with the function $f_\alpha$, we need to design an analysis procedure that maps each trial’s observed data into an *aggregate* “score,” effectively $\langle \Phi, f_\alpha\rangle$. One scheme relies on defining $\widehat{A}_\alpha \;=\;\frac{1}{N}\,\sum_{i=1}^N f_\alpha\bigl(\widehat{\theta}_i\bigr),$ where $\widehat{\theta}_i$ is an *inferred* or *proxy* estimate of the phase from trial $i$. Or more simply, perform an integral transform on the repeated time-series so that the final numeric result is “the distribution’s action on $f_\alpha$.” To capitalize on the QM non-commutative operators insights, we can define *operators* $\widehat{F}_\alpha$ and $\widehat{G}_\beta$ whose matrix elements in a suitable “kime Hilbert space” do *not commute*. Then, $\widehat{F}_\alpha$ and $\widehat{G}_\beta$ represent two measurement families of $\theta$.\ In practice, these are realized by *distinct data analyses or distinct experimental set-ups* that yield functionals of $\Phi$. The specific details are more subtle, because standard quantum commutation $[\hat{X},\hat{P}]\neq 0$ is a strict operator statement in a Hilbert space, whereas here we have distribution-function “commutation.” But the spirit is that we define *distinct expansions or test function families* that cannot be *simultaneously* diagonal in the “kime phase.” Let's explore some parallels between the *quantum tomography* and the *kime-phase tomography* . In *quantum tomography*, we measure the state $\lvert\psi\rangle$ in many bases to reconstruct the *Wigner function* or the *density matrix.* Whereas, in *kime tomography*, we measure $\Phi(\theta;t)$ in many “test function expansions” or “operators” to reconstruct the phase distribution $\Phi$.\ Examples of kime-tomography bases include *Family A* (“Wavelet Test Functions) $\{f_\alpha(\theta)\}$ and *Family B* (“Fourier / Harmonic Functions) $\{g_\beta(\theta)\}$. In repeated sets of trials (each trial = one random $\theta_i$), we evaluate partial integrals (like sample means) that approximate $\int f_\alpha(\theta)\,\Phi(\theta)\,d\theta, \quad \int g_\beta(\theta)\,\Phi(\theta)\,d\theta.$ By combining these measurement outcomes across multiple $\alpha,\beta$ (assuming having enough repeated runs), we invert or fit $\Phi$ subject to consistency with all the observed integrals. If the sets $\{f_\alpha\}$ and $\{g_\beta\}$ have sufficient “resolution” (like a complete basis in $\theta$ space), we may recover $\Phi$. The *“non-commuting”* aspect reflects the need for separate runs to measure $\{f_\alpha\}$ vs. $\{g_\beta\}$, just as QM requires separate runs to measure position vs. momentum. In spacekime representaiton, *each measurement* is the distribution’s action on one family of test functions. Since each trial (or repeated measurement) gives one “realization” of $\theta$, we gather a large data set from many trials to estimate the integral $\int f(\theta)\,\Phi(\theta)\,d\theta$. *Non-commutation* arises if these families of test functions are not simultaneously diagonalizable or do not form a single “common basis.” In practice, that means we cannot glean the integrals for *both* sets from a *single* run. Instead, we record separate sets of repeated experiments, just as in QM we can’t measure $\hat{X}$ and $\hat{P}$ simultaneously on the same wavefunction copy. *Kime-phase tomography* combines the results from multiple “incompatible” families to gain partial (or complete) information about the distribution $\Phi(\theta)$. This is analogous to how measuring a quantum state in multiple bases yields full wavefunction tomography, including phase information. This *Approach 1* defines a kime-measurement protocol analogous to quantum mechanical phase retrieval by interpreting $\Phi(\theta)$ as a measure to be tested by distinct families of test functions (like measuring different observables), gather repeated data sets for each family separately, and use the results to reconstruct or constrain $\Phi(\theta)$ more fully than any single “basis” alone would allow. ### Strategy 1 Example This example demonstrates the *concept* of *kime-phase tomography* using two families of test functions (*wavelet-like* and *Fourier-like*) to probe a phase distribution $\Phi(\theta)$. The simulation involves: 1. Defining a “kime-phase distribution” $\Phi(\theta)$ supported on $\theta \in [-\pi,\pi]$; in this case, a simple mixture of *two Normal (or [von Mises](https://en.wikipedia.org/wiki/Von_Mises_distribution)) distributions*.\ 2. Simulating repeated measurements, where each measurement corresponds to drawing a random phase $\theta_i \sim \Phi,\ 1\leq i\leq N=200$ and generating a synthetic time-series.\ 3. Adopting two *measurement settings* - the “*Wavelet test basis*” and the “*Fourier test basis.*”. Half the measurements use the "Wavelet" family, and the other half use the "Fourier" family, mimicking non-commuting bases. For simplicity, the *wavelet family* uses $f_\alpha(\theta)=e^{-(\alpha\theta)^2} \times \cos(\alpha\theta)$; this is not a real wavelet, just a "toy wavelet-like function". Similarly, the Fourier family basis uses $g_\beta(\theta) = \cos(\beta \theta)$ and $\alpha,\beta \in \{1,2,3\}$. 4. Within each basis and for each measurement, we estimate an integral $\int f_\alpha(\theta)\,\Phi(\theta)\,d\theta$ or $\int g_\beta(\theta)\,\Phi(\theta)\,d\theta$ by collecting a “score” from the data, but from discrete random draws of $\theta$. This involves quick numerical approximations to $\mathbb{E}[f_\alpha(\theta)] = \int f_\alpha(\theta) \Phi(\theta) d\theta$ and $\mathbb{E}[g_\beta(\theta)] = \int g_\beta(\theta) \Phi(\theta) d\theta.$ 5. By combining or ensembling all the measured integrals across these two “*non-commuting*” families, we can solve for, fit, or approximate the kime-phase distribution $\Phi(\theta)$. For instance, one estimate of the kime-phase distribution can be obtained by $$\Phi_{est}(\theta) = \sum_{\alpha \in A} c_\alpha \times f_\alpha(\theta) + \sum_{\beta \in B} d_\beta \times g_\beta(\theta).$$ **Note**: This example is *proof-of-concept.* In a real fMRI or time-series context, we would need to design more sophisticated measurement operators (*test functions*) that link the phase $\theta$ to the measured signals. However, this simulation shows *in principle* the essential logic of **using multiple incompatible measurement families** to pin down (estimate) the phase distribution. ```{r kime-phase-tomography, echo=TRUE, message=FALSE, warning=FALSE} library(plotly) # ------------------------------------------------------- # 1. Define a hidden "kime-phase distribution" Phi(theta) # We'll do a simple mixture of two von Mises as a stand-in. # Range of theta: [-pi, pi]. # ------------------------------------------------------- dphi <- function(theta) { # mixture of two peaks for demonstration w <- 0.6 # we'll use circular normal approximation via dnorm with wrap-around # for simplicity, ignoring normalizing constants p1 <- dnorm(theta, mean=-1.0, sd=0.7) p2 <- dnorm(theta, mean=+1.5, sd=0.4) val <- w*p1 + (1-w)*p2 # we do a naive re-normalization across [-pi, pi] # in real code, better to do a proper normalization or use dvonmises return(val) } # Create a function to sample from this distribution by rejection or simple discretization samplePhi <- function(n) { # discrete approach: make a grid, accumulate, sample ngrid <- 2000 th_grid <- seq(-pi, pi, length.out=ngrid) pdf_vals <- dphi(th_grid) # normalize pdf_vals <- pdf_vals / sum(pdf_vals) # cumulative cdf_vals <- cumsum(pdf_vals) # random draws rU <- runif(n) # invert idx <- findInterval(rU, cdf_vals) return(th_grid[pmax(1, pmin(idx, ngrid))]) } # ------------------------------------------------------- # 2. Simulate repeated measurements: # In reality, each measurement i yields one random draw theta_i. # Then we get some time-series f_i(t). But for demonstration: # We'll just store the phase. We'll "measure" integrals in different # bases as if we had operators f_alpha or g_beta. # ------------------------------------------------------- N <- 200 # total repeated measurements theta_samples <- samplePhi(N) # We'll say half the measurements use the "Wavelet" family, # the other half use the "Fourier" family, mimicking non-commuting bases set_seed <- 123 idx_wav <- sample.int(N, size = N/2) idx_four <- setdiff(seq_len(N), idx_wav) # ------------------------------------------------------- # 3. Define test function families: # Family A: "Wavelet-like" => f_alpha(theta) # Family B: "Fourier-like" => g_beta(theta) # # We'll keep it small: alpha, beta in {1,2,3} as a toy. # ------------------------------------------------------- f_wavelet <- function(theta, alpha) { # Simple "Mexican-hat"-like or Morlet-like in theta-space # We'll do a naive shape, e.g. e^(-(alpha*theta)^2) * cos(alpha*theta) # Not a real wavelet, just a "toy wavelet-like function" return(exp(-0.5*(alpha*theta)^2)*cos(alpha*theta)) } g_fourier <- function(theta, beta) { # Standard Fourier basis = sin / cos. Let's do cos for simplicity: return(cos(beta*theta)) } alphas <- c(1,2,3) betas <- c(1,2,3) # We'll measure integrals for alpha in alphas or beta in betas. # ------------------------------------------------------- # 4. "Measure" partial integrals: # Each measurement i in wavelet basis => pick random alpha from {1,2,3}? # Or measure them all? # For simplicity, we assume each measurement i yields f_alpha(theta_i) # for alpha in a random subset or one alpha per measurement. # # Similarly for the Fourier basis g_beta. # # We'll store the resulting "observed integrals" in arrays. # ------------------------------------------------------- set.seed(999) obs_wavelet <- data.frame(alpha = integer(), val = numeric()) obs_fourier <- data.frame(beta = integer(), val = numeric()) for(i in idx_wav) { # pick alpha randomly from {1,2,3} to measure alpha_i <- sample(alphas, size=1) val_i <- f_wavelet(theta_samples[i], alpha=alpha_i) obs_wavelet <- rbind(obs_wavelet, data.frame(alpha=alpha_i, val=val_i)) } for(i in idx_four) { # pick beta randomly from {1,2,3} to measure beta_i <- sample(betas, size=1) val_i <- g_fourier(theta_samples[i], beta=beta_i) obs_fourier <- rbind(obs_fourier, data.frame(beta=beta_i, val=val_i)) } # We now have naive approximations to # E[f_alpha(\theta)] = \int f_alpha(\theta) Phi(\theta) d\theta # but from discrete random draws of theta. # ------------------------------------------------------- # 5. Estimate these integrals: # We'll do a simple average of obs vals for each alpha or beta # => This approximates \int f_alpha(\theta)\Phi(\theta)d\theta # ------------------------------------------------------- wavelet_est <- aggregate(val ~ alpha, data=obs_wavelet, FUN=mean) fourier_est <- aggregate(val ~ beta, data=obs_fourier, FUN=mean) wavelet_est fourier_est # wavelet_est$val[wavelet_est$alpha==alpha] is the numerical estimate # of int f_alpha(\theta) Phi(\theta) dtheta # ------------------------------------------------------- # 6. Representing Phi in a small parametric basis and solve: # We'll do a linear combination of the same wavelet & Fourier basis: # Phi_est(\theta) = sum_{alpha in A} c_alpha * f_alpha(theta) # + sum_{beta in B} d_beta * g_beta(theta). # # We'll choose alpha,beta in {1,2,3} => total of 6 unknown coefficients. # Then we match the integrals we measure with the predicted integrals # from the param. # ------------------------------------------------------- # Let alphaSet = {1,2,3}, betaSet={1,2,3} # We define: phi_est(\theta) = sum_{a} c_a f_wavelet(...) + sum_{b} d_b g_fourier(...) # The integral of f_wavelet wrt phi_est is: # \int f_wavelet_a(\theta) [ sum_{a'} c_{a'} f_wavelet_{a'}(\theta) # + sum_{b'} d_{b'} g_fourier_{b'}(\theta) ] dtheta # We'll discretize the integral on a fine grid in theta. Let's define a function that, # given (c_a, d_b), returns predicted integrals for each alpha or beta. theta_grid <- seq(-pi, pi, length.out=600) dtheta <- theta_grid[2]-theta_grid[1] # Precompute a large matrix of [f_wavelet_a(theta_grid), g_fourier_b(theta_grid)] for each alpha,beta fw_mat <- sapply(alphas, function(a) f_wavelet(theta_grid, a)) colnames(fw_mat) <- paste0("fw_a", alphas) gf_mat <- sapply(betas, function(b) g_fourier(theta_grid, b)) colnames(gf_mat) <- paste0("gf_b", betas) # We'll define a function to produce integrals int f_wavelet(alpha)*phi_est dtheta # and int g_fourier(beta)*phi_est dtheta, for a proposed vector of coefficients c & d: predictIntegrals <- function(coeffs_cd) { # Suppose length(coeffs_cd) = 6 => c_1..c_3, d_1..d_3 c_vec <- coeffs_cd[1:3] d_vec <- coeffs_cd[4:6] # Reconstruct phi_est on the grid: phi_est_grid <- fw_mat %*% c_vec + gf_mat %*% d_vec # Now compute the integral of f_wavelet_a * phi_est # That is (on the grid): # int( f_wavelet_a(theta_grid) * phi_est_grid ) dtheta # We'll produce a vector of length=3 for alpha=1..3: int_fw <- numeric(length(alphas)) for(i in seq_along(alphas)) { int_fw[i] <- sum( fw_mat[,i] * phi_est_grid ) * dtheta } # Similarly for the g_fourier(b): int_gf <- numeric(length(betas)) for(j in seq_along(betas)) { int_gf[j] <- sum( gf_mat[,j] * phi_est_grid ) * dtheta } return(list(wave=int_fw, four=int_gf)) } # We'll define a cost function: the squared difference between predicted integrals # and the measured integrals from wavelet_est, fourier_est wave_meas <- wavelet_est$val names(wave_meas) <- as.character(wavelet_est$alpha) four_meas <- fourier_est$val names(four_meas) <- as.character(fourier_est$beta) costFun <- function(par_cd) { pred <- predictIntegrals(par_cd) # match alpha=1..3 in order res_w <- pred$wave - wave_meas res_f <- pred$four - four_meas sum(res_w^2 + res_f^2) } # ------------------------------------------------------- # 7. Solve or minimize costFun: # ------------------------------------------------------- init_guess <- rep(0,6) fit <- optim(init_guess, costFun, method="BFGS") fit$par cat("Optimized cost:", fit$value, "\n") # fit$par => c_1..c_3, d_1..d_3 final_pred <- predictIntegrals(fit$par) # Compare final_pred vs wave_meas, four_meas data.frame( alpha=alphas, pred=round(final_pred$wave,3), obs=round(wave_meas,3) ) data.frame( beta=betas, pred=round(final_pred$four,3), obs=round(four_meas,3) ) # ------------------------------------------------------- # 8. Plot the reconstructed Phi vs. the true distribution # ------------------------------------------------------- # Reconstruct Phi_est on a grid coeffs_cd <- fit$par c_vec <- coeffs_cd[1:3] d_vec <- coeffs_cd[4:6] phi_est_grid <- fw_mat %*% c_vec + gf_mat %*% d_vec # True distribution on the same grid phi_true_grid <- sapply(theta_grid, dphi) # normalize each for plotting phi_true_grid <- phi_true_grid / (sum(phi_true_grid)*dtheta) phi_est_grid <- phi_est_grid / (sum(phi_est_grid)*dtheta) # might be negative parts, but let's see # par(mfrow=c(1,1)) # plot(theta_grid, phi_true_grid, type="l", col="blue", lwd=2, # ylim=range(c(phi_true_grid, phi_est_grid)), # xlab=expression(theta), ylab="density", main="Kime-phase tomography") # lines(theta_grid, phi_est_grid, col="red", lwd=2) # legend("topright", legend=c("True Phi","Estimated Phi"), col=c("blue","red"), lty=1, lwd=2) p <- plot_ly() %>% add_trace(x = theta_grid, y = phi_true_grid, type = 'scatter', ### True Phi mode = 'lines', line = list(color = 'blue', width = 2), name = 'True Phi') %>% add_trace(x = theta_grid, y = phi_est_grid, #### Estimated Phi type = 'scatter', mode = 'lines', line = list(color = 'red', width = 2), name = 'Estimated Phi') %>% layout(title = "Kime-phase tomography", xaxis = list(title = expression(theta)), yaxis = list(title = "Density"), legend = list(x = 0.75, y = 0.9)) p ``` In this simulation, the **hidden distribution** $\Phi(\theta)$ is a mixture of “circular Gaussians” (or an approximate approach with normal distributions truncated to $[- \pi,\pi]$). Simulating *repeated measurements* reflects sampling $\theta_i \sim \Phi$. In a real scenario, each $\theta_i$ would come from an *observable time series*, but for demonstration we directly sample $\theta_i$. We used *two measurement families*, a *wavelet-like* test functions $f_\alpha(\theta)$ and a *Fourier-like* test functions $g_\beta(\theta)$, which “non-commuting” in the sense that measuring $\theta_i$ with wavelet functions is separate from measuring with Fourier functions. For the wavelet or Fourier basis, we get $\{f_\alpha(\theta_i)\}$ or $\{g_\beta(\theta_j)\}$ by averaging over many $\theta_i$ to approximate $\int f_\alpha(\theta)\,\Phi(\theta)\,d\theta$. We define a parametric form $\Phi_{\text{est}}(\theta) = \sum c_\alpha\,f_\alpha(\theta) + \sum d_\beta\,g_\beta(\theta)$ and fit the coefficients by matching the measured integrals to the predicted integrals from $\Phi_{\text{est}}$.\ When the the underying phase distribution is complicated, or if we only have a small set of test functions, the reconstruction may be rough. Increasing the number of test functions (i.e., bigger sets $\alpha,\beta$) and the number of repeated measurements $N$ yields a better approximation, akin to how quantum state tomography improves with more measurement bases and more samples. In practice, to extend this approach to actual fMRI, or other longitudinal time-series data, requires defining how each measurement setting (wavelet or Fourier basis) is *implemented* as an operation on the actual time series. For example, *wavelet basis* means convolving the measured fMRI signal (in time) with a certain wavelet, then extracting a phase-dependent amplitude. For the *Fourier basis*, we can use FFT over a certain window that yields the coefficient at frequency $\omega$.\ Then, we *map* each measurement outcome to a “test function value” that depends on the hidden $\theta$. In real data, $\theta$ is not directly known, so we might do a specialized inference step or calibration. Ensembling (or aggregating) these measurement outcomes over many trials or repeated runs estimates the integrals $\int f_\alpha(\theta)\,\Phi(\theta)\,d\theta$. In the experiment above, the *estimated* phase distribution $\widehat{\Phi}$ (6-modes) roughly resembles the *true* distribution (bi-modal), which has two well-localized modes. In our experiment, the *estimated* phase distribution uses only a small basis of *six total functions*, $f_{1,2,3}$ wavelet-like and $g_{1,2,3}$ Fourier-like. This can be too rigid or produce “ringing” artifacts, adding more wavelet scales or higher Fourier frequencies can improve the approximation. Also, the balance between wavelet and Fourier terms can be adjusted, so that the distribution can better capture both peaked modes and smooth tails. In principle, using a bigger basis or a more flexible family supports smoother approximations with fewer spurious “wiggles.” Some wavelet or polynomial bases might be more natural for capturing multi-modal phase priors with rapidly decaying tails. For instance, *Legendre polynomials* or *higher B-splines* on $[-\pi, \pi]$ might yield fewer spurious oscillations than a small set of cosines. When the naive least-squares fit yields *overfitting* (excessive high-frequency structure), a *penalized regression* approach can help. For instance, *Tikhonov or Ridge penalty*\ $\min_{c,d} \left [\underbrace{\sum(\text{residuals}^2)}_{fidelity} + \underbrace{\lambda\,\bigl(\|c\|^2 + \|d\|^2\bigr)}_{regularizer}\right ],$ for some tuning $\lambda>0$ shrinks large oscillatory coefficients in the wavelet or Fourier expansions. This encourages fewer basis functions with large coefficients, simplifying the distribution’s shape. Assuming the phase distribution as a function in a *reproducing kernel Hilbert space (RKHS)*, then the solution is found by penalizing large second derivatives or local curvature. Also, in practice we can adopt a parametric mixture family closer to the true phase distribution. For instance, if the *true* distribution is a *two-component mixture of von Mises* ($vM$) (or circular Gaussians), it might be more direct to *parameterize* the estimate with the *same family*, e.g. $\Phi(\theta) \;=\; w\,vM(\theta;\mu_1,\kappa_1) + (1-w)\,vM(\theta;\mu_2,\kappa_2).$ Then fit the parameters $(w,\mu_1,\kappa_1,\mu_2,\kappa_2)$ by matching the measured integrals or the direct samples. This strategy is analogous to *modeling the distribution with a mixture of 2–3 ‘peaks’*, then solve for mixture weights and location/scale. It naturally yields a shape that has *only two modes* - no artificial high-frequency bumps. When dealing with only a small number of repeated measurements (draws of $\theta$) in each “basis,” then *sampling noise* can create spurious detail *increase* $N$ (the number of repeated draws) and/or *measure more integrals* in each basis or in additional “incompatible” bases. With more data constraints, the solution is less likely to overfit small random fluctuations. #### *Renormalization* of the Estimated Kime-phase Distibution In the *Approach 1* simulation, we see that the estimated (recovered) kime-phase distribution $\hat{\Phi}(t)$ may have negative values, which is unphysical/unrealistic. This can happen due to using a simple linear combination of wavelet-like and Fourier-like basis functions. Of course, this is an oversimplified protocol and *renormalization* of the reconstructed distribution $\hat{\Phi}(t)$ may address the issue by enforcing *positivity* and *unitarity* for interpreting $\Phi(\theta)$ as a *probability density.* In our simulation, we defined a *parametric estimate* of the form $$\Phi_{\text{est}}(\theta) = \sum_{\alpha} c_{\alpha} \, f_{\alpha}(\theta) + \sum_{\beta} d_{\beta} \, g_{\beta}(\theta),$$ where the basis functions $f_{\alpha}$ and $g_{\beta}$ typically take *positive* and *negative* values (e.g., wavelets, cosines). Thus, the linear combination can easily go negative at some $\theta$ if the coefficients $\{c_{\alpha}, d_{\beta}\}$ produce partial cancellations or overshoot. there are multiple *renormalization* strategies that can be used to address this issue and enforce kime-phase density properties on the reconstructed $\Phi_{\text{est}}(\theta)$. For instance, using unconstrained least-squares merely *minimizes squared errors* in matching integrals $\int f_{\alpha}\,\Phi(\theta)\,d\theta$ or $\int g_{\beta}\,\Phi(\theta)\,d\theta$, without imposing *any positivity constraint*. So, there is no mechanism to prevent negative portions in $\Phi_{\text{est}}$. A naive renormalization by *rescaling* the final estimate, $$\Phi_{\text{est}}(\theta) \leftarrow\; \frac{\Phi_{\text{est}}(\theta)} {\int \Phi_{\text{est}}(\theta)\,d\theta}$$ ensures *unitarity*, the integral is $1$, yet, the function can still have negative regions. Exponential parametrization provides a straightforward fix to represent $\Phi_{\text{est}}$ in an *exponential family* form, e.g., $\Phi_{\text{est}}(\theta) = \exp\Bigl(\,\sum_{\alpha} c_{\alpha}\,f_{\alpha}(\theta)\Bigr),$ or in a typical *Fourier or wavelet basis*, $$\Phi_{\text{est}}(\theta) = \exp\Bigl(\, c_{0} + \sum_{k=1}^K \bigl[c_k\,\cos(k\,\theta) \;+\; s_k\,\sin(k\,\theta)]\Bigr).$$ Then automatically, $\Phi_{\text{est}}(\theta)$ is *positive*, but still need to *renormalize* by dividing by its integral. Of course, this over-complexifies the problem, as the integrals $\int f_{\alpha}(\theta)\,\Phi_{\text{est}}(\theta)\,d\theta$ may no longer be linear in the coefficients $c_{\alpha}$. Another strategy is to choose basis functions that are *themselves* nonnegative, e.g. B-splines with compact support, or certain radial basis expansions.\ Then a linear combination with *nonnegative coefficients* ensures $\Phi_{\text{est}}(\theta)\ge 0$. In that case, we need to solve a *constrained* problem: $c_{\alpha}\ge 0$, which can be done with quadratic programming or more advanced optimization tools. Post-processing *clampping* is a simpler (though less principled) fix requiring computing the unconstrained estimate $\Phi_{\text{est}}(\theta)$, setting all negative values to $0$, and renormalizing over $\theta$.\ The “clamp-and-normalize” approach might not preserve integral constraints as precisely as a fully constrained method, but it ensures a nonnegative final result. Alternatively, we can normalize the distribution once we get a function $\Phi_{\text{est}}(\theta) \ge 0$, by $$\Phi_{\text{est}}(\theta) \leftarrow\; \frac{\Phi_{\text{est}}(\theta)} {\int_{-\pi}^{\pi} \Phi_{\text{est}}(\theta)\,d\theta}.$$ In practice, we can normalize over the grid `theta_grid` the functional values of `phi_est_grid`, by `phi_est_int <- sum(phi_est_grid) * dtheta` and `phi_est_grid <- phi_est_grid / phi_est_int`, assuming `phi_est_grid >= 0` so the integral is well-defined and positive.\ Then `phi_est_grid` forms a valid probability density function (pdf) over $[- \pi,\pi]$. #### Opportunities for Enhancement - *Strengthen the Analogy to Quantum Non-Commutativity*: In QM, non-commuting observables (e.g., position/momentum) enforce a fundamental limit on simultaneous measurements. In spacekime, "incompatibility" is loosely defined as non-orthogonal test functions (e.g., wavelet vs. Fourier). However, this does not directly replicate the operator-based non-commutativity of QM. The *no-overlap condition* $\int f_\alpha(\theta)\,g_\beta(\theta)\,d\theta \;\neq\; 0$ is insufficient to guarantee complementary information or a no-go theorem similar to the uncertainty principle. We can consider formalizing incompatibility using operator theory or signal-processing uncertainty principles (e.g., time-frequency trade-offs in wavelets vs. Fourier), and defining a commutator-like metric for test function families to quantify their mutual informativeness. - *Lack of Uniqueness and Stability*: The linear combination $\Phi_{\text{est}}(\theta) = \sum c_\alpha\,f_\alpha(\theta) + \sum d_\beta\,g_\beta(\theta)$ assumes $\Phi$ lies in the span of the chosen basis. If the test functions do not form a complete basis or are ill-conditioned, the inversion may be non-unique or unstable. We can explore using overcomplete dictionaries (e.g., frames) or adaptive basis selection and regularizing the inversion (e.g., Tikhonov, sparsity penalties) to handle noise and undersampling. - *Ambiguity in Phase Interpretation*: The kime-phase $\theta$ is treated as a random draw, but its physical meaning and connection to time-series data (e.g., fMRI) needs to be explicated. Unlike QM phases, which arise from wavefunction structure, $\theta$ does not have a direct operational definition. Is it feasible to ground $\theta$ as a measurable quantity (e.g., phase coherence in oscillatory signals) and define test functions as transformations of observable data (e.g., wavelet coefficients)? - *Simplistic Data Generation*: The simulation directly samples $\theta\sim\Phi_{true}(t)$, bypassing the critical step of inferring $\theta$ from the simulated fMRI time-series data. Perhaps explicate the mapping between signals and $\theta$, introducing noise and bias. Is it feasible to simulate time-series signals (e.g., synthetic fMRI) with known $\theta$ dependence? Can we use signal-processing methods (e.g., Hilbert transform, wavelet ridges) to estimate $\theta$ from data? - *Inefficient Data Usage*: Splitting trials into wavelet/Fourier groups mimics QM measurement collapse but wastes data. In reality, each $\theta_i$ could be analyzed with multiple test functions. Perhaps allow all test functions to act on each trial’s data to improve the statistical power. How to model the correlations between test function outputs. - *Operator-Theoretic Framework*: How to represent test functions as operators in a Hilbert space, with commutators defining incompatibility? Perhpas use spectral theory to formalize "measurements." Can we derive bounds on joint resolution of $\Phi$ using wavelet/Fourier families, similar to time-frequency uncertainty. - *Algorithmic Improvements*: Shoud we add $L_1/L_2$ penalties to the cost function to suppress spurious oscillations. Alternatively, see if Bayesian inference incorporating priors over $\Phi$ (e.g., smoothness, sparsity) can help quantify uncertainty in estimates. ### Analogies between QM Non-Commutativity and SK Representation with Incompatible Functional Bases In quantum mechanics (QM), we have an algebra of operators on a Hilbert space where two observables $\hat{A}$, $\hat{B}$ *do not commute* if $[\hat{A}, \hat{B}] \neq 0$. Non-commutation implies fundamental limits on simultaneous measurability (Heisenberg uncertainty). Mathematically, the *commutator* $[\hat{A}, \hat{B}] = \hat{A}\hat{B} - \hat{B}\hat{A}$ captures the idea that measuring $\hat{A}$ can disturb the outcome distribution of $\hat{B}$, or vice versa. By contrast, in the “kime-phase distribution” we have a distribution $\Phi(\theta)$ and a set of test functions $\{f_\alpha\}$, each giving a real number $\langle \Phi, f_\alpha\rangle$. “Incompatibility” reflects “families of test functions that cannot be measured on the *same* sample realization of $\theta$.” But *that* notion of “cannot measure both simultaneously” is not guaranteed by an operator algebra or a fundamental commutator—rather, it’s due to how we label draws from $\theta \sim \Phi$ and choose to do *different integral transforms* on separate runs. Hence, bridging from “nonorthogonal test functions” to a rigorous “operator non-commutation” requires new structural definitions. - Turning Test Functions into Operators: One way to introduce an operator-theoretic perspective is to treat each test function $f(\theta)$ as either: 1. *A multiplication operator*: $\hat{F} : \Psi(\theta) \;\mapsto\; f(\theta) \,\Psi(\theta)$ for $\Psi(\theta)$ in some function space (like $L^2$ over $\theta\in [-\pi,\pi]$), or 2. *An integral transform operator*: $\hat{F}(\Psi)(\theta) = \int K_f(\theta, \theta')\,\Psi(\theta')\,d\theta'.$ In both cases, each “test function” becomes an operator $\hat{F}$ and we can attempt to *compute commutators* $[\hat{F}, \hat{G}]$ for distinct test functions $f$ and $g$. If $\hat{F}$ is multiplication by $f(\theta)$ and $\hat{G}$ is multiplication by $g(\theta)$, then $[\hat{F}, \hat{G}] = 0$ trivially, because multiplication by $fg$ is commutative. So that alone does *not* replicate quantum non-commutation. When $\hat{F}$ or $\hat{G}$ includes derivatives, convolution kernels, or partial Fourier transforms, one can get non-zero commutators. - Reframe an “Uncertainty Relation” in $(\theta)$-Space: Even if we define a nontrivial operator algebra, we want an *operational meaning*, e.g., “it's impossible to measure both $\hat{F}$ and $\hat{G}$ simultaneously to arbitrary precision.” In QM, that arises from the Hilbert-space structure and the Born rule. In the kime-phase approach, we can enforce a “Hilbert norm” $\|\Psi\|$ on distributions or wavefunctions in $\theta$ and let “measurement” correspond to computing $\langle \Psi, \hat{F}\Psi\rangle$. When$[\hat{F}, \hat{G}] \neq 0$, we can explore a Cauchy–Schwarz–like inequality that implies $\Delta F\,\Delta G \ge \text{constant}.$ - Using Information-Theoretic Measures of Complementarity: Another route to formalizing “incompatibility” is to define how *informative* each measurement family is about $\theta$. Then one might show that if two families of test functions are “maximally” complementary, measuring one family leaves little or no information about the other. This can look like an uncertainty principle in the sense $H(\theta \mid \text{Family A data}) + H(\theta \mid \text{Family B data}) \;\ge\; \text{some bound,}$ where $H$ is an entropy measure, and “Family A data” is the set of integrals with the wavelet basis. If the two families produce highly overlapping expansions, we can glean both from the same data. But if the expansions are “complementary” in a strong sense, it might be impossible to gain both sets of integrals from a single sample. That suggests a trade-off reminiscent of quantum complementarity. - *Define a Kime Hilbert Space*: Let the “states” be (equivalence classes of) wavefunctions $\Psi(\theta)$ or distributions $\Phi(\theta)$. Impose an inner product $\langle \Psi_1, \Psi_2\rangle$. *Map Each “Test Function Family” to Self-Adjoint Operators* and instead of just computing $\int f(\theta)\,\Phi(\theta)\,d\theta$, define an operator $\hat{F}$ on the Hilbert space. Possibly $\hat{F}$ includes differentiation or convolution so that $\hat{F}\hat{G} \neq \hat{G}\hat{F}$. Then, evaluate $[\hat{F}, \hat{G}]$ for two such operators. If it is nonzero, investigate the consequences for “observables.” Specify how “observing $\hat{F}$” changes (or not) the state $\Psi$. It may not be nontrivial to get a truly quantum-like “collapse.” Skipping the collapse or disturbance postulate, we might only get partial analogies to quantum complementarity, not a full “non-commuting measurement” principle. Given that $[\hat{F},\hat{G}]\neq 0$, check if there is a limit or bound on the simultaneous determination of the integral values of $f(\theta)$ and $g(\theta)$. Possibly define an “uncertainty relation” in terms of standard deviations or entropies of these integral measurements. ### Refined Simulation Experiment (Approach 1) In this revised simulation, we use a Bayesian framework starting with some prior phase model, e.g., Laplace phase distribution, and using evidence-based iterative refinement to estimate the Bayesian posterior phase distribution, still using a “forward model” $$x(t) = \mathcal{M}\bigl(t,\theta; \,\boldsymbol{\phi}\bigr) +\epsilon(t),$$ where $\theta$ is the randomly drawn “kime-phase” from the currently estimated Bayesian posterior distribution, still $\theta(t)$ if time-varying, $\mathcal{M}$ is a known or hypothesized function that generates the fMRI time-series from $\theta$, $\boldsymbol{\phi}$ stands for additional model parameters (e.g., amplitude, baseline, shape, HRF in fMRI, etc.), and $\epsilon(t)$ is noise, e.g., Gaussian or realistic Rician fMRI noise with correlation in time. We can simulate many repeated fMRI runs or “trials,” each with a newly sampled $\theta_i$. For each run $i$, we produce $x_i(t) = \mathcal{M}\bigl(t,\theta_i; \,\boldsymbol{\phi}\bigr) + \epsilon_i(t).$ We will explore specific functions $\mathcal{M}$ that generate the fMRI time-series from the randomly drawn $\theta$, and update the entire simulation experiment. The revised simulation protocol involves the following steps, which are repeated until convergence to a stable estimate for the distribution $\Phi(\theta)$. 1. *Prior*: Start with a prior distribution over $\theta$, for instance $\Phi_0(\theta)$.\ 2. *Forward Simulation*: For each trial $i$, sample $\theta_i\sim \Phi_{n}(\theta)$, the *current* posterior or prior, depending on the iteration strategy, generate a noisy signal $x_i(t)$ via $\mathcal{M}(t,\theta_i;\boldsymbol{\phi}) + \epsilon_i(t)$.\ 3. *Inference*: Given $x_i(t)$, update the belief about $\theta_i$ and possibly about $\boldsymbol{\phi}$, producing a posterior $\Phi_{n+1}(\theta)$. Let's start with an *initial* prior, e.g. a **Laplace** or **Gaussian** distribution over $\theta$, $\Phi_0(\theta) \;\propto\; \exp\bigl(-|\theta|/b\bigr)\quad\text{(Laplace)},$ or $\Phi_0(\theta) \;\propto\; \exp\bigl(-\tfrac{\theta^2}{2\sigma^2}\bigr),$ restricted to $[- \pi,\pi]$. For time-varying $\theta(t)$, we can specify a small random walk prior or a correlated prior in time. To explicate the forward model, in each “trial” $i$, we draw $\theta_i \sim \Phi_n(\theta)$, generate the synthetic fMRI time series $x_i(t) = \mathcal{M}\bigl(t,\theta_i;\,\boldsymbol{\phi}\bigr) + \epsilon_i(t)$, and assume the noise $\epsilon_i(t)$ is Gaussian or Rician with temporal correlation. Again, $\boldsymbol{\phi}$ represents the hyerparameter vector containing amplitude, baseline, HRF shape, etc. Given the observed $x_i(t)$, the Bayesian posterior distribution, $\Phi_{n+1}(\theta)$, is updated from the prior $\Phi_n(\theta)$ by $$p(\theta \mid x_i) \propto \bigl[\,p(x_i\mid\theta)\bigr] \Phi_n(\theta).$$ We can combine data from multiple independent trials $\{x_i(t)\}_{i=1}^N$, $$p(\theta_1,\ldots,\theta_N \mid \{x_i\}) \propto\; \prod_{i=1}^N p(x_i \mid \theta_i)\,\Phi_n(\theta_i).$$ In the special case when all trials share a *single* $\theta$, which is less likely in typical repeated-measure scenarios, there would be a single posterior. We can use *MCMC* or a *variational* method to approximate the posterior distribution.\ Alternatively, a *particle filter* approach may keep a sample of “particles” for $\theta$, propagate them forward, weigh them by the likelihood of $x_i$, and resample. Yet, another option is to use *EM-like* iterative scheme when $\theta$ is high-dimensional or time-varying to perform partial expectation (E-step) updates. At the end, we get the posterior $\Phi_{n+1}(\theta)$. In the next iteration, we sample $\theta_i$ from that new distribution, generate new data, etc. We can consider the final $\Phi_{\mathrm{post}}(\theta)$ as the best estimate of the underlying “phase distribution.” Let's explore two alternative “forward model” functions $\mathcal{M}$ for simulating fMRI time-series from a random phase $\theta$. - *Model A*: Sinusoidal + HRF Convolution: $x(t) = \Bigl(A\,\sin\!\bigl(\omega \, t + \theta\bigr)\Bigr) \ast\; \mathrm{HRF}(t) + \text{drift}(t) + \epsilon(t),$ where the *sinusoidal* $A\,\sin(\omega t + \theta)$ is a canonical oscillation with unknown phase $\theta$, the *convolution* $(\cdot)\ast\mathrm{HRF}$ is a typical fMRI hemodynamic response function (e.g. a gamma function or canonical double-gamma), the *drift(t)* is a low-frequency polynomial or spline to mimic slow drift in real BOLD signals, $\epsilon(t)$ is noise, e.g., *Gaussian* $\epsilon\sim N(0,\sigma^2)$, *Rician*, or *AR(1)* correlated noise, and the *parameter vector* $\boldsymbol{\phi} = (A,\omega,\mathrm{HRF\_params}, \dots)$. Because $\theta$ modulates the *phase* of the sinusoid, each trial can shift the timing of the BOLD response in a consistent manner, where we can interpret $\theta$ as a “neuronal phase offset”. - *Model B*: Event-Related or Block-Design Onset Times: $x(t) = \Bigl(\sum_{k=1}^K \mathbf{1}\bigl[t\in (t_k+\theta,\,t_k+\theta+\Delta)\bigr]\Bigr) \ast\;\mathrm{HRF}(t) + \text{drift}(t) + \epsilon(t),$ where we assume an fMRI/BOLD signal with *block design* or *event design* with onsets $\{t_k\}$, the unknown phase $\theta$ *shifts all onsets* by $\theta$ units in time, so that the entire block or event structure is advanced or delayed, convolution with the HRF produce the BOLD-like signal, and the additional amplitude or baseline parameters and noise are similar to *Model A* above. Assuming that the “phase” is an offset in event timing across repeated runs, $\theta$ shifting the entire design is realistic. Each trial might place the block/event boundary earlier or later in time; effectively we measure how well the data fits that shift. the next experiment demonstrates the Bayesian approach in `R` and `Stan`. Specifically, we show an *end-to-end* simulation in a *Bayesian* framework for unknown kime-phases $\{\theta_i\}$. In this experiment, we define a forward model, a sinusoidal input at unknown phase $\theta$ convolved with a simple $HRF + noise$, simulate repeated runs, each with a fresh $\theta_i\sim \text{Laplace}(0,b)$, and infer those phases in `Stan` via MCMC. Running this simulation in a fresh `R` session requires the package `rstan` to be installed and a `C++` toolchain configured, see the [rstan docs](https://mc-stan.org/users/interfaces/rstan)). The entire protocol involves *Stan model compilation* and running MCMC, and produces a data frame comparing the *true* vs. *posterior mean* kime-phases. the simulation can be enhanced by changing the shape of the HRF, or by letting the HRF parameters be unknown and included in the Stan model, by modifying the noise distribution (e.g. Rician or correlated), and by updating the number of runs, sampling intervals, amplitude priors, etc. In a real fMRI context, for dealing with multi-voxel data, multi-parameter block designs, or a richer prior on $\theta$, this *end-to-end protocol* may need to be refined. ```{r approachA_Refined1, echo=TRUE, message=FALSE, warning=FALSE} ############################################################ # 1. ENV SETUP ############################################################ # Make sure rstan and plot packages arer installed: # install.packages(c("rstan", "ggplot2")) library(rstan) # for Stan interface library(ggplot2) # for basic plotting # rstan config for quicker examples rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) ############################################################ # 2. DEFINE A SIMPLE HRF AND HELPER FUNCTIONS ############################################################ # A tiny "canonical HRF": single gamma-like shape # for demonstration. Typically using a more realistic double-gamma, etc. make_hrf <- function(length_out=32, peak=6, scale=1) { # We'll define a gamma PDF shape ~ dgamma(t, shape=?, rate=?) tseq <- seq(0, length_out-1, by=1) # shape=peak, rate=peak? => quick guess # this can be refined hrf_values <- dgamma(tseq, shape=peak, rate=peak/scale) hrf_values / max(hrf_values) # normalize peak to 1 } # Laplace sampler for phases rLaplace <- function(n, mu=0, b=1) { # For each uniform draw u in (-1/2,1/2), transform u <- runif(n, min=-0.5, max=0.5) sapply(u, function(ui) mu - b*sign(ui)*log1p(-2*abs(ui))) } ############################################################ # 3. SIMULATE DATA (FORWARD MODEL) ############################################################ simulate_data <- function(R, N, A=1, freq=0.33, noise_sd=0.2, b=1.0, t_max=30) { # R = # of runs (repeated measurements) # N = # of time points per run # freq => frequency in cycles/sec (just for demonstration) # b => Laplace scale for phase draws # t_max => total length in time # Return a list with: # t_vals: the time vector # x_mat: a matrix (R x N) of noisy signals # theta_true: the actual phases used # time axis t_vals <- seq(0, t_max, length.out=N) # define HRF hrf_vec <- make_hrf(length_out=32, peak=6, scale=1) # sample phases from Laplace theta_true <- rLaplace(R, 0, b) x_mat <- matrix(NA, nrow=R, ncol=N) # convolve each run's sinusoid with HRF for(r in 1:R) { # create raw sinusoid with phase shift raw_signal <- A * sin(2*pi*freq * t_vals + theta_true[r]) # do discrete convolve conv_full <- convolve(raw_signal, rev(hrf_vec), type="open") # match length N conv_full <- conv_full[1:N] # add noise noisy <- conv_full + rnorm(N, 0, noise_sd) x_mat[r, ] <- noisy } list(t_vals=t_vals, x_mat=x_mat, theta_true=theta_true) } ############################################################ # 4. STAN MODEL CODE ############################################################ # We'll assume each run has an unknown phase theta[r]. # We'll do a single amplitude A shared across runs, or let each run have its own amplitude. # We'll keep the HRF fixed (like a known shape). # # The Stan model tries to recover theta[r] by matching the data x[r,] with # x_hat[r, n] = A*sin(2*pi*freq * t_vals[n] + theta[r]) convolved with known HRF # Then a normal likelihood with std dev sigma_noise. stan_model_code <- " data { int R; // number of runs int N; // number of time points per run real t_vals[N]; // time vector matrix[R, N] x_mat; // observed data int len_hrf; vector[len_hrf] hrf_vec; // known HRF shape real freq; // known frequency } parameters { real A; // amplitude real sigma_noise; // noise std vector[R] theta_raw; // unconstrained phases } transformed parameters { // We'll map theta_raw into [-pi, pi] range just for demonstration // We can do a tanh approach or we can do mod vector[R] theta; for(r in 1:R) { // let's do a tanh approach: [-1,1] => [-pi, pi] theta[r] = pi() * tanh(theta_raw[r]); } } model { // Priors A ~ normal(1, 1); // guess amplitude near 1 sigma_noise ~ normal(0.2, 0.1) T[0,]; theta_raw ~ normal(0, 1); // broad for phase // likelihood for(r in 1:R) { // Build raw sinusoid vector[N] raw_signal_r; for(n in 1:N) { raw_signal_r[n] = A * sin(2*pi()*freq * t_vals[n] + theta[r]); } // convolve with hrf => we do a discrete approach // need a new vector conv_full, length N // We'll do naive loop convolution vector[N] conv_full = rep_vector(0, N); for(i in 1:N) { // i => index in conv result // sum over j => raw_signal_r[j]*hrf_vec[i-j+1] real tmp_sum = 0; for(j in 1:i) { int k = i-j+1; // index in hrf if(k <= len_hrf) { tmp_sum += raw_signal_r[j]*hrf_vec[k]; } } conv_full[i] = tmp_sum; } // now conv_full is the convolved signal // likelihood x_mat[r] ~ normal(conv_full, sigma_noise); } } " ############################################################ # 5. RUN SIMULATION ############################################################ set.seed(1234) R <- 8 # number of runs N <- 60 # number of timepoints per run freq_used <- 0.33 sim_out <- simulate_data(R=R, N=N, A=1.0, freq=freq_used, noise_sd=0.15, b=1.0, t_max=20) # We store t_vals, x_mat, theta_true t_vals <- sim_out$t_vals x_mat <- sim_out$x_mat theta_true <- sim_out$theta_true hrf_vec <- make_hrf(length_out=32, peak=6, scale=1) ############################################################ # 6. PREPARE DATA FOR STAN + COMPILE MODEL ############################################################ stan_data <- list( R = R, N = N, t_vals = t_vals, x_mat = x_mat, len_hrf = length(hrf_vec), hrf_vec = hrf_vec, freq = freq_used ) # compile stan_model <- stan_model(model_code=stan_model_code) ############################################################ # 7. SAMPLE / FIT MCMC ############################################################ fit <- sampling(stan_model, data=stan_data, iter=2000, chains=2, seed=999, control=list(adapt_delta=0.9)) # print(fit, pars=c("A","sigma_noise")) fit theta_draws <- extract(fit, pars="theta")$theta dim(theta_draws) # should be iterations x R ############################################################ # 8. CHECK PHASE ESTIMATES vs. TRUE ############################################################ # We'll compute the posterior mean for each run's phase theta_est_mean <- apply(theta_draws, 2, mean) theta_est_sd <- apply(theta_draws, 2, sd) df_phases <- data.frame( run_id = 1:R, theta_true = theta_true, theta_est = theta_est_mean, est_sd = theta_est_sd ) cat("\nPosterior mean vs. true phase:\n") df_phases # Quick plot # ggplot(df_phases, aes(x=theta_true, y=theta_est)) + # geom_point(color="blue") + # geom_abline(intercept=0, slope=1, linetype="dashed") + # labs(title="True vs. estimated phase", x="True θ", y="Estimated θ") + # coord_fixed() # Convert ggplot to plotly library(plotly) # Assuming df_phases contains theta_true and theta_est columns # First create the reference line # Get the range of data for the reference line range_min <- min(min(df_phases$theta_true), min(df_phases$theta_est)) range_max <- max(max(df_phases$theta_true), max(df_phases$theta_est)) ref_line <- data.frame(x = c(range_min, range_max), y = c(range_min, range_max)) # Create the plot fig <- plot_ly() %>% # Add scatter points add_trace(data = df_phases, x = ~theta_true, y = ~theta_est, type = 'scatter', mode = 'markers', marker = list(color = 'blue',size = 6), name = 'Phase Values') %>% # Add reference line add_trace(data = ref_line, x = ~x, y = ~y, type = 'scatter', mode = 'lines', line = list(dash = 'dash', color = 'black'), name = 'y = x') %>% # Layout settings layout(title = list(text = "True vs. estimated phase", y = 0.95), xaxis = list(title = "True θ", scaleanchor = "y", # This ensures equal scaling scaleratio = 1 # This maintains aspect ratio ), yaxis = list(title = "Estimated θ"), showlegend = TRUE) fig ``` This *Bayesian* approach ensures each data set informs our knowledge of $\theta$ (and any other parameters). Over repeated runs, we effectively refine $\Phi(\theta)$. It also naturally handles noise, bias, and partial identifiability in $\theta$. This version of the simulation shows incorporating the kime representation in a Bayesian framework starting with a prior *phase distribution*, a forward model for generating signals from that phase, repeated trials, and an inference step that yields a posterior distribution of phases. This may be a more realistic approach than the simpler simulation shown earlier that “directly samples $\theta\sim\Phi$ and treats it as known.” ## Approach 1A: A Hilbert Space and Operator Reformulation In the above kime-phase estimation framework, the heuristic notion of "*incompatibility*" between test function families (e.g., $f_\alpha$ vs. $g_\beta$) based solely on *non-orthogonality*, $\int f_\alpha(\theta) g_\beta(\theta) d\theta \neq 0$, does not guarantee complementary information. This is unlike quantum mechanics (QM), where non-commutativity arises from operator algebra in a *Hilbert space.* Hence, this reformulation addresses that downside by defining a formal operator framework for the kime-test functions in $L^2[-\pi, \pi]$, computing their commutators explicitly, and exploring the conditions under which non-commutativity can be achieved. This may provide a more rigorous basis for *kime-phase tomography* analogous to *QM tomography*. **Note**: In the section below, we reflect on the meaning of $\frac{d}{d\theta}$ as a derivative operator in $L^2[-\pi, \pi]$ over the $\theta$-domain of $\Phi(\theta; t)$, not a stochastic derivative. $\theta$ is a coordinate, not a process, so standard calculus applies. There is no need for Itô calculus, which applies to stochastic processes where a variable (e.g., $\theta(t)$) *evolves randomly over time* $t$, governed by a stochastic differential equation (SDE). In kime, $\theta$ is *sampled independently* at each $t$ from $\Phi(\theta; t)$, with no explicit temporal dynamics (e.g., $d\theta = \mu dt + \sigma dW_t$) -- more precisely, there’s no continuous stochastic process or differential equation defining how $\theta$ evolves from $t$ to $t + dt$. Each $\theta_n(t)$ is a fresh draw, not a step in a trajectory. The independence of samples across $t$ means no explicit temporal correlation (e.g., autocorrelation) is modeled between $\theta_n(t)$ and $\theta_n(t')$. The derivative $\frac{d}{d\theta}$ operates over the $\theta$-domain of the density $\Phi(\theta; t)$ at fixed $t$, not over $t$. Thus, Itô calculus is unnecessary here. The randomness is in the sampling process, not in a continuous evolution of $\theta$. For a fixed measurement $n$, $\theta_n(t)$ and $\theta_n(t')$ are independent draws from $\Phi(\theta; t)$ and $\Phi(\theta; t')$, respectively. There’s no correlation imposed between $\theta_n(t)$ and $\theta_n(t')$ by a temporal process—just new, independent samples at each $t$. This supports the "no explicit temporal dynamics" claim. At the same time, $\Phi(\theta; t)$ varies with $t$ and the ensemble of samples $\{\theta_n(t)\}_{n=1}^N$ at time $t$ has different statistical characteristics, e.g., mean $\mu(t)$, concentration $\kappa(t)$, than at $t'$. This introduces a form of *temporal dependence* in the distribution, *not in the individual* $\theta_n(t)$ trajectories. - *Non-Commutativity*: $[\hat{\Theta}, \hat{G}_\beta] = i \beta I$ holds, providing a QM-like foundation for complementary information. - *Kime Context*: Operators act on $\psi(\theta) = \sqrt{\Phi}$, with randomness handled via sampling for data generation, not operator definition. In the *initial Approach 1* formulation, two families of test functions *wavelet-like* ($f_\alpha(\theta)$) and *Fourier-like* ($g_\beta(\theta)$) were used to probe the phase distribution $\Phi(\theta; t)$. It argues these families are "incompatible" if they are non-orthogonal, implying they should provide complementary information about $\Phi$. However, non-orthogonality only ensures overlap in $L^2[-\pi, \pi]$, not a fundamental limit on simultaneous measurement or complementary resolution, as in QM’s uncertainty principle. In QM, operators $A$ and $B$ are non-commuting ($[A, B] \neq 0$) when they *cannot be simultaneously diagonalized*, leading to $\Delta A \Delta B \geq \frac{1}{2}|\langle [A, B] \rangle|$. The kime framework needs a similar operator-theoretic foundation. In thsi *reformulation of Aproach 1*, we define the operators $\hat{F}_\alpha$ and $\hat{G}_\beta$ in $L^2[-\pi, \pi]$ corresponding to $f_\alpha$ and $g_\beta$, and compute their commutators explicitely. **Definition** [Kime Hilbert Space and Operators]. Consider the Hilbert space, $L^2[-\pi, \pi]$, of square-integrable functions over $\theta \in [-\pi, \pi]$, equipped with the inner product $$\langle \psi_1, \psi_2 \rangle = \int_{-\pi}^\pi \psi_1(\theta) \overline{\psi_2(\theta)} d\theta,$$ and the norm $\|\psi\|^2 = \langle \psi, \psi \rangle$. The phase distribution $\Phi(\theta; t)$ is a *probability density* (non-negative, $\int \Phi d\theta = 1$), and kime test functions acting on $\Phi$ via $\langle \Phi, f \rangle = \int f(\theta) \Phi(\theta) d\theta$. Let's explicate the operators whose expectation values yield these integrals. Consider the simplest case of *multiplication operators* where we define $\hat{F}_\alpha$ and $\hat{G}_\beta$ as multiplication operators - $\hat{F}_\alpha \psi(\theta) = f_\alpha(\theta) \psi(\theta)$, - $\hat{G}_\beta \psi(\theta) = g_\beta(\theta) \psi(\theta)$, where $\psi(\theta) \in L^2[-\pi, \pi]$, and $f_\alpha, g_\beta$ are bounded functions, e.g., $f_\alpha(\theta) = e^{-(\alpha \theta)^2} \cos(\alpha \theta)$, $g_\beta(\theta) = \cos(\beta \theta)$, to ensure $\hat{F}_\alpha, \hat{G}_\beta$ are well-defined operators. Then, we can compute their commutator by explicating each of its two components - $\hat{F}_\alpha \hat{G}_\beta \psi(\theta) = f_\alpha(\theta) (\hat{G}_\beta \psi(\theta)) = f_\alpha(\theta) g_\beta(\theta) \psi(\theta)$, - $\hat{G}_\beta \hat{F}_\alpha \psi(\theta) = g_\beta(\theta) (\hat{F}_\alpha \psi(\theta)) = g_\beta(\theta) f_\alpha(\theta) \psi(\theta)$. Since multiplication is commutative, $f_\alpha(\theta) g_\beta(\theta) = g_\beta(\theta) f_\alpha(\theta)$, $$[\hat{F}_\alpha, \hat{G}_\beta] = \hat{F}_\alpha \hat{G}_\beta - \hat{G}_\beta \hat{F}_\alpha = 0.$$ Therefore, the *multiplication operators*, $\hat{F}_\alpha$ and $\hat{G}_\beta$ commute, implying that they can be simultaneously diagonalized; their eigenfunctions are Dirac deltas $\delta(\theta - \theta_0)$. This does not replicate QM non-commutativity and provides no basis for an uncertainty relation or complementary information. **Example 1 [Commutation of Multiplication and Differentiation]**: The first example of non-commuting operators borrows intuition from QM, where the *position* $\hat{x}=x$ and *momentum* $\hat{p} = -i \frac{d}{d\theta}$ operators are non-commutative, $[\hat{x}, \hat{p}] = i$, which indicates their differential nature. Let's show that the kime *differential* and *multiplication* operators are *non-commutative.* Define - $\hat{F}_\alpha = \hat{\Theta}$, the *multiplication operator* by $\theta$ (QM position-like) $\hat{\Theta} \psi(\theta) = \theta \psi(\theta)$, and - $\hat{G}_\beta = \hat{P} = -i \frac{d}{d\theta}$, the *derivative operator* (QM momentum-like). Therefore, - $\hat{F}_\alpha \hat{G}_\beta \psi(\theta) = \hat{\Theta} (-i \frac{d}{d\theta} \psi(\theta)) = -i \theta \frac{d}{d\theta} \psi(\theta)$, and - $\hat{G}_\beta \hat{F}_\alpha \psi(\theta) = -i \frac{d}{d\theta} (\hat{\Theta} \psi(\theta)) = -i \frac{d}{d\theta} (\theta \psi(\theta)) = -i \left( \psi(\theta) + \theta \frac{d}{d\theta} \psi(\theta) \right)$, Hence, their commutator is $[\hat{F}_\alpha, \hat{G}_\beta] \psi(\theta) = -i \theta \frac{d}{d\theta} \psi(\theta) - \left(-i \psi(\theta) - i \theta \frac{d}{d\theta} \psi(\theta)\right) = i \psi(\theta)$. In operator form, $[\hat{\Theta},\hat{P}]\equiv [\theta, -i \frac{d}{d\theta}] = i I,$ where $I$ is the identity operator. This mirrors the QM canonical commutation relation $[\hat{x}, \hat{p}] = i \hbar$; here, $\hbar = 1$. To explicate the induced *uncertainty relations*, let's first examine the *expectation*. Assuming $\Phi$ is the true kime-phase density, given a state $\psi(\theta) = \sqrt{\Phi(\theta)}$, The *expected values* of the *multiplication* and *derivative* kime operators are - (multiplication) $\langle \hat{F}_\alpha \rangle\equiv \langle \hat{\Theta} \rangle = \int_{-\pi}^\pi \theta \Phi(\theta) d\theta$, - (derivative) $\langle \hat{G}_\beta \rangle\equiv \langle -i \frac{d}{d\theta} \rangle = \int_{-\pi}^\pi \sqrt{\Phi(\theta)} \left(-i \frac{d}{d\theta} \sqrt{\Phi(\theta)}\right) d\theta$, which requires boundary conditions, e.g., periodicity or rapid decay at $\pm \pi$. Then, the *uncertainty relation* is $\Delta \Theta \Delta P \geq \frac{1}{2} |\langle [\hat{\Theta}, \hat{P}] \rangle| = \frac{1}{2},$ where $\hat{P} = -i \frac{d}{d\theta}$, providing a rigorous basis for complementary information. To linking this relation to kime-test functions, instead of using $f_\alpha(\theta) = e^{-(\alpha \theta)^2} \cos(\alpha \theta)$ as we did earlier, now we define the multiplicaiton and derivative operators as - $\hat{F}_\alpha = \hat{\Theta}$, probing the "*position*" of $\theta$, and - $\hat{G}_\beta = -i \beta \frac{d}{d\theta}$, sensing the *scaled momentum*, tied to frequency-like behavior (similar to Fourier modes). The actions $\langle \Phi, \hat{F}_\alpha \Phi \rangle = \int \theta \Phi(\theta)^2 d\theta$ and $\langle \Phi, \hat{G}_\beta \Phi \rangle = -i \beta \int \Phi(\theta) \frac{d}{d\theta} \Phi(\theta) d\theta$ yield *moments* and *derivatives* of $\Phi$, respectively, offering complementary views and independet information about hte kime-phase distribution. **Example 1 [Commutation of Convolution and Multiplication]**: The second example of non-commutation involves the convolution operator for wavelet-like behavior. Multiplication by a kime-test function $f_\alpha(\theta)$ doesn’t capture the spatial structure of wavelets. Defining $\hat{F}_\alpha$ as a convolution operator - $\hat{F}_\alpha \psi(\theta) = \int_{-\pi}^\pi f_\alpha(\theta - \theta') \psi(\theta') d\theta'$ (convolution), - $\hat{G}_\beta \psi(\theta) = \cos(\beta \theta) \psi(\theta)$ (multiplication, Fourier-like). The pair of components of the commutator are - $\hat{F}_\alpha \hat{G}_\beta \psi(\theta) = \int_{-\pi}^\pi f_\alpha(\theta - \theta') \cos(\beta \theta') \psi(\theta') d\theta'$, and - $\hat{G}_\beta \hat{F}_\alpha \psi(\theta) = \cos(\beta \theta) \int_{-\pi}^\pi f_\alpha(\theta - \theta') \psi(\theta') d\theta'$. For this kime-test function, $f_\alpha(\theta) = e^{-(\alpha \theta)^2} \cos(\alpha \theta)$, - $\hat{F}_\alpha \hat{G}_\beta \psi(\theta) = \int_{-\pi}^\pi e^{-(\alpha (\theta - \theta'))^2} \cos(\alpha (\theta - \theta')) \cos(\beta \theta') \psi(\theta') d\theta'$, and - $\hat{G}_\beta \hat{F}_\alpha \psi(\theta) = \cos(\beta \theta) \int_{-\pi}^\pi e^{-(\alpha (\theta - \theta'))^2} \cos(\alpha (\theta - \theta')) \psi(\theta') d\theta'$. Hence, the *(Convolution-Multiplication) Commutator* $[\hat{F}_\alpha, \hat{G}_\beta] \neq 0$ because $\cos(\beta \theta)$ varies with $\theta$, while the convolution kernel shifts and scales differently. The *exact computation* is complex, but numerically for $\psi(\theta) = 1$, $\hat{F}_\alpha \hat{G}_\beta 1 \neq \hat{G}_\beta \hat{F}_\alpha 1$, indicating non-commutativity. This aligns with the wavelet vs. Fourier intuition where *convolution* (local structure) and *multiplication* (global oscillation) probe the kime-phase $\Phi$ *differently.* Let’s implement and validate this in a simulation confirming the non-commutativity and complementarity of the *multiplication* and *differentiation* operators. ```{r approach1A_Part1, warning=FALSE, message=FALSE} library(dplyr) # Define the Hilbert space: theta in [-pi, pi] theta <- seq(-pi, pi, length.out = 1000) dtheta <- theta[2] - theta[1] # Test functions f_alpha <- function(theta, alpha = 1) exp(-(alpha * theta)^2) * cos(alpha * theta) g_beta <- function(theta, beta = 1) cos(beta * theta) # Operators F_alpha <- function(psi, theta, alpha = 1) { # Convolution with f_alpha sapply(theta, function(th) { kern <- f_alpha(th - theta, alpha) sum(kern * psi) * dtheta }) } G_beta <- function(psi, theta, beta = 1) { # Multiplication by g_beta g_beta(theta, beta) * psi } # Derivative operator for comparison D_beta <- function(psi, theta, beta = 1) { -1i * beta * diff(c(psi[length(psi)], psi))[1:length(theta)] / dtheta } # Commutator commutator <- function(op1, op2, psi, theta, p1 = 1, p2 = 1) { op1_op2 <- op1(op2(psi, theta, p2), theta, p1) op2_op1 <- op2(op1(psi, theta, p1), theta, p2) op1_op2 - op2_op1 } # Test with uniform psi psi <- rep(1, length(theta)) / sqrt(2 * pi) # Normalized comm_FG <- commutator(F_alpha, G_beta, psi, theta, 1, 1) comm_ThetaD <- commutator(function(psi, theta, a) theta * psi, D_beta, psi, theta, 1, 1) # Norms norm_FG <- sqrt(sum(abs(comm_FG)^2) * dtheta) norm_ThetaD <- sqrt(sum(abs(comm_ThetaD)^2) * dtheta) cat("Norm of [F_alpha, G_beta] =", norm_FG, "\n") # Non-zero cat("Norm of [Theta, D_beta] =", norm_ThetaD, "\n") # Non-zero # Probe a distribution Phi <- dnorm(theta, 0, 0.5); Phi <- Phi / sum(Phi * dtheta) exp_F <- sum(F_alpha(Phi, theta, 1) * Phi) * dtheta exp_G <- sum(G_beta(Phi, theta, 1) * Phi) * dtheta cat("Expectation F_alpha =", exp_F, "\n") cat("Expectation G_beta =", exp_G, "\n") ``` Note that $[F_\alpha, G_\beta] \neq 0$ confirms non-commutativity for convolution vs. multiplication and $[\Theta, D_\beta] \neq 0$ matches QM expectations. The expectations $\langle \Phi, \hat{F}_\alpha \Phi \rangle$ and $\langle \Phi, \hat{G}_\beta \Phi \rangle$ differ, suggesting complementary information about $\Phi$. The next simulation reflects the *Approach 1A: Hilbert Space and Operator Reformulation*. This code replaces the *test functions* in *Apprach 1* with *operator-based measurements* (position and momentum), validates their non-commutativity via an uncertainty principle, and reconstructs the phase distribution using both operators. ```{r approach1A_Part2, warning=FALSE, message=FALSE} library(ggplot2) library(dplyr) # 1. Define True Kime-Phase Distribution (Mixture of Von Mises) # ------------------------------------------------------------- dphi <- function(theta, kappa1=5, kappa2=5, mu1=-1, mu2=1.5, w=0.6) { # Mixture of two von Mises distributions vm1 <- exp(kappa1 * cos(theta - mu1)) / (2 * pi * besselI(kappa1, 0)) vm2 <- exp(kappa2 * cos(theta - mu2)) / (2 * pi * besselI(kappa2, 0)) w * vm1 + (1 - w) * vm2 } # 2. Sample Theta from Distribution # --------------------------------- set.seed(123) theta_grid <- seq(-pi, pi, length.out=1000) phi_true <- dphi(theta_grid) theta_samples <- sample(theta_grid, 1000, replace=TRUE, prob=phi_true) # 3. Estimate Phi using Kernel Density Estimation (KDE) # ----------------------------------------------------- bw <- 0.3 # Bandwidth kde <- density(theta_samples, from=-pi, to=pi, bw=bw, n=1000) phi_est <- approx(kde$x, kde$y, theta_grid)$y phi_est <- phi_est / sum(phi_est * diff(theta_grid)[1]) # Normalize # 4. Compute Position (Theta) Variance # ------------------------------------ mean_theta <- sum(theta_grid * phi_est * diff(theta_grid)[1]) var_theta <- sum((theta_grid - mean_theta)^2 * phi_est * diff(theta_grid)[1]) # 5. Compute Momentum (Fourier) Variance # -------------------------------------- # Wavefunction: psi(theta) = sqrt(phi_est) psi <- sqrt(phi_est) # Fourier transform to momentum space psi_fft <- fft(psi) k <- 2 * pi * seq(-0.5, 0.5, length.out=1000) # Momentum grid # Compute |psi_fft|^2 (momentum probability density) momentum_prob <- abs(psi_fft)^2 momentum_prob <- momentum_prob / sum(momentum_prob) # Normalize # Momentum variance mean_k <- sum(k * momentum_prob) var_k <- sum((k - mean_k)^2 * momentum_prob) # 6. Validate Uncertainty Principle # --------------------------------- # Theoretical minimum: 0.5 (for [Θ, P] = i) uncertainty_product <- var_theta * var_k cat(sprintf("Uncertainty Product: %.3f (≥ 0.5? %s)\n", uncertainty_product, ifelse(uncertainty_product >= 0.5, "Yes", "No"))) # 7. Reconstruct Phi Using Gaussian Assumption (Example) # ------------------------------------------------------ # Assume Gaussian from mean_theta and var_theta (for illustration) phi_gaussian <- dnorm(theta_grid, mean_theta, sqrt(var_theta)) phi_gaussian <- phi_gaussian / sum(phi_gaussian * diff(theta_grid)[1]) # 8. Plot Results # --------------- df <- data.frame( theta = theta_grid, True = phi_true, KDE = phi_est, Gaussian_Fit = phi_gaussian ) # #### PLOT # ggplot(df, aes(x = theta)) + # geom_line(aes(y = True, color = "True Distribution"), linewidth=1) + # geom_line(aes(y = KDE, color = "KDE Estimate"), linetype="dashed") + # geom_line(aes(y = Gaussian_Fit, color = "Gaussian Fit"), linetype="dotted") + # labs(title = "Kime-Phase Distribution Estimation", # subtitle = sprintf("Uncertainty Product: %.2f", uncertainty_product), # x = expression(theta), y = "Density") + # scale_color_manual(values = c("True Distribution" = "blue", # "KDE Estimate" = "red", # "Gaussian Fit" = "green")) + # theme_minimal() library(plotly) # Create plot plot_ly(df, x = ~theta) %>% # True Distribution line add_trace( y = ~True, name = "True Distribution", type = "scatter", mode = "lines", line = list( color = "blue", width = 2 ) ) %>% # KDE Estimate line add_trace( y = ~KDE, name = "KDE Estimate", type = "scatter", mode = "lines", line = list( color = "red", width = 2, dash = "dash" ) ) %>% # Gaussian Fit line add_trace( y = ~Gaussian_Fit, name = "Gaussian Fit", type = "scatter", mode = "lines", line = list( color = "green", width = 2, dash = "dot" ) ) %>% # Layout configuration layout( title = list( text = paste0( "Time-Phase Distribution Estimation
", "Uncertainty Product: ", sprintf("%.2f", uncertainty_product), "" ), font = list(size = 20) ), xaxis = list( title = "θ", titlefont = list(size = 14), gridcolor = "lightgray", zerolinecolor = "gray" ), yaxis = list( title = "Density", titlefont = list(size = 14), gridcolor = "lightgray", zerolinecolor = "gray" ), legend = list( x = 1.1, y = 0.9, bordercolor = "lightgray", borderwidth = 1 ), paper_bgcolor = "white", plot_bgcolor = "white", margin = list(t = 100), showlegend = TRUE, hovermode = "x unified" ) %>% # Add modebar buttons config( modeBarButtonsToAdd = list( "hoverclosest", "hovercompare" ), displaylogo = FALSE ) ``` The *position operator* ($\Theta$) directly computes the mean and variance of $\theta$ and the *momentum operator* ($P$) uses the Fourier transform to compute momentum variance, avoiding numerical differentiation. The *uncertainty principle validation* is realized by computing the product of variances $\Delta \Theta \Delta P$ and checking against the theoretical lower bound of $0.5$. The kime-phase *distribution reconstruction* uses kernel density estimation (KDE) to approximate $\Phi$ from samples. the simulation includes a Gaussian fit based on position variance for comparison (simplified example). The resulting plot shows the *true distribution*, *KDE estimate*, and a *Gaussian fit* constrained by the position variance. The uncertainty product is printed, validating the non-commutativity of $\Theta$ and $P$. This (revised) *Approach* $1A$ more rigorously embeds the phase estimation within a Hilbert space framework, leveraging Fourier duality to implement momentum measurements and explicitly testing quantum-inspired uncertainty relations. In this reformulation, $\hat{F}_\alpha$ is *convolution* with $f_\alpha(\theta)$ and captures local phase structure, whereas $\hat{G}_\beta$ is *multiplication* by $g_\beta(\theta)$ and captures global oscillatory behavior. Alternatively, we can use $\hat{\Theta}$ and $\hat{D}_\beta = -i \beta \frac{d}{d\theta}$ for a position-momentum analogy. Non-commutativity ensures that measuring $\hat{F}_\alpha$ disturbs $\hat{G}_\beta$ outcomes, as standard QM. For example, $\hat{F}_\alpha$ emphasizes local peaks in $\Phi$, while $\hat{G}_\beta$ resolves periodic modes. The uncertainty relation $\Delta F_\alpha \Delta G_\beta \geq \frac{1}{2} |\langle [F_\alpha, G_\beta] \rangle|$ can be numerically validated, replacing the heuristic non-orthogonality condition. ### The Meaning of $\frac{d}{d\theta}$ in Kime How do we define a derivative operator when $\theta$ is a random variable rather than a deterministic coordinate? Do we need Itô calculus, distributional derivatives, or another framework to clarify the role of $\theta$ and $\frac{d}{d\theta}$ and explicate the derivative operator? In QM, $\hat{x}$ and $\hat{p}$ operate on wavefunctions $\psi(x)$ in $L^2(\mathbb{R})$, where $x$ is a continuous spatial variable. The derivative $\frac{d}{dx}$ is well-defined as an operator on differentiable functions. In the kime context, $\kappa = t e^{i\theta}$ introduces $\theta$ as a phase, randomly sampled from $\Phi(\theta; t)$ for each measurement at time $t$. *Approach 1 probes* $\Phi(\theta; t)$ via *test functions*, e.g., $f_\alpha(\theta)$, computing expectations like $\langle \Phi, f_\alpha \rangle = \int f_\alpha(\theta) \Phi(\theta; t) d\theta$. In Approach 1A reformulation, we defined operators on a Hilbert space, e.g., $L^2[-\pi, \pi]$, acting on functions $\psi(\theta)$, where $\psi(\theta)$ represent $\sqrt{\Phi(\theta; t)}$ or another state-like object. When $\theta$ is a random draw, $\frac{d}{d\theta}$ *seems* ill-defined without a continuous trajectory. Yet, $\Phi(\theta; t)$ is a density over $\theta$, suggesting a functional space where $\theta$ is a variable of integration, not a stochastic process evolving in time (as in Itô calculus). Since $\Phi(\theta; t)$ is a probability density over $\theta \in [-\pi, \pi]$ at fixed $t$, we treat $\theta$ as a coordinate in a function space, not a random variable evolving dynamically. Thus, - *Hilbert Space*: $L^2[-\pi, \pi]$, where functions $\psi(\theta)$ are square-integrable over $\theta$, and $\Phi(\theta; t) = |\psi(\theta)|^2$ (analogous to a probability density in QM). - *State*: $\psi(\theta)$ could be $\sqrt{\Phi(\theta; t)}$ (assuming $\Phi$ is real and non-negative), representing the "state" of the kime-phase distribution at time $t$. - *Operators*: Define $\hat{F}_\alpha$ and $\hat{G}_\beta$ as operators on $\psi(\theta)$, with expectations computed against $\Phi(\theta; t)$. The *differentiation operator* $\hat{G}_\beta = -i \frac{d}{d\theta}$ in $L^2[-\pi, \pi]$, where $\theta$ is the independent variable of the function space, not a random draw per se. The random sampling occurs when we generate data (e.g., $\theta_n(t) \sim \Phi$), but the operators act on the density or its square root. Thus, the *multiplication operator* $\hat{\Theta} \psi(\theta) = \theta \psi(\theta)$ and the *differentiation operator* $\hat{G}_\beta \psi(\theta) = -i \beta \frac{d}{d\theta} \psi(\theta)$, where $\beta$ scales the "momentum" (frequency-like) contribution. *Boundary Conditions*: Since $\theta \in [-\pi, \pi]$ is a circular domain (phase wraps around), we impose periodic boundary conditions: $\psi(-\pi) = \psi(\pi)$, ensuring $\hat{G}_\beta$ is well-defined on differentiable functions in $L^2[-\pi, \pi]$. *Commutator*:$[\hat{\Theta}, \hat{G}_\beta] \psi(\theta) = -i \beta \theta \frac{d}{d\theta} \psi(\theta) - (-i \beta \psi(\theta) - i \beta \theta \frac{d}{d\theta} \psi(\theta)) = i \beta \psi(\theta)$, where - $\hat{\Theta} \hat{G}_\beta \psi(\theta) = \theta (-i \beta \frac{d}{d\theta} \psi(\theta)) = -i \beta \theta \frac{d}{d\theta} \psi(\theta)$, and - $\hat{G}_\beta \hat{\Theta} \psi(\theta) = -i \beta \frac{d}{d\theta} (\theta \psi(\theta)) = -i \beta \left( \psi(\theta) + \theta \frac{d}{d\theta} \psi(\theta) \right)$. Thus, $[\hat{\Theta}, \hat{G}_\beta] = i \beta I$ confirms non-commutativity, with a commutator proportional to $\beta$. *Itô calculus does not apply* to stochastic processes where a variable, $\theta(t)$, evolves randomly over time $t$, governed by a stochastic differential equation (SDE). In kime, $\theta$ is sampled independently at each $t$ from $\Phi(\theta; t)$, with no explicit temporal dynamics (e.g., $d\theta = \mu dt + \sigma dW_t$). And the derivative $\frac{d}{d\theta}$ operates over the $\theta$-domain of the density $\Phi(\theta; t)$ at fixed $t$, not over $t$. Thus, Itô calculus is unnecessary here. The randomness is in the sampling process, not in a continuous evolution of $\theta$. A *distributional derivative* (e.g., weak derivative) would be needed if $\psi(\theta)$ or $\Phi(\theta; t)$ were not differentiable (e.g., included Dirac deltas). For simplicity, assume $\Phi(\theta; t)$ is smooth (e.g., *von Mises*), so the standard derivative suffices. If $\Phi$ were less regular, we’d define $\hat{G}_\beta$ in the sense of distributions, $\langle \hat{G}_\beta \psi, \phi \rangle = -i \beta \int_{-\pi}^\pi \psi(\theta) \frac{d}{d\theta} \phi(\theta) d\theta$, for test functions $\phi \in C^\infty[-\pi, \pi]$. Hence, the *reformulated Approach 1 with non-commuting operators* is based on 1. *Hilbert Space*: $L^2[-\pi, \pi]$, with $\psi(\theta) = \sqrt{\Phi(\theta; t)}$. 2. *Operators*: - *Multiplication operator*: $\hat{F}_\alpha = \hat{\Theta}$, probing the "position" (mean) of $\theta$, and - *Differentiation operator*: $\hat{G}_\beta = -i \beta \frac{d}{d\theta}$, probing frequency-like variations. 3. *Measurement*: In the state, $\psi(\theta) = \sqrt{\Phi(\theta; t)}$, we compute the expectations $\langle \psi | \hat{F}_\alpha | \psi \rangle = \int \theta \Phi(\theta; t) d\theta$, $\langle \psi | \hat{G}_\beta | \psi \rangle = -i \beta \int \psi(\theta) \frac{d}{d\theta} \psi(\theta) d\theta$. 4. *Kime-Phase Recovery/Reconstruction*: Use multiple $\beta$ values to estimate $\Phi(\theta; t)$ via moments or Fourier components. 5. *Expectations*: $\langle \hat{\Theta} \rangle = \int_{-\pi}^\pi \theta \Phi(\theta; t) d\theta$ and $$\langle \hat{G}_\beta \rangle = -i \beta \int_{-\pi}^\pi \sqrt{\Phi(\theta; t)} \frac{d}{d\theta} \sqrt{\Phi(\theta; t)} d\theta = -i \beta \int_{-\pi}^\pi \frac{1}{2} \frac{\frac{d}{d\theta} \Phi(\theta; t)}{\sqrt{\Phi(\theta; t)}} d\theta.$$ Mind that the statement, "*no explicit temporal dynamics*" refers to the absence of a stochastic differential equation (SDE) or a continuous-time process governing $\theta$’s evolution over $t$. For example, in a stochastic process like a Wiener process or an Itô diffusion, $\theta(t)$ would evolve according to something like $d\theta = \mu dt + \sigma dW_t$, where $W_t$ is Brownian motion, implying a smooth, temporally correlated trajectory with increments dependent on previous values. Such a model introduces *temporal dependence* through a differential structure, where $\theta(t + dt)$ depends on $\theta(t)$ via drift ($\mu$) and diffusion ($\sigma$). In contrast, the Approach 1A kime framework treats $\theta$ as follows - At each fixed time $t$, a value $\theta_n(t)$ is sampled independently from the distribution $\Phi(\theta; t)$ for each measurement $n = 1, 2, \ldots, N$. - The next time point $t' = t + \Delta t$ yields a new sample $\theta_n(t')$ from $\Phi(\theta; t')$, with no explicit rule (like an SDE) linking $\theta_n(t)$ to $\theta_n(t')$. The "no explicit temporal dynamics" claim hinges on the idea that $\theta_n(t)$ and $\theta_n(t')$ are *independent (IID) draws*, lacking a mechanistic model of temporal evolution, e.g., no Markovian transition or differential equation. However, $\Phi(\theta; t)$ is time-dependent, meaning the distribution from which $\theta$ is sampled changes with $t$. For instance, in the simulation, $\mu(t) = 2\pi \sin(0.1 t)$ and $\kappa(t) = 5 (1 + 0.5 \sin(2\pi t / T))$ define a von Mises distribution $\Phi(\theta; t) = \text{vM}(\mu(t), \kappa(t))$, where both *mean* and *concentration* parameters vary over time. This implies that the statistical properties of $\theta$, e.g., mean, variance, evolve with $t$, even if the samples at different $t$ are drawn independently. This time-dependence of $\Phi(\theta; t)$ introduces a subtlety; while there’s no *dynamic process* linking $\theta(t)$ to $\theta(t')$, like an SDE, the distribution’s parameters are functions of $t$, suggesting an *implicit* temporal structure through the changing $\Phi$. *Statistical independence across time*: For a fixed measurement $n$, $\theta_n(t)$ and $\theta_n(t')$ are independent draws from $\Phi(\theta; t)$ and $\Phi(\theta; t')$, respectively. There’s no correlation imposed between $\theta_n(t)$ and $\theta_n(t')$ by a temporal process—just new, independent samples at each $t$. This supports the "no explicit temporal dynamics" claim. At the same time, $\Phi(\theta; t)$ varies with $t$ and the ensemble of samples $\{\theta_n(t)\}_{n=1}^N$ at time $t$ has different statistical characteristics, e.g., mean $\mu(t)$, concentration $\kappa(t)$, than at $t'$. This introduces a form of *temporal dependence* in the distribution, *not in the individual* $\theta_n(t)$ trajectories. There’s no continuous stochastic process or differential equation defining how $\theta$ evolves from $t$ to $t + dt$. Each $\theta_n(t)$ is a fresh draw, not a step in a trajectory. The independence of samples across $t$ means no explicit temporal correlation, e.g., autocorrelation, is modeled between $\theta_n(t)$ and $\theta_n(t')$. Hence, in kime, $\theta$ is sampled independently at each $t$ from a time-dependent distribution $\Phi(\theta; t)$, with no explicit stochastic temporal dynamics, e.g., $d\theta = \mu dt + \sigma dW_t$, linking successive samples. However, the distribution $\Phi(\theta; t)$ itself varies with $t$, introducing an implicit temporal dependence in the statistical properties of $\theta$. Let's also clarify the kime-operators and $\frac{d}{d\theta}$ in the reformulation of Approach 1A, The derivative $\hat{G}_\beta = -i \beta \frac{d}{d\theta}$ operates on functions $\psi(\theta) = \sqrt{\Phi(\theta; t)}$ in $L^2[-\pi, \pi]$ at a fixed $t$. Here, $\theta$ is a coordinate in the phase domain, not a time-evolving variable. The operator probes the spatial (phase) structure of $\Phi(\theta; t)$, not its temporal evolution. Below is an updated simulation with time-dependent $\Phi$ refining the prior simulation to explicitly express the $\Phi(\theta; t)$’s time-dependence. ```{r approach1A_Part3, warning=FALSE, message=FALSE} library(circular) library(plotly) set.seed(123) # Parameters T <- 100 N <- 1000 t <- seq(0, 10, length.out = T) theta_grid <- seq(-pi, pi, length.out = 1000) dtheta <- theta_grid[2] - theta_grid[1] # True time-dependent von Mises distribution mu_true <- 2 * pi * sin(0.1 * t) kappa_true <- 5 * (1 + 0.5 * sin(2 * pi * t / T)) # Operators Theta_op <- function(psi, theta) theta * psi G_beta_op <- function(psi, theta, beta = 1) { -1i * beta * diff(c(psi[length(psi)], psi))[1:length(theta)] / dtheta } # Forward model generate_signal <- function(t, theta) sin(2 * pi * 0.5 * t + theta) + rnorm(length(t), 0, 0.1) # Simulate data theta_true <- matrix(NA, N, T) x <- matrix(NA, N, T) for (n in 1:N) { for (tt in 1:T) { theta_true[n, tt] <- rvonmises(1, mu = mu_true[tt], kappa = kappa_true[tt]) } x[n, ] <- generate_signal(t, theta_true[n, ]) } # Manual Hilbert transform hilbert_manual <- function(x) { n <- length(x) fft_x <- fft(x) h <- numeric(n); h[1] <- 1; h[2:(n/2)] <- 2; h[(n/2 + 1):n] <- 0 fft(h * fft_x, inverse = TRUE) / n } # Phase inference theta_est <- t(apply(x, 1, function(row) Arg(hilbert_manual(row)))) # Moment estimation m1_est <- colMeans(exp(1i * theta_est)) mu_hat <- Arg(m1_est) R <- abs(m1_est) kappa_hat <- sapply(R, function(r) if (r >= 0.999) 100 else if (r <= 0.001) 0 else uniroot(function(k) besselI(k, 1) / besselI(k, 0) - r, c(0, 100))$root) # Operator expectations at sample t Phi_t <- dvonmises(theta_grid, mu = mu_true[50], kappa = kappa_true[50]) Phi_t <- Phi_t / sum(Phi_t * dtheta) psi_t <- sqrt(Phi_t) exp_Theta <- sum(Theta_op(psi_t, theta_grid) * psi_t) * dtheta exp_G <- sum(G_beta_op(psi_t, theta_grid, 1) * psi_t) * dtheta cat("At t[50]: =", exp_Theta, ", =", exp_G, "\n") # Plot df <- data.frame(t = t, mu_true = mu_true %% (2 * pi) - pi, kappa_true = kappa_true, mu_hat = mu_hat %% (2 * pi) - pi, kappa_hat = kappa_hat) fig_mu <- plot_ly(df) %>% add_trace(x = ~t, y = ~mu_true, type = "scatter", mode = "lines", name = "True μ") %>% add_trace(x = ~t, y = ~mu_hat, type = "scatter", mode = "lines", name = "Est μ") %>% layout(title = "Mean Phase", yaxis = list(title = "μ(t)")) fig_kappa <- plot_ly(df) %>% add_trace(x = ~t, y = ~kappa_true, type = "scatter", mode = "lines", name = "True κ") %>% add_trace(x = ~t, y = ~kappa_hat, type = "scatter", mode = "lines", name = "Est κ") %>% layout(title = "Concentration", yaxis = list(title = "κ(t)")) fig_mu fig_kappa ``` This reflects $\Phi(\theta; t)$’s time-dependence while maintaining independent sampling, aligning with the Approach 1A refined formulation. ## Approach 1B: KPT via Test Function Action KPT probes \(\Phi(\theta; t)\) by evaluating its action on families of test functions, analogous to measuring a quantum state in different bases. In quantum mechanics, a state \(|\psi\rangle\) is measured in bases like position (\(|x\rangle\)) and momentum (\(|p\rangle\)), with non-commuting operators \(\hat{X}\) and \(\hat{P}\), yielding distributions \(|\langle x | \psi \rangle|^2\) and \(|\langle p | \psi \rangle|^2\). These measurements constrain the wavefunction, including its phase, but require repeated trials due to the uncertainty principle. In the kime framework, we treat \(\Phi(\theta; t)\) as a distribution on \([-\pi, \pi)\), with \(\int_{-\pi}^{\pi} \Phi(\theta; t) d\theta = 1\). We probe \(\Phi\) via its action on test functions \(f(\theta)\): \[ \langle \Phi, f \rangle = \int_{-\pi}^{\pi} f(\theta) \Phi(\theta; t) d\theta, \] which is analogous to a quantum expectation value \(\langle \hat{A} \rangle_\psi\). To capture complementary information about \(\Phi\), we use two families of test functions - {\it{Family A (Fourier-like)}}: \(f_\alpha(\theta) = sin(\alpha \theta), \(\alpha \in \{1, 2, 3,4,5\}\), modified to be periodic on \([-\pi, \pi)\). - {\it{Family B (Fourier-like)}}: \(g_\beta(\theta) = \cos(\beta \theta)\), \(\beta \in \{1, 2, 3,4,5\}\), naturally periodic. These families are chosen to capture local (wavelet) and global (Fourier) features of \(\Phi\). While quantum tomography relies on non-commuting operators, here we define “incompatibility” operationally: measurements with \(f_\alpha\) and \(g_\beta\) are performed on separate samples, mimicking the need for separate runs in quantum mechanics due to non-commutativity. Given \(N\) repeated measurements at time \(t\), each yielding a sample \(\theta_i \sim \Phi(\theta; t)\), we estimate the integrals: \[ \mathbb{E}[f_\alpha(\theta)] = \int_{-\pi}^{\pi} f_\alpha(\theta) \Phi(\theta; t) d\theta \approx \frac{1}{N_f} \sum_{i=1}^{N_f} f_\alpha(\theta_i), \] \[ \mathbb{E}[g_\beta(\theta)] = \int_{-\pi}^{\pi} g_\beta(\theta) \Phi(\theta; t) d\theta \approx \frac{1}{N_g} \sum_{i=1}^{N_g} g_\beta(\theta_i), \] where \(N_f + N_g = N\), and the samples are split between the two families (e.g., \(N_f = N_g = N/2\)). We reconstruct \(\Phi\) by parameterizing: \[ \Phi_{\text{est}}(\theta) = \sum_{\alpha} c_\alpha f_\alpha(\theta) + \sum_{\beta} d_\beta g_\beta(\theta), \] and fitting \(c_\alpha, d_\beta\) to match the measured integrals. To ensure \(\Phi_{\text{est}} \geq 0\), we use a constrained optimization: \[ \min_{c_\alpha, d_\beta \geq 0} \sum_{\alpha} \left( \frac{1}{N_f} \sum_{i=1}^{N_f} f_\alpha(\theta_i) - \int_{-\pi}^{\pi} f_\alpha(\theta) \Phi_{\text{est}}(\theta) d\theta \right)^2 + \sum_{\beta} \left( \frac{1}{N_g} \sum_{i=1}^{N_g} g_\beta(\theta_i) - \int_{-\pi}^{\pi} g_\beta(\theta) \Phi_{\text{est}}(\theta) d\theta \right)^2 + \lambda \left( \sum_{\alpha} c_\alpha^2 + \sum_{\beta} d_\beta^2 \right), \] where \(\lambda > 0\) is a Tikhonov penalty to prevent overfitting. After optimization, we normalize: \[ \Phi_{\text{est}}(\theta) \leftarrow \frac{\Phi_{\text{est}}(\theta)}{\int_{-\pi}^{\pi} \Phi_{\text{est}}(\theta) d\theta}. \] This ensures \(\Phi_{\text{est}}\) is a valid probability density. To formalize incompatibility, we define operators on \(L^2([-\pi, \pi])\): - \(\hat{F}_\alpha \psi(\theta) = f_\alpha(\theta) \psi(\theta)\), a multiplication operator. - \(\hat{G}_\beta \psi(\theta) = -i \beta \frac{d}{d\theta} (g_\beta(\theta) \psi(\theta))\), a differential operator. The commutator \([\hat{F}_\alpha, \hat{G}_\beta] \neq 0\), as: \[ [\hat{F}_\alpha, \hat{G}_\beta] \psi = f_\alpha (-i \beta) \frac{d}{d\theta} (g_\beta \psi) + i \beta \frac{d}{d\theta} (f_\alpha g_\beta \psi) \neq 0, \] due to the derivative terms. This non-commutativity justifies separate measurements, as in quantum tomography. A simulation demonstrates KPT with \(N = 200\) samples at a fixed \(t\), where the true \(\Phi(\theta; t)\) is a mixture of two von Mises distributions: \[ \Phi(\theta; t) = 0.6 vM(\theta; 0, 5) + 0.4 vM(\theta; \pi, 5). \] We use \(N_f = N_g = 100\), with \(\alpha, \beta \in \{1, 2, 3\}\), and fit \(\Phi_{\text{est}}\) using the above optimization. Figure \ref{fig:Approach1.Fig1} shows that \(\Phi_{\text{est}}\) (red) approximates the true \(\Phi\) (blue), capturing the bimodal structure, though with some ringing due to the limited basis size. For practical applications like fMRI, we propose extracting \(\theta\) from oscillatory signals (e.g., via Hilbert transform) and applying wavelet or Fourier transforms to compute test function integrals. Future work should address the time-dependence of \(\Phi(\theta; t)\) and optimize the choice of test functions (e.g., using B-splines for smoother reconstructions). ```{r approach1B_TestFuncAct, echo=TRUE, message=FALSE, warning=FALSE} # Load required libraries library(circular) # For von Mises distribution library(plotly) # For 3D surface plotting library(htmlwidgets) # For saving as SVG library(nloptr) # For constrained optimization # Set seed for reproducibility set.seed(123) # Von Mises mixture density function von_mises_mixture <- function(theta, w, mu1, mu2, k1, k2) { # vm1 <- dvonmises(theta, mu1, k1) / (2 * pi * besselI(k1, 0)) # vm2 <- dvonmises(theta, mu2, k2) / (2 * pi * besselI(k2, 0)) vm1 <- dvonmises(theta, mu1, k1) vm2 <- dvonmises(theta, mu2, k2) w * vm1 + (1 - w) * vm2 } # Define time-dependent parameters for Phi(theta; t) w_t <- function(t) 0.5 + 0.3 * sin(2 * pi * t / 50) mu1_t <- function(t) (0.5 * t) %% (2 * pi) - pi mu2_t <- function(t) (0.5 * t + pi) %% (2 * pi) - pi k1_t <- function(t) 5 + 2 * sin(2 * pi * t / 30) k2_t <- function(t) 5 - 2 * sin(2 * pi * t / 30) # True Phi(theta; t) phi_true <- function(theta, t) { w <- w_t(t) mu1 <- mu1_t(t) mu2 <- mu2_t(t) k1 <- k1_t(t) k2 <- k2_t(t) von_mises_mixture(theta, w, mu1, mu2, k1, k2) } # Test function families # Fourier-like (sine functions) f_alpha <- function(theta, alpha, sigma = 0.5) { sin(alpha * theta) } # Fourier-like (cosine functions) g_beta <- function(theta, beta) { cos(beta * theta) } # Generate test functions for a range of alpha and beta alphas <- 1:5 betas <- 1:5 test_functions <- list( f = lapply(alphas, function(a) function(theta) f_alpha(theta, a)), g = lapply(betas, function(b) function(theta) g_beta(theta, b)) ) # Simulate measurements and reconstruct Phi simulate_kpt <- function(t_values, theta_grid, N = 500, lambda = 0.01) { N_f <- N_g <- N / 2 # Split samples between families phi_true_vals <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) phi_est_vals <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) mse <- numeric(length(t_values)) for (i in seq_along(t_values)) { t <- t_values[i] # True Phi at time t phi_true_vals[i, ] <- phi_true(theta_grid, t) # Simulate samples theta_i ~ Phi(theta; t) w <- w_t(t) mu1 <- mu1_t(t) mu2 <- mu2_t(t) k1 <- k1_t(t) k2 <- k2_t(t) theta_samples <- c( rvonmises(N_f, mu1, k1), # Samples for f_alpha rvonmises(N_g, mu2, k2) # Samples for g_beta ) # Adjust samples based on mixture weight mix <- runif(N) theta_samples <- ifelse(mix < w, theta_samples[1:N_f], theta_samples[(N_f+1):N]) # Compute sample averages for each test function f_means <- sapply(test_functions$f, function(f) mean(sapply(theta_samples[1:N_f], f))) g_means <- sapply(test_functions$g, function(g) mean(sapply(theta_samples[(N_f+1):N], g))) # Define the objective function for optimization objective <- function(coeffs) { c_alpha <- coeffs[1:length(alphas)] d_beta <- coeffs[(length(alphas)+1):length(coeffs)] phi_est <- function(theta) { sum_f <- sum(sapply(seq_along(alphas), function(a) c_alpha[a] * f_alpha(theta, alphas[a]))) sum_g <- sum(sapply(seq_along(betas), function(b) d_beta[b] * g_beta(theta, betas[b]))) sum_f + sum_g } # Compute integrals numerically f_integrals <- sapply(test_functions$f, function(f) { mean(sapply(theta_grid, function(theta) f(theta) * phi_est(theta))) }) g_integrals <- sapply(test_functions$g, function(g) { mean(sapply(theta_grid, function(theta) g(theta) * phi_est(theta))) }) # Sum of squared residuals + Tikhonov penalty sum((f_means - f_integrals)^2) + sum((g_means - g_integrals)^2) + lambda * (sum(c_alpha^2) + sum(d_beta^2)) } # Constraints: coefficients >= 0 opts <- list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1e-6, "maxeval" = 1000) init_coeffs <- rep(0.1, length(alphas) + length(betas)) res <- nloptr( x0 = init_coeffs, eval_f = objective, lb = rep(0, length(init_coeffs)), opts = opts ) # Reconstruct Phi_est coeffs <- res$solution c_alpha <- coeffs[1:length(alphas)] d_beta <- coeffs[(length(alphas)+1):length(coeffs)] phi_est <- function(theta) { sum_f <- sum(sapply(seq_along(alphas), function(a) c_alpha[a] * f_alpha(theta, alphas[a]))) sum_g <- sum(sapply(seq_along(betas), function(b) d_beta[b] * g_beta(theta, betas[b]))) sum_f + sum_g } phi_est_vals[i, ] <- sapply(theta_grid, phi_est) # Ensure positivity and normalize phi_est_vals[i, ] <- pmax(phi_est_vals[i, ], 0) # integral <- mean(phi_est_vals[i, ]) * (2 * pi) / length(theta_grid) # Assume equidistant theta_grid integral <- sum(phi_est_vals[i, ])*(theta_grid[2] - theta_grid[1]) phi_est_vals[i, ] <- phi_est_vals[i, ] / integral # Compute MSE mse[i] <- mean((phi_true_vals[i, ] - phi_est_vals[i, ])^2) } return(list(phi_true = phi_true_vals, phi_est = phi_est_vals, mse = mse)) } # simulate_kpt <- function(t_values, theta_grid, N = 500, lambda = 0.01) { # # Ensure N is even for simplicity, or add handling for odd N if needed # if (N %% 2 != 0) { # warning("N is odd, adjusting N to be even: ", N + 1) # N <- N + 1 # } # N_f <- N_g <- N / 2 # # # --- Add checks for required objects --- # required_funcs <- c("f_alpha", "g_beta", "w_t", "mu1_t", "mu2_t", "k1_t", "k2_t", "phi_true") # missing_funcs <- required_funcs[!sapply(required_funcs, exists, mode = "function")] # if (length(missing_funcs) > 0) stop("Missing required functions: ", paste(missing_funcs, collapse=", ")) # # required_objs <- c("alphas", "betas", "test_functions") # missing_objs <- required_objs[!sapply(required_objs, exists)] # if (length(missing_objs) > 0) stop("Missing required objects: ", paste(missing_objs, collapse=", ")) # # if (!is.list(test_functions) || !all(c("f", "g") %in% names(test_functions))) stop("test_functions must be a list with names 'f' and 'g'") # if (!is.vector(alphas) || !is.numeric(alphas)) stop("alphas must be a numeric vector") # if (!is.vector(betas) || !is.numeric(betas)) stop("betas must be a numeric vector") # if (length(test_functions$f) != length(alphas)) stop("Length of test_functions$f must match length of alphas") # if (length(test_functions$g) != length(betas)) stop("Length of test_functions$g must match length of betas") # # --- End checks --- # # # phi_true_vals <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) # phi_est_vals <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) # mse <- numeric(length(t_values)) # # for (i in seq_along(t_values)) { # t <- t_values[i] # cat("\nProcessing t =", t, "(iteration", i, "of", length(t_values), ")\n") # DEBUG PROGRESS # # # True Phi at time t # phi_true_vals[i, ] <- phi_true(theta_grid, t) # # # Simulate samples theta_i ~ Phi(theta; t) # # w <- w_t(t) # Not strictly needed for sampling if calculating means separately # mu1 <- mu1_t(t) # mu2 <- mu2_t(t) # k1 <- k1_t(t) # k2 <- k2_t(t) # # # --- DEBUG: Check mu and kappa values --- # cat(" Parameters: mu1=", mu1, " k1=", k1, " mu2=", mu2, " k2=", k2, "\n") # if (!all(sapply(c(mu1, k1, mu2, k2), is.numeric)) || any(is.na(c(mu1, k1, mu2, k2)))) { # stop("Non-numeric or NA mu/kappa values detected at t=", t) # } # if (k1 < 0 || k2 < 0) { # stop("Negative kappa (k) detected at t=", t) # } # # --- End DEBUG --- # # # Generate samples # theta_samples_f <- rvonmises(N_f, mu = mu1, kappa = k1) # theta_samples_g <- rvonmises(N_g, mu = mu2, kappa = k2) # # # --- DEBUG: Check generated samples --- # if (any(!is.numeric(theta_samples_f)) || any(is.na(theta_samples_f))) { # warning("Non-numeric or NA samples generated for f at t=", t) # print(summary(theta_samples_f)) # stop("Stopping due to invalid samples for f") # } # if (any(!is.numeric(theta_samples_g)) || any(is.na(theta_samples_g))) { # warning("Non-numeric or NA samples generated for g at t=", t) # print(summary(theta_samples_g)) # stop("Stopping due to invalid samples for g") # } # cat(" Samples generated successfully.\n") # # --- End DEBUG --- # # # # Compute sample averages # cat(" Calculating f_means...\n") # f_means <- sapply(seq_along(test_functions$f), function(idx) { # f <- test_functions$f[[idx]] # current_alpha <- alphas[idx] # if (!is.numeric(current_alpha) || is.na(current_alpha)) { # stop(paste("Non-numeric or NA alpha detected at index", idx, "value:", current_alpha)) # } # # # --- DEBUG: Check inside inner sapply for f_means --- # inner_results <- sapply(theta_samples_f, function(theta_val) { # if (!is.numeric(theta_val) || is.na(theta_val)) { # stop("Non-numeric or NA theta value passed to test function f: ", theta_val) # } # # Use tryCatch to pinpoint the exact failing call # res <- tryCatch({ # f(theta_val) # Apply the test function # }, error = function(e) { # cat(" ERROR occurred in test_functions$f[[", idx, "]] (alpha=", current_alpha, ") with theta=", theta_val, "\n", sep="") # print(e) # Print the specific error (should be the non-numeric arg error) # stop("Error during f(theta) execution.") # Stop execution here to investigate # }) # if (!is.numeric(res) || is.na(res)) { # cat(" WARNING: Non-numeric or NA result from test_functions$f[[", idx, "]] (alpha=", current_alpha, ") for theta=", theta_val, " -> result:", res, "\n", sep="") # # Decide whether to stop or allow NAs # # stop("Non-numeric result from test function f") # return(NA) # Allow NA propagation for now # } # return(res) # }) # # --- End DEBUG --- # mean(inner_results, na.rm = TRUE) # }) # cat(" f_means calculated.\n") # # # cat(" Calculating g_means...\n") # g_means <- sapply(seq_along(test_functions$g), function(idx) { # g <- test_functions$g[[idx]] # current_beta <- betas[idx] # if (!is.numeric(current_beta) || is.na(current_beta)) { # stop(paste("Non-numeric or NA beta detected at index", idx, "value:", current_beta)) # } # # # --- DEBUG: Check inside inner sapply for g_means --- # inner_results <- sapply(theta_samples_g, function(theta_val) { # if (!is.numeric(theta_val) || is.na(theta_val)) { # stop("Non-numeric or NA theta value passed to test function g: ", theta_val) # } # res <- tryCatch({ # g(theta_val) # Apply the test function # }, error = function(e) { # cat(" ERROR occurred in test_functions$g[[", idx, "]] (beta=", current_beta, ") with theta=", theta_val, "\n", sep="") # print(e) # stop("Error during g(theta) execution.") # }) # if (!is.numeric(res) || is.na(res)) { # cat(" WARNING: Non-numeric or NA result from test_functions$g[[", idx, "]] (beta=", current_beta, ") for theta=", theta_val, " -> result:", res, "\n", sep="") # # stop("Non-numeric result from test function g") # return(NA) # } # return(res) # }) # # --- End DEBUG --- # mean(inner_results, na.rm = TRUE) # }) # cat(" g_means calculated.\n") # # # --- Rest of the function (Objective, nloptr, Reconstruction, MSE) --- # # Add similar checks inside the objective function if needed, especially around phi_est calls # # # Define the objective function for optimization # objective <- function(coeffs) { # # ... (objective function code as before) ... # # Add checks within phi_est if necessary # phi_est <- function(theta) { # # Ensure theta is numeric # if(!is.numeric(theta)) stop("Non-numeric theta passed to phi_est") # sum_f <- sum(sapply(seq_along(alphas), function(a_idx) { # if(!is.numeric(c_alpha[a_idx])) stop("Non-numeric c_alpha") # if(!is.numeric(alphas[a_idx])) stop("Non-numeric alpha in phi_est") # res_f <- f_alpha(theta, alphas[a_idx]) # if(!is.numeric(res_f)) stop("Non-numeric result from f_alpha in phi_est") # c_alpha[a_idx] * res_f # })) # sum_g <- sum(sapply(seq_along(betas), function(b_idx) { # if(!is.numeric(d_beta[b_idx])) stop("Non-numeric d_beta") # if(!is.numeric(betas[b_idx])) stop("Non-numeric beta in phi_est") # res_g <- g_beta(theta, betas[b_idx]) # if(!is.numeric(res_g)) stop("Non-numeric result from g_beta in phi_est") # d_beta[b_idx] * res_g # })) # if(!is.numeric(sum_f + sum_g)) stop("Non-numeric result from phi_est") # return(sum_f + sum_g) # } # # ... (rest of objective function) ... # } # # # ... (nloptr call, reconstruction, normalization, MSE) ... # cat(" Optimization and reconstruction complete for t =", t, "\n") # } # End of loop over t_values # # return(list(phi_true = phi_true_vals, phi_est = phi_est_vals, mse = mse)) # } # Run the simulation t_values <- seq(0, 100, by = 10) theta_grid <- seq(-pi, pi, length.out = 200) results <- simulate_kpt(t_values, theta_grid, N = 500, lambda = 0.01) # Quantitative summaries cat("Mean Squared Errors at each t:\n") for (i in seq_along(t_values)) { cat(sprintf("t = %.1f: MSE = %.6f\n", t_values[i], results$mse[i])) } cat(sprintf("Average MSE: %.6f\n", mean(results$mse))) # Plot the results using plotly # True Phase surface fig <- plot_ly() %>% add_surface(x = theta_grid, y = t_values, z = results$phi_true, colorscale = "Viridis", showscale = FALSE, name = "true Phi(theta; t)", scene = 'scene1') %>% layout(scene = list(aspectmode = "cube")) # Estimated Phi surface fig <- fig %>% add_surface(x = theta_grid, y = t_values, z = ( results$phi_est), # approx normalize the Phi_est type = "surface", colorscale = "Magma", name = "Estimated Phi(theta; t)", scene = 'scene2', showscale = FALSE) %>% layout(scene2 = list(aspectmode = "cube")) # Difference True - Est Phase surface fig <- fig %>% add_surface(x = theta_grid, y = t_values, z = (results$phi_true - (results$phi_est)), type = "surface", colorscale = "Magma", name = "Difference True - Estimated Phase", scene = 'scene3', showscale = FALSE) %>% layout(scene3 = list(aspectmode = "cube")) fig <- fig %>% layout(title = "True (left), Estimated (middle), and Difference Kime-Phase Surfaces", showlegend = TRUE, legend = list(x = 0.85, y = 0.1), annotations = list(list(showarrow = FALSE, text = '(true)', xref = 'paper', yref = 'paper', x = 0.15, y = 0.95, xanchor = 'center', yanchor = 'bottom', font = list(size = 16)), list(showarrow = FALSE, text = '(estimated)', xref = 'paper', yref = 'paper', x = 0.5, y = 0.95, xanchor = 'center', yanchor = 'bottom', font = list(size = 16)), list(showarrow = FALSE, text = '(difference)', xref = 'paper', yref = 'paper', x = 0.85, y = 0.95, xanchor = 'center', yanchor = 'bottom', font = list(size = 16)) ) ) # Add synchronization fig <- fig %>% onRender( "function(el) { el.on('plotly_relayout', function(d) { const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key)); if (camera.length) { const scenes = Object.keys(el.layout).filter((key) => /^scene\\d*/.test(key)); const new_layout = {}; scenes.forEach(key => { new_layout[key] = {...el.layout[key], camera: {...d[camera]}}; }); Plotly.relayout(el, new_layout); } }); }" ) fig ``` ## Approach 2: Kime-Phase Estimation via Non-Commuting Operators Approach 2 estimates the kime-phase \(\phi(t) = \theta(t)\) from a time-series signal \(s(t) = A(t) e^{i\phi(t)}\), where \(\phi(t) \sim \Phi(\theta; t)\), and subsequently reconstruct the kime-phase distribution \(\Phi(\theta; t)\). This aligns with the goal of Kime-Phase Tomography (KPT) to estimate \(\Phi(\theta; t)\), but unlike Approach 1, which operates in the \(\theta\)-domain, Approach 2 works directly with time-series signals in the time domain, making it more practical for applications like EEG or fMRI. This approach is based on - Hilbert Space: \(L^2([0, T])\), the space of square-integrable functions over a finite time interval \(t \in [0, T]\), with inner product $\langle f, g \rangle = \int_0^T f(t) \overline{g(t)} dt.$ This is appropriate for finite-duration signals in practical applications. - Signal Model: The signal is \(s(t) = A(t) e^{i\phi(t)}\), where \(A(t) \geq 0\) is the amplitude, and \(\phi(t) \in [-\pi, \pi)\) is the kime-phase, sampled from a time-dependent distribution \(\Phi(\theta; t)\). - Kime Operators: Three linear operators \(K_1, K_2, K_3\) that probe the signal in time, frequency, and scale domains, respectively. These operators are designed to be non-commuting, providing complementary information about \(\phi(t)\). 1. The Time-Domain Operator (\(K_1\)): Multiplication by time, analogous to the position operator in quantum mechanics $K_1 s(t) = t s(t).$ This operator emphasizes the temporal structure of the signal. 2. Frequency-Domain Operator (\(K_2\)): The derivative operator, analogous to the momentum operator $K_2 s(t) = -i \frac{d}{dt} s(t).$ This operator probes the frequency content of the signal, as the derivative of \(e^{i\omega t}\) yields \(-i \omega e^{i\omega t}\). 3. Scale-Domain Operator (\(K_3\)): Convolution with a Morlet wavelet at an optimal scale $K_3 s(t) = \int_0^T \psi(t - t'; a^*) s(t') dt',$ where \(\psi(t; a) = \frac{1}{\sqrt{a}} e^{-(t/a)^2/2} \cos\left(\omega_0 \frac{t}{a}\right) \cdot \mathbb{I}_{[-a, a]}(t)\), \(\omega_0 = 5\) (to balance time-frequency resolution), and \(\mathbb{I}_{[-a, a]}(t)\) ensures compact support to avoid boundary effects. The optimal scale $a^*$ is $a^* = \arg\max_a \left| \int_0^T \psi(t - t'; a) s(t') dt' \right|.$ This operator captures multi-scale features of the signal. To ensure complementary information, let's compute the pairwise commutators $$K_1 K_2 s(t) = t \left(-i \frac{d}{dt} s(t)\right),$$ $$K_2 K_1 s(t) = -i \frac{d}{dt} (t s(t)) = -i s(t) - i t \frac{d}{dt} s(t),$$ $$[K_1, K_2] s(t) = K_1 K_2 s(t) - K_2 K_1 s(t) = -i t \frac{d}{dt} s(t) + i s(t) + i t \frac{d}{dt} s(t) = i s(t),$$ Hence, $[K_1, K_2] = i I,$ where \(I\) is the identity operator. This mirrors the quantum mechanical relation \([\hat{x}, \hat{p}] = i \hbar\) (with \(\hbar = 1\)). The uncertainty relation is $\Delta K_1 \Delta K_2 \geq \frac{1}{2} |\langle [K_1, K_2] \rangle| = \frac{1}{2},$ where \(\Delta K_j = \sqrt{\langle K_j^2 \rangle - \langle K_j \rangle^2}\), and expectations are computed over the signal \(s(t)\), $\langle K_j \rangle = \int_0^T s(t)^* (K_j s(t)) dt.$ The commutotars of \([K_1, K_3]\) and \([K_2, K_3]\) are generally non-trivial due to the convolution operation. For example, $K_1 K_3 s(t) = t \int_0^T \psi(t - t'; a^*) s(t') dt',$ $K_3 K_1 s(t) = \int_0^T \psi(t - t'; a^*) (t' s(t')) dt',$ and $[K_1, K_3] s(t) = \int_0^T \psi(t - t'; a^*) (t - t') s(t') dt',$ which is non-zero since \(t \neq t'\) in general. Similarly, \([K_2, K_3] \neq 0\), reflecting time-scale and frequency-scale trade-offs. We’ll validate these numerically in the simulation. Phase estimation allows inferring the phase from each operator’s output - \(K_1\): \(\phi_1(t) = \arg(K_1 s(t)) = \arg(t s(t))\). Since \(t \geq 0\), this simplifies to \(\arg(s(t))\) adjusted for the sign of \(t\), but in practice, we use the analytic signal directly. - \(K_2\): \(\phi_2(t) = \arg(K_2 s(t)) = \arg\left(-i \frac{d}{dt} s(t)\right)\). - \(K_3\): \(\phi_3(t) = \arg(K_3 s(t))\), from the wavelet transform at scale $a^*$. KPT is based on combining the phase estimates using weights based on circular dispersion, which is more robust for phase data than variance $\text{Disp}(\phi_j(t)) = 1 - \left| \frac{1}{W} \sum_{k} e^{i \phi_j(t_k)} \right|,$ where the sum is over a time window \(W\) around \(t\). The weights are $w_j(t) = \frac{1}{\text{Disp}(\phi_j(t)) + \epsilon},$ with \(\epsilon = 10^{-6}\) to avoid division by zero. The ensemble phase is $\hat{\phi}(t) = \arg\left( \sum_{j=1}^3 w_j(t) e^{i\phi_j(t)} \right).$ Kime-phase distribution estimation, \(\Phi(\theta; t)\), from \(\hat{\phi}(t)\), uses kernel density estimation (KDE) with a von Mises kernel $$\Phi_{\text{est}}(\theta; t) = \frac{1}{W} \sum_{k} \frac{1}{2\pi I_0(\kappa)} e^{\kappa \cos(\theta - \hat{\phi}(t_k))},$$ where the sum is over a time window \(W\) around \(t\), \(\kappa\) is the concentration parameter (e.g., \(\kappa = 10\)), and \(I_0(\kappa)\) is the modified Bessel function of the first kind. Phase distribution unitarity is ensured via normalization $\int_{-\pi}^\pi \Phi_{\text{est}}(\theta; t) d\theta = 1.$ For regularization, we can project the signal into a Reproducing Kernel Hilbert Space (RKHS) with a Gaussian kernel $K(t, t') = \exp\left(-\frac{(t - t')^2}{2\sigma^2}\right), \quad \sigma = 1.$ In the RKHS, the operators become - \(K_1[s](t) = t \langle s, K(\cdot, t) \rangle_{\mathcal{H}} = t \int_0^T s(t') K(t', t) dt'\), - \(K_2[s](t) = -i \frac{d}{dt} \langle s, K(\cdot, t) \rangle_{\mathcal{H}} = -i \int_0^T s(t') \frac{\partial}{\partial t} K(t', t) dt'\), where \(\frac{\partial}{\partial t} K(t', t) = -\frac{t' - t}{\sigma^2} K(t', t)\), and - $K_3[s](t) = \int_0^T \psi(t - t'; a^*) \langle s, K(\cdot, t') \rangle_{\mathcal{H}} dt'$. In practice, it may be simpler to employ a direct signal processing approach, however, the RKHS framework ensures smoothness and regularization, which can be implemented by smoothing the signal before applying the operators. The following Approach 2 simulation is based on - Time tomain: \(t \in [0, 30]\), with \(T = 10,000\) points (\(f_s \approx 333.33 \, \text{Hz}\)), - Amplitude: \(A(t) = 1 + 0.5 \sin(2\pi \cdot 0.1 t)\), - Phase: $\phi(t) = \begin{cases} 2\pi \cdot 0.3 t & t < 5 \\ 2\pi \cdot 0.5 t + \pi/2 & t \geq 5 \end{cases}$, with frequencies \(0.3 \, \text{Hz}\) (pre-jump) and \(0.5 \, \text{Hz}\) (post-jump), and a \(\pi/2\) jump at \(t = 5\), and - Signal: \(s(t) = A(t) e^{i\phi(t)}\), with noise: \(x(t) = \text{Re}(s(t)) + \mathcal{N}(0, 0.05)\). The algorithm has the folloing steps 1. Generate the signal and compute the analytic signal using the Hilbert transform. 2. Apply operators \(K_1, K_2, K_3\) to estimate \(\phi_1(t), \phi_2(t), \phi_3(t)\). 3. Ensemble to obtain \(\hat{\phi}(t)\). 4. Estimate \(\Phi(\theta; t)\) using von Mises KDE over time windows. 5. Visualize \(\phi(t)\) vs. \(\hat{\phi}(t)\) over time and \(\Phi(\theta; t)\) vs. \(\Phi_{\text{est}}(\theta; t)\) as 3D surfaces. 6. Validate the uncertainty relation \(\Delta K_1 \Delta K_2 \geq \frac{1}{2}\). ```{r approach2_Revised, echo=TRUE, message=FALSE, warning=FALSE} # Load required libraries library(circular) # For von Mises kernel library(plotly) # For 3D surface plotting library(htmlwidgets) # For saving as SVG # Set seed for reproducibility set.seed(123) # Manual Hilbert transform function with zero-padding hilbert <- function(x, pad_factor = 2) { n <- length(x) # Zero-pad the signal to reduce edge effects pad_length <- n * (pad_factor - 1) x_padded <- c(x, rep(0, pad_length)) n_padded <- length(x_padded) # Compute the FFT of the padded signal X <- fft(x_padded) # Frequency indices for the padded signal if (n_padded %% 2 == 0) { pos_freq <- 2:(n_padded/2) neg_freq <- (n_padded/2+2):n_padded nyquist <- n_padded/2 + 1 } else { pos_freq <- 2:((n_padded+1)/2) neg_freq <- ((n_padded+1)/2+1):n_padded nyquist <- NULL } # Create the Hilbert transform multiplier h <- rep(0, n_padded) h[1] <- 1 if (length(pos_freq) > 0) h[pos_freq] <- 2 if (!is.null(nyquist)) h[nyquist] <- 1 if (length(neg_freq) > 0) h[neg_freq] <- 0 # Apply the multiplier and compute the inverse FFT X_analytic <- X * h s_analytic_padded <- fft(X_analytic, inverse = TRUE) / n_padded # Trim the padded region to return to original length s_analytic <- s_analytic_padded[1:n] return(s_analytic) } # Generate the signal T <- 10000 t <- seq(0, 30, length.out = T) dt <- t[2] - t[1] A <- 1 + 0.5 * sin(2 * pi * 0.1 * t) phi_true <- sapply(t, function(ti) { if (ti < 5) 2 * pi * 0.3 * ti else 2 * pi * 0.5 * ti + pi/2 }) s <- A * exp(1i * phi_true) x <- Re(s) + rnorm(T, 0, 0.05) # Add noise # Compute the analytic signal with zero-padding s_analytic <- hilbert(x, pad_factor = 2) # Operator K1: Multiplication by t K1_s <- t * s_analytic phi_1 <- Arg(s_analytic) %% (2 * pi) - pi # Operator K2: Derivative -i d/dt using centered difference K2_s <- rep(0+0i, T) K2_s[1] <- -1i * (s_analytic[2] - s_analytic[1]) / dt for (i in 2:(T-1)) { K2_s[i] <- -1i * (s_analytic[i+1] - s_analytic[i-1]) / (2 * dt) } K2_s[T] <- -1i * (s_analytic[T] - s_analytic[T-1]) / dt phi_2 <- Arg(K2_s) %% (2 * pi) - pi # Operator K3: Simplified wavelet convolution with fixed scale morlet_wavelet <- function(t, a = 2, omega0 = 5) { (1/sqrt(a)) * exp(-(t/a)^2/2) * cos(omega0 * t/a) * (abs(t) <= a) } K3_s <- rep(0+0i, T) for (i in 1:T) { tau <- t - t[i] psi <- morlet_wavelet(tau) K3_s[i] <- sum(psi * s_analytic) * dt } phi_3 <- Arg(K3_s) %% (2 * pi) - pi # Compute circular dispersion for weights circular_dispersion <- function(phi, window = 100) { n <- length(phi) disp <- numeric(n) for (i in 1:n) { idx <- max(1, i - window/2):min(n, i + window/2) mean_vec <- mean(exp(1i * phi[idx])) disp[i] <- 1 - abs(mean_vec) } disp } disp_1 <- circular_dispersion(phi_1) disp_2 <- circular_dispersion(phi_2) disp_3 <- circular_dispersion(phi_3) # Ensemble phase estimates epsilon <- 1e-6 w_1 <- 1 / (disp_1 + epsilon) w_2 <- 1 / (disp_2 + epsilon) w_3 <- 1 / (disp_3 + epsilon) phi_hat <- Arg(w_1 * exp(1i * phi_1) + w_2 * exp(1i * phi_2) + w_3 * exp(1i * phi_3)) phi_hat <- phi_hat %% (2 * pi) - pi # Empirical shift correction using cross-correlation # Compute cross-correlation between phi_hat and phi_true phi_true_wrapped <- phi_true %% (2 * pi) - pi cross_corr <- ccf(phi_hat, phi_true_wrapped, lag.max = 1000, plot = FALSE) lag <- cross_corr$lag[which.max(abs(cross_corr$acf))] shift_seconds <- lag * dt cat(sprintf("Detected time shift: %.3f seconds\n", shift_seconds)) # Apply the shift to phi_hat shift_steps <- round(abs(lag)) if (lag > 0) { # phi_hat is delayed, shift it back phi_hat_shifted <- rep(0, T) phi_hat_shifted[1:(T-shift_steps)] <- phi_hat[(shift_steps+1):T] phi_hat_shifted[(T-shift_steps+1):T] <- phi_hat[T-shift_steps] } else { # phi_hat is advanced, shift it forward phi_hat_shifted <- rep(0, T) phi_hat_shifted[(shift_steps+1):T] <- phi_hat[1:(T-shift_steps)] phi_hat_shifted[1:shift_steps] <- phi_hat[1] } phi_hat <- phi_hat_shifted # Estimate Phi(theta; t) using von Mises KDE theta_grid <- seq(-pi, pi, length.out = 200) t_values <- seq(0, 30, by = 3) window <- 500 # Time window for KDE kappa <- 10 # Concentration parameter phi_true_dist <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) phi_est_dist <- matrix(NA, nrow = length(t_values), ncol = length(theta_grid)) for (i in seq_along(t_values)) { t_center <- t_values[i] idx <- which(abs(t - t_center) <= window * dt / 2) # True distribution phi_window <- (phi_true[idx] %% (2 * pi)) - pi kde_true <- sapply(theta_grid, function(theta) { mean(dvonmises(theta - phi_window, mu = 0, kappa = kappa)) }) phi_true_dist[i, ] <- kde_true / sum(kde_true * (theta_grid[2] - theta_grid[1])) # Estimated distribution phi_hat_window <- phi_hat[idx] kde_est <- sapply(theta_grid, function(theta) { mean(dvonmises(theta - phi_hat_window, mu = 0, kappa = kappa)) }) phi_est_dist[i, ] <- kde_est / sum(kde_est * (theta_grid[2] - theta_grid[1])) } # Compute uncertainty relation Delta K1 * Delta K2 delta_K1 <- sd(Re(K1_s)) delta_K2 <- sd(Re(K2_s)) uncertainty_product <- delta_K1 * delta_K2 cat(sprintf("Uncertainty Product Delta K1 * Delta K2: %.6f (should be >= 0.5)\n", uncertainty_product)) # Plot 1: Phase trajectory fig_phase <- plot_ly() %>% add_trace(x = t, y = phi_true %% (2 * pi) - pi, type = "scatter", mode = "lines", name = "True Phase", line = list(color = "black")) %>% add_trace(x = t, y = phi_hat, type = "scatter", mode = "lines", name = "Estimated Phase", line = list(color = "red")) %>% layout( title = "Approach 2: True vs Estimated Phase (Shift Corrected)", xaxis = list(title = "Time (s)"), yaxis = list(title = "Phase (radians)") ) # Plot 2: Distribution Phi(theta; t) fig_true_dist <- plot_ly( x = theta_grid, y = t_values, z = phi_true_dist, type = "surface", colorscale = "Viridis", name = "True Phi", showscale = TRUE ) fig_est_dist <- plot_ly( x = theta_grid, y = t_values, z = phi_est_dist +1, # offset the Phi_est to see better the differences type = "surface", colorscale = "Magma", name = "Estimated Phi", showscale = TRUE ) fig_dist <- subplot(fig_true_dist, fig_est_dist, nrows = 1, shareX = TRUE, shareY = TRUE) %>% layout( title = "Approach 2: True vs Estimated Phi(theta; t) (Shift Corrected)", scene = list( xaxis = list(title = "θ (radians)"), yaxis = list(title = "t (seconds)"), zaxis = list(title = "Phi(θ; t)") ), scene2 = list( xaxis = list(title = "θ (radians)"), yaxis = list(title = "t (seconds)"), zaxis = list(title = "Phi_est(θ; t)") ) ) # Save plots # htmlwidgets::saveWidget(fig_phase, "approach2_phase_shift_corrected.html", selfcontained = TRUE) # htmlwidgets::saveWidget(fig_dist, "approach2_distribution_shift_corrected.html", selfcontained = TRUE) # Display plots fig_phase fig_dist # Compute MSE for the distribution mse <- mean((phi_true_dist - phi_est_dist)^2) cat(sprintf("Mean Squared Error for Phi(theta; t): %.6f\n", mse)) ``` This approach 2 estimates the kime-phase \(\phi(t)\) using non-commuting operators in time, frequency, and scale domains. It reconstructs \(\Phi(\theta; t)\) using von Mises KDE, aligning with KPT’s goal and the quantitative metrics and graphical results of the numerical simulation show the reliability in practice. ## Approach 4: A Unified Framework for Kime-Phase Tomography (KPT) Capitalizing on the strengths and advantages of the prior *Approaches 1A, 2A*, and *3A*, *Approach 4* proposes a *Unified Framework for Kime-Phase Tomography (KPT)*. ### Overview The Kime-Phase Tomography (KPT) framework relies on the Hilbert spaces for time and phase (as $\mathcal{H}_t$ and $\mathcal{H}_\theta$) and the combined kime space $\mathcal{K} = \mathcal{H}_t \otimes \mathcal{H}_\theta$ to define the central *kime-operators* (time-domain, frequency-domain, scale-domain, phase-domain, and RKHS projection). Each of the primary operators $\{K_1,K_2,K_3\}$ is specified in functional-analytic terms, with their action on signals. The wavelet operator $K_3$ is defined by a proper integral transform and the frequency operator $K_2 = -\,i \,\tfrac{d}{dt}$ ties to the commutation relation with $K_1$. The kime-commutator calculations (e.g.\,[$K_1,K_2$] $=\, i\,\mathcal{I}$) parallel the standard QM derivations for time-frequency operators, and the extension to the wavelet (scale) operator preserves the same notion of “complementary” phase information across multiple bases. Using the completeness of Fourier series on $[- \pi,\pi]$, we argue that *sufficiently many trigonometric moments uniquely determine the phase PDF*. This links multi-basis measurements and *identifiability* of $\Phi(\theta; t)$. The RKHS projections represent $\phi(t)$ via $\arg(\mathcal{P}_K[s](t))$, which is consistent with standard reproducing-kernel arguments where the kernel localizes the signal and captures amplitude-phase structure. Our numerical simulations involve the special case of *von Mises-distributed phases* and simulated fMRI signals. A synthetic neural phase signal $\theta(t)$ is sampled from a *von Mises* distribution $\mathrm{vM}\left (\overbrace{\mu(t)}^{location}, \overbrace{\kappa(t)}^{concentration}\right)$. In the simulation, we convolve the neural-phase-driven “complex time” (kime) signal with a canonical or double-$\gamma$ HRF, to simulate realistic BOLD response, and adds noise. The *primary operators* include the *time-domain operator*, $K_1[s](t) = t \cdot s(t)$, the *frequency-domain operator*, $K_2[s](t) = -\,i\,\frac{d}{dt} s(t)$ approximated by finite-difference derivatives, and the *scale-domain operator*, $K_3[s](t)$, via the continuous wavelet transform, utilizing a Morlet wavelet to locate the dominant scale. The *RKHS projection* is done before these operators for smoothing, as shown in Theorem 3. For each operator $K_j$, we use the repeated measurements (across multiple simulated runs/subjects) to estimate the phase using hte complex-argument function, $\phi_{j,n}(t) = \arg\bigl(K_j[s_n](t)\bigr)$. Then, we obtain an ensemble phase-estimate as an uncertainty-weighted combination of these phase estimates across the non-commuting bases $\{K_1, K_2, K_3\}$ to estimate the overall circular moment $\bar{m}(t)$. Then, a *von Mises* parameter solver recovers $\hat{\mu}(t)$ and $\hat{\kappa}(t)$ by matching the ratio $I_{1}(\kappa)/I_{0}(\kappa)$ to $|\bar{m}(t)|$. This follows the moment-based derivation in Theorem 4. For verification, we compute the circular MAE for phase and RMSE for the concentration by comparing the recovered parameters $\hat{\mu}(t), \hat{\kappa}(t)$ against the known ground truth. Empirical results should confirm that the multi-basis combination recovers $\mu(t)$ and $\kappa(t)$ to within reasonable accuracy. This von Mises phase prior experiment *demonstrates* the more general principle that non-commuting operators capture complementary phase information. Some *potential limitations* of the empirical validation include the reliance on von Mises parametric form, as the algorithmic performance may degrade for multi-modal or heavy-tailed distributions. Also the data requirements of $N > 200$ trials for stable recovery may be impractical for real fMRI where a typical repeats may be in the range of $N\in[10,30]$. Finally, the choice of the RKHS kernel (Gaussian) and bandwidth may not be optimal and could affect the resulting kime-phase recovery. ### Mathematical KPT Framework **Definition 1** (*Kime-Domain Signal Space*). Let $\mathcal{H}_t = L^2(\mathbb{R})$ be the Hilbert space of square-integrable complex-valued functions on the time domain, with inner product: $$\langle f, g \rangle_{\mathcal{H}_t} = \int_{\mathbb{R}} f(t) \overline{g(t)} \, dt$$ **Definition 2** (*Phase-Domain Space*). Let $\mathcal{H}_\theta = L^2([-\pi, \pi])$ be the Hilbert space of square-integrable functions on the phase domain, with inner product $$\langle \psi, \phi \rangle_{\mathcal{H}_\theta} = \int_{-\pi}^{\pi} \psi(\theta) \overline{\phi(\theta)} \, d\theta$$ equipped with periodic boundary conditions $\psi(-\pi) = \psi(\pi)$. **Definition 3** (*Kime Space*). The kime space $\mathcal{K}$ is defined as the tensor product $\mathcal{H}_t \otimes \mathcal{H}_\theta$, representing signals in both time and phase domains. **Definition 4** (*Reproducing Kernel Hilbert Space, RKHS*). The RKHS $\mathcal{R}_K$ is a subspace of $\mathcal{H}_t$ with reproducing kernel $K: \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{C}$ satisfying - For any $t \in \mathbb{R}$, $K(\cdot, t) \in \mathcal{R}_K$, and - For any $f \in \mathcal{R}_K$ and $t \in \mathbb{R}$, $f(t) = \langle f, K(\cdot, t) \rangle_{\mathcal{R}_K}$ **Definition 5** (*Kime-Phase Distribution*). A kime-phase distribution $\Phi(\theta; t)$ is a time-dependent probability density function on $[-\pi, \pi]$ satisfying $$\Phi(\theta; t) \geq 0, \quad \int_{-\pi}^{\pi} \Phi(\theta; t) \, d\theta = 1 \quad \forall t \in \mathbb{R}$$ **Definition 6** (*Complex Kime*). For each time $t$, the complex kime is defined as $\kappa(t) = t e^{i\theta(t)}$, where $\theta(t) \sim \Phi(\cdot; t)$. **Definition 7** (*Primary Kime Operators*). The primary kime operators include - *Time-domain operator*: $K_1: \mathcal{H}_t \rightarrow \mathcal{H}_t$ defined by $K_1[s](t) = t \cdot s(t)$. - *Frequency-Domain Operator*: $K_2: \mathcal{H}_t \rightarrow \mathcal{H}_t$ defined by $K_2[s](t) = -i \frac{d}{dt}s(t)$, and - *Scale-Domain Operator*: For a mother wavelet $\psi \in \mathcal{H}_t$, let $W_\psi[s](a,b) = \langle s, \psi_{a,b} \rangle_{\mathcal{H}_t}$ be the continuous wavelet transform with $\psi_{a,b}(t) = \frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right)$. Then, the scale-domain operator $K_3: \mathcal{H}_t \rightarrow \mathcal{H}_t$ is $$K_3[s](t) = \int_{\mathbb{R}} \int_{\mathbb{R}_+} W_\psi[s](a,b) \, \frac{1}{\sqrt{a}} \psi\left(\frac{t-b}{a}\right) \, \frac{da \, db}{a^2}.$$ - *Phase-Domain Operators*: In $\mathcal{H}_\theta$, we can define a pair of phase-domain operators: - *Position operator*: $\Theta[\phi](\theta) = \theta \cdot \phi(\theta)$, and - *Momentum operator*: $P[\phi](\theta) = -i\frac{d}{d\theta}\phi(\theta)$. - *RKHS Projection Operator*: Given a kernel $K$, the RKHS operator $\mathcal{P}_K: \mathcal{H}_t \rightarrow \mathcal{R}_K$ is defined by $$\mathcal{P}_K[s](t) = \int_{\mathbb{R}} s(\tau) K(t, \tau) \, d\tau .$$ **Definition 8** (Observable Signal). An observable kime-signal $s(t)$ with amplitude $A(t)$ and phase $\phi(t)$ is defined as $s(t) = A(t)e^{i\phi(t)}$, where $\phi(t)$ is sampled from distribution $\Phi(\cdot; t)$. **Definition 9** (fMRI BOLD Signal Model). In fMRI, the observed BOLD signal $x(t)$ can be modeled as $x(t) = \int_{\mathbb{R}} h(t-\tau) s(\tau) \, d\tau + \epsilon(t)$, where $h(t)$ is the *hemodynamic response function* and $\epsilon(t)$ is *noise.* Having these basic definitions, we will explore some of the mathematical results that underpin the kime-operator framework for kime-phase recovery using repeated measurement observations of a controlled experiment, such as repeated fMRI runs (within and across participants) in an event-related block design. **Theorem 1** (*Time-Frequency Commutation*). The operators $K_1$ (time-domain operator) and $K_2$ (frequency-domain operator) are *incompatible*, i.e., they have a non-trivial commutator, $[K_1, K_2] = K_1K_2 - K_2K_1 = i\mathcal{I}$, where $\mathcal{I}$ is the identity operator on $\mathcal{H}_t$. This indicates that the phase-reconstructions corresponding to this pair of kime-operators differentially probe the kime-phase and jointly, they recover complementary phase information. *Proof:* For any $s \in \mathcal{H}_t$ $$\begin{align} K_1K_2[s](t) &= t \cdot \left (-i\frac{d}{dt}s(t)\right ) = -it\frac{ds}{dt}\\ K_2K_1[s](t) &= -i\frac{d}{dt}(ts(t)) = -is(t) - it\frac{ds}{dt} \end{align} .$$ Therefore, $$\begin{align} [K_1, K_2][s](t) &= K_1K_2[s](t) - K_2K_1[s](t)\\ &= -it\frac{ds}{dt} - \left(-is(t) - it\frac{ds}{dt}\right)= is(t)= i\mathcal{I}[s](t) . \end{align}$$ Thus, $[K_1, K_2] = i\mathcal{I}$. $\square$ **Theorem 2** (*Uncertainty Relation*). Given a signal $s \in \mathcal{H}_t$ $$\Delta K_1 \cdot \Delta K_2 \geq \frac{1}{2}|\langle [K_1, K_2] \rangle| = \frac{1}{2},$$ where $\Delta K_j = \sqrt{\langle K_j^2 \rangle - \langle K_j \rangle^2}$ for $j=1,2$ and expectations are with respect to $s$. *Proof:* Assuming $\|s\|=1$ without loss of generality. Using the Cauchy-Schwarz inequality and the commutation relation from *Theorem 1*, $$\Delta K_1 \cdot \Delta K_2 \geq \frac{1}{2}|\langle [K_1, K_2] \rangle| = \frac{1}{2}|\langle i\mathcal{I} \rangle|= \frac{1}{2}|i\langle s, s \rangle|= \frac{1}{2}|i| \cdot \|s\|^2 = \frac{1}{2} \cdot 1 \cdot 1= \frac{1}{2}.\ \square$$ **Theorem 3** (*RKHS Representation*). Given a signal $s \in \mathcal{H}_t$ and a reproducing kernel $K$, the phase function $\phi(t)$ can be represented as the [complex argument](https://en.wikipedia.org/wiki/Argument_(complex_analysis)), $\arg()$, of the RKHS projection operator, $\mathcal{P}_K : \mathcal{H}_t \to \mathbb{C} \ni \mathcal{P}_K[s](t)$, i.e., $\phi(t) = \arg(\mathcal{P}_K[s](t))$. *Proof:* Let $s(t) = A(t)e^{i\phi(t)}$. Then, $$\mathcal{P}_K[s](t) = \int_{\mathbb{R}} s(\tau) K(t, \tau) \, d\tau = \int_{\mathbb{R}} A(\tau)e^{i\phi(\tau)} K(t, \tau) \, d\tau .$$ By the reproducing property and assuming a sufficiently localized kernel $$\mathcal{P}_K[s](t) \approx A(t)e^{i\phi(t)} \int_{\mathbb{R}} K(t, \tau) \, d\tau = C \cdot A(t)e^{i\phi(t)},$$ where $C$ is a non-zero constant. Therefore, $$\arg(\mathcal{P}_K[s](t)) = \arg(C \cdot A(t)e^{i\phi(t)}) = \arg(e^{i\phi(t)}) = \phi(t).\ \square$$ **Theorem 4** (*Phase Recovery from Multiple Bases*). Given repeated observations in multiple non-commuting bases defined by operators $K_1, K_2, K_3$, the kime-phase distribution $\Phi(\theta; t)$ can be *uniquely determined* if the observations are sufficient. *Proof:* Let $\phi_j(t) = \arg(K_j[s](t))$ for $j=1,2,3$ be the phase observations in each basis and consider *von Mises distribution* as an example, $\Phi(\theta; t) = \text{vM}(\mu(t), \kappa(t))$. Then, the first circular moment determines $\mu(t)$ and $\kappa(t)$ $$\mathbb{E}[e^{i\theta}] = \int_{-\pi}^{\pi} e^{i\theta} \Phi(\theta; t) \, d\theta = e^{i\mu(t)} \frac{I_1(\kappa(t))}{I_0(\kappa(t))}.$$ Since we have non-commuting observations $$\mu(t) = \arg\left(\sum_{j=1}^3 w_j(t) e^{i\phi_j(t)}\right),$$ where $w_j(t)$ are weights inversely proportional to variance. The concentration parameter $\kappa(t)$ is determined by solving the system $$\frac{I_1(\kappa(t))}{I_0(\kappa(t))} = \left|\sum_{j=1}^3 w_j(t) e^{i\phi_j(t)}\right|,$$ which uniquely determines $\Phi(\theta; t)$ when the observations provide *sufficient information* across bases. $\square$ The above *Theorem 4* (Phase Recovery from Multiple Bases) makes a strong assumption about von Mises phase distribution, which is restrictive, but the result can be generalized. **Theorem 5** (*Generalized Phase Recovery from Multiple Bases*). Given a kime-phase distribution $\Phi(\theta; t)$ and assuming sufficient observations in multiple non-commuting bases defined by operators $K_1, K_2, K_3$, the phase distribution can be uniquely determined under certain uniqueness conditions 1. *Trigonometric Moment Identifiability*: A circular distribution is uniquely determined by its complete set of trigonometric moments $\{\alpha_k, \beta_k\}_{k=1}^{\infty}$ where $\alpha_k = \mathbb{E}[\cos(k\theta)]$ and $\beta_k = \mathbb{E}[\sin(k\theta)]$, 2. *Information Complementarity*: The operators $K_1, K_2, K_3$ must provide complementary information about different moments of the phase distribution, and 3. *Sufficiency Condition*: The observations must constrain enough trigonometric moments to uniquely specify $\Phi(\theta; t)$ within the class of distributions being considered. **Proof**: The proof proceeds in the following steps. Step 1: reframe the problem in terms of the moment problem. Let $\mathcal{P}([-\pi, \pi])$ be the space of probability distributions on $[-\pi, \pi]$. For any $\Phi \in \mathcal{P}([-\pi, \pi])$, define its trigonometric moments: $$m_k(\Phi) = \int_{-\pi}^{\pi} e^{ik\theta} \Phi(\theta) d\theta, \quad k \in \mathbb{Z}.$$ Note that $m_0(\Phi) = 1$ (normalization), $m_{-k}(\Phi) = \overline{m_k(\Phi)}$ (conjugate symmetry), and $|m_k(\Phi)| \leq 1$ (bounded). Step 2: To establish the connection between moments and distribution uniqueness requires the following lemma, which states that distributions are completely described by their moments, see the note on MGF's at the end of this section. **Lemma 1:** Two distributions $\Phi_1, \Phi_2 \in \mathcal{P}([-\pi, \pi])$ are identical if and only if $m_k(\Phi_1) = m_k(\Phi_2)$ for all $k \in \mathbb{Z}$. *Proof of Lemma 1:* This follows directly from the uniqueness of Fourier series representation. If $\Phi_1 \neq \Phi_2$, then their difference $\Phi_1 - \Phi_2$ has a non-zero Fourier coefficient, meaning some moment $m_k(\Phi_1) \neq m_k(\Phi_2)$. Conversely, if all moments match, the distributions must be identical almost everywhere. Step 3: Next we connect operator measurements to moment estimation. For each operator $K_j$, we observe phases $\{\phi_{j,n}(t)\}_{n=1}^N$ across $N$ repeated measurements at time $t$. These provide empirical estimates of certain functions of the moments $$\hat{m}_{j,k}(t) = \frac{1}{N}\sum_{n=1}^N e^{ik\phi_{j,n}(t)}$$ However, due to the non-commutativity of the operators, these empirical moments capture different aspects of the underlying distribution $\Phi(\theta; t)$, which requires another lemma. **Lemma 2:** Under the action of operator $K_j$, the expected value of $\hat{m}_{j,k}(t)$ is a function $F_j$ of the true moments $\{m_l(\Phi)\}_{l \in \mathbb{Z}}$ $$\mathbb{E}[\hat{m}_{j,k}(t)] = F_j(k, \{m_l(\Phi)\}_{l \in \mathbb{Z}}),$$ where $F_j$ is a continuous function that depends on the specific operator $K_j$. *Proof of Lemma 2:* Each operator $K_j$ transforms the input signal according to its action, resulting in a phase that depends on the underlying distribution in a specific way. The exact form of $F_j$ depends on the operator definition, but it is a continuous function of the moments due to the continuous nature of the operators. Step 4: To show that non-commuting operators provide complementary information we need a third lemma. **Lemma 3:** If operators $K_i$ and $K_j$ do not commute (i.e., $[K_i, K_j] \neq 0$), then there exist moments $m_l(\Phi)$ such that $$\frac{\partial F_i(k, \{m_l(\Phi)\})}{\partial m_l(\Phi)} \neq \frac{\partial F_j(k, \{m_l(\Phi)\})}{\partial m_l(\Phi)}.$$ In other words, the operators have different sensitivities to certain moments of the distribution. *Proof of Lemma 3:* This follows from the definition of non-commutativity. If two operators commute, they share eigenfunctions and can be simultaneously diagonalized. Non-commuting operators have different eigenbases, which means they respond differently to certain input patterns—specifically, they have different sensitivities to certain frequency components or moments of the input distribution. Step 5: To establish the conditions for unique determination, let's define the *total information* from all operators as the set $$\mathcal{I}(\Phi) = \{F_j(k, \{m_l(\Phi)\}) : j \in \{1,2,3\}, k \in \mathbb{Z}\}.$$ The theorem's main result is that the distribution $\Phi(\theta; t)$ is uniquely determined by the observations *if and only if* the mapping $\Phi \mapsto \mathcal{I}(\Phi)$ is injective on the space of distributions being considered (a *necessary and sufficient* condition). Let's show both directions. 1. ($\Longrightarrow$) If $\Phi \mapsto \mathcal{I}(\Phi)$ is injective, then different distributions produce different observable patterns, allowing unique identification. 2. ($\Longleftarrow$) Conversely, suppose $\Phi \mapsto \mathcal{I}(\Phi)$ is not injective. Then there exist distinct distributions $\Phi_1 \neq \Phi_2$ such that $\mathcal{I}(\Phi_1) = \mathcal{I}(\Phi_2)$. This means all observations would be identical, making it impossible to distinguish between $\Phi_1$ and $\Phi_2$. Step 6: Let's specialize to practical cases. For parametric families with finite-dimensional parameterization (e.g., von Mises or mixtures thereof), the injectivity condition reduces to a more tractable form shown in the following two corrolaries. **Corollary 1:** For a parametric family $\{\Phi_\theta : \theta \in \Theta \subset \mathbb{R}^d\}$ where $\Theta$ is a $d$-dimensional parameter space, unique determination is possible if the Jacobian matrix $$J(\theta) = \left[ \frac{\partial}{\partial \theta_i} F_j(k, \{m_l(\Phi_\theta)\}) \right]_{(j,k),i}$$ has full column rank for all $\theta \in \Theta$. This is a direct result of the inverse function theorem from multivariable calculus. If the Jacobian has full column rank, the mapping is locally injective around each point in the parameter space. If this holds globally, the mapping is globally injective. **Corollary 2:** For a von Mises distribution $\text{vM}(\mu, \kappa)$, observations from operators $K_1$ and $K_2$ are sufficient for unique determination if they provide independent constraints on the first trigonometric moment $m_1 = e^{i\mu}\frac{I_1(\kappa)}{I_0(\kappa)}$. The von Mises distribution is completely determined by its first trigonometric moment $m_1$. If the operators $K_1$ and $K_2$ provide independent constraints on $m_1$ (which follows from their non-commutativity), the parameters $\mu$ and $\kappa$ can be uniquely determined. Therefore, the *generalized theorem* establishes *necessary* and *sufficient* conditions for unique phase distribution determination using multiple non-commuting operators. The key insight is that non-commuting operators provide complementary views of the underlying distribution, helping constrain it more effectively than any single operator could. The degree to which uniqueness is possible depends on both the complexity of the true distribution and the information content of the operators used. $\square$ **Notes**: Due to the non-commutativity of $K_1, K_2, K_3$, these operators capture different aspects of the phase distribution - $K_1$ (time-domain): Sensitive to localized phase concentrations - $K_2$ (frequency-domain): Sensitive to phase transitions - $K_3$ (scale-domain): Sensitive to multi-scale phase patterns. Assuming sufficient diversity of measurements, the combined information constrains the characteristic function $\varphi(k)$ enough to uniquely determine $\Phi(\theta; t)$. Let's consider some specific families of phase distributions. 1. **Parametric Families**: For distributions with finite-dimensional parameterizations (like von Mises, wrapped Cauchy, wrapped normal), we only need to constrain finitely many moments to uniquely determine the distribution. 2. **Mixtures of Parametric Distributions**: For mixtures of $M$ parametric components, uniqueness requires constraining $2M-1$ complex moments $\{\varphi(1), \varphi(2), ..., \varphi(2M-1)\}$. 3. **Non-parametric Distributions**: Any distribution can be approximated arbitrarily well by constraining a sufficiently large number of trigonometric moments, with the approximation error decreasing as more moments are constrained. The theorem is practically useful even without knowing the true distribution family, as the multi-basis approach constrains more trigonometric moments than single-basis approaches. For most neurophysiological signals, the phase distribution is well-approximated by a low-order mixture model (typically 1-3 components), which requires only a small number of moments to constrain. The non-commuting bases provide complementary information about different moments, improving identification even with finite data. This generalized theorem provides a broader foundation for KPT, allowing application to a wider range of phase distributions beyond von Mises. For practical applications, we can start with the von Mises assumption (which often works well) but can extend to mixture models or non-parametric approaches when data suggests more complex distributions. ### The Moment Problem in Circular Distributions The question of whether a distribution is uniquely determined by its moments is known as the "*moment problem*" in probability theory. For circular distributions on $[-\pi, \pi]$, we need to consider the trigonometric moments rather than ordinary moments. For a circular distribution $\Phi(\theta)$ on $[-\pi, \pi]$, the trigonometric moments are $$m_k = \mathbb{E}[e^{ik\theta}] = \int_{-\pi}^{\pi} e^{ik\theta} \Phi(\theta) d\theta = \alpha_k + i\beta_k,$$ where $\alpha_k = \mathbb{E}[\cos(k\theta)]$ and $\beta_k = \mathbb{E}[\sin(k\theta)]$. Unlike ordinary moments on the real line, which can sometimes fail to characterize a distribution uniquely (the famous "moment problem"), trigonometric moments on a bounded circular domain always uniquely determine the distribution, due to the completeness of the trigonometric functions. The set $\{e^{ik\theta}\}_{k \in \mathbb{Z}}$ forms a *complete orthogonal basis* for $L^2[-\pi, \pi]$. By the Fourier theory, any square-integrable function on $[-\pi, \pi]$ can be expressed as $f(\theta) = \sum_{k=-\infty}^{\infty} c_k e^{ik\theta}$, where $c_k$ are the Fourier coefficients. For a probability density $\Phi(\theta)$, these coefficients are precisely related to the trigonometric moments, $c_k = \frac{1}{2\pi} \overline{m_k}$. Thus, knowing all trigonometric moments $\{m_k\}_{k \in \mathbb{Z}}$ is equivalent to knowing the exact Fourier series of $\Phi(\theta)$, which uniquely determines the distribution. In the circular setting, the *characteristic function* serves the role that the MGF does for distributions on the real line. The characteristic function of a circular distribution is $$\varphi(t) = \mathbb{E}[e^{it\theta}] = \int_{-\pi}^{\pi} e^{it\theta} \Phi(\theta) d\theta .$$ This is precisely the sequence of trigonometric moments when evaluated at integer values $\varphi(k) = m_k, \quad \text{for } k \in \mathbb{Z}$. The characteristic function directly encodes all trigonometric moments and has these important properties 1. **Uniqueness**: If two distributions have the same characteristic function, they are identical. 2. **Inversion Formula**: The distribution can be recovered from its characteristic function using the inversion formula $$\Phi(\theta) = \frac{1}{2\pi} \sum_{k=-\infty}^{\infty} \overline{m_k} e^{ik\theta}.$$ 3. **Convolutions**: The characteristic function of a sum of independent random variables is the product of their characteristic functions. ### Computational Algorithm 1: Unified Kime-Phase Tomography Let's explicate the *Unified Kime-Phase Tomography* algorithm in the special case of using *von Mises phase prior distribution*. **Input:** BOLD time series $\{x_n(t)\}_{n=1}^N$ from $N$ repeated measurements. **Output:** Estimated kime-phase distribution $\hat{\Phi}(\theta; t)$. 1. **Preprocessing:** First apply bandpass filtering to each $x_n(t)$ to remove noise and then compute analytic signal via Hilbert transform, $s_n(t) = x_n(t) + i\mathcal{H}[x_n(t)]$. 2. **Multi-basis Measurement:** Obtain the three phase recovery estimates using each o fhte kime-operators (in each of the 3 bases). - *Time-domain*: $\phi_{1,n}(t) = \arg(K_1[s_n](t))$, - *Frequency-domain*: $\phi_{2,n}(t) = \arg(K_2[s_n](t))$, and - *Scale-domain*: $\phi_{3,n}(t) = \arg(K_3[s_n](t))$. 3. **Basis Uncertainty Quantification:** Compute the circular variance in each basis $V_j(t) = 1 - \left|\frac{1}{N}\sum_{n=1}^N e^{i\phi_{j,n}(t)}\right|,\ j\in\{1,2,3\}$ and estimate the corresponding weights, $w_j(t) = \frac{1/V_j(t)}{\sum_{k=1}^3 1/V_k(t)}$. 4. **Phase Distribution Estimation:** To obtain an *ensemble* phase reconstruction, we first compute the *raw moment* $m_j(t) = \frac{1}{N}\sum_{n=1}^N e^{i\phi_{j,n}(t)}$ and then estimate the *weighted moment* $\bar{m}(t) = \sum_{j=1}^3 w_j(t) m_j(t)$. With these, we can approximate the von Mises parameters, $\hat{\mu}(t) = \arg(\bar{m}(t))$ (*mean* or *location*) and $\hat{\kappa}(t)$ (*concentration*) by solving $\frac{I_1(\hat{\kappa}(t))}{I_0(\hat{\kappa}(t))} = |\bar{m}(t)|$. 5. **Return** The result of this algorithm is an estimated phase recovery distribution $\hat{\Phi}(\theta; t) = \text{vM}(\hat{\mu}(t), \hat{\kappa}(t))$. ### Algorithm 2: Robust Phase Jump Detection To handle more intricate kime-phase prior distributions, such as mixtures of various distributions, we need to extend and modify *Alforithm 1* as follows. **Input:** Estimated phase series $\hat{\phi}(t)$. **Output:** Detected phase jumps $\{(t_k, \Delta\phi_k)\}$. 1. Compute phase derivative: $\dot{\phi}(t) = \frac{d\hat{\phi}(t)}{dt}$ 2. Compute median absolute deviation: $\text{MAD} = \text{median}(|\dot{\phi}(t) - \text{median}(\dot{\phi}(t))|)$ 3. Set threshold: $\tau = \text{median}(\dot{\phi}(t)) + 5 \cdot \text{MAD}$ 4. For each time $t$, where $|\dot{\phi}(t)| > \tau$, record the jump, $(t, \hat{\phi}(t+\delta) - \hat{\phi}(t-\delta))$, where $\delta$ is a small time window 5. Merge consecutive jumps within a threshold window 6. Return list of jumps $\{(t_k, \Delta\phi_k)\}$. ### Algorithm 3: RKHS-Based Phase Recovery To capitalize on *Approach 3A* (above), we can also refine the algorithm to allow for RKHS kernel estimation. This RKHS-based algorithm is more refined implementation since it directly incorporates the kernel-based smoothing and projection emphasized in the kime operator-theoretic framework. The earlier Algorithm 1 is valid, but it represents essentially a simpler “baseline” version (no RKHS). The recovered phase $\phi_{1}^{K}(t)=\arg\bigl(t \cdot s_n^K(t)\bigr)$, is consistent with keeping phases in $\bigl[-\pi,\pi\bigr)$, since the complex-argument function $\arg(\cdot)$ returns the principal argument in that exact interval. Note that *Algorithm 2* may be more accurate than *Algorithm 1*, since it adds the extra step $s_n^K(t) \;=\;\mathcal{P}_K\bigl[s_n\bigr](t)$ to project each subject’s/time series signal into the chosen reproducing-kernel Hilbert space (RKHS). This is a smoothing or *regularization* necessary for real fMRI data where the noise level is significant, e.g., $SNR < 0.04$. Theorem 3 shows that $\phi(t)\approx\arg\bigl(\mathcal{P}_K[s](t)\bigr)$ where the kernel projection represents the mathematical device to appropriately localize the signal and reduce the noise. While, *Algorithm 1* uses the same multi-basis logic—time, frequency, scale operators, it omits the kernel-based projection (regularization) step. The simpler *Algorithm 1* may still work well in synthetic or high-SNR settings, it may be suboptimal for noisy signals needing more advanced smoothing. This is reflected in teh following *Algorithm 3* included below. **Input:** BOLD time series $\{x_n(t)\}_{n=1}^N$ and kernel $K$. **Output:** Estimated phase $\hat{\phi}(t)$. 1. Project each signal to RKHS: $s_n^K(t) = \mathcal{P}_K[s_n](t)$ 2. Compute time-domain phase: $\phi_1^K(t) = \arg(t \cdot s_n^K(t))$ 3. Compute frequency-domain phase using STFT, $\phi_2^K(t) = \arg(\mathcal{F}[s_n^K \cdot w](t, \omega_{\text{max}}))$, where $\omega_{\text{max}}$ is the frequency with maximum power 4. Compute scale-domain phase using CWT: $\phi_3^K(t) = \arg(W_\psi[s_n^K](a_{\text{max}}, t))$, where $a_{\text{max}}$ is the scale with maximum power 5. Compute weighted average: $\hat{\phi}(t) = \arg\left(\sum_{j=1}^3 w_j(t) e^{i\phi_j^K(t)}\right)$ 6. Return $\hat{\phi}(t)$. Note that the (simulated) true phases should be wrapped to $[-\pi, \pi)$ to match the RKHS-recovered estimates. The recovered phases, e.g., `phi_K1`, `phi_K2`, `phi_K3`, use complex-argument function $\arg()$ to restrict the phase values to $[-\pi, \pi)$, e.g., `phi_K1[n,] <- Arg(K1_operator(s_rkhs, t))`. In our experiment, the true phase is generated via `mu_func = function(t) 2*pi*0.1*t`, which exceeds \([-\pi, \pi)\) as \(t\) grows. However, `rvonmises` internally wraps the mean parameter `mu` to \([-\pi, \pi)\), so the *sampled phases* are correct. The stored `true_mean` (ground truth) needs to also be wrapped to avoid a mismatch when comparing true phases to their corresponding estimates. ```{r approach4_KPT_Demo1, warning=FALSE, message=FALSE} # Complete Implementation of Unified Kime-Phase Tomography # with Application to fMRI Data Analysis library(signal) library(circular) library(plotly) library(dplyr) library(RColorBrewer) # -------------------------------------------------------------- # Core Kime-Phase Tomography Functions # -------------------------------------------------------------- #' Compute analytic signal via Hilbert transform compute_analytic_signal <- function(x) { n <- length(x) X <- fft(x) # Create Hilbert transform filter h <- numeric(n) h[1] <- 1 # DC component h[2:(n/2)] <- 2 # Positive frequencies h[(n/2 + 1):n] <- 0 # Negative frequencies # Apply filter and compute inverse FFT analytic <- fft(h * X, inverse = TRUE) / n return(analytic) } #' Create hemodynamic response function create_hrf <- function(t_hrf, type = "canonical") { if (type == "gamma") { hrf <- dgamma(t_hrf, shape = 6, scale = 1) } else if (type == "double-gamma") { hrf <- dgamma(t_hrf, shape = 6, scale = 1) - 0.1 * dgamma(t_hrf, shape = 16, scale = 1) } else { # canonical hrf <- (t_hrf/1.5)^3 * exp(-(t_hrf/1.5)) / gamma(4) } return(hrf / max(hrf)) # Normalize } # PLOT the HDR function # Create time points for the HRF t_hrf <- seq(0, 20, by = 0.1) # Calculate the three HRF types hrf_canonical <- create_hrf(t_hrf, type = "canonical") hrf_gamma <- create_hrf(t_hrf, type = "gamma") hrf_double_gamma <- create_hrf(t_hrf, type = "double-gamma") # Create the interactive plot plot_ly() %>% add_trace(x = t_hrf, y = hrf_canonical, type = "scatter", mode = "lines", name = "Canonical HRF", line = list(color = "#3366CC", width = 2), hoverinfo = "text", text = ~paste("Time:", round(t_hrf, 1), "s
", "Amplitude:", round(hrf_canonical, 3))) %>% add_trace(x = t_hrf, y = hrf_gamma, type = "scatter", mode = "lines", name = "Gamma HRF", line = list(color = "#DC3912", width = 2, dash = "dash"), hoverinfo = "text", text = ~paste("Time:", round(t_hrf, 1), "s
", "Amplitude:", round(hrf_gamma, 3))) %>% add_trace(x = t_hrf, y = hrf_double_gamma, type = "scatter", mode = "lines", name = "Double-Gamma HRF", line = list(color = "#109618", width = 2, dash = "dot"), hoverinfo = "text", text = ~paste("Time:", round(t_hrf, 1), "s
", "Amplitude:", round(hrf_double_gamma, 3))) %>% layout(title = list(text = "Hemodynamic Response Function (HRF) Models",font = list(size = 18)), xaxis = list(title = "Time (seconds)", titlefont = list(size = 14), zeroline = TRUE, zerolinecolor = "#CCCCCC", gridcolor = "rgba(0,0,0,0.1)"), yaxis = list(title = "Normalized Amplitude", titlefont = list(size = 14), zeroline = TRUE, zerolinecolor = "#CCCCCC", gridcolor = "rgba(0,0,0,0.1)", range = c(-0.2, 1.05) # Allow space for negative values in double-gamma ), legend = list(x = 0.3, y = 0.99, xanchor = "left", yanchor = "top", bgcolor = "rgba(255,255,255,0.7)", bordercolor = "#CCCCCC", borderwidth = 1), shapes = list(# Add a horizontal line at y=0 list(type = "line", x0 = 0, x1 = 20, y0 = 0, y1 = 0, line=list(color="#CCCCCC", width = 1))), annotations = list(# Add annotations explaining each model list(x = 5, y = 0.9, text = "Canonical: (t/1.5)³ exp(-t/1.5) / Γ(4)", showarrow = FALSE, font = list(color = "#3366CC")), list(x = 5, y = 0.8, text = "Gamma: Γ(shape=6, scale=1)", showarrow = FALSE, font = list(color = "#DC3912")), list(x = 7, y = 0.7, text = "Double-Gamma: Γ(shape=6, scale=1) - 0.1×Γ(shape=16, scale=1)", showarrow = FALSE, font = list(color = "#109618"))), hovermode = "closest", hoverlabel = list(bgcolor = "white",font = list(size = 12)), paper_bgcolor = "white", plot_bgcolor = "white", margin = list(t=80, r=50, b=60, l=60)) %>% # Add buttons to highlight individual curves config(modeBarButtonsToAdd = list( list(name = "Show All", icon = list(path = "M22 12c0 5.5-4.5 10-10 10S2 17.5 2 12 6.5 2 12 2s10 4.5 10 10zm-2 0c0-4.4-3.6-8-8-8s-8 3.6-8 8 3.6 8 8 8 8-3.6 8-8zm-8 4c2.2 0 4-1.8 4-4s-1.8-4-4-4-4 1.8-4 4 1.8 4 4 4z"), click = htmlwidgets::JS("function(gd) { Plotly.restyle(gd, 'visible', true); }"))), displaylogo = FALSE ) # Define complex time-varying phase with two jumps complex_mu <- function(t) { ifelse(t < 10, 0.2*t, ifelse(t < 20, 0.2*10 + 0.4*(t-10) + pi/2, 0.2*10 + 0.4*10 + pi/2 + 0.3*(t-20) + pi/4)) } # Time-varying concentration parameter complex_kappa <- function(t) { 5 + 2*sin(2*pi*t/30) + t/10 } #' Generate synthetic fMRI data with known phase properties generate_synthetic_fmri <- function(N = 10, T = 500, tmax = 30, mu_func = complex_mu, # function(t) 2*pi*0.1*t, kappa_func = complex_kappa, # function(t) 5 + 2*sin(2*pi*t/tmax), # +t/10, noise_sd = 0.1, hrf_type = "canonical") { # Time vector t <- seq(0, tmax, length.out = T) dt <- t[2] - t[1] # HRF convolution kernel t_hrf <- seq(0, 20, length.out = 100) hrf <- create_hrf(t_hrf, type = hrf_type) # Generate signals for each repetition true_phases <- matrix(NA, nrow = N, ncol = T) bold_signals <- matrix(NA, nrow = N, ncol = T) for (n in 1:N) { # Generate true phase from von Mises distribution at each time point true_mean <- sapply(t, mu_func) true_kappa <- sapply(t, kappa_func) # Sample phases from von Mises distribution phases <- numeric(T) for (i in 1:T) { phases[i] <- rvonmises(1, mu = true_mean[i], kappa = true_kappa[i]) } # Create complex signal neural_signal <- exp(1i * phases) # Convolve with HRF to get BOLD signal bold <- Re(convolve(neural_signal, rev(hrf), type = "open"))[1:T] # Add noise bold_noisy <- bold + rnorm(T, 0, noise_sd) # Store results true_phases[n, ] <- phases bold_signals[n, ] <- bold_noisy } # Modify true_mean generation: true_mean <- sapply(t, mu_func) %% (2*pi) # Wrap to [0, 2π) true_mean <- ifelse(true_mean > pi, true_mean - 2*pi, true_mean) # Convert to [-π, π) # Return results with ground truth return(list( t = t, bold = bold_signals, true_phases = true_phases, true_mean = true_mean, true_kappa = true_kappa, hrf = hrf )) } #' Kime operator K1 (time-domain) K1_operator <- function(s, t) { return(t * s) } #' Kime operator K2 (frequency-domain) K2_operator <- function(s, t) { dt <- t[2] - t[1] # Approximate derivative using central differences ds <- c(0, diff(s)) / dt # Return -i * d/dt return(-1i * ds) } #' Compute continuous wavelet transform compute_cwt <- function(s, t, scales) { # Define Morlet wavelet morlet <- function(t, scale, w0 = 6) { norm <- (pi^(-0.25)) / sqrt(scale) return(norm * exp(1i * w0 * t / scale) * exp(-0.5 * (t / scale)^2)) } n <- length(t) wt <- matrix(0, length(scales), n) # Compute wavelet transform for each scale for (i in seq_along(scales)) { scale <- scales[i] # Create wavelet at this scale t_kernel <- seq(-5 * scale, 5 * scale, length.out = min(201, n)) psi <- morlet(t_kernel, scale) # Convolve signal with wavelet conv_full <- convolve(s, rev(psi), type = "open") # Extract relevant part of convolution start_idx <- floor(length(conv_full) / 2) - floor(n / 2) + 1 wt[i, ] <- conv_full[start_idx:(start_idx + n - 1)] } return(wt) } #' Kime operator K3 (scale-domain) K3_operator <- function(s, t, scales = NULL) { # Set default scales if not provided if (is.null(scales)) { dt <- t[2] - t[1] fs <- 1/dt f_min <- 0.01 f_max <- fs/10 scales <- 1 / (2 * pi * seq(f_min, f_max, length.out = 32)) } # Compute wavelet transform wt <- compute_cwt(s, t, scales) # For each time point, find scale with maximum power dominant <- numeric(length(t)) for (i in 1:ncol(wt)) { dominant[i] <- wt[which.max(abs(wt[, i])), i] } return(dominant) } #' Project signal to RKHS rkhs_projection <- function(s, t, sigma = 0.5) { dt <- t[2] - t[1] n <- length(t) projected <- numeric(n) # Define Gaussian kernel K <- function(t1, t2) exp(-((t1 - t2)^2) / (2 * sigma^2)) # Apply kernel to signal for (i in 1:n) { kernel_vals <- sapply(t, function(tj) K(t[i], tj)) projected[i] <- sum(s * kernel_vals) * dt } return(projected) } #' Estimate von Mises parameters from complex-valued data estimate_vonmises <- function(z) { # Compute resultant length R <- abs(mean(z)) # Compute mean direction mu <- Arg(mean(z)) # Estimate concentration parameter if (R >= 0.999) { kappa <- 100 # High concentration } else if (R <= 0.001) { kappa <- 0 # Uniform distribution } else { # Solve for kappa using Bessel function ratio kappa_est <- function(kappa) { besselI(kappa, 1) / besselI(kappa, 0) - R } # Solve numerically with appropriate bounds kappa <- tryCatch({ uniroot(kappa_est, interval = c(0, 100))$root }, error = function(e) { # Fallback approximation if (R < 0.53) { 2 * R + R^3 + 5 * R^5 / 6 } else { 1 / (2 * (1 - R) - (1 - R)^2 - (1 - R)^3) } }) } return(list(mu = mu, kappa = kappa)) } #' Detect phase jumps in a phase time series detect_phase_jumps <- function(phase, t, threshold_mult = 5) { # Unwrap phase for continuous derivative unwrapped <- phase for (i in 2:length(phase)) { diff_phase <- phase[i] - phase[i-1] if (abs(diff_phase) > pi) { unwrapped[i:length(phase)] <- unwrapped[i:length(phase)] - sign(diff_phase) * 2 * pi } } # Compute numerical derivative dt <- t[2] - t[1] phase_deriv <- c(0, diff(unwrapped) / dt) # Compute median and MAD for robust threshold med_deriv <- median(phase_deriv, na.rm = TRUE) mad_deriv <- median(abs(phase_deriv - med_deriv), na.rm = TRUE) # Set threshold threshold <- med_deriv + threshold_mult * mad_deriv # Find jumps jump_indices <- which(abs(phase_deriv) > threshold) # Merge consecutive jumps if (length(jump_indices) > 0) { merged_indices <- jump_indices[1] for (i in 2:length(jump_indices)) { if (jump_indices[i] - jump_indices[i-1] > 3) { merged_indices <- c(merged_indices, jump_indices[i]) } } # Create output data frame jumps <- data.frame( time = t[merged_indices], magnitude = phase_deriv[merged_indices] ) return(jumps) } else { return(data.frame(time = numeric(0), magnitude = numeric(0))) } } #' Unified Kime-Phase Tomography algorithm unified_kpt <- function(bold_signals, t, preprocess = TRUE, kernel_sigma = 0.5) { N <- nrow(bold_signals) T <- ncol(bold_signals) dt <- t[2] - t[1] # Store results phi_K1 <- matrix(NA, N, T) phi_K2 <- matrix(NA, N, T) phi_K3 <- matrix(NA, N, T) # Step 1: Preprocessing analytic_signals <- matrix(complex(real = 0, imaginary = 0), N, T) for (n in 1:N) { signal <- bold_signals[n,] if (preprocess) { # Bandpass filter (0.01-0.1 Hz) nyquist <- 1/(2*dt) lowcut <- 0.01/nyquist highcut <- 0.1/nyquist butter_filter <- signal::butter(2, c(lowcut, highcut), type = "pass") filtered <- signal::filtfilt(butter_filter, signal) # Get analytic signal via Hilbert transform analytic_signals[n,] <- compute_analytic_signal(filtered) } else { # Skip filtering analytic_signals[n,] <- compute_analytic_signal(signal) } # Step 2: Multi-basis measurement # Apply RKHS projection if requested if (kernel_sigma > 0) { s_rkhs <- rkhs_projection(analytic_signals[n,], t, kernel_sigma) # Apply operators in each domain phi_K1[n,] <- Arg(K1_operator(s_rkhs, t)) phi_K2[n,] <- Arg(K2_operator(s_rkhs, t)) phi_K3[n,] <- Arg(K3_operator(s_rkhs, t)) } else { # Apply operators directly without RKHS projection phi_K1[n,] <- Arg(K1_operator(analytic_signals[n,], t)) phi_K2[n,] <- Arg(K2_operator(analytic_signals[n,], t)) phi_K3[n,] <- Arg(K3_operator(analytic_signals[n,], t)) } } # Step 3: Basis uncertainty quantification # Compute circular variance for each basis and time point v_K1 <- v_K2 <- v_K3 <- numeric(T) for (i in 1:T) { v_K1[i] <- 1 - abs(mean(exp(1i * phi_K1[,i]))) v_K2[i] <- 1 - abs(mean(exp(1i * phi_K2[,i]))) v_K3[i] <- 1 - abs(mean(exp(1i * phi_K3[,i]))) } # Calculate weights inversely proportional to variance w1 <- w2 <- w3 <- numeric(T) for (i in 1:T) { # Avoid division by zero with a small epsilon eps <- 1e-10 total_weight <- 1/(v_K1[i] + eps) + 1/(v_K2[i] + eps) + 1/(v_K3[i] + eps) w1[i] <- 1/(v_K1[i] + eps) / total_weight w2[i] <- 1/(v_K2[i] + eps) / total_weight w3[i] <- 1/(v_K3[i] + eps) / total_weight } # Step 4: Phase distribution estimation mu_est <- kappa_est <- numeric(T) for (i in 1:T) { # Compute moments m1 <- mean(exp(1i * phi_K1[,i])) m2 <- mean(exp(1i * phi_K2[,i])) m3 <- mean(exp(1i * phi_K3[,i])) # Weighted combination m_combined <- w1[i] * m1 + w2[i] * m2 + w3[i] * m3 # Estimate von Mises parameters vm_params <- estimate_vonmises(c(m_combined)) mu_est[i] <- vm_params$mu kappa_est[i] <- vm_params$kappa } # Detect phase jumps jumps <- detect_phase_jumps(mu_est, t) return(list( t = t, mu = mu_est, kappa = kappa_est, phi_K1 = phi_K1, phi_K2 = phi_K2, phi_K3 = phi_K3, weights = list(w1 = w1, w2 = w2, w3 = w3), variances = list(v_K1 = v_K1, v_K2 = v_K2, v_K3 = v_K3), phase_jumps = jumps )) } #' Evaluate KPT performance evaluate_kpt <- function(kpt, ground_truth) { # Compute circular mean absolute error for phase cmae <- function(true, est) { # Ensure phases are in [-pi, pi] true <- ((true + pi) %% (2*pi)) - pi est <- ((est + pi) %% (2*pi)) - pi # Compute circular difference diff <- abs(((true - est + pi) %% (2*pi)) - pi) return(mean(diff)) } # Compute RMSE for concentration kappa_rmse <- sqrt(mean((ground_truth$true_kappa - kpt$kappa)^2)) # Compute circular MAE for phase phase_cmae <- cmae(ground_truth$true_mean, kpt$mu) # Compute phase jump detection accuracy true_jumps <- detect_phase_jumps(ground_truth$true_mean, kpt$t) # Match detected jumps to true jumps jump_detection <- list( true_jumps = true_jumps, detected_jumps = kpt$phase_jumps ) # Return metrics return(list( phase_cmae = phase_cmae, kappa_rmse = kappa_rmse, jump_detection = jump_detection )) } #' Visualize Kime-Phase estimation results visualize_kpt_results <- function(kpt, ground_truth = NULL) { # Create basic data frame df <- data.frame( t = kpt$t, mu_est = kpt$mu, kappa_est = kpt$kappa ) # Add ground truth if available has_truth <- !is.null(ground_truth) if (has_truth) { df$mu_true <- ground_truth$true_mean df$kappa_true <- ground_truth$true_kappa } # Phase plot p1 <- plot_ly(df, x = ~t) %>% add_trace(y = ~mu_est, type = "scatter", mode = "lines", name = "Estimated Mean Phase", line = list(color = "blue", width = 2)) if (has_truth) { p1 <- p1 %>% add_trace(y = ~mu_true, type = "scatter", mode = "lines", name = "True Mean Phase µ(t)", line = list(color = "black", dash = "dot")) } p1 <- p1 %>% layout( title = "Estimated Mean Phase µ(t)", xaxis = list(title = "Time"), yaxis = list(title = "Phase (radians)") ) # Concentration plot p2 <- plot_ly(df, x = ~t) %>% add_trace(y = ~kappa_est, type = "scatter", mode = "lines", name = "Estimated Concentration", line = list(color = "red", width = 2)) if (has_truth) { p2 <- p2 %>% add_trace(y = ~kappa_true, type = "scatter", mode = "lines", name = "True Concentration κ(t)", line = list(color = "black", dash = "dot")) } p2 <- p2 %>% layout( title = "Estimated Concentration κ(t)", xaxis = list(title = "Time"), yaxis = list(title = "Concentration Parameter", range = list(0.0, 40.0)) ) # Create weights plot p3 <- plot_ly(x = kpt$t) %>% add_trace(y = kpt$weights$w1, type = "scatter", mode = "lines", name = "Time-Domain (K₁)", line = list(color = "blue")) %>% add_trace(y = kpt$weights$w2, type = "scatter", mode = "lines", name = "Frequency-Domain (K₂)", line = list(color = "red")) %>% add_trace(y = kpt$weights$w3, type = "scatter", mode = "lines", name = "Scale-Domain (K₃)", line = list(color = "green")) %>% layout( title = "Basis Weights Over Time", xaxis = list(title = "Time"), yaxis = list(title = "Weight") ) # Combine plots fig <- subplot(p1, p2, p3, nrows = 3, shareX = TRUE, titleY = TRUE) %>% layout( title = "Unified Kime-Phase Tomography Results", showlegend = TRUE ) return(fig) } # -------------------------------------------------------------- # Main Analysis Script # -------------------------------------------------------------- run_kpt_analysis <- function() { # Generate synthetic fMRI data set.seed(123) fmri_data <- generate_synthetic_fmri( N = 15, # 15 repetitions/subjects T = 600, # 600 timepoints tmax = 30, # 30 seconds mu_func = complex_mu, kappa_func = complex_kappa, noise_sd = 100, # 20 or 0.3, hrf_type = "double-gamma" ) # Define different parameter configurations to test config1 <- list(preprocess = TRUE, kernel_sigma = 0.5) # RKHS with moderate smoothing config2 <- list(preprocess = TRUE, kernel_sigma = 0.0) # No RKHS projection config3 <- list(preprocess = FALSE, kernel_sigma = 0.5) # No preprocessing, with RKHS # Apply KPT with different configurations cat("Running KPT with configuration 1 (RKHS + Preprocess)...\n") kpt1 <- unified_kpt(fmri_data$bold, fmri_data$t, preprocess = config1$preprocess, kernel_sigma = config1$kernel_sigma) cat("Running KPT with configuration 2 (No RKHS + Preprocess)...\n") kpt2 <- unified_kpt(fmri_data$bold, fmri_data$t, preprocess = config2$preprocess, kernel_sigma = config2$kernel_sigma) cat("Running KPT with configuration 3 (RKHS + No Preprocess)...\n") kpt3 <- unified_kpt(fmri_data$bold, fmri_data$t, preprocess = config3$preprocess, kernel_sigma = config3$kernel_sigma) # Evaluate performance metrics1 <- evaluate_kpt(kpt1, fmri_data) metrics2 <- evaluate_kpt(kpt2, fmri_data) metrics3 <- evaluate_kpt(kpt3, fmri_data) # Create comparison table performance_table <- data.frame( Configuration = c("RKHS + Preprocess", "No RKHS + Preprocess", "RKHS + No Preprocess"), Phase_CMAE = c(metrics1$phase_cmae, metrics2$phase_cmae, metrics3$phase_cmae), Kappa_RMSE = c(metrics1$kappa_rmse, metrics2$kappa_rmse, metrics3$kappa_rmse) ) cat("\nPerformance Comparison:\n") performance_table # Visualize results for the best configuration fig <- visualize_kpt_results(kpt1, fmri_data) fig # Analyze basis dominance around jumps # Extract basis weights for analysis basis_weights <- data.frame( t = fmri_data$t, K1_weight = kpt1$weights$w1, K2_weight = kpt1$weights$w2, K3_weight = kpt1$weights$w3 ) # Find time points where each basis dominates K1_dominant <- which(basis_weights$K1_weight > pmax(basis_weights$K2_weight, basis_weights$K3_weight)) K2_dominant <- which(basis_weights$K2_weight > pmax(basis_weights$K1_weight, basis_weights$K3_weight)) K3_dominant <- which(basis_weights$K3_weight > pmax(basis_weights$K1_weight, basis_weights$K2_weight)) # Detect true jumps true_jumps <- detect_phase_jumps(fmri_data$true_mean, fmri_data$t) # Analyze basis dominance around phase jumps jump_basis_dominance <- lapply(true_jumps$time, function(tj) { window <- 1.0 # 1-second window around jump indices <- which(fmri_data$t >= tj - window & fmri_data$t <= tj + window) # Count dominant bases in this region K1_count <- sum(indices %in% K1_dominant) K2_count <- sum(indices %in% K2_dominant) K3_count <- sum(indices %in% K3_dominant) return(data.frame( jump_time = tj, K1_dominant_count = K1_count, K2_dominant_count = K2_count, K3_dominant_count = K3_count, dominant_basis = which.max(c(K1_count, K2_count, K3_count)) )) }) # Combine results jump_basis_analysis <- do.call(rbind, jump_basis_dominance) cat("\nBasis Dominance Around Phase Jumps:\n") jump_basis_analysis # Return results for further analysis return(list( data = fmri_data, kpt1 = kpt1, kpt2 = kpt2, kpt3 = kpt3, fig = fig, metrics = list(metrics1 = metrics1, metrics2 = metrics2, metrics3 = metrics3), performance = performance_table, jump_analysis = jump_basis_analysis )) } # # Run analysis if script is executed directly # if (!exists("source_run")) { # results <- run_kpt_analysis() # } results <- run_kpt_analysis() ## First plot the MULTIPLE REPEATS of the simulated fMRI data # # Simulated fMRI series # fmri_simData <- generate_synthetic_fmri( # N = 15, # 15 repetitions/subjects # T = 600, # 600 timepoints # tmax = 30, # 30 seconds # mu_func = complex_mu, # kappa_func = complex_kappa, # noise_sd = 0.15, # hrf_type = "double-gamma" # ) # Create a color palette n_simulations <- 11 # Total number of simulations # color_palette <- colorRampPalette(brewer.pal(9, "Blues"))(n_simulations) # Option 1: Use a spectral color palette (rainbow-like) # color_palette <- colorRampPalette(brewer.pal(11, "Spectral"))(n_simulations) # # Option 2: Use a diverging color palette color_palette <- colorRampPalette(brewer.pal(11, "RdYlBu"))(n_simulations) # Initialize the plot with x values p <- plot_ly(x = results$data$t) # Add all traces in a single loop for (i in 1:n_simulations) { p <- add_trace(p, y = results$data$bold[i,], type = "scatter", mode = "lines", name = paste0("Sim fMRI (BOLD) Repeat ", i), line = list(color = color_palette[i], width = 0.5), opacity = 0.8, showlegend = TRUE) } # Add layout properties p <- p %>% layout( title = "Multiple Repeats (N=11) of fMRI BOLD Simulations", xaxis = list(title = "Time"), yaxis = list(title = "BOLD Signal"), legend = list(x = 1.02, y = 0.98, xanchor = "left", bordercolor = "#EEEEEE", borderwidth = 1), margin = list(r = 150), # Extra right margin for legend hovermode = "closest" ) # Display the plot p # Then plot the model results results$fig ``` ### Algorithm 4: Brain network fMRI Analytics Let's implement a complete example of *Kime-Phase Tomography (KPT)* applied to *synthetic fMRI data* with non-trivial properties, including phase jumps and time-varying dynamics. We demonstrate how the theoretical foundations translate into practical algorithms for recovering phase information from complex time-dependent signals. The experimental design involves analyzing synthetic fMRI BOLD signals with known ground truth to validate our methodology. The data incorporates realistic features of neuroimaging data - Multiple subjects/repetitions ($N=15$) - Hemodynamic response function convolution - Physiological noise - Time-varying phase dynamics - Abrupt phase transitions (simulating state changes). To complete this *unified Kime-Phase Tomography framework*, we will also consider a more *advanced and comprehensive example* that demonstrates the application to real-world fMRI data analysis problems. This experiment creates a visualization of synthetically generated fMRI data from multiple brain regions with known coupling patterns (ROI network interactions). 1. **Multi-Region Time Series**: BOLD signals from $5$ different brain regions across $5$ subjects (out of the $10$ generated). 2. **Hierarchical Coupling**: The data is generated with a specific *hierarchical ROI-connectivity pattern*. Region 1 drives regions 2 and 3 (coupling strengths $0.7$ and $0.5$). Region 2 drives region 4 (coupling strength $0.6$). And Region 3 drives region 5 (coupling strength $0.8$). 3. **Region-Specific Characteristics**: Each region has a different base oscillation frequency (ranging from $0.05\ Hz$ to $0.13\ Hz$). Regions that are coupled (e.g., Region 2 follows Region 1) show similar temporal patterns, but with individual variability. The coupling is implemented through phase synchronization, where the phase of a driven region is pulled toward the phase of the driving region. All signals incorporate the hemodynamic response function (HRF) typical of BOLD fMRI, which causes temporal smoothing and delays. Each subject has unique noise patterns, creating realistic inter-subject variability while preserving the underlying connectivity structure. And finally, the connectivity between regions is visualized as a heatmap, showing the directed influence from source to target regions. This shows how the KPT framework can be applied to detect phase relationships across brain regions in fMRI data, which is particularly relevant for studying functional connectivity in neuroimaging research. ```{r approach4_KPT_DemoAdvanced_fMRI, warning=FALSE, message=FALSE} # Visualization of fMRI data with 5 subjects library(plotly) library(dplyr) library(RColorBrewer) library(igraph) # Source necessary functions create_hrf <- function(t_hrf, type = "canonical") { if (type == "gamma") { hrf <- dgamma(t_hrf, shape = 6, scale = 1) } else if (type == "double-gamma") { hrf <- dgamma(t_hrf, shape = 6, scale = 1) - 0.1 * dgamma(t_hrf, shape = 16, scale = 1) } else { # canonical hrf <- (t_hrf/1.5)^3 * exp(-(t_hrf/1.5)) / gamma(4) } return(hrf / max(hrf)) # Normalize } # Generate multi-region fMRI data with high variability generate_high_variability_fmri <- function(regions = 5, N = 10, T = 500, noise_level = 0.5, phase_noise_level = 0.5, freq_variability = 0.05, amp_variability = 0.4) { # Create time vector t <- seq(0, 30, length.out = T) dt <- t[2] - t[1] # Create coupling matrix with hierarchical structure coupling_matrix <- matrix(0, regions, regions) # Region 1 drives regions 2 and 3 coupling_matrix[1, 2] <- 0.7 coupling_matrix[1, 3] <- 0.5 # Region 2 drives region 4 coupling_matrix[2, 4] <- 0.6 # Region 3 drives region 5 coupling_matrix[3, 5] <- 0.8 # Generate base phases for each region base_phases <- matrix(NA, regions, T) # Initialize with uncoupled phases - each region has its own base frequency for (r in 1:regions) { freq <- 0.05 + 0.02 * (r-1) # Base frequencies from 0.05 to 0.13 Hz base_phases[r,] <- 2 * pi * freq * t } # Apply coupling (phase synchronization) for (i in 1:5) { # Iterate several times to converge for (r1 in 1:regions) { for (r2 in 1:regions) { if (coupling_matrix[r1, r2] > 0) { # Phase of r2 is pulled toward phase of r1 phase_diff <- base_phases[r1,] - base_phases[r2,] # Wrap to [-pi, pi] phase_diff <- ((phase_diff + pi) %% (2*pi)) - pi # Apply coupling base_phases[r2,] <- base_phases[r2,] + coupling_matrix[r1, r2] * phase_diff } } } } # Generate BOLD signals with higher variability bold_signals <- array(NA, dim = c(N, regions, T)) # HRF for convolution t_hrf <- seq(0, 20, length.out = 100) hrf <- create_hrf(t_hrf, type = "double-gamma") # Subject-specific parameters to increase heterogeneity subject_params <- list() for (n in 1:N) { # Each subject has unique frequency offsets for each region freq_offsets <- rnorm(regions, 0, freq_variability) # Each subject has unique amplitude for each region amplitudes <- runif(regions, 1 - amp_variability, 1 + amp_variability) # Each subject has a unique HRF variation hrf_variation <- runif(1, 0.8, 1.2) # Some subjects might have additional low-frequency drift drift_amplitude <- runif(1, 0, 0.5) drift_frequency <- runif(1, 0.005, 0.015) subject_params[[n]] <- list( freq_offsets = freq_offsets, amplitudes = amplitudes, hrf_variation = hrf_variation, drift_amplitude = drift_amplitude, drift_frequency = drift_frequency ) } for (n in 1:N) { params <- subject_params[[n]] for (r in 1:regions) { # Apply subject-specific frequency variation modified_phase <- base_phases[r,] + cumsum(rep(params$freq_offsets[r], T)) * 2 * pi * dt # Add significant random phase noise phase_noise <- rnorm(T, 0, phase_noise_level) # Create complex signal with subject-specific amplitude neural_signal <- params$amplitudes[r] * exp(1i * (modified_phase + phase_noise)) # Modify HRF for this subject subject_hrf <- hrf * params$hrf_variation # Convolve with HRF bold <- Re(stats::convolve(neural_signal, rev(subject_hrf), type = "open"))[1:T] # Add subject-specific drift drift <- params$drift_amplitude * sin(2 * pi * params$drift_frequency * t) # Add higher measurement noise (lower SNR) bold_noisy <- bold + drift + rnorm(T, 0, noise_level) bold_signals[n, r, ] <- bold_noisy } } return(list( t = t, bold = bold_signals, true_phases = base_phases, coupling_matrix = coupling_matrix, subject_params = subject_params )) } # Generate the data set.seed(123) multiregion_data <- generate_high_variability_fmri( regions = 5, N = 10, T = 300, noise_level = 5.0, # high!!! # for low noise: noise_level = 0.5 phase_noise_level = 0.4, freq_variability = 0.05, amp_variability = 0.4 ) # Extract data for visualization t <- multiregion_data$t bold_signals <- multiregion_data$bold coupling_matrix <- multiregion_data$coupling_matrix true_phases <- multiregion_data$true_phases # Select 5 subjects subjects_to_plot <- c(1, 3, 5, 7, 9) # 1. Create fMRI time series visualization for 5 subjects timeseries_plot <- function() { # Create subplot for each region plots <- list() for (r in 1:5) { # Extract data for this region region_data <- data.frame() for (s_idx in seq_along(subjects_to_plot)) { s <- subjects_to_plot[s_idx] temp_df <- data.frame( time = t, bold = bold_signals[s, r, ], subject = paste("Subject", s) ) region_data <- rbind(region_data, temp_df) } # Create plot for this region p <- plot_ly(region_data, x = ~time, y = ~bold, color = ~subject, type = "scatter", mode = "lines") %>% layout( title = paste("Region", r), xaxis = list(title = ""), yaxis = list(title = "BOLD Signal"), showlegend = (r == 1) # Only show legend for first region ) plots[[r]] <- p } # Combine into a subplot subplot(plots, nrows = 5, shareX = TRUE, titleY = TRUE) %>% layout( title = "fMRI BOLD Signals (5 Selected Subjects)", xaxis = list(title = "Time (s)") ) } # 2. Create coupling matrix heatmap coupling_heatmap <- function() { plot_ly( x = paste("Region", 1:5), y = paste("Region", 1:5), z = t(coupling_matrix), # Transpose to match visualization convention type = "heatmap", colorscale = "Blues", showscale = TRUE, zmin = 0, zmax = 1 ) %>% layout( title = "Region Coupling Matrix (Source → Target)", xaxis = list(title = "Source Region"), yaxis = list(title = "Target Region") ) } # 3. Create network graph visualization # Fixed network graph function network_graph <- function() { # Create graph g <- graph_from_adjacency_matrix( coupling_matrix > 0, mode = "directed", weighted = TRUE ) # Set edge weights for visualization E(g)$width <- E(g)$weight * 5 # Generate layout layout <- layout_with_fr(g) # Create vectors for plotting edges <- get.edgelist(g) # Node coordinates node_x <- layout[,1] node_y <- layout[,2] # Edge coordinates edge_x <- edge_y <- list() for (i in 1:nrow(edges)) { from_idx <- as.numeric(edges[i, 1]) to_idx <- as.numeric(edges[i, 2]) edge_x[[i]] <- c(node_x[from_idx], node_x[to_idx], NA) edge_y[[i]] <- c(node_y[from_idx], node_y[to_idx], NA) } edge_x <- unlist(edge_x) edge_y <- unlist(edge_y) # Create plot fig <- plot_ly() %>% # Add edges add_trace( x = edge_x, y = edge_y, mode = "lines", line = list(width = 2, color = "rgba(150, 150, 150, 0.8)"), hoverinfo = "none", showlegend = FALSE ) %>% # Add nodes add_trace( x = node_x, y = node_y, mode = "markers+text", marker = list( size = 30, color = "rgba(31, 119, 180, 0.8)", line = list(width = 2, color = "rgb(0, 0, 0)") ), text = paste("R", 1:5), textposition = "middle center", textfont = list(color = "white", size = 14), hoverinfo = "text", hovertext = paste("Region", 1:5), showlegend = FALSE ) # Add edge labels for (i in 1:nrow(edges)) { from_idx <- as.numeric(edges[i, 1]) to_idx <- as.numeric(edges[i, 2]) # Find weight weight <- coupling_matrix[from_idx, to_idx] # Calculate midpoint mid_x <- (node_x[from_idx] + node_x[to_idx]) / 2 mid_y <- (node_y[from_idx] + node_y[to_idx]) / 2 # Add label fig <- fig %>% add_annotations( x = mid_x, y = mid_y, text = sprintf("%.1f", weight), showarrow = FALSE, font = list(size = 12) ) } # Complete the layout fig <- fig %>% layout( title = "Brain Region Connectivity Network", showlegend = FALSE, xaxis = list( title = "", showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE ), yaxis = list( title = "", showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE ) ) return(fig) } # Run the function with example data # For this example, we'll create a sample coupling matrix if not already defined if (!exists("coupling_matrix")) { coupling_matrix <- matrix(0, 5, 5) coupling_matrix[1, 2] <- 0.7 coupling_matrix[1, 3] <- 0.5 coupling_matrix[2, 4] <- 0.6 coupling_matrix[3, 5] <- 0.8 } # # Generate and display the network graph # network_plot <- network_graph() # network_plot # 4. Create phase visualization phase_plot <- function() { # Create data frame for phase plotting phase_plot_data <- data.frame() for (r in 1:5) { temp_df <- data.frame( time = t, phase = true_phases[r, ], region = paste("Region", r) ) phase_plot_data <- rbind(phase_plot_data, temp_df) } # Plot the phase patterns plot_ly(phase_plot_data, x = ~time, y = ~phase, color = ~region, type = "scatter", mode = "lines") %>% layout( title = "True Phase Patterns by Region", xaxis = list(title = "Time (s)"), yaxis = list(title = "Phase (radians)"), legend = list(title = list(text = "Region")) ) } # Generate all plots ts_plot <- timeseries_plot() heatmap_plot <- coupling_heatmap() network_plot <- network_graph() phases_plot <- phase_plot() # Display plots ts_plot heatmap_plot network_plot phases_plot ``` Now that we have the simulated fMRI data, we can explore *kime phase tomography.* ```{r approach4_KPT_DemoAdvanced2, warning=FALSE, message=FALSE} # Advanced Applications of Kime-Phase Tomography for Neuroimaging # This script extends the unified KPT framework to: # 1. Multi-region phase coupling analysis # 2. Non-parametric phase distribution modeling # 3. Event-related design applications # 4. Phase transition detection in cognitive tasks library(signal) library(circular) library(plotly) library(dplyr) library(igraph) # For network analysis # ----------------------------------------------------------------------- # 1. Multi-Region Phase Coupling Analysis # ----------------------------------------------------------------------- #' Generate synthetic multi-region fMRI data with known phase coupling #' #' @param regions Number of brain regions #' @param N Number of repetitions/subjects #' @param T Number of timepoints #' @param noise_level fMRI/BOLD signal noise level #' @param phase_noise_level kime-phase noise level #' @param freq_variability frequency variability #' @param amp_variability amplitude variability #' @param coupling_matrix Matrix specifying phase coupling between regions #' @return List containing generated data and ground truth generate_multiregion_fmri <- function(regions = 5, N = 10, T = 500, noise_level = 0.5, phase_noise_level = 0.5, freq_variability = 0.05, amp_variability = 0.4) { # Create time vector t <- seq(0, 30, length.out = T) dt <- t[2] - t[1] # Create coupling matrix with hierarchical structure coupling_matrix <- matrix(0, regions, regions) # Region 1 drives regions 2 and 3 coupling_matrix[1, 2] <- 0.7 coupling_matrix[1, 3] <- 0.5 # Region 2 drives region 4 coupling_matrix[2, 4] <- 0.6 # Region 3 drives region 5 coupling_matrix[3, 5] <- 0.8 # Generate base phases for each region base_phases <- matrix(NA, regions, T) # Initialize with uncoupled phases - each region has its own base frequency for (r in 1:regions) { freq <- 0.05 + 0.02 * (r-1) # Base frequencies from 0.05 to 0.13 Hz base_phases[r,] <- 2 * pi * freq * t } # Apply coupling (phase synchronization) for (i in 1:5) { # Iterate several times to converge for (r1 in 1:regions) { for (r2 in 1:regions) { if (coupling_matrix[r1, r2] > 0) { # Phase of r2 is pulled toward phase of r1 phase_diff <- base_phases[r1,] - base_phases[r2,] # Wrap to [-pi, pi] phase_diff <- ((phase_diff + pi) %% (2*pi)) - pi # Apply coupling base_phases[r2,] <- base_phases[r2,] + coupling_matrix[r1, r2] * phase_diff } } } } # Generate BOLD signals with higher variability bold_signals <- array(NA, dim = c(N, regions, T)) # HRF for convolution t_hrf <- seq(0, 20, length.out = 100) hrf <- create_hrf(t_hrf, type = "double-gamma") # Subject-specific parameters to increase heterogeneity subject_params <- list() for (n in 1:N) { # Each subject has unique frequency offsets for each region freq_offsets <- rnorm(regions, 0, freq_variability) # Each subject has unique amplitude for each region amplitudes <- runif(regions, 1 - amp_variability, 1 + amp_variability) # Each subject has a unique HRF variation hrf_variation <- runif(1, 0.8, 1.2) # Some subjects might have additional low-frequency drift drift_amplitude <- runif(1, 0, 0.5) drift_frequency <- runif(1, 0.005, 0.015) subject_params[[n]] <- list( freq_offsets = freq_offsets, amplitudes = amplitudes, hrf_variation = hrf_variation, drift_amplitude = drift_amplitude, drift_frequency = drift_frequency ) } for (n in 1:N) { params <- subject_params[[n]] for (r in 1:regions) { # Apply subject-specific frequency variation modified_phase <- base_phases[r,] + cumsum(rep(params$freq_offsets[r], T)) * 2 * pi * dt # Add significant random phase noise phase_noise <- rnorm(T, 0, phase_noise_level) # Create complex signal with subject-specific amplitude neural_signal <- params$amplitudes[r] * exp(1i * (modified_phase + phase_noise)) # Modify HRF for this subject subject_hrf <- hrf * params$hrf_variation # Convolve with HRF bold <- Re(stats::convolve(neural_signal, rev(subject_hrf), type = "open"))[1:T] # Add subject-specific drift drift <- params$drift_amplitude * sin(2 * pi * params$drift_frequency * t) # Add higher measurement noise (lower SNR) bold_noisy <- bold + drift + rnorm(T, 0, noise_level) bold_signals[n, r, ] <- bold_noisy } } return(list( t = t, bold = bold_signals, true_phases = base_phases, coupling_matrix = coupling_matrix, subject_params = subject_params )) } #' Calculate phase synchronization between two signals #' #' @param phi1 Phase time series of first signal #' @param phi2 Phase time series of second signal #' @return Phase locking value (0-1) phase_locking_value <- function(phi1, phi2) { # Compute phase difference phase_diff <- phi1 - phi2 # Complex phase difference z <- exp(1i * phase_diff) # Phase locking value plv <- abs(mean(z)) return(plv) } #' Apply KPT to multi-region fMRI data #' #' @param bold_data Multi-region BOLD data [subjects × regions × time] #' @param t Time vector #' @return List of KPT results for each region apply_multiregion_kpt <- function(bold_data, t) { N <- dim(bold_data)[1] # Number of subjects regions <- dim(bold_data)[2] # Number of regions # Apply KPT to each region region_results <- list() for (r in 1:regions) { # Extract BOLD data for this region bold_matrix <- matrix(bold_data[, r, ], nrow = N) # Apply KPT cat("Processing region", r, "...\n") kpt_result <- unified_kpt(bold_matrix, t, preprocess = TRUE, kernel_sigma = 0.5) region_results[[r]] <- kpt_result } # Calculate phase coupling between all region pairs plv_matrix <- matrix(0, regions, regions) for (r1 in 1:regions) { for (r2 in 1:regions) { if (r1 != r2) { plv <- phase_locking_value(region_results[[r1]]$mu, region_results[[r2]]$mu) plv_matrix[r1, r2] <- plv } } } return(list( region_results = region_results, plv_matrix = plv_matrix )) } #' Visualize multi-region phase synchronization #' #' @param multiregion_kpt Results from apply_multiregion_kpt #' @param threshold PLV threshold for visualization #' @return Plotly figure visualize_phase_connectivity <- function(multiregion_kpt, threshold = 0.4) { # Extract PLV matrix plv_matrix <- multiregion_kpt$plv_matrix regions <- nrow(plv_matrix) # Create network graph g <- graph_from_adjacency_matrix( plv_matrix > threshold, mode = "directed", weighted = TRUE ) # Set edge weights for visualization E(g)$width <- E(g)$weight * 5 # Layout layout <- layout_with_fr(g) # Create vectors for plotting edges <- get.edgelist(g) edge_weights <- E(g)$width # Node coordinates node_x <- layout[,1] node_y <- layout[,2] # Edge coordinates edge_x <- edge_y <- list() for (i in 1:nrow(edges)) { from_idx <- as.numeric(edges[i, 1]) to_idx <- as.numeric(edges[i, 2]) edge_x[[i]] <- c(node_x[from_idx], node_x[to_idx], NA) edge_y[[i]] <- c(node_y[from_idx], node_y[to_idx], NA) } edge_x <- unlist(edge_x) edge_y <- unlist(edge_y) # Create plot fig <- plot_ly() %>% # Add edges add_trace( x = edge_x, y = edge_y, mode = "lines", line = list(width = 1, color = "rgba(150, 150, 150, 0.8)"), hoverinfo = "none", showlegend = FALSE ) %>% # Add nodes add_trace( x = node_x, y = node_y, mode = "markers", marker = list( size = 15, color = "rgba(31, 119, 180, 0.8)", line = list(width = 2, color = "rgb(0, 0, 0)") ), text = paste("Region", 1:regions), hoverinfo = "text", showlegend = FALSE ) %>% layout( title = "Phase Synchronization Network", showlegend = FALSE, xaxis = list( title = "", showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE ), yaxis = list( title = "", showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE ) ) return(fig) } # ----------------------------------------------------------------------- # 2. Non-Parametric Phase Distribution Modeling # ----------------------------------------------------------------------- #' Estimate non-parametric phase distribution #' #' @param phases Matrix of phase observations [subjects × time] #' @param t Time vector #' @param bw Kernel bandwidth #' @return List with distribution estimates estimate_nonparametric_phase_dist <- function(phases, t, bw = 0.3) { T <- ncol(phases) N <- nrow(phases) # Prepare results density_estimates <- list() # For each time point, estimate circular density for (i in 1:T) { # Extract phases at this time point theta <- phases[, i] # Generate grid for density evaluation grid <- seq(-pi, pi, length.out = 100) # Compute von Mises kernel density density_vals <- numeric(length(grid)) for (j in seq_along(grid)) { # Sum of von Mises kernels centered at each observation kernel_sum <- sum(exp(bw * cos(grid[j] - theta))) density_vals[j] <- kernel_sum } # Normalize density_vals <- density_vals / (sum(density_vals) * (grid[2] - grid[1])) density_estimates[[i]] <- list( grid = grid, density = density_vals ) } return(list( t = t, density_estimates = density_estimates )) } #' Visualize non-parametric phase distribution over time #' #' @param nonparam_dist Results from estimate_nonparametric_phase_dist #' @param time_indices Which time points to visualize #' @return Plotly figure visualize_nonparametric_dist <- function(nonparam_dist, time_indices = NULL) { t <- nonparam_dist$t densities <- nonparam_dist$density_estimates # If time indices not specified, choose evenly spaced points if (is.null(time_indices)) { n_points <- 5 time_indices <- round(seq(1, length(t), length.out = n_points)) } # Create data frame for plotting plot_data <- data.frame() for (idx in time_indices) { temp_df <- data.frame( theta = densities[[idx]]$grid, density = densities[[idx]]$density, time = t[idx] ) plot_data <- rbind(plot_data, temp_df) } # Convert to factor for discrete color scale plot_data$time_factor <- factor(plot_data$time) # Create plot fig <- plot_ly(plot_data, x = ~theta, y = ~density, color = ~time_factor, type = "scatter", mode = "lines") %>% layout( title = "Non-Parametric Phase Distribution Over Time", xaxis = list(title = "Phase (radians)"), yaxis = list(title = "Density"), legend = list(title = list(text = "Time (s)")) ) return(fig) } # ----------------------------------------------------------------------- # 3. Event-Related Design Analysis # ----------------------------------------------------------------------- #' Generate event-related fMRI data with phase reset #' #' @param N Number of subjects #' @param T Number of timepoints #' @param tmax Maximum time #' @param event_times Vector of event times #' @param reset_magnitude Magnitude of phase reset #' @return List containing generated data generate_event_related_fmri <- function(N = 10, T = 600, tmax = 60, event_times = c(10, 30, 50), reset_magnitude = pi/2) { # Time vector t <- seq(0, tmax, length.out = T) dt <- t[2] - t[1] # Find indices of events event_indices <- sapply(event_times, function(et) which.min(abs(t - et))) # Generate signals bold_signals <- matrix(NA, nrow = N, ncol = T) true_phases <- matrix(NA, nrow = N, ncol = T) # HRF t_hrf <- seq(0, 20, length.out = 100) hrf <- create_hrf(t_hrf, type = "double-gamma") for (n in 1:N) { # Base frequency varies slightly by subject freq <- 0.08 + 0.01 * rnorm(1) # Initialize continuous phase phase <- numeric(T) phase[1] <- 2 * pi * runif(1) # Random initial phase # Evolve phase with resets at events for (i in 2:T) { # Check if this is an event onset is_event <- i %in% event_indices if (is_event) { # Apply phase reset reset_direction <- sample(c(-1, 1), 1) # Random direction phase[i] <- phase[i-1] + reset_direction * reset_magnitude } else { # Normal phase evolution phase[i] <- phase[i-1] + 2 * pi * freq * dt } } # Wrap to [-pi, pi] phase <- ((phase + pi) %% (2*pi)) - pi # Create neural signal neural_signal <- exp(1i * phase) # Add amplitude modulation at events amplitude <- rep(1, T) for (idx in event_indices) { # Amplitude rises after event window <- idx:(min(idx + 100, T)) amplitude[window] <- amplitude[window] + 0.5 * exp(-(1:length(window))/20) } # Apply amplitude neural_signal <- amplitude * neural_signal # Convolve with HRF bold <- Re(convolve(neural_signal, rev(hrf), type = "open"))[1:T] # Add noise bold_signals[n, ] <- bold + rnorm(T, 0, 0.15) true_phases[n, ] <- phase } return(list( t = t, bold = bold_signals, true_phases = true_phases, event_times = event_times )) } #' Analyze phase reset in event-related data #' #' @param kpt_result Result from unified_kpt #' @param event_times Vector of event times #' @param window_pre Pre-event window in seconds #' @param window_post Post-event window in seconds #' @return List with phase reset analysis analyze_phase_reset <- function(kpt_result, event_times, window_pre = 2, window_post = 5) { t <- kpt_result$t dt <- t[2] - t[1] # Find indices of events event_indices <- sapply(event_times, function(et) which.min(abs(t - et))) # Calculate window sizes in samples pre_samples <- round(window_pre / dt) post_samples <- round(window_post / dt) # For each event, extract phase before and after reset_analysis <- list() for (i in seq_along(event_indices)) { idx <- event_indices[i] # Extract time windows pre_idx <- max(1, idx - pre_samples):idx post_idx <- idx:(min(length(t), idx + post_samples)) # Calculate phase changes phase_pre <- kpt_result$mu[pre_idx] phase_post <- kpt_result$mu[post_idx] # Calculate phase derivative phase_deriv_pre <- c(0, diff(phase_pre)) phase_deriv_post <- c(0, diff(phase_post)) # Unwrap phase for consistency phase_pre_unwrap <- phase_pre for (j in 2:length(phase_pre_unwrap)) { diff_phase <- phase_pre_unwrap[j] - phase_pre_unwrap[j-1] if (abs(diff_phase) > pi) { phase_pre_unwrap[j:length(phase_pre_unwrap)] <- phase_pre_unwrap[j:length(phase_pre_unwrap)] - sign(diff_phase) * 2 * pi } } phase_post_unwrap <- phase_post for (j in 2:length(phase_post_unwrap)) { diff_phase <- phase_post_unwrap[j] - phase_post_unwrap[j-1] if (abs(diff_phase) > pi) { phase_post_unwrap[j:length(phase_post_unwrap)] <- phase_post_unwrap[j:length(phase_post_unwrap)] - sign(diff_phase) * 2 * pi } } # Calculate reset magnitude reset_magnitude <- phase_post_unwrap[1] - phase_pre_unwrap[length(phase_pre_unwrap)] # Calculate pre/post consistency pre_consistency <- abs(mean(exp(1i * phase_pre))) post_consistency <- abs(mean(exp(1i * phase_post))) # Store results reset_analysis[[i]] <- list( event_time = event_times[i], phase_pre = phase_pre, phase_post = phase_post, phase_deriv_pre = phase_deriv_pre, phase_deriv_post = phase_deriv_post, reset_magnitude = reset_magnitude, pre_consistency = pre_consistency, post_consistency = post_consistency ) } return(list( t = t, events = event_times, reset_analysis = reset_analysis )) } #' Visualize phase reset analysis #' #' @param reset_analysis Results from analyze_phase_reset #' @return Plotly figure visualize_phase_reset <- function(reset_analysis) { t <- reset_analysis$t dt <- t[2] - t[1] events <- reset_analysis$events # Create data frame for plotting plot_data <- data.frame() for (i in seq_along(reset_analysis$reset_analysis)) { event_data <- reset_analysis$reset_analysis[[i]] event_time <- event_data$event_time # Calculate relative time pre_time <- seq(-length(event_data$phase_pre) + 1, 0) * dt post_time <- seq(0, length(event_data$phase_post) - 1) * dt # Combine pre and post data time_rel <- c(pre_time, post_time[-1]) phase <- c(event_data$phase_pre, event_data$phase_post[-1]) phase_deriv <- c(event_data$phase_deriv_pre, event_data$phase_deriv_post[-1]) temp_df <- data.frame( time_rel = time_rel, phase = phase, phase_deriv = phase_deriv, event_id = factor(i) ) plot_data <- rbind(plot_data, temp_df) } # Create plot fig <- plot_ly() %>% add_trace(data = plot_data, x = ~time_rel, y = ~phase, color = ~event_id, type = "scatter", mode = "lines", name = "Phase", line = list(width = 2)) %>% add_trace(data = plot_data, x = ~time_rel, y = ~phase_deriv, color = ~event_id, type = "scatter", mode = "lines", name = "Phase Derivative", line = list(width = 1, dash = "dash"), visible = "legendonly") %>% add_trace(x = c(0, 0), y = c(-pi, pi), type = "scatter", mode = "lines", line = list(color = "black", width = 1, dash = "dot"), showlegend = FALSE) %>% layout( title = "Phase Reset Analysis Around Events", xaxis = list(title = "Time Relative to Event (s)"), yaxis = list(title = "Phase (radians)"), legend = list(title = list(text = "Event")) ) return(fig) } # ----------------------------------------------------------------------- # 4. Cognitive State Detection # ----------------------------------------------------------------------- #' Detect cognitive states from phase dynamics #' #' @param kpt_result Result from unified_kpt #' @param window_size Window size for state detection (seconds) #' @param k Number of states to identify #' @return List with state detection results detect_cognitive_states <- function(kpt_result, window_size = 3, k = 3) { t <- kpt_result$t dt <- t[2] - t[1] # Convert window size to samples window_samples <- round(window_size / dt) # Create feature matrix # For each window, extract: # 1. Mean phase derivative # 2. Phase variability # 3. Concentration parameter # Number of windows n_windows <- length(t) - window_samples + 1 # Feature matrix features <- matrix(0, n_windows, 3) for (i in 1:n_windows) { window_idx <- i:(i + window_samples - 1) phase_window <- kpt_result$mu[window_idx] kappa_window <- kpt_result$kappa[window_idx] # Unwrap phase for derivative calculation phase_unwrap <- phase_window for (j in 2:length(phase_unwrap)) { diff_phase <- phase_unwrap[j] - phase_unwrap[j-1] if (abs(diff_phase) > pi) { phase_unwrap[j:length(phase_unwrap)] <- phase_unwrap[j:length(phase_unwrap)] - sign(diff_phase) * 2 * pi } } # Calculate phase derivative phase_deriv <- diff(phase_unwrap) / dt # Features features[i, 1] <- mean(phase_deriv) # Mean phase velocity features[i, 2] <- sd(phase_deriv) # Variability in phase velocity features[i, 3] <- mean(kappa_window) # Mean concentration } # Scale features features_scaled <- scale(features) # K-means clustering kmeans_result <- kmeans(features_scaled, centers = k) # Assign states to each time point # For simplicity, we'll assign the state of the window starting at each point states <- rep(NA, length(t)) for (i in 1:n_windows) { states[i] <- kmeans_result$cluster[i] } # Fill in remaining points with the last state states[(n_windows+1):length(t)] <- states[n_windows] # Return results return(list( t = t, states = states, features = features, kmeans = kmeans_result )) } #' Visualize cognitive states #' #' @param state_detection Results from detect_cognitive_states #' @param kpt_result Result from unified_kpt #' @return Plotly figure visualize_cognitive_states <- function(state_detection, kpt_result) { t <- state_detection$t states <- state_detection$states # Create data frame df <- data.frame( t = t, phase = kpt_result$mu, concentration = kpt_result$kappa, state = factor(states) ) # Create color map for states state_colors <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd") # Create plot fig <- plot_ly() %>% add_trace(data = df, x = ~t, y = ~phase, color = ~state, type = "scatter", mode = "lines", line = list(width = 2, shape = "hv"), colors = state_colors[1:length(unique(states))]) %>% layout( title = "Cognitive State Detection from Phase Dynamics", xaxis = list(title = "Time (s)"), yaxis = list(title = "Phase (radians)"), legend = list(title = list(text = "State")) ) return(fig) } # ----------------------------------------------------------------------- # Run Advanced Analysis Example # ----------------------------------------------------------------------- run_advanced_analysis <- function() { # 1. Multi-region analysis cat("Generating multi-region fMRI data...\n") multiregion_data <- generate_multiregion_fmri(regions = 5, N = 10, T = 300, noise_level = 5.0) cat("Applying KPT to multi-region data...\n") multiregion_kpt <- apply_multiregion_kpt(multiregion_data$bold, multiregion_data$t) cat("Visualizing phase connectivity...\n") connectivity_plot <- visualize_phase_connectivity(multiregion_kpt, threshold = 0.4) connectivity_plot # 2. Non-parametric phase distribution cat("Estimating non-parametric phase distribution...\n") # Extract phases from first region phases_matrix <- matrix(0, nrow = 10, ncol = 300) for (n in 1:10) { phases_matrix[n,] <- multiregion_data$true_phases[1,] } nonparam_dist <- estimate_nonparametric_phase_dist(phases_matrix, multiregion_data$t) cat("Visualizing non-parametric distribution...\n") nonparam_plot <- visualize_nonparametric_dist(nonparam_dist) nonparam_plot # 3. Event-related analysis cat("Generating event-related fMRI data...\n") event_data <- generate_event_related_fmri(N = 15, T = 600, tmax = 60, event_times = c(10, 30, 50)) cat("Applying KPT to event-related data...\n") event_kpt <- unified_kpt(event_data$bold, event_data$t, preprocess = TRUE, kernel_sigma = 0.5) cat("Analyzing phase reset...\n") reset_analysis <- analyze_phase_reset(event_kpt, event_data$event_times) cat("Visualizing phase reset...\n") reset_plot <- visualize_phase_reset(reset_analysis) reset_plot # 4. Cognitive state detection cat("Detecting cognitive states...\n") state_detection <- detect_cognitive_states(event_kpt, window_size = 3, k = 3) cat("Visualizing cognitive states...\n") state_plot <- visualize_cognitive_states(state_detection, event_kpt) state_plot # Return results return(list( multiregion = list( data = multiregion_data, kpt = multiregion_kpt, plot = connectivity_plot ), nonparametric = list( dist = nonparam_dist, plot = nonparam_plot ), event_related = list( data = event_data, kpt = event_kpt, reset_analysis = reset_analysis, plot = reset_plot ), cognitive_states = list( detection = state_detection, plot = state_plot ) )) } # # Run if script is executed directly # if (!exists("source_run")) { # results <- run_advanced_analysis() # } #### Report and PLOT results .... cat("Generating multi-region fMRI data...\n") multiregion_data <- generate_multiregion_fmri(regions = 5, N = 10, T = 300) cat("Applying KPT to multi-region data...\n") multiregion_kpt <- apply_multiregion_kpt(multiregion_data$bold, multiregion_data$t) cat("Visualizing phase connectivity...\n") connectivity_plot <- visualize_phase_connectivity(multiregion_kpt, threshold = 0.4) connectivity_plot # 2. Non-parametric phase distribution cat("Estimating non-parametric phase distribution...\n") # Extract phases from first region phases_matrix <- matrix(0, nrow = 10, ncol = 300) for (n in 1:10) { phases_matrix[n,] <- multiregion_data$true_phases[1,] } nonparam_dist <- estimate_nonparametric_phase_dist(phases_matrix, multiregion_data$t) cat("Visualizing non-parametric distribution...\n") nonparam_plot <- visualize_nonparametric_dist(nonparam_dist) nonparam_plot # 3. Event-related analysis cat("Generating event-related fMRI data...\n") event_data <- generate_event_related_fmri(N = 15, T = 600, tmax = 60, event_times = c(10, 30, 50)) cat("Applying KPT to event-related data...\n") event_kpt <- unified_kpt(event_data$bold, event_data$t, preprocess = TRUE, kernel_sigma = 0.5) cat("Analyzing phase reset...\n") reset_analysis <- analyze_phase_reset(event_kpt, event_data$event_times) cat("Visualizing phase reset...\n") reset_plot <- visualize_phase_reset(reset_analysis) reset_plot # 4. Cognitive state detection cat("Detecting cognitive states...\n") state_detection <- detect_cognitive_states(event_kpt, window_size = 3, k = 3) cat("Visualizing cognitive states...\n") state_plot <- visualize_cognitive_states(state_detection, event_kpt) state_plot library(htmltools) str_output <- capture.output(str(event_kpt)) html_output <- tags$pre(paste(str_output, collapse = "\n")) html_output # str(event_kpt) ``` Notice the object structures of the result of the *KPT model* of the event-related fMRI data, which is an object `event_kpt <- unified_kpt(event_data$bold, event_data$t, preprocess = TRUE, kernel_sigma = 0.5)`.
SOCR Resource Visitor number Web Analytics SOCR Email