Linear Algebra in Room Escape

# coding: utf-8

"""
The following code based on
https://github.com/pmneila/Lights-Out/blob/master/lightsout.py
"""

from operator import add
from itertools import chain, combinations
import numpy as np
from scipy import ndimage

# First, we define operation in Galois Field (https://en.wikipedia.org/wiki/Finite_field)
# Here what we need is an Z/2Z (A mod 2 Galois Fields)
class GF2(object):
    """Galois field GF(2)."""

    def __init__(self, a=0):
        self.value = int(a) & 1

    def __add__(self, rhs):
        return GF2(self.value + GF2(rhs).value)

    def __mul__(self, rhs):
        return GF2(self.value * GF2(rhs).value)

    def __sub__(self, rhs):
        return GF2(self.value - GF2(rhs).value)

    def __div__(self, rhs):
        return GF2(self.value / GF2(rhs).value)

    def __repr__(self):
        return str(self.value)

    def __eq__(self, rhs):
        if isinstance(rhs, GF2):
            return self.value == rhs.value
        return self.value == rhs

    def __le__(self, rhs):
        if isinstance(rhs, GF2):
            return self.value <= rhs.value
        return self.value <= rhs

    def __lt__(self, rhs):
        if isinstance(rhs, GF2):
            return self.value < rhs.value
        return self.value < rhs

    def __int__(self):
        return self.value

    def __long__(self):
        return self.value

# Encapsulate operation for vectorization computation
GF2array = np.vectorize(GF2)

def gjel(A):
    """Gauss-Jordan elimination."""
    nulldim = 0
    for i in xrange(len(A)):
        pivot = A[i:,i].argmax() + i
        if A[pivot,i] == 0:
            nulldim = len(A) - i
            break
        row = A[pivot] / A[pivot,i]
        A[pivot] = A[i]
        A[i] = row

        for j in xrange(len(A)):
            if j == i:
                continue
            A[j] -= row*A[j,i]
    return A, nulldim

def GF2inv(A):
    """Inversion(逆) and eigenvectors(特征向量) of the null-space of a GF2 matrix."""
    n = len(A)
    assert n == A.shape[1], "Matrix must be square"

    A = np.hstack([A, np.eye(n)])

    B, nulldim = gjel(GF2array(A))

    inverse = np.int_(B[-n:, -n:])
    E = B[:n, :n]
    null_vectors = []
    if nulldim > 0:
        null_vectors = E[:, -nulldim:]
        null_vectors[-nulldim:, :] = GF2array(np.eye(nulldim))
        null_vectors = np.int_(null_vectors.T)

    print "inverse of matrix:"
    print inverse
    print "eigenvectors of matrix:"
    print null_vectors
    return inverse, null_vectors

def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

def lightsoutbase(n):
    # base of the lights out problem of size (n, n)
    a = np.eye(n*n)
    a = np.reshape(a, (n*n, n, n))

    # construct the A_{ij} Matrix
    a = np.array(map(ndimage.binary_dilation, a))
    lo_base = np.reshape(a, (n*n, n*n))
    return lo_base

class LightsOut(object):
    """Lights-Out solver."""
    
    def __init__(self, size=5):
        self.n = size
        self.base = lightsoutbase(self.n)
        self.invbase, self.null_vectors = GF2inv(self.base)
    
    def solve(self, b):
        b = np.asarray(b)
        assert b.shape[0] == b.shape[1] == self.n, "incompatible shape"
        
        if not self.issolvable(b):
            raise ValueError, "The given setup is not solvable"
        
        # Find the base solution.
        first = np.dot(self.invbase, b.ravel()) & 1
        
        # Given a solution, we can find more valid solutions
        # adding any combination of the null vectors.
        # Find the solution with the minimum number of 1's.
        solutions = [(first + reduce(add, nvs, 0))&1 for nvs in powerset(self.null_vectors)]
        final = min(solutions, key=lambda x: x.sum())
        return np.reshape(final, (self.n,self.n))
    
    def issolvable(self, b):
        """Determine if the given configuration is solvable.
        
        A configuration is solvable if it is orthogonal to
        the null vectors of the base.
        """
        b = np.asarray(b)
        assert b.shape[0] == b.shape[1] == self.n, "incompatible shape"
        b = b.ravel()
        p = map(lambda x: np.dot(x,b)&1, self.null_vectors)
        return not any(p)

if __name__ == "__main__":

    b = np.array([[1,0,1],
                  [0,0,1],
                  [1,0,0]])
    print "\nlights status:"
    print b

    lo = LightsOut(3)

    bsol = lo.solve(b)

    print "\nThe solution of"
    print b
    print "is"
    print bsol