grid_zero.py

import random
import numpy as np
#from cvxpy import *

def p(M):
    """Method to pretty print our matrix. """
    for i in range(len(M)):
        print(' '.join([str(num) for num in M[i]]))

def generate_M(N=3, max_touches=2):
    """Generate a solvable matrix. """

    M = [[0]*N for _ in range(N)] # Our matrix.

    # NOTE: Some touches can be rewritten as a set of other touches (if the game board length is odd).
    # Example: in the case of a 3x3 game board, touching the lights
    # in the first and second diagram produces equivalent game states.
    # 0 x 0
    # 0 x 0
    # x 0 0
    # and
    # 0 0 0
    # 0 0 0
    # 0 0 x
    # are equivalent (x represents the square that you touch).
    # This means that, although unlikely, we can accidentally make a solvable matrix
    # (when the game board lenght is odd) that needs fewer touches than we anticipate.

    touches = set() # Unique touches applied to our matrix.
    for _ in range(max_touches):
        i, j = random.randint(0, N-1), random.randint(0, N-1) # Pick a square at random.
        while (i, j) in touches: # Make sure the touch hasn't been applied before.
            i, j = random.randint(0, N-1), random.randint(0, N-1)
        touches.add((i, j))
        touch(M, i, j)
    return M

def touch(M, i, j):
    # Flip the values on the horizontal and the vertical.
    for a in range(len(M)):
        M[i][a] = 1 - M[i][a]
    for b in range(len(M)):
        M[b][j] = 1 - M[b][j]
    M[i][j] = 1 - M[i][j]

def make_A(n=2):
    eye = np.identity(n)
    ones = np.ones(n**2).reshape((n,n))
    # Ones matrix along the diagonal and identities everywhere else.
    A = [[eye if i !=j else ones for i in range(n)] for j in range(n)]
    return np.bmat(A) # Puts together blocks of matrices.

def _find_row_with_1_in_column(A, pivot_row, pivot_column):
    for row in range(pivot_row, len(A)):
        if A[row, pivot_column] == 1:
            return row
    return -1

# Quick unit test.
A = make_A(2)
assert _find_row_with_1_in_column(A, 0, 0) == 0
assert _find_row_with_1_in_column(A, 2, 1) == 3

def solve_Ax_eq_b_mod_p(A0, b0, p=2):
    """Solves Ax=b mod p using Gaussian elimination.

    Gaussian elimination does swaps of rows and adds one row to another.
    We augment A with b and perform all operations that were done on A on b also. """

    A = A0.copy()
    b = b0.copy()

    N = len(A)

    pivot = (0, 0)  # Start at top-left corner.
    while pivot != (N, N):
        pivot_row, pivot_column = pivot

        # We want to find the row where we have a 1 in this column and swap this columns.
        # So if pivot is (1, 1) and the matrix is
        # 1 1 1
        # 0 0 0
        # 0 1 0.
        # We want to swap rows 1 and 2:
        # 1 1 1
        # 0 1 0
        # 0 0 0.
        row_with_1_in_column = _find_row_with_1_in_column(A, pivot_row, pivot_column)
        if row_with_1_in_column != pivot_row:
            # The following swapping DOES NOT work without copy.
            A[pivot_row], A[row_with_1_in_column] = A[row_with_1_in_column], A[pivot_row].copy()

            # Augmented operation on b.
            b[pivot_row], b[row_with_1_in_column] = b[row_with_1_in_column], b[pivot_row].copy()

        # Now we can use this 1 in this column to cancel out everything above and below (mod p).
        # Go from top to bottom (except this row) and add the inverse of current row.
        for row in range(N):
            if row != pivot_row and A[row, pivot_column] == 1:
                A[row] = (A[row] - A[pivot_row]) % p

                # Augmented operation on b.
                b[row] = (b[row] - b[pivot_row]) % p

        # Move pivot down the diagonal.
        pivot = pivot_row + 1, pivot_column + 1

    if (np.isclose(A0.dot(b) % p, b0)).all():
        # At this point, if A is invertible, then A becomes the identity matrix and
        # so b becomes the solution mod p. The solution is unique.
        if not N & 1: # If N is even...
            return A, b
        # Otherwise we might have gotten lucky and x is in the span of the rows of A.
        # What remains is to find the solution with the least-norm (since these are not unique).
        else:
            # TODO: Here is where we need to return the least-norm solution.
            # Right now b is a solution, but not guaranteed to be the least-norm solution.
            return A, b

            # An attempt to solve the problem using the library cvxpy (doesn't work right now).
            """
            x = Variable(len(b0))
            objective = Minimize(sum(x))
            constraints = [A0 * x == b0, x >= 0, x <= 1]
            prob = Problem(objective, constraints)

            result = prob.solve()
            print('status: %s' % prob.status)
            print('result: %f' % result)
            return A, x.value
            """
    else:
        # Because we are explicitly constructing solvable cases,
        # for the time being, we shouldn't hit this else statement (or something is wrong).
        raise AssertionError('Encountered an unsolvable case. ')

def answer(M):
    n = len(M)
    A = make_A(n)
    b = np.array(M).transpose().reshape(n**2)
    A1, x = solve_Ax_eq_b_mod_p(A, b)
    return A1, sum(x)

if __name__ == '__main__':
    for _ in range(5):
        N = 6 # Length of the side of the matrix.
        num_touches = 4
        M = generate_M(N, num_touches)
        p(M)
        _, x = answer(M)
        print('Required %i touches. ' % x)