Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
linear_algebra_game.py
# Supporting code for this blog article:
# https://thirld.com/blog/2018/05/17/linear-algebra-solves-a-childhood-mystery/
#
# License: public domain.
import itertools
import numpy as np
def pretty_print_matrix(name, M):
"""
Prints a numpy array, replacing zeros with spaces.
"""
print "{}:".format(name)
print str(M).replace('0', ' ')
print
def make_action_matrix(board_size):
"""
Returns the action matrix for a board of size board_size x board_size.
M[:, a] has 1s for the cells flipped by action a.
"""
affected_cell_rel_coords = [(0, 0), (-1, 0), (1, 0), (0, -1), (0, 1)]
num_cells = num_actions = board_size * board_size
M = np.zeros((num_cells, num_actions), dtype=np.int)
for r, c in itertools.product(range(board_size), repeat=2):
action_id = r * board_size + c
for dr, dc in affected_cell_rel_coords:
ar, ac = r + dr, c + dc
if 0 <= ar < board_size and 0 <= ac < board_size:
affected_cell_id = ar * board_size + ac
M[affected_cell_id, action_id] = 1
return M
def get_reduced_row_echelon_form_z2(M):
"""
Computes the reduced row echelon form of M in Z2.
Returns (K, R), where R is the reduced row echelon form of M, and K is the
transformation matrix, s.t. R = np.dot(K, M).
This is a "paper and pencil" algorithm; it is polynomial but not fast.
"""
num_rows, num_cols = M.shape
# Holds the transformation matrix s.t. R = np.dot(K, M).
K = np.eye(num_rows, dtype=np.int)
# Holds the reduced row echelon matrix.
R = M.copy()
# Returns the leading index of row r of R, or num_rows if the row is all 0.
def L(r):
for i, val in enumerate(R[r, :]):
if val == 1:
return i
return num_rows
# Adds row a to row r in R and K.
def add_row_to(a, r):
R[r, :] = (R[r, :] + R[a, :]) % 2
K[r, :] = (K[r, :] + K[a, :]) % 2
# Checks that R = np.dot(K, M) holds true.
def check_invariants():
assert (np.dot(K, M) % 2 == R).all()
# Step 1: Convert to row echelon form (this is not unique).
for r in xrange(num_rows):
while L(r) < num_rows:
for r_to_add in xrange(num_rows):
if r_to_add != r and L(r_to_add) == L(r):
add_row_to(r_to_add, r)
check_invariants()
break
else:
break # No more reductions possible on this row.
row_indices = np.argsort([L(r) for r in xrange(num_rows)])
R = R[row_indices, :]
K = K[row_indices, :]
check_invariants()
# Step 2: Convert to reduced row echelon form (this is unique).
for r in xrange(num_rows):
for c in xrange(L(r) + 1, num_cols):
if R[r, c] != 1:
continue
for r_to_add in xrange(num_rows):
if r_to_add != r and L(r_to_add) == c:
add_row_to(r_to_add, r)
check_invariants()
break
return R, K
def z2_vector_to_int(arr):
"""
Returns an int where the i-th bit is on iff arr[i] is 1.
"""
result = 0
for i in xrange(len(arr)):
assert arr[i] in (0, 1)
result |= arr[i] << i
return result
def int_to_z2_vector(arr_int, size):
"""
Inverse of the above: Converts an int to an array in Z2 of the given size.
"""
result = np.zeros(size, dtype=np.int)
for i in xrange(size):
result[i] = (arr_int & (1 << i)) > 0
return result
def check_all_boards_bfs(M):
"""
Checks all boards for solvability using breadth-first search.
M is the action matrix. Returns an array is_solvable, where is_solvable[b]
is True iff the board b is solvable. Each board is represented as an int,
where each bit represents the status of a cell.
"""
# Represent each action as an int.
num_cells, num_actions = M.shape
actions = []
for c in xrange(num_actions):
actions.append(z2_vector_to_int(M[:, c]))
# is_solvable[b] is True iff board b (represented as an int) is solvable.
num_boards = 2 ** num_cells
is_solvable = np.zeros(num_boards, dtype=np.bool)
# Breadth-first search starting from the empty board.
is_solvable[0] = True
boards_to_expand = [0]
num_iters = 0
while len(boards_to_expand):
num_iters += 1
if num_iters % 100000 == 0:
print "num_iters={}".format(num_iters)
board = boards_to_expand.pop()
assert is_solvable[board]
for action in actions:
new_board = board ^ action
if not is_solvable[new_board]:
is_solvable[new_board] = True
boards_to_expand.append(new_board)
return is_solvable
def get_span_of_null_space_z2(R):
"""
Computes the null space of R in Z2.
R must be in reduced row echelon form. Returns a matrix N, whose columns
span null(R). This is a "paper and pencil" algorithm.
"""
num_rows, num_cols = R.shape
# Returns the leading index of row r of R, or num_rows if the row is all 0.
def L(r):
for i, val in enumerate(R[r, :]):
if val == 1:
return i
return num_rows
# Find the leading index for each row.
leading_indices = [L(r) for r in xrange(num_rows)]
# Find the column indices that do not have leading ones.
free_var_col_indices = [
c for c in xrange(num_cols) if c not in leading_indices]
# Construct the span matrix.
span = np.zeros((num_cols, len(free_var_col_indices)), dtype=np.int)
for c in xrange(num_cols):
try:
free_var_span_index = free_var_col_indices.index(c)
# This is a free variable. Express it as equal to itself.
span[c, free_var_span_index] = 1
except ValueError:
# This is not a free variable. Express it as a sum of the free
# variables.
row_index = leading_indices.index(c)
for i, col_index in enumerate(free_var_col_indices):
span[c, i] = R[row_index, col_index]
return span
def check_all_boards_checksums(c1, c2):
"""
Checks all boards for solvability using the checksum approach.
c1 and c2 are the checksums described in the article. Returns an array
is_solvable, where is_solvable[b] is True iff the board b is solvable.
"""
assert len(c1) == len(c2) > 0
num_boards = 2 ** len(c1)
c1_int = z2_vector_to_int(c1)
c2_int = z2_vector_to_int(c2)
is_solvable = np.zeros(num_boards, dtype=np.bool)
for b_int in xrange(num_boards):
if b_int % 100000 == 0:
print "b_int={}".format(b_int)
if (bin(b_int & c1_int).count('1') % 2 == 0 and
bin(b_int & c2_int).count('1') % 2 == 0):
is_solvable[b_int] = True
return is_solvable
if __name__ == "__main__":
np.set_printoptions(linewidth=200)
# Part 1: Show the action matrix.
M = make_action_matrix(board_size=5)
pretty_print_matrix("M", M)
# Part 2: Show the reduced row echelon form and the transformation matrix.
R, K = get_reduced_row_echelon_form_z2(M)
pretty_print_matrix("R", R)
pretty_print_matrix("K", K)
assert (np.dot(K, M) % 2 == R).all()
# # Part 3: Check the claim that 75% of boards are unsolvable.
is_solvable = check_all_boards_bfs(M)
num_solvable_boards = np.sum(is_solvable)
num_total_boards = len(is_solvable)
print "{} of {} boards ({}%) are solvable".format(
num_solvable_boards, num_total_boards,
100.0 * num_solvable_boards / num_total_boards)
print
# Part 4: Show the span of the null space of the action matrix.
null_span = get_span_of_null_space_z2(R)
pretty_print_matrix("null_span", null_span)
# There are only 3 values in the null space. Check that they are solvable.
assert null_span.shape == (25, 2)
n1 = null_span[:, 0]
n2 = null_span[:, 1]
n3 = (n1 + n2) % 2
assert (np.dot(M, n1) % 2 == 0).all()
assert (np.dot(M, n2) % 2 == 0).all()
assert (np.dot(M, n3) % 2 == 0).all()
assert is_solvable[z2_vector_to_int(n1)]
assert is_solvable[z2_vector_to_int(n2)]
assert is_solvable[z2_vector_to_int(n3)]
# Part 5: Verify the "checksum method" for checking whether a board is
# solvable. The c1 and c2 in the article are in fact members of null_span.
c1 = n2
c2 = n3
pretty_print_matrix("c1", c1.reshape((5, 5)))
pretty_print_matrix("c2", c2.reshape((5, 5)))
print
is_solvable_alt = check_all_boards_checksums(c1, c2)
assert (is_solvable == is_solvable_alt).all()
print
# Some quick tests.
num_cells = M.shape[0]
num_boards = 2 ** num_cells
for value in (0, 1, 2, 345, 67890, num_boards - 1):
assert z2_vector_to_int(int_to_z2_vector(value, num_cells)) == value
# Part 6: Verify the solution method described in the article.
# TODO(cberzan): It should be possible to make this loop (and the
# check_all_boards_* functions) a lot faster by using PyPy and/or by
# replacing some of the numpy operations with bit operations.
for b_int in xrange(num_boards):
if b_int % 100000 == 0:
print "b_int={}".format(b_int)
b = int_to_z2_vector(b_int, num_cells)
d = np.dot(K, b) % 2
if d[-1] == 0 and d[-2] == 0:
assert is_solvable[b_int]
assert (np.dot(M, d) % 2 == b).all()
else:
assert not is_solvable[b_int]
print
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment