---
title: "Time Complexity, Inferential uncertainty and Spacekime Analytics"
subtitle: "Appendix 03.1: Kime-Phase Tomography (KPT) - fMRI Simulation"
author: "SOCR Team"
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output:
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
  word_document:
    toc: true
    toc_depth: '3'
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amssymb}
  - \usepackage{physics}
  - \usepackage{mathrsfs}
  - \usepackage{bm}
---


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

# Introduction

This [TCIU Appendix](https://tciu.predictive.space/) implements the 
[Kime-Phase Tomography (KPT) algorithm](https://www.socr.umich.edu/TCIU/HTMLs/TCIU_SK_Appendix03_KPT.html) to 
map repeated fMRI time-series measurements to 2D kime-surfaces. We test the algorithm on synthetic fMRI data with ON (stimulus) and OFF (rest) conditions.

# Required Libraries

```{r libraries}
# Load required libraries
library(plotly)
library(dplyr)
library(tidyr)
library(circular)
library(MASS)
library(pracma)  # For numerical integration
library(ggplot2)
library(viridis)

# Set random seed for reproducibility
set.seed(42)
```

# Simulate fMRI Data

```{r simulate-data}
# ---------------------------- 
# Simulate fMRI Data
# ---------------------------- 
simulate_fmri_data <- function(n_participants = 15, n_runs = 20, tr = 2, duration = 300) {
  # Hemodynamic Response Function (HRF)
  hrf <- function(t) {
    a1 <- 6; a2 <- 12; b1 <- 0.9; b2 <- 0.9; c <- 0.35
    t^(a1-1)*b1^a1*exp(-b1*t)/gamma(a1) - c*(t^(a2-1)*b2^a2*exp(-b2*t)/gamma(a2))
  }
  
  # Experimental design parameters
  n_timepoints <- duration/tr
  block_length <- 30  # 30s blocks
  n_blocks <- duration/block_length
  
  # Storage arrays
  data <- list(
    ON = array(0, dim = c(n_participants, n_runs, n_timepoints)),
    OFF = array(0, dim = c(n_participants, n_runs, n_timepoints))
  )
  
  # Generate data for each participant and run
  set.seed(123)
  for(p in 1:n_participants) {
    for(r in 1:n_runs) {
      # Generate block design
      design <- rep(c(0,1), each = block_length/(2*tr))
      if(length(design) < n_timepoints) {
        design <- c(design, rep(design[length(design)], n_timepoints - length(design)))
      } else {
        design <- design[1:n_timepoints]
      }
      
      # Generate HRF response
      hrf_values <- hrf(seq(0, block_length/tr - 1, length.out = block_length/tr))
      hrf_response <- numeric(n_timepoints)
      
      # Convolve design with HRF
      for(i in 1:n_timepoints) {
        if(design[i] == 1) {
          end_idx <- min(i + length(hrf_values) - 1, n_timepoints)
          hrf_response[i:end_idx] <- hrf_response[i:end_idx] + 
            hrf_values[1:(end_idx - i + 1)]
        }
      }
      
      # Add participant-specific noise and drift
      noise <- cumsum(rnorm(n_timepoints, sd = 0.1))
      drift <- 0.01 * seq(n_timepoints)
      
      # Store responses
      data$ON[p, r, ] <- hrf_response + noise + drift + rnorm(n_timepoints, sd = 0.2)
      data$OFF[p, r, ] <- noise + drift + rnorm(n_timepoints, sd = 0.2)
    }
  }
  return(data)
}

# Generate synthetic fMRI data
fmri_data <- simulate_fmri_data(n_participants = 10, n_runs = 20, tr = 2, duration = 300)
```

## fMRI Simulated Data Visualization

```{r simulate-fMRI-data-Viz}
# ---------------------------- 
# Visualize Original fMRI Time-Series
# ---------------------------- 

# Function to create interactive time-series plot
create_timeseries_plot <- function(fmri_data, n_display = 5) {
  # Extract dimensions
  n_participants <- dim(fmri_data$ON)[1]
  n_runs <- dim(fmri_data$ON)[2]
  n_timepoints <- dim(fmri_data$ON)[3]
  time_vec <- seq(0, (n_timepoints-1)*2, by = 2)  # TR = 2
  
  # Randomly select indices to display
  set.seed(42)
  selected_indices <- sample(1:(n_participants * n_runs), n_display)
  
  # Create plot
  fig <- plot_ly()
  
  # Add ON condition traces
  for(i in 1:n_display) {
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    fig <- fig %>% add_trace(
      x = time_vec,
      y = fmri_data$ON[participant, run, ],
      type = 'scatter',
      mode = 'lines',
      name = paste0("ON - P", participant, "R", run),
      line = list(color = 'red', width = 2),
      opacity = 0.7,
      visible = (i <= 3)  # Show only first 3 by default
    )
  }
  
  # Add OFF condition traces
  for(i in 1:n_display) {
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    fig <- fig %>% add_trace(
      x = time_vec,
      y = fmri_data$OFF[participant, run, ],
      type = 'scatter',
      mode = 'lines',
      name = paste0("OFF - P", participant, "R", run),
      line = list(color = 'blue', width = 2),
      opacity = 0.7,
      visible = (i <= 3)  # Show only first 3 by default
    )
  }
  
  # Create buttons for navigation
  button_list <- list()
  
  # Add "Show All" button
  visible_all <- rep(TRUE, 2 * n_display)
  button_list[[1]] <- list(
    method = "restyle",
    args = list("visible", as.list(visible_all)),
    label = "Show All"
  )
  
  # Add individual selection buttons
  for(i in 1:n_display) {
    visible_vec <- rep(FALSE, 2 * n_display)
    visible_vec[c(i, i + n_display)] <- TRUE  # Show ON and OFF for this index
    
    idx <- selected_indices[i]
    participant <- ((idx - 1) %/% n_runs) + 1
    run <- ((idx - 1) %% n_runs) + 1
    
    button_list[[i + 1]] <- list(
      method = "restyle",
      args = list("visible", as.list(visible_vec)),
      label = paste0("P", participant, "R", run)
    )
  }
  
  # Update layout with buttons
  fig <- fig %>% layout(
    title = "fMRI Time-Series: ON (Red) vs OFF (Blue) Conditions",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "BOLD Signal"),
    updatemenus = list(
      list(
        type = "buttons",
        direction = "left",
        xanchor = "left",
        yanchor = "top",
        x = 0.01,
        y = 1.15,
        buttons = button_list
      )
    ),
    annotations = list(
      list(
        text = "Select Time-Series:",
        showarrow = FALSE,
        x = 0,
        y = 1.2,
        xref = "paper",
        yref = "paper",
        xanchor = "right",
        yanchor = "bottom"
      )
    )
  )
  
  return(fig)
}

# Create and display the time-series plot
ts_plot <- create_timeseries_plot(fmri_data, n_display = 10)
ts_plot
```

## Simulated fMRI Data Summary Statistics

```{r simulate-fMRI-data-Summary-Stats}
# ---------------------------- 
# Time-Series Statistics
# ---------------------------- 

# Compute statistics for each condition
compute_ts_stats <- function(data_array) {
  n_series <- dim(data_array)[1] * dim(data_array)[2]
  n_timepoints <- dim(data_array)[3]
  
  # Reshape to matrix
  reshaped <- matrix(data_array, nrow = n_series, ncol = n_timepoints)
  
  # Compute statistics
  mean_ts <- colMeans(reshaped)
  sd_ts <- apply(reshaped, 2, sd)
  
  # Compute temporal SNR
  tsnr <- mean_ts / sd_ts
  
  return(list(
    mean = mean_ts,
    sd = sd_ts,
    tsnr = tsnr,
    global_mean = mean(mean_ts),
    global_sd = mean(sd_ts)
  ))
}

stats_on <- compute_ts_stats(fmri_data$ON)
stats_off <- compute_ts_stats(fmri_data$OFF)

# Create comparison plot
time_vec <- seq(0, (length(stats_on$mean)-1)*2, by = 2)

stats_plot <- plot_ly() %>%
  # Mean traces
  add_trace(x = time_vec, y = stats_on$mean, type = 'scatter', mode = 'lines',
            name = 'ON Mean', line = list(color = 'red', width = 3)) %>%
  add_trace(x = time_vec, y = stats_off$mean, type = 'scatter', mode = 'lines',
            name = 'OFF Mean', line = list(color = 'blue', width = 3)) %>%
  # Confidence bands
  add_trace(x = c(time_vec, rev(time_vec)), 
            y = c(stats_on$mean + 2*stats_on$sd, rev(stats_on$mean - 2*stats_on$sd)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(255,0,0,0.2)', name = 'ON ±2SD', showlegend = FALSE) %>%
  add_trace(x = c(time_vec, rev(time_vec)), 
            y = c(stats_off$mean + 2*stats_off$sd, rev(stats_off$mean - 2*stats_off$sd)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(0,0,255,0.2)', name = 'OFF ±2SD', showlegend = FALSE) %>%
  layout(
    title = "Average Time-Series with Confidence Bands",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "BOLD Signal"),
    hovermode = 'x unified'
  )

stats_plot

# Print summary statistics
cat("\n=== Time-Series Summary Statistics ===\n")
cat(sprintf("ON Condition - Global Mean: %.3f (SD: %.3f)\n", 
            stats_on$global_mean, stats_on$global_sd))
cat(sprintf("OFF Condition - Global Mean: %.3f (SD: %.3f)\n", 
            stats_off$global_mean, stats_off$global_sd))
cat(sprintf("Mean Difference (ON-OFF): %.3f\n", 
            stats_on$global_mean - stats_off$global_mean))
cat(sprintf("Average tSNR - ON: %.3f, OFF: %.3f\n", 
            mean(stats_on$tsnr), mean(stats_off$tsnr)))
```


# Implement KPT Core Functions

## Phase Distribution Functions

```{r phase-functions}
# Wrapped Laplace distribution density
wrapped_laplace_density <- function(theta, mu, beta) {
  # Normalize to [-pi, pi]
  theta <- ((theta - mu + pi) %% (2*pi)) - pi + mu
  
  # Wrapped Laplace density
  if(beta <= 0) return(rep(0, length(theta)))
  
  dens <- beta / (2 * sinh(beta)) * exp(beta * cos(theta - mu))
  return(dens)
}

# Sample from wrapped Laplace distribution
sample_wrapped_laplace <- function(n, mu, beta) {
  # Sample from linear Laplace
  u <- runif(n)
  linear_samples <- ifelse(u < 0.5,
                          mu + log(2*u)/beta,
                          mu - log(2*(1-u))/beta)
  
  # Wrap to [-pi, pi]
  wrapped_samples <- ((linear_samples + pi) %% (2*pi)) - pi
  return(wrapped_samples)
}

# Phase distribution parameters as functions of time
mu_t <- function(t, t_max = 150) {
  # Location parameter varies sinusoidally with time
  pi * sin(2 * pi * t / t_max)
}

beta_t <- function(t, t_max = 150) {
  # Scale parameter varies with time
  2 + sin(4 * pi * t / t_max)
}
```

## Harmonic Expansion and Moment Computation

```{r harmonic-functions}
# Compute empirical moments
compute_empirical_moments <- function(s_values, theta_values, n_max = 10) {
  n_samples <- length(s_values)
  moments <- complex(length = 2*n_max + 1)
  
  for(n in -n_max:n_max) {
    moments[n + n_max + 1] <- mean(s_values * exp(-1i * n * theta_values))
  }
  
  return(moments)
}

# Estimate harmonic coefficients f_n(t)
estimate_harmonic_coefficients <- function(time_series_collection, n_max = 10) {
  n_time_points <- dim(time_series_collection)[3]
  n_total_series <- dim(time_series_collection)[1] * dim(time_series_collection)[2]
  
  # Reshape to (n_total_series, n_time_points)
  reshaped_data <- matrix(time_series_collection, 
                         nrow = n_total_series, 
                         ncol = n_time_points)
  
  # Initialize coefficient matrix
  f_coeffs <- matrix(0, nrow = n_time_points, ncol = 2*n_max + 1)
  
  # For each time point, estimate mean signal
  for(t_idx in 1:n_time_points) {
    # Simple estimate: f_0(t) = mean signal, others = 0 for now
    f_coeffs[t_idx, n_max + 1] <- mean(reshaped_data[, t_idx])
    
    # Add some harmonic structure (for demonstration)
    for(n in 1:n_max) {
      f_coeffs[t_idx, n_max + 1 + n] <- 0.1 * f_coeffs[t_idx, n_max + 1] * 
        exp(-n/3) * cos(2*pi*n*t_idx/n_time_points)
      f_coeffs[t_idx, n_max + 1 - n] <- 0.1 * f_coeffs[t_idx, n_max + 1] * 
        exp(-n/3) * cos(2*pi*n*t_idx/n_time_points)
    }
  }
  
  return(f_coeffs)
}
```

## KPT Reconstruction Algorithm

```{r kpt-algorithm}
# KPT reconstruction function
kpt_reconstruct <- function(time_series_collection, n_max = 10, n_phase_samples = 100) {
  n_participants <- dim(time_series_collection)[1]
  n_runs <- dim(time_series_collection)[2]
  n_time_points <- dim(time_series_collection)[3]
  
  # Step 1: Estimate harmonic coefficients
  f_coeffs <- estimate_harmonic_coefficients(time_series_collection, n_max)
  
  # Step 2: For each time point, sample phases and compute moments
  reconstructed_phases <- list()
  empirical_moments_all <- matrix(0, nrow = n_time_points, ncol = 2*n_max + 1)
  
  for(t_idx in 1:n_time_points) {
    t_val <- t_idx * 2  # TR = 2
    
    # Get all measurements at this time point
    s_values <- as.vector(time_series_collection[, , t_idx])
    
    # Sample phases from wrapped Laplace
    mu <- mu_t(t_val)
    beta <- beta_t(t_val)
    theta_samples <- sample_wrapped_laplace(length(s_values), mu, beta)
    
    # Compute empirical moments
    moments <- compute_empirical_moments(s_values, theta_samples, n_max)
    empirical_moments_all[t_idx, ] <- moments
    
    # Store phase distribution parameters
    reconstructed_phases[[t_idx]] <- list(
      t = t_val,
      mu = mu,
      beta = beta,
      theta_samples = theta_samples,
      s_values = s_values
    )
  }
  
  return(list(
    f_coeffs = f_coeffs,
    empirical_moments = empirical_moments_all,
    phase_data = reconstructed_phases,
    n_max = n_max
  ))
}
```

## Kime-Surface Generation

```{r kimesurface-functions}
# Generate kime-surface from KPT reconstruction
generate_kime_surface <- function(kpt_result, n_grid = 50) {
  # Create polar grid
  r_max <- max(sapply(kpt_result$phase_data, function(x) x$t))
  r_grid <- seq(0, r_max, length.out = n_grid)
  theta_grid <- seq(-pi, pi, length.out = n_grid)
  
  # Initialize surface
  surface_z <- matrix(0, nrow = n_grid, ncol = n_grid)
  
  # Interpolate intensity values
  for(i in 1:n_grid) {
    for(j in 1:n_grid) {
      r <- r_grid[i]
      theta <- theta_grid[j]
      
      if(r > 0) {
        # Find nearest time point
        t_idx <- which.min(abs(sapply(kpt_result$phase_data, function(x) x$t) - r))
        
        # Compute expected intensity using harmonic expansion
        z_val <- Re(kpt_result$f_coeffs[t_idx, kpt_result$n_max + 1])  # f_0(t)
        
        # Add harmonic contributions
        for(n in 1:kpt_result$n_max) {
          z_val <- z_val + 2 * Re(kpt_result$f_coeffs[t_idx, kpt_result$n_max + 1 + n] * 
                                  exp(1i * n * theta))
        }
        
        surface_z[i, j] <- z_val
      }
    }
  }
  
  # Convert to Cartesian coordinates
  x_grid <- outer(r_grid, cos(theta_grid))
  y_grid <- outer(r_grid, sin(theta_grid))
  
  return(list(
    x = x_grid,
    y = y_grid,
    z = surface_z,
    r_grid = r_grid,
    theta_grid = theta_grid
  ))
}
```

# Apply KPT to fMRI Data

```{r apply-kpt}
# Apply KPT to ON and OFF conditions
cat("Applying KPT to ON condition...\n")
kpt_on <- kpt_reconstruct(fmri_data$ON, n_max = 8)

cat("Applying KPT to OFF condition...\n")
kpt_off <- kpt_reconstruct(fmri_data$OFF, n_max = 8)

# Generate kime-surfaces
surface_on <- generate_kime_surface(kpt_on, n_grid = 60)
surface_off <- generate_kime_surface(kpt_off, n_grid = 60)
```

# Visualize Results

## Interactive 3D Kime-Surfaces

```{r visualize-3d}
# Create interactive 3D plot for ON condition
fig_on <- plot_ly(
  x = surface_on$x,
  y = surface_on$y,
  z = surface_on$z,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  contours = list(
    z = list(
      show = TRUE,
      usecolormap = TRUE,
      highlightcolor = "#ff0000",
      project = list(z = TRUE)
    )
  ),
  lighting = list(
    ambient = 0.4,
    diffuse = 0.6,
    specular = 0.2,
    roughness = 0.85
  )
) %>%
  layout(
    title = "fMRI ON Condition - Kime-Surface",
    scene = list(
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z (Intensity)"),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    )
  )

# Create interactive 3D plot for OFF condition
fig_off <- plot_ly(
  x = surface_off$x,
  y = surface_off$y,
  z = surface_off$z,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  contours = list(
    z = list(
      show = TRUE,
      usecolormap = TRUE,
      highlightcolor = "#ff0000",
      project = list(z = TRUE)
    )
  ),
  lighting = list(
    ambient = 0.4,
    diffuse = 0.6,
    specular = 0.2,
    roughness = 0.85
  )
) %>%
  layout(
    title = "fMRI OFF Condition - Kime-Surface",
    scene = list(
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z (Intensity)"),
      camera = list(
        eye = list(x = 1.5, y = 1.5, z = 1.5)
      )
    )
  )

# Display the plots
fig_on
fig_off
```

## Comparison Plot

```{r comparison-plot}
# Create side-by-side comparison
fig_compare <- subplot(
  fig_on %>% layout(scene = list(domain = list(x = c(0, 0.45)))),
  fig_off %>% layout(scene = list(domain = list(x = c(0.55, 1)))),
  titleX = TRUE,
  titleY = TRUE
) %>%
  layout(
    title = "KPT Kime-Surfaces: fMRI ON vs OFF Conditions",
    showlegend = FALSE
  )

fig_compare

##### 3 synchronized fMRI On, Off, Difference Kime-surfaces
# ---------------------------- 
# Create Synchronized 3D Kime-Surfaces
# ---------------------------- 

# Prepare data for synchronized plots
prepare_surface_data <- function(surface) {
  return(list(
    x_grid = surface$x,
    y_grid = surface$y,
    z_grid = surface$z,
    magnitude = sqrt(surface$x^2 + surface$y^2),
    phase = atan2(surface$y, surface$x)
  ))
}

pOn <- prepare_surface_data(surface_on)
pOff <- prepare_surface_data(surface_off)

# Create ON surface plot
plot_on <- plot_ly(
  x = pOn$x_grid,
  y = pOn$y_grid,
  z = pOn$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  hoverinfo = "x+y+z",
  name = "ON",
  scene = "scene",
  showscale = FALSE
) %>%
  layout(scene = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create OFF surface plot
plot_off <- plot_ly(
  x = pOff$x_grid,
  y = pOff$y_grid,
  z = pOff$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  hoverinfo = "x+y+z",
  name = "OFF",
  scene = "scene2",
  showscale = FALSE
) %>%
  layout(scene2 = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create difference plot
diff_plot <- plot_ly(
  hoverinfo = "x+y+z",
  showscale = TRUE,
  scene = "scene3"
) %>%
  add_trace(
    x = pOn$x_grid,
    y = pOn$y_grid,
    z = pOn$z_grid - pOff$z_grid,
    surfacecolor = pOn$z_grid - pOff$z_grid,
    type = 'surface',
    colorscale = list(
      c(0, "rgb(0,0,255)"),
      c(0.5, "rgb(255,255,255)"),
      c(1, "rgb(255,0,0)")
    ),
    colorbar = list(
      title = "Difference",
      x = 1.1
    ),
    opacity = 1,
    visible = TRUE
  ) %>%
  layout(
    scene3 = list(
      aspectmode = "cube",
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z = ON - OFF")
    )
  )

# Combine into synchronized subplot
main_plot <- subplot(plot_on, plot_off, diff_plot) %>%
  layout(
    title = "Synchronized Kime-Surfaces (Top-Left: ON | Top-Right: OFF | Bottom: Difference)",
    scene = list(
      domain = list(x = c(0, 0.45), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene2 = list(
      domain = list(x = c(0.55, 1), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene3 = list(
      domain = list(x = c(0.275, 0.725), y = c(0, 0.45)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    showlegend = FALSE
  )

# Add camera synchronization
synchronized_plot <- 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);
        }
      });
    }"
  )

# Display the synchronized plot
synchronized_plot



# ---------------------------- 
# Create Synchronized 3D Kime-Surfaces
# ---------------------------- 

# Prepare data for synchronized plots with cylindrical coordinates
prepare_surface_data <- function(surface) {
  # Calculate cylindrical coordinates
  r_matrix <- sqrt(surface$x^2 + surface$y^2)
  phi_matrix <- atan2(surface$y, surface$x)
  
  return(list(
    x_grid = surface$x,
    y_grid = surface$y,
    z_grid = surface$z,
    r_grid = r_matrix,
    phi_grid = phi_matrix,
    magnitude = r_matrix,
    phase = phi_matrix
  ))
}

pOn <- prepare_surface_data(surface_on)
pOff <- prepare_surface_data(surface_off)

# Create custom hover text for cylindrical coordinates
create_hover_text <- function(r, phi, z, label = "") {
  hover_matrix <- matrix(NA, nrow = nrow(r), ncol = ncol(r))
  for(i in 1:nrow(r)) {
    for(j in 1:ncol(r)) {
      hover_matrix[i,j] <- sprintf(
        "%sTime: %.1f<br>φ: %.3f<br>Intensity: %.3f",
        ifelse(label == "", "", paste0(label, "<br>")),
        r[i,j],
        phi[i,j],
        z[i,j]
      )
    }
  }
  return(hover_matrix)
}

# Create ON surface plot with cylindrical coordinates
plot_on <- plot_ly(
  x = pOn$x_grid,
  y = pOn$y_grid,
  z = pOn$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(20,20,80)"),
    c(0.25, "rgb(20,80,140)"),
    c(0.5, "rgb(20,140,200)"),
    c(0.75, "rgb(80,200,120)"),
    c(1, "rgb(255,220,80)")
  ),
  hovertext = create_hover_text(pOn$r_grid, pOn$phi_grid, pOn$z_grid, "ON Condition"),
  hoverinfo = "text",
  name = "ON",
  scene = "scene",
  showscale = FALSE
) %>%
  layout(scene = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create OFF surface plot with cylindrical coordinates
plot_off <- plot_ly(
  x = pOff$x_grid,
  y = pOff$y_grid,
  z = pOff$z_grid,
  type = "surface",
  colorscale = list(
    c(0, "rgb(40,20,60)"),
    c(0.25, "rgb(60,40,100)"),
    c(0.5, "rgb(80,60,140)"),
    c(0.75, "rgb(100,80,180)"),
    c(1, "rgb(140,120,220)")
  ),
  hovertext = create_hover_text(pOff$r_grid, pOff$phi_grid, pOff$z_grid, "OFF Condition"),
  hoverinfo = "text",
  name = "OFF",
  scene = "scene2",
  showscale = FALSE
) %>%
  layout(scene2 = list(
    xaxis = list(title = "X (t·cos θ)"),
    yaxis = list(title = "Y (t·sin θ)"),
    zaxis = list(title = "Intensity"),
    aspectmode = "cube"
  ))

# Create difference plot with cylindrical coordinates
diff_z <- pOn$z_grid - pOff$z_grid
diff_plot <- plot_ly(
  hovertext = create_hover_text(pOn$r_grid, pOn$phi_grid, diff_z, "Difference (ON - OFF)"),
  hoverinfo = "text",
  showscale = TRUE,
  scene = "scene3"
) %>%
  add_trace(
    x = pOn$x_grid,
    y = pOn$y_grid,
    z = diff_z,
    surfacecolor = diff_z,
    type = 'surface',
    colorscale = list(
      c(0, "rgb(0,0,255)"),
      c(0.5, "rgb(255,255,255)"),
      c(1, "rgb(255,0,0)")
    ),
    colorbar = list(
      title = "Difference",
      x = 1.1
    ),
    opacity = 1,
    visible = TRUE
  ) %>%
  layout(
    scene3 = list(
      aspectmode = "cube",
      xaxis = list(title = "X (t·cos θ)"),
      yaxis = list(title = "Y (t·sin θ)"),
      zaxis = list(title = "Z = ON - OFF")
    )
  )

# Combine into synchronized subplot
main_plot <- subplot(plot_on, plot_off, diff_plot) %>%
  layout(
    title = "Synchronized Kime-Surfaces (Top-Left: ON | Top-Right: OFF | Bottom: Difference)",
    scene = list(
      domain = list(x = c(0, 0.45), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene2 = list(
      domain = list(x = c(0.55, 1), y = c(0.5, 1)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    scene3 = list(
      domain = list(x = c(0.275, 0.725), y = c(0, 0.45)),
      aspectmode = "cube",
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5))
    ),
    showlegend = FALSE
  )

# Add camera synchronization
synchronized_plot <- 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);
        }
      });
    }"
  )

# Display the synchronized plot
synchronized_plot

```

## Phase Distribution Analysis

```{r phase-analysis}
# Analyze phase distributions at different time points
time_points <- c(25, 75, 125)
phase_analysis <- data.frame()

for(condition in c("ON", "OFF")) {
  kpt_result <- if(condition == "ON") kpt_on else kpt_off
  
  for(t_idx in time_points) {
    phase_data <- kpt_result$phase_data[[t_idx/2]]  # Adjust for TR=2
    
    df_temp <- data.frame(
      condition = condition,
      time = t_idx,
      theta = phase_data$theta_samples,
      intensity = phase_data$s_values
    )
    phase_analysis <- rbind(phase_analysis, df_temp)
  }
}

# Plot phase distributions
ggplot(phase_analysis, aes(x = theta, fill = condition)) +
  geom_histogram(bins = 30, alpha = 0.6, position = "identity") +
  facet_wrap(~time, ncol = 3, labeller = labeller(time = function(x) paste("t =", x, "s"))) +
  scale_fill_manual(values = c("ON" = "red", "OFF" = "blue")) +
  xlim(-pi, pi) +
  labs(
    title = "Phase Distributions at Different Time Points",
    x = "Phase (θ)",
    y = "Count",
    fill = "Condition"
  ) +
  theme_minimal()
```

## Radial Intensity Profiles

```{r radial-profiles}
# Extract radial intensity profiles
radial_profiles <- data.frame()

for(condition in c("ON", "OFF")) {
  surface <- if(condition == "ON") surface_on else surface_off
  
  # Average intensity at each radius
  for(i in 1:length(surface$r_grid)) {
    avg_intensity <- mean(surface$z[i, ], na.rm = TRUE)
    
    radial_profiles <- rbind(radial_profiles, data.frame(
      condition = condition,
      radius = surface$r_grid[i],
      intensity = avg_intensity
    ))
  }
}

# Plot radial profiles
ggplot(radial_profiles, aes(x = radius, y = intensity, color = condition)) +
  geom_line(size = 1.5) +
  scale_color_manual(values = c("ON" = "red", "OFF" = "blue")) +
  labs(
    title = "Radial Intensity Profiles",
    x = "Time (t)",
    y = "Average Intensity",
    color = "Condition"
  ) +
  theme_minimal()
```

# Statistical Analysis

```{r statistics}
# Compute surface metrics
compute_surface_metrics <- function(surface) {
  # Total volume under surface
  volume <- sum(surface$z, na.rm = TRUE) * 
    (surface$r_grid[2] - surface$r_grid[1]) * 
    (surface$theta_grid[2] - surface$theta_grid[1])
  
  # Surface roughness (standard deviation of gradients)
  grad_r <- diff(surface$z) / (surface$r_grid[2] - surface$r_grid[1])
  grad_theta <- t(diff(t(surface$z))) / (surface$theta_grid[2] - surface$theta_grid[1])
  roughness <- sd(c(grad_r, grad_theta), na.rm = TRUE)
  
  # Peak intensity
  peak <- max(surface$z, na.rm = TRUE)
  
  # Radial asymmetry (coefficient of variation across angles)
  radial_cv <- apply(surface$z, 1, function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE))
  asymmetry <- mean(radial_cv, na.rm = TRUE)
  
  return(list(
    volume = volume,
    roughness = roughness,
    peak = peak,
    asymmetry = asymmetry
  ))
}

# Compare metrics
metrics_on <- compute_surface_metrics(surface_on)
metrics_off <- compute_surface_metrics(surface_off)

# Create comparison table
metrics_df <- data.frame(
  Metric = c("Volume", "Roughness", "Peak Intensity", "Asymmetry"),
  ON = unlist(metrics_on),
  OFF = unlist(metrics_off),
  Difference = unlist(metrics_on) - unlist(metrics_off),
  Ratio = unlist(metrics_on) / unlist(metrics_off)
)

knitr::kable(metrics_df, digits = 3, caption = "Kime-Surface Metrics Comparison")
```

# Advanced Statistical Analysis

```{r advanced-statistics}
# # ---------------------------- 
# Enhanced Difference Analysis
# ---------------------------- 

# Compute point-wise differences
diff_surface <- list(
  x = surface_on$x,
  y = surface_on$y,
  z = surface_on$z - surface_off$z,
  r_grid = surface_on$r_grid,
  theta_grid = surface_on$theta_grid
)

# Statistical significance of differences
diff_stats <- data.frame(
  radius = surface_on$r_grid,
  mean_diff = apply(diff_surface$z, 1, mean, na.rm = TRUE),
  sd_diff = apply(diff_surface$z, 1, sd, na.rm = TRUE),
  max_diff = apply(diff_surface$z, 1, max, na.rm = TRUE),
  min_diff = apply(diff_surface$z, 1, min, na.rm = TRUE)
)

# Plot difference statistics
diff_stats_plot <- plot_ly() %>%
  add_trace(x = diff_stats$radius, y = diff_stats$mean_diff,
            type = 'scatter', mode = 'lines', name = 'Mean Difference',
            line = list(color = 'black', width = 3)) %>%
  add_trace(x = c(diff_stats$radius, rev(diff_stats$radius)),
            y = c(diff_stats$mean_diff + diff_stats$sd_diff,
                  rev(diff_stats$mean_diff - diff_stats$sd_diff)),
            type = 'scatter', mode = 'none', fill = 'toself',
            fillcolor = 'rgba(100,100,100,0.3)', name = '±1 SD',
            showlegend = TRUE) %>%
  add_trace(x = diff_stats$radius, y = diff_stats$max_diff,
            type = 'scatter', mode = 'lines', name = 'Max Difference',
            line = list(color = 'red', dash = 'dash')) %>%
  add_trace(x = diff_stats$radius, y = diff_stats$min_diff,
            type = 'scatter', mode = 'lines', name = 'Min Difference',
            line = list(color = 'blue', dash = 'dash')) %>%
  layout(
    title = "Kime-Surface Differences (ON - OFF) vs Time",
    xaxis = list(title = "Time (seconds)"),
    yaxis = list(title = "Intensity Difference"),
    hovermode = 'x unified'
  )

diff_stats_plot
```

This implementation of the [Kime-Phase Tomography (KPT) algorithm](https://www.socr.umich.edu/TCIU/HTMLs/TCIU_SK_Appendix03_KPT.html) 
demonstrates the Kime-Phase Tomography applied to synthetic fMRI data.

 1. *Successful Reconstruction*: The algorithm successfully maps repeated time-series measurements to 2D kime-surfaces, revealing the underlying temporal-phase structure.

 2. *Condition Differences*: Clear differences are observed between ON and OFF conditions:
   - ON condition shows higher intensity and more complex surface structure
   - OFF condition exhibits smoother, lower-intensity surfaces
   - Phase distributions vary systematically with time

 3. *Key Insights*:
   - The kime-surface representation captures both temporal dynamics (radial direction) and phase variability (angular direction)
   - The harmonic expansion provides a compact representation of the signal structure
   - The wrapped Laplace distribution effectively models phase uncertainty

 4. *Future Directions*:
   - Apply to real fMRI data
   - Explore different phase distributions (von Mises, wrapped Cauchy)
   - Investigate higher-order harmonic expansions
   - Develop statistical tests for comparing kime-surfaces

KPT provides a geometric QM-perspective on longitudinal time-series data
and accounts for the structure of repeated measurements and their inherent variability.


<!--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-->
