River Crossing Solution

Benjamin Soltoff
January 22, 2016

require(dplyr)

## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require(readr)

## Loading required package: readr

require(magrittr)

## Loading required package: magrittr

require(igraph)

## Loading required package: igraph
## 
## Attaching package: 'igraph'
## 
## The following object is masked from 'package:magrittr':
## 
##     %>%
## 
## The following objects are masked from 'package:dplyr':
## 
##     %>%, as_data_frame, groups, union
## 
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## 
## The following object is masked from 'package:base':
## 
##     union

rm(list = ls())

# number of simulations
n <- 100000

# import matrix with node structure
nodes <- read_csv("river_cross.csv") %>%
  # replace missing values with 0 - no connection
  mutate_each(funs(ifelse(is.na(.), 0, .)))

# prepare data for graph functions - set NA to zero to indicate no direct edge
mat_nodes <- function(nodes){
  nms <- nodes$node
  mat <- nodes[, -1] %>%
    as.matrix
  rownames(mat) <- nms
  colnames(mat) <- nms
  
  return(mat)
}

mat <- mat_nodes(nodes)
mat

##         start island1 island2 island3 island4 island5 island6 end
## start       0       1       1       1       0       0       0   0
## island1     1       0       1       0       1       0       0   0
## island2     1       1       0       1       0       1       0   0
## island3     1       0       1       0       0       0       1   0
## island4     0       1       0       0       0       1       0   1
## island5     0       0       1       0       1       0       1   1
## island6     0       0       0       1       0       1       0   1
## end         0       0       0       0       1       1       1   0

# create graph from adjacency matrix
g <- graph.adjacency(mat, weighted=TRUE)

# Get all path distances
s.paths <- shortest.paths(g, algorithm = "dijkstra")
s.paths <- get.shortest.paths(g, from = 1, to = 8)
s.paths

## $vpath
## $vpath[[1]]
## + 4/8 vertices, named:
## [1] start   island1 island4 end    
## 
## 
## $epath
## NULL
## 
## $predecessors
## NULL
## 
## $inbound_edges
## NULL

# Now what happens if each bridge has 50% chance of breaking?
## serial process
# system.time({
#   trials <- replicate(n, nodes %>%
#                         mutate_each(funs(ifelse(. == 1, rbinom(1, 1, .5), .)), -start, -end) %>%
#                         mat_nodes %>%
#                         graph.adjacency(., weighted = TRUE) %>%
#                         shortest.paths(., algorithm = "dijkstra") %>%
#                         .[1, nrow(.)] %>%
#                         is.finite(.))
# })
# system.time({
#   trials <- replicate(n, nodes %>%
#                         mutate_each(funs(ifelse(. == 1, rbinom(1, 1, .5), .)), -start, -end) %>%
#                         mat_nodes %>%
#                         graph.adjacency(., weighted = TRUE) %>%
#                         get.shortest.paths(., from = 1, to = 8) %>%
#                         .$vpath %>%
#                         lapply(length) %>%
#                         unlist)
# })

## parallel process
system.time({
  require(parallel)
  
  trials <- mclapply(1:n, function(x) nodes %>%
                       # if a bridge exists, 50% chance it breaks
                       mutate_each(funs(ifelse(. == 1, rbinom(1, 1, .5), .)), -start, -end) %>%
                       # convert to matrix
                       mat_nodes %>%
                       # convert to igraph
                       graph.adjacency(., weighted = TRUE) %>%
                       # calculate the shortest distance between nodes
                       shortest.paths(., algorithm = "dijkstra") %>%
                       # get the distance between start and end islands
                       .[1, nrow(.)] %>%
                       # if not finite (Inf), then no path exists
                       is.finite(.),
                     mc.cores = 6) %>%
    unlist
})

## Loading required package: parallel

##    user  system elapsed 
## 712.190   6.318 150.406

# proportion of trials where traveller can still get from start to end island
sum(trials) / length(trials)

## [1] 0.98368