#This scripts are based on the supplemental material of the paper:
# Peña-Araya, V., Pietriga, E., & Bezerianos, A. (2019). A Comparison of Visualizations for Identifying Correlation over Space and Time. IEEE transactions on visualization and computer graphics.

library(plyr)
library(ggplot2)

source("CI-Functions-Bonferroni.R")

allwm <- read.csv("comparison.csv")
data <- subset(allwm, isTraining==0)
data$feature <- revalue(data$feature, c("Center"="Center", "Source"="Repelling Node", "Sink"="Attracting Node", "SourceVortex"="Repelling Focus", "SinkVortex"="Attracting Focus", "Saddle"="Saddle"))
data_time<-subset(data, trialCorrectCount!=0)



features <- c("Center", "Repelling Node", "Attracting Node", "Repelling Focus", "Attracting Focus", "Saddle")
data_time$meanTime = data_time$time/data_time$trialCorrectCount

summary_time <- summarySEwithin(data_time, measurevar="meanTime", withinvars=c("technique", "feature"), idvar="participant")
ggplot(summary_time, aes(x=feature, y=meanTime, fill=technique)) + geom_bar(position=position_dodge(), stat="identity", colour="gray50") + theme_bw(base_size = 30) + geom_errorbar(aes(ymin=meanTime-ci, ymax=meanTime+ci), width=.2, position=position_dodge(.9)) + ylim(-0.5, 16) + scale_fill_manual(values=cbPalette)

summary_error <- summarySEwithin(data, measurevar="errorCount", withinvars=c("technique", "feature"), idvar="participant")
ggplot(summary_error, aes(x=feature, y=errorCount, fill=technique)) + geom_bar(position=position_dodge(), stat="identity", colour="gray50") + theme_bw(base_size = 30) + geom_errorbar(aes(ymin=errorCount-ci, ymax=errorCount+ci), width=.2, position=position_dodge(.9)) + ylim(-0.5, 2) + scale_fill_manual(values=cbPalette)

summary_error_all <- summarySEwithin(data, measurevar="errorCount", withinvars=c("technique"), idvar="participant")
ggplot(summary_error_all, aes(x=technique, y=errorCount, fill=technique)) + geom_bar(stat="identity", colour="gray50") + theme_bw(base_size = 30) + geom_errorbar(aes(ymin=errorCount-ci, ymax=errorCount+ci), width=.2, position=position_dodge(.9)) + ylim(0, 1) + scale_fill_manual(values=cbPalette)


aggregated_time <- ddply(data_time,
                            c("participant","technique", "feature"),
                            summarise,
                            mean_time = mean(meanTime)
                            
)

# order for the transpose
elements <- aggregated_time
elements <- elements [ order(elements$participant, elements$technique), ]
statstable_time <- ddply(elements,
                         c("participant","technique"),
                         summarise,
                         time=mean(mean_time)
)
elements <- statstable_time

#
elements <- reshape(elements, timevar="technique", idvar=c("participant"), direction="wide")
colnames(elements) <- gsub("time.", "", colnames(elements))


# drop columns with N/A
elements <- na.omit(elements)

data <- elements

techniqueA <- bootstrapMeanCI(data$FROLIC)
techniqueB <- bootstrapMeanCI(data$IBFV)
techniqueC <- bootstrapMeanCI(data$PARTICLES)
techniqueD <- bootstrapMeanCI(data$MIX)

analysisData <- c()

