---
title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)"
subtitle: "Structured Big Data Analytics Case-Study (UKBB)"
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: hide
  # word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
```

```{r,warning=FALSE,message=FALSE,echo=FALSE}
library(TCIU)
library(ggplot2)
library(dplyr)
library(doParallel)
library(plotly)
library(rpart)
library(caret)
library(rattle)
library(mice)
library(foreach)
```

In this case study, we will look at another interesting example of a large structured tabular dataset [UK Biobank (UKBB) data](http://www.ukbiobank.ac.uk/data-showcase). The goal of this study is to compare the data representations as well as the subsequent data analytics between classical time-series and kime-series.

# Data Preparation
A previous investigation [Predictive Big Data Analytics using the UK Biobank Data](http://doi.org/10.1038/s41598-019-41634-y) based on $7,614$ clinical information, phenotypic features, and neuroimaging data of $9,914$ UKBB subjects, reported the twenty most salient derived imaging biomarkers associated with mental health conditions (e.g., depression, anxiety). By jointly representing and modeling the significant clinical and demographic variables along with the derived salient neuroimaging features, the researchers predicted the presence and progression of depression and other mental health disorders in the cohort of UKBB participating volunteers. We will explore the effects of kime-direction on the findings based on the same data and similar analytical methods. To streamline the demonstration, enable efficient calculations, and facilitate direct interpretation, we will transform the data into a tight computable object of dimensions $9,914\times 106$ (participants × features) that can be processed, interpreted and visualized more intuitively. An interactive demonstration of this tensor data showing linear and non-linear dimensionality reduction is available [online](https://socr.umich.edu/HTML5/SOCR_TensorBoard_UKBB).

```{r, warning=FALSE, message=FALSE, echo=FALSE}
UKBB_data <- get(load("UKBB_data_cluster_label.Rdata"))
UKBB_Colnames <- colnames(UKBB_data)

# Extract 106 derived columns from dimensionality reduction process.
# 1. top-50 derived NI biomarkers 
load("NI_Biomarkers_namelist.RData")
top50_NI_Biomarkers <- NI_Biomarkers_namelist

# 2. extract main clinical features 
## binary
top25_BinaryClinical_Biomarkers <- c("X1200.0.0", "X1200.2.0", "X1170.0.0", "X1190.2.0", "X1170.2.0", "X2080.0.0", "X6138.2.2", "X20117.0.0", "X6138.0.2", "X2877.0.0", "X20117.2.0", "X2877.2.0", "X1190.0.0", "X4968.2.0", "X1249.2.0", "X1190.1.0", "X1170.1.0", "X2080.2.0", "X4292.2.0", "X2050.0.0", "X1628.0.0", "X1200.1.0", "X20018.2.0", "X4292.0.0", "X3446.0.0")
## polytomous
top31_PolytomousClinical_Biomarkers <- c("X31.0.0", "X22001.0.0", "X1950.0.0", "X1950.2.0", "X1980.0.0", "X2040.2.0", "X1980.2.0", "X2030.0.0", "X2090.0.0", "X2040.0.0", "X1618.2.0", "X1618.0.0", "X1210.0.0", "X2030.2.0", "X2000.0.0", "X1930.0.0", "X2090.2.0", "X2000.2.0", "X1210.2.0", "X1618.1.0", "X4653.2.0", "X1970.2.0", "X1970.0.0", "X1980.1.0", "X1930.2.0", "X4598.2.0", "X4598.0.0", "X4653.0.0", "X2090.1.0", "X2040.1.0", "X4631.2.0")

# 3. Construct the Computable data object including all salient predictors
ColNameList <- c(top50_NI_Biomarkers, top25_BinaryClinical_Biomarkers,top31_PolytomousClinical_Biomarkers)
col.index <- which(colnames(UKBB_data) %in% ColNameList)
tight106_UKBB_data <- UKBB_data[ , col.index]
dim(tight106_UKBB_data)
```

The UKBB archive contains incomplete records. We employed multiple imputation using chained equations (MICE) to obtain complete instances of the data and avoid issues with missing observations. 
```{r, warning=FALSE, message=FALSE, eval = FALSE}
# This chunk is not set up to run due to long runtime; instead we load the pretrained imputation results.
# MICE imputation (With parallel core computing)
number_cores <- detectCores() - 2
clustUKBB <- makeCluster(number_cores)
clusterSetRNGStream(clustUKBB, 1234)
registerDoParallel(clustUKBB)

