In [1]:
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
In [2]:
def is_proper_coloring(colors, G):
    assert len(colors) == 4
    for node, other_node in G.edges:
        if colors[node] == colors[other_node]:
            return False
    return True
In [3]:
G = nx.Graph()
G.add_edge(1, 2)
G.add_edge(2, 3)
G.add_edge(3, 4)
G.add_edge(4, 1)
available_colors = ('blue', 'red', 'green', 'yellow')
In [4]:
pos = nx.circular_layout(G)
In [5]:
current_colors = {1: 'blue', 2: 'red', 3: 'yellow', 4: 'red'}
nx.draw(G, pos=pos, node_color=[current_colors[n] for n in G.nodes], with_labels=True)

is_proper_coloring(current_colors, G)
Out[5]:
True

Glauber dynamics with uniform distribution on proper colorings as the stationary distribution¶

In [ ]:
In [6]:
states_visited = []
T_steps = 30000
nx.draw(G, pos=pos, node_color=[current_colors[n] for n in G.nodes], with_labels=True)
plt.show()
for t in range(T_steps):
    if (t & (t-1) == 0): print(t, end="...")
    n_t = np.random.choice(G.nodes)
    allowed_colors = [c for c in available_colors
                        if  is_proper_coloring({**current_colors, n_t: c}, G)]
    c_t = np.random.choice(allowed_colors)
    #print("allowed for node", n_t, ":", allowed_colors)
    current_colors[n_t] = c_t
    states_visited.append(tuple(current_colors[n] for n in G.nodes))
    if t< 10:
        print('random node', n_t)
        print("allowed_colors", allowed_colors)
        print('random color', c_t)
        nx.draw(G, pos=pos, node_color=[current_colors[n] for n in G.nodes], with_labels=True)
        plt.show()

0...random node 2
allowed_colors ['red', 'green']
random color green

1...random node 4
allowed_colors ['red', 'green']
random color green

2...random node 2
allowed_colors ['red', 'green']
random color red

random node 3
allowed_colors ['blue', 'yellow']
random color yellow

4...random node 3
allowed_colors ['blue', 'yellow']
random color yellow

random node 4
allowed_colors ['red', 'green']
random color red

random node 1
allowed_colors ['blue', 'green', 'yellow']
random color green

random node 1
allowed_colors ['blue', 'green', 'yellow']
random color yellow

8...random node 1
allowed_colors ['blue', 'green', 'yellow']
random color green

random node 2
allowed_colors ['blue', 'red']
random color red

