---
title: "TCIU: The Laplace Transform of Longitudinal Data: Analytic Duality of Time-series and Kimesurfaces"
output: 
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 5
    toc_float:
      collapsed: false
      smooth_scroll: true
---

<!-- 
output: word_document
-->

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(warn = -1)
require(pracma)
require(cubature)
require(plotly)
library(doParallel)
ncor = 8 # number of cores for parallel computing
```

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): R^+ \longrightarrow C$, the **Laplace transform** is the function of a complex frequency argument, $F(z)={\mathcal {L}}(f)(z):C \longrightarrow 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):R^+ \longrightarrow 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 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 in 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).$$


# *Discrete* Laplace Transforms

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

```{r error=FALSE, message=FALSE, warning=FALSE}
# Numerical Laplace Transform (LT) and its Inverse (ILT)

################## 1. (forward) Laplace Transform
#' @title Laplace Transform
#' @description numerically compute the Laplace Transform
#'
#' @param FUNCT Input function f(t)
#' @param z complex domain value to evaluate the F(z)=LT(FUNCT)(z)
#' 
#' @examples
#' f <- function(t) { t }; z= 1+1i; LT(f, z); F <- function (z) { 1/z^2 }; F(z)
#'
LT <- function(FUNCT, z){
  FUNCT <- match.fun(FUNCT) # input R-domain FUNCT should be interpreted as a function
  # define the integrand function (complex-valued)
  # Note that cubic numeric integration is used (Unified Cubature Integration Interface),
  # Integral limits are exact [0, Inf)
  integrand <- function(t) {   return(FUNCT(t) * exp(-z*t)) }
  return(cubintegrate(f = function(t) Re(integrand(t)), 
                      lower = 0, upper = Inf, method = "pcubature")$integral +
           1i * cubintegrate(f = function(t) Im(integrand(t)), 
                             lower = 0, upper = Inf, method = "pcubature")$integral)
}
```

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

```{r error=FALSE, message=FALSE, warning=FALSE}
### 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
f <- function(t) { exp(-5*t) }  # test function, F(z)=L(f)(z)=1/(z+5)
LT(f, z); 1/(z+5)
```

## *Discrete Inverse* Laplace Transform

### Basic functions

Let's start by defining some basic functions that transform between *Cartesian* and *Polar* coordinates.

```{r error=FALSE, message=FALSE, warning=FALSE}
################ 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 eturns cartesian coordinate y from a pair of polar coordinates (r, phi)
y.rphi <- function(r, phi){
  return(r*sin(phi))
}
```

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

```{r error=FALSE, message=FALSE, warning=FALSE}
################## 3. 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 co-ordinates (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, pp531-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 countour
#' @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 countour
bromwichPathDerivative=function (phi, gamma = 0.5) {
  return((0+1i) * (bromwichContour(phi, gamma)^2 + (gamma * tan(phi)/cos(phi))^2)^0.5)
}
```

### Define the discrete Inverse Laplace Transform (ILT)

Next, we will define the discrete `ILT()` function, which inverts the Laplace Transform, $\mathcal{L}^{-1}$.

```{r error=FALSE, message=FALSE, warning=FALSE}
################## 4. Inverse Laplace Transform #######################################################
#' @title Inverse Laplace Transform
#' @description Numerical inverse of the Laplace Transform (ILT)
#'
#' @param FUNCT the function F(z), typically a Laplace-Transform of a function f(t)
#' @param t time domain value to evaluate the ILT(F)(t)
#' @param nterms number of terms to use in the numerical inversion (odd number)
#' @param m width of the contour path in C; too small values may lead to singularities on the negative x-axis; 
#' too large valued may lead to numerical instability for large positive x-axis
#' @param fail_val value to return in event of failure to converge
#' @param gamma value on the positive x-axis for the vertical line representing the countour
#' @msg Boolean to show/hide warnings
#'
ILT <- function(FUNCT, t, nterms = 31L, m=1, gamma=0.5, fail_val = complex(0+0i), msg=TRUE){
  FUNCT <- match.fun(FUNCT)           # input C-valued FUNCT should be interpreted as a function
  if(identical(t, 0)) { t <- 10^-50 } # accurate for very low values of t>0
  n_attempts <- 10
  
  for (attempt in 1:n_attempts) {  # repeated numerical integration attempts
    dphi <- 2*pi/nterms
    phi <- -pi + (1:nterms - 1/2)*dphi
    z <- optimContour(phi, m, t)*cos(phi) + optimContour(phi, m, t)*sin(phi)*1i
    # Plot optimal contour   # plot(Re(z), Im(z))
    L <- vapply(z, FUNCT, 1i) # Specify the rerurn Function type must be a single C-value!
    # L <- sapply(z, FUNCT)  # Funct returns array or list, not a single C-value!
    
    if(any(!is.finite(L))){
      if(msg==TRUE & identical(attempt, as.integer(n_attempts))) { 
        # if original method fails, construct ILT ob bromwich contour
        sprintf("Laplace Transform inversion failed after %s attempts.\n", n_attempts)
        # Note: Try alternative Bromwich contour in C, which may not be very stable
        nterms = 1000L
        tot <- 0
        dphi <- pi/nterms
        for (n in 1:nterms) {
          phi <- -pi/2 + (n - 1/2) * dphi
          ds <- dphi * bromwichPathDerivative(phi, gamma)
          x <- x.rphi(bromwichContour(phi, gamma), phi)
          y <- y.rphi(bromwichContour(phi, gamma), phi)
          tot <- tot + exp((x + y*(0+1i)) * t) * FUNCT(x + y*(0+1i)) * ds
        }
        return(tot/(2 * pi * (0+1i)))
      }
      m <- m*2   # contour path increment
      nterms <- round(nterms*1.378, 0) # random irregular number to avoid sampling same point again
      next
    }
    em <- exp(z*t)
    ds_dphi <- vapply(phi, analyticalPathDerivative, 1i, m, t)
    break
  }
  return(sum(em*L*ds_dphi)*dphi/(2i*pi))
}
```

Let's test the `ILT()` function using some [simple functions](https://en.wikipedia.org/wiki/Laplace_transform#Table_of_selected_Laplace_transforms). 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}.$$

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

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

<!-- https://www.wolframalpha.com/input/?i=integral+%28sin%28r%29cos%28t-r%29dr+r+from+0+to+t%29+ -->

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}
library(ggplot2)
library(plotly)

### 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: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,
              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",
                            aspectratio = list(x=1, y=1, z=0.5),
                            zaxis = list(range = c(-5,5))
              )
          ) # 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
```

