This Spacekime TCIU Learning Module presents the core elements of spacekime analytics including:

Import of repeated measurement longitudinal data,
Numeric (stitching) and analytic (Laplace) kimesurface reconstruction from time-series data,
Forward prediction modeling extrapolating the process behavior beyond the observed time-span $[0,T]$,
Group comparison discrimination between cohorts based on the structure and properties of their corresponding kimesurfaces. For instance, statistically quantify the differences between two or more groups,
Unsupervised clustering and classification of individuals, traits, and other latent characteristics of cases included in the study,
Construct low-dimensional visual representations of large repeated measurement data across multiple individuals as pooled kimesurfaces (parameterized 2D manifolds),
Statistical comparison, topological quantification, and analytical inference using kimesurface representations of repeated-measurement longitudinal data.

1 Preliminary setup

TCIU and other R package dependencies …

library(TCIU)
library(DT)
library("R.matlab")
# library(AnalyzeFMRI) # https://cran.r-project.org/web/packages/AnalyzeFMRI/index.html
library(ggfortify)
library(plotly)
library(manipulateWidget)
library(transport)
library(shapes)

2 Longitudinal Data Import

In this case, we are just loading some exemplary fMRI data, which is available here.

# pathname <- file.path("./test_data", "subj1_run1_complex_all.mat")
mat1 = readMat("./test_data/subj1_run1_complex_all.mat")
bigim1 = mat1$bigim[,64:1,,]
bigim1_mod = Mod(bigim1) # Modulus
smoothmod = GaussSmoothArray(bigim1_mod,sigma = diag(3,3))
# dim(bigim1) = 64 64 40
# bigim1 contains the complex image space
# dimensions are 64x*64y*40z*160t, corresponding to x,y,z,time
load("./test_data/mask.rda") # load the 3D nifti data of the mask
load("./test_data/mask_label.rda") # load the 3D nifti data of the mask with labels
load("./test_data/mask_dict.rda") # load the label index and label name
label_index = mask_dict$index
label_name = as.character(mask_dict$name)
label_mask = mask_label
load("./test_data/hemody.rda") # load the hemodynamic contour of the brain

3 Time-series graphs

3.1 Interactive time-series visualization

The TCIU function fmri_time_series() is used to create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data. We can also add a reference plotly object to the plot. This function is based on the GTSplot function from package TSplotly.

3.1.1 Example fMRI(x=4, y=42, z=33, t)

# reference_plot = sample[[4]]
fmri_time_series(bigim1, c(44,42,33), is.4d = TRUE) #, ref = reference_plot)

4 Kime-series/kimesurfaces (spacekime analytics protocol)

The following examples connect to several ongoing spacekime analytics, including kimesurface analyticity, kimesurface regularization, and distribution-time dynamics.

The first example uses a time-series simulation to illustrate how to transform the fMRI time-series at a fixed voxel location into a kimesurface (kime-series).

Notes:

Validate all steps in this time-series to kimesurfaces transformation protocol of simulated data, and finalize this 3D visualization.
Try it with real fMRI data at brain voxel locations associated with the finger-tapping task or musical genre study.

4.1 Pseudo-code

Randomly generate $8$ $phi=\phi$ kime-phases for each of the $10$ time radii. This yields an $8\times 10$ array (phi_8_vec) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$, (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc.
Optionally, order all kime-phases (rows) from small to large for each column.
Fix the $\nu=(x,y,z)$ voxel location and extract the fMRI time-course $fMRI_{\nu}(t), \forall 1\leq t\leq 160$.
For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into $80$ ON (Activation) and $80$ OFF (rest) states, or sub-series.
Construct a data-frame with $160$ rows and $4$ columns; time ($1:10$), phases ($8$), states ($2$), and fMRI_value (Complex or Real intensity).
Convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.
Check for visual (graphical) and numeric differences between the fMRI intensities during the ON vs. OFF states
Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have a large range.
Generate the plot_ly text labels that will be shown on mouse hover (pop-up dialogues) over each kimesurface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. The labels are $21\times21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kimesurface (ON vs. OFF stimuli).
Generate the $21\times21$ kime-domain Cartesian grid by polar-transforming the polar coordinates $(t,\phi)$ into Cartesian counterparts $(x,y)$.
Interpolate the fMRI intensities (natively anchored at $(t,\phi)$) to $fMRI(x,y), \forall -11\leq x,y \leq 11, \sqrt{x^2+y^2}\leq 10$.
Use plot_ly to display in 3D the kime-series as 2D manifolds (kimesurfaces) over the Cartesian domain.

4.2 Function main step: Time-series to kimesurfaces Mapping

4.2.1 Generate the kime-phases

Randomly generate $8$ $phi=\phi$ kime-phases for each of the $10$ time radii. This yields an $8\times 10$ array (phi_8_vec) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$; (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes; (3) use an $AR(p)$ autoregressive model with memory (affected by prior samples at earlier $p$ times), etc.

Optionally, order all kime-phases (rows) from small to large for each column.

# plot Laplacian
# ggfortify::ggdistribution(extraDistr::dlaplace, seq(-pi, pi, 0.01), m=0, s=0.5)

t <- seq(-pi, pi, 0.01)
f_t <- extraDistr::dlaplace(t, mu=0, sigma=0.5)
plot_ly(x=t, y=f_t, type="scatter", mode="lines", name="Laplace Distribution") %>%
  layout(title="Laplace Distribution: Kime-Phase Prior")

# randomly generate 8 phi kime-phases for each of the 10 time radii
phi_8_vec <- matrix(NA, ncol=10, nrow = 8)
for (t in 1:10) { 
  # for a given t, generate 8 new phases
  set.seed(t);
  phi_8_vec[ ,t] <-
    extraDistr::rlaplace(8, mu=0, sigma=0.5)
  # rank-order the phases for consistency
  # within the same foliation leaf
  phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t])
  # force phases in [-pi: pi)
  for (i in 1:8) {
    if (phi_8_vec[i,t] < -pi) 
      phi_8_vec[i,t] <- -pi
    if (phi_8_vec[i,t] >= pi) 
      phi_8_vec[i,t] <- pi
  }
}
# order all kime-phases (rows) from small to large for each column
# phi_8_vec_ordered <- apply(phi_8_vec, 2, sort)

4.2.2 Structural Data Preprocessing

Here should be included any study-specific data preprocessing. In the case of fMRI series, we may need to split off the individual repeated measurement fMRI time-series from the master (single) fMRI time-series according to the specific event-related fMRI design.

For simplicity, consider a simulated binary stimulus paradigm (e.g., activation (ON) and rest (OFF) event-related design). We can split the $160$-timepoint fMRI series into $80$ ON (Activation) and $80$ OFF (rest) states, or sub-series, consisting of $8$ repeats, each of length $10$ time points, where each time point corresponds to about $2.5\ s$ of clock time.

fMRI_ON<-bigim1_mod[40,42,33,][c(rep(TRUE,10),rep(FALSE,10))]
fMRI_OFF<-bigim1_mod[40,42,33,][c(rep(FALSE,10),rep(TRUE,10))]

4.2.3 Intensity Data Preprocessing

In practice, some spatial smoothing (blurring) the 1D time-series or their 2D array (tensor representations as 2D matrices ($row=repeat$, $column=time$) to reduce noise and make the data more natural (low-pass filtering, avoiding high-pitch noise). Sometimes, we may also need to temper the magnitude of the raw time-series intensities, which can have a large range.

# Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates
library(spatstat)

matrix_ON <- matrix(0, nrow = 21, ncol = 21) 
matrix_OFF <- matrix(0, nrow = 21, ncol = 21) 
for (t in 1:10) {
  for (p in 1:8) {
    x = 11+t*cos(phi_8_vec[p,t])
    y = 11+t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- fMRI_ON[(p-1)*10 +t] # What is matrix_ON?
    matrix_OFF[x,y] <- fMRI_OFF[(p-1)*10 +t]
  }
}
# smooth/blur the matrices
# matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
# matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))
matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))

4.2.4 Generate `plotly` labels

Generate the plot_ly text labels that will be shown over the graph, upon mouse hovering (pop-up dialogues) over each kimesurface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. The labels are $21\times 21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kimesurface (ON vs. OFF stimuli).

# fix the plot_ly Text Labels
x <- vector()
y <- vector()
i <- 1
for (t in 1:10) {
  for (p in 1:8) {
    x[i] = 11+t*cos(phi_8_vec[p,t])
    y[i] = 11+t*sin(phi_8_vec[p,t])
    i <- i+1
  }
}

# hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext)
# custom_txt <- matrix(NA, nrow=21, ncol=21)
# hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext)
# custom_txtOFF <- matrix(NA, nrow=21, ncol=21)
# 
# for (x in 1:21) {
#    for (y in 1:21) {
#      t = sqrt((x-11)^2 + (y-11)^2)
#      p = atan2(y-11, x-11)
#      custom_txt[x,y] <- paste(' fMRI: ', round(hoverText[(x-1)*21+y, 3], 3),
#                     '\n time: ', round(t, 0),
#                     '\n phi: ', round(p, 2)) 
#      custom_txtOFF[x,y] <- paste(' fMRI: ', round(hoverTextOFF[(x-1)*21+y, 3], 3),
#                     '\n time: ', round(t, 0),
#                     '\n phi: ', round(p, 2)) 
#    }
# }

4.2.5 Cartesian space interpolation

Interpolate the fMRI intensities, natively anchored at polar (kime) coordinates) $\left (\underbrace{t}_{time},\underbrace{\phi}_{repeat} \right)$, into Cartesian coordinates $fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10$. Note that for specificity, we hard-coded polar-grid parameterization to time $t\in\{1,2,3, \cdots,10\},\ |T|=10$ and phase $\phi=seq\left (from=-\pi, to=\pi, step=\frac{2\pi}{20}\right )\in\{-3.1415927, -2.8274334, \cdots, 0.0,\cdots, 2.8274334, 3.1415927 \},\ |\Phi|=21$.

# xx2 <- 11 + c(-10:10) %o% cos(seq(-pi, pi, 2*pi/20))
# yy2 <- 11 + c(-10:10)  %o% sin(seq(-pi, pi, 2*pi/20))
xx2 <- 11 + seq(0,10,1/2) %o% cos(seq(-pi, pi, 2*pi/20))
yy2 <- 11 + seq(0,10,1/2)  %o% sin(seq(-pi, pi, 2*pi/20))
#zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21)
zz2 <- matrix_ON_smooth
ww2 <- matrix_OFF_smooth
dd2 <- matrix_ON_smooth-matrix_OFF_smooth

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
f <- list(family = "Courier New, monospace", size = 18, color = "black")
x <- list(title = "k1", titlefont = f)
y <- list(title = "k2", titlefont = f)
z <- list(title = "fMRI Kime-series", titlefont = f)

# matrix_ON_smooth,matrix_OFF_smooth is indexed in Cartesian coordinates
# Reinterpolate the fMRI intensities according to Polar grid
library(akima)
ON_transformed <- matrix(0, nrow = 21, ncol = 21)
OFF_transformed <- matrix(0, nrow = 21, ncol = 21)
ON_OFF_transformed <- matrix(0, nrow = 21, ncol = 21)

cart_x <- as.vector(rep(seq(1,21),21))
cart_y <- as.vector(sort(rep(seq(1,21),21)))
cart_z_on <- as.vector(matrix_ON_smooth)
cart_z_diff <- as.vector(matrix_ON_smooth-matrix_OFF_smooth)
cart_z_off <- as.vector(matrix_OFF_smooth)

for(i in 1:21){
  for(j in 1:21){
    #interpolations
    int_res_on <- interp(cart_x,cart_y,cart_z_on,xx2[i,j],yy2[i,j])
    int_res_off <- interp(cart_x,cart_y,cart_z_off,xx2[i,j],yy2[i,j])
    int_res_diff <- interp(cart_x,cart_y,cart_z_diff,xx2[i,j],yy2[i,j])
    #insert data
    ON_transformed[i,j] = as.numeric(int_res_on$z)
    OFF_transformed[i,j] = as.numeric(int_res_off$z)
    ON_OFF_transformed[i,j] = as.numeric(int_res_diff$z)
    # if None set to 0
    if(is.na(ON_transformed[i,j])){
      ON_transformed[i,j] = 0
    }
    if(is.na(ON_OFF_transformed[i,j])){
      ON_OFF_transformed[i,j] = 0
    }
    if(is.na(OFF_transformed[i,j])){
      ON_OFF_transformed[i,j] = 0
    }
  }
}
# plot_ly(x = xx2,y = yy2,z = ON_transformed) %>% add_surface()

4.2.6 Cartesian representation

Generate the $21\times 21$ kime-domain Cartesian grid hovering texts.

# hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext)
custom_txt <- matrix(NA, nrow=21, ncol=21)
# hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext)
custom_txtOFF <- matrix(NA, nrow=21, ncol=21)
custom_txt_DIFF <- matrix(NA, nrow=21, ncol=21)
custom_txt_ON_OFF <- matrix(NA, nrow=21, ncol=21)