16...32...64...128...256...512...1024...2048...4096...8192...16384...
In [7]:
from collections import Counter
counts = Counter(states_visited)
counts
Out[7]:
Counter({('blue', 'green', 'yellow', 'red'): 318,
         ('blue', 'green', 'yellow', 'green'): 390,
         ('blue', 'red', 'yellow', 'green'): 401,
         ('blue', 'red', 'yellow', 'red'): 375,
         ('green', 'red', 'yellow', 'red'): 392,
         ('yellow', 'red', 'yellow', 'red'): 357,
         ('green', 'blue', 'yellow', 'red'): 367,
         ('green', 'blue', 'green', 'red'): 371,
         ('green', 'red', 'yellow', 'blue'): 368,
         ('yellow', 'blue', 'yellow', 'red'): 281,
         ('yellow', 'blue', 'green', 'red'): 323,
         ('yellow', 'blue', 'green', 'blue'): 323,
         ('yellow', 'blue', 'red', 'blue'): 333,
         ('yellow', 'blue', 'yellow', 'blue'): 315,
         ('yellow', 'green', 'yellow', 'blue'): 323,
         ('yellow', 'green', 'red', 'blue'): 339,
         ('green', 'blue', 'red', 'blue'): 305,
         ('yellow', 'red', 'green', 'blue'): 445,
         ('yellow', 'red', 'yellow', 'green'): 393,
         ('yellow', 'red', 'blue', 'green'): 416,
         ('blue', 'red', 'blue', 'green'): 413,
         ('blue', 'yellow', 'blue', 'green'): 408,
         ('blue', 'yellow', 'red', 'green'): 403,
         ('blue', 'yellow', 'red', 'yellow'): 352,
         ('blue', 'green', 'red', 'green'): 386,
         ('red', 'green', 'red', 'green'): 432,
         ('red', 'green', 'red', 'blue'): 372,
         ('red', 'blue', 'red', 'blue'): 375,
         ('red', 'blue', 'red', 'yellow'): 404,
         ('red', 'blue', 'green', 'yellow'): 363,
         ('red', 'green', 'red', 'yellow'): 420,
         ('red', 'yellow', 'red', 'yellow'): 341,
         ('red', 'yellow', 'blue', 'yellow'): 347,
         ('red', 'green', 'blue', 'yellow'): 416,
         ('red', 'green', 'blue', 'green'): 365,
         ('red', 'yellow', 'blue', 'green'): 374,
         ('red', 'yellow', 'red', 'green'): 386,
         ('red', 'green', 'yellow', 'green'): 402,
         ('blue', 'green', 'blue', 'yellow'): 335,
         ('blue', 'green', 'blue', 'green'): 384,
         ('blue', 'green', 'blue', 'red'): 353,
         ('yellow', 'green', 'blue', 'red'): 375,
         ('yellow', 'red', 'blue', 'red'): 404,
         ('blue', 'red', 'blue', 'red'): 365,
         ('blue', 'red', 'green', 'red'): 401,
         ('blue', 'yellow', 'green', 'red'): 400,
         ('green', 'yellow', 'green', 'red'): 304,
         ('green', 'yellow', 'blue', 'red'): 280,
         ('green', 'yellow', 'blue', 'yellow'): 311,
         ('red', 'blue', 'yellow', 'green'): 366,
         ('blue', 'yellow', 'blue', 'yellow'): 347,
         ('blue', 'yellow', 'green', 'yellow'): 350,
         ('green', 'yellow', 'red', 'yellow'): 331,
         ('green', 'yellow', 'red', 'blue'): 234,
         ('red', 'yellow', 'green', 'yellow'): 366,
         ('red', 'yellow', 'green', 'blue'): 323,
         ('red', 'blue', 'green', 'blue'): 352,
         ('yellow', 'red', 'yellow', 'blue'): 387,
         ('green', 'red', 'green', 'blue'): 378,
         ('green', 'blue', 'green', 'blue'): 304,
         ('green', 'blue', 'green', 'yellow'): 374,
         ('green', 'yellow', 'green', 'yellow'): 297,
         ('green', 'yellow', 'green', 'blue'): 266,
         ('green', 'blue', 'red', 'yellow'): 350,
         ('red', 'green', 'yellow', 'blue'): 361,
         ('yellow', 'green', 'yellow', 'green'): 407,
         ('yellow', 'blue', 'yellow', 'green'): 372,
         ('red', 'blue', 'yellow', 'blue'): 336,
         ('red', 'yellow', 'red', 'blue'): 305,
         ('green', 'blue', 'yellow', 'blue'): 307,
         ('yellow', 'red', 'green', 'red'): 386,
         ('green', 'red', 'green', 'red'): 331,
         ('green', 'red', 'blue', 'red'): 317,
         ('green', 'red', 'blue', 'yellow'): 309,
         ('yellow', 'green', 'yellow', 'red'): 336,
         ('yellow', 'green', 'blue', 'green'): 396,
         ('red', 'blue', 'red', 'green'): 433,
         ('yellow', 'blue', 'red', 'green'): 390,
         ('yellow', 'green', 'red', 'green'): 315,
         ('blue', 'red', 'blue', 'yellow'): 327,
         ('blue', 'red', 'green', 'yellow'): 359,
         ('blue', 'yellow', 'blue', 'red'): 304,
         ('blue', 'green', 'red', 'yellow'): 372,
         ('green', 'red', 'green', 'yellow'): 306})
In [8]:
# proportion of time points spent at each coloring
pi_estimated = np.array(list(counts.values()))/T_steps
pi_estimated
Out[8]:
array([0.0106    , 0.013     , 0.01336667, 0.0125    , 0.01306667,
       0.0119    , 0.01223333, 0.01236667, 0.01226667, 0.00936667,
       0.01076667, 0.01076667, 0.0111    , 0.0105    , 0.01076667,
       0.0113    , 0.01016667, 0.01483333, 0.0131    , 0.01386667,
       0.01376667, 0.0136    , 0.01343333, 0.01173333, 0.01286667,
       0.0144    , 0.0124    , 0.0125    , 0.01346667, 0.0121    ,
       0.014     , 0.01136667, 0.01156667, 0.01386667, 0.01216667,
       0.01246667, 0.01286667, 0.0134    , 0.01116667, 0.0128    ,
       0.01176667, 0.0125    , 0.01346667, 0.01216667, 0.01336667,
       0.01333333, 0.01013333, 0.00933333, 0.01036667, 0.0122    ,
       0.01156667, 0.01166667, 0.01103333, 0.0078    , 0.0122    ,
       0.01076667, 0.01173333, 0.0129    , 0.0126    , 0.01013333,
       0.01246667, 0.0099    , 0.00886667, 0.01166667, 0.01203333,
       0.01356667, 0.0124    , 0.0112    , 0.01016667, 0.01023333,
       0.01286667, 0.01103333, 0.01056667, 0.0103    , 0.0112    ,
       0.0132    , 0.01443333, 0.013     , 0.0105    , 0.0109    ,
       0.01196667, 0.01013333, 0.0124    , 0.0102    ])
In [9]:
plt.bar([a for a in range(len(pi_estimated))], pi_estimated)
Out[9]:
<BarContainer object of 84 artists>
In [10]:
for coloring in [0,1,2,3,4,5]:
    print(pi_estimated[coloring]**-1)

94.33962264150944
76.92307692307692
74.81296758104737
80.0
76.53061224489795
84.03361344537815
In [11]:
np.average(pi_estimated)**-1
Out[11]:
84.0
In [12]:
# # Compare with the exact formula:
# https://math.stackexchange.com/a/3393053/484640
r = len(available_colors)
r*(r-1)*(r**2 - 3*r + 3)
Out[12]:
84
In [ ]: