Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

### Introduction

Advent of code 2018, day 14, was a tough one to do in R. Here is a link. I did this problem a lot of different ways, trying to find a nice way to do it. In the end, I did it using vectors in base R, data.table and finally using the standard template library in Rcpp. I had never used Rcpp really before, or even programmed in C++, so this was a bit of a challenge! I do know C pretty well, but not great, so that was some help. I’ll try to document some of the learning steps that I needed to go through.

### First Star

Here’s a very brief recap of the puzzle. We are building a vector that can get arbitrarily long by adding the values at two positions of the vector (and doing some manipulation), then we change the two positions by moving them to the right by between 1 and 10 spaces. If the new position falls off the end of the vector, then we wrap around to the beginning.

My puzzle input was 380621, so the first job was to find the value of the ten recipes that follow the 380,621st recipe. OK, that doesn’t seem to hard to do in R.

This was essentially my first attempt.

make_recipe <- function(sc, e1, e2) {
s <- sc[e1] + sc[e2]
if(floor(s/10) == 0)
return(s)
return(c(floor(s/10), s%%10))
}

my_mod <- function(x, m) {
ifelse(x%%m == 0, m, x%%m)
}

sc <- c(3, 7) #first two recipes have values of 3 and 7.
e1 <- 1       #this is the location of elf 1's recipe
e2 <- 2       #this is the location of elf 2's recipe

sc <- c(sc, make_recipe(sc, e1, e2))   #add the value of the new recipe to the end of the vector.
a <- length(sc)                        #find number of recipes so we can move the elves
e1 <- my_mod(e1 + sc[e1] + 1, a)       #move elf 1
e2 <- my_mod(e2 + sc[e2] + 1, a)       #move elf 2

system.time(
while(a < 400000) {
sc <- c(sc, make_recipe(sc, e1, e2))
a <- length(sc)
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)
#  if(a %% 10000 == 0)
#    print(a)  #Just so we can monitor progress!
} #This takes a pretty long time on my computer, but it does eventually finish.
)
sc[380622:380631] #First star

That wasn’t too, too bad, but let’s see if we can’t speed it up. The elapsed time was 365.496, which seems excessive! I think the big slow-down was that I was building the vector one or two elements at a time. In the version below, I decided to build up a bunch of recipes until I needed to wrap the elf’s position back around to the start of the recipe vector. Then, I would mostly be adding vectors to smallish vectors, and occasionally combining two longer vectors. This should be faster.

make_recipe <- function(sc, e1, e2) {
s <- sc[e1] + sc[e2]
if(floor(s/10) == 0)
return(s)
return(c(floor(s/10), s%%10))
}

my_mod <- function(x, m) {
ifelse(x%%m == 0, m, x%%m)
}

sc <- c(3, 7) #first two recipes have values of 3 and 7.
e1 <- 1       #this is the location of elf 1's recipe
e2 <- 2       #this is the location of elf 2's recipe

sc <- c(sc, make_recipe(sc, e1, e2))   #add the value of the new recipe to the end of the vector.
a <- length(sc)                        #find number of recipes so we can move the elves
e1 <- my_mod(e1 + sc[e1] + 1, a)       #move elf 1
e2 <- my_mod(e2 + sc[e2] + 1, a)       #move elf 2

builder <- integer(0)
system.time(
while(a < 400000) {
if(e1 < length(sc) - 20 && e2 < length(sc) - 20) {  # I am being careful here. I don't want the elves to wrap without
builder <- c(builder, make_recipe(sc, e1, e2))    # knowing about it!
} else {
if(length(builder) > 0) {                         # Here, we need to either just add the new recipes to the end,
sc <- c(sc, builder)                            # or add all of the built recipes.
builder <- integer(0)
}
sc <- c(sc, make_recipe(sc, e1, e2))
}
a <- length(sc)
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)
# if(length(builder) == 0)                             # again, just measuring progress.
#    print(a)
}
)
##    user  system elapsed
##  20.647   9.976  32.488
sc[380622:380631]
##   6 9 8 5 1 0 3 1 2 2

Woo-hoo! That was much better. The elapsed time was 29 seconds, a better than 10x speed-up. There are several things I might try to speed this up further, and those are going to be looked at once I get to the second part of the challenge.

