Import packages

 library(RColorBrewer)
 library(ggplot2)
 library(scales)
 library(gridExtra)
 library(mefa)

Helper functions

A function to return integer valued breaks for chart scales

integer_breaks <- function(n = 5, ...) {
  breaker <- pretty_breaks(n, ...)
  function(x) {
     breaks <- breaker(x)
     breaks[breaks == floor(breaks)]
  }
}

A function to load the gradient data

loadGradientData <- function(year, range, hosts) {
  ### Load the data
  input_file_name <- paste(year, "VCM_gradientAnalysis_", range, "m_", hosts, "hosts_0.csv", sep="")
  input_file_path <- file.path(paste("../../results/useCase1/1A/",year, sep=""),  input_file_name)
  data <- read.csv(input_file_path, head=TRUE,sep=",",stringsAsFactors=TRUE)

  ### Replace values where no gradient exists
  data <- subset(data, Gradient != "")
  data <- droplevels(data)
  data[data == -1] <- NA
  data$Gradient <- droplevels(data$Gradient)
  data
}

A function to load the runner passing time data.

loadRunnerData <- function(year) {
  ### Load the data
  input_file_path <- file.path(paste("../../results/useCase1/1A/",year, sep=""),  "checkpointPasses.csv")
  data <- read.csv(input_file_path, head=TRUE,sep=",",stringsAsFactors=FALSE)
  # Transpose the data to have checkpoints as columns
  data <- t(data)
  colnames(data) <- data[1,]
  data[-1,]
}

A function to crate the basic plot template (scales, axes, and lables):

generatePlot <- function(gradientData, year, range) {
  ### Construct the plot ommitting NA
  plot <- ggplot(na.omit(gradientData)) 
  
  ### Add axis scales
  plot <- plot +  coord_cartesian(xlim = c(0, 36000), ylim = c(0,120))  
  
  ### Add axis ticks
  plot = plot + scale_x_continuous(breaks=c(3600,7200,10800,14400,18000,21600,25200,28800,32400), labels=c("08:00", "09:00","10:00","11:00","12:00","13:00","14:00","15:00","16:00")) +
 scale_y_continuous(breaks = integer_breaks())
  
  ### Add plot and axis labels. Set legend to inclue all levels
  plot + ggtitle(paste("Gradient Analysis - range: ", range,", year: ", year, sep="")) + theme(panel.background = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1, vjust=1, size=8), axis.text.y = element_text(size=8), panel.margin = unit(0.7, "lines")) +  scale_fill_discrete(drop=FALSE)
}

Adds the gradient stability bar chart to the plot:

withGradientStability <- function(chart) {
  chart + geom_bar(aes(x = Time, y = Hops, fill = Gradient, width=20),  stat = "identity")
}

Shades the temporal area where runners are between checkpoints:

withRunnerPassingRange <- function(chart, runnerTimes) {
    ranges <- data.frame(xmin=numeric(9), xmax=numeric(9), ymin=0, ymax=Inf, Gradient=character(9), stringsAsFactors=FALSE)

  for(checkpoint in seq(1:9)){
    ranges$xmin[checkpoint] <-  as.numeric(as.character(min(runnerTimes[,checkpoint])))
    ranges$xmax[checkpoint] <-  as.numeric(as.character(max(runnerTimes[,checkpoint + 1])))
    ranges$Gradient[checkpoint] <- paste("C", checkpoint - 1, "-C", checkpoint, sep="") 
  } 
  ### Convert Gradient column to factor.
  ranges$Gradient <- as.factor(ranges$Gradient)
  ### Add shading to chart representing when people pass.
  chart + geom_rect(data=ranges, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), alpha = 0.8, fill = 'grey')
}

Adds horizontal lines to show when a transmission window is active for a given transmission speed

withTransmissionWindow <- function(chart, summary, speeds) {
    # Copy frame, repeating for number of speeds, add a level indicator
    frameSize <- nrow(summary)
    newdata <- rep(summary, length(speeds))
    newdata$level <- rep(1:length(speeds), each=frameSize)
    # Add a column representing the transmission speed
    newdata$transmissionSpeed <- rep(speeds, each=frameSize)
    # Add a column representing the transmission time
    newdata$transmissionTime <- with(newdata, transmissionSpeed*round(averageHops,0))
    # Add a column representing the effective start time
    newdata$effectiveStartTime <- with(newdata, startTime + transmissionTime)
    # Add a column representing the effective end time
    newdata$effectiveEndTime <- with(newdata, endTime - transmissionTime)
    
    
    # Take the subset of valid transmission windows and chart them
    window <- subset(newdata, effectiveStartTime < effectiveEndTime)
    if(nrow(window) > 0){
        chart <- chart + geom_segment(data=window, aes(x=effectiveStartTime, xend=effectiveEndTime, y=20+(10*level), yend=20+(10*level)), show_guide=FALSE) 
    }
chart
}

Facet the chart by gradient:

### Create a faceted chart by gradient.
withGradientFacets <- function(chart) {
   chart + facet_grid( Gradient ~ ., drop=FALSE)
}

Calculate when stable gradients exist:

calculateStableGradients <- function(gradientData) {
      
    # Data for each gradient  
    Gradient <- character()
    startTime <- numeric()
    endTime <- numeric()
    totalHops <- numeric()
   
    allLevels <- levels(gradientData$Gradient) 
    for(aLevel in allLevels){
        data <- subset(gradientData, gradientData$Gradient == aLevel)
        tempStartTime <- 0
        tempEndTime <- 0
        tempHops <- 0
        gradientActive <- FALSE
        
        # Iterate through data
        for(row in 1:nrow(data)){
            # Extract the hop value from the row
            hopValue <- data[row,3]
            
            if(gradientActive){
                if(is.na(hopValue)){
                    tempEndTime <- data[row-1,1]
                    # commit totals and reset
                    Gradient <- c(Gradient, aLevel)
                    startTime <- c(startTime, tempStartTime)
                    endTime <- c(endTime, tempEndTime)
                    totalHops <- c(totalHops, tempHops)
                    tempHops <- 0
                    gradientActive <- FALSE
                }else{
                    tempHops <- tempHops + hopValue
                }
            }else{
                if(is.na(hopValue)){
                    next
                }else{
                    tempStartTime <- data[row,1]
                    tempHops <- tempHops + hopValue
                    gradientActive <- TRUE
                }
            }
        }
        
        # Catch case where active gradient as last step
        if(gradientActive){
              tempEndTime <- data[nrow(data),1]
              # commit totals
              Gradient <- c(Gradient, aLevel)
              startTime <- c(startTime, tempStartTime)
              endTime <- c(endTime, tempEndTime)
              totalHops <- c(totalHops, tempHops)
              tempHops <- 0
              gradientActive <- FALSE
        }
    }
    
    # Create data frame from rows
    frame <- data.frame(Gradient, startTime, endTime, totalHops)
    frame$Gradient <- as.factor(frame$Gradient)
    # Add columns for steps and average hops
    frame$stepCount <- with(frame, 1 + (endTime - startTime) / 20)
    frame$elapsedSeconds <- with(frame, stepCount * 20 - 1)
    frame$averageHops <- with(frame, totalHops/stepCount, 0)
    frame
}

Analysis

Define the dataset years, transmission ranges, host counts, and gradient names and then loop through all posibilities to create the corresponding plots:

years <- c("2013", "2014")
ranges <- c(10,50,100,200,300)
speeds <- c(1,2,5,10,20)
hosts <- c("2013"=8332, "2014"=7069)

for(year in years){
  runnerTimes <- loadRunnerData(year)
  for(range in ranges){
      gradientData <- loadGradientData(year, range, hosts[year])
      gradientSummary <- calculateStableGradients(gradientData)
      plot <- generatePlot(gradientData, year, range)
      plot <- withRunnerPassingRange(plot, runnerTimes)
      plot <- withGradientStability(plot)
      plot <- withTransmissionWindow(plot, gradientSummary, speeds)
      plot <- withGradientFacets(plot)
      print(plot)
   } 
}