imp_tight106_UKBB_data <-
  foreach(no = 1:number_cores, 
          .combine = ibind, 
          .export = "tight106_UKBB_data",
          .packages = "mice") %dopar% {mice(tight106_UKBB_data, m=2,maxit=3, printFlag=T, seed=1234, method = 'cart')}

comp_imp_tight106_UKBB_data <- as.matrix(complete(imp_tight106_UKBB_data), 
                                         dimnames = list(NULL, colnames(tight106_UKBB_data)))
stopCluster(clustUKBB)
save(comp_imp_tight106_UKBB_data, file = "origin_after_imputation.RData")
```


# Modeling Structure
To assess the quality of kime-series based analytical inference, we will focus on predicting a specific clinical phenotype – Ever depressed for a whole week 1, i.e. column "X4598.2.0". A number of model-based and model-free methods can be used for supervised and unsupervised, retrodictive and predictive strategies for binary, categorical, or continuous regression, clustering and classification. We will demonstrate a Random Forest classification approach for forecasting this binary clinical depression phenotype. 

Other preprocessing steps were necessary prior to the main data analytics. To construct a kime-series, we randomly splitted the UKBB archive into 11 time epochs. We used all 11 epochs for phase estimation, and we iteratively used 10 epochs for modeling training and 1 epoch for model testing. Many alternative kime-direction aggregator functions can be employed, for this demonstration, we used nil-phase, averaging and smoothing splines with different degrees of freedom separately to estimate the unknown phase-angles.

Note that in later analysis, we found that the first 50 columns of neuroimaging data do not contribute much to this specific classification task. So in the following analysis, the building blocks of our model contain only elements from the last 56 clinical features (54 of them as predictors, 1 as response variable, 1 discarded due to close relationship with the response variable). 

```{r, warning=FALSE, message=FALSE}
load("origin_after_imputation.RData")
# randomly splitted the UKBB archive into 11 time epochs
set.seed(1234)
suffle_rid <- sample(1:nrow(comp_imp_tight106_UKBB_data))
comp_imp_tight106_UKBB_data <- comp_imp_tight106_UKBB_data[suffle_rid[1:9900],]
col_name <- colnames(comp_imp_tight106_UKBB_data)
dim(comp_imp_tight106_UKBB_data) <- c(11, 900, dim(comp_imp_tight106_UKBB_data)[2])
dim(comp_imp_tight106_UKBB_data)
```

# Fourier Transform
We implemented Fourier Transform to shift the 54 predictors from time domain to frequency domain. 
```{r, warning=FALSE, message=FALSE}
# fft of all 54 predictors
set.seed(1234)
FT_epochs_tight106_UKBB <- array(complex(), c(11, 900, 54))
mag_FT_epochs_tight106_UKBB <- array(NA, c(11, 900, 54))
phase_FT_epochs_tight106_UKBB <- array(NA, c(11, 900, 54))