### Second Star

For the second star, I need to keep going until I see the pattern 380621. That could take a while. There are some really cool theorems about waiting for patterns like this; I think they follow from Blackwell’s Theorem. For example, the expected time before observing HTHHTHTT is less than the expected time before observing TTTTTTTT. Anyway, I noticed that this builder idea wasn’t working well with longer vectors, because the builder itself was getting way too big. So, I just made a maximum builder size. Let’s compare the speed of that versus the speed of the previous algorithm in the first star case.

sc <- c(3, 7) #first two recipes have values of 3 and 7.
e1 <- 1       #this is the location of elf 1's recipe
e2 <- 2       #this is the location of elf 2's recipe

sc <- c(sc, make_recipe(sc, e1, e2))   #add the value of the new recipe to the end of the vector.
a <- length(sc)                        #find number of recipes so we can move the elves
e1 <- my_mod(e1 + sc[e1] + 1, a)       #move elf 1
e2 <- my_mod(e2 + sc[e2] + 1, a)       #move elf 2

builder <- integer(0)
system.time(
while(a < 400000) {
if(e1 < length(sc) - 20 && e2 < length(sc) - 20 && length(builder) < 1000) {  #ONLY CHANGE: Max builder size is 1000
builder <- c(builder, make_recipe(sc, e1, e2))
} else {
if(length(builder) > 0) {                         # Here, we need to either just add the new recipes to the end,
sc <- c(sc, builder)                            # or add all of the built recipes.
builder <- integer(0)
}
sc <- c(sc, make_recipe(sc, e1, e2))
}
a <- length(sc)
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)
#if(length(builder) == 0)                             # again, just measuring progress.
#  print(a)
}
)
##    user  system elapsed
##   4.950   0.904   5.905
sc[380622:380631]
##   6 9 8 5 1 0 3 1 2 2

That’s about 6 seconds. Things are looking up! Note that is a 60x speed up over the original algorithm! I happen to know that I am going to need to build a vector of length about 40 million for the second star, which would mean that we would need 600 seconds, even if there is no slowdown as the vectors get longer (which I seriously doubt).

library(data.table)
sc <- rep(0L, 500000) #assume max size is 500,000
sc <- data.table(sc)

make_recipe_data_table <- function(sce1, sce2) {
s <- sce1 + sce2  #sc$sc[e1] + sc$sc[e2]
if(floor(s/10) == 0)
return(s)
return(c(floor(s/10), s%%10))
} #thought this might be faster. not sure why, now?

my_mod_data_table <- function(x, m) {
ifelse(x%%m == 0, m, x%%m)
}

cur_pos <- 1L
set(sc, i = c(1L, 2L), j = 1L, value = c(3L, 7L))  #note that to keep things faster, we have to specify that everything is an integer. Don't want to force data.table to cast, as that slowed things way down.
e1 <- 1L
e2 <- 2L
cur_pos <- cur_pos + 2L

sce1 <- sc$sc[e1] sce2 <- sc$sc[e2]

rows <- cur_pos
} else{
rows <- c(cur_pos, cur_pos + 1L)
}

set(sc, i = rows, j = 1L, value = add_vec)
a <- cur_pos - 1L
e1 <- as.integer(my_mod_data_table(e1 + sc$sc[e1] + 1, a)) e2 <- as.integer(my_mod_data_table(e2 + sc$sc[e2] + 1, a))

continue <- TRUE
system.time(
while(a < 400000) {
sces <- sc$sc[c(e1, e2)] add_vec <- as.integer(make_recipe_data_table(sces, sces)) if(length(add_vec) == 1L) { rows <- cur_pos } else{ rows <- c(cur_pos, cur_pos + 1L) } set(sc, i = rows, j = 1L, value = add_vec) cur_pos <- cur_pos + length(add_vec) a <- cur_pos - 1L e1 <- my_mod_data_table(e1 + sces + 1L, a) e2 <- my_mod_data_table(e2 + sces + 1L, a) #if(i%%1000000L == 0) { # print(a) #ss <- str_flatten(as.character(sc$sc))
#continue <- !str_detect(ss, "380621")
#print(continue)
#  }
}
)
##    user  system elapsed
##   9.490   0.057   9.615

