--- title: "Laplace Transform and Kimesurface Transform of TCIU Analytics" author: "

SOCR Team

" date: "`r format(Sys.time(),'%m/%d/%Y')`" 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: show vignette: > %\VignetteIndexEntry{Laplace Transform and Kimesurface Transform of TCIU Analytics} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- # Background The Laplace transform 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). In this vignette, we introduce the use of inv_kimesurface_transform() and kimesurface_transform() functions in our package. We demonstrate how to do Laplace transform and kimesurface transform with our functions. ```{r warning=FALSE, message=FALSE} require(TCIU) require(doParallel) require(cubature) require(oro.nifti) require(magrittr) require(plotly) require(ggplot2) ``` # Laplace Transform (LT) and inverse Laplace Transform (ILT) ## discrete LT and analytical form of LT Here we first apply the discrete Laplace Transform (LT) function in our package on the sine function with the domain of $[0 : 2\pi]$ to see whether it has the same 2D surface function as the analytic form of LT of sine. The analytic form is $1/(1+z^2)$. ```{r eval=FALSE, fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9} # For this part of code, we comment it out and import the output plot already generated before to reduce time # But this part of code can be run successfully. If you are interested, you can try it on your computer! # discrete Laplace Transform of sine range_limit = 2 x2 = seq(from = 0, to = range_limit, length.out = 50)[2:50] # drop the first row to avoid real part value of 0 y2 = seq(from = 0, to = range_limit, length.out = 50)[2:50] # drop the first column 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 ncors = 2 # please choose the ncors according to the number of cores your PC has # it is better that you increase the number of cores used for parallel computing if your computer allows cl <- makeCluster(ncors) registerDoParallel(cl) F = list() for (i in 1:length(x2) ){ F[[i]] = foreach(j = 1:length(y2), .export='cubintegrate', .packages='cubature') %dopar% { TCIU::LT(FUNCT=f_sin, complex(real=x2[i], imaginary = y2[j])) } } stopCluster(cl) F_vec = lapply(F, unlist) z2_grid = unlist(do.call(rbind, F_vec)) # explicit form of Laplace Transform of sine laplace_sine = function(p) { 1/(p^2 + 1) } # Exact Laplace transform of sin(x), continuous function XY = expand.grid(X=x2, Y=y2) 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 # make the two plots in the same plot to compare lt_ilt_plotly = plot_ly(hoverinfo="none", showscale = FALSE)%>% add_trace(z=Re(sine_z)-1, type="surface", surfacecolor=Im(sine_z)) %>% add_trace(z = Re(z2_grid), type="surface", opacity=0.7, surfacecolor=Im(z2_grid) )%>% layout(title = "Laplace Transform, LT(sin()), Height=Re(LT(sin())), Color=Re(LT(sin())) \n Contrast Exact (Continuous) vs. Approximate (Discrete) Laplace Transform", showlegend = FALSE) lt_ilt_plotly ``` ```{r fig.align="center", warning=FALSE, message=FALSE} sample_save[[1]] ``` From the plot, we can easily tell that our discrete LT function generate the same 2D surface as the analytic LT function does. ## ILT on discrete LT of sine Here we first apply the discrete Laplace Transform (LT) function in our package on the sine function with the domain of $[0, 2\pi]$, and then use the inverse Laplace Transform (ILT) function in our package to prove the Laplace Transformation has been converted back to sine. ```{r warning=FALSE, message=FALSE, fig.width = 9, fig.align = "center", eval=FALSE} # For this part of code, we comment it out and import the output plot already generated before to reduce time # But this part of code can be run successfully. If you are interested, you can try it on your computer! # discrete Laplace Transform of sine f_sin = function(t) { sin(t) } lt_sine = function(z) TCIU::LT(f_sin, z) # inverse Laplace Transform on the lt_sine # using parallel computing to speed up code tvalsn <- seq(0, pi*2, length.out = 20) cl <- makeCluster(ncors) registerDoParallel(cl) sinvalsn <- foreach(t=1:length(tvalsn), .export='cubintegrate', .packages='cubature') %dopar% { TCIU::ILT(FUNCT=lt_sine, t=tvalsn[t]) } stopCluster(cl) sinvalsn = unlist(sinvalsn) # make the plot of the result from ILT # to see whether it still looks like sine sinvalsn_df2 <- as.data.frame(cbind(Re=Re(sinvalsn),Im=Im(sinvalsn), Sin=sin(tvalsn), time_points=tvalsn)) lt_ilt_sine = 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(sine)")) + 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)) lt_ilt_sine ``` ```{r, fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9} sample_save[[2]] ``` # Kimesurface transform and inverse Kimesurface transform In this part, we still apply the Kimesurface transform on our sine function, and then check whether the inverse Kimesurface transform can convert it back to our original form of sine. Kimesurface transform mainly does a Laplace Transform and then applies a height and surface transformation to have a better view of the plot. ## Kimesurface transform of sine function Below, we apply the Kimesurface transform on the sine function, and make the plot of the 2D function using the magnitude of the complex value as y. ```{r warning=FALSE, message=FALSE, eval=FALSE} x = seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50]; # do Kimesurface transform on sine function z2_grid = kimesurface_transform(FUNCT = function(t) { sin(t) }, real_x = x, img_y = y) # make the plot after Kimesurface transformation surf_color <- atan2(Im(z2_grid), Re(z2_grid)) colorscale = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x))) 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 ``` ```{r fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9} sample_save[[3]] ``` Also, you can try to apply the Kimesurface transformation on our fMRI data for a brain. Each voxel of the brain contains a set of time series data that can be analyzed using our transform. However, due to the number of cores limited for parallel computing in the vignette, we did not run this chuck. You can try it if you are interested. ```{r warning=FALSE, message=FALSE, eval = FALSE, fig.align = "center", fig.width = 9} # load the fMRI data 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 # extract a set of time series data from a voxel of fMRI data xA_fMRI_1D_x20_y20_z11 <- fMRIVolume[20, 20, 11, ] # get the smooth function of this time series data # Instead of using the extremely noisy fMRI data and avoiding integration problems, # smooth "f" and use the **smooth version, f** time_points <- seq(0+0.001, 2*pi, length.out = 180) f <- smooth.spline(ceiling((180*time_points)/(2*pi)), xA_fMRI_1D_x20_y20_z11, df = 10)$y # Define the f(t)=smooth(fMRI)(t) signal as a function of real time 0 2*pi) { return ( 0 ) } else { return ( f[ceiling((180*t)/(2*pi))] ) } } # do the Kimesurface transform ncors = 8 # please choose the ncors according to the number of cores your PC has x = seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50]; # do Kimesurface transform on sine function z2_grid_fmri = kimesurface_transform(FUNCT = fmri_funct, glb_para="f", real_x = x, img_y = y, parallel_computing = TRUE, ncor=ncors) # make the plot of function after Kimesurface transformation surf_color <- atan2(Im(z2_grid_fmri), Re(z2_grid_fmri)) colorscale = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x))) magnitude <- (sqrt( Re(z2_grid_fmri)^2+ Im(z2_grid_fmri)^2)) p_fmri <- 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_fmri ``` ```{r fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9} sample_save[[4]] ``` ## Inverse Kimesurface Transform After the Kimesurface transformation, we apply the inverse Kimesurface transformation, and see that we can get a curve that captures the most trend of the sine function. ```{r warning=FALSE, message=FALSE, fig.width = 9, fig.align = "center", eval=FALSE} time_points <- seq(0+0.001, 2*pi, length.out = 180) inv_data = inv_kimesurface_transform(time_points, z2_grid) inv_data = inv_kimesurface_transform(time_points, z2_grid,num_length = 23, m=1, msg=TRUE) time_Intensities_ILT_df2 <- as.data.frame(cbind(Re=scale(Re(inv_data$Smooth_Reconstruction)), Im=scale(Re(inv_data$Raw_Reconstruction)), fMRI=scale(Re(sin(time_points))), time_points=time_points)) 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") pppp<-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= bquote("Comparison between" ~ "f(t)=Smooth(Sin)(t)" ~ "and Smooth(ILT(LT(Sin)))(t); Range [" ~ 0 ~":"~ 2*pi~"]")) ``` ```{r fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9} sample_save[[5]] ```