for (i in 1:11) {
  FT_epochs_tight106_UKBB[i, , ] <- fft(comp_imp_tight106_UKBB_data[i, , c(51:94,97:106)])
  # col96 <X4598.2.0>: response variable; col95 <X4598.0.0>: highly related to the response variable
  # fourier transform
  X2 <- FT_epochs_tight106_UKBB[i, , ]
  mag_FT_epochs_tight106_UKBB[i, , ] <- sqrt(Re(X2)^2+Im(X2)^2); 
  phase_FT_epochs_tight106_UKBB[i, , ] <- atan2(Im(X2), Re(X2)); 
}
dim(phase_FT_epochs_tight106_UKBB)
```

# Phase Estimation
## Nil Phase Estimation
We applied Nil-Phase estimation as a trivial baseline. The idea is to estimate phase with 0 at each position of the fft transformed data. We'll show its model performance compared to true-phase and other phase estimation methods in the next section.
```{r, warning=FALSE, message=FALSE}
# Nil Phase
set.seed(1234)
ift_NilPhase <- array(NA, c(11, 900, 54))
for(i in 1:11){
  Real_Nil <- mag_FT_epochs_tight106_UKBB[i, , ] * cos(0)
  Imaginary_Nil <- mag_FT_epochs_tight106_UKBB[i, , ] * sin(0)
  ift_NilPhase[i,,] <- Re(fft(Real_Nil+1i*Imaginary_Nil, inverse = T)/length(mag_FT_epochs_tight106_UKBB[i,,]))
}
```

## Average Phase Estimation
By this method, we average over the 11 estimated phase for each position of the $900\times 54$ matrix.
```{r,warning=FALSE, message=FALSE}
# Average Phase
avgPhase_FT_epochs_tight106_UKBB <- apply(phase_FT_epochs_tight106_UKBB, c(2,3), mean)
set.seed(1234)
ift_AvgPhase <- array(NA, c(11, 900, 54))
for(i in 1:11){
  Real_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * cos(avgPhase_FT_epochs_tight106_UKBB)
  Imaginary_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * sin(avgPhase_FT_epochs_tight106_UKBB)
  ift_AvgPhase[i,,] <- Re(fft(Real_Avg+1i*Imaginary_Avg, inverse = T)/length(mag_FT_epochs_tight106_UKBB[i,,]))
}
```

## Smoothing Spline Phase Estimation
We fit a smoothing spline to the 11 estimated phase values for each position of the $900\times 54$ matrix. We'll see later that different degrees of freedom of the smoothing spline function will result in different prediction accuracy. The average phase and nil phase estimation can actually be seen as special cases of the smoothing spline estimation with degrees of freedom df = 1 and df = 0 respectively.

We visualize this estimation using one biomarker column "X31.0.0" as an example:
```{r,warning=FALSE,message=FALSE}
plot(1:11,phase_FT_epochs_tight106_UKBB[,10,1],'b',main = "phase_FT_epochs_tight106_UKBB[,10,1]",ylab = "fft phase",xlab = "epochs 1-11",lty = 1,lwd = 1)
lines(1:11,smooth.spline(1:11,phase_FT_epochs_tight106_UKBB[,10,1],df = 10)$y,col = 'red')
lines(1:11,smooth.spline(1:11,phase_FT_epochs_tight106_UKBB[,10,1],df = 5)$y,col = 'blue')
legend("topleft",legend = c("df = 10","df = 5"),lty = 1:1,col = c("red","blue"),cex = 0.8)
```

# Inverse Fourier Transform
After getting the phase estimation, we need to reverse back into the time domain to conduct further classification tasks. The ability to reflect the structure of the original dataset after inverse fft represents the adequacy of the phase estimation process.

We show the results of smoothing spline phase estimation with df = 10 as an example:

```{r,warning=FALSE,message=FALSE}
# inverse FFT example - smoothing spline df = 10
set.seed(1234)
ift_splinePhase <- array(NA, c(11, 900, 54))
smooth_spline_value <- array(NA, c(11, 900, 54))
df <- 10

for(i in 1:900){
  for(j in 1:54){
    smooth_spline_value[,i,j] <- smooth.spline(1:11,phase_FT_epochs_tight106_UKBB[,i,j],df = df)$y
  }
}

for(i in 1:11){
  Real_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * cos(smooth_spline_value[i, , ])
  Imaginary_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * sin(smooth_spline_value[i, , ])
  ift_splinePhase[i,,] <- Re(fft(Real_Avg+1i*Imaginary_Avg, inverse = T)/length(mag_FT_epochs_tight106_UKBB[i,,]))
}
```

*Original Data*
```{r,warning=FALSE,message=FALSE}
head(comp_imp_tight106_UKBB_data[1, ,52],100)
```

*Raw Inverse FFT*
```{r,warning=FALSE,message=FALSE}
head(ift_splinePhase[1,,2],100)
```

*Round Inverse FFT*
```{r,warning=FALSE,message=FALSE}
head(round(ift_splinePhase[1,,2]),100)
```

# Inverse FFT based Random Forest
Then, we compared the model performance of a classification task based on original time-series data and transferred data utilizing a different phase estimation approach.

## Original Data based Random Forest (Unknown True-Phase)
We first build the model based on original data. For comparison across different phase estimation methods, we use the 4th epoch as the test set and the other 10 epochs as the training set for all methods. We have also implemented cross-validation to check performance with each 11 epoch as the testing set and found that the model performance remains similar. So for simplicity, we will only show the result of one train-test epoch setting.

```{r,warning=FALSE,message=FALSE}
library(randomForest)
rf_train_data_ift <- c()
label_train <- c()
test_epoch_idx <- 4
for(i in setdiff(1:11,test_epoch_idx)){
  rf_train_data_ift <- rbind(rf_train_data_ift,comp_imp_tight106_UKBB_data[i,,c(51:94,97:106)])
  label_train <- c(label_train,comp_imp_tight106_UKBB_data[i, ,96])
}
rf_train_data_ift <- as.data.frame(rf_train_data_ift)
rf_train_label_ift <- as.factor(label_train)

