SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
This TCIU Section in Chapter 6 extends previous work to show alternative strategies to transform repeated data measurements of longitudinal data as 2D parametric manifolds (kimesurfaces).
There are alternative strategies to map series of repeated longitudinal measurements into kimesurfaces, i.e., 2D manifolds parameterized by complex time \(\kappa = t e^{i\theta}\). Types of such transformations include analytical, numerical, and algorithmic approaches. The selection of an appropriate transform of series of repeated longitudinal measurements into kimesurfaces generally depends on the specific data characteristics including data density, noise levels, and interpretability needs. The resulting kimesurface can effectively unify repeated time-series into a structured 2D manifold, enabling novel subsequent AI modeling and analytics to support gaining intuition and insights into the phenomenon’s temporal-phase variability.
There are multiple complementary approaches for mapping repeated time-series (longitudinal data) into a 2D kimesurface, i.e., a manifold parametrized by the kime magnitude and phase \((t,\theta)\), where \(t = |\kappa|\), the usual notion of time, and \(\theta\) is an added phase index indexing the repeated samples.
Analytical transforms of time-series to kimesurfaces provide interpretability but require mathematical tractability.
One idea is to mathematically transform the time-series data \(f(t)\) into a complex function \(F(\kappa)\) with \(\kappa = t e^{i\theta}\). For instance, consider the generalized Laplace transform, \(\mathcal{L}:\{f\in L^1(\mathbb{R}^+)\} \longrightarrow \{F: \mathbb{C}\overbrace{\longrightarrow}^{analytic}\mathbb{C}\}\), which is defined over \(\kappa = t e^{i\theta}\in \mathbb{C}\) by
\[F(\kappa) = \mathcal{L}\{f(t)\}(\kappa) = \int_0^\infty f(t) e^{-\kappa t} dt = \int_0^\infty f(t) e^{-t^2 e^{i\theta}} dt,\]
where \(\theta\) as a parameter. This maps time-series data to a complex surface, with \(\theta\) indexing the repeated longitudinal trials. The output of the Laplace transform is a complex analytic function with suitable decay, which has a convergent power series expansion whose coefficients decompose the function \(F\) into its (statistical) moments. Recall that the standard Laplace transform (LT) of \(f(t)\) is \(\mathcal{L}\{f\}(s) = \int_0^\infty f(t)\,e^{-st}dt\), for \(\mathrm{Re}(s) > 0\). The LT generalizes to
\[F(\kappa)\;=\;\int_0^\infty f(t)\,
e^{-\kappa\,t}\,dt
\quad\text{or}\quad F(\kappa)\;=\;\int_0^\infty f(\tau)\,\kappa^{-
\tau}\,d\tau,\] where \(\kappa =
t\,e^{i\theta}\) might be treated as the transform
variable.
However, this requires ensuring the integral convergence and
appropriate interpretation of how \(\kappa\mapsto t e^{i\theta}\) factors
in.
Alternatively, we can employ the Fourier series expansion for each fixed \(t\), expand observations across \(\theta\) into the Fourier coefficients, \(c_k(t)\), \(f(t, \theta) = \sum_{k=-\infty}^\infty c_k(t) e^{ik\theta},\) where \(c_k(t)\) are time-dependent coefficients estimated from repeated trials.
In practice, having a repeated dataset \(\{f_n(t)\}\), requires transforming each one \(F_n(\kappa) = \int_0^T f_n(t)\,\varphi(\kappa,t)\,dt\) using some transform kernel \(\varphi(\kappa,t)\). Then, if the phase \(\theta\sim\Phi(t)\) is the experiment index, we gather all \(\{F_n(\kappa)\}_{n-1}^N\) and interpret them as sampling phase angles in the complex plane. We can also consider a Fourier transform in time, but then replace the real frequency \(\omega\) by the phase angle \(\theta\).
Another alternative is to consider a Hilbert transform approach to define an analytic continuation of the real time series, thereby obtaining a phase function for each sample. In many of these transforms, a purely real time variable is vital as input. Extending time to kime, \(t\to \kappa = t\, e^{i\theta}\) can complicate standard theorems on convergence, hence we need to validate the rigorous integral transforms into a complex domain.
Many examples of Laplace-transformed time-series are shown in the TCIU Chapter 6 Tutorial.
First plot the simulated On (stimulus) and Off (rest) fMRI time-series at a fixed voxel location \((44,42,33)\in\mathbb{R}^3\), along with the averaged (pooled) On and Off signal over all \(8\) repeats in the single run (epoch) of \(160\) time-points.
# Extract just the On kime-series at voxel (44,42,33), each time-series has 8 repeats!
library(TCIU)
library("R.matlab")
library(plotly)
# Data import
# 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
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.
# create the LT
NuLT = function(datax, datay, inputz, k = 3, fitwarning = FALSE,
mirror = FALSE, range = 2*pi) {
datax = as.numeric(datax)
datay = as.numeric(datay)
n = length(datax)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*range
if(mirror){
x1 = c(x1,rev(2*range - x1))/2
n = 2*n
datay = c(datay, rev(datay))
plot(x1, datay)
}
#generate the coefficients in indefinite integral of t^n*exp(-zt)
coef = 1;
coefm = as.matrix(coef)
for(i in 1:k){
coefm = cbind(coefm,0)
coef = c(coef*i,1)
coefm = rbind(coefm,coef)
}
# these coefficients ordered by ^0, ^1, ^2, ... in column format
# compute 1, z, z^2...,z^k
zz = cbind(1,inputz)
zt = inputz
for (i in 2:k){
zt = zt*inputz
zz = cbind(zz,zt)
}
zd = zt*inputz
# compute 1, x, x^2...,x^k
tx = x1;
xm = cbind(1,x1)
for (i in 2:k){
tx = tx*x1
xm = cbind(xm,tx)
}
# sum over intervals
result = 0*inputz
ii = 1
while(ii+k <= n) {
A = xm[seq(ii,ii+k),c(0:k+1)]
b = datay[seq(ii,ii+k)]
# polyfit might be faster when using polynomial basis, while matrix inverse, `solve()`,
# is the more general approach that works for any function basis
polyc = as.numeric(solve(A,b))
#ordered by ^0, ^1, ^2, ... in column format
# Enter a new function variable qualityCheck=FALSE
# check fit quality; this step can be skipped for speed/efficiency
# if (qualityCheck) { .... }
if (fitwarning){
xx = seq(A[1,2],A[k+1,2],length.out = 100);
yy = polyval(rev(polyc),xx)
if(max(abs(yy-mean(b)))>2*max(abs(b-mean(b)))){
print(c("Warning: Poor Fit at ",ii,", Largest Deviation is",
max(abs(yy-mean(b)))))
print(c("Spline Polynomial is", polyc),3);
#print(c(polyval(rev(polyc),A[,2]),b))
plot(xx, yy, main="Polynomial fit", ylab="", type="l", col="blue")
lines(A[,2],b, col="red")
legend("topleft",c("fit","data"),fill=c("blue","red"))
print(" ")
}
}
# Use vector/matrix operations to avoid looping,
# some of the expressions look weird
# May need to actually compare the efficiency/speed of
# vector based vs. standard numeric calculations
m1 = t(t(polyc*coefm)*A[1,])
m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1]
m2 = t(t(polyc*coefm)*A[k+1,])
m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
exp(-inputz*A[k+1,2])*colSums(t(zz)*m22))/zd
result = result+intgl
ii=ii+k
}
# Computations over the last interval
if(ii < n){
nk = n-ii;
A = xm[seq(ii,ii+nk),c(0:nk+1)]
b = datay[seq(ii,ii+nk)]
nc = as.numeric(solve(A,b))
nc = c(nc,seq(0,0,length.out = k-nk))
A = xm[seq(ii,ii+nk),]
m1 = t(t(nc*coefm)*A[1,])
m11 = as.numeric(tapply(m1, col(m1)-row(m1), sum))[0:k+1]
m2 = t(t(nc*coefm)*A[nk+1,])
m22 = as.numeric(tapply(m2, col(m2)-row(m2), sum))[0:k+1]
# cc = colSums(coefm*polyc)
intgl = (exp(-inputz*A[1,2])*colSums(t(zz)*m11)-
exp(-inputz*A[nk+1,2])*colSums(t(zz)*m22))/zd
result = result+intgl
}
#offset = 0.05*pi
#result = result+datay[n]*(exp(-2*pi*inputz)-exp(-(2*pi+offset)*inputz))/inputz
return(result)
}
Let’s test the discrete LT using the \(\sin(x),\cos(x)\) function.
datax = seq(1,200)
n = length(datax)
x1 = n/(n+0.5)*((datax-min(datax))/(max(datax)-min(datax)))*(2*pi)
datay = cos(x1) # test function!!!
Lout = 61
range_limit <- 20
x2 <- seq(from = 0.1, to = range_limit, length.out = Lout)
y2 <- seq(from = 0, to = range_limit, length.out = Lout)
x2_g = x2 %o% seq(1,1, length.out = Lout)
y2_g = seq(1,1, length.out = Lout)%o%y2
z2_grid = x2 %o% y2
argz = as.vector(x2_g + 1i*y2_g)
LTz = NuLT(x1, datay, argz)
rec1 = matrix(LTz,nrow = Lout)
recm = abs(rec1)
recr = Re(rec1)
reci = Im(rec1)
ph = reci/recm
surf_color <- atan2(reci,recr)
# colorscale = cbind(seq(-pi, pi, by=1/10, rainbow(11)))
p <- plot_ly(hoverinfo="none", showscale = TRUE) %>%
add_trace(x = x2_g, y = y2_g, z = surf_color,
# F-magnitude or Re(F), # z = Im(z2_grid), # Real or Imaginary part of F(t)
surfacecolor=log(recm), # colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title =
"Truncated Laplace transform of cos(x), Color = Mag(Z)",
showlegend = TRUE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0),
xaxis=list(title= "Real=Re(z)"), yaxis=list(title = "Imaginary=Im(z)"),
zaxis=list(title= "Z = tan-1(Im(L[f](z))/Re(L[f](z))))")) ) # 1:1:1 aspect ratio
p