---
title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)"
subtitle: "Basic TCIU Protocol for Predictive Spacekime Analytics using Repeated-Measurement Longitudinal Data"
author: "SOCR Team"
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output: 
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,cache = FALSE, warning = FALSE)
```

This [Spacekime TCIU Learning Module](https://tciu.predictive.space/) 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.


 
# Preliminary setup

TCIU and other `R` package dependencies ...

```{r setup1, message=FALSE, error=FALSE, warning=FALSE}
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)
```

# Longitudinal Data Import

In this case, we are just loading some exemplary fMRI data, which is [available here](https://drive.google.com/file/d/1Z5CvfotR032J37vHWf0_5ooAmDRpVVl8/view?usp=sharing).

```{r dataImport, message=FALSE, error=FALSE, warning=FALSE}
# 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
```

Alternatively, we can simulate synthetic ON/OFF fMRI data that can be saved as
a pair of arrays ON (stimulus) and OFF (rest) conditions.

```{r simDataGen, eval=FALSE, message=FALSE, error=FALSE, warning=FALSE}
library(plotly)
library(dplyr)

# 1. Simulate On/Off fMRI time-series data
simulate_fmri_data <- function(n_time = 80, n_runs = 5, subjects = 20) {
  set.seed(123)
  
  # Example: 'ON' data
  arr_on <- array(
    rnorm(n_time * n_runs * subjects, mean=0, sd=1),
    dim = c(n_time, subjects, n_runs)
  )
  # 'OFF' could be something different. For demonstration, let's shift it
  arr_off <- array(
    rnorm(n_time * n_runs * subjects, mean=1, sd=1),
    dim = c(n_time, subjects, n_runs)
  )
  
  list(ON = arr_on, OFF = arr_off)
}

# 2. Display the On and Off series
# This function converts the 3D array data into a tidy data frame format
# Calculates the mean time series for each condition (ON/OFF)
# Creates an interactive plot with:
# 
# Individual run data shown with low opacity
# Mean condition data shown with thicker lines
# Interactive tooltips with details on hover
# Legend distinguishing between conditions and runs
# Appropriate axis labels and title
# Function to create interactive plotly visualization of fMRI data
plot_fmri_timeseries <- function(fmri_data, subjects = c(1), title="On/Off fMRI Time-Series Data") {
  # Extract data from the simulation
  on_data <- fmri_data$ON[ , subjects, ]
  off_data <- fmri_data$OFF[ , subjects, ]
  
  # fix the dimensions
  on_data  <- array(on_data, dim = c(dim(on_data)[1], length(subjects), dim(on_data)[2]))
  off_data <- array(off_data, dim = c(dim(off_data)[1], length(subjects), dim(off_data)[2]))
    
  n_time <- dim(on_data)[1]
  n_runs <- dim(on_data)[3]
  
  # Create a data frame for plotting
  plot_data <- data.frame()
  
  # Process ON data
  for (run in 1:n_runs) {
    temp_df <- data.frame(Time = 1:n_time, Value = on_data[, 1, run],
                          Condition = "ON", Run = paste("Run", run))
    plot_data <- rbind(plot_data, temp_df)
  }
  
  # Process OFF data
  for (run in 1:n_runs) {
    temp_df <- data.frame(Time = 1:n_time, Value = off_data[, 1, run],
                          Condition = "OFF", Run = paste("Run", run))
    plot_data <- rbind(plot_data, temp_df)
  }
  
  # Calculate mean values for each condition across runs
  mean_data <- plot_data %>%
    group_by(Time, Condition) %>%
    summarize(Mean_Value = mean(Value), .groups = 'drop')
  
  # Create interactive plot
  p <- plot_ly() %>%
    # Add individual run traces with lower opacity
    add_trace(data = plot_data, x = ~Time, y = ~Value, color = ~Condition,
               type = 'scatter', mode = 'lines', linetype = ~Run, opacity = 0.3,
               showlegend = TRUE, hoverinfo = 'text', text = ~paste('Condition:',
                  Condition, '<br>Run:', Run, '<br>Time:', Time, '<br>Value:', round(Value, 3))
    ) %>%
    # Add mean traces with higher opacity and wider lines
    add_trace(data = mean_data, x = ~Time, y = ~Mean_Value, color = ~Condition,
               type = 'scatter', mode = 'lines', line = list(width = 4),
               name = ~paste(Condition, "Mean"), hoverinfo = 'text',
               text = ~paste('Condition:', Condition, 'Mean', 
                   '<br>Time:', Time, '<br>Value:', round(Mean_Value, 3))
    ) %>%
    layout(title = title,
      xaxis = list(title = "Time Point"), yaxis = list(title = "Signal Value"),
      legend = list(orientation = 'h', y=-0.1), hovermode = "closest")
  return(p)
}

