# filter heatmap colors  ------------------------------------------------------------

filter_heatmap_colors <- function(Labels.data.file, remove.areas) {
  tryCatch({
    flog.debug("Loading data", name = log.name)
    
    remove.areas <- str_replace_all(remove.areas, "_", " ")
    
    Labels.data <- read_delim(file = Labels.data.file, delim = ";")
    #Labels.data$fMRIName <- gsub(pattern = " ", replacement = "_", x = Labels.data$fMRIName)
    
    Labels.data <- Labels.data %>%
      filter(!(fMRIName %in% remove.areas))
    
    return(Labels.data)
    
  }, error = function(e) {
    flog.error("ERROR!", name = log.name)
    flog.error(e, name = log.name)
  })
}


# filter heatmap colors  ------------------------------------------------------------

process_heatmap_colors <- function(Labels.data, Heatmap.Level.color) {
  tryCatch({
    flog.debug("Process heatmap", name = log.name)
    
    # Level color is taken from config file
    Heatmap.brains.sorted <- left_join(Labels.data, 
                                       Heatmap.Level.color, 
                                       by = "Parental_brain_region")  %>% 
      mutate(Parental_brain_region = factor(Parental_brain_region, 
                                            levels = Heatmap.brain.levels)) %>%
      arrange(desc(Hemisphere), 
              Parental_brain_region, Color)
    
    return(Heatmap.brains.sorted)
    
  }, error = function(e) {
    flog.error("ERROR!", name = log.name)
    flog.error(e, name = log.name)
  })
}

# Functions needed for data processing

# Matlab files ------------------------------------------------------------

# load matlab files and run correlations
loadMatCorr <- function(file, raw.data.folder = "./Data/MRI Raw signal data/",
                        Treatment, Impulsivity, Matlab.brain.labels, 
                        remove.areas = c(paste('Cortical_subplate_L'),paste('Cortical_subplate_R'),
                                         paste('Corpus_callosum_L'),paste('Corpus_callosum_R'),
                                         paste('lateral_ventricle_L'),paste('lateral_ventricle_R'),
                                         paste('third_ventricle_L'),paste('third_ventricle_R'),
                                         paste('fourth_ventricle_L'),paste('fourth_ventricle_R'),
                                         paste('white_Matter_L'),paste('white_Matter_R'),
                                         paste('Spinal_cord_L'),paste('Spinal_cord_R'))){
  
  flog.debug(paste("Processing file:", file, sep = " "), name = log.name)
  Data.temp <- readMat(con = paste0(raw.data.folder, "/",
                                    Treatment, "/",
                                    Impulsivity, "/",
                                    file))
  
  Data.temp <- data.frame(Data.temp$SeedSeries)
  
  # change column names to brain atlas names
  names(Data.temp) <- Matlab.brain.labels
  
  # remove brain areas (white matter mainly)
  excl_reg <- -match(remove.areas, names(Data.temp))
  Data.temp <- Data.temp[,excl_reg]
  
  # run correlations
  correlation <- rcorr(as.matrix(Data.temp))$r
  
  Data.gathered <- data.frame(From = rownames(correlation)[row(correlation)[upper.tri(correlation, diag = FALSE)]], 
                              To = colnames(correlation)[col(correlation)[upper.tri(correlation, diag = FALSE)]], 
                              correlation.rho = correlation[upper.tri(correlation, diag = FALSE)])
  
  # set IDs
  Data.gathered$Impulsivity <- Impulsivity
  Data.gathered$File <- file
  Data.gathered$Treatment <- Treatment
  
  return(Data.gathered)
}

# Two step BH correction --------------------------------------------------

two.step.BH <- function(data, method = "TSBH", alpha = 0.05){
  
  # this uses the multtest package to correct for multiple comparisons
  
  stats <- mt.rawp2adjp(rawp = data, proc = method, alpha = alpha)
  
  # the p values are sorted with increasing value, so need to resort to original order
  p.ordered <- stats$adjp[order(stats$index),] 
  
  return(p.ordered[,2]) # the first column is the unadjusted p values
  
}

# HvsL t.test -------------------------------------------------------------
HvsL.tTest <- function(data){

  
  if (all(data$FisherZ == 0)) {
    
    pval <- data.frame(p.value = NA,
                       t.value = NA)
    
  } else {
    
    pval.data <- t.test(FisherZ ~ Impulsivity, data = data, paired = FALSE, var.equal = TRUE)
    
    pval <- data.frame(p.value = pval.data$p.value,
                       t.value = as.numeric(pval.data$statistic))
  }
  
  return(pval)
}


# Sum data for bar plot ---------------------------------------------------

sum_data <- function(Data, single.hemisphere = "TRUE", filter = "p.value < 0.05"){
  
  if (single.hemisphere) {
    Data <- Data %>% 
      mutate(From = gsub(pattern = " L| R", replacement = "", From))
  }
  
  Sum <- Data %>% 
    filter_(filter)  %>%
    mutate(Direction = if_else(t.value < 0, "Negative", "Positive")) %>%
    group_by(Treatment, From, Direction) %>% 
    summarise(N = n()) %>% 
    mutate(N = if_else(Direction == "Negative", -N, N)) %>% ungroup()
  
  return(Sum)
}

# Smirnov test ------------------------------------------------------------

KSmirnov.test <- function(data){
  
  if (all(data$FisherZ == 0)) {
    
    pval <- data.frame(p.value = NA,
                       t.value = NA)
    
  } else {
    
    data <- data %>% select(From, To, Impulsivity, FisherZ) %>% 
      spread(key = Impulsivity, value = FisherZ)
    
    pval.data <- suppressWarnings(ks.test(x = data$high, y = data$low, alternative = "two.sided"))
    
    pval <- data.frame(KS.p.value = pval.data$p.value,
                       KS.D.value = as.numeric(pval.data$statistic))
  }
  
  return(pval)
}

# Two-way RM ANOVA --------------------------------------------------------

TwowayRM_ANOVA <- function(Data){
  model <- aov(FisherZ ~ Impulsivity * Treatment + Error(File/Treatment), data = Data) #>Treatment/File
  summary <- summary(model)

  Impulsivity <- summary$`Error: File`[[1]]["Impulsivity",] %>% mutate(Stat = row.names(.))
  Treatment <- summary$`Error: File:Treatment`[[1]]["Treatment",] %>% mutate(Stat = row.names(.))
  Interaction <- summary$`Error: File:Treatment`[[1]]["Impulsivity:Treatment",] %>% mutate(Stat = row.names(.))
  
  ANOVA.temp <- bind_rows(Impulsivity, 
                          Treatment,
                          Interaction) 

  return(ANOVA.temp)
}
