4-color-world.R

# 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())