# 2.1 Generate simulated data
fmri_data <- simulate_fmri_data(n_time = 80, n_runs = 5, subjects = 20)
# 2.2 Create interactive plot
plot_fmri_timeseries(fmri_data=fmri_data, subjects = c(1), title="Simulated On/Off fMRI Time-Series Data")


# 3. The triple Kimesurface visualization function
visualize_kimesurfaces_tripleREV <- function(fmri_data) {
  t_min <- 1;     t_max <- dim(fmri_data$ON)[1]
  th_min<- 1;     th_max<- dim(fmri_data$ON)[2]*dim(fmri_data$ON)[3]
  
  t_seq     <- c(t_min:t_max)
  theta_seq <- c(th_min:th_max)
  
  # flatten the repeats within subject over runs (5) and across subjects (20) == th_max
  # into a single dimensions, random sampling for theta,
  reshapedOn <- matrix(fmri_data$ON, nrow = th_max, ncol = t_max)
  reshapedOff <- matrix(fmri_data$OFF, nrow = th_max, ncol = t_max)
  
  # 2D surface in (t,theta) space
  fig <- plot_ly() %>%
    # First: ON
    add_surface(
      z = reshapedOn,
      colorscale = "Viridis",
      showscale  = FALSE,
      name       = "ON",
      scene      = "scene1"
    ) %>%
    layout(
      scene = list(
        domain     = list(x = c(0, 0.33), y = c(0, 1)),
        aspectmode = "cube",
        xaxis = list(title = "time (t)"),
        yaxis = list(title = "Runs * Participants (θ)"),
        zaxis = list(title = "Intensity")
      )
    ) %>%
    # Second: OFF
    add_surface(
      z = reshapedOff,
      colorscale = "Viridis",
      showscale  = FALSE,
      name       = "OFF",
      scene      = "scene2"
    ) %>%
    layout(
      scene2 = list(
        domain     = list(x = c(0.33, 0.66), y = c(0, 1)),
        aspectmode = "cube",
        xaxis = list(title = "time (t)"),
        yaxis = list(title = "Runs * Participants (θ)"),
        zaxis = list(title = "Intensity")
      )
    ) %>%
    # Third: Overlay
    add_surface(
      z = reshapedOn,
      colorscale = "Viridis",
      opacity    = 0.5,
      showscale  = FALSE,
      name       = "ON (overlay)",
      scene      = "scene3"
    ) %>%
    add_surface(
      z = reshapedOff,
      colorscale = "Portland",
      opacity    = 0.8,
      showscale  = FALSE,
      name       = "OFF (overlay)",
      scene      = "scene3"
    ) %>%
    layout(
      scene3 = list(
        domain     = list(x = c(0.66, 1), y = c(0, 1)),
        aspectmode = "cube",
        xaxis = list(title = "time (t)"),
        yaxis = list(title = "Runs * Participants (θ)"),
        zaxis = list(title = "Intensity")
      )
    ) %>%
    layout(
      title      = "Kimesurface Visualization of fMRI ON and OFF",
      showlegend = TRUE,
      legend     = list(x=0.85, y=0.1),
      annotations= list(
        list(showarrow=FALSE, text='(ON)',  xref='paper', yref='paper',
             x=0.15, y=0.95, xanchor='center', yanchor='bottom',
             font=list(size=16)),
        list(showarrow=FALSE, text='(OFF)', xref='paper', yref='paper',
             x=0.5,  y=0.95, xanchor='center', yanchor='bottom',
             font=list(size=16)),
        list(showarrow=FALSE, text='(OVERLAY)',
             xref='paper', yref='paper', x=0.85, y=0.95,
             xanchor='center', yanchor='bottom',
             font=list(size=16))
      )
    ) %>%
    # Sync the cameras
    onRender("
      function(el) {
        el.on('plotly_relayout', function(d) {
          const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
          if (camera.length) {
            const scenes = Object.keys(el.layout).filter((key) => /^scene\\d*/.test(key));
            const new_layout = {};
            scenes.forEach(key => {
              new_layout[key] = {...el.layout[key], camera: {...d[camera]}};
            });
            Plotly.relayout(el, new_layout);
          }
        });
      }
    ")
  
  fig
}

# 4. Display the kime-surfaces
fmri_data <- simulate_fmri_data()
fig <- visualize_kimesurfaces_tripleREV(fmri_data)
fig
```

# Time-series graphs

## 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`.


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

```{r simpleFMRI_Plot, fig.height=10, fig.width=10, message=FALSE, error=FALSE, warning=FALSE}
# reference_plot = sample[[4]]
fmri_time_series(bigim1, c(44,42,33), is.4d = TRUE) #, ref = reference_plot)
```

# Kime-series/kimesurfaces (spacekime analytics protocol)

The following examples connect to several ongoing spacekime analytics, including
[kimesurface analyticity](https://www.socr.umich.edu/TCIU/HTMLs/Chapter3_Radon_NikodymDerivatives.html), [kimesurface regularization](https://www.socr.umich.edu/TCIU/HTMLs/Chapter4_Laplace_Transform_Timeseries_Kimesurfaces.html), and 
[distribution-time dynamics](https://www.socr.umich.edu/TCIU/HTMLs/Chapter3_ReproducingKernelHilbertSpaces.html).

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*.
 
 - [TCIU Chapter 3 Kime-Phase Tomography (KPT)](https://www.socr.umich.edu/TCIU/HTMLs/Chapter3_Kime_Phase_Problem.html) and [TCIU Appendix on Mapping Repeated Measurement Longitudinal Data (Time-series) to 2D Manifolds](https://www.socr.umich.edu/TCIU/HTMLs/Chapter6_TCIU_MappingLongitudinalTimeseries_2_Kimesurfaces.html) provide mechanisms for transforming longitudinal processes from time-courses to kime-surfaces (manifolds).

## 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.


## Function main step: Time-series to kimesurfaces Mapping

### 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.

```{r kimePhaseLaplacePriorPlot, message=FALSE, error=FALSE, warning=FALSE}
# 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)
```

### 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.

```{r fMRI_On_Off_Simulation, message=FALSE, error=FALSE, warning=FALSE}
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))]
```




### 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.

```{r intensityPreProcessing, message=FALSE, error=FALSE, warning=FALSE}
# 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))
```

### 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*).

```{r plot_lyHoverLabels, message=FALSE, error=FALSE, warning=FALSE}
# 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
  }
}
```



```{r plot_lyCartesianHoverText, message=FALSE, error=FALSE, warning=FALSE}
# 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)) 
#    }
# }
```

### 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$.

```{r polar2CartesianCoordTransform, message=FALSE, error=FALSE, warning=FALSE}
# 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)
```



```{r reinterpolateCartesian, message=FALSE, error=FALSE, warning=FALSE}
# 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()
```
### Cartesian representation

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

```{r plot_lyCartesianHoverText_Mod, message=FALSE, error=FALSE, warning=FALSE}
# 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))
   }
}
```

```{r ON_surface_plot, message=FALSE, error=FALSE, warning=FALSE}
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)
    # 
```




```{r OFF_surface_plot, message=FALSE, error=FALSE, warning=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()
```

```{r raw_difference, message=FALSE, error=FALSE, warning=FALSE}
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))
```



```{r ON_OFF_surface_plot, message=FALSE, error=FALSE, warning=FALSE}
#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.

### 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.

```{r longDataFrameMapping, message=FALSE, error=FALSE, warning=FALSE}
# datatable(fmri_kimesurface(bigim1_mod,c(44,42,33))[[1]])
```

### 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*.

```{r plot_lyKimeSurfaceDisplay, message=FALSE, error=FALSE, warning=FALSE}
# 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
```

## 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](https://en.wikipedia.org/wiki/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](https://www.socr.umich.edu/TCIU/HTMLs/Chapter3_Radon_NikodymDerivatives.html), [kimesurface regularization](https://www.socr.umich.edu/TCIU/HTMLs/Chapter4_Laplace_Transform_Timeseries_Kimesurfaces.html), and 
[distribution-time dynamics](https://www.socr.umich.edu/TCIU/HTMLs/Chapter3_ReproducingKernelHilbertSpaces.html).

```{r kimePhase_AR1_LaplacePriorPlot, message=FALSE, error=FALSE, warning=FALSE}
# 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*, `r attributes(phi_8_vec_AR1[[1]]$coef)$names[1]`, 
$\alpha_1(r)=$ `r phi_8_vec_AR1[[1]]$coef[1]`, and 
 - *Mean*, `r attributes(phi_8_vec_AR1[[1]]$coef)$names[2]`, $\mu(r)=$ `r phi_8_vec_AR1[[1]]$coef[2]`.
 - The corresponding $AR(1)$ model residuals were $\{\varepsilon_t(r)\}_{t=1}^{10} =\{$ `r as.vector(phi_8_vec_AR1[[1]]$residuals)` $\}$.

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.

```{r intensityPreProcessing1, message=FALSE, error=FALSE, warning=FALSE}
# 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))
```

### 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*).

```{r plot_lyHoverLabels1, message=FALSE, error=FALSE, warning=FALSE}
# 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
  }
}
```


### 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$.

```{r polar2CartesianCoordTransform1, message=FALSE, error=FALSE, warning=FALSE}
# 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)
```



```{r reinterpolateCartesian1, message=FALSE, error=FALSE, warning=FALSE}
# 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()
```

### Cartesian representation

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

```{r plot_lyCartesianHoverText_Mod_AR1, message=FALSE, error=FALSE, warning=FALSE}
# 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))
   }
}
```

```{r ON_surface_plot1, message=FALSE, error=FALSE, warning=FALSE}
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()
```


```{r OFF_surface_plot1, message=FALSE, error=FALSE, warning=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 = "Kime-Phase AR(1) Model: OFF kimesurface at a fixed voxel location",
           scene = list(xaxis = x, yaxis = y, zaxis = z)) %>% hide_colorbar()
```

```{r raw_difference1, message=FALSE, error=FALSE, warning=FALSE}
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()
```



```{r ON_OFF_surface_plot1, message=FALSE, error=FALSE, warning=FALSE}
#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()
```


## Analytical Time-series to kimesurface Transformation (*Laplace*)

### 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.

```{r analyticKimesurfaceTransformLaplace, message=FALSE, error=FALSE, warning=FALSE}
# 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*.

### Discrete LT

```{r discreteLT, message=FALSE, error=FALSE, warning=FALSE}
# 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.



```{r testLT_cosine, message=FALSE, error=FALSE, warning=FALSE}
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}\ .$$

```{r ground_truthLT_cosine, message=FALSE, error=FALSE, warning=FALSE}
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.

```{r test_cosine_numeric, message=FALSE, error=FALSE, warning=FALSE}
print(paste0("The average difference between the ground truth and the numerical integration is ",
             sum(abs(ground_truth -rec1))/length(ground_truth)))
```


```{r testLT_sine, message=FALSE, error=FALSE, warning=FALSE}
# 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
```

```{r}
# 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](https://socr.umich.edu/DSPA2/DSPA2_notes/DSPA_Appendix_10_Spatiotemporal_Interpolation.html) to find out how to regularize either 
[regularly (equally-spaced) sampled longitudinal data](https://socr.umich.edu/DSPA2/DSPA2_notes/DSPA_Appendix_10_Spatiotemporal_Interpolation.html#21_Regularly_spaced_sample) or 
[*irregularly* (unequally-spaced) sampled longitudinal data](https://socr.umich.edu/DSPA2/DSPA2_notes/DSPA_Appendix_10_Spatiotemporal_Interpolation.html#22_Irregularly_spaced_sample).

```{r testLT_sine1, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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);
#         }
#       });
#     }")
```




```{r co-rotated_plot, message=FALSE, error=FALSE, warning=FALSE}

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);
        }
      });
    }")