for (xdir in 1:21) {
   for (ydir in 1:21) {
     custom_txt[xdir,ydir] <- paste(' fMRI: ', round(ON_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txtOFF[xdir,ydir] <- paste(' fMRI: ', round(OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txt_DIFF[xdir,ydir] <- paste(' fMRI: ', round(ON_transformed[xdir,ydir]-OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txt_ON_OFF[xdir,ydir] <- paste(' fMRI: ', round(ON_OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
   }
}

plot_ly(x = xx2,y = yy2,z = ON_transformed, type = "surface", colors=c("#FFFFFF","#0000FF"),
          text = custom_txt, hoverinfo = "text", showlegend = FALSE) %>%
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines",
              line = list(width = 10, color="red"), name="Space(x)",
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", title = "ON kimesurface/Kime-Series at a fixed voxel location",scene = list(xaxis = x, yaxis = y, zaxis = z), showlegend = FALSE)

plot_ly(x = ~xx2, y = ~yy2, z = ~OFF_transformed, type = "surface",colors=c("#FFFFFF","#0000FF"),   # scatterpolar
          text = custom_txtOFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", title = "OFF kimesurface/Kime-Series at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = z), showlegend = FALSE)

# x <- as.vector(rep(seq(0.5,20.5),21))
# y <- as.vector(sort(rep(seq(0.5,20.5),21)))
# z <- as.vector(matrix_ON_smooth)
# intper <- interp(x,y,z)
# plot_ly(x = intper$x,y = intper$y,z = matrix_ON_smooth) %>% add_surface()

plot_ly(x = ~xx2, y = ~yy2, z = ~ON_transformed-OFF_transformed, type = "surface",colors=c("#FFFF00","#FFFFFF","#0000FF"),   # scatterpolar
          text = custom_txt_DIFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", title = "Raw ON-OFF difference kimesurface/Kime-Series at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = z), showlegend = FALSE)

# x <- as.vector(rep(seq(0.5,20.5),21))

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
zd <- list(title = "fMRI Kime-ON/OFF difference", titlefont = f)
dd2 <- matrix_ON_smooth-matrix_OFF_smooth
# dd2scale<-fmri_split_ab_bl(dd2)
plot_ly(x = ~xx2, y = ~yy2, z = ~ON_OFF_transformed, type = "surface", #colors = dd2scale, 
        colors=c("#FFFF00","#FFFFFF","#0000FF"),  # scatterpolar
          text = custom_txt_ON_OFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=-0.15:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", title = "Difference for ON & OFF kimesurface/Kime-Series at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = zd), showlegend = FALSE)

The last two kimesurface plots are almost identical as the interpolation process is fairly robust.

4.2.7 Generate a long data-frame

Convert the entire dataset into a long data-frame, i.e., construct a data-frame (DF) with $\underbrace{160}_{total}=\underbrace{2}_{binary\\ On/Off}\times \underbrace{8}_{repeats}\times \underbrace{10}_{timepoints}$ rows and $4$ columns representing time ($1:10$), phases ($8$), states ($2$), and fMRI_value (Complex or Real time-series intensity). Then, transform the long DF representing the fMRI_ON and fMRI_OFF time-sources from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.

# datatable(fmri_kimesurface(bigim1_mod,c(44,42,33))[[1]])

4.2.8 Display kimesurfaces

Use plot_ly to display in 3D the kime-series as 2D manifolds (kimesurfaces) over the Cartesian domain. In this specific binary case-study, we demonstrate the following $3$ kimesurface reconstructions:

the On kimesurface
the Off kimesurface, and
the difference On - Off kimesurface.

# fmri_kimesurface(bigim1_mod,c(44,42,33))[[2]]  # On kimesurface
# fmri_kimesurface(bigim1_mod,c(44,42,33))[[3]]  # Off kimesurface
# fmri_kimesurface(bigim1_mod,c(44,42,33))[[4]]  # Difference On - Off kimesurface

4.3 Autoregressive Kime-Phase Modeling

Above we demonstrated transforming repeated sampling longitudinal processes (time-series) as 2D manifolds (kimesurfaces) using a kime-phase prior distribution, such as the truncated Laplace phase distribution.

An alternative numerical protocols utilizes autoregressive model, $AR(p)$. For instance, here is an $AR(1)$ kime-phase model $\hat{\theta}_{t}=\alpha_1 \theta_{t-1}+ \varepsilon _{t}$, where $\varepsilon_{t}\sim N(0,\sigma_{\varepsilon}^2)$ (white noise) or $\varepsilon_{t}\sim \operatorname {Laplace}(\mu=0,scale=\sigma_{\varepsilon}^2)$ (Laplace phase), $\theta_{t}\equiv \hat{\theta}_{t} {\text{ mod }} 2\pi \sim \Phi_{[-\pi,\pi)}(t)$, and the (time-dependent) kime-phase distribution $\Phi_{[-\pi,\pi)}(t)$ is compactly supported over the (periodic) interval $[-\pi,\pi)$.

When $|\alpha_1 |<1$,This $AR(1)$ model is weak-stationary and represents the output of a stable filter with white noise input. In this case, $\forall\ t$, the mean $\mu=\mathbb{E} (\theta_{t})$ is constant, due to the model’s weak stationarity. Since, $\mathbb{E} (\theta_{t})=\alpha_1 \mathbb{E} (\theta_{t-1}) + \mathbb{E} (\varepsilon_{t}),$ $\Longrightarrow \mu =\alpha \mu +0,$. therefore, the (stationary) mean is trivial, $\mu =0$, and the (stationary) variance is constant, ${\text{var}}(\theta_{t})=\mathbb{E} (\theta_{t}^{2}) - \mu ^{2}={\frac {\sigma _{\varepsilon}^{2}}{1-\alpha ^{2}}}.$ This stationarity is due to the fact that ${\text {var}}(\theta_{t})=\alpha ^{2}{\text{var}}(\theta_{t-1})+\sigma_{\varepsilon }^{2}$.
However, if $\alpha =1$, the variance of $\theta_t$ depends on time lag $t$. Hence, the variance of the series diverges to infinity as $t\to\infty$, implying that the model is non-stationary.

In the more general case, the $AR(p)$ autoregressive model is \[{\hat{\theta}}_{t}= \sum_{i=1}^{p}\alpha_{i} \theta_{t-i} + \varepsilon _{t}, \ \ \theta_{t}\equiv \hat{\theta}_{t} {\text{ mod }} 2\pi \sim \Phi_{[-\pi,\pi)}(t)\]

where the coefficients $\alpha_{i}, \ i\in\{1, 2, \cdots, p\}$}} are hyper-parameters that can be set, estimated, or optimized to control the covariance function of the model. The $AR(p)$ autoregressive model may also be expressed using the backshift operator $B$

\[\hat{\theta}_{t} = \sum_{i=1}^{p}\alpha_{i}B^{i}\theta_{t} + \varepsilon _{t}, \ \ \theta_{t}\equiv \hat{\theta}_{t} {\text{ mod }} 2\pi \sim \Phi_{[-\pi,\pi)}(t).\]

Rearranging the terms on both hands sides, $\alpha [B]X_{t} = \varepsilon_{t}$, where $\alpha(z)\equiv 1 - \sum_{i=1}^{p} \alpha_{i}z^{i}$, each $z_{i}\in\mathbb{C}$ satisfies $|z_{i}|>1$, and the backshift operator $B^{i}\theta_{t}=\theta_{t-i},\ \forall\ i\in\mathbb{N}$.

Let’s run the same fMRI kimesurface reconstruction we did above (using the Laplace kime-phase prior distribution model) on the $AR(p=1)$ kime-phase temporal dynamics $t\in\{1,2,3,\cdots,10\}$, and visually compare the fMRI On, Off, and On-Off difference kimesurface reconstructions. As expected these kimesurafces look smoother and analytic.

These examples connect to several ongoing spacekime analytics R&D projects, including kimesurface analyticity, kimesurface regularization, and distribution-time dynamics.

# AR(1) model for generating 8 phi kime-phases for each of the 10 time radii
# phi_8_vec_AR1 <- matrix(NA, ncol=10, nrow = 8)
library(tidyr)

# Initialize all the phases: all repeats: 1<=i<=8, and times: 1<=t<=10
for (t in 1:10) { 
    phi_8_vec[ , t] <- extraDistr::rlaplace(8, mu=0, sigma=0.5) 
}

# AR(1) Modeling for all 8 fMRI repeats; for ARIMA(p,q,r), update order=c(p,q,r)
phi_8_vec_AR1 <- list()
for (i in 1:8) {
    phi_8_vec_AR1[[i]] <- arima(phi_8_vec[i , ], order=c(1,0,0))
}

# AI modeling and forcing all phases in [-pi: pi), modulo or via wrap around correction
for (i in 1:8) { 
  for (t in 2:10) { # fitting AR(p=1) model starts with t=1+p=2 
    # Use AR(1) models to adjust the phases
    phi_8_vec[i, t] <- phi_8_vec_AR1[[i]]$coef[2] + 
      phi_8_vec_AR1[[i]]$coef[1] * (phi_8_vec[i, t-1] - phi_8_vec_AR1[[i]]$coef[2]) 
      + phi_8_vec_AR1[[i]]$residuals[t]  ## Should the random error (residuals) be added???
    
    # force all phases in [-pi: pi), e.g., phi_8_vec_AR1[i,t] %% pi  (wrap around)
    if (phi_8_vec[i,t] < -pi) 
      phi_8_vec[i,t] <- phi_8_vec[i,t] %% pi
    if (phi_8_vec[i,t] >= pi) 
      phi_8_vec[i,t] <- ((phi_8_vec[i,t] + pi) %% pi) -pi
  }
}

phi_8_vec_AR1_df <- as.data.frame(t(phi_8_vec))
phi_8_vec_AR1_df$time <- c(1:10)
phi_8_vec_AR1_df_long <- phi_8_vec_AR1_df %>%
  pivot_longer(!time, names_to = "repeats", values_to = "values")
  
plot_ly(data = phi_8_vec_AR1_df_long,x = ~time, y=~values, type="scatter", mode="lines", 
        color = ~repeats, name= ~repeats) %>%
  layout(title="AR(1) Kime-Phase Prior Model based on Laplace Noise", 
         legend = list(title="Repeats", orientation = 'h', y=-0.2))

In the example above, the $AR(1)$ model for a fixed fMRI sequence repeat index $r$ is

\[fMRI_t(r) - \mu(r) = \alpha_1(r)(fMRI_{t-1}(r) - \mu(r)) + \varepsilon_t(r)\ .\] In this experiment, for the first fMRI repeat index, $r=1$, the estimated $AR(p=1)$ model coefficients and residuals were:

Slope, ar1, $\alpha_1(r)=$ 0.2878037, and
Mean, intercept, $\mu(r)=$ 0.359729.
The corresponding $AR(1)$ model residuals were $\{\varepsilon_t(r)\}_{t=1}^{10} =\{$ -0.5487502, -0.0093747, 0.8677366, 1.2673865, -0.1304475, 0.3686666, -1.3133557, 0.2258199, -0.3095218, -0.2290063 $\}$.

Observation: As expected, the $AR(p=1)$ (longitudinal) model of the kime-phases stabilizes with time. This may be important when kimesurface analyticity, smoothness, or other properties of the time-series to manifold transformations may be desirable. Let’s also illustrate the resulting On, Off, and On-Off difference kimesurface reconstructions. These ($AR(p=1)$ regularized) kimesurfaces can be contrasted against their earlier counterparts, which used the raw Laplace kime-phase samples without $AR(1)$ modeling.

# Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates
library(spatstat)

matrix_ON <- matrix(0, nrow = 21, ncol = 21) 
matrix_OFF <- matrix(0, nrow = 21, ncol = 21) 
for (t in 1:10) {
  for (p in 1:8) {
    x = 11+t*cos(phi_8_vec[p,t])
    y = 11+t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- fMRI_ON[(p-1)*10 +t] # What is matrix_ON?
    matrix_OFF[x,y] <- fMRI_OFF[(p-1)*10 +t]
  }
}
# smooth/blur the matrices
matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))

4.3.1 Generate `plotly` labels

# fix the plot_ly Text Labels
x <- vector()
y <- vector()
i <- 1
for (t in 1:10) {
  for (p in 1:8) {
    x[i] = 11+t*cos(phi_8_vec[p,t])
    y[i] = 11+t*sin(phi_8_vec[p,t])
    i <- i+1
  }
}

4.3.2 Cartesian space interpolation

# xx2 <- 11 + c(-10:10) %o% cos(seq(-pi, pi, 2*pi/20))
# yy2 <- 11 + c(-10:10)  %o% sin(seq(-pi, pi, 2*pi/20))
xx2 <- 11 + seq(0,10,1/2) %o% cos(seq(-pi, pi, 2*pi/20))
yy2 <- 11 + seq(0,10,1/2)  %o% sin(seq(-pi, pi, 2*pi/20))
#zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21)
zz2 <- matrix_ON_smooth
ww2 <- matrix_OFF_smooth
dd2 <- matrix_ON_smooth-matrix_OFF_smooth

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
f <- list(family = "Courier New, monospace", size = 18, color = "black")
x <- list(title = "k1", titlefont = f)
y <- list(title = "k2", titlefont = f)
z <- list(title = "AR(1) fMRI Kime-series", titlefont = f)

# matrix_ON_smooth,matrix_OFF_smooth is indexed in Cartesian coordinates
# Reinterpolate the fMRI intensities according to Polar grid
library(akima)
ON_transformed <- matrix(0, nrow = 21, ncol = 21)
OFF_transformed <- matrix(0, nrow = 21, ncol = 21)
ON_OFF_transformed <- matrix(0, nrow = 21, ncol = 21)

cart_x <- as.vector(rep(seq(1,21),21))
cart_y <- as.vector(sort(rep(seq(1,21),21)))
cart_z_on <- as.vector(matrix_ON_smooth)
cart_z_diff <- as.vector(matrix_ON_smooth-matrix_OFF_smooth)
cart_z_off <- as.vector(matrix_OFF_smooth)

for(i in 1:21){
  for(j in 1:21){
    #interpolations
    int_res_on <- interp(cart_x,cart_y,cart_z_on,xx2[i,j],yy2[i,j])
    int_res_off <- interp(cart_x,cart_y,cart_z_off,xx2[i,j],yy2[i,j])
    int_res_diff <- interp(cart_x,cart_y,cart_z_diff,xx2[i,j],yy2[i,j])
    #insert data
    ON_transformed[i,j] = as.numeric(int_res_on$z)
    OFF_transformed[i,j] = as.numeric(int_res_off$z)
    ON_OFF_transformed[i,j] = as.numeric(int_res_diff$z)
    # if None set to 0
    if(is.na(ON_transformed[i,j])){
      ON_transformed[i,j] = 0
    }
    if(is.na(ON_OFF_transformed[i,j])){
      ON_OFF_transformed[i,j] = 0
    }
    if(is.na(OFF_transformed[i,j])){
      ON_OFF_transformed[i,j] = 0
    }
  }
}
# plot_ly(x = xx2,y = yy2,z = ON_transformed) %>% add_surface()

4.3.3 Cartesian representation

Generate the $21\times 21$ kime-domain Cartesian grid hovering texts.

# hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext)
custom_txt <- matrix(NA, nrow=21, ncol=21)
# hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext)
custom_txtOFF <- matrix(NA, nrow=21, ncol=21)
custom_txt_DIFF <- matrix(NA, nrow=21, ncol=21)
custom_txt_ON_OFF <- matrix(NA, nrow=21, ncol=21)

for (xdir in 1:21) {
   for (ydir in 1:21) {
     custom_txt[xdir,ydir] <- paste(' AR(p=1) fMRI: ', round(ON_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txtOFF[xdir,ydir] <- paste(' fMRI: ', round(OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txt_DIFF[xdir,ydir] <- paste(' fMRI: ', round(ON_transformed[xdir,ydir]-OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
     custom_txt_ON_OFF[xdir,ydir] <- paste(' fMRI: ', round(ON_OFF_transformed[xdir,ydir], 3),
                    '\n time: ', round((xdir-1)/2, 0),
                    '\n phi: ', round(-pi+pi/10*(ydir-1), 2))
   }
}

plot_ly(x = xx2,y = yy2,z = ON_transformed, type = "surface", colors=c("#FFFFFF","#0000FF"),
          text = custom_txt, hoverinfo = "text", showlegend = FALSE) %>%
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines",
              line = list(width = 10, color="red"), name="Space(x)",
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode="turntable", 
           title="Kime-Phase AR(1) Model:\n ON kimesurface at a fixed voxel location",
           scene=list(xaxis=x, yaxis = y, zaxis = z)) %>% hide_colorbar()

plot_ly(x = ~xx2, y = ~yy2, z = ~OFF_transformed, type = "surface",colors=c("#FFFFFF","#0000FF"),   # scatterpolar
          text = custom_txtOFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", 
           title = "Kime-Phase AR(1) Model: OFF kimesurface at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = z)) %>% hide_colorbar()

plot_ly(x = ~xx2, y = ~yy2, z = ~ON_transformed-OFF_transformed, type = "surface",colors=c("#FFFF00","#FFFFFF","#0000FF"),   # scatterpolar
          text = custom_txt_DIFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=0:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", 
           title = "Kime-Phase AR(1) Model:\n Raw ON-OFF difference kimesurface at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = z)) %>% hide_colorbar()

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
zd <- list(title = "AR(1) fMRI Kime-ON/OFF difference", titlefont = f)
dd2 <- matrix_ON_smooth-matrix_OFF_smooth
# dd2scale<-fmri_split_ab_bl(dd2)
plot_ly(x = ~xx2, y = ~yy2, z = ~ON_OFF_transformed, type = "surface", #colors = dd2scale, 
        colors=c("#FFFF00","#FFFFFF","#0000FF"),  # scatterpolar
          text = custom_txt_ON_OFF, hoverinfo = "text", showlegend = FALSE) %>% 
    #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"),
    #          type="surface", text = custom_txtOFF, hoverinfo = "text",
    #          opacity=0.3, showscale = FALSE, showlegend = FALSE) %>%
    # trace the main Z-axis
    add_trace(x=11, y=11, z=-0.15:0.15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Space(x)", 
              hoverinfo="none", showlegend = FALSE) %>%
    layout(dragmode = "turntable", 
           title = "Kime-Phase AR(1) Model:\n Difference for ON & OFF kimesurfaces at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = zd)) %>% hide_colorbar()

4.4 Analytical Time-series to kimesurface Transformation (Laplace)

4.4.1 Data Preparation

First plot the simulated On (stimulus) and Off (rest) fMRI time-series at a fixed voxel location $(44,42,33)\in\mathbb{R}^3$, along with the averaged (pooled) On and Off signal over all $8$ repeats in the single run (epoch) of $160$ time-points.

# Extract just the On kime-series at voxel (44,42,33), each time-series has 8 repeats!
onFMRISeries <- 
  bigim1_mod[44,42,33, c(1:10, 21:30, 41:50, 61:70, 81:90, 101:110, 121:130, 141:150)]
# the corresponding Off kime-series at voxel (44,42,33) will be the temporal complement
offFMRISeries <- 
  bigim1_mod[44,42,33, -c(1:10, 21:30, 41:50, 61:70, 81:90, 101:110, 121:130, 141:150)]

t_indx <- seq(1, 80, 1) # index of time (not actual time value)
f_On <- onFMRISeries
f_Off <- offFMRISeries

vline <- function(x=0, color = "lightgray") {
  list(type="line", y0=0, y1=1, yref="paper", name="time break points", 
       opacity = 0.5, x0=x, x1=x, line=list(color=color))
}

hline <- function(y = 0, color = "blue") {
  list(type="line", x0=0, x1=1, xref="paper", name="intensity break points",
       opacity = 0.3, y0=y, y1=y, line=list(color=color))
}

plot_ly() %>%
  # On fMRI time-series
  add_trace(x=t_indx, y=f_On, type="scatter", mode="lines", name="On time-series") %>%
  # On fMRI time-series
  add_trace(x=t_indx, y=f_Off, type="scatter", mode="lines", name="Off time-series") %>%
  # Repeated measurement break points
  # add_trace(x=c(10,20,30,40,50,60,70,80), y=f_Off, type="scatter", mode="lines", name="Off time-series") %>%
  layout(title="3D fMRI Simulation: On & Off Time-series at Voxel(44,42,33)",
         shapes = list(
           vline(10),vline(20),vline(30),vline(40),vline(50),vline(60),vline(70),vline(80)),
         legend = list(orientation='h', y=-0.2))

# Compute and plot against each other the Average On and Average Off time-series
seriesAvg <- function(f=0, period=0) {
  tempAvg <- rep(0,period)  # initialize avg signal
  for (i in c(1:period)) {
    counter =0  #  empirically count the number of repeated samples (here it's 8)
    for (j in c(1:length(f))) {
      if (j %% period == i-1) {  # e.g., 159 %% 10 # [1] 9
        tempAvg[i] = tempAvg[i] + f[j]
        counter = counter + 1
      }
    }
    tempAvg[i] = tempAvg[i]/counter
  }
  return(tempAvg)
}

period <- 10
onAvg  <- seriesAvg(f=f_On,  period=period) 
offAvg <- seriesAvg(f=f_Off, period=period) 

plot_ly() %>%
  # Average On fMRI time-series
  add_trace(x=c(1:period), y=onAvg, type="scatter", mode="lines", name="On Average") %>%
  # On fMRI time-series
  add_trace(x=c(1:period), y=offAvg, type="scatter", mode="lines", name="Off Average") %>%
  # Repeated measurement break points
  layout(title=
  "3D fMRI Simulation: On & Off Time-series\n (averaged over all 8 repeats) at Voxel (44,42,33)",
  legend = list(orientation='h', y=-0.2))

Next, we’ll define and apply the Laplace Transform (LT) and its inverse (ILT) and use them to show the analytical kimesurface reconstruction.

4.4.2 Discrete LT

# create the LT
NuLT = function(datax, datay, inputz, k = 3, fitwarning = FALSE, 
                mirror = FALSE, range = 2*pi) {
  
  datax = as.numeric(datax)
  datay = as.numeric(datay)

  n = length(datax)
  x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*range
  
  if(mirror){
    x1 = c(x1,rev(2*range - x1))/2
    n = 2*n
    datay = c(datay, rev(datay))
    plot(x1, datay)
  }
  
  #generate the coefficients in indefinite integral of t^n*exp(-zt)
  coef = 1;
  coefm = as.matrix(coef)
  for(i in 1:k){
    coefm = cbind(coefm,0)
    coef = c(coef*i,1)
    coefm = rbind(coefm,coef)
  }
  # these coefficients ordered by ^0, ^1, ^2, ... in column format
  
  # compute 1, z, z^2...,z^k
  zz = cbind(1,inputz)
  zt = inputz
  for (i in 2:k){
    zt = zt*inputz
    zz = cbind(zz,zt)
  }
  zd = zt*inputz
  
  # compute 1, x, x^2...,x^k
  tx = x1;
  xm = cbind(1,x1)
  for (i in 2:k){
    tx = tx*x1
    xm = cbind(xm,tx)
  }
  
  
  # sum over intervals
  result = 0*inputz
  ii = 1
  while(ii+k <= n) {
    A = xm[seq(ii,ii+k),c(0:k+1)]
    b = datay[seq(ii,ii+k)]
    # polyfit might be faster when using polynomial basis, while matrix inverse, `solve()`,
    # is the more general approach that works for any function basis
    polyc = as.numeric(solve(A,b))

  
    #ordered by ^0, ^1, ^2, ... in column format
    
    # Enter a new function variable qualityCheck=FALSE
    # check fit quality; this step can be skipped for speed/efficiency
    # if (qualityCheck) { .... }
    
    if (fitwarning){
      xx = seq(A[1,2],A[k+1,2],length.out = 100);
      yy = polyval(rev(polyc),xx)
      if(max(abs(yy-mean(b)))>2*max(abs(b-mean(b)))){
        print(c("Warning: Poor Fit at ",ii,", Largest Deviation is",
                max(abs(yy-mean(b)))))
        print(c("Spline Polynomial is", polyc),3);
        #print(c(polyval(rev(polyc),A[,2]),b))
        plot(xx, yy, main="Polynomial fit", ylab="", type="l", col="blue")
        lines(A[,2],b, col="red")
        legend("topleft",c("fit","data"),fill=c("blue","red"))
        print(" ")
      }
    }
    
    # Use vector/matrix operations to avoid looping, 
    # some of the expressions look weird
    # May need to actually compare the efficiency/speed of
    # vector based vs. standard numeric calculations
    
    m1 = t(t(polyc*coefm)*A[1,])
    m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1] 
    
    m2 = t(t(polyc*coefm)*A[k+1,])
    m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
    
    intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
               exp(-inputz*A[k+1,2])*colSums(t(zz)*m22))/zd
    result = result+intgl
    ii=ii+k
  }
  
  # Computations over the last interval
  if(ii < n){
    nk = n-ii;
    A = xm[seq(ii,ii+nk),c(0:nk+1)]
    b = datay[seq(ii,ii+nk)]
    nc = as.numeric(solve(A,b))
    nc = c(nc,seq(0,0,length.out = k-nk))
    
    A = xm[seq(ii,ii+nk),]
    m1 = t(t(nc*coefm)*A[1,])
    m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1]
    
    m2 = t(t(nc*coefm)*A[nk+1,])
    m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
    
    # cc = colSums(coefm*polyc)
    intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
               exp(-inputz*A[nk+1,2])*colSums(t(zz)*m22))/zd
    result = result+intgl
  }
  #offset = 0.05*pi
  #result = result+datay[n]*(exp(-2*pi*inputz)-exp(-(2*pi+offset)*inputz))/inputz
  return(result)
}

Let’s test the discrete LT using the $\sin(x),\cos(x)$ function.

datax = seq(1,200)

n = length(datax)

x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)

datay = cos(x1)  # test function!!!

Lout = 61
range_limit <- 20
x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
y2 <- seq(from = 0, to = range_limit, length.out = Lout)

x2_g = x2 %o% seq(1,1, length.out = Lout)    
y2_g =  seq(1,1, length.out = Lout)%o%y2 
z2_grid = x2 %o% y2

argz = as.vector(x2_g + 1i*y2_g)
LTz = NuLT(x1, datay, argz)
rec1 = matrix(LTz,nrow = Lout)

recm = abs(rec1)
recr = Re(rec1)
reci = Im(rec1)
ph = reci/recm
surf_color <- atan2(reci,recr)
# colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
p <-  plot_ly(hoverinfo="none", showscale = TRUE) %>%
    add_trace(x = x2_g, y = y2_g, z = surf_color, 
    # F-magnitude or Re(F),   # z = Im(z2_grid),  # Real or Imaginary part of F(t)
          surfacecolor=log(recm), # colorscale=colorscale, #Phase-based color
          type = 'surface', opacity=1, visible=T) %>%
  layout(title = 
    "Truncated Laplace transform of cos(x), Color = Mag(Z)", 
         showlegend = TRUE,
         scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0), 
              xaxis=list(title= "Real=Re(z)"), yaxis=list(title = "Imaginary=Im(z)"),
              zaxis=list(title= "Z = tan-1(Im(L[f](z))/Re(L[f](z))))")) ) # 1:1:1 aspect ratio
p

Compare the above numerical calculation of the Laplace Transform (LT) of $f(x)=\cos(x)$ against the below kimesurface of the exact analytical (ground truth integration) of the LT of $f(x)=\cos(x)$, which is

\[F(z)\equiv \mathcal{L}(f)(z)\equiv \int_0^{2\pi}f(t)e^{-zt}dt = \frac{z-e^{-2\pi z}z}{1+z^2}\ .\]

z_val <- (x2_g+y2_g*1i)
ground_truth <- (z_val - exp(-2*pi*z_val)*z_val)/(1+z_val^2)
surf_ground_truth <- atan2(Im(ground_truth),Re(ground_truth))
plot_ly(hoverinfo="none", showscale = TRUE) %>%add_trace(x = x2_g, y = y2_g, z = surf_ground_truth,surfacecolor=log(abs(ground_truth)),type = 'surface', opacity=1, visible=T)%>%layout(title = 
    "Ground Truth for the truncated Laplace transform of cos(x), Color = Mag(Z)",scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0), 
              xaxis=list(title= "Real=Re(z)"), yaxis=list(title = "Imaginary=Im(z)"),
              zaxis=list(title= "Z = tan-1(Im(L[f](z))/Re(L[f](z)))")) )