I tried builder lengths of 500 and 2000 as well, and the 1000 length seemed to be the fastest.

sc <- c(3, 7) #first two recipes have values of 3 and 7.
e1 <- 1       #this is the location of elf 1's recipe
e2 <- 2       #this is the location of elf 2's recipe

sc <- c(sc, make_recipe(sc, e1, e2))   #add the value of the new recipe to the end of the vector.
a <- length(sc)                        #find number of recipes so we can move the elves
e1 <- my_mod(e1 + sc[e1] + 1, a)       #move elf 1
e2 <- my_mod(e2 + sc[e2] + 1, a)       #move elf 2

builder <- integer(0)
system.time(
while(a < 400000) {
if(e1 < length(sc) - 20 && e2 < length(sc) - 20 && length(builder) < 1000) {  #ONLY CHANGE: Max builder size is 1000
builder <- c(builder, make_recipe(sc, e1, e2))
} else {
if(length(builder) > 0) {                         # Here, we need to either just add the new recipes to the end,
sc <- c(sc, builder)                            # or add all of the built recipes.
builder <- integer(0)
}
sc <- c(sc, make_recipe(sc, e1, e2))
}
a <- length(sc)
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)
#if(length(builder) == 0)                             # again, just measuring progress.
#  print(a)
}
)
##    user  system elapsed
##   5.190   0.875   6.145

For the second star, we had to look until we got the pattern 380621. I had a hard time with this one, basically because I wasn’t patient enough, and I thought my code must have a bug.

My idea was that I would pretty much do what I did above, but I would have a maximum builder size so that part doesn’t start slowing down too much. I feel like I could have had a builder for my builder, but I didn’t want to go crazy. But wait, first, I wanted to know what parts of my previous algorithm were slowing me down the most.

make_recipe_2 <- function(sce1, sce2) {
s <- sce1 + sce2  #sc$sc[e1] + sc$sc[e2]
if(floor(s/10) == 0)
return(s)
return(c(floor(s/10), s%%10))
}

sc <- rep(1, 4000000)
microbenchmark::microbenchmark(
make_recipe(sc, 1000, 200000),
c(sc, 1),
{sce1 <- sc; sce2 <- sc; make_recipe_2(sce1, sce2)},
{d <- sample(4000000, 2); sc[d]},
{d <- sample(4000000, 2); sc[d]; sc[d]}
)
## Unit: microseconds
##                                                                          expr
##                                                  make_recipe(sc, 1000, 2e+05)
##                                                                      c(sc, 1)
##  {     sce1 <- sc     sce2 <- sc[2e+05]     make_recipe_2(sce1, sce2) }
##                                       {     d <- sample(4e+06, 2)     sc[d] }
##                       {     d <- sample(4e+06, 2)     sc[d]     sc[d] }
##        min         lq        mean     median         uq       max neval
##      1.034     3.9760    11.31215    12.1425    17.0285    39.012   100
##  10014.234 13465.5705 17120.52131 17175.6540 18713.6340 66800.367   100
##      1.644     2.8445   110.69945    13.3150    15.6845  9928.314   100
##   2314.042  2600.5335  4998.77946  2997.7615  8734.1630 13446.779   100
##   2318.337  2608.5985  5566.62225  2905.7695  8793.2450 61288.794   100
##  cld
##  a
##    c
##  a
##   b
##   b

I also wondered whether it would be faster to not pass sc to make_recipe, but oddly it doesn’t seem like it is. At any rate, adding those values to the end of the long vectors is for sure what is slowing things down. Let’s see how I do:

sc <- c(3, 7)
e1 <- 1
e2 <- 2

a <- length(sc)
sc <- c(sc, make_recipe(sc, e1, e2))
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)

