Map Coloring with OMPR and HiGHS

Author

notesofdabbler

Published

September 10, 2022

Vertex coloring problem in a graph is to find the minimum number of colors such that no two connected vertices have the same color. A special case is map coloring problem where the vertices are states/countries in a map and edges are pairs of states/countries that are adjacent to each other. Here we illustrate the map coloring problem with OMPR modeling language and HiGHS solver. Specifically, we want to find the minimum number of colors and the color of each state in US so that no two adjacent states have the same color. Only portions of R code are shown here. The .qmd file is in this location

Code
# load libraries and source functions
library(dplyr)
library(glue)
library(ompr)
library(ompr.highs)
library(highs)
library(purrr)
library(ggplot2)
library(maps)
library(mapdata)

theme_set(theme_light())

State Adjacency Data

We use the state adjacency data from Gregg Lind’s site and process it in a format that we can use for OMPR model.

Code
#
# US states adjacency data from
# https://writeonly.wordpress.com/2009/03/20/adjacency-list-of-states-of-the-united-states-us/
#

x = readLines("data/US_state_adjacency.txt")

edgelist = list()
e = 0
for (i in 1:length(x)) {
  stloc = strsplit(x[i], ",")[[1]]
  ns = length(stloc)
  if(ns > 1) {
    for (i in 2:ns) {
      e = e + 1
      edgelist[[e]] = c(stloc[1], stloc[i])   
    }
  }
}

edge_df = as.data.frame(do.call(rbind, edgelist))
names(edge_df) = c("from", "to")

nodes_df = bind_rows(edge_df %>% distinct(from) %>% rename(state = from),
          edge_df %>% distinct(to) %>% rename(state = to)) %>% distinct(state)
nodes_df = nodes_df %>% mutate(id = seq(1, nrow(nodes_df)))

edge_df = inner_join(edge_df, nodes_df %>% rename(from = state), by = "from") %>% rename(fromid = id)
edge_df = inner_join(edge_df, nodes_df %>% rename(to = state), by = "to") %>% rename(toid = id)

edge_df = edge_df %>% filter(fromid < toid)

There are many mathematical formulations of this problem. We use a basic formulation for the example here.

\[ \begin{array}{llr} \min & \sum_{c=1}^Cy_c & (a)\\ & \sum_{c=1}^Cx_{ic} = 1, \;i=1,2,\ldots,N & (b)\\ & x_{ic} + x_{jc} \leq y_c, \; \mbox{when }i, j \mbox{ are adjacent} & (c)\\ & x_{ic} \; binary \\ & y_c \; binary \end{array} \]

where:

  • \(y_c=1\) if color \(c\) is used, \(x_{ic}=1\) if state \(i\) is colored with color \(c\).
  • Objective (a) is to minimize the number of colors used
  • Constraint (b) ensures that each state gets some color
  • Constraint (c) ensures that if state \(i\) and \(j\) are adjacent, they don’t get the same color.
# OMPR model

ns = nrow(nodes_df)
nc = 4
edge_str = edge_df %>% mutate(edge_str = glue("{fromid}_{toid}")) %>% pull(edge_str)

mdl = MIPModel()
mdl = mdl %>% add_variable(x[i, c], i = 1:ns, c = 1:nc, type = "integer", lb = 0, ub = 1)
mdl = mdl %>% add_variable(y[c], c = 1:nc, type = "integer", lb = 0, ub = 1)
mdl = mdl %>% set_objective(sum_over(y[c], c=1:nc), sense = "min")
mdl = mdl %>% add_constraint(sum_over(x[i, c], c = 1:nc) == 1, i = 1:ns)
mdl = mdl %>% add_constraint(x[i, c] + x[j, c] <= y[c], i = 1:ns, j = 1:ns, c = 1:nc, glue("{i}_{j}") %in% edge_str)

Solving the model shows that 4 colors are sufficient to color the US map (Hawaii and Alaska are excluded).

Code
# solve model
sol = mdl %>% solve_model(highs_optimizer())


# Get solution and visualize in a map
sol[["status"]]
[1] "optimal"
Code
sol[["objective_value"]]
[1] 4
Code
soldf = sol %>% get_solution(x[i, c]) %>% filter(value == 1) %>%
  rename(id = i, colid = c)

soldf = inner_join(soldf, nodes_df, by = "id")

colmap = c("1" = "red", "2" = "blue", "3" = "green", "4" = "yellow")
soldf = soldf %>% mutate(stcol = colmap[colid])

The solution is visualized in a map using code from this location

Code
#
# https://jtr13.github.io/cc19/different-ways-of-plotting-u-s-map-in-r.html
#

state <- map_data("state")
state_name_abb = tibble(stname = tolower(state.name), stabb = state.abb)

state = inner_join(state, state_name_abb %>% rename(region = stname), by = "region")
state = inner_join(state, soldf %>% select(state, stcol) %>% rename(stabb = state), by = "stabb")

ggplot(data=state, aes(x=long, y=lat, fill=stcol, group=group)) + 
  geom_polygon(color = "white") + 
  guides(fill="none") + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) + 
  ggtitle('U.S. Map with States') + 
  coord_fixed(1.3)

Code
sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] mapdata_2.3.0         maps_3.4.0            ggplot2_3.3.6        
[4] purrr_0.3.4           highs_0.1-2           ompr.highs_0.0.1.9000
[7] ompr_1.0.2            glue_1.6.2            dplyr_1.0.10         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9        pillar_1.8.1      compiler_3.6.3    tools_3.6.3      
 [5] digest_0.6.29     gtable_0.3.1      jsonlite_1.8.0    evaluate_0.16    
 [9] lifecycle_1.0.1   tibble_3.1.8      checkmate_2.1.0   lattice_0.20-40  
[13] pkgconfig_2.0.3   rlang_1.0.5       Matrix_1.2-18     cli_3.3.0        
[17] yaml_2.3.5        xfun_0.32         fastmap_1.1.0     listcomp_0.4.1   
[21] withr_2.5.0       stringr_1.4.1     knitr_1.40        generics_0.1.3   
[25] vctrs_0.4.1       htmlwidgets_1.5.4 grid_3.6.3        tidyselect_1.1.2 
[29] data.table_1.14.2 R6_2.5.1          fansi_1.0.3       rmarkdown_2.16   
[33] farver_2.1.1      magrittr_2.0.3    ellipsis_0.3.2    scales_1.2.1     
[37] backports_1.4.1   htmltools_0.5.3   colorspace_2.0-3  labeling_0.4.2   
[41] utf8_1.2.2        stringi_1.7.8     munsell_0.5.0     lazyeval_0.2.2