| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
Spacekime Analytics (Time Complexity and Inferential Uncertainty)
TCIU: The Laplace Transform of Longitudinal Data: Analytic Duality of Time-series and Kimesurfaces
SOCR Team
07/03/2024
The Laplace Transform, \(\mathcal{L}\), allows us to examine the relations between the space-time and space-kime representations of longitudinal data. The Fourier transformation is a linear operator that maps complex-valued functions of real variables (e.g., space, time domains) to complex valued functions of other real variables (e.g., frequency domain). The Laplace transform is similar, however, it sends complex-valued functions of positive real variables (e.g., time) to complex-valued functions defined on complex variables (e.g., kime).
1 Continuous Forward and Inverse Laplace Transforms (LT)
The forward and inverse (continuous) Laplace transforms are defined below.
- For a given function function (of time) \(f(t): \mathbb{R}^+ \longrightarrow \mathbb{C}\), the Laplace transform is the function of a complex frequency argument, \(F(z)={\mathcal {L}}(f)(z):\mathbb{C} \longrightarrow \mathbb{C}\):
\[{\mathcal {L}}(f)(z)=F(z)=\int_{0}^{\infty} {f(t)e^{-z t} dt}.\]
- For a given function function of a complex frequency argument, \(F(z)\), the Inverse Laplace transform (ILT) is the function of a positive real (time-like) argument \(f(t)={\mathcal {L}}^{-1}(F)(t):\mathbb{R}^+ \longrightarrow \mathbb{C}\), which is defined in terms of a complex path integral (a.k.a. Bromwich integral or Fourier–Mellin integral):
\[f(t)={\mathcal {L}}^{-1}(F)(t)={\frac {1}{2\pi i}} \lim_{T\to \infty } \int_{\gamma -iT}^{\gamma +iT} { e^{z t}F(z) dz},\]
where the parameter \(\gamma\in \mathbb{R}\) is chosen so that the entire complex contour path of the integral is inside of the region of convergence of \(F(z)\).
The Laplace transform plays an interesting role as an expected value in the field of probability and statistics. If \(X\) is a random variable, then its Laplace transform, i.e., the LT of its probability density function \(f_X\), is given by the expectation of an exponential:
\[{\mathcal {L}}(X)={\mathcal {L}}(f)(z)= E\left ( e^{-zX}\right ).\]
Note that when \(z=−t\in R\), then the \({\mathcal {L}}(X)=E\left ( e^{tX}\right )\) is just the moment generating function of \(X\).
Another application of LT involves recovering the cumulative distribution function of a continuous random variable \(X\), \(F_{X}(x)=\int_{y<x}{f_X(y) dy}\):
\[ F_{X}(x)={\mathcal {L}}^{-1} \left \{ {\frac {1}{z}} E \left[e^{-zX}\right]\right \}\!(x)= {\mathcal {L}}^{-1} \left \{ {\frac {1}{z}} {\mathcal {L}}(f)(z) \right \}\!(x).\]
2 Discrete Forward Laplace Transforms
The discrete Forward and Inverse Laplace Transforms are
computed by numerically estimating the corresponding continuous
transformation integrals. The numerical approximation of the LT is
defined by the LT() method below.
Now we can test and compare the discrete and continuous LT (\(\mathcal{L}\)) on a pair of functions that have closed-form analytical expressions.
\[f_1(t)=\mathcal{L}^{-1}(F_1)(t)=t\ \underset{ILT}{\overset{LT}{\longleftrightarrow}} \ F_1(z)=\mathcal{L}(f_1)(z)=\frac{1}{z^2}.\] And \[ f_2(t)=\mathcal{L}^{-1}(F_2)(t)=e^{-5t} \ \underset{ILT}{\overset{LT}{\longleftrightarrow}}\ F_2(z)=\mathcal{L}(f_2)(z)=\frac{1}{z+5}.\]
### 2. Tests the LT
z= 1+1i # Complex-domain value
f <- function(t) { t } # test function, F(z)=L(f)(z)=1/z^2
LT(f, z); 1/z^2## [1] 0-0.5i## [1] 0-0.5i## [1] 0.1621622-0.027027i## [1] 0.1621622-0.027027i3 Discrete Inverse Laplace Transform
3.1 Basic functions
Let’s start by defining some basic functions that transform between Cartesian and Polar coordinates.
################ Cartesian to Polar Coordinate Transform ##############################################
#' Cartesian to Polar Coordinate Transform
#' @title Cartesian to Polar Coordinate Transform
#' @description In r.xy Returns polar coordinate r from a pair of Cartesian coordinates (x, y)
#' @param x x co-ordinate
#' @param y y co-ordinate
#' @return r or phi respectively from x and y
r.xy <- function(x, y){
return(sqrt(x^2 + y^2))
}
#' @description In r.xy Returns polar coordinate phi from a pair of Cartesian coordinates (x, y).
phi.xy <- function(x, y){
return(atan2(y, x))
}
################ Polar to Cartesian Coordinate Transform ##############################################
#' Polar to Cartesian Coordinate Transform
#' @title Polar to Cartesian Coordinate Transform
#' @description In x.rphi Returns cartesian coordinate x from a pair of polar coordinates (r, phi)
#'
#' @param r distance from origin (Time)
#' @param phi phase
#' @return x or y respectively from r and phi
x.rphi <- function(r, phi){
return(r*cos(phi))
}
#' @description In y.rphi turns cartesian coordinate y from a pair of polar coordinates (r, phi)
y.rphi <- function(r, phi){
return(r*sin(phi))
}3.2 Contour paths in \(C\)
Next we will define alternative contour paths in the complex plane,
\(C\), the optimal contour
optimContour(), the Bromwich contour, and their
corresponding numerical derivatives along contour paths,
analyticalPathDerivative() and
bromwichPathDerivative() functions.
################## Optimum contour path in Complex plane (needed for ILT) ##########################
# Uses Evans, Chung (2000) method
#'
#' @title Optimum Contour Path
#' @description The optimum contour path in polar coordinates (r,phi)
#'
#' @param phi phi (polar phase) value
#' @param m contour width; small values may lead to singularities on negative x-axis;
#' large values may lead to instabilities on the positive x-axis
#' @param t time (R^+ domain) variable, interdependent with the contour width, m
#'
#' @references Evans, Chung, 2000: Laplace transform inversions using optimal contours
#' in the complex plane, International Journal of Computer Mathematics, v73, pp 531-543.
#'
optimContour <- function(phi, m=1, t=5){
if(identical(t, 0)){ t <- 5 } # reset t=0 to avoid singularities
return(m*phi/(t*sin(phi)))
}
#' Numerical derivative of contour path length with respect to phi
#' return(((m*phi/(t*sin(phi)))^2 + (1/sin(phi) - phi/(tan(phi)*sin(phi)))^2)^0.5)
#' s as magnitude: real
#' s as dx and dy components: complex
#' @param phi phi (polar phase) value
#' @param m contour width; small values may lead to singularities on negative x-axis;
#' large values may lead to instabilities on the positive x-axis
#' @param t time value
analyticalPathDerivative <- function(phi, m=1, t=5) {
if(identical(t, 0)) { t <- 100 } # reset t = 0 to avoid singularities
if(identical(phi, 0)) { dx_dphi <- 0 } # avoid small phi value singularities
else{ dx_dphi <- (m/t)*((sin(phi) - phi*cos(phi))*cos(phi)/(sin(phi))^2 - phi) }
dy_dphi <- (m/t)
return( complex(real=dx_dphi, imaginary = dy_dphi) )
}
#' @title Bromwich contour path
#' @description Bromwich contour path - a straight vertical line through x=gamma
#' @param phi phi (polar phase) value
#' @param gamma value on the positive x-axis for the vertical line representing the contour
#' @param t time value
bromwichContour = function(phi, gamma = 0.5) {
return(gamma/cos(phi))
}
#' Numerical derivative of Bromwich contour
#' @param phi phi (polar phase) value
#' @param gamma value on the positive x-axis for the vertical line representing the contour
bromwichPathDerivative=function (phi, gamma = 0.5) {
return((0+1i) * (bromwichContour(phi, gamma)^2 + (gamma * tan(phi)/cos(phi))^2)^0.5)
}3.3 Define the discrete Inverse Laplace Transform (ILT)
Next, we will define the discrete ILT() function, which
inverts the Laplace Transform, \(\mathcal{L}^{-1}\).
Let’s test the ILT() function using some simple
functions. For instance, we can use the same 2 functions we defined
above:
\[f_1(t)=\mathcal{L}^{-1}(F_1)(t)=t\ \underset{ILT}{\overset{LT}{\longleftrightarrow}} \ F_1(z)=\mathcal{L}(f_1)(z)=\frac{1}{z^2},\] and \[ f_2(t)=\mathcal{L}^{-1}(F_2)(t)=e^{-5t} \ \underset{ILT}{\overset{LT}{\longleftrightarrow}}\ F_2(z)=\mathcal{L}(f_2)(z)=\frac{1}{z+5}.\]
# Tests the ILT
t = 1/5 # R-domain value
f <- function(z) { 1/(z^2) } # test function, F(z)=L(f)(z)=1/z^2 ==> f(t)=t
ILT(f, t); t## [1] 0.2+0i## [1] 0.2f <- function(z) { 1/(z+5) } # test function, F(z)=L(f)(z)=1/(z+5) ==> f(t)=e^-{5t}
ILT(f, t); exp(-5*t)## [1] 0.3678794+0i## [1] 0.3678794For longitudinal data, the Laplace transform provides one explicit mapping of the duality between functions of time (time-series) and functions of kime (kimesurfaces). The next sections illustrate that dichotomy.
4 Regularizing Kimesurfaces by Stereographic Projection of the Laplace Transform
As infinities and singularities may introduce substantive analytical challenges associated with the spacekime analytics based on the Laplace transformed time-series, various regularization schemes can be introduced to annihilate such singularities over the complex domain. One example of a Laplace transform regularization utilizes stereographic projection, which removes the singularities \(F(z)\) by compactifying the range space.
The surface of the unit sphere \(S^2\subseteq \mathbb{R}^3\) is the set comprised of all 3D Cartesian coordinates \({\bf{x}}=(x, y, z)\in\mathbb{R}^3\) satisfying the equation of a sphere \(||{\bf{x}}||^2=x^2 + y^2 + z^2 = 1\). Set the north pole of the sphere to be \(N = (0, 0, 1)\) and \(M=S^2 \setminus \{N\}\) represent all points on the sphere other than the north pole \(N\). The canonical axial (transferse) plane \(A=\{(x, y, z)\in\mathbb{R}^3 \ |\ z = 0\}\) runs through the center of the sphere and splits it into 2 hemispheres with the equator \(E=S^2\cap A\equiv \{(x, y, z)\in\mathbb{R}^3 \ |\ z = 0\ \& \ x^2+y^2=1\}\).
For each point on the sphere other than the north pole, \(\forall P\in M\), \(\exists! \ \vec{NP}\), a line through \(N\) and \(P\), such that \(\vec{NP} \cap A=P'\) is the stereographic projection of \(P\) onto the axial plane \(A\).
For each point on the sphere, \((x, y, z)\in M\), the bijective stereographic maps \(\mathcal{S}:M\to A\), \(\mathcal{S}(x,y,z)=(X,Y)\), and \(\mathcal{S}^{-1}:A\to M\), \(\mathcal{S}^{-1}(X,Y)=(x,y,z)\), such \(\mathcal{S}^{-1}\circ \underbrace{\mathcal{S}(x,y,z)}_{(X,Y)} = (x,y,z)\) and \(\mathcal{S}\circ \underbrace{\mathcal{S}^{-1}(X,Y)}_{(x,y,z)} = (X,Y)\) are inverse of each other:
\[\mathcal{S}(x,y,z)\equiv \left({\frac {x}{1-z}},{\frac {y}{1-z}}\right)=(X,Y)\in A, \ \forall (x,y,z)\in M,\ z\not=1,\] \[\mathcal{S}^{-1}(X,Y)\equiv \left({\frac {2X}{1+X^{2}+Y^{2}}},{\frac {2Y}{1+X^{2}+Y^{2}}}, {\frac {-1+X^{2}+Y^{2}}{1+X^{2}+Y^{2}}}\right)=(x,y,z)\in M.\]
Remember that we are trying to regularize the range \(\mathbb{C}\) of the Laplace transform values, \(F(z)\in\mathbb{C}\), not its kime domain, \(z\in\mathbb{C}\).
We can also consider the 3D spherical coordinates parameterizing the sphere by the zenith angle, \(0\leq \phi \leq \pi\) and the azimuth phase, \(0 \leq \theta\leq 2\pi\), \(M=\{(r\equiv 1, \phi, \theta)\ |\ 0\leq \phi \leq \pi, 0 \leq \theta\leq 2\pi\}\).
And similarly parameterize the 2D axial projection plane using polar coordinates \(A=\{(R, \Theta)\ |\ R\geq 0, 0 \leq \Theta\leq 2\pi\}\). Then the bijective stereographic maps are
\[\mathcal{S}(\phi,\theta)\equiv \left({\frac {\sin \phi }{1-\cos \phi }},\theta \right)=\left(\cot {\frac {\phi }{2}},\theta \right)=(R,\Theta )\in A,\] \[\mathcal{S}^{-1}(R,\Theta)\equiv \left(2\arctan \left (\frac {1}{R}\right ), \Theta \right)= (\phi ,\theta ) \in M.\]
The wellposedness of these maps requires specifying that \(\phi=\pi \iff R = 0\). Let’s derive the analytical expression of the inverse stereographic map in polar coordinates.
The 2D schematic below shows the projection triangle consisting of \(A\) is the center of the sphere, the points \(A\) and \(B\) are on the axial projection complex plane, and the pair of points \(N\) (North Pole) and \(D\) that are on the sphere.
circle_seg <- function(from, to,scaling){
by <- (to-from)/4
t <- seq(from*pi/180, to*pi/180 , by*pi/180)
x0 <- 0
y0 <- 0
r <- 0.1*scaling
x <- x0 + r*cos(t)
y<- y0 + r*sin(t)
lists <- list(x, y)
return (lists)
}# set up the circle segment
#return five point to interpolate the circle segment
#from: starting phase, to: ending phase
scaling <- 3
lista <- circle_seg(-90,37,3)
listb <- circle_seg(0,37,3)
listc <- circle_seg(90,37,3)
listd <- circle_seg(-90,-53,3)
liste <- circle_seg(0,-53,3)
listf <- circle_seg(90,-53,3)
width=780
height=400
library(plotly)
x <- c(0,2,0.8,0.4,0)
y <- c(0,0,0.6,0.8,1)
xp <- c(0,0,0.3,1/3,0.6)
yp <- c(0,1,0.1,0,-0.8)
data <- data.frame(x, y)
t <- list(family = "Arial", size = 14, color = toRGB("grey50"))
fig1 <- plot_ly(data, x = ~x, y = ~y, type = 'scatter', mode = 'lines',
text = c("A","B","D","E","N"),width=width,height=height)
fig1 <- fig1 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,0.6), x = c(0,0.8), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.8), x = c(0,0.4), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace(y = c(0.6), x = c(0.8), mode = "c",
text = c("D"), textposition = "top right") %>%
add_trace( x=lista[[1]], y = lista[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.07986,0.07986,0.07986,0.08986,0.09986)*scaling-0.01,
y = c(0.04,0.05,0.06,0.06,0.06)*scaling-0.01,
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.4),y=c(-0.2),text=TeX("\\varphi"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue')))) %>% config(mathjax = "cdn")
fig2 <- plot_ly(data, x = ~x, y = ~y, type = 'scatter', mode = 'lines',
text = c("A","B","D","E","N"),width=width,height=height)
fig2 <- fig2 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,0.6), x = c(0,0.8), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.8), x = c(0,0.4), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace(y = c(0.6), x = c(0.8), mode = "c",
text = c("D"), textposition = "top right") %>%
add_trace( x=~listb[[1]], y = ~listb[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.07986,0.07986,0.07986,0.08986,0.09986)*scaling-0.01,
y = c(0.04,0.05,0.06,0.06,0.06)*scaling-0.01,
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.5),y=c(0.15),text=TeX("\\varphi '"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue')))) %>% config(mathjax = "cdn")
fig3 <- plot_ly(data, x = ~x, y = ~y, type = 'scatter', mode = 'lines',
text = c("A","B","D","E","N"),width=width,height=height)
fig3 <- fig3 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,0.6), x = c(0,0.8), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.8), x = c(0,0.4), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace(y = c(0.6), x = c(0.8), mode = "c",
text = c("D"), textposition = "top right") %>%
add_trace( x=~listc[[1]], y = ~listc[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.07986,0.07986,0.07986,0.06986,0.05986)*scaling,
y = c(0.08,0.07,0.06,0.06,0.06)*scaling,
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.2),y=c(0.4),text=TeX("\\varphi ''"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue')))) %>% config(mathjax = "cdn")
data <- data.frame(xp, yp)
t <- list(family = "Arial", size = 14, color = toRGB("grey50"))
fig4 <- plot_ly(data, x = ~xp, y = ~yp, type = 'scatter', mode = 'lines',
text = c("A","N","E","B","D"),width=width,height=height)
fig4 <- fig4 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,-0.8), x = c(0,0.6), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.1), x = c(0,0.3), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace( x=listd[[1]], y = listd[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.181,0.135,0.09,0.16,0.17),
y = c(-0.24,-0.245,-0.25,-0.3,-0.27),
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.1),y=c(-0.4),text=TeX("\\varphi"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue'))))%>% config(mathjax = "cdn")
fig5 <- plot_ly(data, x = ~xp, y = ~yp, type = 'scatter', mode = 'lines',
text = c("A","N","E","B","D"),width=width,height=height)
fig5 <- fig5 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,-0.8), x = c(0,0.6), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.1), x = c(0,0.3), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace( x=liste[[1]], y = liste[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.181,0.231,0.281,0.181,0.181),
y = c(-0.24,-0.24,-0.24,-0.19,-0.14),
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.4),y=c(-0.2),text=TeX("\\varphi'"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue')))) %>% config(mathjax = "cdn")
fig6 <- plot_ly(data, x = ~xp, y = ~yp, type = 'scatter', mode = 'lines',
text = c("A","N","E","B","D"),width=width,height=height)
fig6 <- fig6 %>% add_markers() %>%
add_text(textfont = t, textposition = "top right") %>%
add_lines(y = c(0,-0.8), x = c(0,0.6), line = list(color = "grey"),
inherit = FALSE, showlegend = FALSE) %>%
add_lines(y = c(0,0.1), x = c(0,0.3), line = list(color = "grey", dash='dot'),
inherit = FALSE, showlegend = FALSE) %>%
add_trace( x=listf[[1]], y = listf[[2]], line = list(shape = "spline",color = '#e763fa'),mode="lines") %>%
#set up the triangle
add_trace(
x = c(0.181,0.231,0.281,0.181,0.181),
y = c(-0.24,-0.24,-0.24,-0.19,-0.14),
fill = 'toself',
fillcolor = '#e763fa',
line = list(
color = '#e763fa'
),
text = "Points only",
mode="lines"
)%>%
add_annotations(x=c(0.4),y=c(-0.2),text=TeX("\\varphi''"),showarrow=FALSE,font=list(size=24)) %>%
layout(xaxis= list(showticklabels = FALSE, title="R", range=c(-1.1,2.2)),
yaxis = list(title="", range=c(-1.2,1.2)), showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1)),
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue'))))%>% config(mathjax = "cdn")
fig <- subplot(fig1,fig2,fig3,fig4,fig5,fig6,nrows=2) %>% config(mathjax = "cdn") %>%
layout(annotations = list(
list(x = 0.1 , y = 1.05, text = TeX("A1:\ \\varphi\\in (0,\\pi)"),
showarrow = F, xref='paper', yref='paper'),
list(x = 0.5 , y = 1.05, text = TeX("A2:\ \\varphi'\\in (-\\frac{\\pi}{2},\\frac{\\pi}{2})"),
showarrow = F, xref='paper', yref='paper'),
list(x = 0.9 , y = 1.05, text = TeX("A3:\ \\varphi''\\in (0,\\pi)"),
showarrow = F, xref='paper', yref='paper')))
figWe present three equivalent definitions of the spherical angle using representation of \(R\) varying the starting axis of the spherical zenith angle, \(\varphi\in (0,\pi)\) (left most column A1), \(\varphi'\in (-\frac{\pi}{2},\frac{\pi}{2})\) (middle column A2), \(\phi''\in (0,\pi)\)(right column A3), where \[\varphi' = \varphi - \frac{\pi}{2}, \varphi''=\pi -\varphi\] By definition then, the \(\varphi'\) on the top of column A2 is positive and on the bottom is negative. We show that \(\varphi = 2\text{arctan}(R)\), \(\varphi' = 2\text{arcsin}\left (\frac{R^2-1}{R^2+1}\right )\), \(\varphi'' = 2\text{arctan}(\frac{1}{R})\) and they are equivalent.
Formula for A1: The derivation for column A1 is straightforward as \(\varphi =2\angle ANB=2\arctan (\frac{R}{1})=2\arctan (R)\)
Formula for A3: \(\varphi''=\pi -\varphi=\pi -2\arctan (R)=2(\frac{\pi}{2}-\arctan (R))=2{\text{arccot}}(R)=2\arctan(\frac{1}{R})\)
Formula for A2: Considering \(AB=R, AD=AN=1, AE\perp BN, AN \perp AB\), We need to explicate \(\phi=\angle DAB\) in terms of \(R\). Let \(\angle NAE=\angle DAE = \beta\). Set \(AE=x\). First we observed the similarity relation \(\triangle EAB \sim \triangle ACB\). Then \(\frac{AE}{AC}=\frac{x}{1}=\frac{R}{\sqrt{R^2+1}}=\frac{AB}{BC}\).
Now \[\begin{equation} \sin(\phi)= \cos(2\beta) = 2\cos^2(\beta)-1=2x^2-1 = \frac{R^2-1}{R^2+1} \end{equation}\ .\]
Therefore, \(\phi=\arcsin\left (\frac{R^2-1}{R^2+1}\right )\).
- Finally, the A1 and A2 equivalence can be verified from the identity \(2\arctan(R)=\arcsin\left (\frac{R^2-1}{R^2+1}\right )+\frac{\pi}{2}\).
# generate the DF with lines
rangeX <- rangeY <- 5 # x-axis range for animation
len <- 20 # number of animation frames
t <- seq(from=-rangeX, to=rangeX, length.out = len) # equidistance P1x coordinates
# Define the projection point P1
P1x <- t; P1y <- 0
# define the *corresponding* point on the circle/sphere
Px <- (2*P1x)/(1+P1x^2+P1y^2)
Py <- (-1+P1x^2 + P1y^2)/(1+P1x^2 + P1y^2)
lineX <- P1x
lineY <- rep(P1y, len)
for (i in 1:length(P1x)) {
if (abs(P1x[i]) < 1) {
lineX[i] <- Px[i]
lineY[i] <- Py[i]
}
}
df <- data.frame(Nx=0, Ny=1, P1x=P1x, P1y=P1y, Px=Px, Py=Py,
lineX=lineX, lineY=lineY, frame=c(1:length(P1x)))
# render the animation
library(plotly)
df %>% plot_ly(x=~Px, y = ~Py, frame=~frame, type='scatter', name="3D Point on Sphere",
mode='lines+markers', marker=list(color='green', size=15), showlegend=F) %>%
add_segments(x = ~Nx, y = ~Ny, xend = ~lineX[frame], yend = ~lineY[frame],
showlegend = F, name="stereographic map") %>%
add_trace(x = ~P1x, y = ~P1y, frame=~frame, showlegend = F,
marker=list(color='red', size=15), name="2D Projection") %>%
layout(title="Bijective Stereographic Projection Map",
shapes = list(list(type = 'circle', name="Sphere cross-section",
xref = 'x', x0 = -1, x1 = 1, yref = 'y', y0 = -1, y1 = 1,
line = list(color = 'blue'))),
xaxis = list(title="Axial plane (X,Y)", range=c(-rangeX,rangeX)),
yaxis = list(title="Azimuthal Axis (Z)", range=c(-rangeY,rangeY)))5 Laplace Transform Applications
5.1 ILT of Kimesurfaces \((F(z))\) to Time-series \((f(t))\)
5.1.1 Example 1: \(F(z)=F_1(z)+F_2(z)\times F_3(z) +F_4(z)\)
Suppose we want to apply the ILT (\(\mathcal{L}^{-1}\) ) to reconstruct a time-series, \(\hat{f}(t)=\mathcal{L}^{-1}(F)(t)\), corresponding to a given a composite kimesurface: \[F(z)=\mathcal{L}(f)=\underbrace{\frac{1}{z+1} }_{F_1(z)=\mathcal{L}\left( f_1(t)=e^{-t}\right )} + \underbrace{\frac{1}{z^2 + 1}}_{F_2(z)=\mathcal{L}\left( f_2(t)=\sin(t)\right )} \times \underbrace{\frac{z}{z^2 + 1}}_{F_3(z)=\mathcal{L}\left( f_3(t)=\cos(t)\right )} + \underbrace{\frac{1}{z^2}}_{F_4(z)=\mathcal{L}\left( f_4(t)=t\right )} .\]
Mind the decomposition of the compound (time and kime) functions in terms of their building blocks, i.e., simpler functions. Using the linearity and convolution-to-product properties of the Laplace transform, we have:
\[F(z)=F_1(z)+F_2(z)\times F_3(z) +F_4(z)\] Therefore, \[f(t)=\mathcal{L}^{-1}\left ( F \right ) = \mathcal{L}^{-1}\left ( F_1+F_2\times F_3 +F_4) \right )=\] \[\mathcal{L}^{-1}( F_1) + \left ( \underbrace{\mathcal{L}^{-1}(F_2) * \mathcal{L}^{-1} (F_3)}_{\text{convolution}}\right )(t) + \mathcal{L}^{-1}(F_4)=\]
\[ \mathcal{L}^{-1}(\mathcal{L}(f_1))(t) + \left ( \mathcal{L}^{-1}(\mathcal{L}(f_2)) * \mathcal{L}^{-1}(\mathcal{L}(f_3))\right )(t) + \mathcal{L}^{-1}(\mathcal{L}(f_4))(t).\]
Finally,
\[f(t)=\mathcal{L}^{-1}(F)(t)=f_1(t) + \left ( f_2 * f_3\right )(t) + f_4(t)=e^{-t}+ \int_{0}^{t}{\sin(\tau) \times \cos(t-\tau)d\tau} +t=t+e^{-t}+\frac{t\sin(t)}{2}. \]
Let’s validate this analytical derivation by using the discrete LT.
### Analytical-Model-based ILT: F(z) = 1/(z+1) + (1/(z^2 + 1))*(z/(z^2 + 1)) + 1/(z^2)
fMRI_Kimesurface <- function (z) {
return ( 1/(z+1) + (1/(z^2 + 1))*(z/(z^2 + 1)) + 1/(z^2) ) # F(z)
}
# Test F=fMRI_kimesurface
# fMRI_Kimesurface(0.5+ 0.5i)
######### Display the actual kimesurface: F(z)
x2 <- seq(from = -2, to = 2, length.out = 200)
y2 <- seq(from = -2, to = 2, length.out = 200)
z2_grid = array(complex(), dim=c(length(x2), length(y2)))
for (i in 1:dim(z2_grid)[1]) {
for (j in 1:dim(z2_grid)[2]) {
z2_grid[i,j] = fMRI_Kimesurface(complex(real=x2[i], imaginary = y2[j]))
}
}
image(Im(z2_grid), axes=FALSE); contour(Im(z2_grid), add=T, col="red", lwd=2, drawlabels=FALSE, axes=FALSE);
contour(Re(z2_grid), add=T, col="blue", lwd=2, drawlabels=F,
main = "Im(f) Image with Re(f) and Im(f) Contours in red and blue", axes=F);
title("Re(f) Image with Re(f) and Im(f) Contours in red and blue", font = 5)
Recall that kimesurfaces are defined over complex-time (Complex-pane)
and have complex-value ranges. The plot below shows the original
kimesurface as a 2D manifold whose height (vertical dimension)
is represented by the Real part of the kimesurface intensity,
Re(fMRI_Kimesurface), and the surface color
encodes the Imaginary part of the kimesurface intensity,
Im(fMRI_Kimesurface).
# Surface height= Re(fMRI_Kimesurface); Surface color = Im(fMRI_Kimesurface)
z3 <- Im(z2_grid)
surf_color <- z3
for (i in 1:dim(z2_grid)[1]) {
for (j in 1:dim(z2_grid)[2]) {
if (z3[i,j]>0) surf_color[i,j] <- min(z3[i,j], 5)
else surf_color[i,j] <- max(z3[i,j], -5)
}
}
colorscale = cbind(seq(0, 1, by=1/(length(x2) - 1)), rainbow(length(x2)))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = Re(z2_grid), # z = Im(z2_grid), # Real or Imaginary part of f(t)
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T #,
# contour=list(x = list(highlight = FALSE),
# y = list(highlight = FALSE),
# z = list( highlight = TRUE, highlightcolor = "blue"),
# color="#000", width=15, lwd=10,
# opacity=1.0, hoverinfo="none")
) %>%
layout(title =
paste0("Kime-Surface, Height=Re(F), Color=Im(F) \n",
"F=LT(f)=1/(z+1) + (1/(z^2 + 1))*(z/(z^2 + 1)) + 1/(z^2)"),
showlegend = FALSE,
scene = list(aspectmode = "manual", zaxis = list(range = c(-5,5)),
aspectratio = list(x=1, y=1, z=0.5)) #,
# xaxis = list(range = c(min(x2), max(x2))),
# yaxis = list(range = c(min(x2), max(x2)))
# )
) # 1:1:1 aspect ratio
pWe can also render the kimesurface as a manifold whose height (vertical dimension) is represented by the Magnitude and color encodes the Phase of the kimesurface.
magnitude <- sqrt((Re(z2_grid))^2+(Im(z2_grid))^2)
x5 <- Re(z2_grid); y5 <- Im(z2_grid)
phase <- atan2(y5,x5)
surf_color <- phase
for (i in 1:200) {
for (j in 1:200) {
if (phase[i,j]>0) surf_color[i,j] <- min(phase[i,j], 5)
else surf_color[i,j] <- max(phase[i,j], -5)
}
}
colorscale = cbind(seq(0, 1, by=1/(length(x2) - 1)), rainbow(length(x2)))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = magnitude, # z = Magnitude
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T,
contour=list(x = list(highlight = FALSE),
y = list(highlight = FALSE),
z = list( highlight = TRUE, highlightcolor = "blue"),
color="#000", width=15, lwd=10,
opacity=1.0, hoverinfo="none")
) %>%
layout(title =
paste0("Kime-Surface, Height=Magnitude(F), Color=Phase(F) \n",
"F=LT(f)=1/(z+1) + (1/(z^2 + 1))*(z/(z^2 + 1)) + 1/(z^2)"),
showlegend = FALSE,
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1, z=0.5),
zaxis = list(range = c(-2,5))
)
) # 1:1:1 aspect ratio
p5.1.1.1 Regularize the Kimesurface
Recall the inverse stereographic projection mapping, which needs to be applied to the values in the range of the Laplace transform.
- Cartesian coordinate inverse stereographic projection mapping:
\[\mathcal{S}^{-1}\left (\underbrace{X}_{Re(z2\_grid)}, \underbrace{Y}_{Im(z2\_grid)} \right )\equiv \left({\frac {2X}{1+X^{2}+Y^{2}}},{\frac {2Y}{1+X^{2}+Y^{2}}}, {\frac {-1+X^{2}+Y^{2}}{1+X^{2}+Y^{2}}}\right)=(x,y,z)\in M.\]
- Spherical coordinate inverse stereographic projection mapping:
\[\mathcal{S}^{-1}\left (\underbrace{\underbrace{R}_{LT\ Magnitude,\ ||(z2\_grid)||}, \overbrace{\Theta}^{LT\ Phase,\ \arctan\left (\frac{Re(z2\_grid)}{Im(z2\_grid)}\right )} }_{2D\ polar\ coordinates} \right )\equiv \\ \left(1, 2\arctan \left (\frac {1}{R}\right ),\Theta \right)= \underbrace{\ \ \left (\overbrace{1}^{sphere\ radius}, \overbrace{\phi}^{zenith\ angle} , \overbrace{\theta}^{azimuth\ } \right )\ \ }_{3D\ spherical\ coordinates} \in M.\]
\[M=\{(r\equiv 1, \phi, \theta)\ |\ 0\leq \phi \leq \pi,\ 0 \leq \theta\leq 2\pi\}.\]
# regularize kimesurface by computing the Stereographic Projection of the LT
# generate the DF with lines
magnitude <- sqrt((Re(z2_grid))^2+(Im(z2_grid))^2)
x5 <- Re(z2_grid); y5 <- Im(z2_grid)
phase <- atan2(y5,x5)
# spherical coordinates str(phi1); str(theta1) # num [1:200, 1:200]
phi1 <- 2*atan2(1, magnitude) # zenith
theta1 <- phase # azimuth
# simple 2D plots of the zenith and azimuth values of the regularized LT
# Mind the tempering of the LT -- no singularities!
#image(phi1); hist(phi1) # Zenith
plot_ly(z=~phi1, type="heatmap") %>%
layout(title="Heatmap of Zenith Angle (phi)") %>% hide_colorbar()plot_ly(x=~as.vector(phi1), type = "histogram", name = "Zenith Angle",
histnorm = "probability") %>%
# add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,
# fill = "tozeroy", name = "Normal Density") %>%
layout(title='Zenith Angle Histogram',
xaxis = list(title = "Zenith"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) # image(theta1); hist(theta1) # Azimuth; mind the symmetry in [-\pi, \pi]
plot_ly(z=~theta1, type="heatmap") %>%
layout(title="Heatmap of Azimuth Angle (theta)") %>% hide_colorbar()plot_ly(x=~as.vector(theta1), type = "histogram", name = "Azimuth Angle",
histnorm = "probability") %>%
# add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,
# fill = "tozeroy", name = "Normal Density") %>%
layout(title='Azimuth Angle Histogram',
xaxis = list(title = "Azimuth"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) # Spherical to Cartesian coordinate mapping
# plot the regularized kimesurface
surf_color <- theta1 # azimuth
x2new <- seq(from = -pi, to = pi, length.out = 200)
colorscale = cbind(seq(0, 1, by=1/(length(x2new) - 1)), rainbow(length(x2new)))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = phi1, # z = zenith
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color (azimuth)
type = 'surface', opacity=1, visible=T #,
# contour=list(x = list(highlight = FALSE),
# y = list(highlight = FALSE),
# z = list( highlight = TRUE, highlightcolor = "blue"),
# color="#000", width=15, lwd=10,
# opacity=1.0, hoverinfo="none")
) %>%
layout(title =
paste0("Regularized Kime-Surface (LT)\n Height=Zenith(F), Color=Azimuth(F) \n",
"F=LT(f)=1/(z+1) + (1/(z^2 + 1))*(z/(z^2 + 1)) + 1/(z^2)"),
showlegend = FALSE,
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1, z=0.5),
zaxis = list(range = c(-2,5))
)
) # 1:1:1 aspect ratio
p5.1.1.2 Inverse Laplace Transform
Let’s now reconstruct and plot the reconstructed 1D function \(f(t)=\mathcal{L}^{-1}(F)(t)\).
This uses the raw (unregularized) kimesurface, but we can also use the regularized Laplace transform (kimesurface).
##### f(t): 1D version #############################################
time_points <- seq(0+0.001, 2*pi, length.out = 160)
# Apply ILT to recover f=ILT(F) and evaluate f(time_points)
# using parallel computing to improve coding speed
cl <- makeCluster(4)
registerDoParallel(cl)
f_result <- foreach(t=1:length(time_points)) %dopar% {
TCIU::ILT(FUNCT=fMRI_Kimesurface, t=time_points[t],
nterms = 31L, m = 1, fail_val = complex(1))
}
stopCluster(cl)
f <- array(complex(1), dim = length(time_points)) # length(f)
f[1:length(f)] = unlist(f_result)
# Convert f=ILT(F) to a data-frame for ggplot
fMRI_time_Intensities_ILT_df <- as.data.frame(cbind(Re=Re(f),Im=Im(f),time_points=time_points))
# ggplot(fMRI_time_Intensities_ILT_df, aes(x=Re(time_points))) +
# geom_line(aes(y = Im, colour="Imaginary"), lwd=3)+
# geom_line(aes(y = Re, colour = "Real"), lwd=3) +
# scale_color_manual("Index",
# breaks=c("Imaginary", "Real"),
# values = c("steelblue", "darkred")) +
# xlab("Time") + ylab("fMRI Image Intensities (f)") +
# # ggtitle("Real (Re) and Imaginary (Im) parts \n of the original fMRI Time-series, f=ILT(F)") +
# labs(title =
# "ILT Reconstructed fMRI Time-series, f=ILT(F)") +
# scale_y_continuous() +
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5))
plot_ly(fMRI_time_Intensities_ILT_df, x=~Re(time_points), y=~Im,
type="scatter", mode="lines", name="Imaginary") %>%
add_trace(x=~Re(time_points), y=~Re, type="scatter", mode="lines",
name="Real") %>%
layout(title="ILT Reconstructed fMRI Time-series, f=ILT(F)",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="fMRI Image Intensities (f)"))5.1.2 Example 2: \(F(z)=|z|\)
Another example is finding the time-series Laplace dual to the kimesurface \[F(z)=|z| .\]
### F(z)= ABS(z) ==> f(t) ~= -1/t
F <- function (z) {
return (abs(z))
}
neg_reciprocal <- array(complex(1), dim = length(time_points));
for (t in 1:length(time_points)) {
neg_reciprocal[t] <- ILT(FUNCT=F, t=time_points[t], nterms = 31L, m = 1, fail_val = complex(1))
}
# Convert f=ILT(F) to a data-frame for ggplot
fMRI_time_Intensities_ILT_df <- as.data.frame(cbind(Re=Re(neg_reciprocal),
Im=Im(neg_reciprocal),time_points=time_points))
# ggplot(fMRI_time_Intensities_ILT_df, aes(x=time_points), ylim=c(-10,10)) +
# geom_line(aes(y = Im, color="Imaginary"), linetype=1, lwd=3)+
# geom_line(aes(y = Re, color = "Real"), lwd=2) +
# ggtitle("Original Time-series, f=ILT(F), F=|z|") +
# xlab("Time") + ylab("Intensities (f)") +
# scale_color_manual(name="Index",
# breaks=c("Imaginary", "Real"),
# values = c("steelblue", "darkred"))+
# scale_y_continuous(limits = c(-5, 1)) +
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5))
plot_ly(fMRI_time_Intensities_ILT_df, x=~time_points, y=~Re, type="scatter", mode="lines",
name="Real") %>%
add_trace(x=~time_points, y=~Im, type="scatter", mode="lines",
name="Imaginary") %>%
layout(title="Original Time-series, f=ILT(F), F=|z|",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity", range=c(-5,0.1)))5.1.3 Example 3: \(F(z)=\frac{w}{z^2+w^2}\)
\[f(t)=\mathcal{L}^{-1}(F_1)(t)=\sin(w t)\ \underset{ILT}{\overset{LT}{\longleftrightarrow}} \ F_1(z)=\mathcal{L}(f_1)(z)=\frac{w}{z^2+w^2}.\]
We can choose \(w=2\) without loss of generosity.
###### f(t)=sin(w t)*u(t) <--> F(z)= w/(z^2+w^2)
# https://en.wikipedia.org/wiki/Laplace_transform
################## Analytical-Model-based ILT: F(z) = 2/(z^2+2^2); w=2
fMRI_Kimesurface <- function (z) {
return (2/(z^2+4))
}
# Test F=fMRI_kimesurface
# fMRI_Kimesurface(1+ 1i) # 1.04-0.78i
# fMRI_Kimesurface(0.5+ 0.5i) # 0.4923077-0.0615385i
######################## Display the actual kimesurface: F(z)
x2 <- seq(from = -2, to = 2, length.out = 200)
y2 <- seq(from = -2, to = 2, length.out = 200)
z2_grid = array(complex(), dim=c(length(x2), length(y2)))
for (i in 1:dim(z2_grid)[1]) {
for (j in 1:dim(z2_grid)[2]) {
z2_grid[i,j] = fMRI_Kimesurface(complex(real=x2[i], imaginary = y2[j]))
}
}
z3 <- Im(z2_grid)
surf_color <- z3
for (i in 1:200) {
for (j in 1:200) {
if (z3[i,j]>0) {surf_color[i,j] <- min(z3[i,j], 5) }
else{surf_color[i,j] <- max(z3[i,j], -5)}
}
}
colorscale = cbind(seq(0, 1, by=1/(length(x2) - 1)), rainbow(length(x2)))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = Re(z2_grid), # z = Im(z2_grid), # Real or Imaginary part of f(t)
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title = "Kime-Surface, Height=Re(F), Color=Im(F) \n F=LT(f)=2/(z^2+4)", showlegend = FALSE,
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1, z=0.5),
zaxis = list(range = c(-2,5))
)
) # 1:1:1 aspect ratio
p##### f(t): 1D version #############################################
time_points <- seq(0+0.001, 2*pi, length.out = 160)
f <- array(NA, dim = length(time_points));
# using parallel computing to speed up the code
cl <- makeCluster(ncor)
registerDoParallel(cl)
ILT_F_result <-
foreach(t = 1:length(time_points)) %dopar% {
TCIU::ILT(FUNCT=fMRI_Kimesurface, t=time_points[t], nterms = 31L, m = 1, fail_val = complex(1))
}
stopCluster(cl)
ILT_F <- array(NA, dim = length(time_points)) # length(f)
ILT_F[1:length(ILT_F)] = unlist(ILT_F_result)
for (t in 1:length(time_points)) {
f[t] = sin(2*time_points[t]) # f(t)=sin(w t)*u(t)
}
# Convert f=ILT(F) to a data-frame for ggplot
fMRI_time_Intensities_ILT_df <- as.data.frame(cbind(time_points=time_points, Re_f=Re(f),Im_f=Im(f),
Re_ILT_F=Re(ILT_F),Im_ILT_F=Im(ILT_F)))
# ggplot(fMRI_time_Intensities_ILT_df, aes(x=time_points)) +
# geom_line(aes(y = Re_f, color = "Sin"), lwd=1, lty="dashed") +
# geom_line(aes(y = Re_ILT_F+0.01, color = "Real f'"), lwd=2, lty="dotted") +
# geom_line(aes(y = Im_ILT_F, color="Imaginary f'"), linetype=1, lwd=2) +
# scale_color_manual("Index",
# values = c("Sin"="darkred","Real f'"="darkgreen", "Imaginary f'"="steelblue")) +
# # ggtitle("Real (Re) and Imaginary (Im) parts \n of the original fMRI Time-series, f=ILT(F)")
# labs(title = "Original fMRI Time-series f(t)=sin(2t) and \n Reconstructed f'(t)=ILT(F)=ILT(LT(f))",
# subtitle = bquote("F" ~ "=" ~ "LT(f)" ~ "=" ~ 2/(z^2+2^2))) +#"F=LT(f)=2/(z^2+2^2)")
# xlab("Time") + ylab("fMRI Image Intensities (f and f')") +
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5),
# plot.subtitle = element_text(hjust = 0.5))
plot_ly(fMRI_time_Intensities_ILT_df, x=~time_points, y=~Re_f, type="scatter", mode="lines",
name="Sine") %>%
add_trace(x=~time_points, y=~Re_ILT_F+0.02, type="scatter", mode="markers",
name="Real f") %>%
add_trace(x=~time_points, y=~Im_ILT_F, type="scatter", mode="lines",
name="Imaginary f'") %>%
layout(title="Original fMRI Time-series f(t)=sin(2t) and \n Reconstructed f'(t)=ILT(F)=ILT(LT(f))",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))5.2 Laplace Transformation of Time-series \((f(t))\) to Kimesurfaces \((F(z))\)
Next, we will demonstrate the reverse operation - finding the Laplace-dual kimesurface corresponding to a given time-series. Starting with the time-series:
5.2.1 Discrete Representation of \(\mathcal{L}^{-1}(\mathcal{L}(\sin(x)))\)
Let’s explore the duality of the LT, (\(\mathcal{L}^{-1}(\mathcal{L})=I?\)), using \(f(t)=\sin(t)\).
#######################################################################
###### ILT on discrete laplace transform of sine (analytical form) ####
#######################################################################
f_sin <- function(t) { sin(t) }
# Define the LT(sin) as a C-valued function
lt_func = function(z) TCIU::LT(f_sin, z)# LT(FUNCT, z), not LT(z, FUNCT) # discrete Laplace Transform of sine
tvalsn <- seq(0, pi*2, length.out = 20)
# using parallel computing to improve coding speed
cl <- makeCluster(ncor)
registerDoParallel(cl)
sinvalsn <- foreach(t=1:length(tvalsn),
.export='cubintegrate',
.packages='cubature') %dopar% {
TCIU::ILT(FUNCT=lt_func, t=tvalsn[t])
}
stopCluster(cl)
sinvalsn = unlist(sinvalsn)
sinvalsn_df <- as.data.frame(cbind(Re=Re(sinvalsn),Im=Im(sinvalsn), Sin=sin(tvalsn), time_points=tvalsn))
# ggplot(sinvalsn_df, aes(x=time_points))+
# geom_line(aes(y=Re, color="Real"), linetype=1, lwd=2) +
# geom_line(aes(y = Sin, color="Sin"), linetype=2, lwd=1) +
# scale_color_manual(name="Index",
# values = c("Real"="steelblue", "Sin"="darkred"))+
# labs(title = "Original fMRI Time-series f(t)=sin(t) and \n Reconstructed f'(t)=ILT(F)=ILT( discrete LT(f))",
# subtitle = bquote("F" ~ "=" ~ "discrete LT(f)")) + #"F=LT(f)=2/(z^2+2^2)"
# xlab("Time") + ylab("fMRI Image Intensities (f and f')") +
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5),
# plot.subtitle = element_text(hjust = 0.5))
plot_ly(sinvalsn_df, x=~time_points, y=~Re, type="scatter", mode="lines",
name="Real") %>%
add_trace(x=~time_points, y=~Sin, type="scatter", mode="lines",
name="Sine") %>%
add_trace(x=~time_points, y=~Im, type="scatter", mode="lines",
name="Imaginary") %>%
layout(title="Original fMRI Time-series f(t)=sin(t) and \n Reconstructed f'(t)=ILT(F)=ILT( discrete LT(f))",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))5.3 Foliating the Laplace Transform to Synthetically Generate Simulated Time-series
\[\begin{array}{ccc} & LT & \\ \overbrace{f(t):\mathbb{R}^+\to \mathbb{C}}^{spacetime\ signal} & \xrightarrow{\hspace{2cm}} & \overbrace{F(\kappa):\mathbb{C}\to \mathbb{C}}^{LT(signal)} \\ & & \\ \overset{Samples}{\theta=\theta_o} \Bigg\updownarrow & & \Bigg\downarrow \overset{Foliate}{\overset{\kappa=te^{i\theta}}{\theta=\theta_o}}\\ & ILT & \\ \underbrace{g(t|\theta_o):\mathbb{R}^+\to \mathbb{C}}_{SK\ simulated\ data} & \xleftarrow{\hspace{2cm}} & \underbrace{F(\kappa=te^{i\theta_o}):\mathbb{C}\to \mathbb{C}}_{kimesurface\ leaf} \\ \end{array}\]
See the TCIU manifold foliation notes.
#######################################################################
###### ILT on discrete Laplace transform of sine (analytical form) ####
#######################################################################
f_sin <- function(t) { sin(t) }
kimesurfaceFoliation <- function (fun=f_sin, n=50, theta=pi/4, centerOffset=10) {
# Compute the LT(sin) Foliation (at \theta=\pi/2) a C-valued function
# n=100 # sample size
z <- complex(length.out = n, real = 0, imaginary = 0) # initialize the kimesurface foliation domain
# theta <- pi/4 # foliation along a fixed kime-phase (theta= pi/4)
theta <- rep(theta, n) # initialize the kime-phase of hte foliation plane
off <- centerOffset # offset the start of hte foliation plane, if needed
for (i in c(1:n)) { # linear grid over foliation plane at \theta=\pi/4
the = theta[i] # in general, theta ~ \Phi (random sample from the phase distribution)
radius = (i*2*pi)/n # time restricted to [0, 2\pi]
z[i] = complex(real = off + radius*cos(the), imaginary = off + radius*sin(the))
}
# Generic LT function - kimesurface, defined over each kime in the complex domain
lt_func = function(z) TCIU::LT(fun, z)# LT(FUNCT, z), not LT(z, FUNCT) # discrete Laplace Transform of sine
# Define the leaf-generating function intersecting the foliating plane with the kimesurface
lt_func_Leaf = function(z) {
r = abs(z) # kime-magnitude
# phi = atan2(Im(z), Re(z)) # kime-phase
return (lt_func(complex(real=r*cos(theta), imaginary=r*sin(theta))))
}
tvalsn <- seq(0, pi*2, length.out = n)
# using parallel computing to improve coding speed
cl <- makeCluster(detectCores()-3)
registerDoParallel(cl)
sinvalsn <- foreach(t=1:length(tvalsn),
.export='cubintegrate',
.packages='cubature') %dopar% {
TCIU::ILT(FUNCT=lt_func_Leaf, t=tvalsn[t])
}
stopCluster(cl)
sinvalsn = unlist(sinvalsn)
reMin <- min(Re(sinvalsn))
imMin <- min(Im(sinvalsn))
sinvalsn_df <- as.data.frame(cbind(Re=log(2+abs(reMin)+Re(sinvalsn)),
Im=log(2+abs(imMin)+Im(sinvalsn)), Sin=fun(tvalsn), time_points=tvalsn))
return(sinvalsn_df)
}
# Sample 1: \theta_1 = pi/4
# sinvalsn_df1 <- kimesurfaceFoliation(fun=f_sin, n=50, theta=pi/4, centerOffset=10)
sinvalsn_df1 <- kimesurfaceFoliation(fun=f_sin, n=50, theta=pi/(2.5), centerOffset=10)
# Sample 2: \theta_2 = pi/3
sinvalsn_df2 <- kimesurfaceFoliation(fun=f_sin, n=50, theta=pi/3, centerOffset=10)
# ggplot(sinvalsn_df2, aes(x=time_points))+
# geom_line(aes(y=Re, color="Real"), linetype=1, lwd=2) +
# geom_line(aes(y = Sin, color="Sin"), linetype=2, lwd=1) +
# scale_color_manual(name="Index",
# values = c("Real"="steelblue", "Sin"="darkred"))+
# labs(title = "Original fMRI Time-series f(t)=sin(t) and \n Reconstructed f'(t)=ILT(F)=ILT( discrete LT(f))",
# subtitle = bquote("F" ~ "=" ~ "discrete LT(f)")) + #"F=LT(f)=2/(z^2+2^2)"
# xlab("Time") + ylab("fMRI Image Intensities (f and f')") +
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5),
# plot.subtitle = element_text(hjust = 0.5))
# Show Real and Imaginary timeseries
plot_ly(x=sinvalsn_df1$time_points,y=sinvalsn_df1$Sin, type="scatter", mode="lines", name="Sine") %>%
# first sample, corresponding to \theta= pi/2
add_trace(x=sinvalsn_df1$time_points, y=sinvalsn_df1$Re, type="scatter", mode="lines",
name="Re(theta=pi/4)") %>%
add_trace(x=sinvalsn_df1$time_points, y=sinvalsn_df1$Im, type="scatter", mode="lines",
name="Im(theta=pi/4)") %>%
# second sample, corresponding to \theta= pi/3
add_trace(x=sinvalsn_df2$time_points, y=sinvalsn_df2$Re, type="scatter", mode="lines",
name="Re(theta=pi/3)") %>%
add_trace(x=sinvalsn_df2$time_points, y=sinvalsn_df2$Im, type="scatter", mode="lines",
name="Im(theta=pi/3)") %>%
layout(title="Original Time-series f(t)=sin(t) and Random Time-series Sampling by \n Inverting the Foliated Kimesurface, g(t)=ILT(F(t,theta))",
legend = list(orientation = 'h', y=-0.2),
xaxis=list(title="Time"), yaxis=list(title="Real and Imaginaty parts\n Foliation ILT(Leaf Function) Intensity"))# Plot only the Magnitudes
plot_ly(x=sinvalsn_df1$time_points,y=sinvalsn_df1$Sin, type="scatter", mode="lines", name="Sine") %>%
# first sample, corresponding to \theta= pi/2
add_trace(x=sinvalsn_df1$time_points, y=sqrt((sinvalsn_df1$Re)^2+(sinvalsn_df1$Im)^2),
type="scatter", mode="lines", name="MAG(theta=pi/4)") %>%
# second sample, corresponding to \theta= pi/3
add_trace(x=sinvalsn_df2$time_points, y=sqrt((sinvalsn_df2$Re)^2+(sinvalsn_df2$Im)^2),
type="scatter", mode="lines", name="MAG(theta=pi/3)") %>%
layout(title="Original Time-series f(t)=sin(t) and Random Time-series Sampling by \n Inverting the Foliated Kimesurface, g(t)=ILT(F(t,theta))",
legend = list(orientation = 'h', y=-0.2),
xaxis=list(title="Time"), yaxis=list(title="Magnitude of Foliated ILT(Leaf Function) Intensity"))6 ILT of the explicit definition of the kimesurface function \(F(z)=\mathcal{L}(\sin)(z)\)
Next, we will compare the original time-signal \(f(t)=\sin(t)\) against a spline-smoothed \(f_1(t)=\mathcal{L}^{-1}(\mathcal{L}(\sin))(t)\) over the range \(t\in [0 : 2\pi]\).
# reduce the grid-resolution from 200*200 down to (50-1)*(50-1)
range_limit = 2
x2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# drop the first row to avoid real part value of 0
y2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# drop the first column to to avoid imaginary part value of 0
# Recompute the LT(sin) discretized on lower-res grid
z2_grid = array(dim=c(length(x2), length(y2)))# x2 %o% y2
f_sin <- function(t) { sin(t) }
# kime surface transform
# use parallel computing to speed up code
cl <- makeCluster(ncor)
registerDoParallel(cl)
F = list()
for (i in 1:length(x2) ){
F[[i]] =
foreach(j = 1:length(y2),
.export='cubintegrate',
.packages='cubature') %dopar% {
F_result = TCIU::LT(FUNCT=f_sin, complex(real=x2[i], imaginary = y2[j]))# LT(FUNCT, z), not LT(z, FUNCT)
mag = log(sqrt( Re(F_result)^2+ Im(F_result)^2))
# log-transform the magnitude to temper the kimesurface amplitude
phase = atan2(Im(F_result), Re(F_result))
mag * exp(1i*phase)#z2_grid[i,j] =
}
}
stopCluster(cl)
F_vec = lapply(F, unlist)
z2_grid = unlist(do.call(rbind, F_vec))
#### define Kimesurface_fun #####
Kimesurface_fun <- function (z, array_2D) {
# array_2D <- z2_grid
# convert z in C to Cartesian (x,y) coordinates
x1 <- ceiling(Re(z))-1; # if (x1<2 || x1>dim(array_2D)[1]) x1 <- 2
y1 <- ceiling(Im(z))-1; # if (y1<2 || y1>dim(array_2D)[2]) y1 <- 2
# if exceed the domain use the default 1
if(!is.na(x1)){
if((x1 < 1) || (x1 > dim(array_2D)[1])){ x1 <- 1 }
}
if(!is.na(y1)){
if((y1 < 1) || (y1 > dim(array_2D)[2])){ y1 <- 1 }
}
# Exponentiate to Invert the prior (LT) log-transform of the kimesurface magnitude
val1 = complex(length.out=1, real=Re(array_2D[x1, y1]), imaginary = Im(array_2D[x1, y1]))
mag = exp(sqrt( Re(val1)^2+ Im(val1)^2))
# mag = sqrt( Re(val1)^2+ Im(val1)^2)
phase = atan2(Im(val1), Re(val1))
value <- complex(real=Re(mag * exp(1i*phase)), imaginary = Im(mag * exp(1i*phase)))
return ( value )
}
# Kimesurface_fun(5+5i, z2_grid)
##### Time-domain grid (regular equidistant positive real break points)
time_points <- seq(0+0.001, 2*pi, length.out = 160)
f2 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f2[t] <- ILT(FUNCT=function(z) Kimesurface_fun(z, array_2D=z2_grid),
t= time_points[t], nterms = 31L,
m = 1, gamma=0.5, fail_val = complex(1), msg=TRUE)
}
# interpolate f(t) to 20 samples in [0 : 2*pi]. Note that 20~2*PI^2
tvalsn <- seq(0, pi*2, length.out = 20)
f3 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f3[t] <- f2[ceiling(t/20)]
}
f4 <- smooth.spline(time_points, Re(f3), spar = 1)$y; # plot(f4)
time_Intensities_ILT_df2 <- as.data.frame(cbind(Re=Re(f4), Im=Re(f3),
Sin=sin(time_points),
time_points=time_points))
min_range <- range(Re(f4))[1]; max_range <- range(Re(f4))[2]
time_Intensities_ILT_df2$Re <- time_Intensities_ILT_df2$Re/(8*max_range/9)
time_Intensities_ILT_df2$Im <- time_Intensities_ILT_df2$Im/(8*max_range/9)
colnames(time_Intensities_ILT_df2) <-
c("Smooth Reconstruction", "Raw Reconstruction", "Original sin()", "time_points")
df <- reshape2::melt(time_Intensities_ILT_df2, id.var = "time_points")
# ggplot(df, aes(x = time_points, y = value, colour = variable)) +
# geom_line(linetype=1, lwd=3)+
# ylab("Function Intensity") + xlab("Time") +
# theme(legend.position="top")+
# labs(title=
# "Comparison between f(t)=sin(t) and SplineSmooth(ILT(LT(sin)))(t); Range [" ~ 0 ~":"~ 2*pi~"]")
plot_ly(df, x=~time_points, y=~value, type="scatter", mode="lines",
color =~variable, name=~variable) %>%
layout(title="Comparison between f(t)=sin(t) and SplineSmooth(ILT(LT(sin)))(t); Range=[0,2pi]",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))7 Compare the analytical and discrete Laplace transforms
Using the function \(f(t)=\sin(t)\) and its closed-form LT, \(F(z)=\mathcal{L}(f)(z)=\frac{1}{z^2 + 1}\), we will compare the analytical and discrete kime-surface reconstructions.
# range_limit = 2
# x2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# # drop the first row to avoid real part value of 0
# y2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# # drop the first column to to avoid imaginary part value of 0
#############################################################################
### Explicit (continuous-function) form of laplace transformation of sine ###
#############################################################################
XY = expand.grid(X=x2,Y=y2)
# XY
laplace_sine = function(p) { 1/(p^2 + 1) } # Exact laplace transform of sin(x), continuous function
complex_xy = mapply(complex, real=XY$X,imaginary=XY$Y)
sine_z =laplace_sine(complex_xy)
dim(sine_z) = c(length(x2), length(y2));
# dim(sine_z) # [1] 49 49
# Log-transform the Intensity of the Kime surface (to make a better-looking surface plot, not a big deal)
mag_sine_z = log(sqrt( Re(sine_z)^2+ Im(sine_z)^2))
# log-transform the magnitude to temper the kimesurface amplitude
phase_sine_z = atan2(Im(sine_z), Re(sine_z))
sine_z_new = mag_sine_z * exp(1i*phase_sine_z)
# draw the surface whose heights are based on real, imaginary,
# magnitude and phase of the both results respectively
# the two kinds of result are from closed-formed LT and discrete LT of sine
plot_ly(hoverinfo="none", showscale = FALSE)%>%
add_trace(z=Re(sine_z_new)-0.5, type="surface", surfacecolor=phase_sine_z) %>%
add_trace(z = Re(z2_grid), type="surface", opacity=0.7, surfacecolor=Im(z2_grid) )%>%
layout(title =
"Kime-Surface, LT(sin()), Height=Re(LT(sin())), Color=Re(LT(sin())) \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform", showlegend = FALSE)plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z=Im(sine_z_new)-0.5, type="surface") %>%
add_trace(z = Im(z2_grid), type="surface", opacity=0.7)%>%
layout(title =
"Kime-Surface, LT(sin()), Height=Im(LT(sin())) \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform", showlegend = FALSE)mag_sine_z_new = sqrt(Re(sine_z_new)^2 + Im(sine_z_new)^2)
mag_z2_grid = sqrt( Re(z2_grid)^2+ Im(z2_grid)^2)
plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = mag_sine_z_new-0.5, type="surface") %>%
add_trace(z = mag_z2_grid, type="surface", opacity=0.7) %>%
layout(title =
"Kime-Surface, LT(sin()), Height=Magnitude(LT(sin())) \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform", showlegend = FALSE)phase_sine_z_new = atan2(Im(sine_z_new), Re(sine_z_new))
phase_z2_grid = atan2(Im(z2_grid), Re(z2_grid))
plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = phase_sine_z_new-0.5, type="surface") %>%
add_trace(z = phase_z2_grid, type="surface", opacity=0.7) %>%
layout(title =
"Kime-Surface, LT(sin()), Height=Phase(LT(sin())) \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform", showlegend = FALSE)plot_ly(hoverinfo="none", showscale = FALSE)%>%
add_trace(z=mag_sine_z_new-0.5, surfacecolor=phase_sine_z_new, type="surface") %>%
add_trace(z=mag_z2_grid, surfacecolor=phase_z2_grid, type="surface", opacity=0.7) %>%
layout(title =
"Kime-Surface, LT(sin()), Height=Phase(LT(sin())), Color=Phase \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform",
showlegend = FALSE,
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1, z=1.0),
zaxis = list(range = c(-2,3))
)
)In these 3D displays, we intentionally offset vertically the two kimesurfaces to illustrate their identical geometric, topological, and shape characteristics. This validates that the theoretical and the pragmatic implementation of the Laplace Transform agree on the example function, \(f(t)=sin(t)\).
8 Compare the analytical and discrete Inverse Laplace Transforms
We can also similarly validate the agreement between the exact (continuous) theoretical and the approximate (discrete) Inverse Laplace Transforms using the kimesurface, \(F(z)=\mathcal{L}(f)(z)=\frac{1}{z^2 + 1}\), corresponding to the original function \(f(t)=sin(t)\).
# Reconstruct and plot the 1D function f(t) = ILT(LT(sin(x))(z))(t) and compare to sin(t)
################ ILT on analytic LT of sine ###########################################################
##### Time-domain grid (regular equidistant positive real break points)
time_points <- seq(0+0.001, 2*pi, length.out = 160)
f20 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f20[t] <- ILT(FUNCT=function(z) Kimesurface_fun(z, array_2D=sine_z_new),
t= time_points[t], nterms = 31L,
m = 1, gamma=0.5, fail_val = complex(1), msg=TRUE)
}
# interpolate f(t) to 20 samples in [0 : 2*pi]. Note that 20~2*PI^2
tvalsn <- seq(0, pi*2, length.out = 20)
f30 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f30[t] <- f20[ceiling(t/20)]
}
f40 <- smooth.spline(time_points, Re(f30), spar = 1)$y; # plot(f4)
################ ILT on discrete LT of sine ###########################################################
##### Time-domain grid (regular equidistant positive real break points)
time_points <- seq(0+0.001, 2*pi, length.out = 160)
f2 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f2[t] <- ILT(FUNCT=function(z) Kimesurface_fun(z, array_2D=z2_grid),
t= time_points[t], nterms = 31L,
m = 1, gamma=0.5, fail_val = complex(1), msg=TRUE)
}
# interpolate f(t) to 20 samples in [0 : 2*pi]. Note that 20~2*PI^2
tvalsn <- seq(0, pi*2, length.out = 20)
f3 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f3[t] <- f2[ceiling(t/20)]
}
f4 <- smooth.spline(time_points, Re(f3), spar = 1)$y; # plot(f4)
############## make the comparison plot ###############################################################
time_Intensities_ILT_df2 <- as.data.frame(cbind(Re_discrete=Re(f4), Re_analytic=Re(f40),
Sin=sin(time_points),
time_points=time_points))
min_range_Re_discrete <- range(Re(f4))[1]; max_range_Re_discrete <- range(Re(f4))[2]
time_Intensities_ILT_df2$Re_discrete <- time_Intensities_ILT_df2$Re_discrete/(8*max_range_Re_discrete/9)
min_range_Re_analytic <- range(Re(f40))[1]; max_range_Re_analytic <- range(Re(f40))[2]
time_Intensities_ILT_df2$Re_analytic <- time_Intensities_ILT_df2$Re_analytic/(8*max_range_Re_analytic/9)
# ggplot(time_Intensities_ILT_df2, aes(x=time_points))+
# geom_line(aes(y=Re_discrete, color="Smooth Reconstruction of Discrete LT"), lwd=1, lty="dashed") +
# geom_line(aes(y=Re_analytic, color="Analytic Reconstruction of Analytic LT"), lwd=2,lty="dotted") +
# geom_line(aes(y = Sin, color="Original Sin"), lwd=2) +
# scale_color_manual(name="Index",
# values = c("Smooth Reconstruction of Discrete LT"="#F8766D",
# "Analytic Reconstruction of Analytic LT"="#00BA38",
# "Original Sin"="#619CFF"))+
# labs(title = "Comparison among f(t)=sin(t), SplineSmooth(ILT(LT(sin)))(t), \n
# and Analytic ILT(LT(sin))(t); Range [" ~ 0 ~":"~ 2*pi~"]") +
# xlab("Time") + ylab("Function Intensity") +
# theme(legend.title = element_text(size=14, color = "black", face="bold"),
# legend.position = "top",
# panel.grid.minor.y = element_blank(),
# panel.grid.major.y = element_blank(),
# plot.title = element_text(hjust = 0.5),
# plot.subtitle = element_text(hjust = 0.5))
plot_ly(time_Intensities_ILT_df2, x=~time_points, y=~Re_discrete, type="scatter", mode="lines", name="Smooth Reconstruction of Discrete LT") %>%
add_trace(x=~time_points, y=~Re_analytic, type="scatter", mode="markers",
name="Analytic Reconstruction of Analytic LT") %>%
add_trace(x=~time_points, y=~sin(time_points), type="scatter", mode="lines",
name="Original Sine Function") %>%
layout(title="Comparison among f(t)=sin(t), SplineSmooth(ILT(LT(sin)))(t), \n
and Analytic ILT(LT(sin))(t); Range=[0,2pi]",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))9 fMRI example
Let’s demonstrate the Laplace dual of a single time-series representing the real-valued intensity of a functional magnetic resonance imaging (fMRI) data at a fixed spatial location (single voxel in the brain). The entire fMRI data is a 4D hypervolume with intensities encoding the blood oxygenation level dependence at a specific spacetime location \((x,y,z,t)\). For simplicity, we are only focusing on one fixed spatial voxel location \((x_o,y_o,z_o)\). Start by loading and plotting the original fMRI data, \(f_{(x_o,y_o,z_o)}(t)\).
As the fMRI series is extremely noisy (almost nowhere differentiable), we need to apply preprocess filtering (smoothing) to ensure the LT is well-defined, i.e., the fMRI time-series is integrable against the exponential kernel.
9.1 Load the fMRI data
fMRIURL <- "https://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE)
# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec
fMRIVolDims <- dim(fMRIVolume);
# fMRIVolDims
time_dim <- fMRIVolDims[4];
# time_dim ## 180
# 2. extract the time-course of 1D mid-axial slice (3D) hypervolume
xA_fMRI_1D_x20_y20_z11 <- fMRIVolume[20, 20, 11, ]; # length(xA_fMRI_1D_x20_y20_z11) # 180
# hist(xA_fMRI_1D_x20_y20_z11)
# Now, combine your two 1D timeseries into one dataframe for joint hist plotting as densities.
# First make a new column in each that will be
# a variable to identify where they came from later.
times <- c(1:180)
# Smooth noisy f(t)=fMRI
f <- array(complex(1), dim = length(times)); # length(f)
f <- smooth.spline(times, xA_fMRI_1D_x20_y20_z11, df = 10)$y # smooth f (fMRI)
# Construct DF for ggplot
xA_df <- as.data.frame(cbind(times, xA_fMRI_1D_x20_y20_z11, f))
colnames(xA_df) <- c("time", "xA", "smooth")
# ggplot(xA_df, aes(x=times)) +
# geom_point(aes(y = scale(xA)), colour = "darkred", size = 3)+
# geom_line(aes(y = scale(xA)), color="steelblue", lty=1, lwd=1)+
# geom_line(aes(y = scale(smooth)), color = "green", lwd=2) + ### , lty="dashed") +
# ggtitle("Original fMRI Time-series f(t) and Smoothed curve") +
# xlab("Time") + ylab("fMRI Image Intensities (f and f')") +
# scale_color_manual(name="Index",labels=c("Real", "Imaginary"),
# values = c("Re"="darkred", "Im"="steelblue"))+
# theme_grey(base_size = 16) +
# theme(legend.title = element_text(size=14, color = "salmon", face="bold"),
# legend.position="top",
# axis.text.x = element_blank(),
# axis.text.y = element_blank(),
# axis.ticks = element_blank())
plot_ly(xA_df, x=~time, y=~scale(xA), type="scatter", mode="markers+lines", name="fMRI") %>%
add_trace(x=~time, y=~scale(smooth), type="scatter", mode="markers+lines",
name="Smooth fMRI") %>%
layout(title="Original fMRI Time-series f(t) and Smoothed Interpolation",
legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))9.2 Compute the LT \(\mathcal{L}(fMRI)(z)\)
Next, we will compute and display in 3D the Laplace-dual, \(F_{(x_o,y_o,z_o)}(z)=\mathcal{L}(f_{(x_o,y_o,z_o)}(z)\) and display the kimesurface as a complex-valued, complex-time, analytic (holomorphic) function.
Warning: This \(\mathcal{L}(fMRI)(z)\) kimesurface calculation is intensive. To ensure near-real-time computation, below we only show a low-resolution \(49\times49\) kimesurface representation where the kime-manifold height and color represent the kimesurface magnitude and phase, respectively.
time_points <- seq(0+0.001, 2*pi, length.out = 180)
f <- array(complex(1), dim = length(time_points)); # length(f)
# Instead of using the extremely noisy fMRI data, avoid integration problems,
# smooth "f" and use the **smooth version, f**
f <- smooth.spline(ceiling((180*time_points)/(2*pi)), xA_fMRI_1D_x20_y20_z11, df = 10)$y # smooth f (fMRI)
# Define the f(t)=smooth(fMRI)(t) signal as a function of real time 0<t<=2*pi
f_funct <- function(t) {
# f <- match.fun(f) # input R-valued FUNCT should be interpreted as f(t)=fMRI(t) time-series function
if (t < 0+0.001 || t > 2*pi) { return ( 0 ) # sprintf("Out of Range ...")
} else {
return ( f[ceiling((180*t)/(2*pi))] )
}
}
# Evaluate the LT at z, F(z)=LT(f)(z)
z= 1+1i;
# LT(f_funct, z) # Complex-domain value
######################## Display the actual kimesurface: F(z)
range_limit = 2
x2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# drop the first row to avoid real part value of 0
y2 = seq(from = 0, to = range_limit, length.out = range_limit*25)[2:50]
# drop the first column to to avoid imaginary part value of 0
#
# # Recompute the LT(sin) discretized on lower-res grid
# z2_grid = array(dim=c(length(x2), length(y2)))# x2 %o% y2
cl <- makeCluster(ncor)
registerDoParallel(cl)
F = list()
for (i in 1:length(x2) ){
F[[i]] =
foreach(j = 1:length(y2),
.export='cubintegrate',
.packages='cubature') %dopar% {
F_result = TCIU::LT(f_funct, complex(real=x2[i], imaginary = y2[j]))
mag = log(sqrt( Re(F_result)^2+ Im(F_result)^2))
# log-transform the magnitude to temper the kimesurface amplitude
phase = atan2(Im(F_result), Re(F_result))
mag * exp(1i*phase)#z2_grid[i,j] =
}
}
stopCluster(cl)
F_vec = lapply(F, unlist)
z2_grid = unlist(do.call(rbind, F_vec))
surf_color <- atan2(Im(z2_grid), Re(z2_grid))
colorscale = cbind(seq(0, 1, by=1/(length(x2) - 1)), rainbow(length(x2)))
magnitude <- (sqrt( Re(z2_grid)^2+ Im(z2_grid)^2))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = magnitude,
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title = "fMRI Kime-Surface, F=LT(fMRI) \n Height=Mag(F), Color=Phase(F)", showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0) ) ) # 1:1:1 aspect ratio
p9.3 Regularize the fMRI Kimesurface
Recall the earlier discussion of the inverse stereographic projection mapping. In this case, it is applied to the values in the range of the fMRI Laplace transform.
# regularize kimesurface by computing the Stereographic Projection of the LT
magnitude <- sqrt((Re(z2_grid))^2+(Im(z2_grid))^2)
x5 <- Re(z2_grid); y5 <- Im(z2_grid)
phase <- atan2(y5,x5)
# spherical coordinates str(phi1); str(theta1) # num [1:200, 1:200]
phi1 <- 2*atan2(1, magnitude) # zenith
theta1 <- phase # azimuth
# simple 2D plots of the zenith and azimuth values of the regularized LT
# Mind the tempering of the LT -- no singularities!
# image(phi1); hist(phi1)
plot_ly(z=~phi1, type="heatmap") %>%
layout(title="Heatmap of Zenith Angle (theta)") %>% hide_colorbar()plot_ly(x=~as.vector(phi1), type = "histogram", name = "Zenith Angle",
histnorm = "probability") %>%
# add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,
# fill = "tozeroy", name = "Normal Density") %>%
layout(title='Zenith Angle Histogram',
xaxis = list(title = "Zenith"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) # image(theta1); hist(theta1)
plot_ly(z=~theta1, type="heatmap") %>%
layout(title="Heatmap of Azimuth Angle (theta)") %>% hide_colorbar()plot_ly(x=~as.vector(theta1), type = "histogram", name = "Azimuth Angle",
histnorm = "probability") %>%
# add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,
# fill = "tozeroy", name = "Normal Density") %>%
layout(title='Azimuth Angle Histogram',
xaxis = list(title = "Azimuth"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) # Spherical to Cartesian coordinate mapping
# plot the regularized kimesurface
surf_color <- theta1 # azimuth
x2new <- seq(from = -pi, to = pi, length.out = 200)
colorscale = cbind(seq(0, 1, by=1/(length(x2new) - 1)), rainbow(length(x2new)))
p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
add_trace(z = phi1, # z = zenith
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color (azimuth)
type = 'surface', opacity=1, visible=T #,
# contour=list(x = list(highlight = FALSE),
# y = list(highlight = FALSE),
# z = list( highlight = TRUE, highlightcolor = "blue"),
# color="#000", width=15, lwd=10,
# opacity=1.0, hoverinfo="none")
) %>%
layout(title =
paste0("Regularized fMRI Kime-Surface, F=LT(fMRI) \n",
"Height=Zenith(F), Color=Azimuth(F) \n"),
showlegend = FALSE,
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1, z=0.5),
zaxis = list(range = c(-2,5))
)
) # 1:1:1 aspect ratio
p9.4 Invert the LT \(\mathcal{L}^{-1}(\mathcal{L}(fMRI))(t)\)
Finally, we can invert the fMRI kimesurface back in the time-domain and compare it to the original fMRI time-series.
##### Time-domain grid (regular equidistant positive real break points)
# using parallel computing to speed up code
cl <- makeCluster(ncor)
registerDoParallel(cl)
f2_reuslt <-
foreach (t = 1:length(time_points)) %dopar% {
TCIU::ILT(FUNCT=function(z) Kimesurface_fun(z, array_2D=z2_grid),
t= time_points[t], nterms = 31L,
m = 1, gamma=0.5, fail_val = complex(1), msg=TRUE)
}
stopCluster(cl)
###
f2 <- array(complex(1), dim = length(time_points))
for(t in 1:length(time_points)){
f2[t] = ILT(FUNCT=function(z) Kimesurface_fun(z, array_2D=z2_grid),
t= time_points[t], nterms = 31L,
m = 1, gamma=0.5, fail_val = complex(1), msg=TRUE)}
###
f2 <- array(complex(1), dim = length(time_points))# length(f), f=ILT(F)
f2[1:length(f2)] = unlist(f2_reuslt)
# interpolate f(t) to 20 samples in [0 : 2*pi]. Note that 20~2*PI^2
tvalsn <- seq(0, pi*2, length.out = 20)
f3 <- array(complex(1), dim = length(time_points)); # length(f), f=ILT(F)
for (t in 1:length(time_points)) {
f3[t] <- f2[ceiling(t/20)]
}
f4 <- smooth.spline(time_points, Re(f3), spar = 1)$y; # plot(f4)
# scale all 3 signals to get congruence ...
# time_Intensities_ILT_df2 <- as.data.frame(cbind(Re=scale(Re(f4)), Im=scale(Re(f3)),
# fMRI=scale(Re(f_funct(time_points))),
# time_points=time_points))
myFunct <- time_points
for (i in 1:length(time_points)) {
myFunct[i] <- f_funct(time_points[i])
}
time_Intensities_ILT_df2 <- as.data.frame(cbind(Re=scale(Re(f4)), Im=scale(Re(f3)),
fMRI=scale(Re(myFunct)),
time_points=time_points))
min_range <- range(Re(f4))[1]; max_range <- range(Re(f4))[2]
colnames(time_Intensities_ILT_df2) <-
c("Smooth Reconstruction", "Raw Reconstruction", "Original fMRI", "time_points")
df <- reshape2::melt(time_Intensities_ILT_df2, id.var = "time_points")
# ggplot(df, aes(x = time_points, y = value, colour = variable)) +
# geom_line(linetype=1, lwd=3)+
# ylab("Function Intensity") + xlab("Time") +
# theme(legend.position="top")+
# labs(title="Comparison between f(t)=fMRI(t) and SplineSmooth(ILT(LT(fMRI)))(t); Range [" ~ 0 ~":"~ 2*pi~"]")
plot_ly(df, x =~time_points, y = ~value, color = ~variable, type="scatter",
mode="lines", name=~variable) %>%
layout(title="Comparison between f(t)=fMRI(t) and SplineSmooth(ILT(LT(fMRI)))(t);
Range [0,2pi]", legend = list(orientation = 'h'),
xaxis=list(title="Time"), yaxis=list(title="Function Intensity"))The observed discrepancy between the smoothed original fMRI time-series and the raw and smoothed reconstructions, ILT(LT(fMRI)), is due to the over-simplified coarse representation of the kimesurface and the discrete transform approximation to the continuous Laplace transform (forward and inverse).
This example illustrates the foundation of spacekime analytics, where we derive statistical inference by lifting the observed time-series into kime-surfaces, which contain additional geometric and topological information that can then be modeled, interrogated, and analyzed to obtain robust prediction or forecasting.