library(plyr)
library(dplyr)
library(MPDiR)
library(quickpsy)
library(fitdistrplus)
library(ggplot2)

source("1-read-and-calculate_errors_time_1.R") # to read global values


# set directory where script is
sourceDir <- dirname (rstudioapi::getActiveDocumentContext()$path) 
defaultpath <- sourceDir

#remove(list = ls())
print(defaultpath)
setwd(defaultpath)

filename <- "results/aggregated_per_percentage.csv"
errors_table <- read.csv(filename)

filename_analysis <- "results/error_analysis"

### CI analysis

source("CI-Functions-Bonferroni.R")
mydata <- errors_table

# deal wih missing values
# mydata$mean_abs_error[is.na(mydata$mean_abs_error)] <- 0 

# order for the transpose
elements <- mydata
elements <- elements [ order(elements$participant_id, elements$speed_duration), ]
#
myvars <- c("participant_id", "speed_duration", "mean_abs_error")
elements <- elements [myvars]
#aggregating all cases per participant
statstable <- ddply(elements,
                    c("participant_id","speed_duration"),
                    summarise,
                    mean_abs_error=mean(mean_abs_error)
)
elements <- statstable

#
elements <- reshape(elements, timevar="speed_duration", idvar=c("participant_id"), direction="wide")
colnames(elements) <- gsub("mean_abs_error.", "", colnames(elements))

#############
# all elements analysis
#############

data <- elements

techniqueA <- bootstrapMeanCI(data[,speed_duration[1]]) #30-fast-irregular-movement
techniqueB <- bootstrapMeanCI(data[,speed_duration[2]]) #30-fast-linear-movement
techniqueC <- bootstrapMeanCI(data[,speed_duration[3]]) #15-slow-irregular-movement
techniqueD <- bootstrapMeanCI(data[,speed_duration[4]]) #15-slow-linear-movement


analysisData <- c()
analysisData$name <- c(speed_duration[4],speed_duration[3],speed_duration[2],speed_duration[1])
analysisData$pointEstimate <- c(techniqueD[1],techniqueC[1], techniqueB[1], techniqueA[1])
analysisData$ci.max <- c(techniqueD[3],techniqueC[3], techniqueB[3], techniqueA[3])
analysisData$ci.min <- c(techniqueD[2],techniqueC[2], techniqueB[2], techniqueA[2])

datatoprint <- data.frame(factor(analysisData$name),analysisData$pointEstimate, analysisData$ci.min, analysisData$ci.max)
colnames(datatoprint) <- c("speed", "mean_error", "lowerBound_CI", "upperBound_CI ")
tmp <- paste(filename_analysis, "means.csv", collapse="_")
write.csv(datatoprint, file=tmp)

## Code that should print chart, there is some issue
colnames(datatoprint) <- c("technique", "mean_time", "lowerBound_CI", "upperBound_CI ") #We use the name mean_time for the value of the mean even though it's not a time, it's just to parse the data for the plot
barChart(datatoprint,analysisData$name ,nbTechs = 4, ymin = 0, ymax = 6, "", "Avg. Error. Error Bars, Bootstrap 95% CIs")
ggsave("plots/error-means.pdf",device = pdf, width=6, height=2)

# CIs with adapted alpha value for multiple comparisons

diffAB = bootstrapMeanCI_corr(data[,speed_duration[1]] - data[,speed_duration[2]], 4) # 30-fast-irregular-movement - 30-fast-linear-movement
diffCD = bootstrapMeanCI_corr(data[,speed_duration[3]] - data[,speed_duration[4]], 4) # 15-fast-irregular-movement - 15-fast-linear-movement
diffAC = bootstrapMeanCI_corr(data[,speed_duration[1]] - data[,speed_duration[3]], 4) # 30-fast-irregular-movement - 15-fast-irregular-movement
diffBD = bootstrapMeanCI_corr(data[,speed_duration[2]] - data[,speed_duration[4]], 4) # 30-fast-linear-mmovement - 15-fast-linear-movement

analysisData <- c()
analysisData$name <- c(paste(speed_duration[2], speed_duration[4], sep=" - "), # 30-fast-linear-movement - 15-fast-linear-movement
                       paste(speed_duration[1], speed_duration[3], sep=" - "), # 30-fast-irregular-movement - 15-fast-irregular-movement
                       paste(speed_duration[3], speed_duration[4], sep=" - "), # 15-fast-irregular-movement - 15-fast-linear-movement
                       paste(speed_duration[1], speed_duration[2], sep=" - ")) # 30-fast-irregular-movement - 30-fast-linear-movement
analysisData$pointEstimate <- c(diffBD[1],diffAC[1], diffCD[1], diffAB[1])
analysisData$ci.max <- c(diffBD[3],diffAC[3], diffCD[3], diffAB[3])
analysisData$ci.min <- c(diffBD[2],diffAC[2], diffCD[2], diffAB[2])
analysisData$level <- c(diffBD[4],diffAC[4], diffCD[4], diffAB[4])
analysisData$ci_corr.max <- c(diffBD[6],diffAC[6], diffCD[6], diffAB[6])
analysisData$ci_corr.min <- c(diffBD[5],diffAC[5], diffCD[5], diffAB[5])

datatoprint <- data.frame(factor(analysisData$name), analysisData$pointEstimate, analysisData$ci.max, analysisData$ci.min, analysisData$level, analysisData$ci_corr.max, analysisData$ci_corr.min)
colnames(datatoprint) <- c("current_speed", "mean_error", "upperBound_CI","lowerBound_CI", "corrected_CI", "upperBound_CI_corr", "lowerBound_CI_corr") #We use the name mean_time for the value of the mean even though it's not a time, it's just to parse the data for the plot
tmp <- paste(filename_analysis, "diffs.csv", collapse="_")
write.csv(datatoprint,file=tmp)

## Code that should print chart
colnames(datatoprint) <- c("technique", "mean_time", "lowerBound_CI", "upperBound_CI", "corrected_CI","lowerBound_CI_corr","upperBound_CI_corr") #We use the name mean_time for the value of the mean even though it's not a time, it's just to parse the data for the plot
barChart_corr(datatoprint, analysisData$name, nbTechs = 4, ymin = -0.5, ymax = 2.5, "", "")
ggsave("plots/error-diffs.pdf",device = pdf, width=8, height=2)