inst/doc/tciu-LT-kimesurface.R

In TCIU: Spacekime Analytics, Time Complexity and Inferential Uncertainty

## ----warning=FALSE, message=FALSE---------------------------------------------
require(TCIU)
require(doParallel)
require(cubature)
require(oro.nifti)
require(magrittr)
require(plotly)
require(ggplot2)

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

## ----fig.align="center", warning=FALSE, message=FALSE-------------------------
sample_save[[1]]

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

## ---- fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9-------
sample_save[[2]]

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

## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9--------
sample_save[[3]]

## ----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<t<=2*pi
#  fmri_funct <- function(t) {
#    if (t < 0+0.001 || t > 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

## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9--------
sample_save[[4]]

## ----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~"]"))

## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9--------
sample_save[[5]]