## LATERAL VISUAL OCCLUSION DOES NOT CHANGE WALKING TRAJECTORIES

# NOTE: for target-heading angle data, positive values indicate that the participant is headed to the RIGHT of the target

rm(list=ls(all=TRUE)) # clear workspace

plot_width = 6 # inches
plot_width_path = 4 # inches
plot_height = 6 # inches

require(ggplot2)
require(zoo)

## Define a couple of functions
stderr <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      stderr = stderr(x[[col]]))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- plyr::rename(data_sum, c("mean" = varname))
  return(data_sum)
}

filepath_source <<- getwd()
filepath_data <<- file.path(filepath_source,'Data')
filepath_output <<- file.path(filepath_source,'Output')
dir.create(filepath_output, showWarnings = FALSE)

filename <- list.files(path = filepath_data,pattern = "\\.csv$") # find all CSV files in directory

participant_id <- sub("(.*?)_.*", "\\1", filename)
trial_code <- sub(".*?_(.*)\\(.*", "\\1", filename)
iteration <- sub(".*?_.*\\((.*)\\).*", "\\1", filename)
route_name <- vector(mode = "character", length = length(filename))
files_with_obstacles <- grepl(".*\\)_(.*?)_segmented.csv", filename) # TODO: This can't be the best way to code this (and the next) line
route_name[files_with_obstacles] <- sub(".*\\)_(.*?)_segmented.csv", "\\1", filename)[files_with_obstacles]

obstacle_locations <<- data.frame(obstacle1_x = c(0,0,0),obstacle1_z = c(3,3,5),obstacle2_x = c(-0.75,0.75,-0.75),obstacle2_z = c(5,5,2),obstacle3_x = c(0.5,-0.5,0.5),obstacle3_z = c(2,2,3)) # define obstacle locations as a global variable

filelist <- data.frame(filename, participant_id, trial_code, iteration, route_name) # create data.frame using the vectors created above
filelist[] <- lapply(filelist, as.character) # convert factors to character # TODO: Why is this step necessary? Can't they be made to START as characters?
filelist$dframe <- lapply(file.path(filepath_data,filename), function(x) read.csv(x, header=TRUE)) # load a CSV file to each member of the list
filelist$iteration <- as.numeric(filelist$iteration) # coerce to numeric

paradigm_environments <- data.frame(corridor = c('0','3','6','9'),obstacles = c('2','5','8','11'),open = c('1','4','7','10'),NoHMD=c('b','m')) # these constants define which type of experimental environment each paradigm code refers to
paradigm_viewing_conditions <- data.frame(binocular = c('0','1','2','b',NA,NA),monocular = c('3','4','5','m',NA),head_referenced_loss = c('6','7','8',NA,NA,NA),trunk_referenced_loss = c('9','10','11',NA,NA,NA),baseline1='baseline1',baseline2='baseline2',baseline3='baseline3',baseleft='baseleft',baseright='baseright') # these constants define which type of experimental environment each paradigm code refers to

filelist$collision <- unlist(lapply(filelist$route_name,function(x) grepl("collision",x,ignore.case=TRUE))) # note walks that led to a collision

target_location <- c(0,7) # def
for(i in 1:nrow(filelist)){
  filelist$dframe[[i]]$head_body_yaw_diff <- filelist$dframe[[i]]$head_yaw - filelist$dframe[[i]]$body_yaw # calculate difference between head and body angle for each frame
    filelist$dframe[[i]]$head_yaw_wrt_target <- filelist$dframe[[i]]$head_yaw - ((atan2(target_location[2] - filelist$dframe[[i]]$head_z,target_location[1] - filelist$dframe[[i]]$head_x) * 180) / pi - 90) # this tells us, in degrees, the extent to which the head was directly facing the target throughout the walk
  
  filelist$range[[i]] <- max(filelist$dframe[[i]]$head_z,na.rm=TRUE) - min(filelist$dframe[[i]]$head_z,na.rm=TRUE) # calculate total Z range of this dataset
  
  for(j in 1:nrow(filelist$dframe[[i]])){ # TODO: Try doing this with an apply()-like function
    filelist$dframe[[i]]$segment_mid_distance_from_target[j] <- mean(c(filelist$dframe[[i]]$segment_start_distance_from_target[j],filelist$dframe[[i]]$segment_end_distance_from_target[j]))
  }
  
  # From trial number, determine the type of paradigm this was
  if(is.element(filelist$trial_code[i],paradigm_environments$corridor)){
    filelist$environment_type[i] <- 'Corridor'
  } else if(is.element(filelist$trial_code[i],paradigm_environments$obstacles)){
    filelist$environment_type[i] <- 'Obstacles'
  } else if(is.element(filelist$trial_code[i],paradigm_environments$open)){
    filelist$environment_type[i] <- 'Open'
  } else if(is.element(filelist$trial_code[i],paradigm_environments$NoHMD)){
    filelist$environment_type[i] <- 'NoHMD'
  } else{
    filelist$environment_type[i] <- ''  
  }
  
  if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$binocular)){
    filelist$viewing_conditions[i] <- 'Binocular'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$monocular)){
    filelist$viewing_conditions[i] <- 'Monocular'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$head_referenced_loss)){
    filelist$viewing_conditions[i] <- 'HeadReferenced'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$trunk_referenced_loss)){
    filelist$viewing_conditions[i] <- 'TrunkReferenced'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$baseline1)){
    filelist$viewing_conditions[i] <- 'Baseline1'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$baseline2)){
    filelist$viewing_conditions[i] <- 'Baseline2'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$baseline3)){
    filelist$viewing_conditions[i] <- 'Baseline3'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$baseleft)){
    filelist$viewing_conditions[i] <- 'BaseLeft'
  } else if(is.element(filelist$trial_code[i],paradigm_viewing_conditions$baseright)){
    filelist$viewing_conditions[i] <- 'BaseRight'
  } else{
    filelist$viewing_conditions[i] <- 'unknown_viewing_conditions'
  }
}

## FUNCTION DEFINITIONS