Let's now reconstruct and plot the 1D function $f(t)=\mathcal{L}^{-1}(F)(t)$.

```{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% { 
  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))
```

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



### 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 lost 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) # [1] 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% {
    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))

```


## 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 see 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 functione 
lt_func = function(z) 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% { 
  ILT(FUNCT=lt_func, t=tvalsn[t])
}
stopCluster(cl)
sinvalsn = unlist(sinvalsn)
sinvalsn_df2 <- as.data.frame(cbind(Re=Re(sinvalsn),Im=Im(sinvalsn), Sin=sin(tvalsn), time_points=tvalsn))
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))
```

# ILT on 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 = 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-domail 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 : 2Pi]")
```


# 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-funciton) 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 heigh 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 the 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-domail 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-domail 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 : 2Pi]") +
  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))

```


# fMRI example

Let's demonstrate the Laplace dual of 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}
library(EBImage)
library(brainR)
# library(RNifti)
library(spatstat) 

fMRIURL <- "http://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-corse 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 # smoth 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")

library(ggplot2)
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())
```

## 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 problmes, 
# 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 # smoth 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

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

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 = 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
```

## 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-domail 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% {
    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))# 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 congruency ...
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))
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)=Smooth(fMRI)(t) and Smooth(ILT(LT(fMRI)))(t);\n Range [0 : 2Pi]")

```

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, analyzed to obtain robust prediction or forecasting.


<!--html_preserve-->
<div>
    	<footer><center>
			<a href="http://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img src="http://counter.digits.net/?counter=SOCR"
	 			align="middle" border="0" height="20" hspace="4" vspace="2" width="60">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="http://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="http://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

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