Tea & Stats

Increases

To count the number of times elements are increasing in a vector it’s as simple as

depths <- c(199, 200, 208, 210, 200, 207, 240, 269, 260, 263)
sum(diff(depths) > 0)

## [1] 7

for which I defined a function called increases.

Rolling sum

For part two, we first want to calculate the three-depth moving sum, then we count the increases as in part one. There are plenty of solutions in external R packages for getting lagged (and leading) vectors, for instance dplyr::lag() and dplyr::lead():

depths + dplyr::lead(depths) + dplyr::lead(depths, 2)

##  [1] 607 618 618 617 647 716 769 792  NA  NA

Or you could even calculate the rolling sum using a pre-made solution in zoo (Z’s Ordered Observations, a time-series analysis package).

zoo::rollsum(depths, 3)

## [1] 607 618 618 617 647 716 769 792

To avoid loading any external packages at this stage, I defined my own base R function called rolling_sum(), which uses tail and head with negative lengths to omit the first and last elements of the vector:

head(depths, -2) + head(tail(depths, -1), -1) + tail(depths, -2)

## [1] 607 618 618 617 647 716 769 792

As David Schoch points out, you can just use the lag argument of diff to make this entire puzzle into a one-liner:

sapply(c(1, 3), \(lag) sum(diff(depths, lag) > 0))

## [1] 7 5

Depth sum

Read in the input as a two-column data frame using read.table(). I gave mine nice column names, cmd and value, but this isn’t essential.

Then take advantage of the fact that TRUE == 1 and FALSE == 0 to make a mathematical ifelse-type statement for the horizontal and vertical movements. In my R package, this is implemented as a function called dive():

x <- (cmd == 'forward') * value
y <- ((cmd == 'down') - (cmd == 'up')) * value
sum(x) * sum(y)

Cumulative depth sum

Part two is much like part one, but now y represents (change in) aim, and (change in) depth is derived from that. Don’t forget the function cumsum(), which can save you writing a loop! Here is the body of my function dive2():

x <- (cmd == 'forward') * value
y <- ((cmd == 'down') - (cmd == 'up')) * value
depth <- cumsum(y) * x
sum(x) * sum(depth)

Power consumption

There are a few different ways you could approach part one, but my approach was first to read in the data as a data frame of binary integers using the function read.fwf(). Then, find the most common value in each column using the base function colMeans() and rounding the result.

According to the instructions, in the event of a tie you should take 1 to be the most common digit. Although this is familiar to real life—0.5 rounds up to 1—computers don’t work this way: R rounds to even instead (see ?round). Because zero is even, that means round(0.5) yields 0. To get around this, add 1 before rounding, then subtract it again.

My function power_consumption(), which once again takes advantage of TRUE being equivalent to 1 and FALSE to 0:

common <- round(colMeans(x) + 1) - 1
binary_to_int(common) * binary_to_int(!common)

To convert a vector of binary digits to decimal, I use the following utility function:

binary_to_int <- function(x) {
  sum(x * 2 ^ rev(seq_along(x) - 1))
}

However, if using a string representation then there’s a handy function in base R called strtoi() that you could also use for this (thanks to Riinu Pius for that tip).

Life support

Part two finds the common digits in a successively decreasing set of binary numbers. A loop is appropriate here, since we can halt once there is only one number left. As this loop will only run (at most) 12 times in total, it shouldn’t be too slow in R.

Function life_support():

life_support <- function(x) {
  oxygen <- co2 <- x
  for (j in 1:ncol(x)) {
    if (nrow(oxygen) > 1) {
      common <- most_common(oxygen)
      oxygen <- oxygen[oxygen[, j] == common[j], ]
    }
    if (nrow(co2) > 1) {
      common <- most_common(co2)
      co2 <- co2[co2[, j] != common[j], ]
    }
  }
  binary_to_int(oxygen) * binary_to_int(co2)
}

There might be cleverer ways of doing this.

Bingo

This is one of those problems where half the battle is working out which data structure to use. I wrote a function read_draws() that reads in the first line of the file to get the drawn numbers, then separately reads in the remainder of the file to get the bingo cards stacked as a data frame. Later we take advantage of the fact that the bingo cards are square to split the data frame into a list of matrices.

read_draws <- function(file) {
  draws <- scan(file, sep = ',', nlines = 1, quiet = TRUE)
  cards <- read.table(file, skip = 1)
  list(draws = draws, cards = cards)
}

As numbers are called out, I replace them in the dataset with NAs. Then the helper score_card() counts the number of NAs in each row and column. If there are not enough, we return zero, else we calculate the score.

score_card <- function(mat, draw) {
  marked <- is.na(mat)
  if (all(c(rowMeans(marked), colMeans(marked)) != 1))
    return(0)
  sum(mat, na.rm = TRUE) * draw
}

Then we put it all together, looping through the draws, replacing numbers with NA and halting as soon as someone wins. Function play_bingo() is defined as follows, using just base R commands:

play_bingo <- function(draws, cards) {
  size <- ncol(cards)
  ncards <- nrow(cards) / size
  ids <- rep(1:ncards, each = size)

  for (d in draws) {
    cards[cards == d] <- NA
    score <- sapply(split(cards, ids), score_card, draw = d)
    if (any(score > 0))
      return(score[score > 0])
  }
}

Last caller

Part two is very similar, but we throw away each winning bingo card as we go to avoid redundant computation, eventually returning the score when there is only one left. Here is function play_bingo2(), which uses the same two utility functions:

