obtaining_hurst_exponent_execution_period

library(readxl)
library(knitr)
library(dplyr)
library(ggplot2)
library(HKprocess)
library(parallel)
library(openxlsx)

setwd("C:/Users/Marcel/Desktop/R Scripts Thesis/Sequence Execution Period")
fractal_df <-read_excel("movementperiod_go_until_last_step+65_merged.xlsx")
fractal_df$Block <- factor(fractal_df$Block, levels = as.character(1:10))
fractal_df$Subject <- as.factor(fractal_df$Subject)
fractal_df$Trial <- factor(fractal_df$Trial, levels = as.character(1:48))
fractal_df$Sequence <- as.factor(fractal_df$Sequence)
str(fractal_df)

# Calculate the minimum trial length
trial_lengths <- fractal_df %>% group_by(Subject, Block, Trial) %>% summarise(Trial_Length = n())
min_length <- min(trial_lengths$Trial_Length)

# Filter each trial to the minimum length
fractal_df <- fractal_df %>% group_by(Subject, Block, Trial) %>% slice(1:min_length)

# Check trial lengths again
trial_lengths <- fractal_df %>% group_by(Subject, Block, Trial) %>% summarise(Trial_Length = n())

# Splitting the dataset into the learning blocks and the 9th and 10th block seperately to seperate the computation load into more manageable chunks
Sequence1 <- fractal_df %>% 
  filter(Sequence == 1)

Sequence2 <- fractal_df %>%
  filter(Sequence == 2)

# Function to calculate Hurst exponent for each time series
calculate_hurst_exponent_hk <- function(time_series, sample_size = 50) {
    hurst_samples <- inferH(time_series, n = sample_size)
    hurst_estimate <- median(hurst_samples)
    return(hurst_estimate)
}

# Define the columns you want to iterate over
columns_to_analyze <- c("COM_acceleration_X", "COM_acceleration_Y", "COM_acceleration_Z")

# Set RNG kind suitable for parallel computations. It ensures that the RNG can produce multiple independent streams of random numbers.
RNGkind("L'Ecuyer-CMRG")

# Parallel Processing Setup
num_cores <- detectCores() - 1  # Use one less core than available
cl <- makeCluster(num_cores)

# Initialize different Random Number Generator streams on each worker node. The iseed parameter ensures reproducibility by making the same independent RNG streams get allocated to the individual worker nodes. Then each time the inferH method is called within a worker node, it uses the current RNG state, generates the necessary numbers, and then advances the RNG state.This makes sure that every call to calculate the hurst exponent (using the inferH function) uses a different set of random numbers.
clusterSetRNGStream(cl, iseed = 12345)

# Load Necessary Libraries on Worker Nodes
clusterEvalQ(cl, {
  library(dplyr)
  library(HKprocess)
})

# Export Necessary Variables and Functions to the Cluster
clusterExport(cl, varlist = c("Sequence1", "Sequence2", "calculate_hurst_exponent_hk", "columns_to_analyze"))

# Function to Calculate Hurst Exponents in Parallel for Each Group
parallel_hurst_calculation <- function(sequence, time_series_column) {
  # Split the Data into Groups Manually
  grouped_data <- sequence %>% group_by(Subject, Block, Trial) %>% group_split()
  
  results <- parLapply(cl, grouped_data, function(group) {
    hurst_value <- calculate_hurst_exponent_hk(group[[time_series_column]])
    result <- data.frame(Subject = unique(group$Subject),
                         Block = unique(group$Block),
                         Trial = unique(group$Trial),
                         Hurst = hurst_value)
    return(result)
  })
  return(do.call(rbind, results))
}

# Perform Parallel Computation
results_parallel <- parallel_hurst_calculation(Sequence1, columns_to_analyze[1])
write.xlsx(results_parallel, file = "ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_X.xlsx")

results_parallel2 <- parallel_hurst_calculation(Sequence1, columns_to_analyze[2])
write.xlsx(results_parallel2, file = "ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_Y.xlsx")

results_parallel3 <- parallel_hurst_calculation(Sequence1, columns_to_analyze[3])
write.xlsx(results_parallel3, file = "ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_Z.xlsx")

