Enumeration of spanning subgraphs of the NxN grid

enumeratemazes.py

# Problem:
# How many possible mazes are there on an NxN grid? Where a maze is a
# collection of links between grid nodes such that each node has 1-4 links
# and has a path to every other node.
#   -- @relsqui, https://twitter.com/relsqui/status/970555868390465536
#
# This counting approach implements a method similar to that described in
# "Spanning Trees in Grid Graphs", Paul Raff, https://arxiv.org/pdf/0809.2551.pdf
#
# The key idea is to create a transition matrix counting the number of
# ways of going from one partition to another, where the partition
# specifies which points on the edge of the grid are connected to each other.
# Once we do that for all K-partitions, we can easily compute the
# number of "mazes" for all KxN grids.
#
# The resulting sequence for NxN grids
# 1, 5, 431, 555195, 10286937043, 2692324030864335
# is not in OEIS, but there was a reference to it in a Japanese set of slides
# referring to Tutte polynomials, which may provide a better way of computing
# the sequence.
# http://www-imai.is.s.u-tokyo.ac.jp/~imai/lecture/AdvancedAlgorithmsReport2_20140709.pdf
# https://en.wikipedia.org/wiki/Tutte_polynomial
# The official name for these "mazes" appears to be a "spanning subgraph",
# i.e., the subgraph spans all nodes but has the same number of components.

def partition( collection ):
    """Generator that yields all partitions of "collection" into sets,
    represented as lists or lists."""
    if len( collection ) == 1:
        yield [ collection ]
    else:
        first = collection[0]
        rest = collection[1:]
        for p in partition( rest ):
            # First element could be in any of the sets returned
            for i in xrange( len( p ) ):
                yield p[:i] + [ [first] + p[i] ] + p[i+1:]
            yield [[first]] + p

def partitions( k ):
    return list( sorted( p ) for p in partition( range( 1, k+1 ) ) )
    
def isConsecutiveRange( nums ):
    """Assume the numbers are in order?"""
    for (n, x) in enumerate( nums, nums[0] ):
        if n != x:
            return False
    return True
    
def startingVector( k, partitions ):
    """Count the number of ways to connect the k vertices in a 1xk grid
    for each of the numbered partitions.

    This is 1 if all the numbered verticies are adjacent, and zero otherwise.
    """
    def possibleInK( p ):
        for s in p:
            if not isConsecutiveRange( s ):
                return 0
        return 1
    
    return [ possibleInK( p ) for p in partitions ]

import itertools
import numpy as np

def allSubsets( collection ):
    for cardinality in xrange( len( collection ) + 1 ):
        for s in itertools.combinations( collection, cardinality ):
            yield s
                
def adjacencyMatrices( k, startPartition ):
    # We will use 1..k for the end and k+1..2k for the previous vertices
    a = np.zeros( ( 2 * k + 1, 2 * k + 1 ), np.int8 )
    # Fill in the startPosition adjacencies
    for p in startPartition:
        if len( p ) > 0:
            for (i,j) in zip( p, p[1:] ):
                a[k+i,k+j] = 1
                a[k+j,k+i] = 1
    # Now we'll iterate over all horizontal and vertical edges
    while True:
        yield a
        # Increment verticals
        carry = True
        for v in xrange( 1, k ):
            if a[v,v+1] == 0:
                a[v,v+1] = 1
                a[v+1,v] = 1
                carry = False
                break
            else:
                # Carry to next "digit"
                a[v,v+1] = 0
                a[v+1,v] = 0
        if carry:
            for h in xrange( 1, k+1 ):
                if a[h,k+h] == 0:
                    a[h,k+h] = 1
                    a[k+h,h] = 1
                    carry = False
                    break
                else:
                    a[h,k+h] = 0
                    a[k+h,h] = 0
        if carry:
            # All done
            break

import networkx as nx
import matplotlib.pyplot as plt
import random

def drawAdj( k, a ):
    graph = nx.convert_matrix.from_numpy_matrix( a )
    positions = [ (100,100) ] + \
        [ ( 200 + random.randint(-10,10), 100 * i ) for i in xrange( 1, k+1 ) ] + \
        [ ( 100 + random.randint(-10,10), 100 * i ) for i in xrange( 1, k+1 ) ]
    positions[0] = positions[1]    
    nx.draw( graph, positions, node_size=5 )
    print a
    for c in nx.connected_components( graph ):
        print c    
    plt.show()

def isGraphCompatible( k, a, partitions ):
    """Returns the index of the partition this graph is compatible with,
    or None if not compatible."""
    graph = nx.convert_matrix.from_numpy_matrix( a )
    # FIXME: this conversion is probably expensive, and
    # many graphs could be eliminated right away rather than continuing
    # to enumerate all their connected components.
    ccs = list( nx.connected_components( graph ) )
    #print "Matrix", a, "components", ccs
    
    # Each of the connected components we want has elements from 1..k, but
    # may also have larger items.  Any connected component with just k+1 and
    # higher means we have lost overall connectivity.
    part = []
    for c in ccs:
        c2 = filter( lambda x : x <= k,
                     iter( c ) )
        if len( c2 ) == 0:
            return None
        c3 = sorted( c2 )
        if c3 == [ 0 ]:
            continue
        part.append( c3 )

    part.sort()
    #print "part", part
    try:
        return partitions.index( part )
    except ValueError:
        return None
                
def buildCountingMatrix( k ):
    pp = partitions( k )
    numPart = len( pp )
    count = np.zeros( ( numPart, numPart ), np.int64 )
    for i in xrange( numPart ):
        for a in adjacencyMatrices( k, pp[i] ):
            j = isGraphCompatible( k, a, pp )
            if j != None:
                count[i,j] += 1
    return ( startingVector( k, pp ), count )

def nxnTable():
    for k in xrange( 1, 100 ):
        print k,
        (v, a) = buildCountingMatrix( k )
        ap = np.linalg.matrix_power( a, k - 1 )
        print np.matmul( v, ap )[0]