---
title: "A directed graph puzzle (Bodycombe 2020)"
subtitle: "Burger run"
author: "Andrew J. Sims"
date: "18th June 2020"
output: 
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 5
    fig_caption: true
    df_print: kable
vignette: >
  %\VignetteIndexEntry{A directed graph puzzle (Bodycombe 2020)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: "REFERENCES.bib"
csl: "nature-no-et-al.csl"
---

```{r}
#| purl = FALSE,
#| include = FALSE
knitr::opts_chunk$set(
  collapse = TRUE,
  echo = FALSE,
  comment = "#>"
)
```

```{r}
#| purl = FALSE
#nolint start
```

```{r}
library(rdecision)
```

```{r}
#| purl = FALSE
#nolint end
```

# Introduction
This puzzle was published in *New Scientist* in June 2020 [@bodycombe2020].
It is a practical example of a problem in graph theory. This vignette explains 
how the puzzle can be solved with `redecison`.

# The puzzle
Three friends agree to drive from A to B via the shortest road possible
(driving down or right at all times). They are hungry, so also want to drive
through a Big Burger restaurant, marked in red. They are arguing about how many
shortest routes will pass through exactly one Big Burger. Xenia: "I reckon there
are 10." Yolanda: "I'd say more like 20." Zara: "No you're both wrong, I bet 
there are more than 50." Who is right, or closest to right?

```{r}
#| results = "hide",
#| fig.keep = "last",
#| fig.align = "center"
# new page
grid::grid.newpage()
# functions to transform coordinates and distances in graph space (0:300)
# to grid space (cm)
fig.size <- dev.size("cm")
scale <- max(300.0 / fig.size[[1L]], 300.0 / fig.size[[2L]])
gx <- function(x) {
  xcm <- fig.size[[1L]] / 2.0 + (x - 150.0) / scale
  return(xcm)
}
gy <- function(y) {
  ycm <- fig.size[[2L]] / 2.0 + (y - 150.0) / scale
  return(ycm)
}
gd <- function(d) {
  dcm <- d / scale
  return(dcm)
}
# grid
for (x in seq(50L, 250L, 50L)) {
  grid::grid.move.to(
    x = grid::unit(gx(x), "cm"), y = grid::unit(gy(50.0), "cm")
  )
  grid::grid.line.to(
    x = grid::unit(gx(x), "cm"), y = grid::unit(gy(250.0), "cm"),
    gp = grid::gpar(lwd = 2.0)
  )
}
for (y in seq(50L, 250L, 50L)) {
  grid::grid.move.to(
    x = grid::unit(gx(50.0), "cm"), y = grid::unit(gy(y), "cm")
  )
  grid::grid.line.to(
    x = grid::unit(gx(250.0), "cm"), y = grid::unit(gy(y), "cm"),
    gp = grid::gpar(lwd = 2.0)
  )
}
grid::grid.text(
  label = "A", x = grid::unit(gx(45.0), "cm"), y = grid::unit(gy(255.0), "cm"),
  gp = grid::gpar(fontsize = 14.0)
)
grid::grid.text(
  label = "B", x = grid::unit(gx(255.0), "cm"), y = grid::unit(gy(45.0), "cm"),
  gp = grid::gpar(fontsize = 14.0)
)
# restaurants
BB <- data.frame(
  x0 = c(150.0, 100.0, 210.0, 160.0, 250.0, 110.0, 50.0),
  y0 = c(60.0, 110.0, 100.0, 150.0, 160.0, 200.0, 210.0),
  x1 = c(150.0, 100.0, 240.0, 190.0, 250.0, 140.0, 50.0),
  y1 = c(90.0, 140.0, 100.0, 150.0, 190.0, 200.0, 240.0)
)
apply(BB, MARGIN = 1L, function(r) {
  grid::grid.move.to(
    x = grid::unit(gx(r[["x0"]]), "cm"), y = grid::unit(gy(r[["y0"]]), "cm")
  )
  grid::grid.line.to(
    x = grid::unit(gx(r[["x1"]]), "cm"),
    y = grid::unit(gy(r[["y1"]]), "cm"),
    gp = grid::gpar(col = "red", lwd = 6.0, lend = "square")
  )
})
```

# Constructing the graph
The grid has 25 nodes and 40 edges (20 horizontal and 20 vertical). These form
a directed graph because it is allowed to drive down or right only. Seven of 
the edges are defined as "Big Burger" edges. Because it is not possible 
to find a path from any node which revisits that node, the graph is
acyclic (a directed acyclic graph, DAG).

Although it possible to construct the graph by creating 25 node 
objects explicitly, it is more compact to create a list of vertices
in a loop construct. Indices $i = [1 .. 5]$ and $j = [1 .. 5]$ are
used to identify grid intersections in the vertical and horizontal directions
respectively. Each node is labelled as $N_{i,j}$ and the index of node
$N_{i,j}$ in the list is $5(i-1)+j$. Similarly, the 40 edges (arrows) are
constructed more compactly in a list, with horizontal edges being labelled
$H_{i,j}$ (the horizontal edge joining node $N_{i,j}$ to node $N_{i,j+1}$) and
the vertical edges similarly as $V_{i,j}$.

```{r}
#| construct-graph,
#| echo = TRUE
# node index function
idx <- function(i, j) {
  return(5L * (i - 1L) + j)
}
# create vertices
N <- vector(mode = "list", length = 5L * 4L)
for (i in seq(5L)) {
  for (j in seq(5L)) {
    N[[idx(i, j)]] <- Node$new(paste0("N", i, j))
  }
}
# create edges
H <- vector(mode = "list", length = 5L * 4L)
ie <- 1L
for (i in seq(5L)) {
  for (j in seq(4L)) {
    a <- Arrow$new(
      N[[idx(i, j)]], N[[idx(i, j + 1L)]], paste0("H", i, j)
    )
    H[[ie]] <- a
    ie <- ie + 1L
  }
}
V <- vector(mode = "list", length = 4L * 5L)
ie <- 1L
for (i in seq(4L)) {
  for (j in seq(5L)) {
    a <- Arrow$new(
      N[[idx(i, j)]], N[[idx(i + 1L, j)]], paste0("V", i, j)
    )
    V[[ie]] <- a
    ie <- ie + 1L
  }
}
# create graph
G <- Digraph$new(V = N, A = c(V, H))
```

```{r}
#| purl = FALSE
# test that graph properties are as expected
stopifnot(
  G$is_simple(),
  !G$is_connected(),
  G$is_weakly_connected(),
  !G$is_tree(),
  !G$is_polytree(),
  G$is_acyclic()
)
```

# Finding the paths
Method `paths` finds all possible paths between any two nodes, where a *path*
is defined as a sequence of distinct and adjacent nodes. Because the
restaurants are specific edges, each path is converted to a *walk*, which
is a path defined as sequence of connected, non-repeating edges. 

In this case, the number of restaurants traversed by each path is counted
by comparing the label associated with each edge in each path with the labels
of the edges which contain a restaurant.

Note that although we cannot guarantee that node A is saved *within* the graph
at index 1 and node B is saved at index 25, we do know that A and B are saved
at indices 1 and 25 in the local list `V`.

```{r}
#| findpaths,
#| echo = TRUE,
#| results = "markdown"
# get all paths from A to B
A <- N[[1L]]
B <- N[[25L]]
P <- G$paths(A, B)
# convert paths to walks
W <- lapply(P, FUN = G$walk)
# count and tabulate how many special edges each walk traverses
BB <- c("V11", "H22", "V25", "H33", "V32", "H44", "V43")
nw <- vapply(W, FUN.VALUE = 1L, FUN = function(w) {
  lv <- vapply(w, FUN.VALUE = TRUE, FUN = function(e) e$label() %in% BB)
  return(sum(lv))
})
# tabulate
ct <- as.data.frame(table(nw))
```

```{r}
#| purl = FALSE
# test that 23 paths traverse one special edge
stopifnot(
  all.equal(ct[[which(ct[, "nw"] == 1L), "Freq"]], 23L)
)
```

# Solution found by `rdecision`
The number of paths which pass through exactly one Big Burger 
is `r ct$Freq[ct$nw == 1L]`. In total there are `r sum(ct$Freq)` paths from A
to B, with the number of restaurants $n$, traversed by each path as follows:

```{r}
#| results = "markdown"
names(ct) <- c("n", "frequency")
ct
```

# Provided solution
Yolanda's estimate is closest - there are 23 shortest routes from A to B
that pass through exactly one Big Burger. One way to solve this kind of 
puzzle is to systematically work from A and keep track of how many
ways there are of reaching each point. With this problem, you should
keep a separate count of how many ways there are of reaching each point
after (a) zero or (b) one Big Burger visits. For line segments that contain
a Big Burger, (b) becomes equal to (a) then becomes equal to 0 with the
old value for (b) effectively discarded.

# References