combine_all_participants_together_and_plot_target_heading_angle <- function(filelist,environment_type,viewing_conditions,iteration){
  # Function to take all trials of a given type and combine into one big data.frame, then plot the mean target-heading angle
  indices_of_dframes_from_this_paradigm <- which(is.element(filelist$environment_type,environment_type) & is.element(filelist$viewing_conditions,viewing_conditions) & is.element(filelist$iteration,iteration)) # find the trials we are interested in. By using is.element(), one can use multiple selections; e.g to average all iterations, send c(1,2,3)
  combined_dframe <- filelist$dframe[[1]][0,] # create new dframe with same columns as filelist
  for(i in 1:length(indices_of_dframes_from_this_paradigm)){
    combined_dframe <- rbind(combined_dframe,filelist$dframe[[indices_of_dframes_from_this_paradigm[i]]]) 
  }
  
  print(paste0("Mean target-heading angle: ",round(mean(combined_dframe$target_heading_angle,na.rm=TRUE),2),'°'))
  print(paste0("Stderr target-heading angle: ",round(stderr(combined_dframe$target_heading_angle),2),'°'))
  
  dframe_summary <- data_summary(combined_dframe,varname = "target_heading_angle", groupnames = "segment_mid_distance_from_target") # calculate mean and stderr values
  
  p <- ggplot(data = combined_dframe,aes(x = segment_mid_distance_from_target,y = target_heading_angle))
  p <- p + geom_point()
  p <- p + geom_errorbar(data = dframe_summary,aes(ymin=target_heading_angle - stderr,ymax=target_heading_angle + stderr), width=0.1, colour = 'red') # error bars showing stderr
  p <- p + xlab("Distance from target (m)") + ylab("Target-heading angle (°)")
  p <- p + scale_x_reverse() # reverse X axis
  #p <- p + ggtitle(paste0(environment_type,' ',viewing_conditions,' (n = ',length(indices_of_dframes_from_this_paradigm),')')) # uncomment this line to include n = <...>
  p <- p + ggtitle(paste0(environment_type,' ',viewing_conditions))
  p <- p + theme_classic()
  ggsave(plot = p, file.path(filepath_output,paste(viewing_conditions[1],environment_type,"target-heading-angle.png",sep='_')), type = "cairo-png",width=plot_width,height=plot_height)
}

get_participant_average_dframes <- function(filelist,environment_type,viewing_conditions,iteration,include_collisions = TRUE){
  # Combine trials to get mean values
  if(include_collisions){
    indices_of_dframes_from_this_paradigm <- is.element(filelist$environment_type,environment_type) & is.element(filelist$viewing_conditions,viewing_conditions) & is.element(filelist$iteration,iteration) # find the trials we are interested in. By using is.element(), one can use multiple selections; e.g to average all iterations, send c(1,2,3)
  }else{ # if we are NOT to include walks with an obstacle collision
    indices_of_dframes_from_this_paradigm <- is.element(filelist$environment_type,environment_type) & is.element(filelist$viewing_conditions,viewing_conditions) & is.element(filelist$iteration,iteration) & !filelist$collision
  }
  aggregated_participant_totals <- filelist[0,] # create new data.frame with same columns as filelist
  
  i = 0 # init
  for(participant_id in unique(filelist[indices_of_dframes_from_this_paradigm,]$participant_id)){ # find runs by the same participant and collapse into mean dframes. This way, each participant contributes equally to the overall statistics (eg if one participant has a particularly curved walk, this avoids them unduly affecting the results)
    i = i + 1 # keep count of how many participants we've got in this paradigm
    indices_of_files_to_use <- which(indices_of_dframes_from_this_paradigm & filelist$participant_id == participant_id) # these are the numbers of the CSV files which contain data pertaining to the paradigm (and route) we are currently interested in AND the participant we are interested in
    dframes_from_this_paradigm_and_participant <- filelist$dframe[indices_of_files_to_use]
    dframe_mean_for_this_participant_doing_this_paradigm <- setNames(as.data.frame(`dim<-`(apply(do.call(cbind,lapply(dframes_from_this_paradigm_and_participant, c, recursive=TRUE)), 1, mean,na.rm=TRUE), dim(dframes_from_this_paradigm_and_participant[[1]]))), colnames(dframes_from_this_paradigm_and_participant[[1]])) # this combines the dataframes into a single frame containing the mean average values for this participant, ignoring NAs
    
    aggregated_participant_totals <- rbind(aggregated_participant_totals,filelist[indices_of_files_to_use[1],]) # add the general info about this paradigm. Since all three iterations should be the same, we just take the data from the first one.
    aggregated_participant_totals$dframe[[i]] <- dframe_mean_for_this_participant_doing_this_paradigm # replace the $dframe entry with the averaged dataframe
  }
  aggregated_participant_totals$filename <- NULL
  aggregated_participant_totals$trial_code <- NULL
  aggregated_participant_totals$route_name <- NULL
  if(length(iteration) == 0){
    aggregated_participant_totals$iteration <- NULL
  }else{
    aggregated_participant_totals$iteration <- as.numeric(aggregated_participant_totals$iteration)
  }
  return(aggregated_participant_totals)
}

