--- title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)" subtitle: "TCIU: The Laplace Transform of Longitudinal Data: Analytic Duality of Time-series and Kimesurfaces" author: "SOCR Team " date: "`r format(Sys.time(),'%m/%d/%Y')`" # header-includes: # - \usepackage{tikzcd} output: html_document: theme: spacelab highlight: tango includes: before_body: TCIU_header.html toc: true number_sections: true toc_depth: 2 toc_float: collapsed: false smooth_scroll: true code_folding: hide --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warings = FALSE) library(pracma) library(cubature) library(plotly) library(ggplot2) library(doParallel) library(TCIU) library(EBImage) library(brainR) library(spatstat) # number of cores for parallel computing ncor = 8 ``` The [Laplace Transform](https://en.wikipedia.org/wiki/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). # *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 f(t)=t ILT(f, t); t f <- 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) ``` For 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. # 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](https://en.wikipedia.org/wiki/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. ```{r, warning=F, message=F, error=F} 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) } ``` ```{r pressure, warning=F, message=F, error=F} # 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'))) fig ``` We 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}$. ```{r stereographic2D, message=F, error=F, warning=F} # 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))) ``` # Laplace Transform Applications ## ILT of Kimesurfaces $(F(z))$ to Time-series $(f(t))$ ### 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. ```{r warning=FALSE, error=FALSE, message=FALSE} ### 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)`. ```{r warning=FALSE, error=FALSE, message=FALSE} # 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 p ``` We 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. ```{r warning=FALSE, error=FALSE, message=FALSE} 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 p ``` #### 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\}.$$ ```{r regularizedKimesurface1, warning=FALSE, error=FALSE, message=FALSE} # 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 p ``` #### 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). ```{r warning=FALSE, error=FALSE, message=FALSE} ##### 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)")) ``` ### Example 2: $F(z)=|z|$ Another example is finding the time-series Laplace dual to the kimesurface $$F(z)=|z| .$$ ```{r warning=F, error=F, message=F} ### 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))) ``` ### 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. ```{r warning=FALSE, error=FALSE, message=FALSE} ###### 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 ``` ```{r warning=FALSE, error=FALSE, message=FALSE} ##### 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")) ``` ## 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: ### 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)$. ```{r warning=FALSE, error=FALSE, message=FALSE} ####################################################################### ###### 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")) ``` ## 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](https://www.socr.umich.edu/TCIU/HTMLs/Chapter2_ComplexWavefunctions_ComplexTime.html#1_Manifold_Foliation). ```{r LT_Foliation, warning=FALSE, error=FALSE, message=FALSE} ####################################################################### ###### ILT on discrete Laplace transform of sine (analytical form) #### ####################################################################### f_sin <- function(t) { sin(t) } f_sinNoise <- function(t) { sin(t) + rnorm(1, mean=0, sd=0.1) } 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 the 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")) ``` # 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]$. ```{r warning=FALSE, error=FALSE, message=FALSE} # 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")) ``` # 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. ```{r warning=FALSE, error=FALSE, message=FALSE} # 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)$. # 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)$. ```{r eval=TRUE, echo=TRUE, warning=FALSE, error=FALSE, message=FALSE} # 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")) ``` # 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. ## Load the fMRI data ```{r eval=TRUE, warning=FALSE, error=FALSE, message=FALSE} 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")) ``` ## 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. ```{r eval=TRUE, warning=FALSE, error=FALSE, message=FALSE} 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 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 p ``` ## 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. ```{r regularizedFMRIKimesurface, warning=FALSE, error=FALSE, message=FALSE} # 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 p ``` ## 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. ```{r eval=TRUE, echo=TRUE, warning=FALSE, error=FALSE, message=FALSE} ##### 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](https://spacekime.org/), 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.
SOCR Resource Visitor number Web Analytics SOCR Email