Colour Maps using Integer Programming

Operations Research has a lot of cool and practical applications such as Supply Chain Optimization (my favourite!), crew scheduling, vehicle routing, and portfolio optimization.

This article discusses one of these cool Operations Research applications. We will solve a map colouring puzzle using a technique known as Integer Programming (IP). The puzzle states that given a map with bordered entities (e.g. countries, states, counties) what is the minimum number of colours required to colour each entity such that no two adjacent entities have the same colour.

This map colouring problem corresponds to a famous problem in Operations Research known as the vertex colouring problem. In the vertex colouring problem, we are given an undirected graph and asked to assign colours to its vertices such that no two connected vertices have the same colour and of course use the least number of colours possible. The correspondence between map colouring and vertex colouring is very obvious, which makes it a great example of how “business problems” map to “generic Operations Research problems” as discussed in my previous article The Big Picture of Operations Research.

The general process is depicted in the figure below.

• Choose an uncoloured map
• Model the map as an undirected graph
• Colour the undirected graph using your favourite algorithm
• Convert the coloured undirected graph back to the corresponding map

I wrote two Python packages, vertex-color-solver and map-grapher, that will help us carry out this process in about 10 lines of code. Let’s walk through the code.

First you need to install the packages by following the instructions provided on their Github repositories.

After you install the packages successfully, import the following classes and functions

from vertex_colorer.solver import Solver
from map_grapher.core import load_map_data, plot_map, color_map,\
save_to_file, build_adjmat

The functions loaded from map_grapher allow you to select a map, model it as an undirected graph, and use a coloured undirected graph to colour your map. The package currently supports all county maps of states in the United States as well as a national map.

In this demo, I will choose the Arizona as an example. Let’s load the map

state = load_map_data('Arizona')

the state variable holds map data loaded from a geojson file. Next we will plot the map using the plot_map function

figure = plot_map(state, map_title='Arizona map')

If you are using a Jupyter Notebook, you can view the map directly from the figure, if you are on the command line you would want to save the figure to a file

save_to_file (figure, 'arizona_black_and_white_map.png')

Nice! Now the next step is to model this map as an undirected graph. We can do this using the following code

matrix = build_adjmat(state)

build_adjmat function takes map data and converts it to an adjacency matrix, which is the data structure we will use to represent an undirected graph. The adjacency matrix looks like this

[[1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
 [0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1],
 [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1],
 [1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0],
 [1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1],
 [0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0],
 [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0],
 [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
 [0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0],
 [0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1],
 [1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0],
 [0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1],
 [0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1]]

Each element represents an edge between two vertices of the undirected graph. If there is an edge (i.e. the two vertices are connected, or from the business problem’s perspective, the corresponding map entities are adjacent) then the value is 1, otherwise it is 0.

We can now pass this matrix to our solver

solver = Solver(matrix)

The solver is instantiated and holds the undirected graph. To colour the graph, you need to pick an approach to use. Here I am using the term approach instead of algorithm, because they are not quite the same. For example, the IP solver we will use in this demo models the undirected graph itself as an Integer Program then uses a Branching algorithm (among other things) to solve the model, then convert it back to the undirected graph. Alternatively, we could operate on the undirected graph directly using a Dynamic Programming algorithm or a greedy heuristic.

Let’s ask the solver to colour the undirected graph using IP

solver.solve('IP')

and voila! we are done. The solution can be accessed using the solution attribute.

solver.solution

this will give us the following data structure

[[0, 0], [1, 1], [2, 3], [3, 1], [4, 2], [5, 0], [6, 2], [7, 3], [8, 3], [9, 3], [10, 3], [11, 0], [12, 2], [13, 1], [14, 2]]]

each pair in this list holds two values, the first is a vertex identifier (an index) and the second is the colour assigned.

Finally we pass this solution up to the business level problem using the color_map function

colored_map = color_map(state, solver.solution)

and we follow the same procedure to plot the map and save it to a file

fig = plot_map(colored_map, map_title='map of Arizona')
save_to_file (fig, 'arizona_colored_map.png')

Here is the entire script

#!/usr/bin/env python

from vertex_colorer.solver import Solver
from map_grapher.core import load_map_data, plot_map, color_map,\
save_to_file, build_adjmat

SELECTED_MAP = 'Arizona'

if __name__=="__main__":
    state = load_map_data(SELECTED_MAP)
    figure = plot_map(state, SELECTED_MAP+'_bw')
    matrix = build_adjmat(state)
    solver = Solver(matrix)
    solver.solve('IP')
    colored_map = color_map(state, solver.solution)
    fig = plot_map(colored_map, map_title=SELECTED_MAP+'_colored')
    save_to_file (fig, 'arizona_colored_map.png')

You should be able to use the code yourself and colour some maps. Notice how depending on the state you choose to colour, the solver runtime varies. After experimenting with few states myself, I recommend you start with ones that do not have that many counties and work your way up. You will notice that runtime is not a linear function of the number of counties (vertices) but rather how dense the connections are. You will see that when you tryto run California (Good luck!).

So far, I have only implemented the IP solver. Next stop is to implement other solvers such as Constraint Programming, Dynamic Programming and greedy heuristics and see how they compare with the IP. The full code and more detailed explanation are available on this gist. Enjoy and feel free to write to me if you have questions!