get_walk_curvature <- function(real_walk){
  # Define a function that we will run optimise() on:
  generate_equiangle_path <- function(constant_angular_offset,real_walk,target_position){
    model_walk <- matrix(NA,nrow(real_walk),ncol(real_walk)) # init
    model_walk[,2] <- real_walk[,2] # the Z-axis column will be all the real_walk Z positions, as this is where we will be modelling
    #model_walk[,2] <- sort(real_walk[,2],na.last=TRUE) # the Z-axis column will be all the real_walk Z positions, as this is where we will be modelling
    
    previous_non_na_position <- 1 # keep track of the non-NA positions, as we draw our model walk from this
    model_walk[1,] <- 0
    for(i in 2:nrow(model_walk)){ # our first position is the same as the real data. We work from there with the model
      if(!is.na(model_walk[i,2])){ # only model at points where there are data from the real walk
        previous_distance_from_target <- target_position - model_walk[previous_non_na_position,]
        angular_displacement_from_target <- atan((previous_distance_from_target[1]/previous_distance_from_target[2])) / (pi/180)
        this_step_angular_displacement <- angular_displacement_from_target + constant_angular_offset
        model_walk[i,1] <- model_walk[previous_non_na_position,1] + (model_walk[i,2] - model_walk[previous_non_na_position,2]) * tan(this_step_angular_displacement * (pi/180))
        previous_non_na_position <- i
      }
    }
    model_walk[1,] <- NA # clear the first point, as this is (0,0) and was placed deliberately above
    return(model_walk)
  }
  find_best_fit_equiangle_walk <- function(constant_angular_offset,real_walk,target_position){
    model_walk = generate_equiangle_path(constant_angular_offset,real_walk,target_position)
    sum_of_residuals = sum(abs(real_walk[,1] - model_walk[,1]),na.rm=TRUE)
    return(sum_of_residuals)
  }
  target_position <- c(0,7) # def
  real_walk <- real_walk[order(real_walk[,2]),] # order all the real_walk data by the Z axis (if we're working with multiple walks combined, this is necessary)
  
  constant_angular_offset = optimise(find_best_fit_equiangle_walk,c(-60,60),real_walk,target_position)$minimum # find the constant angular offset that best models the walk data. ±60? are just sensible bounds for the final value
  print(paste('Constant angular offset:',round(constant_angular_offset,2)))
  model_walk = generate_equiangle_path(constant_angular_offset,real_walk,target_position)
  return(model_walk)
}

generate_walk_path_plots <- function(filelist,environment_type,viewing_conditions,iteration,route_name = "",plot_mean_walk_path_even_in_obstacle_conditions = FALSE){
  # Plot walk paths (average and/or individual walks)
  
  plot_individual_walks <- TRUE
  
  if(route_name == 'all'){
    indices_of_dframes_from_this_paradigm <- which(is.element(filelist$environment_type,environment_type) & is.element(filelist$viewing_conditions,viewing_conditions) & is.element(filelist$iteration,iteration)) # find the trials we are interested in. By using is.element(), one can use multiple selections; e.g to average all iterations, send c(1,2,3)
  }else{
    indices_of_dframes_from_this_paradigm <- which(is.element(filelist$environment_type,environment_type) & is.element(filelist$viewing_conditions,viewing_conditions) & is.element(filelist$iteration,iteration)  & is.element(filelist$route_name,route_name)) # find the trials we are interested in. By using is.element(), one can use multiple selections; e.g to average all iterations, send c(1,2,3)
  }
  print(paste('n =',length(indices_of_dframes_from_this_paradigm),'for this paradigm')) # report n
  if(length(indices_of_dframes_from_this_paradigm) > 0){ # only continue if there is something to plot
    
    all_walks_x <- vector() # init
    all_walks_z <- vector() # init
    for(walkindex in indices_of_dframes_from_this_paradigm){
      all_walks_x <- append(all_walks_x,filelist$dframe[[walkindex]]$head_x)
      all_walks_z <- append(all_walks_z,filelist$dframe[[walkindex]]$head_z)
    }
    
    all_walks <- matrix(NA,nrow=length(all_walks_x),ncol=2) # init
    all_walks[,1] <- all_walks_x
    all_walks[,2] <- all_walks_z
    
    # TODO: There must surely be a more efficient way to code the above?
    
    model_walk <- get_walk_curvature(all_walks) # get the model walk that best approximates all the walks together
    model_walk <- data.frame(head_x = model_walk[,1],head_z = model_walk[,2]) # convert to data.frame so ggplot2 can handle
    
    p <- ggplot()
    p <- p + xlab("X position (m)") + ylab("Z position (m)")
    p <- p + xlim(-2,2)
    p <- p + ylim(0,7)
    if(tolower(environment_type) == "corridor"){
      p <- p + geom_segment(aes(x = -2, y = 0, xend = -2, yend = 7), size = 3) # add corridor wall line
      p <- p + geom_segment(aes(x = 2, y = 0, xend = 2, yend = 7), size = 3) # add corridor wall line
    }else if(tolower(environment_type) == "obstacles"){
      p <- p + geom_point(aes(x = obstacle_locations$obstacle1_x[iteration], y = obstacle_locations$obstacle1_z[iteration]), shape = 1, size = 5) # add obstacle marker
      p <- p + geom_point(aes(x = obstacle_locations$obstacle2_x[iteration], y = obstacle_locations$obstacle2_z[iteration]), shape = 1, size = 5) # add obstacle marker
      p <- p + geom_point(aes(x = obstacle_locations$obstacle3_x[iteration], y = obstacle_locations$obstacle3_z[iteration]), shape = 1, size = 5) # add obstacle marker
    }
    p <- p + geom_point(aes(x = 0, y = 7), shape = 16, size = 10) # add walk end marker
    #p <- p + ggtitle(paste0(environment_type,' ',viewing_conditions,' (n = ',length(indices_of_dframes_from_this_paradigm),')')) # uncomment this to include n = <...>
    if(environment_type == 'NoHMD'){
      p <- p + ggtitle(viewing_conditions) # for No HMD condition, we don't need to mention the environment type
    }else{
      p <- p + ggtitle(paste0(environment_type,' ',viewing_conditions)) 
    }
    p <- p + coord_fixed() # fix aspect ratio
    p <- p + theme_bw()
    if(plot_individual_walks){ # can overlay the individual walks if desired:
      furthest_position_right <- list() # init
      furthest_position_left <- list() # init
      for(i in indices_of_dframes_from_this_paradigm){
        p <- p + geom_path(data = filelist$dframe[[i]],aes(x = as.numeric(head_x),y = as.numeric(head_z)),alpha=0.2,size = 2,lineend='round',colour="black")#colour=i%%21+1) # TODO: Why is it necessary to use as.numeric here?
        furthest_position_right <- append(furthest_position_right,max(filelist$dframe[[i]]$head_x))
        furthest_position_left <- append(furthest_position_left,min(filelist$dframe[[i]]$head_x))
      }
    }
    
    if(tolower(environment_type) != "obstacles" | plot_mean_walk_path_even_in_obstacle_conditions){ # if NOT an obstacle paradigm
      if(plot_mean_walk_path_even_in_obstacle_conditions){
        mean_plot_colour <- "red"
      }else{
        mean_plot_colour <- "white"
      }
      
      #p <- p + geom_vline(xintercept = max(unlist(furthest_position_right),na.rm=TRUE),size = 0.1) # add line showing rightmost position
      #p <- p + geom_vline(xintercept = min(unlist(furthest_position_left),na.rm=TRUE),size = 0.1) # add line showing leftmost position
      #p <- p + geom_point(data = model_walk,aes(x = head_x,y=head_z),colour = 'red') # add the calculated model walk
      
      # Add mean and stderr walk path on top:
      dframe_mean <- setNames(as.data.frame(`dim<-`(apply(do.call(cbind,lapply(filelist$dframe[indices_of_dframes_from_this_paradigm], c, recursive=TRUE)), 1, mean,na.rm=FALSE), dim(filelist$dframe[[1]]))), colnames(filelist$dframe[[1]])) # this combines the dataframes into a single frame containing the mean average values for each, ignoring NAs
      for(column_name in names(dframe_mean)){
        names(dframe_mean)[names(dframe_mean) == column_name] <- paste0(column_name,"_mean") # change name of each column, so we can differentiate from the stderr data 
      }
      
      dframe_stderr <- setNames(as.data.frame(`dim<-`(apply(do.call(cbind,lapply(filelist$dframe[indices_of_dframes_from_this_paradigm], c, recursive=TRUE)), 1, stderr), dim(filelist$dframe[[1]]))), colnames(filelist$dframe[[1]])) # this combines the dataframes into a single frame containing the stderr average values for each, ignoring NAs
      for(column_name in names(dframe_stderr)){
        names(dframe_stderr)[names(dframe_stderr) == column_name] <- paste0(column_name,"_stderr") # change name of each column, so we can differentiate from the stderr data 
      }
      
      dframe_combined <- cbind(dframe_mean,dframe_stderr) # combine dframes
      
      p <- p + geom_point(data = dframe_combined,aes(x = head_x_mean,y=head_z_mean),colour=mean_plot_colour) # add the mean walk position 
      p <- p + geom_errorbarh(data = dframe_combined,aes(x = head_x_mean,xmin = head_x_mean - head_x_stderr,xmax = head_x_mean + head_x_stderr,height=0.1, y = head_z_mean),colour = mean_plot_colour) # error bars showing stderr
      
      ggsave(plot = p, file.path(filepath_output,paste(environment_type,viewing_conditions,as.character(iteration[1]),route_name,"walk_path.png",sep='_')),type = "cairo-png",width=plot_width_path,height=plot_height)
    }else{
      dframe_combined <- NULL # this object not used here
      ggsave(plot = p, file.path(filepath_output,paste(environment_type,iteration,viewing_conditions,route_name,"walk_path.png",sep='_')),type = "cairo-png",width=plot_width_path,height=plot_height)
    }
    
    if(tolower(environment_type) != "obstacles"){
      # Produce another version with zoomed in view:
      p <- p + xlim(-0.5,0.5) + ylim(3,4)
      p <- p + coord_cartesian()
      p <- p + geom_vline(xintercept=0,linetype="dotted") + geom_hline(yintercept = 3.5,linetype="dotted")
      p <- p + ggtitle(NULL)
      p <- p + theme_void()
      ggsave(plot = p, file.path(filepath_output,paste(environment_type,viewing_conditions,"walk_path-zoomed.png",sep='_')),type = "cairo-png",width=plot_width_path,height=1)
    }
  }else{
    dframe_combined <- NULL # this object not used here
  }
  return(dframe_combined)
}

library(reshape2)
library(effsize) # used to get Hedge's g

calibrate = TRUE
interpolate_data = FALSE

for(i in 1:nrow(filelist)){
  # POST-HOC CALIBRATION:
  if(calibrate){
    if(is.element(filelist$trial_code[i],paradigm_environments$NoHMD)){
      # Apply calibration factor to all non-HMD trials:
      filelist$dframe[[i]]$head_x <-  filelist$dframe[[i]]$head_x + 0.17 # apply 17cm linear correction factor to X axis due to known calibration error in the room/cameras
      filelist$dframe[[i]]$body_x <-  filelist$dframe[[i]]$body_x + 0.17 # as above
    } else{
      # Apply calibration factor to all HMD trials:
      filelist$dframe[[i]]$head_x <-  filelist$dframe[[i]]$head_x + 0.07 # apply 7cm linear correction factor to X axis due to known calibration error in the room/cameras
      filelist$dframe[[i]]$body_x <-  filelist$dframe[[i]]$body_x + 0.07 # as above
    }
  }
  # CUBIC SPLINE INTERPOLATION:
  if(interpolate_data & sum(!is.na(filelist$dframe[[i]]$target_heading_angle)) >= 2){ # as long as some data are present that can be interpolated (i.e. at least two points)
    print(paste0(mean(unlist(lapply(filelist$dframe,function(x) sum(is.na(x$head_z))/length(x$head_z)))),'% of the data will be interpolated'))
    print("Interpolating with cubic splines...")
    filelist$dframe[[i]][c("target_heading_angle","head_yaw","head_body_yaw_diff","head_z","head_x")] <- lapply(filelist$dframe[[i]][c("target_heading_angle","head_yaw","head_body_yaw_diff","head_z","head_x")], function(x) na.spline(x) + 0*na.approx(x,na.rm=FALSE)) # interpolate NAs using cubic splines
  }
}

###############################################
# EXPERIMENT 1: No HMD binocular vs. monocular
###############################################

# Determine total change in mean target-heading angle across all trials, separated into binocular and monocular conditions
aggregated_participant_totals.binocular <- get_participant_average_dframes(filelist,'NoHMD','Binocular',1:3)
aggregated_participant_totals.binocular$mean_target_heading_angle <- sapply(aggregated_participant_totals.binocular$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
aggregated_participant_totals.monocular <- get_participant_average_dframes(filelist,'NoHMD','Monocular',1:3)
aggregated_participant_totals.monocular$mean_target_heading_angle <- sapply(aggregated_participant_totals.monocular$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
print(paste("Target-heading angle",round(mean(aggregated_participant_totals.monocular$mean_target_heading_angle,na.rm=TRUE) - mean(aggregated_participant_totals.binocular$mean_target_heading_angle,na.rm=TRUE),2),"to the right under monocular viewing as compared to binocular"))

# Determine whether there is an order effect for target-heading angle in the monocular and binocular no HMD trials:
aggregated_participant_totals.experiment1 <- NULL # init
for(viewing_conditions in c('Binocular','Monocular')){
  for(iteration in c(1:3)){
    print(paste(viewing_conditions,iteration))
    aggregated_participant_totals <- get_participant_average_dframes(filelist,'NoHMD',viewing_conditions,iteration) # NOTE: In this instance, there is only one file for each participant. However, we call this function as an easy way of pulling out the necessary data
    aggregated_participant_totals$iteration <- iteration
    aggregated_participant_totals$mean_target_heading_angle <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
    
    aggregated_participant_totals$mean_total_speed <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$total_speed,na.rm=TRUE))
    
    print(paste("Mean target-heading angle:",round(mean(aggregated_participant_totals$mean_target_heading_angle,na.rm=TRUE),2)))
    print(paste("Stderr target-heading angle:",round(stderr(aggregated_participant_totals$mean_target_heading_angle),2)))
    
    aggregated_participant_totals.experiment1 <- rbind(aggregated_participant_totals.experiment1,aggregated_participant_totals) # add all conditions into a single dataframe
  }
  
  # Now do the above for all binocular or monocular trials combined (as we report the overall SD in the paper):
  aggregated_participant_totals <- get_participant_average_dframes(filelist,'NoHMD',viewing_conditions,1:3)
  aggregated_participant_totals$mean_target_heading_angle <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
  
  print(paste("Overall",viewing_conditions,"mean target-heading angle:",round(mean(aggregated_participant_totals$mean_target_heading_angle,na.rm=TRUE),2)))
  print(paste("Overall",viewing_conditions,"stderr target-heading angle:",round(stderr(aggregated_participant_totals$mean_target_heading_angle),2)))
}

aggregated_participant_totals.experiment1$viewing_conditions <- as.factor(aggregated_participant_totals.experiment1$viewing_conditions) # convert to factor
aggregated_participant_totals.experiment1$iteration <- as.factor(aggregated_participant_totals.experiment1$iteration) # convert to factor

aggregated_participant_totals.experiment1.reshaped <- dcast(aggregated_participant_totals.experiment1, participant_id ~ iteration + viewing_conditions,value.var = 'mean_target_heading_angle')
write.table(aggregated_participant_totals.experiment1.reshaped,'forJASP_experiment1.csv',row.names = FALSE,na="")

aggregated_participant_totals.experiment1.reshaped_total_speed <- dcast(aggregated_participant_totals.experiment1, participant_id ~ iteration + viewing_conditions,value.var = 'mean_total_speed')
write.table(aggregated_participant_totals.experiment1.reshaped_total_speed,'forJASP_experiment1_totalspeed.csv',row.names = FALSE,na="")

# Plot mean walk paths for both paradigms:
for(viewing_conditions in c('Binocular','Monocular')){
  generate_walk_path_plots(filelist,'NoHMD',viewing_conditions,1:3)
  combine_all_participants_together_and_plot_target_heading_angle(filelist,'NoHMD',viewing_conditions,1:3) # get target-heading plots
}

## BOOTSTRAPPING

# CHECK FOR MAIN EFFECT OF BINOCULAR VS MONOCULAR VIEWING:
print('Monocular vs binocular')
g <- cohen.d(aggregated_participant_totals.monocular[,'mean_target_heading_angle'],aggregated_participant_totals.binocular[,'mean_target_heading_angle'],hedges.correction=TRUE,na.rm=TRUE,paired=TRUE) # calculate Hedge's g (rather than Cohen's d as sample size < 20)
print(paste('Hedge\'s g:',round(g$estimate,2),'// 95% CIs:',round(g$conf.int[[1]],2),'to',round(g$conf.int[[2]],2))) # this is when each condition is compared to Binocular viewing

# CHECK FOR ORDER EFFECT:
for(iteration in c('1','2','3')){
  print(paste('',iteration,'',sep = ' *** ')) # print name of the group we are tesing this time
  g <- cohen.d(aggregated_participant_totals.experiment1.reshaped[,paste0(iteration,'_Monocular')],aggregated_participant_totals.experiment1.reshaped[,paste0(iteration,'_Binocular')],hedges.correction=TRUE,na.rm=TRUE,paired=TRUE) # calculate Hedge's g (rather than Cohen's d as sample size < 20)
  print(paste('Hedge\'s g:',round(g$estimate,2),'// 95% CIs:',round(g$conf.int[[1]],2),'to',round(g$conf.int[[2]],2)))
}

##################################
# EXPERIMENT 2(A): HMD walk paths
##################################
# Determine total change in mean target-heading angle across all open evironment trials, separated into viewing conditions
for(viewing_conditions in c('Binocular','Monocular','HeadReferenced','TrunkReferenced')){
  aggregated_participant_totals <- get_participant_average_dframes(filelist,c('Corridor','Open'),viewing_conditions,1:3)
  aggregated_participant_totals$mean_target_heading_angle <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
  print(paste("Mean target-heading angle for all trials under",viewing_conditions,"conditions =",round(mean(aggregated_participant_totals$mean_target_heading_angle,na.rm=TRUE),3)))
}

