Applying the 4 color theorem to list of neighbor polygons stocked in a graph array

I want to apply the 4 color theorem ie. each neighbor polygons on a map should have a different color. The theorem state that only 4 colors is needed for any kind of map.

As input I have an array of polygon containing id and color id and a graph array of adjacent polygons. As output I want to associate a color id between 0-3 (or at maximum 0-4) to each polygons ensuring that adjacent polygons have different color id.

I have the following code (from this question) that returns a different color id for each neighboor, but it does not ensure that the minimum number of colors is used.

#array of input polygon tuples (id, color_id) color_id is None at beginning
input_polygons = [(0, None), (1, None), (2, None), (3, None), (4, None), (5, None), (6, None), (7, None), (8, None)]

#graph array of neighbors ids
adjacent_polygons_graph = [[0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [1, 0], [2, 0], [2, 3], [3, 0], [3, 2], [3, 6], [3, 8], [4, 0], [4, 5], [4, 6], [5, 0], [5, 4], [6, 3], [6, 4], [6, 7], [7, 6], [8, 3]]

colored = []
for polygon in input_polygons:
    nbrs = [second for first, second in adjacent_polygons_graph if first == polygon[0]]
    used_colors = []
    for nbr in nbrs:
        used_colors += [second for first, second in colored if first == nbr]
    polygon[1]=[color for color in range(10) if color not in used_colors][0]
    colored.append(polygon)

I would like the script to reduce at maximum the number of colors (4 or 5) using brute force loop or other method.

Here is the graph coloring code extracted from the code in my math.stackexchange answer. That code is admittedly a little complicated because it uses Knuth's Algorithm X for two tasks: solving the tromino packing problem, and then solving the coloring problem for each packing problem solution.

This code finds 24 3-color solutions or 756 4-color for your sample data in less than 1 second. There are actually more solutions, but I arbitrarily set the colors for the first 2 nodes to reduce the total number of solutions, and to speed up the searching algorithm a little.

''' Knuth's Algorithm X for the exact cover problem, 
    using dicts instead of doubly linked circular lists.
    Written by Ali Assaf

    From http://www.cs.mcgill.ca/~aassaf9/python/algorithm_x.html

    and http://www.cs.mcgill.ca/~aassaf9/python/sudoku.txt

    Python 2 / 3 compatible

    Graph colouring version
'''

from __future__ import print_function
from itertools import product

def solve(X, Y, solution):
    if not X:
        yield list(solution)
    else:
        c = min(X, key=lambda c: len(X[c]))
        Xc = list(X[c])

        for r in Xc:
            solution.append(r)
            cols = select(X, Y, r)
            for s in solve(X, Y, solution):
                yield s
            deselect(X, Y, r, cols)
            solution.pop()

def select(X, Y, r):
    cols = []
    for j in Y[r]:
        for i in X[j]:
            for k in Y[i]:
                if k != j:
                    X[k].remove(i)
        cols.append(X.pop(j))
    return cols

def deselect(X, Y, r, cols):
    for j in reversed(Y[r]):
        X[j] = cols.pop()
        for i in X[j]:
            for k in Y[i]:
                if k != j:
                    X[k].add(i)

#Invert subset collection
def exact_cover(X, Y):
    newX = dict((j, set()) for j in X)
    for i, row in Y.items():
        for j in row:
            newX[j].add(i)
    return newX

def colour_map(nodes, edges, ncolours=4):
    colours = range(ncolours)

    #The edges that meet each node
    node_edges = dict((n, set()) for n in nodes)
    for e in edges:
        n0, n1 = e
        node_edges[n0].add(e)
        node_edges[n1].add(e)

    for n in nodes:
        node_edges[n] = list(node_edges[n])

    #Set to cover
    coloured_edges = list(product(colours, edges))
    X = nodes + coloured_edges

    #Subsets to cover X with
    Y = dict()
    #Primary rows
    for n in nodes:
        ne = node_edges[n]
        for c in colours:
            Y[(n, c)] = [n] + [(c, e) for e in ne]

    #Dummy rows
    for i, ce in enumerate(coloured_edges):
        Y[i] = [ce]

    X = exact_cover(X, Y)

    #Set first two nodes
    partial = [(nodes[0], 0), (nodes[1], 1)]
    for s in partial:
        select(X, Y, s)

    for s in solve(X, Y, []):
        s = partial + [u for u in s if not isinstance(u, int)]
        s.sort()
        yield s

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

input_polygons = [
    (0, None), (1, None), (2, None), (3, None), (4, None), 
    (5, None), (6, None), (7, None), (8, None),
]

#graph array of neighbors ids
adjacent_polygons_graph = [
    (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (1, 0), (2, 0), (2, 3), 
    (3, 0), (3, 2), (3, 6), (3, 8), (4, 0), (4, 5), (4, 6), 
    (5, 0), (5, 4), (6, 3), (6, 4), (6, 7), (7, 6), (8, 3),
]

# Extract the nodes list 
nodes = [t[0] for t in input_polygons]

# Get an iterator of all solutions with 3 colours
all_solutions = colour_map(nodes, adjacent_polygons_graph, ncolours=3)

# Print all solutions
for count, solution in enumerate(all_solutions, start=1):
    print('%2d: %s' % (count, solution))

output

 1: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 1), (7, 0), (8, 1)]
 2: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 1), (7, 0), (8, 0)]
 3: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 1), (7, 2), (8, 1)]
 4: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 1), (7, 2), (8, 0)]
 5: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 0), (7, 2), (8, 0)]
 6: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 0), (7, 1), (8, 0)]
 7: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 0), (7, 2), (8, 1)]
 8: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 2), (5, 1), (6, 0), (7, 1), (8, 1)]
 9: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 1), (5, 2), (6, 0), (7, 1), (8, 1)]