play_bingo2 <- function(draws, cards) {
  size <- ncol(cards)

  for (d in draws) {
    ncards <- nrow(cards) / size
    ids <- rep(1:ncards, each = size)
    cards[cards == d] <- NA
    score <- sapply(split(cards, ids), score_card, draw = d)
    if (any(score > 0)) {
      if (ncards == 1)
        return(score[score > 0])
      cards <- cards[ids %in% which(score == 0), ]
    }
  }
}

Further optimizations are possible. For example: as written, we calculate every intermediate winner’s score, but we only really need to do it for the first (part 1) and last (part 2) winners.

Also, we could draw more than one number at a time, as we know that nobody’s going to win until at least the fifth draw (for 5×5 cards) and from there, increment according to the minimum number of unmarked numbers on any row or column.

I didn’t bother implementing either of these, as it already runs quickly enough.

Median

While it’s possible to brute-force this puzzle by simply calculating the fuel requirement at every single point (within the range of the inputs), you can do it about 200× faster by treating it as an optimization problem.

The total fuel required for any potential position is

x <- scan('input.txt', sep = ',')
f <- function(pos) sum(abs(x - pos))

where x are the initial locations of the crabs. Then run it through optimize(), and round to the nearest integer position:

sol <- optimize(f, range(x))$minimum
f(round(sol))

However, there is an even faster analytical solution!

sol <- median(x)

Thanks to Claire Little for pointing this out.

Mean

Part two just has a slightly different function to optimize. Using the formula for the sum of an arithmetic progression:

f2 <- function(pos) {
  n <- abs(x - pos)
  sum(n / 2 * (1 + n))
}

Then we can simply minimize this function as before.

sol <- optimize(f2, range(x))$minimum
f2(round(sol))

However, there’s a shortcut for this part as well! Calculate the mean of the initial positions, and work out which of the two nearest integers gives the minimum result:

min(
  f2(floor(mean(x))),
  f2(ceiling(mean(x)))
)

Thanks to Jonatan Pallesen. This is about 5 times faster than my optimizer.

And here is what the functions look like for my input dataset:

Unique digits

Read in the data and then the first part is just a one-liner:

input <- do.call(rbind, strsplit(readLines(input_file(8)), '[^a-z]+'))

count_unique <- function(x) {
  sum(nchar(x[, -(1:10)]) %in% c(2, 3, 4, 7))
}

Segment matching

I really wanted to solve part two using graph theory, by representing the puzzle as a maximum bipartite matching problem. However, I couldn’t quite get this to work. My final solution is instead just a lot of leg work.

Essentially you solve the problem by hand and then encode the process programmatically. Recognize that some digits have segments in common, or not in common, and use this to eliminate the possibilities. I stored the solutions in a named vector, which I was able to use to look up the digits found so far.

The function setdiff() comes in useful.

contains <- function(strings, letters) {
  vapply(strsplit(strings, ''),
         function(s) all(strsplit(letters, '')[[1]] %in% s),
         logical(1))
}

output_value <- function(vec) {
  segments <- c('abcefg', 'cf', 'acdeg', 'acdfg', 'bcdf',
                'abdfg', 'abdefg', 'acf', 'abcdefg', 'abcdfg')
  nchars <- setNames(nchar(segments), 0:9)

  # Sort the strings
  vec <- sapply(strsplit(vec, ''), function(d) paste(sort(d), collapse = ''))
  sgn <- head(vec, 10)
  out <- tail(vec, 4)

  # Store the known values
  digits <- setNames(character(10), 0:9)
  unique <- c('1', '4', '7', '8')
  digits[unique] <- sgn[match(nchars[unique], nchar(sgn))]

  # Remaining digits have 5 or 6 segments:
  sgn <- setdiff(sgn, digits)
  digits['3'] <- sgn[nchar(sgn) == 5 & contains(sgn, digits['1'])]
  digits['6'] <- sgn[nchar(sgn) == 6 & !contains(sgn, digits['1'])]
  sgn <- setdiff(sgn, digits)
  digits['0'] <- sgn[nchar(sgn) == 6 & !contains(sgn, digits['4'])]
  sgn <- setdiff(sgn, digits)
  digits['9'] <- sgn[nchar(sgn) == 6]
  sgn <- setdiff(sgn, digits)
  digits['2'] <- sgn[
    contains(sgn, do.call(setdiff,
                          unname(strsplit(digits[c('8', '6')], ''))))
  ]
  digits['5'] <- setdiff(sgn, digits)

  # Combine four output digits:
  as.numeric(paste(match(out, digits) - 1, collapse = ''))
}

Lowest points

You can find all the lowest points with a one-liner:

lowest <- function(h) {
  h < cbind(h, Inf)[, -1] &          # right
    h < rbind(h, Inf)[-1, ] &        # down
    h < cbind(Inf, h[, -ncol(h)]) &  # left
    h < rbind(Inf, h[-nrow(h), ])    # up
}

Then do sum(h[lowest(h)]) to get the result, where h is a numeric matrix of the input data.

Basins

The second part is harder and doesn’t immediately lead from the first. Initially I thought of replacing each lowest point with Inf, then finding the new lowest points and repeating the process until all the basins are found. However, the basins are simply all points where the height is < 9, so you can find the basins in a single step.

The tricky part is labelling them separately, so you can count up their respective sizes.

The boring way of doing this is just to loop over the indices and label the points that neighbour already-labelled ones (starting with the lowest points as the initial labels), doing several passes until everything (except the 9s) is labelled.