# Determine total change in mean target-heading angle across all trials, separated into binocular and monocular conditions
aggregated_participant_totals.binocular <- get_participant_average_dframes(filelist,c('Corridor','Open'),'Binocular',1:3)
aggregated_participant_totals.binocular$mean_target_heading_angle <- sapply(aggregated_participant_totals.binocular$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
aggregated_participant_totals.monocular <- get_participant_average_dframes(filelist,c('Corridor','Open'),'Monocular',1:3)
aggregated_participant_totals.monocular$mean_target_heading_angle <- sapply(aggregated_participant_totals.monocular$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
print(paste("Target-heading angle",round(mean(aggregated_participant_totals.monocular$mean_target_heading_angle,na.rm=TRUE) - mean(aggregated_participant_totals.binocular$mean_target_heading_angle,na.rm=TRUE),2),"to the right under monocular viewing as compared to binocular"))

g <- cohen.d(aggregated_participant_totals.monocular[,'mean_target_heading_angle'],aggregated_participant_totals.binocular[,'mean_target_heading_angle'],hedges.correction=TRUE,na.rm=TRUE,paired=TRUE) # calculate Hedge's g (rather than Cohen's d as sample size < 20)
print(paste('Hedge\'s g:',round(g$estimate,2),'// 95% CIs:',round(g$conf.int[[1]],2),'to',round(g$conf.int[[2]],2))) # this is when each condition is compared to Binocular viewing

# Get participant-averaged dframes for each paradigm type:
aggregated_participant_totals.experiment2 <- NULL # init
for(environment_type in c('Corridor','Open')){
  for(viewing_conditions in c('Binocular','Monocular','HeadReferenced','TrunkReferenced')){
    combine_all_participants_together_and_plot_target_heading_angle(filelist,environment_type,viewing_conditions,1:3)
    generate_walk_path_plots(filelist,environment_type,viewing_conditions,1:3)
    aggregated_participant_totals <- get_participant_average_dframes(filelist,environment_type,viewing_conditions,c(1:3))
    aggregated_participant_totals$mean_target_heading_angle <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$target_heading_angle,na.rm=TRUE)) # grab mean target-heading angle from each dframe
    aggregated_participant_totals$constant_angular_offset <- sapply(aggregated_participant_totals$dframe,function(x) (x$target_heading_angle[1])) # grab constant angular offset value. Since these are always the same across all rows, just take the first value from each
    
    aggregated_participant_totals$mean_head_yaw <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$head_yaw,na.rm=TRUE))
    aggregated_participant_totals$mean_body_yaw <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$body_yaw,na.rm=TRUE))
    aggregated_participant_totals$mean_head_body_yaw_diff <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$head_body_yaw_diff,na.rm=TRUE))
    
    aggregated_participant_totals$mean_total_speed <- sapply(aggregated_participant_totals$dframe,function(x) mean(x$total_speed,na.rm=TRUE))
    
    print(paste0("Mean target-heading angle: ",round(mean(aggregated_participant_totals$mean_target_heading_angle,na.rm=TRUE),2),'°'))
    print(paste0("Stderr target-heading angle: ",round(stderr(aggregated_participant_totals$mean_target_heading_angle),2),'°'))
    print(paste0("Mean constant angular offset: ",round(mean(aggregated_participant_totals$constant_angular_offset,na.rm=TRUE),2),'°'))
    print(paste0("Stderr constant angular offset: ",round(stderr(aggregated_participant_totals$constant_angular_offset),2),'°'))
    print(paste0("Mean head yaw: ",round(mean(aggregated_participant_totals$mean_head_yaw,na.rm=TRUE),2),'°'))
    print(paste0("Stderr head yaw: ",round(stderr(aggregated_participant_totals$mean_head_yaw),2),'°'))
    print(paste0("Mean body yaw: ",round(mean(aggregated_participant_totals$mean_body_yaw,na.rm=TRUE),2),'°'))
    print(paste0("Stderr body yaw: ",round(stderr(aggregated_participant_totals$mean_body_yaw),2),'°'))
    print(paste0("Mean head/body yaw difference: ",round(mean(aggregated_participant_totals$mean_head_body_yaw_diff,na.rm=TRUE),2),'°'))
    print(paste0("Stderr head/body yaw difference: ",round(stderr(aggregated_participant_totals$mean_head_body_yaw_diff),2),'°'))
    
    aggregated_participant_totals.experiment2 <- rbind(aggregated_participant_totals.experiment2,aggregated_participant_totals) # add all conditions into a single dataframe
  }
}

# Draw individual walk paths for obstacle conditions:
for(viewing_conditions in c('Binocular','Monocular','HeadReferenced','TrunkReferenced')){
  for(iteration in 1:3){
    generate_walk_path_plots(filelist,'Obstacles',viewing_conditions,iteration,'all')
  }
}

aggregated_participant_totals.experiment2.reshaped.mean_target_heading_angle <- dcast(aggregated_participant_totals.experiment2,participant_id ~ viewing_conditions + environment_type,value.var = 'mean_target_heading_angle')
write.table(aggregated_participant_totals.experiment2.reshaped.mean_target_heading_angle,'forJASP_experiment2_mean_target_heading_angle.csv',row.names = FALSE,na="")
aggregated_participant_totals.experiment2.reshaped.mean_head_yaw <- dcast(aggregated_participant_totals.experiment2,participant_id ~ viewing_conditions + environment_type,value.var = 'mean_head_yaw')
write.table(aggregated_participant_totals.experiment2.reshaped.mean_head_yaw,'forJASP_experiment2_mean_head_yaw_wrt_target.csv',row.names = FALSE,na="")

# Check that walking speed remained the same across trials:
aggregated_participant_totals.experiment2.reshaped.walking_speed <- dcast(aggregated_participant_totals.experiment2,participant_id ~ viewing_conditions + environment_type,value.var = 'mean_total_speed')
write.table(aggregated_participant_totals.experiment2.reshaped.walking_speed,'forJASP_experiment2_mean_total_speed.csv',row.names = FALSE,na="")

