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)
# devtools::install_github("dirkschumacher/ompr@milp")
# CC by
map_data <- rgdal::readOGR("https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_50m_admin_0_countries.geojson", "OGRGeoJSON")
# this gives us 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)
library(ROI)
library(ROI.plugin.cbc)
library(ompr.roi)
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
system.time({
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))
result <- solve_model(model, with_ROI("cbc", presolve = 1))
})
model
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 = "ADMIN") %>%
left_join(select(color_data, ADMIN, color),
by = c("id" = "ADMIN")) %>%
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_map("gilbert") +
hrbrthemes::theme_ipsum() +
viridis::scale_fill_viridis(discrete = TRUE, option = "D") +
ggtitle("The world in 4 colors",
subtitle = "Find the minimal number of colors, such that no two adjacent countries share the same color") +
theme(axis.text = element_blank())
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.