Trials can be represented with a graph

All potential sequences of trials will be represented as a graph. The graphs nodes will correspond to conditions, and edges between the nodes will correspond to allowable transitions. To represent these graphs, I’ll use the DiagrammeR package.

# library(DiagrammeR)
library(magrittr)
# library(dplyr)

In the graph of my experiment, there will be three nodes for each of the three conditions (Figure 1).

nodes <- DiagrammeR::create_node_df(
  n = 3,
  label = c("zero","low","high"))

DiagrammeR::create_graph(nodes_df = nodes) %>%
  DiagrammeR::render_graph(layout = "tree")

In my experiment, I want trials to go from low to high, zero to high, or high to high, high to low, and high to zero (Figure 2). Including only these five types of transitions means excluding a few of the possible edges that could be in the graph. For example, I do not want any zero contrast trials to follow any other zero contrast trials, nor do I want a low contrast trial to follow a zero contrast trial.

edges <- DiagrammeR::create_edge_df(
  from = c(1,2,3,3,3),
  to = c(3,3,3,1,2)) 

DiagrammeR::create_graph(
  nodes_df = nodes,
  edges_df = edges)  %>%
  DiagrammeR::render_graph(
    layout = "tree")

Constructing a sequence of trials will correspond to walking along the edges, from node to node. That walk will be Eulerian if each edge is be visited exactly once. With so few edges, it’s easy enough to visualize an Eulerian walk through the edges. One possible Eulerian walk (a cycle⁴, even) is shown in Figure 3.

edges_labelled <- DiagrammeR::create_edge_df(
  from = c(3,2,3,3,1),
  to = c(2,3,3,1,3),
  label = as.character(1:5)) 

DiagrammeR::create_graph(
  nodes_df = nodes,
  edges_df = edges_labelled) %>%
  DiagrammeR::render_graph(
    layout = "tree")

The cycle in Figure 3 implies a workable sequence of six trials, but the use of this Eulerian conceptualization will be how it automates creating much longer sequences. For example, to achieve 21 trials the edges could be replicated four times. Figure 4 shows the graph with replicated edges, and already it looks too messy to traverse by sight. A real experiment will involve hundreds of trials, meaning that we’d like a way to automatically traverse an Eulerian circuit.

edges_messy <- DiagrammeR::create_edge_df(
  from = rep(c(3,2,3,3,1), each=4),
  to = rep(c(2,3,3,1,3), each=4)) 

DiagrammeR::create_graph(
  nodes_df = nodes,
  edges_df = edges_messy) %>%
  DiagrammeR::render_graph(layout = "tree")

Hierholzer’s algorithm automates Eulerian cycles

Fortunately, there exists and algorithm for making Eulerian cycles that is both simple to implement and quick to run. First, here is a helper function to replicate edges, replicate_edges.

#' replicate_edges
#'
#' @param edge_df output of DiagrammeR::create_edge_df (will only need columns `to` and `from`)
#' @param n_reps integer number of times that the edges should be replicated
#'
#' @return replicated edge dataframe
replicate_edges <- function(edge_df, n_reps){
  replicate(n_reps, edge_df, simplify = FALSE) %>%
    dplyr::bind_rows() %>%
    dplyr::mutate(id = 1:dplyr::n())
}

The next function will generate the Eulerian circuit, walk_circuit. It will take in an edge dataframe (possibly replicated) and output a vector containing the nodes listed in the order that they were reached. Again, I won’t spend too long explaining why this works. But the basic idea is to traverse the edges, deleting edges as you walk along them. You’ll eventually reach a dead-end. If there are still more edges, then backtrack until you can travel along an edge that will result in a different dead-end. Save a list of the nodes that were traveled while backtracking, and these nodes will contain the circuit.

#' walk_circuit
#'
#' @param edge_df edge dataframes 
#' @param curr_v vertex at which to start the circuit
#'
#' @return vector consisting of Eulerian circuit along edges
#'
#' @details modified python script from https://gregorulm.com/finding-an-eulerian-path/.
walk_circuit <- function(edge_df, curr_v){
  
  # helpful to have the edges stored by node
  adj <- edge_df %>%
    dplyr::group_split(from)
  
  # vector to store final circuit 
  circuit <- c() 
  # Maintain a stack to keep vertices 
  # start from given node
  curr_path <- curr_v 
  while (length(curr_path)){
    # If there's a remaining edge 
    if (nrow(adj[[curr_v]])){
      # Push the vertex 
      curr_path <- c(curr_path,curr_v)   
      # Find the next vertex using an edge 
      next_v_ind <- sample.int(nrow(adj[[curr_v]]), size=1)
      next_v <- adj[[curr_v]]$to[next_v_ind]
      # and remove that edge 
      adj[[curr_v]] <- adj[[curr_v]][-next_v_ind,]
      # Move to next vertex 
      curr_v <- next_v 
    } else{ # back-track to find remaining circuit 
      circuit <- c(circuit, curr_v)   
      # Back-tracking 
      curr_v <- tail(curr_path, n = 1) 
      curr_path <- head(curr_path, n = -1)   
    }
  }   
  return(rev(circuit))
}

Now replicate the edges twice and go for and Eulerian tour.

edges_twice <- replicate_edges(edges, 2)
walk_circuit(edges_twice, 3)

##  [1] 3 1 3 2 3 1 3 3 2 3 3

This sequence is small enough that it’s feasible to verify the Eulerian property by hand, but it’ll be nice to have automate the checking. That is the purpose of this next function, check_blocking.

#' check_blocking
#'
#' @param circuit output of walk_circuit
#' @param nodes nodes_df, output of DiagrammeR::create_node_df. Used to label which nodes were visited during the walk
#'
#' @return tbl containing the counts of each transition type contained in the circuit. 
#' If all went well, the counts should be equal
check_blocking <- function(circuit, nodes){
  tibble::tibble(contrast = circuit, .name_repair = "check_unique") %>%
    dplyr::mutate(
      trial = 1:dplyr::n(),
      contrast = nodes$label[contrast],
      last_contrast = dplyr::lag(contrast)) %>%
    dplyr::filter(trial > 1) %>%
    dplyr::group_by(contrast, last_contrast) %>%
    dplyr::summarise(n = dplyr::n(), .groups = "drop")
}

Now, generate a sequence of 101 trials,

edges_large <- replicate_edges(edges, n_reps = 20)
circuit <- walk_circuit(edges_large, 3)
circuit

##   [1] 3 2 3 3 1 3 2 3 2 3 3 1 3 1 3 3 1 3 3 2 3 1 3 2 3 2 3 1 3 3 1 3 3 2 3 3 3
##  [38] 1 3 1 3 1 3 3 2 3 2 3 3 2 3 2 3 1 3 1 3 3 2 3 1 3 1 3 3 2 3 1 3 3 3 3 3 2
##  [75] 3 3 1 3 2 3 3 2 3 1 3 1 3 3 1 3 2 3 2 3 1 3 3 2 3 2 3

and check that each transition happened equally often

check_blocking(circuit, nodes) %>%
  knitr::kable()