for(dependent_variable in c('mean_target_heading_angle','mean_head_yaw')){ # for both THA and head yaw...
  data <- switch(dependent_variable,
                 mean_target_heading_angle = aggregated_participant_totals.experiment2.reshaped.mean_target_heading_angle,
                 mean_head_yaw = aggregated_participant_totals.experiment2.reshaped.mean_head_yaw)
  print(paste('',dependent_variable,'',sep = ' ***** '))
  for(viewing_conditions in c('Monocular','HeadReferenced','TrunkReferenced')){
    for(environment in c('Corridor','Open')){
      print('')
      print(paste('*',viewing_conditions,environment,'*')) # print name of the group we are tesing this time
      mean_shift_from_binocular <- mean(data[,paste0(viewing_conditions,'_',environment)],na.rm=TRUE) - mean(data[,paste0('Binocular_',environment)],na.rm=TRUE)
      print(paste('Mean shift (°): ',round(mean_shift_from_binocular,2))) # positive means head is further RIGHT of the target than under binocular conditions
      
      g <- cohen.d(data[,paste0(viewing_conditions,'_',environment)],data[,paste0('Binocular_',environment)],hedges.correction=TRUE,na.rm=TRUE,paired=TRUE) # calculate Hedge's g (rather than Cohen's d as sample size < 20)
      print(paste('Hedge\'s g:',round(g$estimate,2),'// 95% CIs:',round(g$conf.int[[1]],2),'to',round(g$conf.int[[2]],2)))
    }
  }
}

aggregated_participant_totals.experiment2$viewing_conditions <- factor(aggregated_participant_totals.experiment2$viewing_conditions,levels = c("Binocular","Monocular","HeadReferenced","TrunkReferenced")) # set as factor so we order by levels

p <- ggplot(aggregated_participant_totals.experiment2,aes(x = viewing_conditions,y = mean_head_body_yaw_diff))
p <- p + geom_boxplot()
p <- p + theme_bw()
p <- p + scale_x_discrete("Viewing conditions",labels = c("Binocular"="Binocular","Monocular"="Monocular","HeadReferenced"="Head-Referenced Loss","TrunkReferenced"="Trunk-Referenced Loss"))
p <- p + xlab("Viewing conditions") + ylab("Head yaw (°)")
ggsave(file.path(filepath_output,"head-yaw-wrt-target.png"),type = "cairo-png",width=7,height=4)

print("*** Mean head yaw wrt target ***")
for(conditions in c('Binocular','HeadReferenced','TrunkReferenced','Monocular')){
  print(conditions)
  print(mean(subset(aggregated_participant_totals.experiment2,viewing_conditions==conditions)$mean_head_yaw))
}

#############################
# EXPERIMENT 2(B): Obstacles
#############################

all_obstacle_data <- NULL # init
for(viewing_conditions in c('Binocular','HeadReferenced','TrunkReferenced','Monocular')){
  for(iteration in c(1:3)){
    generate_walk_path_plots(filelist,'Obstacles',viewing_conditions,iteration)
    aggregated_participant_totals <- get_participant_average_dframes(filelist,'Obstacles',viewing_conditions,iteration,FALSE) # FALSE indicates we will not be including walks with collisions
    
    aggregated_participant_totals$closest_obstacle_on_left_side <- lapply(aggregated_participant_totals$dframe,function(x) x$closest_obstacle_on_left_side[1])
    aggregated_participant_totals$closest_obstacle_on_right_side <- lapply(aggregated_participant_totals$dframe,function(x) x$closest_obstacle_on_right_side[1])
    for(i in 1:nrow(aggregated_participant_totals)){ # TODO: Surely this can be done with an apply()-like function?
      aggregated_participant_totals$obstacle_passing_distance_ratio[i] <- aggregated_participant_totals$closest_obstacle_on_left_side[[i]] / aggregated_participant_totals$closest_obstacle_on_right_side[[i]]
      aggregated_participant_totals$mean_obstacle_distance[i] <- mean(c(aggregated_participant_totals$closest_obstacle_on_left_side[[i]],aggregated_participant_totals$closest_obstacle_on_right_side[[i]]),na.rm=TRUE)
    }
    
    print(paste(viewing_conditions,iteration,'L:R obstacle passing ratio:',round(mean(aggregated_participant_totals$obstacle_passing_distance_ratio,na.rm=TRUE),2),': 1'))
    print(paste('mean obstacle distance:',round(mean(aggregated_participant_totals$mean_obstacle_distance,na.rm=TRUE),2),'m'))
    
    all_obstacle_data <- rbind(all_obstacle_data,aggregated_participant_totals[c('participant_id','viewing_conditions','iteration','closest_obstacle_on_left_side','closest_obstacle_on_right_side','obstacle_passing_distance_ratio','mean_obstacle_distance')])
  }
}

# Plot mean walk paths for obstacle sets 1 and 2:
for(iteration in c(1:3)){
  for(viewing_conditions in c('Binocular','HeadReferenced','TrunkReferenced','Monocular')){
    print(paste("Obstacle set",iteration,":",viewing_conditions))
    for(route_name in c("route1","route2","route3","route4")){
      print(route_name)
      generate_walk_path_plots(filelist,'Obstacles',viewing_conditions,iteration,route_name,TRUE)
    }
  }
}