basins <- function(h) {
  l <- lowest(h)
  h[] <- ifelse(h < 9, NA, Inf)
  h[l] <- 1:sum(l)
  while (anyNA(h)) {
    for (i in 1:nrow(h)) for (j in 1:ncol(h)) {
      if (is.na(h[i, j])) {
        nbrs <- h[cbind(c(max(i - 1, 1), min(i + 1, nrow(h)), i, i),
                        c(j, j, max(j - 1, 1), min(j + 1, ncol(h))))]
        if (any(is.finite(nbrs)))
          h[i, j] <- nbrs[is.finite(nbrs)][1]
      }
    }
  }
  sizes <- table(h[is.finite(h)])
  head(sort(sizes, decreasing = TRUE), 3)
}

To vectorize this in the same way as part one, we define a new binary (infix) operator %c%, analogous to dplyr::coalesce(). What this does is replace an NA value (a basin not yet assigned a label) with its finite neighbour, while leaving Infs (marking basin edges) alone.

"%c%" <- function(x, y) {
  ifelse(is.infinite(x), x,
         ifelse(!is.na(x), x,
                ifelse(!is.infinite(y), y, x)))
}

Then the new function for part two is as follows. It is five times faster to run than the nested loop above.

basins2 <- function(h) {
  l <- lowest(h)
  h[] <- ifelse(h < 9, NA, Inf)
  h[l] <- 1:sum(l)
  while(anyNA(h)) {
    h <- h %c%
      cbind(h, NA)[, -1] %c%        # right
      rbind(h, NA)[-1, ] %c%        # down
      cbind(NA, h[, -ncol(h)]) %c%  # left
      rbind(NA, h[-nrow(h), ])      # up
  }
  sizes <- table(h[is.finite(h)])
  head(sort(sizes, decreasing = TRUE), 3)
}

You can also formulate this as an image analysis problem, effectively treating each basin as an area of similar colour to select, or you can treat it as a network theory problem and apply the igraph package to find graph components.

Corrupt characters

Whilst it’s probably possible to do the first part with some very fancy recursive regular expressions, I don’t know how to use them.

Instead, my method of finding unmatched brackets is simply to search for empty pairs of brackets and successively strip them from the string. Keep doing this until the strings stop changing. Then, get the first closing bracket (if any), using regmatches(). These are the illegal characters.

My function syntax_score() is implemented as follows:

lines <- readLines('input.txt')
old <- ''
while (!identical(old, lines -> old))
  lines <- gsub(r'(\(\)|<>|\{\}|\[\])', '', lines)
illegals <- regmatches(lines, regexpr(r'(\)|>|\}|\])', lines))

The syntax score is calculated using a named vector as a lookup table.

illegal_score <- c(')' = 3, ']' = 57, '}' = 1197, '>' = 25137)
sum(illegal_score[illegals])

Autocomplete

Part two starts the same, but instead of extracting the illegal characters we just throw away those lines that contain them.

illegals <- grep(r'(\)|>|\}|\])', lines)
chars <- strsplit(lines[-illegals], '')

From here, we can calculate the scores using a Reduce operation (from right to left) with another lookup table. The final answer is the median score.

complete_score <- c('(' = 1, '[' = 2, '{' = 3, '<' = 4)
scores <- sapply(chars, Reduce, init = 0, right = TRUE,
                 f = \(c, s) 5 * s + complete_score[c])
median(scores)

The function autocomplete() wraps it all together.

Convoluted octopuses

The process of updating the energy levels can be described using a convolution matrix. It’s easy—like on Day 11 last year—to use a ready-made solution from an image analysis package for this, namely OpenImageR::convolution().

The convolution matrix, or kernel is \[\begin{bmatrix}1 & 1 & 1 \\1 & 0 & 1 \\1 & 1 & 1\end{bmatrix},\] to be applied on an indicator matrix of ‘flashing’ octopuses, and then added to the result. In R,

kernel <- matrix(1:9 != 5, 3, 3)

So we define a function, step1(), that applies a single step of the energy level updating process. Since each octopus can only flash once in a given step, we keep track of those that have already flashed, as well as those currently flashing. A short while() loop repeats until no more octopuses flash.

step1 <- function(x) {
  x <- x + 1
  flashing <- flashed <- x == 10
  while (any(flashing)) {
    x <- x + OpenImageR::convolution(flashing, kernel)
    flashing <- x > 9 & !flashed
    flashed <- flashing | flashed
  }
  x[x > 9] <- 0
  x
}

However, there is a base R alternative to OpenImageR::convolution() that you can substitute in, for negligible speed penalty (despite that package being written in C).

add_neighbours <- function(x) {
  I <- nrow(x)
  J <- ncol(x)
  cbind(x[, -1], 0) +  # Right
  rbind(x[-1, ], 0) +  # Down
  cbind(0, x[, -J]) +  # Left
  rbind(0, x[-I, ]) +  # Up
  rbind(cbind(x[-1, -1], 0), 0) + # SE
  rbind(0, cbind(x[-I, -1], 0)) + # NE
  rbind(cbind(0, x[-1, -J]), 0) + # SW
  rbind(0, cbind(0, x[-I, -J]))   # NW
}

Counting flashes

We then solve both parts one and two with a single function that repeatedly applies the updating step and counts the flashes (number of zeros) each time.

By default, it simply completes iter steps and then returns cumulative total number of flashes.

In part2 mode, the function will terminate early if it encounters a step where all the octopuses are flashing and return the iteration number, otherwise it will generate a warning.

count_flashes <- function(x, iter = 100, part2 = FALSE) {
  count <- 0
  for (i in 1:iter) {
    x <- step1(x)
    nflashes <- sum(x == 0)
    if (part2 & nflashes == prod(dim(x)))
      return(i)
    count <- count + nflashes
  }
  if (part2)
    warning('No synchronization detected in ', iter, ' steps')
  count
}

The whole thing runs in about 20 milliseconds on my input dataset, or 60 milliseconds using the base R convolution function.

