chessboard.py

from fractions import Fraction
from copy import deepcopy

def places():
    for r in range(1, 8+1):
        for c in range(1, 8+1):
            yield (r, c)

def neighbors(b1, b2):
    r1, c1 = b1
    r2, c2 = b2
    if (r1, c1) == (r2, c2):
        return False
    return abs(r1-r2) <= 1 and abs(c1-c2) <= 1

def transitions():
    num_neighbors = {}
    for b1 in places():
        num_neighbors[b1] = len(list(b2 for b2 in places() if neighbors(b1, b2)))
    trans = {}
    for b1 in places():
        for b2 in places():
            if neighbors(b1, b2):
                trans[(b1,b2)] = Fraction(1, num_neighbors[b1])
            else:
                trans[(b1,b2)] = Fraction(0)
    return trans

def transition_matrix(trans):
    P = []
    for b1 in places():
        row = []
        for b2 in places():
            row.append(trans[(b1,b2)])
        P.append(row)
    return P

def verify_transition_matrix(P):
    n = len(P)
    for r in range(n):
        row_sum = Fraction(0)
        for c in range(n):
            row_sum += P[r][c]
        assert row_sum == 1, f'every row should sum to 1, but row {r} sums to {row_sum}'

def prepare_for_solving(P):
    n = len(P)
    # transpose
    for r in range(n):
        for c in range(r+1, n):
            P[r][c], P[c][r] = P[c][r], P[r][c]
    # subtract 1 from the diagonal
    for r in range(n):
        P[r][r] -= 1
    # add the results column (zeros)
    for r in range(n):
        P[r].append(0)
    # add the sum-check row (sum of all the results must be 1)
    P.append([1 for _ in range(n+1)])

def gaussian_eliminate(P):
    def swap_rows(r0, r1):
        n = len(P)
        for c in range(n):
            P[r0][c], P[r1][c] = P[r1][c], P[r0][c]

    def multiply_row(r, k):
        n = len(P)
        for c in range(n):
            P[r][c] *= k

    def add_row_multiple(r0, r1, k):
        n = len(P)
        for c in range(n):
            P[r1][c] += P[r0][c] * k
 
    n = len(P)
    for r in range(n):
        if P[r][r] == 0:
            for r1 in range(r+1, n):
                if P[r1][r] != 0:
                    swap_rows(r, r1)
                    break
            else:
                continue
        multiply_row(r, 1/P[r][r])
        for r1 in range(n):
            if r == r1:
                continue
            add_row_multiple(r, r1, -P[r1][r])

def extract_result(P):
    n = len(P)-1
    return [P[r][n] for r in range(n)]

def multiply(vec, P):
    n = len(P)
    res = [Fraction(0) for _ in range(n)]
    for j in range(n):
        for i in range(n):
            res[j] += vec[i] * P[i][j]
    return res

P = transition_matrix(transitions())
verify_transition_matrix(P)
S = deepcopy(P)
prepare_for_solving(S)
gaussian_eliminate(S)
result = extract_result(S)
verify = multiply(result, P)

print('result:', ' '.join(map(str, result)))
print()
print('verify:', ' '.join(map(str, verify)))
print()

print('Stationary distribution visualization on the chess board:')
board = [[None for _ in range(8)] for _ in range(8)]
for (i, (r, c)) in enumerate(places()):
    board[r-1][c-1] = format(str(result[i]), '>6')
for row in board:
    print(' '.join(row))