linear_algebra_game.py

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