Random Chutes and Ladders Board Generator (June 2016 Community Challenge)

Question

This is my submission for the June 2016 Community Challenge. It is R code that takes as input the number of squares on the board, the desired number of chutes and ladders, and the desired total sum of the chute and ladder deltas. It then samples uniformly at random from all feasible boards with these properties to return a board; each feasible board fitting the input requirements is given exactly the same probability of being selected.

To do this, the code defines a mathematical program that has a binary decision variable for each valid chute or ladder. The constraints are used to ensure a feasible board is returned (no more than one chute/ladder starting at a square, a chute cannot be chained to a ladder, and a ladder cannot be chained to a chute) and that the input requirements are met (desired number of chutes and ladders and desired sum of deltas). Each chute and ladder is assigned a random weight and the optimization model seeks the feasible board with minimum weight; this ensures we are selecting uniformly at random from all feasible boards.

library(lpSolve)

random.game.board <- function(squares=100, num.chute=9, num.ladder=9, delta=-50) {
  if (squares < 3) stop("Game board must have at least 3 squares")
  all.paired.squares <- expand.grid(from=seq_len(squares), to=seq_len(squares))
  chutes <- subset(all.paired.squares, from < squares & to < from)
  ladders <- subset(all.paired.squares, from > 1 & to > from)
  mod <- lp(objective.in = runif(nrow(chutes) + nrow(ladders)),
            const.mat =
              rbind(t(sapply(seq_len(squares), function(x) {
                      # Square x begins no more than one chute or ladder
                      as.numeric(c(chutes$from == x, ladders$from == x))
                    })),
                    t(sapply(seq_len(squares), function(x) {
                      # Square x does not begin a chute and end a ladder
                      as.numeric(c(chutes$from == x, ladders$to == x))
                    })),
                    t(sapply(seq_len(squares), function(x) {
                      # Square x does not end a chute and begin a ladder
                      as.numeric(c(chutes$to == x, ladders$from == x))
                    })),
                    # Number of chutes
                    rep(1:0, c(nrow(chutes), nrow(ladders))),
                    # Number of ladders
                    rep(0:1, c(nrow(chutes), nrow(ladders))),
                    # Delta of chutes and ladders
                    c(chutes$to - chutes$from, ladders$to - ladders$from)),
            const.dir=rep(c("<=", "="), c(3*squares, 3)),
            const.rhs=c(rep(1, 3*squares), num.chute, num.ladder, delta),
            all.bin=TRUE)
  if (mod$status != 0) stop("No feasible game boards")}
  board <- rbind(chutes, ladders)[mod$solution > 0.999,]
  board$delta <- board$to - board$from
  board <- board[order(board$from),]
  row.names(board) <- NULL
  board
}

This is an optimization problem with a huge number of binary variables (9702), but it is so loosely constrained that the open-source lpSolve package can solve it to optimality in about 1 second on my computer with 100 squares and 9 chutes and ladders, yielding a randomly selected board:

set.seed(144)
random.game.board()
#    from  to delta
# 1     3  75    72
# 2    10  43    33
# 3    12  38    26
# 4    18  41    23
# 5    26   4   -22
# 6    34  13   -21
# 7    50  69    19
# 8    53  25   -28
# 9    56  71    15
# 10   60  59    -1
# 11   64  22   -42
# 12   65  29   -36
# 13   66  85    19
# 14   72  84    12
# 15   74  23   -51
# 16   82  44   -38
# 17   89  57   -32
# 18   98 100     2

I would appreciate comments on any aspects of my code, though I would especially be interested in comments on:

1. Whether more efficient approaches exist to sample uniformly across all feasible boards
2. Approaches to vectorize the construction of the constraint matrix
3. The error handling
4. Whether I could condense the three lines at the end of the function that order the rows, restore default row names, and return.

Answer 1

After looking through this code with fresh eyes, I identified a way to improve the construction of the constraint matrix; following the guidance from the meta site, I'm posting this as a self-answer instead of editing the question.

The first class of constraints in my optimization model have a row for each square of the board and a column for each variable (which corresponds to a chute or ladder). The value in row i and column j is 1 if chute/ladder j starts on square i and 0 otherwise. I originally achieved this with the following code:

t(sapply(seq_len(squares), function(x) {
  # Square x begins no more than one chute or ladder
  as.numeric(c(chutes$from == x, ladders$from == x))
}))

Basically, I looped through each square using sapply and obtained a 1/0 vector for whether each chute/ladder began on that square. Since sapply arranges the returned results by column, I transposed the result with t.

This is actually a task much better suited for the outer function:

outer(seq_len(squares), c(chutes$from, ladders$from), "==")

This function associates each square number (seq_len(squares)) with each row and the start square for each variable (c(chutes$from, ladders$from)) with each column, and it returns a matrix of whether the row and column are equal. In addition to being shorter and easier to read, this has an additional property of being fully vectorized, which should make it more efficient in situations where we have a lot of squares.

The same approach can be applied for the two subsequent constraint blocks in the constraint matrix. These three constraint blocks can be written as:

# Each square begins no more than one chute or ladder
outer(seq_len(squares), c(chutes$from, ladders$from), "=="),
# Each square does not begin a chute and end a ladder
outer(seq_len(squares), c(chutes$from, ladders$to), "=="),
# Each square does not end a chute and begin a ladder
outer(seq_len(squares), c(chutes$to, ladders$from), "=="),