set.seed(1234)
rf_Phase_ift_original_data <- randomForest(rf_train_label_ift ~ . ,data= rf_train_data_ift)
rf_Phase_ift_original_data
pred_train = predict(rf_Phase_ift_original_data, type="class")
confusionMatrix(pred_train, rf_train_label_ift)
```

## Smoothing Spline Phase Estimation based Random Forest 
Then we show the results for Random Forest based on transformed data using smoothing spline phase estimation with df = 10 as an example. Later we will show the training and testing accuracy trending with degrees of freedom increasing from 2 to 11. We can see that smoothing spline phase estimation with df = 10 can give us a competitive classification accuracy compared to true-phase Random Forest.

```{r,warning=FALSE,message=FALSE}
rf_train_data_ift <- c()
label_train <- c()
test_epoch_idx <- 4
for(i in setdiff(1:11,test_epoch_idx)){
  rf_train_data_ift <- rbind(rf_train_data_ift,ift_splinePhase[i,,])
  label_train <- c(label_train,comp_imp_tight106_UKBB_data[i, ,96])
}
rf_train_data_ift <- as.data.frame(rf_train_data_ift)
rf_train_label_ift <- as.factor(label_train)

set.seed(1234)
rf_Phase_ift_spline <- randomForest(rf_train_label_ift ~ . ,data= rf_train_data_ift)
rf_Phase_ift_spline
pred_train = predict(rf_Phase_ift_spline, type="class")
confusionMatrix(pred_train, rf_train_label_ift)
```

## Average Phase Estimation based Random Forest
We will see that with average phase estimation, our model has bad classification performance. This also indicates the importance of phase information.

```{r,warning=FALSE,message=FALSE}
rf_train_data_ift <- c()
label_train <- c()
test_epoch_idx <- 4
for(i in setdiff(1:11,test_epoch_idx)){
  rf_train_data_ift <- rbind(rf_train_data_ift,ift_AvgPhase[i,,])
  label_train <- c(label_train,comp_imp_tight106_UKBB_data[i, ,96])
}
rf_train_data_ift <- as.data.frame(rf_train_data_ift)
rf_train_label_ift <- as.factor(label_train)

set.seed(1234)
rf_Phase_ift_avg <- randomForest(rf_train_label_ift ~ . ,data= rf_train_data_ift)
rf_Phase_ift_avg
pred_train = predict(rf_Phase_ift_avg, type="class")
confusionMatrix(pred_train, rf_train_label_ift)
```

## Prediction Accuracy Trend w.r.t. Degrees of Freedom in Phase Estimation
Finally, we visualize the training and testing accuracy of Inverse FFT-based Random Forest with different degrees of freedom in the smoothing spline phase estimation process. Note that we include the Nil-Phase estimation and the Average-Phase estimation as special cases with df = 0 and df = 1 respectively. The dotted horizontal line in the plot represents the testing accuracy of original data based Random Forest. This plot provides strong evidence of the importance of correctly estimating kime-phases.

```{r, warning=FALSE,message=FALSE}
# helper function 1: compute smoothing spline with different degrees of freedom
ift_df <- function(df){
  set.seed(1234)
  ift_splinePhase <- array(NA, c(11, 900, 54))
  smooth_spline_value <- array(NA, c(11, 900, 54))
  for(i in 1:900){
    for(j in 1:54){
      smooth_spline_value[,i,j] <- smooth.spline(1:11,phase_FT_epochs_tight106_UKBB[,i,j],df = df)$y
    }
  }
  for(i in 1:11){
    Real_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * cos(smooth_spline_value[i, , ])
    Imaginary_Avg <- mag_FT_epochs_tight106_UKBB[i, , ] * sin(smooth_spline_value[i, , ])
    ift_splinePhase[i,,] <- Re(fft(Real_Avg+1i*Imaginary_Avg, inverse = T)/length(mag_FT_epochs_tight106_UKBB[i,,]))
  }
  
  return(ift_splinePhase)
}