The visual agreement between the exact and numerical calculations of the kimesurface $F(z)$ corresponding to the signal $f(x)=\cos(x)$ can also be checked by computing the average surface difference as shown below.

print(paste0("The average difference between the ground truth and the numerical integration is ",
             sum(abs(ground_truth -rec1))/length(ground_truth)))

## [1] "The average difference between the ground truth and the numerical integration is 0.000156495820765476"

# datax = seq(1,200)
# 
# n = length(datax)
# 
# x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)
# 
# datay = sin(x1)  # test function!!!
# 
# Lout = 61
# range_limit <- 20
# x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
# y2 <- seq(from = 0, to = range_limit, length.out = Lout)
# 
# x2_g = x2 %o% seq(1,1, length.out = Lout)    
# y2_g =  seq(1,1, length.out = Lout)%o%y2 
# z2_grid = x2 %o% y2
# 
# argz = as.vector(x2_g + 1i*y2_g)
# LTz = NuLT(x1, datay, argz)
# rec1 = matrix(LTz,nrow = Lout)
# 
# recm = abs(rec1)
# recr = Re(rec1)
# reci = Im(rec1)
# ph = reci/recm
# surf_color <- atan2(reci,recr)
# # colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
# p <-  plot_ly(hoverinfo="none", showscale = TRUE) %>%
#     add_trace(x = x2_g, y = y2_g, z = surf_color, 
#     # F-magnitude or Re(F),   # z = Im(z2_grid),  # Real or Imaginary part of F(t)
#           surfacecolor=log(recm), # colorscale=colorscale, #Phase-based color
#           type = 'surface', opacity=1, visible=T) %>%
#   layout(title = 
#     "Laplace transform of the sin(x), Color = log(mag)", 
#          showlegend = TRUE,
#          scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0), 
#               xaxis=list(title= "Real=Re(z)"), yaxis=list(title = "Imaginary=Im(z)"),
#               zaxis=list(title= "Z = tan-1(Im(L[f](z))/Re(L[f](z)))")) ) # 1:1:1 aspect ratio
# p

# z_val <- (x2_g+y2_g*1i)
# ground_truth <- (1 - exp(-2*pi*z_val))/(1+z_val^2)
# surf_ground_truth <- atan2(Im(ground_truth),Re(ground_truth))
# plot_ly(hoverinfo="none", showscale = TRUE) %>%add_trace(x = x2_g, y = y2_g, z = surf_ground_truth,surfacecolor=log(abs(ground_truth)),type = 'surface', opacity=1, visible=T)%>%layout(title = 
#     "Ground Truth for the truncated Laplace transform of sin(x), Color = Mag(Z)",scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0), 
#               xaxis=list(title= "Real=Re(z)"), yaxis=list(title = "Imaginary=Im(z)"),
#               zaxis=list(title= "Z = tan-1(Im(L[f](z))/Re(L[f](z)))")) )

Next, apply the discrete LT to the average-On (onAvg) and average-Off signals (offAvg), interpolating from their original size, $n=10$, to a new supersampled size $n=200$, and transforming the time support from $t\in[1:10]$ in increments of $\Delta t=1$, to $t' \in[0,2\pi)$, in increments of $\Delta t'=\frac{n}{n+0.5}\times \frac{1}{(200-1)2\pi}.$

This numerical longitudinal data preprocessing is done purely to establish some homologies in the structure of the LT domain, i.e., the input space signals (time-series), and the LT Range, i.e., the output space manifold (kimesurface). See the DSPA2 signal interpolation appendix to find out how to regularize either regularly (equally-spaced) sampled longitudinal data or irregularly (unequally-spaced) sampled longitudinal data.

datax = seq(1,200)
n = length(datax)
x_old <- c(1:10)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)
# longitudinal data series are: onAvg & offAvg
xnew <- x1
spl_onAvg <- spline(x=x_old, y=onAvg, xmin=min(x_old), xmax=max(x_old), n=200)
spl_offAvg <- spline(x=x_old, y=offAvg, xmin=min(x_old), xmax=max(x_old), n=200)

plot_ly(type="scatter", mode="markers") %>%
  add_trace(x=x_old, y=onAvg, name="Raw ON-Avg", 
            marker=list(size=20, color="lightblue", sizemode="area")) %>%
  add_trace(x=x_old, y=offAvg, name="Raw OFF-Avg",
            marker=list(size=20, color="orange", sizemode="area")) %>%
  add_trace(x=spl_onAvg$x, y=spl_onAvg$y, type="scatter", mode="markers+lines",
            name="Spline ON-Avg Model", marker=list(size=8, color="blue"),
            line = list(color = "blue", width = 4)) %>%
  add_trace(x=spl_offAvg$x, y=spl_offAvg$y, type="scatter", mode="markers+lines",
            name="Spline OFF-Avg Model", marker=list(size=8, color="red"),
            line = list(color = "red", width = 4)) %>%
  layout(title="Spline Modeling of 1D ON and OFF fMRI data\n (averaged over repeated samples)",
         legend = list(orientation = 'h', y=-0.2)) %>%
  hide_colorbar()

# datay = sin(x1)  # test function!!!

# Define a time-series to kimesurface plotting function
plotKimeSurface <- function(datay = sin(x1),index="", 
           title="Laplace transform of the Time-Series, Color = phase(Z)") {
  Lout = 61
  range_limit <- 20
  x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
  y2 <- seq(from = 0, to = range_limit, length.out = Lout)
  
  x2_g = x2 %o% seq(1,1, length.out = Lout)    
  y2_g =  seq(1,1, length.out = Lout)%o%y2 
  # z2_grid = x2 %o% y2
  
  argz = as.vector(x2_g + 1i*y2_g)
  # print(argz)
  LTz = NuLT(x1, datay, argz)
  rec1 = matrix(LTz,nrow = Lout)
  
  recm = abs(rec1)
  recr = Re(rec1)
  reci = Im(rec1)
  ph = reci/recm
  surf_color <- atan2(reci,recr)
  # colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
  p <-  plot_ly(hoverinfo="none", showscale = TRUE,
           scene = paste("scene",as.character(index),sep = "")) %>%
      add_trace(x = x2_g, y = y2_g, z = surf_color, 
      # F-magnitude or Re(F),   # z=Im(z2_grid),  # Real or Imaginary part of F(t)
            surfacecolor=log(recm), # colorscale=colorscale, #Phase-based color
            type = 'surface', opacity=1, visible=T) %>%
      layout(title = title, showlegend = TRUE,scene=list(aspectmode="manual", aspectratio=list(x=1, y=1, z=1.0), 
                xaxis=list(title= "Real"), yaxis=list(title = "Imaginary"),
                zaxis=list(title= "Z = phase")) ) # 1:1:1 aspect ratio)
  
  return(list("plot"=p,"phase"=surf_color,"magnitude"=recm,"x_grid"=x2_g,"y_grid"=y2_g))
}
#paste("scene",as.character(index),sep = "")
# list(aspectmode="manual", aspectratio=list(x=1, y=1, z=1.0), 
#                 xaxis=list(title= "Real"), yaxis=list(title = "Imaginary"),
#                 zaxis=list(title= "Z = log(mag(Lf(x+iy)))")) ) # 1:1:1 aspect ratio
pOn  <- plotKimeSurface(datay=spl_onAvg$y, "",
           title="ON fMRI ((44,42,33): Laplace transform of Avg Time-Series, Color = phase(Z)")