```


## 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.

### 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.

#### 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.

#### 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.

#### 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](https://en.wikipedia.org/wiki/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 the operator $\operatorname {sgn}(\omega)$ captures the sign of the argument, 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.

#### 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\mathbb{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.

#### 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.

### 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*.

#### Kimesurface Reconstructions via Basic Stitching in the Native-space

```{r Timeseries2KimesurfaceIntegralTransform_P1, message=FALSE, error=FALSE, warning=FALSE}
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);
        }
      });
    }")
```


#### 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*.
 
```{r Timeseries2KimesurfaceIntegralTransform_P2, message=FALSE, error=FALSE, warning=FALSE}
# 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.


# 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)](https://cran.r-project.org/web/packages/TCIU/index.html) 
to *customize* the [TCIU functionality (GitHub)](https://github.com/SOCR/TCIU) for 
the application specific needs!

## 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.

```{r example1, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
bigim1_mask<-bigim1_mod
for (i in 1:160) {
  bigim1_mask[,,,i]<-bigim1_mask[,,,i]*mask
}
```

### Manual examples

A plot *without* (spatio-temporal) *mask* restrictions.

```{r example2, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
#' @title interactive graph object of the fMRI image
#' @description fMRI image visualization method, based on package \code{plotly}.
#'
#' @param fmridata a 4D array contains information for the fMRI spacetime image. The data should only contain the magnitude for the fMRI image.
#' @param option The default is 'manually'. If choose 'auto', then this function will lead you to key in the space (x,y,z) parameters and time (time) parameter for this function to generate graphs.
#' @param voxel_location a 3D array indicating the spatial location of the brain. If option is auto, set the voxel_location as NULL.
#' @param time time location for the voxel
#'
#' @details
#' The function \code{fmri_image} is used to create images for front view, side view, and top view of the fMRI image.
#' When providing the 4D array of the fMRI spacetime image and input the x,y,z position of the voxel, 
#' three views of the fMRI image and the time series image of the voxel will be shown.
#'
#' @author SOCR team <\url{http://socr.umich.edu/people/}>
#'
#' @return an interactive graph object of the fMRI image created by \code{plotly}
#' @export
#'
#' @import plotly
#' @importFrom forecast auto.arima
#' @importFrom ggplot2 ggplot geom_hline geom_vline xlim ylim
#' 
#' @examples
#' fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask)
#' fmri_image(fmri_generate$fmri_data, option='manually', voxel_location = c(40,22,33), time = 4)
#' 

fmri_image = function(fmridata,
                      option = "manually",
                      voxel_location = NULL,
                      time = NULL) {
  
  # Check input dimensions
  if (length(dim(fmridata)) != 4) {
    stop("fmridata must be a 4D array")
  }
  
  xdim = dim(fmridata)[1]
  ydim = dim(fmridata)[2]
  zdim = dim(fmridata)[3]
  tdim = dim(fmridata)[4]
  
  if (option == "auto") {
    message("Please key in x, y, z, and time in sequence")
    x = 0
    y = 0
    z = 0
    t = 0
    
    while (x > xdim | x < 1 | x%%1 != 0) {
      x <- scan(n = 1)
      if (x > xdim | x < 1 | x%%1 != 0) {
        message("input x is not in range or is not an integer! Please retype!")
      }
    }
    while (y > ydim | y < 1 | y%%1 != 0) {
      y <- scan(n = 1)
      if (y > ydim | y < 1 | y%%1 != 0) {
        message("input y is not in range or is not an integer! Please retype!")
      }
    }
    while (z > zdim | z < 1 | z%%1 != 0) {
      z <- scan(n = 1)
      if (z > zdim | z < 1 | z%%1 != 0) {
        message("input z is not in range or is not an integer! Please retype!")
      }
    }
    while (t > tdim | t < 1 | t%%1 != 0) {
      t <- scan(n = 1)
      if (t > tdim | t < 1 | t%%1 != 0) {
        message("input time is not in range or is not an integer! Please retype!")
      }
    }
  } else if (option == "manually") {
    # Check manual inputs
    if (is.null(voxel_location) || length(voxel_location) != 3) {
      stop("voxel_location must be a vector of length 3 (x, y, z)")
    }
    if (is.null(time) || length(time) != 1) {
      stop("time must be a single integer")
    }
    
    x = voxel_location[1]
    y = voxel_location[2]
    z = voxel_location[3]
    t = time
    
    # Validate manual inputs
    if (x < 1 || x > xdim || x%%1 != 0) {
      stop(sprintf("x must be an integer between 1 and %d", xdim))
    }
    if (y < 1 || y > ydim || y%%1 != 0) {
      stop(sprintf("y must be an integer between 1 and %d", ydim))
    }
    if (z < 1 || z > zdim || z%%1 != 0) {
      stop(sprintf("z must be an integer between 1 and %d", zdim))
    }
    if (t < 1 || t > tdim || t%%1 != 0) {
      stop(sprintf("time must be an integer between 1 and %d", tdim))
    }
  }
  
  # Time series data processing
  try1 <- fmridata[x, y, z, ]
  
  # Detrend if pracma is available
  if (requireNamespace("pracma", quietly = TRUE)) {
    try1 <- pracma::detrend(try1, bp = seq(21, tdim, by = 20))
  } else {
    warning("Package 'pracma' not available. Skipping detrending.")
  }
  
  tstry1 <- ts(try1)
  
  # Smoothing
  if (requireNamespace("stats", quietly = TRUE)) {
    ksmth <- stats::ksmooth(1:tdim, tstry1, kernel = "normal", bandwidth = 5)
    smth <- stats::smooth(tstry1)
  } else {
    stop("stats package is required for smoothing")
  }
  
  # Create time series plot directly with plotly (replacing GTSplot)
  time_df <- data.frame(
    Time = 1:tdim,
    Original = as.numeric(tstry1),
    Smooth = as.numeric(smth),
    KSmooth = ksmth$y
  )
  
  # Create interactive time series plot
  TScore <- plot_ly(time_df) %>%
    add_trace(x = ~Time, y = ~Original, type = 'scatter', mode = 'lines',
              name = 'Original', line = list(color = '#3399FF')) %>%
    add_trace(x = ~Time, y = ~Smooth, type = 'scatter', mode = 'lines',
              name = 'Smooth', line = list(color = '#66FF33')) %>%
    add_trace(x = ~Time, y = ~KSmooth, type = 'scatter', mode = 'lines',
              name = 'KSmooth', line = list(color = '#FF6666')) %>%
    layout(
      title = paste0("Voxel (", x, ",", y, ",", z, ") Time Series"),
      xaxis = list(title = "Time"),
      yaxis = list(title = "Signal"),
      hovermode = 'x unified'
    )
  
  # Create slice plots
  # z slice (axial)
  zfMRI <- t(fmridata[, , z, t])
  zslice <- plot_ly(z = ~zfMRI, type = "contour") %>%
    layout(
      title = paste0("Z = ", z, " (Axial View)"),
      xaxis = list(title = "X", range = c(0, xdim)),
      yaxis = list(title = "Y", range = c(0, ydim)),
      shapes = list(
        list(
          type = "line",
          x0 = x, x1 = x, y0 = 0, y1 = ydim,
          line = list(color = "red", width = 2)
        ),
        list(
          type = "line",
          x0 = 0, x1 = xdim, y0 = y, y1 = y,
          line = list(color = "red", width = 2)
        )
      )
    )
  
  # x slice (sagittal)
  xfMRI <- t(fmridata[x, , , t])
  xslice <- plot_ly(z = ~xfMRI, type = "contour") %>%
    layout(
      title = paste0("X = ", x, " (Sagittal View)"),
      xaxis = list(title = "Y", range = c(0, ydim)),
      yaxis = list(title = "Z", range = c(0, zdim)),
      shapes = list(
        list(
          type = "line",
          x0 = y, x1 = y, y0 = 0, y1 = zdim,
          line = list(color = "red", width = 2)
        ),
        list(
          type = "line",
          x0 = 0, x1 = ydim, y0 = z, y1 = z,
          line = list(color = "red", width = 2)
        )
      )
    )
  
  # y slice (coronal)
  yfMRI <- t(fmridata[, y, , t])
  yslice <- plot_ly(z = ~yfMRI, type = "contour") %>%
    layout(
      title = paste0("Y = ", y, " (Coronal View)"),
      xaxis = list(title = "X", range = c(0, xdim)),
      yaxis = list(title = "Z", range = c(0, zdim)),
      shapes = list(
        list(
          type = "line",
          x0 = x, x1 = x, y0 = 0, y1 = zdim,
          line = list(color = "red", width = 2)
        ),
        list(
          type = "line",
          x0 = 0, x1 = xdim, y0 = z, y1 = z,
          line = list(color = "red", width = 2)
        )
      )
    )
  
  # Combine all 4 plots
  result <- subplot(
    TScore, 
    zslice, 
    xslice, 
    yslice,
    nrows = 2,
    titleY = TRUE,
    titleX = TRUE
  ) %>%
    layout(
      title = paste0("fMRI Visualization - Voxel (", x, ",", y, ",", z, "), Time = ", t),
      showlegend = TRUE
    )
  
  return(result)
}

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

A plot *with* (spatio-temporal) *mask* restrictions.

```{r example3, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
fmri_image(bigim1_mask, option = "manually", voxel_location = c(40,40,30), time = 4)
```

## Forecasting with time series

The function `fmri_ts_forecast()` fits the [ARIMA models](https://socr.umich.edu/DSPA2/DSPA2_notes/12_LongitudinalDataAnalysis.html) 
for each *voxel* (spatial volume element) location. This function is based on the `TSplot_gen()` function from package `TSplotly`. 

```{r example4, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
fmri_ts_forecast(smoothmod,c(41,44,33))
```

# 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.

## fMRI data simulation

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

```{r fMRI_SImulation1, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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))
```


## 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()`.


### Examples

```{r fMRI_Sim_2, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
# 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.

```{r pValueSummary, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
summary(pval1)
summary(pval2)
summary(pval3)
summary(pval4)
```


# Motor area visualization

## Visualization and comparison of p-value

### 3D visualization for p-value

```{r example5, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
pval1_3d = fmri_3dvisual(pval1, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="Accounting for HRF")

pval1_3d$plot
```

```{r example5_pvals, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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
```

### 2D visualization for p-value

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

Plot *without isocontours*.


```{r example6, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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*.

```{r example7, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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))
}
```

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

### 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.

```{r example8, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
fmri_pval_comparison_3d(list(pval1, pval2), mask, 
    				            list(0.05, 0.05), list("scale_p", "scale_p"), 
    				            multi_pranges=FALSE)
```

### 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.


```{r example9, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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)
```

# 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](https://github.com/SOCR/TCIU) may be used as *templates* to extend this spacekime functionality to *alternative data formats* and *different data structures.* 

## Phase 1: Detect Potential Activated ROI

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

```{r example10, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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
```

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

```{r example11, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
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")
```

## Phase 3: FDR Correction and Spatial Clustering

```{r example12, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8}
phase3_pval = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
                            spatial_cluster.thr = 0.05,
                            spatial_cluster.size = 5, 
                            show_comparison = FALSE)
```
                           
## 3D visualization based on the activated areas by regions

```{r example13, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8, fig.align = "center"}
fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, title = "Phase 1 p-values", rank = "value")
``` 


```{r example14, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8, fig.align = "center"}
fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, 5, title = "pPhase-1 top five p-values", rank = "value")
``` 


```{r example15, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8, fig.align = "center"}
# 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")
``` 


```{r example16, message=FALSE, error=FALSE, warning=FALSE, fig.height=12, fig.width=12, fig.align = "center"}
fmri_3dvisual_region(list(TCIU::phase2_pval,TCIU::phase3_pval), label_mask, label_index,
                    label_name, title = "Contrasting Phase 2 & Phase 3 p-values")
```

# Statistical Inference on Kimesurfaces

## Quantification of distances between kimesurfaces

[This paper by Daubechies and colleagues](https://doi.org/10.1073/pnas.1112822108) 
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,  

```{r WassersteinDistance2Kimesurfaces, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8, fig.align = "center"}
# 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')
```


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

 - [dist: Distance Matrix Computation](https://rdrr.io/r/stats/dist.html)
 - [Distance Sampling](https://distancesampling.org/R/), see the [vignettes](https://distancesampling.org/R/vignettes/Project-migration.Rmd)
 - [dsm: Density Surface Modelling of Distance Sampling Data](https://cran.r-project.org/web/packages/dsm/)
 - [topoDistance: Calculating Topographic Paths and Distances](https://cran.r-project.org/web/packages/topoDistance/), the the [vignettes](https://cran.r-project.org/web/packages/topoDistance/vignettes/)
 - [Procrustes shape or size-and-shape distance between two configurations](https://rdrr.io/cran/shapes/man/procdist.html)

## Kimesurface *Tensor Linear Modeling* (TLM)

Capitalizing on previous work, [Lock's TLM, DOI: 10.1080/10618600.2017.1401544](https://doi.org/10.1080/10618600.2017.1401544) and 
[Spacekime Analytics, DOI: 10.1007/s00521-021-06789-8](https://doi.org/10.1007%2Fs00521-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} \times \cdots \times I_{K} \times P_{1} \times \cdots \times 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](https://doi.org/10.1007/s00521-021-06789-8)
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.

```{r kimesurfaceTLM, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=8, fig.align = "center"}
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](https://www.socr.umich.edu/TCIU/HTMLs/Chapter4_Laplace_Transform_Timeseries_Kimesurfaces.html) 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. 

## *Topological Kimesurface Analyses* (TKA)

[Topological Data Analysis (TDA)](https://en.wikipedia.org/wiki/Topological_data_analysis)
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](https://en.wikipedia.org/wiki/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](https://cran.r-project.org/web/packages/TDA),
[TDAstats](https://cran.r-project.org/web/packages/TDAstats) and 
[tdaverse](https://github.com/tdaverse/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.

### Quantitative TKA

The examples below use the [TDA framework](https://cran.r-project.org/web/packages/TDA/vignettes/article.pdf) 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.

```{r topoKimesurfaceAnalytics, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=12, fig.align = "center"}
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)
ksTreeKDE <- clusterTree(X=ksX, k=2, h=5.3, density="kde", printProgress=TRUE)
# 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)
}
```


### 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](https://doi.org/10.1007/s41468-017-0008-7) 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.

```{r permutationKimesurfaceInference, message=FALSE, error=FALSE, warning=FALSE, fig.height=8, fig.width=12, fig.align = "center"}
# 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)

# TKA: Topological Barcodes for On (left) and Off (right) fMRI Kimesurfaces

# 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')

plot_barcode(persistentHomology_ON) +
  ggtitle("On fMRI Kimesurface Persistent Homology") +
  xlim(c(0, 500))
plot_barcode(persistentHomology_OFF)+
  ggtitle("Off fMRI Kimesurface Persistent Homology") +
  xlim(c(0, 500))

# 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))
print(paste0("On vs. Off fMRI Kimesurface permutation test p-value for 1-cycle is: ",
             onn_off_PermTest[[2]]$pvalue))
```

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=* `r onn_off_PermTest[[2]]$pvalue`.

### 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](https://doi.org/10.1007/s41468-017-0008-7).
 

<!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->