Counting paths

Today’s puzzle is very obviously a graph theory problem, so it would seem intuitive to break out igraph for this. However, the functions in that package are mainly designed for finding the shortest paths from A to B, not a complete list of all possible paths.

Therefore, it’s easier (as far as I can tell) to solve this with your own data structure and a recursive function or two.

Antoine Fabri suggested a rather neat solution of describing the graph with a named list, where the names of the elements describe a node, and the vector contained therein list the other nodes reachable from that point.

We can read in the data as follows:

read_edges <- function(file) {
  edges <- read.table(file, sep = '-', col.names = c('from', 'to'))
  edges <- rbind(edges, setNames(edges, rev(colnames(edges))))
  edges <- edges[edges$to != 'start', ]
  split(edges$to, edges$from)
}

Then we solve part one with a recursive function that, starting from start, traverses a path all the way to end and increments a counter. To avoid revisiting small caves, we discard any edges to vertices with lowercase names, with the aid of tolower().

count_paths <- function(edgelist, node = 'start') {
  if (node == 'end')
    return(1)
  if (!length(edgelist[[node]]))
    return(0)
  if (node == tolower(node))
    edgelist <- lapply(edgelist, \(v) v[node != v])
  sum(vapply(edgelist[[node]], count_paths, e = edgelist,
             FUN.VALUE = numeric(1)))
}

Revisiting caves

Part two is just a slight adaptation. This time we can’t simply delete small caves as soon as we visit them, so we need to keep track of those visited so far, while keeping the graph intact. Once we find ourselves on a node that we have already visited, it means we are visiting it for the second time, so at this point we can delete all the visited nodes and switch over to the function from part one.

count_paths2 <- function(edgelist, node = 'start', visited = NULL) {
  if (node == 'end')
    return(1)
  if (node == tolower(node)) {
    if (node %in% visited) {
      edgelist <- lapply(edgelist, \(v) v[!v %in% visited])
      return(count_paths(edgelist, node))
    }
    visited <- union(visited, node)
  }
  if (!length(edgelist[[node]]))
    return(0)
  sum(vapply(edgelist[[node]], count_paths2, e = edgelist, visited,
             FUN.VALUE = numeric(1)))
}

R is quite slow at recursion, so this second function takes about 6 seconds to run on my full input data. It may be possible to optimize this further, but that might require the use of igraph, which is written in C, rather than a base R solution.

Folding paper

The hardest bit about this puzzle is reading in the data, since we have two different data structures in the same file.

The neatest way of doing this is to find the separator (a blank line) in the middle of the file, then read through the parts of the file before and after this separately. Since there’s only one file to read, you could of course just inspect it in a text editor and find out this line number by hand. But here’s a programmatic method.

n <- which(readLines('input.txt') == '')
paper <- read.table('input.txt', sep = ',', nrows = n - 1,
                    col.names = c('x', 'y'))
folds <- read.table('input.txt', sep = '=', skip = n)

Next, we need a function that’ll fold the paper in half once. My initial solution, which works perfectly well, is to convert the x and y coordinates to a matrix representation, then simply update the matrix using subset operations.

mat <- with(paper, matrix(0, max(x) + 1, max(y) + 1))
mat[data.matrix(paper) + 1] <- 1

Watch out for R’s one-based indexing. We want to fold along the coordinate f, which is at index f + 1, and so the indices either side are f and f + 2.

fold_once <- function(x, dir, f) {
  if (dir == 'fold along y') {
    x[, 1:f] | x[, ncol(x):(f + 2)]
  } else {
    x[1:f, ] | x[nrow(x):(f + 2), ]
  }
}

At this point you might notice that the paper is always being folded exactly in half every time. So the value after the = in the input is actually redundant; it’s always the middle of the matrix. For instance, in the example data, the matrix has dimension 15×11 and we fold along \(y=7\) (i.e. index 8), which is half of 15. The next instruction \(x=5\) (index 6) is half of the other dimension.

From here, loop over the directions and then plot the final matrix using image().

Tidy origami

Why do they give us the values to fold along if it’s always the middle value? Because it’s not actually very efficient to store these sparse coordinates in a matrix. Instead, recognize, as pointed out by Riinu Pius and Antoine Fabri, that you can just keep the coordinates in a long format and fold them as follows:

fold_once <- function(x, d, f) {
  x[, d] <- ifelse(x[, d] >= f, 2 * f - x[, d], x[, d])
  x[!duplicated(x), ]
}

Then loop over all the instructions with:

fold_paper <- function(x, folds, n = nrow(folds)) {
  for (i in 1:n)
    x <- fold_once(x, folds[i, 1], folds[i, 2])
  x
}

And finally you can plot the result using ggplot2.

library(adventofcode2021)
input <- read_origami(input_file(13))
input$n <- nrow(input[[2]])
folded <- do.call(fold_paper, input)

library(ggplot2)
ggplot(folded) +
  aes(x, y) +
  geom_tile(fill = '#f5c966') +
  coord_fixed() +
  scale_y_reverse() +
  theme_void()

Inserting letters

The initial approach to this problem is to perform the insertions and build up an ever-lengthening string. There might be a way to do this using regular expressions, but I used an approach based on table joins.

First, read in the data:

template <- strsplit(readLines('input.txt', n = 1), '')[[1]]
template <- data.frame(first = head(template, -1),
                       second = tail(template, -1))
rules <- read.table('input.txt', skip = 1)
rules <- data.frame(first  = substr(rules[, 1], 1, 1),
                    second = substr(rules[, 1], 2, 2),
                    insert = rules[, 3])

This should give you two data frames that look something like this:

  first second