analysisData$name <- c("MIX","PS","IBFV","OLIC")
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("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

path  = paste0("plots/")
filename = paste0("time_means_task_all")

write.table(datatoprint, paste0(path,"printed_",filename,".txt",seq=""), sep=",",row.names=FALSE)

barChart(datatoprint, analysisData$name, nbTechs = 4, ymin = 0, ymax = 10, mycolor = "steelblue3", "", "")

######
#A OLIC
#B IBFV
#C PS
#D MIX

diffBA = bootstrapMeanCI_corr(data$IBFV - data$FROLIC, 4)
diffBC = bootstrapMeanCI_corr(data$IBFV - data$PARTICLES, 4)
diffAC = bootstrapMeanCI_corr(data$FROLIC - data$PARTICLES, 4)
diffDA = bootstrapMeanCI_corr(data$MIX - data$FROLIC, 4)
diffDB = bootstrapMeanCI_corr(data$MIX - data$IBFV, 4)
diffDC = bootstrapMeanCI_corr(data$MIX - data$PARTICLES, 4)


analysisData <- c()
analysisData$name <- c("IBFV-OLIC","IBFV-PS","OLIC-PS", "MIX-OLIC", "MIX-IBFV", "MIX-PS") # Symbol name has been changed in paper to Glyph
analysisData$pointEstimate <- c(diffBA[1], diffBC[1], diffAC[1], diffDA[1], diffDB[1], diffDC[1])
analysisData$ci.max <- c(diffBA[6], diffBC[6], diffAC[6], diffDA[6], diffDB[6], diffDC[6])
analysisData$ci.min <- c(diffBA[5], diffBC[5], diffAC[5], diffDA[5], diffDB[5], diffDC[5])
analysisData$level <- c(diffBA[4], diffBC[4], diffAC[4], diffDA[4], diffDB[4], diffDC[4])
analysisData$ci_corr.max <- c(diffBA[6], diffBC[6], diffAC[6], diffDA[6], diffDB[6], diffDC[6])
analysisData$ci_corr.min <- c(diffBA[5], diffBC[5], diffAC[5], diffDA[5], diffDB[5], diffDC[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("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
filenamediff = paste(featureName, "timediff", sep="_")

barChart(datatoprint, analysisData$name, nbTechs = 6, ymin = -4, ymax = 4, mycolor = "steelblue3", "", "")



plotFeature<- function(featureName) {
  elements <- aggregated_time
  elements <- subset(elements, feature==featureName)
  elements <- elements [ order(elements$participant, elements$technique), ]
  statstable_time <- ddply(elements,
                           c("participant","technique"),
                           summarise,
                           time=mean(mean_time)
  )
  elements <- statstable_time
  
  print(summary(statstable_time))
  #
  elements <- reshape(elements, timevar="technique", idvar=c("participant"), direction="wide")
  colnames(elements) <- gsub("time.", "", colnames(elements))
  
  
  # drop columns with N/A
  elements <- na.omit(elements)
  
  data <- elements
  
  techniqueA <- bootstrapMeanCI(data$FROLIC)
  techniqueB <- bootstrapMeanCI(data$IBFV)
  techniqueC <- bootstrapMeanCI(data$PARTICLES)
  techniqueD <- bootstrapMeanCI(data$MIX)
  
  analysisData <- c()
  
  analysisData$name <- c("MIX","PARTICLES","IBFV","OLIC")
  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("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
  
  path  = paste0("plots/")
  filename = paste(featureName, "timep2", sep="_")
  

  
  barChart(datatoprint, analysisData$name, nbTechs = 4, ymin = 0, ymax = 15, mycolor = "steelblue3", "", "")
  ggsave(paste(filename, "pdf", sep="."), device = cairo_pdf,  width=3.5, height=2.4)
  
  diffBA = bootstrapMeanCI_corr(data$IBFV - data$FROLIC, 24)
  diffBC = bootstrapMeanCI_corr(data$IBFV - data$PARTICLES, 24)
  diffAC = bootstrapMeanCI_corr(data$FROLIC - data$PARTICLES, 24)
  diffDA = bootstrapMeanCI_corr(data$MIX - data$FROLIC, 24)
  diffDB = bootstrapMeanCI_corr(data$MIX - data$IBFV, 24)
  diffDC = bootstrapMeanCI_corr(data$MIX - data$PARTICLES, 24)
  
  
  analysisData <- c()
  analysisData$name <- c("IBFV-OLIC","IBFV-PS","OLIC-PS", "MIX-OLIC", "MIX-IBFV", "MIX-PS") 
  analysisData$pointEstimate <- c(diffBA[1], diffBC[1], diffAC[1], diffDA[1], diffDB[1], diffDC[1])
  analysisData$ci.max <- c(diffBA[6], diffBC[6], diffAC[6], diffDA[6], diffDB[6], diffDC[6])
  analysisData$ci.min <- c(diffBA[5], diffBC[5], diffAC[5], diffDA[5], diffDB[5], diffDC[5])
  analysisData$level <- c(diffBA[4], diffBC[4], diffAC[4], diffDA[4], diffDB[4], diffDC[4])
  analysisData$ci_corr.max <- c(diffBA[6], diffBC[6], diffAC[6], diffDA[6], diffDB[6], diffDC[6])
  analysisData$ci_corr.min <- c(diffBA[5], diffBC[5], diffAC[5], diffDA[5], diffDB[5], diffDC[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("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
  filenamediff = paste(featureName, "timediffp2", sep="_")
  
  barChart(datatoprint, analysisData$name, nbTechs = 6, ymin = -11, ymax = 11, mycolor = "steelblue3", "", "")
  ggsave(paste(filenamediff, "pdf", sep="."), device = cairo_pdf,  width=3.5, height=2.4)
}



plotFeature("Attracting Node")
plotFeature("Attracting Focus")
plotFeature("Repelling Focus")
plotFeature("Repelling Node")
plotFeature("Saddle")
plotFeature("Center")