pOff <- plotKimeSurface(datay=spl_offAvg$y,2, 
           title="OFF fMRI ((44,42,33): Laplace transform of Avg Time-Series, Color = phase(Z)")

# combineWidgets(pOn, pOff)

# synch both 3D scenes in 1 3D plot
# install.packages("manipulateWidget")
# subplot(pOn, pOff, nrows = 1, shareX = TRUE)

# main_plot <- subplot(pOn, pOff, nrows = 2, shareX = TRUE, margin = 0.06) %>% 
#   layout(title = "Kimesurfaces On (activation) and Off (rest) fMRI voxel (44,42,33)", 
#          scene  = list(domain = list(x = c(0, 0.5), y = c(0.5, 1)), aspectmode = "cube"), 
#          scene2 = list(domain = list(x = c(0.2, 0.7), y = c(0.5, 1)), aspectmode = "cube"))
# 
# main_plot %>% 
#   htmlwidgets::onRender(
#     "function(x, el) {
#       x.on('plotly_relayout', function(d) {
#         const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
#         if (camera.length) {
#           const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
#           const new_layout = {};
#           scenes.forEach(key => {
#             new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
#           });
#           Plotly.relayout(x, new_layout);
#         }
#       });
#     }")

library(plotly)
library(htmlwidgets)

# Generate some example 3D data
# ... (Create your pOn and pOff plots here)
diff_plot <- plot_ly(hoverinfo="none", showscale = TRUE,
           scene = "scene3") %>%
      add_trace(x = pOn$x_grid, y = pOn$y_grid, z = pOn$phase - pOff$phase, 
      # F-magnitude or Re(F),   # z=Im(z2_grid),  # Real or Imaginary part of F(t)
            surfacecolor=pOn$magnitude-pOff$magnitude, # colorscale=colorscale, #Phase-based color
            type = 'surface', opacity=1, visible=T) %>%
      layout(title = title, showlegend = TRUE,scene3=list(aspectmode="manual", aspectratio=list(x=1, y=1, z=1.0), 
                xaxis=list(title= "Real"), yaxis=list(title = "Imaginary"),
                zaxis=list(title= "Z = Phase")) )
# Combine them into a subplot
main_plot <- subplot(pOn$plot, pOff$plot,diff_plot) %>%
  layout(
    title = "kimesurfaces (Left: On; Right: Off; Bottom: Kimesurface Difference (On vs. Off))",
    scene  = list(domain = list(x = c(0, 0.5), y = c(0.5, 1)), aspectmode = "cube"),
    scene2 = list(domain = list(x = c(0.5, 1), y = c(0.5, 1)), aspectmode = "cube"),
    scene3 = list(domain = list(x = c(0.25, 0.75), y = c(0, 0.5)), aspectmode = "cube")
  )
#https://stackoverflow.com/questions/77757102/possible-to-synchronize-3d-scatter-plotly-plots-in-r
main_plot %>%
  htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }")

4.5 Time-series to Kimesurface Mapping using Integral Transforms

Next, we will explore alternative strategies to map repeated measurement time-series into kimesurfaces. These approaches attempt to resolve the issues with mapping the highly irregular (non-analytic/non-holomorphic) kimesurfaces to more locally smooth analytic functions.

4.5.1 Integral Transformation of Stiching-based Kimesurface Recontrucitons

Protocol:

Apply the standard stitching strategy to transform the repeated time-series into a kimesurface using a random sample of kime-phases, $\varphi\sim\Phi_{[-\pi,\pi)}$.
Use the inverse of an integral transform, such as Hilbert, Fourier, and …, to map the irregular ($\mathbb{C}$-valued) kimesurface into a spacetime-like locally smooth (analytic?) manifold.
Explore inference based on these spacetime-like kimesurface reconstructions and contrast these against other statistical/data-analytic/AI tools for prediction, classification based on the original time-series.

Integral transforms are mathematical operators maping functions from one domain to another, often simplifying the analysis of problems by transforming differential equations into algebraic equations or by converting convolutions into products. Each transform has specific properties, applications, and interpretations, making them essential tools in areas like signal processing, physics, and engineering. Below are commontly used integral transforms, along with their interrelations and interpretations of their forward and reverse transforms.

4.5.1.1 Fourier Transform

The Fourier Transform, one of the most widely used integral transforms, decomposes a function into its constituent frequencies, transforming a time-domain signal into the frequency domain (k-space). For a continuous function $f(t)$, the Fourier transform $\mathcal{F}\{f(t)\} = F(\omega)$ is defined by $F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} \, dt .$

The inverse Fourier transform recovers $f(t)$ from $F(\omega)$ using $f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} \, d\omega .$

The forward Fourier Transform (FT) maps a time-domain signal into the frequency domain, where $F(\omega)$ represents the amplitude and phase of the frequencies present in the original signal. Whereas the inverse FT reconstructs the time-domain representation of the signal by summing over all frequency components weighted by their amplitudes and phases.

4.5.1.2 Laplace Transform

The Laplace Transform (LT) is closely related to the Fourier Transform but is typically used for functions defined on $[0, \infty)$, i.e., time-series. It is particularly useful in solving linear differential equations and control theory problems. For a function of time, $f(t)$, the Laplace transform $\mathcal{L}\{f(t)\} = F(s)$ is defined by $F(s) = \int_0^{\infty} f(t) e^{-st} \, dt .$

Then, the inverse Laplace transform, used to recover $f(t)$ from $F(s)$, is $f(t) = \frac{1}{2\pi i} \int_{\gamma - i\infty}^{\gamma + i\infty} F(s) e^{st} \, ds .$ The forward LT sends a time-domain function into the complex $\mathbb{C}$-domain, capturing both growth and oscillatory behavior. The inverse LT similarly recovers the time-domain function from its $z$-domain representation, often involving complex analysis techniques like contour integration.

4.5.1.3 Hilbert Transform

The Hilbert Transform (HT) is an integral transform that provides a phase shift to each frequency component of a signal. It is particularly useful in signal processing for creating analytic signals and in the study of harmonic analysis. For a function $f(t)$, the Hilbert transform $\mathcal{H}\{f(t)\} = \tilde{f}(t)$ is defined by $\tilde{f}(t) = \frac{1}{\pi} \text{P.V.} \int_{-\infty}^{\infty} \frac{f(\tau)}{t - \tau} \, d\tau ,$ where P.V. denotes the Cauchy principal value of the integral.

The HT produces a signal where each frequency component of the original signal is phase-shifted by $90^\circ$. The combination of the original signal and its Hilbert transform creates an analytic signal (complex signal with only positive frequency components). Applying the Hilbert transform twice to a function $u()$ yields $\operatorname {H} {\bigl (}\operatorname {H} (u){\bigr )}(t)=-u(t)$, assuming the integrals defining both iterations converge in a suitable sense. Hence, the inverse HT is $H^{(-1)}=-\operatorname {H}$. This links the Hilbert transform of $u(t)$ to the Fourier transform of $u(t)$. The Hilbert transform is an operator \[\overbrace{\mathcal {F}}^{FT}{\bigl (}\overbrace{\operatorname {H}}^{HT} (u){\bigr )}(\omega )=-i\operatorname {sgn}(\omega )\cdot {\overbrace{\mathcal {F}}^{FT}}(u)(\omega ).\]

As $\operatorname {sgn} = \operatorname {sgn}(2\pi x)$, this relations applies to the three common FT definitions. The HT is typically applied to real-valued signals to extract the imaginary part of an analytic signal.

4.5.1.4 Z-Transform

The Z-transform is the discrete counterpart of the Laplace Transform, typically used for discrete-time signals. It is extensively used in digital signal processing and control theory. For a discrete sequence $f[n]$, the forward Z-transform $\mathcal{Z}\{f[n]\} = F(z)$ is defined by $F(z) = \sum_{n=-\infty}^{\infty} f[n] z^{-n} .$ The inverse Z-transform is $f[n] = \frac{1}{2\pi i} \oint F(z) z^{n-1} \, dz .$

Forward ZT converts a discrete-time sequence into a complex frequency domain where $z\in\athbb{C}$, where as the inverse ZT recovers the original discrete-time sequence from its $z$-domain representation, typically involving residue calculus or partial fraction expansion.

4.5.1.5 Mellin Transform

The Mellin Transform (MT) is another integral transform often used in problems involving scaling, which connects multiplicative convolutions with additive convolutions. For a function $f(x)$, the Mellin transform $\mathcal{M}\{f(x)\} = F(s)$ is $F(s) = \int_0^{\infty} x^{s-1} f(x) \, dx .$ The inverse MT is $f(x) = \frac{1}{2\pi i} \int_{c - i\infty}^{c + i\infty} F(s) x^{-s} \, ds .$ The forward MT converts a function into the Mellin domain, which is useful when analyzing problems that exhibit scaling properties (e.g., in image processing or fractal analysis). Whereas the inverse MT reconstructs the function from its Mellin transform, often using contour integration techniques.

The Fourier transform is a special case of the Laplace transform where the real part of $z$ (in the Laplace domain) is set to zero, i.e., $z = i\omega$. Thus, the Laplace Transform generalizes the Fourier Transform by considering exponential growth/decay in addition to oscillations. Similarly, the Hilbert transform is closely related to the Fourier transform, as it can be viewed as a frequency-domain filter that shifts the phase of the positive frequencies by $-90^\circ$ and the negative frequencies by $+90^\circ$. This makes the Hilbert transform useful for analytic signal generation. The Z-transform can be viewed as a discrete-time version of the Laplace transform, just as the Discrete-Time Fourier transform (DTFT) is the discrete counterpart of the continuous Fourier Transform. When $z = e^{i\omega}$, the Z-transform becomes the DTFT.
Finally, the Mellin transform is related to the Fourier Transform by a logarithmic scaling of the input variable. If a function’s domain is logarithmically transformed, the Mellin transform becomes a Fourier transform in the logarithmic domain.

4.5.2 Fourier Transform Reconstructions

Protocol:

Render the basic stitching kimesurfaces reconstructed in the native space. these will provide the baseline for contrasting the subsequent kime-space kimesurface reconstructions using stitching in the Fourier domain.
Apply the FT to each time-series instance of the repeated fMRI (on and Off) measurement time-series.
K-space stitching together (using random kime-phase draws, $\varphi\sim\Phi_{[-\pi,\pi)}$) all these complex-valued FT series into a kimesurface in the Fourier domain.
IFT the k-space constructed kimesurfaces into the native space and characterize it’s local smoothness, relative to the smoothness of the kimesurface reconstructions obtained via basic stitching in the native space.

4.5.2.1 Kimesurface Reconstructions via Basic Stitching in the Native-space

library(spatstat)
library(plotly)
library(htmlwidgets)

# Extract just the On kime-series at voxel (44,42,33), each time-series includes 10 time-points and has 8 repeats!
# onFMRISeries <- bigim1_mod[44,42,33, c(1:10, 21:30, 41:50, 61:70, 
#                                        81:90, 101:110, 121:130, 141:150)]
# # the corresponding Off kime-series at voxel (44,42,33) will be the temporal complement
# offFMRISeries <- bigim1_mod[44,42,33, -c(1:10, 21:30, 41:50, 61:70, 
#                                          81:90, 101:110, 121:130, 141:150)]

# t_indx <- seq(1, 80, 1) # index of time (not actual time value)
# f_On <- onFMRISeries
# f_Off <- offFMRISeries

# Compute and plot against each other the Average On and Average Off time-series
period <- 10   # repeated measurement period
numSeries <- 8 # number of fMRI timeseries in the entire run (for each of On and Off signals)
extract8fMRI_TimeSeries <- function(f=0, period=0) {
  # Split the fMRI vector into (numSeries=8) individual time-courses
  fMRI_splitSeries <- split(f, ceiling(seq_along(f) / period))
  return(fMRI_splitSeries)
}

# 1. Extract numSeries=8 series for stitching and kimesurface reconstruction
fMRI_ON_Series <- extract8fMRI_TimeSeries(f=f_On, period=10)
fMRI_OFF_Series <- extract8fMRI_TimeSeries(f=f_Off, period=10)

# 1.a. Stitch together (using random kime-phase draws, $\varphi\sim\Phi_{[-\pi,\pi)}$) all these complex-valued FT series into a kimesurface in k-space.
phi_8_vec <- matrix(NA, ncol=period, nrow = numSeries)
for (t in 1:period) { 
  # for a given t, generate 8 new phases
  set.seed(t);
  phi_8_vec[ ,t] <- extraDistr::rlaplace(numSeries, mu=0, sigma=0.5)
  # rank-order the phases for consistency
  # within the same foliation leaf
  phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t])
  # force phases in [-pi: pi)
  for (i in 1:numSeries) {
    if (phi_8_vec[i,t] < -pi) phi_8_vec[i,t] <- -pi
    if (phi_8_vec[i,t] >= pi) phi_8_vec[i,t] <- pi
  }
}

matrix_ON <- matrix(0, nrow = 2*period+1, ncol = 2*period+1) 
matrix_OFF <- matrix(0, nrow = 2*period+1, ncol = 2*period+1) 
for (t in 1:period) {  # POLAR to CARTESIAN coordinates
  for (p in 1:numSeries) {
    x = (period+1) + t*cos(phi_8_vec[p,t])
    y = (period+1) + t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- fMRI_ON_Series[[p]][t]
    matrix_OFF[x,y] <- fMRI_OFF_Series[[p]][t]
  }
}
# smooth/blur the matrices
# matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
# matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))
matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))

# Display the Re() and Im() parts of the complex k-space kimesurfaces
# Load necessary libraries
# Create sequences for x and y axes
x <- seq(1, 2*period+1, length.out = 2*period+1)
y <- seq(1, 2*period+1, length.out = 2*period+1)

# Create the first surface plot
pON_Re <- plot_ly(x=x, y=y, z=Re(matrix_ON), type="surface", 
                  scene='scene1', name="Re(ON)") %>% 
  layout(title = "Re(fMRI-ON Basic Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Re <- plot_ly(x=x, y=y, z=Re(matrix_OFF), type="surface", 
                   scene='scene2', name="Re(OFF)") %>% 
  layout(title = "Re(fMRI-OFF Basic Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pDIFF_Re <- plot_ly(x=x, y=y, z=Re(matrix_ON-matrix_OFF), type="surface",
                    scene='scene3', name="Re(Diff)") %>% 
  layout(title = "Re(fMRI-ON - fMRI-OFF Basic Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
# Scene titles
annotations=list(list(showarrow=FALSE, text='(ON)', 
                      xref='scene',yref='scene',zref='scene', x=0.0,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='(OFF)', 
                      xref='scene2',yref='scene',zref='scene', x=0.35,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='(DIFF)', 
                      xref='scene3',yref='scene',zref='scene', x=0.7,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)))
# Combine the plots
combined_plot <- subplot(pON_Re, pOFF_Re, pDIFF_Re) %>% #, nrows=2,margin=0.05) %>%
  layout(title="Native-space Reconstruction of fMRI ON and OFF Kimesurfaces using Basic Stitching", annotations=annotations)