1     N      N
2     N      C
3     C      B

   first second insert
1      C      H      B
2      H      H      N
3      C      B      H
4      N      H      C
5      H      B      C
6      H      C      B
7      H      N      C
8      N      N      C
9      B      H      H
10     N      C      B
11     N      B      B
12     B      N      B
13     B      B      N
14     B      C      B
15     C      C      N
16     C      N      C

From here, the method is straightforward. Join the template with the rules and produce a new template comprised of the pairs first -> insert and insert -> second.

template <- with(merge(template, rules),
                 data.frame(first = c(first, insert),
                            second = c(insert, second)))

Repeat this process the required number of times and we will have a data frame containing all the pairs of letters.

There’s a problem, though: merge() does not preserve row order, so we’ve lost track of which letters are the first and the last in the sequence. This means counting up the letters in either column first or second is going to undercount by 1. We could just refer back to the original input, since the first and last letters won’t change.

However, there’s another neat solution, which takes advantage of the fact that all the middle characters will appear in both columns, but the first and last values, if different, will appear an odd number of times.

counts <- (table(unlist(template)) + 1) %/% 2

From here, the difference between the minimum and maximum counts is given by:

diff(range(counts))

Scalable insertion

The above approach works for part one, but once we increase the iteration count to 40, the data frame or string representation is far too long to store in memory. It simply doesn’t scale.

Rather than keep the entire sequence and then count up the letters at the end, we can instead store pair frequencies. To do this, we can call dplyr::count() or use a base R approach.

Add a frequency column to the original data frame:

template <- data.frame(first  = head(template, -1),
                       second = tail(template, -1),
                       n = 1)

or if you want to aggregate on initialization, you could use:

template <- as.data.frame(table(template), stringsAsFactors = FALSE)

Then after each round of insertions, aggregate() the counts:

data <- with(merge(template, rules),
             data.frame(first = c(first, insert),
                        second = c(insert, second)))
template <- aggregate(n ~ first + second, data, sum)

Finally, the total letter counts are given by:

tab <- xtabs(n ~ first + second, template)
(rowSums(tab) + colSums(tab) + 1) %/% 2

Versions

This puzzle takes a long time to read. The solution involves recursive functions that read through the packets and subpackets in turn. My approach takes the bits and throws them away as they are processed, until there are none left (or, for sub-packets, if we have finished parsing the requisite number).

We are working between bases here, so you want a couple of utility functions to convert to and from integers.

to_bits <- function(x) {
  ints <- strtoi(strsplit(x, '')[[1]], 16)
  bits <- sapply(ints, \(n) tail(rev(as.numeric(intToBits(n))), 4))
  as.vector(bits)
}

Unfortunately, strtoi is not your friend for the binary-to-integer conversion, because it will almost surely silently introduce integer overflow on your main input data. Here’s one I prepared earlier that uses floating point numbers so it won’t overflow.

binary_to_int <- function(x) {
  sum(x * 2 ^ rev(seq_along(x) - 1))
}

Now here’s the main recursion.

packet_versions <- function(bits, depth = 0, rem_subs = Inf, eval = FALSE) {
  acc <- 0
  while (any(bits > 0) & rem_subs > 0) {
    rem_subs <- rem_subs - 1
    version <- binary_to_int(bits[1:3])
    acc <- acc + version
    type <- binary_to_int(bits[4:6])
    bits <- tail(bits, -6)
    if (type == 4) { # literal value
      n_groups <- which.min(bits[seq(1, length(bits), 5)])
      sub <- head(bits, n_groups * 5)
      sub <- sub[(seq_along(sub) - 1) %% 5 > 0]
      # literal_value <- binary_to_int(sub)
      bits <- tail(bits, -n_groups * 5)
      next
    }
    # Operator mode:
    lentype <- bits[1]
    bits <- bits[-1]
    if (lentype == 0) {
      sub_length <- binary_to_int(bits[1:15])
      bits <- tail(bits, -15)
      sub  <- head(bits, sub_length)
      acc <- acc + packet_versions(sub, depth + 1)[['acc']]
      bits <- tail(bits, -sub_length)
    } else {
      n_subs <- binary_to_int(bits[1:11])
      bits <- tail(bits, -11)
      sub_result <- packet_versions(bits, depth + 1, n_subs)
      acc <- acc + sub_result[['acc']]
      bits <- tail(bits, sub_result[['length']])
    }
  }
  if (depth)
    acc <- c(acc = acc, length = length(bits))
  acc
}

Be careful with length type ID 1, because you need to save the number of bits you’ve parsed and then send this number back to the main function. Similarly, it’s a good idea to keep track of how deep in the recursion you are at any given time.

Evaluation

The second part is tricky. I decided to adapt my function from part one to generate syntactically-valid R expressions, so that final step is simply evaluating machine-generated code.

