
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt


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


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')


pos = nx.circular_layout(G)


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)

True


states_visited = []
T_steps = 100000
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="...")
    # base chain proposal
    n_t = np.random.choice(G.nodes)
    c_t = np.random.choice(available_colors)
    # proposal from the base chain
    proposal = {**current_colors, n_t: c_t}
    # reject/accept the move based on the metropolis rule
    if is_proper_coloring(proposal, G):
        # accept the move
        current_colors = proposal
        print("accept the move")
    else:
        # reject the move
        current_colors = current_colors
        print("reject the move")
    states_visited.append(tuple(current_colors[n] for n in G.nodes))
    if t< 10:
        print('random node', n_t)
        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...accept the move
random node 4
random color red

1...reject the move
random node 1
random color red

2...accept the move
random node 1
random color green

accept the move
random node 2
random color blue

4...accept the move
random node 3
random color yellow

accept the move
random node 3
random color green

accept the move
random node 2
random color red

accept the move
random node 3
random color yellow

8...accept the move
random node 4
random color red

accept the move
random node 2
random color red


from collections import Counter
counts = Counter(states_visited)
counts

Counter({('blue', 'red', 'yellow', 'red'): 1134,
         ('green', 'red', 'yellow', 'red'): 1163,
         ('green', 'blue', 'yellow', 'red'): 1298,
         ('green', 'blue', 'green', 'red'): 1286,
         ('green', 'red', 'green', 'red'): 1221,
         ('yellow', 'blue', 'yellow', 'red'): 1274,
         ('yellow', 'green', 'yellow', 'red'): 1241,
         ('yellow', 'green', 'yellow', 'blue'): 1095,
         ('red', 'green', 'yellow', 'blue'): 1099,
         ('red', 'blue', 'yellow', 'blue'): 1160,
         ('red', 'blue', 'green', 'blue'): 1153,
         ('green', 'blue', 'yellow', 'blue'): 1127,
         ('yellow', 'blue', 'yellow', 'blue'): 1280,
         ('yellow', 'red', 'yellow', 'blue'): 1142,
         ('yellow', 'red', 'yellow', 'red'): 1196,
         ('yellow', 'red', 'yellow', 'green'): 1263,
         ('yellow', 'blue', 'yellow', 'green'): 1285,
         ('red', 'blue', 'yellow', 'green'): 1268,
         ('yellow', 'red', 'blue', 'red'): 1226,
         ('blue', 'red', 'blue', 'red'): 1196,
         ('blue', 'green', 'blue', 'red'): 1223,
         ('yellow', 'green', 'blue', 'red'): 1222,
         ('yellow', 'blue', 'green', 'red'): 1126,
         ('yellow', 'green', 'yellow', 'green'): 1232,
         ('blue', 'red', 'yellow', 'green'): 1310,
         ('yellow', 'blue', 'red', 'green'): 1370,
         ('red', 'blue', 'red', 'green'): 1252,
         ('yellow', 'green', 'red', 'green'): 1237,
         ('blue', 'green', 'yellow', 'green'): 1125,
         ('blue', 'green', 'red', 'green'): 1208,
         ('blue', 'yellow', 'red', 'green'): 1320,
         ('blue', 'green', 'blue', 'green'): 1110,
         ('blue', 'green', 'blue', 'yellow'): 1255,
         ('blue', 'green', 'red', 'yellow'): 1141,
         ('red', 'green', 'red', 'yellow'): 1237,
         ('red', 'green', 'red', 'blue'): 1093,
         ('yellow', 'green', 'red', 'blue'): 1210,
         ('red', 'green', 'red', 'green'): 1193,
         ('red', 'yellow', 'red', 'green'): 1339,
         ('red', 'yellow', 'red', 'blue'): 1202,
         ('red', 'yellow', 'red', 'yellow'): 1250,
         ('red', 'green', 'yellow', 'green'): 1145,
         ('red', 'blue', 'red', 'yellow'): 1326,
         ('red', 'blue', 'red', 'blue'): 1085,
         ('yellow', 'green', 'blue', 'green'): 1293,
         ('red', 'blue', 'green', 'yellow'): 1157,
         ('green', 'blue', 'green', 'yellow'): 1220,
         ('green', 'blue', 'green', 'blue'): 1200,
         ('green', 'red', 'green', 'yellow'): 1222,
         ('green', 'red', 'green', 'blue'): 1203,
         ('green', 'red', 'blue', 'yellow'): 1207,
         ('green', 'red', 'yellow', 'blue'): 1304,
         ('yellow', 'blue', 'red', 'blue'): 1316,
         ('red', 'yellow', 'green', 'yellow'): 1142,
         ('red', 'yellow', 'green', 'blue'): 1111,
         ('green', 'yellow', 'green', 'blue'): 1061,
         ('green', 'yellow', 'red', 'blue'): 1239,
         ('green', 'blue', 'red', 'blue'): 1195,
         ('green', 'yellow', 'green', 'red'): 1081,
         ('green', 'yellow', 'green', 'yellow'): 1129,
         ('green', 'yellow', 'blue', 'yellow'): 1141,
         ('red', 'yellow', 'blue', 'yellow'): 1199,
         ('red', 'green', 'blue', 'yellow'): 1287,
         ('blue', 'red', 'blue', 'yellow'): 1210,
         ('blue', 'yellow', 'blue', 'yellow'): 1133,
         ('red', 'yellow', 'blue', 'green'): 1323,
         ('red', 'green', 'blue', 'green'): 1310,
         ('blue', 'red', 'blue', 'green'): 1244,
         ('yellow', 'red', 'blue', 'green'): 1206,
         ('yellow', 'blue', 'green', 'blue'): 1192,
         ('green', 'blue', 'red', 'yellow'): 1107,
         ('green', 'yellow', 'red', 'yellow'): 1062,
         ('blue', 'yellow', 'red', 'yellow'): 1162,
         ('blue', 'yellow', 'blue', 'green'): 1137,
         ('blue', 'yellow', 'blue', 'red'): 988,
         ('green', 'yellow', 'blue', 'red'): 1094,
         ('green', 'red', 'blue', 'red'): 1069,
         ('yellow', 'red', 'green', 'blue'): 1142,
         ('blue', 'yellow', 'green', 'red'): 912,
         ('blue', 'yellow', 'green', 'yellow'): 1116,
         ('blue', 'red', 'green', 'red'): 1167,
         ('blue', 'green', 'yellow', 'red'): 1143,
         ('yellow', 'red', 'green', 'red'): 1101,
         ('blue', 'red', 'green', 'yellow'): 1157})


