# 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