# helper function 2: compute model accuracy with different test epoch (mainly developed for the cross-validation part not shown)
rf_model_ift <- function(ift_splinePhase,test_epoch_idx){
  # dim(ift_splinePhase) 11,900,54
  rf_train_data_ift <- c()
  label_train <- c()
  for(i in setdiff(1:11,test_epoch_idx)){
    rf_train_data_ift <- rbind(rf_train_data_ift,ift_splinePhase[i,,])
    label_train <- c(label_train,comp_imp_tight106_UKBB_data[i, ,96])
  }
  rf_train_data_ift <- as.data.frame(rf_train_data_ift)
  rf_train_label_ift <- as.factor(label_train)
  # dim(rf_train_data_ift) 9000, 54
  set.seed(1234)
  rf_Phase_ift <- randomForest(rf_train_label_ift ~ . ,data= rf_train_data_ift)
  # OOB_accuracy <- 1- rf_Phase_ift$err.rate[,1][500] # 500 - number of trees
  
  # rf_Phase_ift 
  pred_train = predict(rf_Phase_ift, type="class")
  train_accuracy = sum(pred_train == rf_train_label_ift)/length(rf_train_label_ift)
  
  # test
  rf_test_data_ift <- ift_splinePhase[test_epoch_idx,,]
  rf_test_label_ift <- as.factor(comp_imp_tight106_UKBB_data[test_epoch_idx, ,96])
  
  pred_test = predict(rf_Phase_ift, rf_test_data_ift, type="class")
  test_accuracy = sum(pred_test == rf_test_label_ift)/length(rf_test_label_ift)

  return(c(train_accuracy,test_accuracy))
}
```

```{r,warning=FALSE,message=FALSE,eval = FALSE}
# This chunk is not setup to run due to long runtime; instead we load the pretrained results.
train_test_accuracy <- c()
for(df in 2:11){
  train_test_accuracy <- rbind(train_test_accuracy,c(rf_model_ift(ift_df(df),4),df = df))
}

train_test_accuracy <- as.data.frame(train_test_accuracy)
names(train_test_accuracy) <- c("train_accuracy","test_accuracy","df")

# add results from nil-phase and average phase 
train_test_accuracy <- rbind(train_test_accuracy,c(rf_model_ift(ift_AvgPhase,4),1)) # avg-phase
train_test_accuracy <- rbind(train_test_accuracy,c(rf_model_ift(ift_NilPhase,4),0)) # nil-phase
origin_acc <- rf_model_ift(comp_imp_tight106_UKBB_data[,,c(51:94,97:106)],4)        # original data
save(train_test_accuracy,origin_acc, file = "train_test_accuracy_testepoch4.RData" )
```

```{r, warning=FALSE,message=FALSE}
library(ggplot2)
load("train_test_accuracy_testepoch4.RData")
accuracy_data <- reshape(train_test_accuracy, idvar = "df", varying = list(1:2),v.names = "accuracy", direction = "long") 

names(accuracy_data) <- c("df","group","accuracy")
accuracy_data$group <- factor(accuracy_data$group,levels = c(1,2),labels = c("train","test"))
ggplot(data = accuracy_data, aes(x= df, y = accuracy, group = group)) +
  geom_line(aes(color=group))+
  geom_point(aes(shape=group,color=group)) + 
  scale_x_continuous(breaks = accuracy_data$df) + 
  geom_text(aes(label=round(accuracy,2)),hjust=0,vjust=0.2)+
  ggtitle("Train & Test Accuracy of Inverse FFT-based RF with Different Degree of Freedom") + 
  geom_hline(yintercept = origin_acc[2], linetype="dashed")  
```

# Kimesurface Transformation
We further compare the characteristics of two groups of people with or without suffering from depression in the kime space. For better visualization, we show the kime surface of one feature ("X4631.2.0: Ever unenthusiastic/disinterested for a whole week") that is most correlated with depression as an example. 

```{r, warning=FALSE,message=FALSE}
# data preparation
load("origin_after_imputation.RData")
set.seed(1234)
suffle_rid <- sample(1:nrow(comp_imp_tight106_UKBB_data))
comp_imp_tight106_UKBB_data <- as.data.frame(comp_imp_tight106_UKBB_data[suffle_rid[1:9900],])

depressed_data <- filter(comp_imp_tight106_UKBB_data, X4598.2.0 == 1) # 5007, 106
not_depressed_data <- filter(comp_imp_tight106_UKBB_data, X4598.2.0 == 0)

depressed_idx = which(comp_imp_tight106_UKBB_data$X4598.2.0 == 1)
not_depressed_idx = which(comp_imp_tight106_UKBB_data$X4598.2.0 == 0)

load("ift_splinePhase.RData")
dim(ift_splinePhase) <- c(9900,54)
depressed_data_ift = ift_splinePhase[depressed_idx,]    # dim: 5007, 54
not_depressed_data_ift = ift_splinePhase[not_depressed_idx,]  # dim: 4893, 54