packet_decode <- function(bits, depth = 0, rem_subs = Inf) {
  acc <- NULL
  while (any(bits > 0) & rem_subs > 0) {
    rem_subs <- rem_subs - 1
    type <- binary_to_int(bits[4:6])
    bits <- tail(bits, -6)
    op <- c('sum', 'prod', 'min', 'max', 'c', '`>`', '`<`', '`==`')[type + 1]
    if (type == 4) { # literal value
      n_groups <- which.min(bits[seq(1, length(bits), 5)])
      sub <- head(bits, n_groups * 5)
      sub <- sub[(seq_along(sub) - 1) %% 5 > 0]
      literal_value <- binary_to_int(sub)
      acc$expr <- c(acc$expr, literal_value)
      bits <- tail(bits, -n_groups * 5)
      next
    }
    # Operator mode:
    lentype <- bits[1]
    bits <- bits[-1]
    if (lentype == 0) {
      sub_length <- binary_to_int(bits[1:15])
      bits <- tail(bits, -15)
      sub  <- head(bits, sub_length)
      sub_result <- packet_decode(sub, depth + 1)
      subexpr <- c(op, sub_result[['expr']])
      acc$expr <- c(acc$expr, list(subexpr))
      bits <- tail(bits, -sub_length)
    } else {
      n_subs <- binary_to_int(bits[1:11])
      bits <- tail(bits, -11)
      sub_result <- packet_decode(bits, depth + 1, n_subs)
      subexpr <- c(op, sub_result[['expr']])
      acc$expr <- c(acc$expr, list(subexpr))
      bits <- tail(bits, sub_result[['length']])
    }
  }
  if (depth)
    acc$length <- length(bits)
  acc
}

The structure of the function is essentially the same, and the result is a deeply nested list representing the hierarchical tree of expressions. From here we want to convert it from a tree into some code, and the function combine_expr does this for us:

combine_expr <- function(node, eval = FALSE) {
  if (length(node) == 1)
    return(node)
  expr <- sprintf('%s(%s)', node[1], paste(node[-1], collapse = ', '))
  if (!eval)
    return(expr)
  eval(str2expression(expr))
}

However, it will only work on leaves of the tree. How can we correctly operate on the entire tree and turn it into a single expression? That requires more recursion. Annoyingly, the function rapply doesn’t quite do what’s needed, so we introduce yet another function, called packet_parse:

packet_parse <- function(tree, eval = FALSE) {
  if (!is.list(tree))
    return(combine_expr(tree, eval))
  expr <- packet_parse(unlist(lapply(tree, packet_parse, eval), use.names = F))
  if (eval & is.character(expr))
    return(eval(parse(text = expr)))
  expr
}

You can choose to view the expression, or evaluate it. Here is my final expression generated by packet_parse:

sum(prod(425542, `<`(247, 247)), sum(121, 21236), prod(`>`(sum(11, 12, 11), sum(7, 10, 7)), 32566), prod(`<`(sum(8, 7, 15), sum(6, 11, 10)), 4507180), min(prod(prod(sum(prod(max(prod(min(prod(min(prod(max(sum(sum(max(sum(sum(min(prod(prod(130)))))))))))))))))))), 139930778832, prod(`>`(52, 667118), 10), 602147, max(62199), prod(14849899, `<`(11716, 26963)), prod(4083, `>`(135, 135)), prod(135, 217, 224), 73, prod(sum(13, 4, 9), sum(12, 15, 7), sum(13, 10, 9)), min(194), prod(182, 197, 136, 2, 242), prod(226, 142, 34, 124), max(4025, 186042), min(30059, 126119002), min(9, 260, 162), prod(`<`(4, 4), 28699), prod(1945, `==`(1714, 1714)), prod(7, `<`(1545, 108)), sum(12), prod(200, `>`(31050, 655605)), 3154, prod(3, `<`(64896, 116)), 3055, prod(13), min(48082, 226938, 1175, 68077774919), sum(66, 15, 181, 1380642642, 11831587), prod(241, 59), prod(150, `>`(2742, 113)), 37007908601, max(52444, 11, 13008816, 2935), 20723, 8, prod(5, `>`(6241732, 759708)), sum(prod(15, 7, 4), prod(14, 2, 12), prod(13, 6, 6)), sum(2877, 229333, 655820, 1020971), sum(39581, 2, 14), max(982557, 44, 31), 68, prod(`==`(11530, 3492), 41177), prod(`==`(236, 918711093), 3937), max(903466, 228, 6, 25989131, 4028), 229, min(299875, 10969849, 11481, 2281, 13), prod(55300721, `>`(63, 63)), prod(244, `>`(sum(7, 13, 7), sum(12, 5, 14))), prod(4494263, `==`(sum(4, 15, 4), sum(3, 3, 14))), prod(`<`(45, 3307915), 58514), prod(3596530693, `<`(sum(3, 12, 4), sum(9, 11, 2))))

Which yields the answer 184487454837 when evaluated.

Reduction

I decided (for better or for worse) to convert the square brackets into list(...) expressions, then parse the input as nested lists.

read_sf <- function(input) {
    out <- lapply(input, function(string) {
        expr <- gsub('\\]', ')', gsub('\\[', 'list(', string))
        eval(parse(text = expr))
    })
    if (length(out) == 1) return(unlist(out, FALSE)) else out
}

Then we need to perform the reductions on our list structures. Splitting is easy, but it’s important that we split only the first node we find and then stop. I used a counter called done to ensure the rapply loop terminates early.

sf_split <- function(lst) {
    done <- FALSE
    rapply(lst, \(n) {
        if (n <= 9 | done) return(n)
        done <<- TRUE
        list(floor(n / 2), ceiling(n / 2))
    }, how = 'replace')
}

To explode, it’s a bit trickier because you need several nodes to interact with one another. I wrote some recursive functions to find the depth of the list, to mark the appropriate node with an NA and then to replace it.

I also discovered the relist function, which makes it possible to flatten a list, do some vector operations (handy for getting the indices just before and after our exploded pair) and then return the data to the original nested structure.

depths <- function(lst, depth = 0) {
    if (!is.list(lst)) return(depth)
    unlist(sapply(lst, depths, depth + 1))
}

replaceNA <- function(lst, replacement) {
    if (all(is.na(lst))) return(replacement)
    if (!is.list(lst)) return(lst)
    lapply(lst, replaceNA, replacement)
}