10: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 1), (5, 2), (6, 0), (7, 1), (8, 0)]
11: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 1), (5, 2), (6, 0), (7, 2), (8, 0)]
12: [(0, 0), (1, 1), (2, 1), (3, 2), (4, 1), (5, 2), (6, 0), (7, 2), (8, 1)]
13: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 2), (7, 0), (8, 0)]
14: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 2), (5, 1), (6, 0), (7, 2), (8, 0)]
15: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 0), (7, 2), (8, 0)]
16: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 2), (5, 1), (6, 0), (7, 1), (8, 0)]
17: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 0), (7, 1), (8, 0)]
18: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 2), (7, 1), (8, 0)]
19: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 2), (7, 1), (8, 2)]
20: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 2), (7, 0), (8, 2)]
21: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 2), (5, 1), (6, 0), (7, 2), (8, 2)]
22: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 0), (7, 2), (8, 2)]
23: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 2), (5, 1), (6, 0), (7, 1), (8, 2)]
24: [(0, 0), (1, 1), (2, 2), (3, 1), (4, 1), (5, 2), (6, 0), (7, 1), (8, 2)]

If you just want a single solution, replace the last for loop with

output_polygons = next(all_solutions)
print(output_polygons)

---

FWIW, here's some code that can be used to create a diagram from the graph data.

# Create a Graphviz DOT file from graph data
colors = {
    None: 'white',
    0: 'red', 1: 'yellow', 
    2: 'green', 3: 'blue'
}

def graph_to_dot(nodes, edges, outfile=sys.stdout):
    outfile.write('strict graph test{\n')
    outfile.write('    node [style=filled];\n')

    for n, v in nodes:
        outfile.write('    {} [fillcolor={}];\n'.format(n, colors[v]))   
    outfile.write('\n')   

    for u, v in edges:
        outfile.write('    {} -- {};\n'.format(u, v))
    outfile.write('}\n')

You can call it like this:

# Produce a DOT file of the colored graph
with open('graph.dot', 'w') as f:
    graph_to_dot(output_polygons, adjacent_polygons_graph, f)

It can also be used to produce a DOT file of the input_polygons.

To generate a PNG image file from the DOT file using the Graphviz neato utility, you can do this (in Bash).

neato -Tpng -ograph.png graph.dot

Here's the DOT file for the first solution above:

strict graph test{
    node [style=filled];
    0 [fillcolor=red];
    1 [fillcolor=yellow];
    2 [fillcolor=yellow];
    3 [fillcolor=green];
    4 [fillcolor=green];
    5 [fillcolor=yellow];
    6 [fillcolor=yellow];
    7 [fillcolor=red];
    8 [fillcolor=yellow];

    0 -- 1;
    0 -- 2;
    0 -- 3;
    0 -- 4;
    0 -- 5;
    1 -- 0;
    2 -- 0;
    2 -- 3;
    3 -- 0;
    3 -- 2;
    3 -- 6;
    3 -- 8;
    4 -- 0;
    4 -- 5;
    4 -- 6;
    5 -- 0;
    5 -- 4;
    6 -- 3;
    6 -- 4;
    6 -- 7;
    7 -- 6;
    8 -- 3;
}

And here's the PNG:

Here's a live version, running on the SageMathCell server.