# proportion of time points spent at each coloring
pi_estimated = np.array(list(counts.values()))/T_steps
pi_estimated

array([0.01134, 0.01163, 0.01298, 0.01286, 0.01221, 0.01274, 0.01241,
       0.01095, 0.01099, 0.0116 , 0.01153, 0.01127, 0.0128 , 0.01142,
       0.01196, 0.01263, 0.01285, 0.01268, 0.01226, 0.01196, 0.01223,
       0.01222, 0.01126, 0.01232, 0.0131 , 0.0137 , 0.01252, 0.01237,
       0.01125, 0.01208, 0.0132 , 0.0111 , 0.01255, 0.01141, 0.01237,
       0.01093, 0.0121 , 0.01193, 0.01339, 0.01202, 0.0125 , 0.01145,
       0.01326, 0.01085, 0.01293, 0.01157, 0.0122 , 0.012  , 0.01222,
       0.01203, 0.01207, 0.01304, 0.01316, 0.01142, 0.01111, 0.01061,
       0.01239, 0.01195, 0.01081, 0.01129, 0.01141, 0.01199, 0.01287,
       0.0121 , 0.01133, 0.01323, 0.0131 , 0.01244, 0.01206, 0.01192,
       0.01107, 0.01062, 0.01162, 0.01137, 0.00988, 0.01094, 0.01069,
       0.01142, 0.00912, 0.01116, 0.01167, 0.01143, 0.01101, 0.01157])


plt.bar([a for a in range(len(pi_estimated))], pi_estimated)

<BarContainer object of 84 artists>


for coloring in [0,1,2,3,4,5]:
    print(pi_estimated[coloring]**-1)

88.18342151675485
85.98452278589853
77.04160246533128
77.76049766718506
81.9000819000819
78.49293563579278


np.average(pi_estimated)**-1

84.0


# # Compare with the exact formula:
# https://math.stackexchange.com/a/3393053/484640
r = len(available_colors)
r*(r-1)*(r**2 - 3*r + 3)

84

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