sf_explode <- function(lst) {
    i <- which(depths(lst) > 4)
    L <- i[1]; R <- i[2]
    x <- unlist(lst)
    if (L > 1)
        x[L - 1] <- x[L - 1] + x[L]
    if (R < length(x))
        x[R + 1] <- x[R] + x[R + 1]
    x[c(L, R)] <- NA
    replaceNA(relist(x, lst), 0)
}

Then we put it all together as reduce operations.

sf_reduce <- function(lst) {
    Reduce(\(a, b) list(sf_reduce1(a), b), lst) |> sf_reduce1()
}

sf_reduce1 <- function(lst) {
    if (max(depths(lst)) > 4) lst <- Recall(sf_explode(lst))
    if (any(unlist(lst) > 9)) lst <- Recall(sf_split(lst))
    lst
}

If you try to reduce the whole list in one operation, you’ll definitely get a stack overflow with the larger input datasets. Thus it’s important to reduce from the left until it’s as simple as possible, then combine it with the next line.

My code is quite slow; it might be more sensible to rewrite this as a repeat loop rather than doing all this recursion.

Finally, calculate the magnitude.

magnitude <- function(lst) {
    if (!is.list(lst)) return(lst)
    Reduce(\(a, b) 3 * magnitude(a) + 2 * magnitude(b), lst)
}

Pairs

There might be a clever way of finding the biggest pair but I use brute force:

combos <- t(combn(length(input), 2))
combos <- rbind(combos, combos[, 2:1])
magnitudes <- apply(combos, 1, \(i) magnitude(sf_reduce(input[i])))
max(magnitudes)

Transformers

This was the hardest problem yet. There were a few ideas I had about how to solve the problem quickly, possibly without even having to find any of the actual coordinates. However, I could not work out how to implement them in the end and resorted to the ‘slow’ way: iterating through the scanners one by one and trying the rotations and reflections in turn until they fit into place.

One key insight is that pairwise Euclidean distances within a scanner’s beacons are invariant to translations and rotations. Thus, if we compute within-scanner distance matrix for each one, we know that adjacent scanners must have at least \(12 * (12 - 1) / 2 = 66\) distance pairs in common. This is a good way to set up the initial search strategy.

David Robinson offers a neat solution which finds the rotations and reflections mathematically, rather than by trial and error.

The clever solution, which I suspect I had been looking for but never quite reached, is to think like a statistician should, and find the transformations using linear regression. Thanks to Ryan McCorvie.

Whatever method you go for, read in the data, with judicious handling of NAs to divide up the scanners.

read_scanners <- function(file) {
  input <- read.csv(file, header = FALSE, col.names = c('x', 'y', 'z'))
  input$scanner <- cumsum(is.na(input$z))
  type.convert(na.omit(input), as.is = TRUE)
}

Here is an initial solution, partly inspired by Ashley Baldry. Generate copies of every scanner’s beacons in the various orientations. There are supposed to be 24 different transformations but I ended up with 48.

orientate_scanners <- function(input) {
  mirror <- expand.grid(x = c(1, -1), y = c(1, -1), z = c(1, -1))
  rotate <- expand.grid(x = 1:3, y = 1:3, z = 1:3)
  rotate <- rotate[apply(rotate, 1, \(x) length(unique(x)) == 3), ]
  lapply(split(input, input$scanner), function(s) {
    s <- as.matrix(s)[, 1:3]
    apply(rotate, 1, \(r) s[, r], simplify = F) |>
      lapply(\(m) apply(mirror, 1, \(r) sweep(m, 2, r, '*'), simplify = F)) |>
      unlist(recursive = F, use.names = F)
  })
}

Then loop through the scanners and keep rotating them until the between-scanner distance matrix contains ≥12 common distances, implying that the beacons are in the right orientation and just need translating (by this distance).

locate_scanners <- function(input) {
  scanners <- orientate_scanners(input)
  scanner_coords <- rep(0, 3)
  beacons <- scanners[[1]][[1]]
  scanners <- scanners[-1]
  while (length(scanners)) {
    for (scanner in scanners) {
      done <- FALSE
      for (orientation in scanner) {
        dtbl <- table(dist2(orientation, beacons))
        if (any(dtbl >= 12)) {
          done <- TRUE
          break
        }
      }
      if (!done) next
      scanners <- setdiff(scanners, list(scanner))
      trans <- names(dtbl)[dtbl >= 12]
      trans <- scan(textConnection(trans), quiet = TRUE)
      scanner_coords <- rbind(scanner_coords, trans)
      beacons <- unique(rbind(beacons, sweep(orientation, 2, trans, '+')))
    }
  }
  list(beacons = beacons, scanners = scanner_coords)
}

The helper function dist2 stores the multi-dimensional pairwise distances as a string, and the most common one is then parsed back into a 3-vector using scan.

dist2 <- function(X, Y) {
  apply(X, 1, \(x) apply(Y, 1, \(y) paste(y - x, collapse = ' ')))
}

Linear regression

It’s not particularly satisfactory having a solution that takes over 10 seconds to run for part one, and systematically checking every rotation and reflection seems rather too much like brute force.

The next day, I refactored Ryan McCorvie’s linear regression solution as follows.

read_scanners <- function(file) {
  input <- read.csv(file, header = FALSE, col.names = c('x', 'y', 'z'))
  split(input, cumsum(is.na(input$z))) |>
    lapply('[', -1, ) |>
    lapply(type.convert, as.is = TRUE)
}
scanners <- read_scanners(input_file(19))

# Compute pairwise distances among beacons *within* each scanner.
distances <- lapply(scanners, dist)

# How many such distances does each pair of scanners have in common?
similarity <- combn(distances, 2, \(x) Reduce(intersect, x), simplify = F)