old_sc <- sc
builder <- integer(0)
continue <- TRUE          #we need to run this loop until we see the pattern 380621. Ha.
i <- 1
build_max <- 1000         #at the beginning, we will build up vectors of length 1000 to add to sc. As sc gets longer, the build also gets longer.
#I should have optimized this more, but I played around with it and this seems OK? Note that it is *WAY* faster to get
#to 400,000. Which is good, because otherwise we would have had a long, long wait.
while(continue) {
i <- i + 1
if(e1 < a - 11 && e2 < a - 11 && length(builder) < build_max) {
builder <- c(builder, make_recipe(sc, e1, e2))
} else {
if(length(builder) > 0) {
sc <- c(sc, builder)
builder <- integer(0)
}
sc <- c(sc, make_recipe(sc, e1, e2))
}
a <- length(sc)
e1 <- my_mod(e1 + sc[e1] + 1, a)
e2 <- my_mod(e2 + sc[e2] + 1, a)
if(i%%20000 == 0) {                                  #The str_flatten + str_detect takes a really long time when the vector gets long. So, I kept
print(a)                                           #track of what hadn't yet been checked and just checked them. We do this every 20K iterates, but
sc_cont <- sc                                      #that number was picked out of my ass.
e1_cont <- e1
e2_cont <- e2
a_cont <- a
builder_cont <- builder
l1 <- max(a - 40000, 1)
ss <- stringr::str_flatten(as.character(sc[l1:a]))
continue <- !stringr::str_detect(ss, "380621")
if(a > 500000)
build_max <- 1500
if(a > 1000000)
build_max <- 2000
if(a > 2000000)
build_max <- 2500
if(a > 4000000)
build_max <- 3000
if(a > 8000000)
build_max <- 3500
if(a > 16000000)
build_max <- 4000
}
}

If you are patient enough, the above code will, I think, eventually work. I have never been patient enough. I have heard people say that I should pre-allocate memory for the vector rather than building it on the fly. I tried pre-allocating sc to have 30 million zeroes, but that was way slower. Maybe if I would have only pre-allocated builder, that would have helped. I’m just not sure.

At any rate, I next decided to use data.table, because it seems that building the vector is still what is taking time. I didn’t use much from data.table. I used set, which is the data.table way of doing vec[] <-, and to access previously set values in the vector, I used $ which as I understand it, casts to a regular vector and used []. I looked into using native data.table things, but they were slower when I benchmarked. library(data.table) sc <- rep(0L, 50000000) #assume max size is 50,000,000 sc <- data.table(sc) make_recipe_data_table <- function(sce1, sce2) { s <- sce1 + sce2 #sc$sc[e1] + sc$sc[e2] if(floor(s/10) == 0) return(s) return(c(floor(s/10), s%%10)) } #thought this might be faster. not sure why, now? my_mod_data_table <- function(x, m) { ifelse(x%%m == 0, m, x%%m) } cur_pos <- 1L set(sc, i = c(1L, 2L), j = 1L, value = c(3L, 7L)) #note that to keep things faster, we have to specify that everything is an integer. Don't want to force data.table to cast, as that slowed things way down. e1 <- 1L e2 <- 2L cur_pos <- cur_pos + 2L sce1 <- sc$sc[e1]
sce2 <- sc$sc[e2] add_vec <- as.integer(make_recipe_data_table(sce1, sce2)) if(length(add_vec) == 1L) { rows <- cur_pos } else{ rows <- c(cur_pos, cur_pos + 1L) } set(sc, i = rows, j = 1L, value = add_vec) cur_pos <- cur_pos + length(add_vec) a <- cur_pos - 1L e1 <- as.integer(my_mod_data_table(e1 + sc$sc[e1] + 1, a))
e2 <- as.integer(my_mod_data_table(e2 + sc$sc[e2] + 1, a)) continue <- TRUE for(i in 1:18000000) { sces <- sc$sc[c(e1, e2)]
rows <- cur_pos
} else{
rows <- c(cur_pos, cur_pos + 1L)
}

set(sc, i = rows, j = 1L, value = add_vec)
a <- cur_pos - 1L
e1 <- my_mod_data_table(e1 + sces + 1L, a)
e2 <- my_mod_data_table(e2 + sces + 1L, a)
if(i%%1000000L == 0) {
print(a)
#ss <- str_flatten(as.character(sc\$sc))
#continue <- !str_detect(ss, "380621")
#print(continue)
}
}
##  1300566
##  2604631
##  3908479
##  5212451
##  6515907
##  7819555
##  9122859
##  10426309
##  11729730
##  13032975
##  14336480
##  15640992
##  16944212
##  18247938
##  19551908
##  20855419
##  22158850
##  23462488