# COMBINED PLOTS TO DEMONSTRATE (NON-)EFFECT OF VIEWING CONDITIONS ON WALKING IN OBSTACLE TASKS
for(iteration in 1:3){
  for(route_name in c('route1','route2','route3','route4')){
    aggregate_plot_data <- data.frame(viewing_conditions = character(),iteration = integer(),route_name = character(),dframe_mean = data.frame(),dframe_stderr = data.frame()) # init
    for(viewing_conditions in c('Binocular','HeadReferenced','TrunkReferenced','Monocular')){
      dframe_combined <- generate_walk_path_plots(filelist,'Obstacles',viewing_conditions,iteration,route_name,TRUE)
      
      # Prepare row to be added to data.frame:
      plot_data <- data.frame(viewing_conditions = viewing_conditions,iteration = iteration,route_name = route_name) # init
      if(!is.null(dframe_combined)){
        plot_data$dframe_combined[[1]] <- dframe_combined # we have to add this column manually, or R tries to expand the data.frame to add it to the parent frame
      }else{
        plot_data$dframe_combined[[1]] <- data.frame() # just use a blank dataframe, as this row has no data
      }
      
      aggregate_plot_data <- rbind(aggregate_plot_data,plot_data)
    }
    
    # Produce combined plot:
    p <- ggplot()
    p <- p + xlab("X position (m)") + ylab("Z position (m)")
    p <- p + xlim(-2,2)
    p <- p + ylim(0,7)
    p <- p + geom_point(aes(x = obstacle_locations$obstacle1_x[iteration], y = obstacle_locations$obstacle1_z[iteration]), shape = 1, size = 5) # add obstacle marker
    p <- p + geom_point(aes(x = obstacle_locations$obstacle2_x[iteration], y = obstacle_locations$obstacle2_z[iteration]), shape = 1, size = 5) # add obstacle marker
    p <- p + geom_point(aes(x = obstacle_locations$obstacle3_x[iteration], y = obstacle_locations$obstacle3_z[iteration]), shape = 1, size = 5) # add obstacle marker
    p <- p + coord_fixed() # fix aspect ratio
    p <- p + theme_bw()
    for(i in 1:nrow(aggregate_plot_data)){
      if(aggregate_plot_data$viewing_conditions[i] == "Binocular"){
        plot_colour = "red"
      }else if(aggregate_plot_data$viewing_conditions[i] == "HeadReferenced"){
        plot_colour = "yellow"
      }else if(aggregate_plot_data$viewing_conditions[i] == "TrunkReferenced"){
        plot_colour = "blue"
      }else if(aggregate_plot_data$viewing_conditions[i] == "Monocular"){
        plot_colour = "green"
      }
      if(length(aggregate_plot_data$dframe_combined[[i]]) > 0){ # if there are any data to draw...
        p <- p + geom_point(data = aggregate_plot_data$dframe_combined[[i]],aes(x = head_x_mean,y=head_z_mean),colour=plot_colour,alpha=0.5)
        p <- p + geom_errorbarh(data = aggregate_plot_data$dframe_combined[[i]],aes(x = head_x_mean,xmin = head_x_mean - head_x_stderr,xmax = head_x_mean + head_x_stderr,height=0.1, y = head_z_mean),colour = plot_colour,alpha=0.5) # error bars showing stderr
      }
    }
    p <- p + scale_colour_manual(name = 'Viewing conditions',values =c('Binocular'='red','HeadReferenced'='yellow','TrunkReferenced'='blue','Monocular'='green'))
    ggsave(file.path(filepath_output,paste("combined_obstacle_plots_",iteration,route_name,"_path.png",sep='_')),type = "cairo-png",width=plot_width_path,height=plot_height)
  }
}

all_obstacle_data.reshaped_mean_obstacle_distance <- dcast(all_obstacle_data,participant_id ~ iteration + viewing_conditions,value.var = 'mean_obstacle_distance')
write.table(all_obstacle_data.reshaped_mean_obstacle_distance,'forJASP_obstacles_mean_obstacle_distance.csv',row.names = FALSE,na="")
all_obstacle_data.reshaped_obstacle_passing_distance_ratio <- dcast(all_obstacle_data,participant_id ~ iteration + viewing_conditions,value.var = 'obstacle_passing_distance_ratio')
write.table(all_obstacle_data.reshaped_obstacle_passing_distance_ratio,'forJASP_obstacles_obstacle_passing_distance_ratio.csv',row.names = FALSE,na="")

all_obstacle_data$iteration <- as.factor(all_obstacle_data$iteration)

## BOOTSTRAPPING
for(dependent_variable in c('obstacle_passing_distance_ratio','mean_obstacle_distance')){ # for both THA and head yaw...
  data <- switch(dependent_variable,
                 obstacle_passing_distance_ratio = dcast(all_obstacle_data,participant_id ~ viewing_conditions,fun.aggregate = mean,value.var = 'mean_obstacle_distance'), # aggregate mean of all iterations, as we're not interested in per-iteration values.
                 mean_obstacle_distance = dcast(all_obstacle_data,participant_id ~ viewing_conditions,fun.aggregate = mean,value.var = 'obstacle_passing_distance_ratio')) # aggregate mean of all iterations, as we're not interested in per-iteration values.
  print(paste('',dependent_variable,'',sep = ' ***** '))
  for(viewing_conditions in c('Monocular','HeadReferenced','TrunkReferenced')){
    print(paste('*',viewing_conditions,'*')) # print name of the group we are tesing this time
    
    mean_shift_from_binocular <- mean(data[,viewing_conditions],na.rm=TRUE) - mean(data[,'Binocular'],na.rm=TRUE)
    print(paste('Mean shift:',round(mean_shift_from_binocular,2))) # positive means head is further RIGHT of the target than under binocular conditions
    
    g <- cohen.d(data[,viewing_conditions],data[,'Binocular'],hedges.correction=TRUE,na.rm=TRUE,paired=TRUE) # calculate Hedge's g (rather than Cohen's d as sample size < 20)
    print(paste('Hedge\'s g:',round(g$estimate,2),'// 95% CIs:',round(g$conf.int[[1]],2),'to',round(g$conf.int[[2]],2)))
  }
}

## Get average head angle with respect to target across participants in obstacle-free trials for HeadReferenced viewing conditions:
aggregated_participant_totals <- get_participant_average_dframes(filelist,c('Open','Corridor'),c('HeadReferenced'),1:3)
for(i in 1:nrow(aggregated_participant_totals)){ # TODO: Surely this can be done with an apply()-like function?
  aggregated_participant_totals$head_yaw_wrt_target.unsignedmean[i] <- mean(abs(aggregated_participant_totals$dframe[[i]]$head_yaw_wrt_target))
}
print(paste('Mean unsigned head yaw deviation from pointing towards the target across participants in non-obstacle HeadReferenced trials is',round(mean(aggregated_participant_totals$head_yaw_wrt_target.unsignedmean),2),'deg ±',round(sd(aggregated_participant_totals$head_yaw_wrt_target.unsignedmean),2),'(SD)'))