# Using this information, which scanners appear to overlap?
overlaps <- lengths(similarity) >= 12 * (12 - 1) / 2

# Just keep those pairs.
similarity <- similarity[overlaps > 0]

# But it helps if we can remember which they are.
pairs <- t(combn(length(scanners), 2)[, overlaps > 0])

# Match up those beacons that two scanners have in common.
match_beacons <- function(s1, s2) {
  m1 <- as.matrix(distances[[s1]])
  m2 <- as.matrix(distances[[s2]])
  matching <- c()
  for (j1 in 1:ncol(m1)) {
    done <- FALSE
    for (j2 in 1:ncol(m2)) {
      if (length(intersect(m1[, j1], m2[, j2])) >= 12) {
        matching <- rbind(matching, cbind(s1 = j1, s2 = j2))
        done <- TRUE
        break
      }
    }
    if (done) next
  }
  matching
}

# For two scanners, how can we orientate the 2nd to match the 1st?
align_scanners <- function(s1, s2) {
  match <- match_beacons(s1, s2)

  # Time for a bit of linear regression!
  LHS <- scanners[[s1]][match[, 's1'], , drop = F]
  RHS <- scanners[[s2]][match[, 's2'], , drop = F]
  lm_x <- lm(LHS$x ~ x + y + z, RHS)
  lm_y <- lm(LHS$y ~ x + y + z, RHS)
  lm_z <- lm(LHS$z ~ x + y + z, RHS)
  
  # The predictions gives the beacon coordinates.
  beacon_coords <- as.data.frame(apply(cbind(
    x = predict(lm_x, scanners[[s2]]),
    y = predict(lm_y, scanners[[s2]]),
    z = predict(lm_z, scanners[[s2]])), 2, round))

  # The model intercept is the translation / the scanner position!
  scanner_coords <- apply(cbind(x = coef(lm_x)[1],
                                y = coef(lm_y)[1],
                                z = coef(lm_z)[1]), 2, round)

  # Globally update
  scanners[[s2]] <<- beacon_coords
  beacons <<- unique(rbind(beacons, beacon_coords))
  scanner_locations <<- rbind(scanner_locations, scanner_coords)
}

# Do the search.
found <- 1
to_find <- 2:length(scanners)
scanner_locations <- c(0, 0, 0)
beacons <- scanners[[1]]
while (length(to_find)) {
  overlaps_with_found <- pairs[xor(pairs[, 1] %in% found,
                                   pairs[, 2] %in% found), , drop = F]
  for (i in 1:nrow(overlaps_with_found)) {
    # Let s1 be the already-aligned scanner.
    s1 <- intersect(overlaps_with_found[i, ], found)
    if (length(s1) > 1) next
    s2 <- setdiff(overlaps_with_found[i, ], s1)
    message('Joining ', s1, ' and ', s2)
    align_scanners(s1, s2)
    found <- c(found, s2)
    to_find <- setdiff(to_find, found)
  }
}

As well as being elegant, it’s also far quicker than my initial approach. You can see Ryan’s original code on GitHub.

Manhattan distance

Part two, assuming you’ve been keeping track of the scanner locations or translations, is a one-liner. Distance matrices are a built-in feature of R!

max(dist(result$scanners, 'manhattan'))

Deterministic dice

For part one, we simply have to play the game. Just write a while or repeat loop to keep going until someone achieves a score of 1000. The following code is more verbose than it needs to be:

deterministic_dice <- function(player1, player2) {
  score1 <- score2 <- nrolls <- 0
  dice <- 1:3
  repeat {
    player1 <- (player1 + sum(dice) - 1) %% 10 + 1
    score1 <- score1 + player1
    nrolls <- nrolls + 3
    if (score1 >= 1000) break

    player2 <- (player2 + sum(dice + 3) - 1) %% 10 + 1
    score2 <- score2 + player2
    nrolls <- nrolls + 3
    if (score2 >= 1000) break

    dice <- (dice + 6 - 1) %% 100 + 1
  }
  prod(nrolls, min(score1, score2))
}

Quantum states

Part two calls for recursion, and with some optimization if you want your code to run in any reasonable time frame. This means that at every opportunity, if there is room to ignore or delete universes that are redundant, we should do so.

Each player is still rolling the die three times, but there are only 7 possible unique sums of the values 1, 2 and 3:

rolls <- table(rowSums(expand.grid(1:3, 1:3, 1:3)))
rolls <- mapply(c, sum = as.numeric(names(rolls)), n = rolls, SIMPLIFY = F)

Then, at each set of rolls, we recursively count up the number of wins associated with each player for the sum of those rolls. We only need one function to do this: for player 2’s wins, we can just swap the players and accumulated scores around and pretend that player 2 is now player 1.

dirac_dice <- function(player1, player2) {
  count_wins <- memoise::memoise(
    function(player1, player2, score1, score2) {
      Reduce(
        \(wins, roll) {
          player1 <- (player1 + roll['sum'] - 1) %% 10 + 1
          score1 <- score1 + player1
          if (score1 >= 21) {
            return(wins + c(roll['n'], 0))
          } else
            wins + roll['n'] *
            rev(count_wins(player2, player1, score2, score1))
        },
        rolls,
        init = c(w1 = 0, w2 = 0))
    })
  max(count_wins(player1, player2, 0, 0))
}

To avoid revisiting universes, some caching is useful, so we can look up the value of any function call that has already been made. You could store and look up these results in a named vector, or you could do it the lazy way like I did and use the package memoise, which simply does this for you behind the scenes.

My code takes about two minutes to run for part two. Not quick, but R isn’t very fast at recursive function calls.