results_parallel <- parallel_hurst_calculation(Sequence2, columns_to_analyze[1])
write.xlsx(results_parallel, file = "ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_X.xlsx")

results_parallel2 <- parallel_hurst_calculation(Sequence2, columns_to_analyze[2])
write.xlsx(results_parallel2, file = "ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_Y.xlsx")

results_parallel3 <- parallel_hurst_calculation(Sequence2, columns_to_analyze[3])
write.xlsx(results_parallel3, file = "ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_Z.xlsx")

# Stop the Cluster
stopCluster(cl)

# Set the working directory to the location of your Excel files
setwd("C:/Users/Marcel/Desktop/R Scripts Thesis/Sequence Execution Period")

# Import DF's for Sequence 1 and add 'Sequence' column, rename 'Hurst'
seq1_acceleration_x <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_X.xlsx") %>%
  mutate(Sequence = 1) %>%
  rename(hurst_acceleration_x = Hurst)
seq1_acceleration_y <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_Y.xlsx") %>%
  mutate(Sequence = 1) %>%
  rename(hurst_acceleration_y = Hurst)
seq1_acceleration_z <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence1_COM_acceleration_Z.xlsx") %>%
  mutate(Sequence = 1) %>%
  rename(hurst_acceleration_z = Hurst)


# Import DF's for Sequence 2 and add 'Sequence' column, rename 'Hurst'
seq2_acceleration_x <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_X.xlsx") %>%
  mutate(Sequence = 2) %>%
  rename(hurst_acceleration_x = Hurst)
seq2_acceleration_y <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_Y.xlsx") %>%
  mutate(Sequence = 2) %>%
  rename(hurst_acceleration_y = Hurst)
seq2_acceleration_z <- read_excel("ExecutionPeriod_Hurst_Exponents_Sequence2_COM_acceleration_Z.xlsx") %>%
  mutate(Sequence = 2) %>%
  rename(hurst_acceleration_z = Hurst)

# Combine Sequence DF's
acceleration_x_combined <- bind_rows(seq1_acceleration_x, seq2_acceleration_x)
acceleration_y_combined <- bind_rows(seq1_acceleration_y, seq2_acceleration_y)
acceleration_z_combined <- bind_rows(seq1_acceleration_z, seq2_acceleration_z)

# Merge all combined dataframes into one, starting with acceleration_x_combined
combined_data <- acceleration_x_combined %>%
  full_join(acceleration_y_combined, by = c("Subject", "Block", "Trial", "Sequence")) %>%
  full_join(acceleration_z_combined, by = c("Subject", "Block", "Trial", "Sequence"))

# Reordering columns so that 'Sequence' is right before 'hurst_acceleration_x'
combined_data <- combined_data %>%
  select(Subject, Block, Trial, Sequence, 
         hurst_acceleration_x, hurst_acceleration_y, hurst_acceleration_z)

combined_data$Block <- factor(combined_data$Block, levels = as.character(1:10))
combined_data$Subject <- as.factor(combined_data$Subject)
combined_data$Trial <- factor(combined_data$Trial, levels = as.character(1:48))
combined_data$Sequence <- as.factor(combined_data$Sequence)

# Import RT data to append relevant collumns into final dataframe for movement preparation period for subsequent analysis
rt_data  <-read_excel("C:/Users/Marcel/Desktop/R Scripts Thesis/RT Data/rt_data_cleaned_merged.xlsx")
rt_data$Block <- factor(rt_data$Block, levels = as.character(1:10))
rt_data$Subject <- as.factor(rt_data$Subject)
rt_data$Trial <- factor(rt_data$Trial, levels = as.character(1:48))
rt_data$Sequence <- as.factor(rt_data$Sequence)
str(rt_data)

# Calculate mean and total RT per trial and grab the trial accuracy
rt_stats <- rt_data %>%
  group_by(Subject, Block, Trial) %>%
  summarise(trial_mean_rt = mean(RT, na.rm = TRUE),
            trial_total_rt = sum(RT, na.rm = TRUE),
            trial_acc = max(Trial.ACC)) %>%
  ungroup()

# Merge combined_data with the rt_stats
combined_data <- combined_data %>%
  left_join(rt_stats, by = c("Subject", "Block", "Trial"))

write.xlsx(combined_data, "trial_data_execution.xlsx")