# Display the synchronized plot
combined_plot %>% htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }")

4.5.2.2 Kimesurface Reconstructions via Fourier-space Stitching

Next, we will examine the kimesurface reconstructions via k-space stitching together. This process uses random kime-phase draws, $\varphi\sim\Phi_{[-\pi,\pi)}$, e.g., Laplace phase distribution) coupled with the complex-valued FT series of the original repeated measurement fMRI time-series into a kimesurface in the Fourier domain. We will contrast the IFT of the k-space constructed kimesurfaces back into the native space and characterize it’s local smoothness, relative to the smoothness of the kimesurface reconstructions obtained via basic stitching in the native space.

# 2. Apply the FT to each instance of the repeated measurement time-series.
library(spatstat)

FT_fMRI_ON_Series <- list()
FT_fMRI_OFF_Series <- list()
for (i in c(1:numSeries)) {
  FT_fMRI_ON_Series[[i]] = fft(fMRI_ON_Series[[i]])
  FT_fMRI_OFF_Series[[i]] = fft(fMRI_OFF_Series[[i]])
}

# 3. Stitch together (using random kime-phase draws, $\varphi\sim\Phi_{[-\pi,\pi)}$) all these complex-valued FT series into a kimesurface in k-space.
phi_8_vec <- matrix(NA, ncol=period, nrow = numSeries)
for (t in 1:period) { 
  # for a given t, generate 8 new phases
  set.seed(t);
  phi_8_vec[ ,t] <- extraDistr::rlaplace(numSeries, mu=0, sigma=0.5)
  # rank-order the phases for consistency
  # within the same foliation leaf
  phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t])
  # force phases in [-pi: pi)
  for (i in 1:numSeries) {
    if (phi_8_vec[i,t] < -pi) 
      phi_8_vec[i,t] <- -pi
    if (phi_8_vec[i,t] >= pi) 
      phi_8_vec[i,t] <- pi
  }
}

matrix_ON <- matrix(0, nrow = 2*period+1, ncol = 2*period+1) 
matrix_OFF <- matrix(0, nrow = 2*period+1, ncol = 2*period+1) 
for (t in 1:period) {  # POLAR to CARTESIAN coordinates
  for (p in 1:numSeries) {
    x = (period+1) + t*cos(phi_8_vec[p,t])
    y = (period+1) + t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- FT_fMRI_ON_Series[[p]][t]
    matrix_OFF[x,y] <- FT_fMRI_OFF_Series[[p]][t]
  }
}
# smooth/blur the matrices
# matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
# matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))
matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))

# Display the Re() and Im() parts of the complex k-space kimesurfaces
# Load necessary libraries
# library(plotly)
# library(htmlwidgets)

# Create sequences for x and y axes
x <- seq(1, 2*period+1, length.out = 2*period+1)
y <- seq(1, 2*period+1, length.out = 2*period+1)

# Create the first surface plot
pON_Re <- plot_ly(x=x, y=y, z=Re(matrix_ON), type="surface", 
                  scene='scene1', name="Re(ON)") %>% 
  layout(title = "Re(fMRI-ON Kime Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Re <- plot_ly(x=x, y=y, z=Re(matrix_OFF), type="surface", 
                   scene='scene2', name="Re(OFF)") %>% 
  layout(title = "Re(fMRI-OFF Kime Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pDIFF_Re <- plot_ly(x=x, y=y, z=Re(matrix_ON-matrix_OFF), type="surface",
                    scene='scene3', name="Re(Diff)") %>% 
  layout(title = "Re(fMRI-ON - fMRI-OFF Kime Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pON_Im <- plot_ly(x=x, y=y, z=Im(matrix_ON), type = "surface", 
                  scene='scene4', name="Im(On)") %>% 
  layout(title = "Im(fMRI-ON Kime Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Im <- plot_ly(x=x, y=y, z=Im(matrix_OFF), type = "surface", 
                  scene='scene5', name="Im(Off)") %>% 
  layout(title = "Im(fMRI-OFF Kime Stiching)") %>% hide_colorbar()
pDIFF_Im <- plot_ly(x=x, y=y, z=Im(matrix_ON-matrix_OFF), type = "surface", 
                  scene='scene6', name="Im(Diff)") %>% 
  layout(title = "Im(fMRI-ON - fMRI-OFF Kime Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()

# Scene titles
annotations=list(list(showarrow=FALSE, text='Re(ON)', 
                      xref='scene',yref='scene',zref='scene', x=0.0,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Re(OFF)', 
                      xref='scene2',yref='scene',zref='scene', x=0.35,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Re(DIFF)', 
                      xref='scene3',yref='scene',zref='scene', x=0.7,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(ON)', 
                      xref='scene4',yref='scene',zref='scene', x=0.0,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(OFF)', 
                      xref='scene5',yref='scene',zref='scene', x=0.35,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(DIFF)', 
                      xref='scene6',yref='scene',zref='scene', x=0.7,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)))

# Combine the plots
combined_plot <- subplot(pON_Re, pOFF_Re, pDIFF_Re, 
                         pON_Im, pOFF_Im, pDIFF_Im) %>% #, nrows=2,margin=0.05) %>%
  layout(title="Real (Re) and Imaginary (Im) fMRI-ON vs. fMRI-OFF Kime Stitching Surfaces", 
         scene=list(domain=list(x=c(0, 0.33), y=c(0.5,1)), aspectmode="cube"), 
         scene2=list(domain=list(x=c(0.33, 0.66), y=c(0.5,1)), aspectmode="cube"),
         scene3=list(domain=list(x=c(0.7, 1), y = c(0.5,1)), aspectmode="cube"),
         scene4=list(domain=list(x=c(0, 0.33), y=c(0.0,0.5)), aspectmode="cube"), 
         scene5=list(domain=list(x=c(0.33,0.66), y=c(0.0,0.5)), aspectmode="cube"),
         scene6=list(domain=list(x=c(0.7, 1), y = c(0.0,0.5)), aspectmode="cube"),
         annotations = annotations)

# Display the synchronized plot
combined_plot %>% htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }")

# 4.  IFT the k-space constructed kimesurface into spacetime and 
# characterize it's local smoothness, relative to the smoothness 
# of the original k-space kimesurface.
# matrix_ON, matrix_OFF, and (matrix_ON - matrix_OFF)
IFT_matrix_ON <- fft(matrix_ON, inverse = T)
IFT_matrix_OFF <- fft(matrix_OFF, inverse = T)

# smooth/blur the matrices
# matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
# matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))
Re_IFT_matrix_ON_smooth <- 
  (1/100)*as.matrix(blur(as.im(Re(IFT_matrix_ON)), sigma=2.0))
Re_IFT_matrix_OFF_smooth <- 
  (1/100)*as.matrix(blur(as.im(Re(IFT_matrix_OFF)), sigma=2.0))
Im_IFT_matrix_ON_smooth <- 
  (1/100)*as.matrix(blur(as.im(Im(IFT_matrix_ON)), sigma=2.0))
Im_IFT_matrix_OFF_smooth <- 
  (1/100)*as.matrix(blur(as.im(Im(IFT_matrix_OFF)), sigma=2.0))
Mag_IFT_matrix_ON_smooth <- 
  (1/100)*as.matrix(blur(as.im(Mod(IFT_matrix_ON)), sigma=2.0))
Mag_IFT_matrix_OFF_smooth <- 
  (1/100)*as.matrix(blur(as.im(Mod(IFT_matrix_OFF)), sigma=2.0))

# # Create sequences for x and y axes
# x <- seq(1, 2*period+1, length.out = 2*period+1)
# y <- seq(1, 2*period+1, length.out = 2*period+1)

# Create the first surface plot
pON_Re <- plot_ly(x=x, y=y, z=Re_IFT_matrix_ON_smooth, type="surface", 
                  scene='scene1', name="Re(ON)") %>% 
  layout(title = "Re(fMRI-ON K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Re <- plot_ly(x=x, y=y, z=Re_IFT_matrix_OFF_smooth, type="surface", 
                   scene='scene2', name="Re(OFF)") %>% 
  layout(title = "Re(fMRI-OFF K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pDIFF_Re <- plot_ly(x=x, y=y, 
                z=(Re_IFT_matrix_ON_smooth - Re_IFT_matrix_OFF_smooth),
                type="surface", scene='scene3', name="Re(Diff)") %>% 
  layout(title = "Re(fMRI-ON - fMRI-OFF K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pON_Im <- plot_ly(x=x, y=y, z=Im_IFT_matrix_ON_smooth, type = "surface", 
                  scene='scene4', name="Im(On)") %>% 
  layout(title = "Im(fMRI-ON K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Im <- plot_ly(x=x, y=y, z=Im_IFT_matrix_OFF_smooth, type = "surface", 
                  scene='scene5', name="Im(Off)") %>% 
  layout(title = "Im(fMRI-OFF K-space Stiching)") %>% hide_colorbar()
pDIFF_Im <- plot_ly(x=x,y=y,z=(Im_IFT_matrix_ON_smooth - Im_IFT_matrix_OFF_smooth),
                  type="surface", scene='scene6', name="Im(Diff)") %>% 
  layout(title = "Im(fMRI-ON - fMRI-OFF K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pON_Mag <- plot_ly(x=x, y=y, z=Mag_IFT_matrix_ON_smooth, type = "surface", 
                  scene='scene7', name="Mag(On)") %>% 
  layout(title = "Mag(fMRI-ON K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
pOFF_Mag <- plot_ly(x=x, y=y, z=Mag_IFT_matrix_OFF_smooth, type = "surface", 
                  scene='scene8', name="Mag(Off)") %>% 
  layout(title = "Mag(fMRI-OFF K-space Stiching)") %>% hide_colorbar()
pDIFF_Mag <- plot_ly(x=x,y=y,
                     z=(Mag_IFT_matrix_ON_smooth - Mag_IFT_matrix_OFF_smooth),
                  type="surface", scene='scene9', name="Im(Diff)") %>% 
  layout(title = "Mag(fMRI-ON - fMRI-OFF K-space Stiching)", 
         scene = list(domain=list(x=c(-period, period), y=c(-period, period)),
                      aspectmode = "cube")) %>% hide_colorbar()
# Scene titles
annotations=list(list(showarrow=FALSE, text='Re(ON)', 
                      xref='scene',yref='scene',zref='scene', x=0.0,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Re(OFF)', 
                      xref='scene2',yref='scene',zref='scene', x=0.35,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Re(DIFF)', 
                      xref='scene3',yref='scene',zref='scene', x=0.7,y=0.9,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(ON)', 
                      xref='scene4',yref='scene',zref='scene', x=0.0,y=0.5,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(OFF)', 
                      xref='scene5',yref='scene',zref='scene', x=0.3,y=0.5,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Im(DIFF)', 
                      xref='scene6',yref='scene',zref='scene', x=0.6,y=0.5,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Mag(ON)', 
                      xref='scene7',yref='scene',zref='scene', x=0.0,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Mag(OFF)', 
                      xref='scene8',yref='scene',zref='scene', x=0.35,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)),
                 list(showarrow=FALSE, text='Mag(DIFF)', 
                      xref='scene9',yref='scene',zref='scene', x=0.7,y=0,z=0.1,
                      xanchor='left',yanchor='bottom',font=list(size=16)))

# Combine the plots
combined_plot <- subplot(pON_Re, pOFF_Re, pDIFF_Re, 
                         pON_Im, pOFF_Im, pDIFF_Im, 
                         pON_Mag, pOFF_Mag, pDIFF_Mag) %>%
  layout(title="Real (Re), Imaginary (Im) and Magnitude of\n the fMRI-ON vs. fMRI-OFF Fourier-space Stitching Surfaces", 
         scene=list(domain=list(x=c(0, 0.33), y=c(0.67,1)), aspectmode="cube"), 
         scene2=list(domain=list(x=c(0.33,0.66),y=c(0.67,1)), aspectmode="cube"),
         scene3=list(domain=list(x=c(0.7, 1), y=c(0.67,1)), aspectmode="cube"),
         scene4=list(domain=list(x=c(0, 0.33), y=c(0.33,0.66)), aspectmode="cube"), 
         scene5=list(domain=list(x=c(0.33,0.66),y=c(0.33,0.66)), aspectmode="cube"),
         scene6=list(domain=list(x=c(0.7, 1), y = c(0.33,0.66)), aspectmode="cube"),
         scene7=list(domain=list(x=c(0, 0.33), y=c(0.0,0.33)), aspectmode="cube"), 
         scene8=list(domain=list(x=c(0.33,0.66), y=c(0.0,0.33)), aspectmode="cube"),
         scene9=list(domain=list(x=c(0.7, 1), y = c(0.0,0.33)), aspectmode="cube"),
         annotations = annotations)

# Display the synchronized plot
combined_plot %>% htmlwidgets::onRender(
    "function(x, el) {
      x.on('plotly_relayout', function(d) {
        const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
        if (camera.length) {
          const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
          const new_layout = {};
          scenes.forEach(key => {
            new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
          });
          Plotly.relayout(x, new_layout);
        }
      });
    }")

We can also consider using other integral transforms to obtain alternative kimesurface reconstructions in the native space.

5 Interactive plotly Example

For most of these examples, users may need to review/explore/extend the core functionality in the TCIU-package source-code (CRAN) to customize the TCIU functionality (GitHub) for the application specific needs!

5.1 Plotly method: interactive way

The function fmri_image() is used to create images for front view, side view, and top view of the fMRI image.

bigim1_mask<-bigim1_mod
for (i in 1:160) {
  bigim1_mask[,,,i]<-bigim1_mask[,,,i]*mask
}

5.1.1 Manual examples

A plot without (spatio-temporal) mask restrictions.

fmri_image(bigim1_mod, option = "manually", voxel_location = c(40,40,30), time = 4)

A plot with (spatio-temporal) mask restrictions.

fmri_image(bigim1_mask, option = "manually", voxel_location = c(40,40,30), time = 4)

5.2 Forecasting with time series

The function fmri_ts_forecast() fits the ARIMA models for each voxel (spatial volume element) location. This function is based on the TSplot_gen() function from package TSplotly.

fmri_ts_forecast(smoothmod,c(41,44,33))

6 Motor area detection

If there are concrete spatial locations (regions, voxels, coordinates) that are of specific interest, then one needs to focus on these locations. We’ll demonstrate this using fMRI brain activation (simulated) in the somatosensory (motor) cortex.

6.1 fMRI data simulation

The function fmri_simulate_func() is used to simulate real-valued fMRI data with specified dimensions (locations and extents).

fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
                                   ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                                   dur = c(10, 10, 10, 10, 10, 10, 10, 10))

6.2 Stimulus detection

The integrated function fmri_stimulus_detect() is designed to apply multiple analytical methods for activation detection, in this case in the (human brain) motor area.

Examples of parametric and non-parametric statistical tests already built in include:

“t-test” and “Wilcoxon-test”, which can be applied to all real*-valued fMRI data,
“Hotelling’s T2”, “Wilk’s Lambda” and “gLRT” methods, which can be applied to all complex-valued fMRI data,
“on_off_diff” and “HRF” (accounting for the hemodynamic response function) methods, which can be applied to $4D$ real-valued fMRI data where the first method is calculating on-off difference period polar volume to get p-values and the second method is using hemodynamic response function and linear regression,
The post-hoc stat-mapping filtering can be also applied to the calculated p-values if specified in the parameter or use the function fmri_post_hoc().

6.2.1 Examples

# statistical method HRF needs parameter ons and dur
pval1 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "HRF" , 
                             ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                             dur = c(10, 10, 10, 10, 10, 10, 10, 10) )

# statistical method t-test for real-valued fMRI data
pval2 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "t-test")

# statistical method Wilk's Lambda for complex-valued data
pval3 = fmri_stimulus_detect(fmridata = bigim1, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10)) ], 
                             method = "Wilks-Lambda" )

