The Laplace Transform, \(\mathcal{L}\), allows us to examine the relations between the space-time and space-kime representations of longitudinal data. The Fourier transformation is a linear operator that maps complex-valued functions of real variables (e.g., space, time domains) to complex valued functions of other real variables (e.g., frequency domain). The Laplace transform is similar, however, it sends complex-valued functions of positive real variables (e.g., time) to complex-valued functions defined on complex variables (e.g., kime).

1 Continuous Forward and Inverse Laplace Transforms (LT)

The forward and inverse (continuous) Laplace transforms are defined below.

For a given function function (of time) \(f(t): 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).\]

2 Discrete Laplace Transforms

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

# 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}.\]

### 2. Tests the LT
z= 1+1i # Complex-domain value
f <- function(t) { t }  # test function, F(z)=L(f)(z)=1/z^2
LT(f, z); 1/z^2

## [1] 0-0.5i

## [1] 0-0.5i

f <- function(t) { exp(-5*t) }  # test function, F(z)=L(f)(z)=1/(z+5)
LT(f, z); 1/(z+5)

## [1] 0.1621622-0.027027i

## [1] 0.1621622-0.027027i

2.2 Discrete Inverse Laplace Transform

2.2.1 Basic functions

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

################ Cartesian to Polar Coordinate Transform ##############################################
#' Cartesian to Polar Coordinate Transform
#' @title Cartesian to Polar Coordinate Transform
#' @description In r.xy Returns polar coordinate r from a pair of Cartesian coordinates (x, y)
#'
#' @param x x co-ordinate
#' @param y y co-ordinate
#' @return r or phi respectively from x and y
r.xy <- function(x, y){
  return(sqrt(x^2 + y^2))
}

#' @description In r.xy Returns polar coordinate phi from a pair of Cartesian coordinates (x, y).
phi.xy <- function(x, y){
  return(atan2(y, x)) 
}

################ Polar to Cartesian Coordinate Transform ##############################################
#' Polar to Cartesian Coordinate Transform
#' @title Polar to Cartesian Coordinate Transform
#' @description In x.rphi Returns cartesian coordinate x from a pair of polar coordinates (r, phi)
#'
#' @param r distance from origin (Time)
#' @param phi phase
#' @return x or y respectively from r and phi
x.rphi <- function(r, phi){
  return(r*cos(phi))
}

#' @description In y.rphi eturns cartesian coordinate y from a pair of polar coordinates (r, phi)
y.rphi <- function(r, phi){
  return(r*sin(phi))
}

2.2.2 Contour paths in \(C\)

Next we will define alternative contour paths in the complex plane, \(C\), the optimal contour optimContour(), the Bromwich contour, and their corresponding numerical derivatives along contour paths, analyticalPathDerivative() and bromwichPathDerivative() functions.

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

2.2.3 Define the discrete Inverse Laplace Transform (ILT)

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

################## 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. 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}.\]

### 5. Tests the ILT
t = 1/5 # R-domain value
f <- function(z) { 1/(z^2) }  # test function, F(z)=L(f)(z)=1/z^2 ==> f(t)=t
ILT(f, t); t

## [1] 0.2-0i

## [1] 0.2

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)

## [1] 0.3678794-0i

## [1] 0.3678794

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.

3 Laplace Transform Applications

3.1 ILT of Kimesurfaces \((F(z))\) to Time-series \((f(t))\)

3.1.1 Example 1: \(F(z)=F_1(z)+F_2(z)\times F_3(z) +F_4(z)\)

Suppose we want to apply the ILT (\(\mathcal{L}^{-1}\) ) to reconstruct a time-series, \(\hat{f}(t)=\mathcal{L}^{-1}(F)(t)\), corresponding to a given a composite kimesurface: \[F(z)=\mathcal{L}(f)=\underbrace{\frac{1}{z+1} }_{F_1(z)=\mathcal{L}\left( f_1(t)=e^{-t}\right )} + \underbrace{\frac{1}{z^2 + 1}}_{F_2(z)=\mathcal{L}\left( f_2(t)=\sin(t)\right )} \times \underbrace{\frac{z}{z^2 + 1}}_{F_3(z)=\mathcal{L}\left( f_3(t)=\cos(t)\right )} + \underbrace{\frac{1}{z^2}}_{F_4(z)=\mathcal{L}\left( f_4(t)=t\right )} .\]

Mind the decomposition of the compound (time and kime) functions in terms 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.

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)

## [1] 1.16-2.28i

######### Display the actual kimesurface: F(z)
x2 <- seq(from = -2, to = 2, length.out = 200)
y2 <- seq(from = -2, to = 2, length.out = 200)
z2_grid = array(complex(), dim=c(length(x2), length(y2)))

for (i in 1:dim(z2_grid)[1]) {
  for (j in 1:dim(z2_grid)[2]) {
    z2_grid[i,j] = fMRI_Kimesurface(complex(real=x2[i], imaginary = y2[j]))
  }
}

image(Im(z2_grid), axes=FALSE); contour(Im(z2_grid), add=T, col="red", lwd=2, drawlabels=FALSE, axes=FALSE);
contour(Re(z2_grid), add=T, col="blue", lwd=2, drawlabels=F, 
        main = "Im(f) Image with Re(f) and Im(f) Contours in red and blue", axes=F); 
title("Re(f) Image with Re(f) and Im(f) Contours in red and blue", font = 5)

Recall that kimesurfaces are defined over complex-time (Complex-pane) and have complex-value ranges. The plot below shows the original kimesurface as a 2D manifold whose height (vertical dimension) is represented by the Real part of the kimesurface intensity, Re(fMRI_Kimesurface), and the surface color encodes the Imaginary part of the kimesurface intensity, Im(fMRI_Kimesurface).

# Surface height= Re(fMRI_Kimesurface); Surface color = Im(fMRI_Kimesurface)
z3 <- Im(z2_grid)
surf_color <- z3
for (i in 1: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

TCIU: The Laplace Transform of Longitudinal Data: Analytic Duality of Time-series and Kimesurfaces