4-color-berlin.R

# based on a article from here https://dirkschumacher.github.io/ompr/articles/problem-graph-coloring.html

library(maptools)
library(dplyr)

# CC by
# license here https://github.com/berlinermorgenpost/Berlin-Geodaten
map_data <- rgdal::readOGR("https://github.com/berlinermorgenpost/Berlin-Geodaten/raw/master/berlin_bezirke.geojson", "OGRGeoJSON")

# this gives as an adjancy list
neighbors <- spdep::poly2nb(map_data)

# a helper function that determines if two nodes are adjacent
is_adjacent <- function(i, j) {
  purrr::map2_lgl(i, j, ~ .y %in% neighbors[[.x]])
}

# build an optimization model

library(ompr)

n <- nrow(map_data@data)
max_colors <- 4 # let's try 4 colors

# based on the formulation from here
# http://wwwhome.math.utwente.nl/~uetzm/do/IP-FKS.pdf
model <- MILPModel() %>%

  # 1 iff node i has color k
  add_variable(x[i, k], type = "binary", i = 1:n, k = 1:max_colors) %>%

  # 1 iff color k is used
  add_variable(y[k], type = "binary", k = 1:max_colors) %>%

  # minimize colors
  # multiply by k for symmetrie breaking (signifcant diff. in solution time)
  set_objective(sum_expr(colwise(k) * y[k], k = 1:max_colors), sense = "min") %>%

  # each node is colored
  add_constraint(sum_expr(x[i, k], k = 1:max_colors) == 1, i = 1:n) %>%

  # if a color k is used, set y[k] to 1
  add_constraint(sum_expr(x[i, k], i = 1:n) <= n * y[k], k = 1:max_colors) %>%

  # no adjacent nodes have the same color
  add_constraint(x[i, k] + x[j, k] <= 1, i = 1:n, j = 1:n, k = 1:max_colors, is_adjacent(i, j))

model

library(ROI)
library(ROI.plugin.glpk)
library(ompr.roi)

result <- solve_model(model, with_ROI("glpk", presolve = TRUE, verbose = TRUE))

result

# Last step is to plot the result. First we will get the colors from the optimal solution.

assigned_colors <- get_solution(result, x[i, k]) %>%
  filter(value > 0.9) %>%
  arrange(i)

library(ggplot2)
color_data <- map_data@data
color_data$color <- assigned_colors$k
plot_data_fort <- fortify(map_data, region = "name") %>%
  left_join(select(color_data, name, color),
            by = c("id" = "name")) %>%
  mutate(color = factor(color))

#Now we have everything to plot it:

ggplot(plot_data_fort, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = color)) +
  coord_quickmap() +
  hrbrthemes::theme_ipsum() +
  viridis::scale_fill_viridis(discrete = TRUE, option = "D") +
  ggtitle("Berlin in 4 colors",
          subtitle = "Find the minimal number of colors, such that no two adjacent districts share the same color")