# do the fdr correction and the spatial clustering
# pval4 is the pval1 after the post-hoc processing
pval4 = fmri_post_hoc(pval1, fdr_corr = "fdr",
                                             spatial_cluster.thr = 0.05,
                                             spatial_cluster.size = 5, 
                                             show_comparison = FALSE)

Summarize the corresponding p-values.

summary(pval1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.8987  1.0000  0.8463  1.0000  1.0000

summary(pval2)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000096 0.8598920 1.0000000 0.8579245 1.0000000 1.0000000

summary(pval3)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000025 0.1877401 1.0000000 0.7476717 1.0000000 1.0000000

summary(pval4)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  1.0000  1.0000  0.9971  1.0000  1.0000

7 Motor area visualization

7.1 Visualization and comparison of p-value

7.1.1 3D visualization for p-value

pval1_3d = fmri_3dvisual(pval1, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="Accounting for HRF")

pval1_3d$plot

pval4_3D = fmri_3dvisual(pval4, mask, p_threshold = 0.05, method="scale_p",
                         multi_pranges=TRUE, title="HRF method with post-hoc correction")

pval4_3D$plot

7.1.2 2D visualization for p-value

Generate 2D plots of the 2D p-value (projection images) in sagittal, coronal and axial views.

Plot without isocontours.

for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=NULL, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}

Plot with isocontours.

for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=hemody, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}

7.2 Comparison of performance of different methods on the same fMRI data

7.2.1 3D p-value comparison

The function fmri_pval_comparison_3d() is used to visualize two p-value maps showing their differences in detecting stimulated brain areas in 3D scenes. Since our original fMRI is too big to use here for different statistical tests, in this example, we compare the differences of stimulated parts of two different fMRI data using the same mask.

fmri_pval_comparison_3d(list(pval1, pval2), mask, 
                                list(0.05, 0.05), list("scale_p", "scale_p"), 
                                multi_pranges=FALSE)

7.2.2 2D p-value comparison

The function fmri_pval_comparison_2d() displays the p-values (generated by different statistical tests on the same fMRI data) exposing their difference via 2D plots. For simplicity here we only compare two different 3D arrays of p-values.

fmri_pval_comparison_2d(list(pval2, pval1), 
                               list('t-test', 'HRF'),
                               list(list(35, 33, 22), list(40, 26, 33)), 
                               hemody_data = NULL, 
                               mask = mask, p_threshold = 0.05, 
                               legend_show = FALSE, method = 'scale_p',
                               color_pal = "YlOrRd", multi_pranges=FALSE)

8 Tri-phase ROI-based Spacekime Analytics

Below we demonstrate a three-tier statistical analysis (tri-phase) of regional activation.

Notes:

For large volumes, these calculations may be expensive (i.e., stats mapping may take a significant time)!
The following examples are specific to $4D$ fMRI spacekime analytics. In most situations, these processing steps need to be adapted to the concrete data format and the specific application domain. In practice, some of the TCIU package functions may be used as templates to extend this spacekime functionality to alternative data formats and different data structures.

8.1 Phase 1: Detect Potential Activated ROI

First, identify large local (regional areas) associated with activations/task stimuli.

phase1_pval = fmri_ROI_phase1(bigim1_mod, mask_label, mask_dict,
                              stimulus_idx = c(1:160)[rep(c(TRUE,FALSE),
                                                          c(10,10))])$all_ROI$pval_t

8.2 Phase 2: ROI-Based Tensor-on-Tensor Regression

phase2_pval = fmri_ROI_phase2(fmridata = bigim1_mod, 
                              label_mask = mask_label, label_dict = mask_dict, 
                              stimulus_idx = c(1, 21, 41, 61, 81, 101, 121, 141),
                              stimulus_dur = c(10, 10, 10, 10, 10, 10, 10, 10),
                              rrr_rank = 3, fmri.design_order = 2,
                              fmri.stimulus_TR = 3, method = "t_test")

8.3 Phase 3: FDR Correction and Spatial Clustering

phase3_pval = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
                            spatial_cluster.thr = 0.05,
                            spatial_cluster.size = 5, 
                            show_comparison = FALSE)

8.4 3D visualization based on the activated areas by regions

fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, title = "Phase 1 p-values", rank = "value")

fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, 5, title = "pPhase-1 top five p-values", rank = "value")

# for 3D visualization, user needs to include empty region in the label
label_index = c(0, label_index)
label_name = c("empty", label_name)
fmri_3dvisual_region(TCIU::phase2_pval, label_mask, label_index,
                    label_name, title = "Phase-2 p-values")

fmri_3dvisual_region(list(TCIU::phase2_pval,TCIU::phase3_pval), label_mask, label_index,
                    label_name, title = "Contrasting Phase 2 & Phase 3 p-values")

9 Statistical Inference on Kimesurfaces

9.1 Quantification of distances between kimesurfaces

This paper by Daubechies and colleagues shows one approach to quantify the similarities and differences between pairs of 2D surfaces (embedded in 3D) using their local structures and global information contained in inter-structure geometric relationships:

Conformal Wasserstein Distance (cW)
Conformal Wasserstein Neighborhood Dissimilarity Distance (cWn)
Continuous Procrustes Distance Between Surfaces

Below, we demonstrate one example of computing the Wasserstein distance and the Procrustes distance between a pair of kimesurfaces,

# Wasserstein distance (transport package)

#' Takes a time-series (1D vector), temporal longitudinal data, 
#' and converts it to complex-values kimesurface
#' 
#' @param datay A numerical vector.
#' @returns A list of 7 objects representing 
#'              [[1]] x-grid points of (2D) kimesurface magnitude
#'              [[2]] y-grid points (2D) of kimesurface magnitude
#'              [[3]] LT kimesurface (complex-number)
#'              [[4]] LT kimesurface magnitude
#'              [[5]] LT kimesurface Real-part
#'              [[6]] LT kimesurface Imaginary-part
#'              [[7]] LT kimesurface phase
#' @examples
#' testTS2KS <- timeSeries2KimeSurface(datay = sin(seq(-pi, pi, length.out=200)))
timeSeries2KimeSurface <- function(datay = sin(x1)) {
  Lout = 61
  range_limit <- 20
  x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
  y2 <- seq(from = 0, to = range_limit, length.out = Lout)
  
  x2_g = x2 %o% seq(1,1, length.out = Lout)
  y2_g =  seq(1,1, length.out = Lout) %o% y2

  z2_grid = x2 %o% y2
  
  argz = as.vector(x2_g + 1i*y2_g)
  LTz = NuLT(x1, datay, argz)
  rec1 = matrix(LTz,nrow = Lout)
  
  recm = abs(rec1)
  recr = Re(rec1)
  reci = Im(rec1)
  ph = reci/recm
  surf_color <- atan2(reci,recr)
  # resultKimesurface <- list()
  # resultKimesurface[[1]] <- x2_g  # x-grid points of (2D) kimesurface magnitude
  # resultKimesurface[[2]] <- y2_g  # y-grid points (2D) of kimesurface magnitude
  # resultKimesurface[[3]] <- rec1  # LT kimesurface (complex-number)
  # resultKimesurface[[4]] <- recm  # LT kimesurface magnitude
  # resultKimesurface[[5]] <- recr  # LT kimesurface Real-part
  # resultKimesurface[[6]] <- reci  # LT kimesurface Imaginary-part
  # resultKimesurface[[7]] <- surf_color # LT kimesurface phase
  p2D <- plot_ly(z=matrix(recm, nrow = sqrt(length(recm))), type="heatmap")
  p3D <- plot_ly(x=x2_g, y=y2_g, z=matrix(recm, nrow = sqrt(length(recm))),
                 type="surface")
  #  plot_ly(x=x2_g, y=y2_g, z=matrix(log(recm), nrow = sqrt(length(recm))), type="surface")
  return(list("x_grid"=x2_g, "y_grid"=y2_g, "complexLT_KS"=rec1, "magLT_KS"=recm,
              "realLT_KS"=recr, "imaginaryLT_KS"=reci, "phaseLT_KS"=surf_color,
              "plot_ly2D" = p2D, "plot_ly3D" = p3D))
}

kimesurface_SplOnAvg <- timeSeries2KimeSurface(datay = spl_onAvg$y)
kimesurface_SplOffAvg <- timeSeries2KimeSurface(datay = spl_offAvg$y)

# Compute the Wasserstein Distance Between the On and Off kimesurfaces 
# convert the matrices LT to pgrid
# Input arrays must be prob mass functions (>=0)
wDistanceMagLT_KS_p1 <- wasserstein(pgrid(kimesurface_SplOnAvg$magLT_KS),
            pgrid(kimesurface_SplOffAvg$magLT_KS), p=1, prob=FALSE)
wDistanceMagLT_KS_p2 <- wasserstein(pgrid(kimesurface_SplOnAvg$magLT_KS),
            pgrid(kimesurface_SplOffAvg$magLT_KS), p=2, prob=FALSE)

wDistanceRealLT_KS_p1 <- wasserstein(pgrid(abs(kimesurface_SplOnAvg$realLT_KS)),
            pgrid(abs(kimesurface_SplOffAvg$realLT_KS)), p=1, prob=FALSE)
wDistanceRealLT_KS_p2 <- wasserstein(pgrid(abs(kimesurface_SplOnAvg$realLT_KS)),
            pgrid(abs(kimesurface_SplOffAvg$realLT_KS)), p=2, prob=FALSE)

# Plot p1 and p2 Wasserstein Distances Between the On and Off kimesurfaces
WassersteinDistances <- c("Wasserstein P1 Distance", "Wasserstein P2 Distance")
KimeSurfaceMagnitude <- c(wDistanceMagLT_KS_p1, wDistanceMagLT_KS_p2)
KimeSurfaceRealPart  <- c(wDistanceRealLT_KS_p1, wDistanceRealLT_KS_p2)
KS_DataFrame <- data.frame(WassersteinDistances, KimeSurfaceMagnitude, KimeSurfaceRealPart)

plot_ly(KS_DataFrame, x = ~WassersteinDistances, y = ~KimeSurfaceMagnitude,
               type = 'bar', name = 'kimesurface (Complex) Magnitude') %>% 
  add_trace(y = ~KimeSurfaceRealPart, name = 'kimesurface Real Part') %>% 
  layout(title="Wasserstein Distances Between the kimesurfaces \n corresponding to the (avg) ON and OFF fMRI Stimuli",
         yaxis = list(title = 'Wasserstein Distances'), 
         xaxis=list(title="p-norm (1 or 2)"), barmode = 'group')

### Procrustes distance  (shapes package)
# procdist(x, y,type="full",reflect=FALSE)
# type - string indicating the type of distance; "full" full Procrustes distance, "partial" partial Procrustes distance, "Riemannian" Riemannian shape distance, "sizeandshape" size-and-shape Riemannian/Procrustes distance

distanceMagLT_KS_Procrustes <- procdist(kimesurface_SplOnAvg$magLT_KS,
            kimesurface_SplOffAvg$magLT_KS, type="full")
distanceMagLT_KS_Riemannian <- procdist(kimesurface_SplOnAvg$magLT_KS,
            kimesurface_SplOffAvg$magLT_KS, type="Riemannian")

distanceRealLT_KS_Procrustes <- procdist(kimesurface_SplOnAvg$realLT_KS,
            kimesurface_SplOffAvg$realLT_KS, type="full")
distancerealLT_KS_Riemannian <- procdist(kimesurface_SplOnAvg$realLT_KS,
            kimesurface_SplOffAvg$realLT_KS, type="Riemannian")

ProcrustesDistances <- c("Procrustes Distance", "Riemannian Distance")
KimeSurfaceMagnitude <- c(distanceMagLT_KS_Procrustes, distanceMagLT_KS_Riemannian)
KimeSurfaceRealPart  <- c(distanceRealLT_KS_Procrustes, distancerealLT_KS_Riemannian)

KS_DataFrame <- data.frame(ProcrustesDistances, KimeSurfaceMagnitude, KimeSurfaceRealPart)

plot_ly(KS_DataFrame, x = ~ProcrustesDistances, y = ~KimeSurfaceMagnitude,
               type = 'bar', name = 'kimesurface (Complex) Magnitude') %>% 
  add_trace(y = ~KimeSurfaceRealPart, name = 'kimesurface Real Part') %>% 
  layout(title="Procrustes & Riemannian Distances Between the kimesurfaces \n corresponding to the (avg) ON and OFF fMRI Stimuli",
         yaxis = list(title = 'Topological Shape Distances'), 
         xaxis=list(title="Procrustes & Riemannian Distances"), barmode = 'group')

9.2 `R` functionality and `R` packages for surface distance calculations

9.3 Kimesurface Tensor Linear Modeling (TLM)

Capitalizing on previous work, Lock’s TLM, DOI: 10.1080/10618600.2017.1401544 and Spacekime Analytics, DOI: 10.1007/s00521-021-06789-8, we can perform tensor linear modeling (TLM) on kimesurfaces. As parametric 2D manifolds, kimesurfaces are second-order tensors, or multidimensional arrays extending to higher dimensions the multivariate linear models, where measured scalar (or vector) outcomes $Y_n$) are regressed on two-way (design) matrices $X_{n\times k}$ $(n\ (samples) \times k\ (features))$ representing the observed inputs (predictors or covariates).

The tensor representations of kimesurfaces encoding the repeated measurement longitudinal data are stored as multidimensional arrays with dimensions reflecting, time, repetition, subjects, clinical traits, brain regions, frequencies, and other observable characteristics that may play a role in describing the output (tensor), i.e., complex-valued ($\mathbb{C}$) kimesurface intensity. Suppose we have a number of kimesurfaces, one kimesurface corresponding to one subject, parameterized on a 2D grid $100 (timepoints) \times 100 (repeats)$. The entire collection of kimesurfaces may be a 4D array (4th order tensor) of dimensions

\[\underbrace{\overbrace{\mathcal{X}}^{kimesurface}}_{design\ tensor}= \left (\underbrace{traits}_{phenotypes}, \underbrace{\overbrace{100}^{timepoints}, \overbrace{100}^{repeats}}_{\kappa\in\mathbb{C}}, \overbrace{2}^{f(\kappa)\in\mathbb{C}} \right ).\]

Assume also that the outcome tensor $\mathcal{Y}_{N, Q_1, \cdots, Q_M}$, of dimensions $\left (\underbrace{traits}_{phenotypes}, attributes\right )$, is predicted using the observed kimesurface tensor $\mathcal{X}_{N, P_1, \cdots, P_L}$. Then, the tensor linear model fitting $\mathcal{Y}= \langle \mathcal{X},\ \mathcal{B}\rangle_L + \mathcal{E}$ allows us to estimate the two tensors $\mathcal{B}_{P_1, \cdots, P_L; Q_1, \cdots, Q_M}$ , representing the TLM effect-sizes (tensor coefficients), and $\mathcal{E}_{N, Q_1, \cdots, Q_M}$, denoting the residual error tensor.

The (contracted) tensor product, $\langle \mathcal{U},\ \mathcal{V}\rangle_L$ operates on an order ($K+L$) tensor, representing the TLM_{I_{1} I_{K} P_{1} P_{L}}$ and an order ($L+M$) tensor $\mathcal{V}_{P_{1} \times \cdots \times P_{L} \times Q_{1} \times \cdots \times Q_{M}}$ and returns another tensor of order ($K+M$), contracting the dimensions $\left (P_{1} \times \cdots \times P_{L}\right )$. This is different from the general outer product of the same two tensors $\mathcal{U}$ and $\mathcal{V}$, which yields a different tensor of dimensions $(K +2L+M)$. The tensor product contracts the common $L$ dimensions shared by the two tensors $\mathcal{U}$ abd $\mathcal{V}$:

\[\langle{\mathbb{U}},{\mathbb{V}}\rangle_{L} \left[ {i_{1} , \ldots , i_{K} ;q_{1} , \ldots , q_{M} } \right] = \mathop \sum \limits_{{p_{1} = 1}}^{{P_{1} }} \cdots \mathop \sum \limits_{{p_{L} = 1}}^{{P_{L} }} {\mathbb{U}}\left[ {i_{1} , \ldots , i_{K} ;p_{1} , \ldots , p_{L} } \right]{\mathbb{V}}\left[ {p_{1} , \ldots , p_{L} ;q_{1} , \ldots , q_{M} } \right] \ .\]

In TLM, the first $L$ modes of the effects tensor $\mathcal{B}$ contract the dimensions of the covariate tensor, $\mathcal{X}$, and do not show in the outcome tensor $\mathcal{Y}$. On the other hand, the last $M$ modes of the effects tensor $\mathcal{B}$ expand along the modes of the outcome tensor $\mathcal{Y}$, and are not present in the covariate (data) tensor $\mathcal{X}$. When sufficient data are available to fit the tensor linear model, we can the estimate the effect tensor $\mathcal{B}$ and the residual error tensor $\mathcal{E}$.

The fitted model facilitates prediction, regression, forecasting, and extrapolation of the outcome tensor $\mathcal{Y}$ using new or prospective data, $\mathcal{X}'$. More specifically, For any outcome tensor index-set $(q_1, \cdots, q_M)$ and $\forall\ 1\leq n\leq N$, the TLM prediction of the corresponding tensor outcome is expressed as follows

\[\underbrace{\mathcal{Y}' \left({n,q_{1} , \ldots , q_{M} } \right )}_{TLM\ prediction} \cong \mathop \sum \limits_{{p_{1} = 1}}^{{P_{1} }} \cdots \mathop \sum \limits_{{p_{L} = 1}}^{{P_{L} }} \underbrace{{\mathcal{X}'}\left ({n,p_{1} , \ldots , p_{L} } \right )}_{new\ data\ tensor} \ \underbrace{\widehat{\mathcal{B}} \left ( {p_{1} , \ldots , p_{L} ;q_{1} , \ldots , q_{M} } \right )}_{Estimated\ effects\ tensor} \ .\]

Once the effects-tensor ${\widehat{\mathcal{B}}}$ is estimated, the residual tensor, $\widehat{\mathcal{E}}_{N, Q_1, \cdots, Q_M}$, is computed by

\[\widehat{\mathcal{E}}= \mathcal{Y} - \langle \mathcal{X},\ {\widehat{\mathcal{B}}}\rangle_L\ .\]

Below we show two simple TLM examples. The first one is a simulation involving $\sin()$, the sine function with different frequencies and the second one uses the On and Off fMRI dataset. In both cases, the first step involves transforming the raw longitudinal time-courses to kimesurfaces. Then, the outcomes (sine function period, in the first example, and On/Off fMRI stimuli, in the second example) are predicted by the fitted tensor linear model (TLM). 2D visualizations of the observed vs. predicted outcomes are shown as scatter plots.

These examples are intentionally oversimplified to demonstrate the mechanics of spacekime analytics using TLM. In practice, research investigators need to draft expropriate analytical protocols that reflect the underlying assumptions (e.g., tensor rank decomposition), study design, research hypotheses, expected outcomes, and scientific validity.

library(MultiwayRegression)
# 1. Simulation example

# TO DO: generalize this function to allow inserting any dataY/function input!!!
timeSeries2KimeSurfaceTensor <- function(datay = sin(x1), freq = c(1:10)) {
  Lout = 61
  range_limit <- 20
  x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
  y2 <- seq(from = 0, to = range_limit, length.out = Lout)
  
  x2_g = x2 %o% seq(1,1, length.out = Lout)
  y2_g =  seq(1,1, length.out = Lout) %o% y2

  z2_grid = x2 %o% y2
  
  argz = as.vector(x2_g + 1i*y2_g)
  
  # compute all 10 sine-frequency (kimesurface) simulated samples
  arr3D <- array(0, c(freq[length(freq)], Lout, Lout))
  for (i in freq) {
    LTz = NuLT(x1, sin(i*x1), argz)
    arr3D[i, , ] = matrix(LTz, nrow = Lout)
  }
  return(arr3D)  # Complex-valued!!!!
}

# str(X) #  num [1:10, 1:61, 1:61]
# str(Y) #  num [1:10, 1:2] univariate scalar 1:10 (sine frequency)
ksComplex <- timeSeries2KimeSurfaceTensor(datay = sin(x1), freq = c(1:10))
Y <- cbind(c(1:10), c(1,0,0,0,0,0,0,0,0,0))

# Mind the TLM regularization (Ridge penalty, L2) to prevent singularities
tlModel <- rrr(X=Re(ksComplex), Y=Y, R=1, lambda = 0.01) # Fit rank 1 model with Ridge regularization
Y_pred <- ctprod(A=Re(ksComplex), B=tlModel$B, 2) 

plot_ly(x=Y[ , 1], y=Y_pred[,2], type="scatter", mode="markers+lines") %>%
  layout(title="Simulated Sine Function Kimesurface TLM", 
         xaxis=list(title="Actual Sine Period"), 
         yaxis=list(title="TLM-predicted Period"))

# 2. fMRI Data example
fMRI2KimeSurfaceTensor <- function(datay = sin(x1)) {
  Lout = 61
  range_limit <- 20
  x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
  y2 <- seq(from = 0, to = range_limit, length.out = Lout)
  
  x2_g = x2 %o% seq(1,1, length.out = Lout)
  y2_g =  seq(1,1, length.out = Lout) %o% y2

  z2_grid = x2 %o% y2
  
  argz = as.vector(x2_g + 1i*y2_g)
  
  # compute all 10 sine-frequency (kimesurface) simulated samples
  arr3D <- array(0, c(1, Lout, Lout))
  LTz = NuLT(x1, datay, argz)
  arr3D[1, , ] = matrix(LTz, nrow = Lout)
  
  return(arr3D)  # Complex-valued!!!!
}

ks1 <- fMRI2KimeSurfaceTensor(datay = spl_onAvg$y)
ks2 <- fMRI2KimeSurfaceTensor(datay = spl_offAvg$y)
Lout <-  61
arrKS3D <- array(c(ks1, ks2), c(2, Lout, Lout))

# str(X) #  num [1:2, 1:61, 1:61]
# str(Y) #  num [1:2, 1:2] univariate scalar 1:2 (fMRI ON and OFF stimuli)
Y <- cbind(c(1:2), c(1,0))  # On and Off fMRI stimulation paradigm

# Mind the TLM regularization (Ridge penalty, L2) to prevent singularities
tlModel <- rrr(X=Re(arrKS3D), Y=Y, R=1, lambda = 0.01) # Fit rank 1 model with Ridge regularization
Y_pred <- ctprod(A=Re(arrKS3D), B=tlModel$B, 2) 

plot_ly(x=Y[ , 1], y=Y_pred[,2], type="scatter", mode="markers+lines") %>%
  layout(title="Simple fMRI Kimesurface TLM", 
         xaxis=list(title="Actual On/Off Stimulus"), 
         yaxis=list(title="TLM-predicted On/Off Paradigm"))

Of course, the most interesting applications of spacekime analytics reflect multiple repeated measurement data. For instance, we can use the $\sin(\cdot)$ function to simulate 10, or 100, noisy signals, corrupted with white noise, that represent multiple repeated observations. Then, we can contrast the extrapolation results from extending the predictions forward in time by using various model-based and model-free methods, including spacekime analytics. For instance, the out-of-temporal-scope spacekime analytic predictions may be based on applying the Inverse Laplace Transform to the kimesurface models obtained by first interpolating in the spacekime domain, and then transforming the kimesurface models back into spacetime for validation and contrasting against the predictions of alternative methods.

9.4 Topological Kimesurface Analyses (TKA)

Topological Data Analysis (TDA) facilitates statistical inference, analytics and modeling of point cloud data via extraction and quantification of topological information about the underlying space that the observed data is embedded in.
Examples of such topological characterizations include estimation of distance functions, distance to a measure functions, kNN density estimators, kernel density estimators, and kernel distances. TDA aims to capture salient topological features of the sub- or super-level sets of these functions by quantifying the embedding space persistent homology.

There are a number of algorithms for computing the persistent homology of the Rips filtration, which represents just one for the persistent homology of the sub- and super-level sets of any functions evaluated over a grid of points, e.g, C++ libraries GUDHI, Dionysus, and PHAT, and R packages TDA, TDAstats and tdaverse.

Obtaining reliable estimates of the persistence homology of kimesurfaces facilitates diagrammatic depictions and quantitative analytics of the induced topological characterizations (functions, distances, and distributions). For repeated measurements longitudinal data (such as time-series)Specifically, topological kimesurface analyses (TKA) enables the identification and explication of spatiotemporal structure of the corresponding kimesurface representations.

9.4.1 Quantitative TKA

The examples below use the TDA framework to showcase the basic steps in TKA/TDA and describe some quantitative kimesurface analytical protocols. Examples of topological quantization functions that will be demonstrated include the following:

The distance function is defined $\forall y \in \mathbb{R}^d$ by $\Delta(y) = \inf_{x\in X} |x − y|_{L_2}$.
Given a probability measure $P$ , the distance to measure (DTM) function is defined $\forall y \in \mathbb{R}^d$ by $d_{m_o}(y)=\left ( \frac{1}{m_o} \int_{0}^{m_o} {\left ( G^{-1}_y(u)\right )^r\ du} \right )^{\frac{1}{r}}$, where $G_y(u)= P\left (|X − y|\leq t\right )$, and the tuning parameters (smoothness) $m_o\in (0, 1)$ and $r \in [1,\infty)$ control the (smoothness/analytic) properties of the DTM. As $m_o \uparrow$, the DTM function becomes smoother. The DTM is a smoother version of the distance function. For discrete kimesurfaces, $X = \{x_1, \cdots , x_n\}$, the distance to measure function is defined by $\hat{d}_{m_o}(y)=\left ( \frac{1}{k} \sum_{x_i\in N_k(y)} {\left ( ||x_i - y||^r\right )} \right )^{\frac{1}{r}}$, where $k = \lceil m_o \times n\rceil$ and the $k$-nearest neighbors $y\in N_k(y)$ surround the discrete kimesurfaces, $X = \{x_1, \cdots , x_n\}$.
$\forall y \in \mathbb{R}^d$, the $k$-nearest neighbor density estimator, is defined by $\hat{\delta}_k(y)=\frac{k}{n\ \nu_d\ r_k^d(y)}$, where $\nu_d$ is the volume of the Euclidean $d$-dimensional unit ball and $r_k^d(\cdot)$ is the Euclidean distance from point $y$ to its $k^{th}$ closest neighbor among the points on the kimesurface $X$.
$\forall y \in \mathbb{R}^d$, the Gaussian Kernel Density Estimator (KDE), is defined by $\hat{p}_h(y)=\frac{1}{n(h\sqrt{2\pi})^d} \sum_{i=1}^{n}{e^{-\frac{||y-x_i||_2^2}{2h^2}}}$, where $h$ is the kernel smoothing parameter.
Density Clustering using the function TDA::clusterTree() provides a mechanism to segment, partition, and classify kimesurfaces. Suppose $f$ is the (kimesurface) density of a (kimesurface) probability distribution $P$ and we have obtained/observed a sample of kimesurfaces $X = \{x_1, \cdots , x_n\}\subset \mathbb{R}^d$. Given a threshold value $\lambda \gt 0$, denote the corresponding super-level set of $f$ by $L_f (\lambda) \equiv closure\left (\{x\in \mathbb{R}^s : f (x) \gt \lambda\}\right )$, whose $d$-dimensional subsets are high-density regions.

The high-density clusters of $P$ are the maximally connected subsets of $L_f (\lambda)$. The evolution and the hierarchy of the high-density clusters of $P$ can be tracked by the level-sets dynamics corresponding to $\lambda\in \mathbb{R}^+$. This evolution $0\leftarrow \lambda\rightarrow \infty$ yields the cluster density tree of the distribution $P$, defined as the collection of sets $T\equiv \{L_f(\lambda),\ \forall\lambda \geq 0\}$ that satisfy the (hierarchical-clustering) tree property, i.e., $\forall A, B \in T$, either $A \subset B$, or $B \subset A$, or $A\cap B = \emptyset$. This $\lambda$-tree construction can be augmented by the $\alpha$-tree and $\kappa$-tree constructions, which encode the probability content of each tree branch, rather than the probability density level. Cluster trees are useful in high dimensional data analytics (HDDA) such as feature rich point-cloud datasets and kimesurfaces whose spatiotemporal organization (structure) may be challenging to explicate in the native space.

library(TDA)
# 1. Simulation example
# str(X) #  num [1:10, 1:61, 1:61]
# str(Y) #  num [1:10, 1:2] univariate scalar 1:10 (sine frequency)
# ksComplex <- timeSeries2KimeSurfaceTensor(datay = sin(x1), freq = c(1:10))
# Y <- cbind(c(1:10), c(1,0,0,0,0,0,0,0,0,0))

#  1.1 Plot the kimesurfaces of the sine-functions corresponding to different periods
#  1.1.A Plot the Native Sine functions
plot_ly(x=x1, y=sin(x1), type="scatter", mode="lines", name="Period=1") %>%
  add_trace(y=sin(4*x1), name="Period=4") %>%
  add_trace(y=sin(9*x1), name="Period=9") %>%
  layout(title="Original Sine-functions with Different Periods", 
         legend = list(orientation='h', y=-0.2))

#  1.1.B Plot the Corresponding sine-kimesurfaces
plot_ly(z=Re(ksComplex)[1, , ], type="surface", name="Period=1") %>%
  add_surface(z=1+Re(ksComplex)[4, , ], type="surface", name="Period=4", opacity=0.4) %>%
  add_surface(z=-1+Re(ksComplex)[9, , ], type="surface", name="Period=9", opacity=0.4) %>%
  layout(title="Kimesurfaces of Sine-functions with Different Periods (1,4,9)") %>%
  hide_colorbar()

# 1.2 Compute the distance function for each point on the Grid
# Re(ksComplex)[1, , ]) is the Real part of the kimesurface of sin(x)

# Construct the 2D kimesurface domain grid
Lout = 61
range_limit <- 20
Xseq <- seq(from = 0.1, to = range_limit, length.out = Lout)
Yseq <- seq(from = 0, to = range_limit, length.out = Lout)

x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
y2 <- seq(from = 0, to = range_limit, length.out = Lout)
ksGrid <- expand.grid(x2, y2)

# Unwind the complex-value 2D kimesurface over the 61*61 grid into
# a 2D a matrix (m=61*61) * (d=2), m is the number of points on the kimesurface X
# and d is the dimension of the space, corresponding to 
# the kimesurface Re(KS) and Im(KS) parts
ksVector <- c(ksComplex[1, , ])
ksX <- cbind(R=Re(ksVector), I=Im(ksVector))

ksDistance <- distFct(X = ksX, Grid = ksGrid)
# plot_ly(z=matrix(ksDistance, ncol=length(Yseq), nrow=length(Xseq)), 
#         type="surface", name="Distance Function")

# 1.3 distance to measure (DTM) function
ksDTM <- dtm(X = ksX, Grid = ksGrid, m0 = 0.1)
# plot_ly(z=matrix(ksDTM, ncol=length(Yseq), nrow=length(Xseq)), 
#         type="surface", name="Distance to Measure (DTM)")

# 1.4 k-nearest neighbor density estimator function
# ks_kNN <- knnDE(X = Re(ksComplex)[1, , ], Grid = z2_grid, k = 30)
# ksVector <- c(ksComplex[9, , ])
# ksX <- cbind(R=Re(ksVector), I=Im(ksVector))
ks_kNN <- knnDE(X = ksX, Grid = ksGrid, k = 10)
# plot_ly(z=matrix(ks_kNN, ncol=length(Yseq), nrow=length(Xseq)), 
#         type="surface", name="k-Nearest Neighbor (kNN) Density Estimator")