# divide into 110 epochs
comp_imp_tight106_UKBB_data_3D <- as.matrix(comp_imp_tight106_UKBB_data)
dim(comp_imp_tight106_UKBB_data_3D) <- c(110,90,106)

ift_splinePhase_3D <- ift_splinePhase
dim(ift_splinePhase_3D) <- c(110,90,54)

# col45: (most salient predictor) X4631.2.0 (Ever unenthusiastic/disinterested for a whole week)
depressed_data_ift_avg_col45 <- c()
not_depressed_data_ift_avg_col45 <- c()
for(i in 1:110){
  depressed_idx_epoch = which(comp_imp_tight106_UKBB_data_3D[i,,96] == 1)
  not_depressed_idx_epoch = which(comp_imp_tight106_UKBB_data_3D[i,,96] == 0)
  depressed_data_ift_avg_col45[i] = mean(ift_splinePhase_3D[i,depressed_idx_epoch,45])
  not_depressed_data_ift_avg_col45[i] = mean(ift_splinePhase_3D[i,not_depressed_idx_epoch,45])
}
```

We apply **NuLT** function which achieves Laplace Transform to two time-series respectively and we visualize two kime-surfaces.
```{r, warning=FALSE,message=FALSE}
source('LT_ILT.R')
x2 = seq(from = 1, to = 11, length.out = 50)
y2 = seq(from = -5, to = 5, length.out = 50)
XY = expand.grid(X = x2, Y = y2)       
complex_xy = mapply(complex, real=XY$X,imaginary=XY$Y)
time_points <- seq(0+0.001, 2*pi, length.out = 110)

Y1 = depressed_data_ift_avg_col45
kime_depressed <- NuLT(time_points, Y1, complex_xy, mirror = FALSE, range = 2*pi,k = 3)
dim(kime_depressed) = c(length(x2), length(y2))

Y2 = not_depressed_data_ift_avg_col45
kime_not_depressed <- NuLT(time_points, Y2, complex_xy, mirror = FALSE, range = 2*pi,k = 3)
dim(kime_not_depressed) = c(length(x2), length(y2))
```


```{r, warning=FALSE,message=FALSE}
kime_grid1 <- kime_depressed
kime_grid2 <- kime_not_depressed
phase_depressed <- atan2(Im(kime_grid1), Re(kime_grid1))
magnitude_depressed <- sqrt( Re(kime_grid1)^2+ Im(kime_grid1)^2)
phase_not_depressed <- atan2(Im(kime_grid2), Re(kime_grid2))
magnitude_not_depressed <- sqrt(Re(kime_grid2)^2+ Im(kime_grid2)^2)

plot_ly(hoverinfo="none", showscale = FALSE)%>%
  add_trace(z=magnitude_depressed, surfacecolor=phase_depressed, type="surface",showlegend=TRUE,name = 'Depressed') %>%
  add_trace(z=magnitude_not_depressed, surfacecolor=phase_not_depressed, type="surface", opacity=0.7,showlegend=TRUE,name = 'Not Depressed') %>%
  
  layout(title = "X4631.2.0 Kimesurface Difference: Depressed vs Not Depressed",
         # X4631.2.0 (Ever unenthusiastic/disinterested for a whole week)
         showlegend = TRUE,
         scene = list(aspectmode = "manual",
                      aspectratio = list(x=1, y=1, z=1.0)  
         )
  )
```

# Inverse Kimesurface Transformation
We further recover the time-series from the kime-surface by applying the **ctILT** to the **NuLT** function output.
```{r, warning=FALSE,message=FALSE}
z2_funct_depressed <- function(z) NuLT(time_points,depressed_data_ift_avg_col45, z, k = 3, mirror = FALSE, fitwarning = FALSE, range = 2*pi)
inv_result_depressed <- ctILT(z2_funct_depressed,nnt = 110,tend = 2*pi)

# plot to validate LT and ILT
plot(time_points, inv_result_depressed, col="blue",
     xlim=c(0,2*pi))
lines(time_points,depressed_data_ift_avg_col45, col="red")
legend("top",c("recovered","data"), fill=c("blue","red"))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted")
```

We plot the difference of recovered and original time-series to validate the LT and ILT process.
```{r, warning=FALSE,message=FALSE}
# diff
plot(time_points,inv_result_depressed - depressed_data_ift_avg_col45,lty = 1,lwd = 1,'l',col = 'blue',ylab = "diff", main = 'Difference of IFT recovery and original time-series')
```

<!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://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-->