Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.