# 1.5 Gaussian Kernel Density Estimator (KDE) function
ksKDE <- kde(X = ksX, Grid = ksGrid, h = 2.5)
# plot_ly(z=matrix(ksKDE, ncol=length(Yseq), nrow=length(Xseq)), 
#         type="surface", name="Gaussian Kernel Density Estimator (KDE)")

# Joint 3D surface plot
plot_ly(z=matrix(ksDistance, ncol=length(Yseq), nrow=length(Xseq)), 
        type="surface", name="Distance Function", hovertemplate = paste('(%{x}, %{y})')) %>%
  add_surface(z=-8+matrix(ksDTM, ncol=length(Yseq), nrow=length(Xseq)), 
        name="Distance to Measure (DTM)", opacity=0.4, hovertemplate = paste('(%{x}, %{y})')) %>%
  add_surface(z=8+matrix(ks_kNN, ncol=length(Yseq), nrow=length(Xseq)), 
        name="k-Nearest Neighbor (kNN) Density Estimator", opacity=0.4, 
        hovertemplate = paste('(%{x}, %{y})')) %>%
  add_surface(z=16+matrix(ksKDE, ncol=length(Yseq), nrow=length(Xseq)), 
        name="Gaussian Kernel Density Estimator (KDE)", opacity=0.4, 
        hovertemplate = paste('(%{x}, %{y})')) %>%
  layout(title="Joint 3D Surface Plot (Distance, DTM, kNN & KDE)") %>%
  hide_colorbar()

# 2. Bootstrap Confidence Bands: (1 − \alpha) confidence band for a function 
# using the bootstrap algorithm
ksBootstrapBands <- 
  bootstrapBand(X = ksX, Grid = ksGrid, FUN = dtm, # options: distFct, kde, and dtm 
                B = 100, parallel=TRUE, alpha=0.1, m0=0.1)  # for KDE: h=0.3)

plot_ly(z=matrix(ksBootstrapBands$fun, ncol=length(Yseq), nrow=length(Xseq)), 
        type="surface", name="Distance-to-Measure (DTM) Function", hovertemplate = paste('(%{x}, %{y})')) %>%
  add_surface(z=-5+matrix(ksBootstrapBands$band[,1], ncol=length(Yseq), nrow=length(Xseq)), 
        name="Lower Bootstrap Band", opacity=0.1, hovertemplate = paste('(%{x}, %{y})')) %>%
  add_surface(z=5+matrix(ksBootstrapBands$band[,2], ncol=length(Yseq), nrow=length(Xseq)), 
        name="Upper Bootstrap Band", opacity=0.1, hovertemplate = paste('(%{x}, %{y})')) %>%
  layout(title="Bootstrap Confidence Bands for the DTM Kimesurface") %>%
  hide_colorbar()

# 3. Density Clustering
# The function clusterTree implements Algorithm 1 in Kent, Rinaldo, and Verstynen (2013).
# Let f be the density of the probability distribution P generating the observed sample 
# X = {x1, . . . , xn} ⊂ Rd. For a threshold value λ > 0, the corresponding super-level set 
# of f is Lf (λ) := cl({x ∈ Rs : f (x) > λ}), and its d-dimensional  subsets are 
# high-density regions. The high-density clusters of P are the maximal connected subsets of Lf (λ).
# By considering all the level sets simultaneously (from λ = 0 to λ = ∞), 
# we track the evolution and the hierarchy of the high-density clusters of P.
# This leads to the notion of the cluster density tree of P defined as the collection of sets 
# T := {Lf (λ), λ ≥ 0} satisfying the tree property:
# A, B ∈ T implies that A ⊂ B or B ⊂ A or A ∩ B = ∅. This is the λ-tree. 
# Alternatively, Kent et al. (2013) introduced the α-tree and κ-tree, which facilitate the
# interpretation of the tree by precisely encoding the probability content of each tree branch
# rather than the density level. Cluster trees are particularly useful for high dimensional data,
# whose spatial organization is difficult to represent.

ksTree    <- clusterTree(X=ksX, k=12, density = "knn", 
                         Nlambda=120, printProgress=TRUE)

## 0   10   20   30   40   50   60   70   80   90   100
## |----|----|----|----|----|----|----|----|----|----|
## ***************************************************

ksTreeKDE <- clusterTree(X=ksX, k=2, h=5.3, density="kde", printProgress=TRUE)

## 0   10   20   30   40   50   60   70   80   90   100
## |----|----|----|----|----|----|----|----|----|----|
## ***************************************************

# ksTreeKDE <- clusterTree(X=ksX, k=10, h=0.3, density="kde", printProgress=TRUE)

# Plot the lambda and kappa trees of the k Nearest Neighbor (kNN) density
# estimator, and the kernel density estimator (KDE)
par(mfrow = c(2, 3))

plot(ksTree, type = "lambda", main = "Kimesurface Lambda Tree (kNN)")
plot(ksTree, type = "kappa", main = "Kimesurface Lappa Tree (kNN)")
plot(ksX, pch = 19, cex = 0.6, xlim = c(-0.01, 0.01), ylim = c(-0.02, 0.0),
     main = "Kimesurface kNN Cluster Labels")
for (i in ksTree[["id"]]){
    points(matrix(ksX[ksTree[["DataPoints"]][[i]],],ncol = 2), 
           col=i, pch=19, cex = 0.6)
}

plot(ksTreeKDE, type = "lambda", main = "Kimesurface Lambda Tree (KDE)")
plot(ksTreeKDE, type = "kappa", main = "Kimesurface Kappa Tree (KDE)")
plot(ksX, pch = 19, cex = 0.6, xlim = c(-0.2, 0.2), ylim = c(-0.2, 0.05), 
     main = "Kimesurface KDE Cluster Labels")
for (i in ksTreeKDE[["id"]]){
    points(matrix(ksX[ksTreeKDE[["DataPoints"]][[i]],],ncol = 2), 
           col=i, pch=19, cex = 0.6)
}

9.4.2 Parametric and Nonparametric Statistical Inference on Kimesurfaces

Both classical parametric and nonparametric (e.g., permutation) tests can be used to identify the statistical significance between kimesurfaces. In general, as the kimesurface distributions and the kimesurface topological characterizations are relatively unknown, the parametric assumptions may be invalid and TDA/TKA-based statistical inference on persistent homology is best quantified using nonparametric tests. Suppose we have a pair of (point-clouds) kimesurfaces $X$ and $Y$ (think of the fMRI-ON and fMRI-Off time-series runs). recall that in classical statistical inference, we compare the data-driven estimates of specific population parameters (e.g., mean, median) of each population (e.g., fMRI activation vs. rest) using a null ($H_o$) and an alternative (research) ($H_a$) hypotheses:

\[H_o:\mu_X=\mu_Y\ \ vs.\ \ H_a:\mu_X \not= \mu_Y\ .\]

Suppose a function $T$ returns the persistent homology of a point-cloud (kimesurface). The kimesurface distribution parameters (e.g., medians) corresponding to the pair of fMRI paradigms (On vs. Off), are statistically contrasted using corresponding (proper) proxy measures of the underlying kimesurface distributions (e.g., sample means, medians, variance, etc.). Without a priori parametric assumptions, we can assess the following hypotheses using a nonparametric permutation test to explicate the statistical inference comparing the quantitative TKA descriptions of the pair of kimesurfaces:

\[H_o: T(X)=T(Y)\ \ vs.\ \ H_a:T(X) \not= T(Y)\ .\]

Note that the permutation-test calculations below are computationally expensive. To expedite the compilation (knitting) of the Rmd into HTML5, the following examples use only (homologous) portions of the pair of kimesurfaces. Hence, the actual statistical significance (e.g., p-values) are coarse approximations of the corresponding p-values. In practice, the exact p-values should always be used.

# Recall:
# ks1 <- fMRI2KimeSurfaceTensor(datay = spl_onAvg$y)
# ks2 <- fMRI2KimeSurfaceTensor(datay = spl_offAvg$y)

# Plot (again) the (average) On and Off fMRI kimesurfaces
plot_ly(z=log(500+Re(ks1[1, , ])), name="fMRI On Kimesurface (Real Part)", 
        type="surface", opacity=1.0, colorscale="Jet", 
        # surfacecolor=Im(ks1[1, , ]),  # Color=Im(KS)?
        hovertemplate = paste('(%{x}, %{y})')) %>%
    add_surface(z=log(600+Re(ks2[1, , ])), name="fMRI Off Kimesurface (Real Part)",
                opacity=0.8, # surfacecolor=Im(ks1[1, , ]),
                hovertemplate = paste('(%{x}, %{y})')) %>%
    layout(title="fMRI On and Off Paradigm Kimesurfaces: Height=Log(Re(KS))") %>%
    hide_colorbar()

plot_ly(z=(Re(ks1[1, , ])-Re(ks2[1, , ])), type="surface", colorscale="Jet",
        name="fMRI Diff (On-Off) Kimesurface (Real Part)",
        opacity=1.0, hovertemplate = paste('(%{x}, %{y})')) %>%
  layout(title="fMRI Diff (On-Off) Kimesurface (Real Parts)") %>%
  hide_colorbar()

# Unwind the complex-value 2D kimesurface over the 61*61 grid into
# a 2D a matrix (m=61*61) * (d=2), m is the number of points on the kimesurface X
# and d is the dimension of the space, corresponding to 
# the kimesurface Re(KS) and Im(KS) parts

###################################################################
######### For speed - 
#########  reduce the kimesurfaces from [1:61,1:61] to 1:20, [1:20]
###################################################################

plot_ly(z=(Re(ks1[1, 1:20, 1:20])-Re(ks2[1, 1:20, 1:20])), 
        type="surface", colorscale="Jet",
        name="reduced fMRI Diff (On-Off) Kimesurface (Real Part)",
        opacity=1.0, hovertemplate = paste('(%{x}, %{y})')) %>%
    layout(title="Reduced fMRI Diff (On-Off) Kimesurface (Real Parts)") %>%
    hide_colorbar()

ksVector_ON <- c(ks1[1, 1:20, 1:20])  # complex-values!!!
ksX_ON <- cbind(R=Re(ksVector_ON), I=Im(ksVector_ON))

ksVector_OFF <- c(ks2[1, 1:20, 1:20])  # complex-values!!!
ksX_OFF <- cbind(R=Re(ksVector_OFF), I=Im(ksVector_OFF))

# load TDAstats
library("TDAstats")

# calculate homologies for both kimesurfaces
persistentHomology_ON <- calculate_homology(ksX_ON, dim = 1)
persistentHomology_OFF <- calculate_homology(ksX_OFF, dim = 1)

# plot topological barcodes for both ON and OFF fMRI kimesurfaces
plt_ON <- ggplotly(plot_barcode(persistentHomology_ON) +
  ggtitle("On fMRI Kimesurface Persistent Homology") +
  xlim(c(0, 500))) 
plt_OFF <- ggplotly(plot_barcode(persistentHomology_OFF)+
  ggtitle("Off fMRI Kimesurface Persistent Homology") +
  xlim(c(0, 500)))

subplot(plt_ON, plt_OFF) %>% 
  layout(title = 'TKA: Topological Barcodes for On (left) and Off (right) fMRI Kimesurfaces')

# Stats Inference: Permutation test
# Remove [1:400,] restriction for real studies, this is just for speed of
# compiling the Rmd to HTML notes!
onn_off_PermTest <- permutation_test(ksX_ON, ksX_OFF, iterations=100) #, update=1)

# report the stat-significance, p-values, for 0-cycles and 1-cycles
print(paste0("On vs. Off fMRI Kimesurface permutation test p-value for 0-cycle is: ",
             onn_off_PermTest[[1]]$pvalue))

## [1] "On vs. Off fMRI Kimesurface permutation test p-value for 0-cycle is: 1"

print(paste0("On vs. Off fMRI Kimesurface permutation test p-value for 1-cycle is: ",
             onn_off_PermTest[[2]]$pvalue))

## [1] "On vs. Off fMRI Kimesurface permutation test p-value for 1-cycle is: 0.01"

In this simple example, the 1-cycle of the persistent homology represents a single topological quantification of the fMRI kimesurfaces corresponding to the On and Off stimuli. The statistical significance quantifying the evidence to reject the null hypotheses, $H_o: T(X)=T(Y)\ \ vs.\ \ H_a:T(X) \not= T(Y)$, is represented by the computed p-value= 0.01.

9.4.3 Background on Topological Data Analysis

Simplicial complex representation of point cloud data, persistent homology quantification, and topological data analysis (TDA) rely on some of the following basic concepts.

A $k$-simplex is denoted by $[v_0,v_1, \cdots ,v_k]$ and represents the convex hull of $k+1$ affinely independent points $v_0,v_1, \cdots ,v_k$. For example, the $0$-simplex $[v_0]$ is just the single point $v_0$, whereas the $1$-simplex $[v_0,v_1]$ is the edge (line-segment) between the vertices $v_0$ and $v_1$, and the $2$-simplex $[v_0,v_1,v_2]$ is the triangle bordered by the three $1$-simplices (edges) $[v_0,v_1]$, $[v_1,v_2]$, and $[v_0,v_2]$.
A $j$-simplex $[u_0, u_1, \cdots ,u_j]$ is a face of another $k$-simplex $[v_0,v_1, \cdots ,v_k]$ when $\{u_0, u_1, \cdots ,u_j\}\subset \{v_0,v_1, \cdots ,v_k\}$.
A simplicial complex $K$ is a countable set of simplices where (1) every face of a simplex in $K$ is also in $K$, and (2) when two simplices $\sigma_1,\sigma_2\in K$, their intersection is either empty or a face of both simplices $\sigma_1,\sigma_2$.
For a finite simplicial complex $K$, a simplicial k-chain is a formal linear combination (over $\mathbb{C}$) of $k$-simplices in $K$. The set of $k$-chains forms a vector space $C_k(K)$. For each $k$-simplex, and by a linear extension, for each $k$-chain, the boundary map $\partial_k: C_k(K) \to C_{k-1}(K)$ is defined by

\[\partial_k([v_0, v_1, \cdots ,v_k]) = \sum_{j=0}^k { (-1)^j[v_0,\cdots ,\hat{v_j}, \cdots ,v_k]}\ \overbrace{\underbrace{=}_{-1\equiv +1}}^{\text{in }\mathbb{Z}_2}\ \sum_{j=0}^k [v_0,\cdots ,\hat{v_j}, \cdots ,v_k]\ ,\]

The second equality above follows since $\mathbb{Z}_2$ is the coefficient group. For the boundary map, $\partial_k$, the elements of $B_k(K) = {\text {img}}\ \partial _{k+1}$ are the boundaries and the elements of $Z_k(K)={\text {ker}}\ \partial_k$ are the cycles. For instance, the $1$-cycle represents the elements of $Z_1(K)={\text {ker}}\ \partial_1$, where $\partial_1: C_1(K) \to C_{0}(K)$.

Note that concatenating boundary maps yields elements in the kernel, $\partial_{k+1}\circ \partial_k=0$, i.e., $B_k(K) \subseteq Z_k(K)$. For the simplicial complex $K$, this property suggests the following $k^{th}$ homology group definition

\[H_k(K)\equiv Z_k(K)/B_k(K)\ .\]

The TKA topological barcode plots above may have different ranges for the horizontal and/or vertical axes, which may aid or complicate the direct comparison between the kimesurfaces. The persistent homology of the two kimesurfaces are compared using the Wasserstein metric (Wasserstein-1 is the Earth Mover’s Distance). The interpretation of the magnitude of the Wasserstein metric for a pair of kimesurfaces is difficult without a baseline, i.e., reference distribution. The nonparametric permutation test permutes the points between the ON and OFF fMRI kimesurfaces to derive the statistical significance, reflecting the power to reject the null hypothesis. Additional details about statistical TDA/TKA inference using permutation tests is provided in this article Hypothesis testing for topological data analysis, DOI: 10.1007/s41468-017-0008-7.

Spacekime Analytics (Time Complexity and Inferential Uncertainty)

Basic TCIU Protocol for Predictive Spacekime Analytics using Repeated-Measurement Longitudinal Data

SOCR Team

10/23/2024