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

source("1-read-and-calculate_errors_time.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_ms[1]]) #0-static
techniqueB <- bootstrapMeanCI(data[,speed_duration_ms[2]]) #1-slow
techniqueC <- bootstrapMeanCI(data[,speed_duration_ms[3]]) #2-fast


analysisData <- c()
analysisData$name <- c(speed_duration_ms[3],speed_duration_ms[2],speed_duration_ms[1])
analysisData$pointEstimate <- c(techniqueC[1], techniqueB[1], techniqueA[1])
analysisData$ci.max <- c(techniqueC[3], techniqueB[3], techniqueA[3])
analysisData$ci.min <- c(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 = 3, ymin = 0, ymax = 6, "", "Avg. Error. Error Bars, Bootstrap 95% CIs")
ggsave("plots/error-means.pdf",device = pdf, width=5, height=2)

# CIs with adapted alpha value for multiple comparisons

diffAB = bootstrapMeanCI_corr(data[,speed_duration_ms[1]] - data[,speed_duration_ms[2]], 3) # 0-static - 1-slow
diffCB = bootstrapMeanCI_corr(data[,speed_duration_ms[3]] - data[,speed_duration_ms[2]], 3) # 2-fast - 1-slow
diffCA = bootstrapMeanCI_corr(data[,speed_duration_ms[3]] - data[,speed_duration_ms[1]], 3) # 2-fast - 0-static

analysisData <- c()
analysisData$name <- c(paste(speed_duration_ms[3], speed_duration_ms[1], sep=" - "), # 2-fast - 0-static
                       paste(speed_duration_ms[3], speed_duration_ms[2], sep=" - "), # 2-fast - 1-slow
                       paste(speed_duration_ms[1], speed_duration_ms[2], sep=" - ")) # 0-static - 1-slow
analysisData$pointEstimate <- c(diffCA[1], diffCB[1], diffAB[1])
analysisData$ci.max <- c(diffCA[3], diffCB[3], diffAB[3])
analysisData$ci.min <- c(diffCA[2], diffCB[2], diffAB[2])
analysisData$level <- c(diffCA[4], diffCB[4], diffAB[4])
analysisData$ci_corr.max <- c(diffCA[6], diffCB[6], diffAB[6])
analysisData$ci_corr.min <- c(diffCA[5], diffCB[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 = 3, ymin = -5, ymax = 6.5, "", "")
ggsave("plots/error-diffs.pdf",device = pdf, width=5, height=2)