Generate all the free polyominoes who's width and height is no larger than 8

Question

I'm working on a puzzle game and I need all the free polyominoes that can fit into an 8 by 8 grid. My code uses Redelmeier's Alogorithm for enumerating them, checks them for rotations and reflections, and makes sure it doesn't go outside the 8 by 8 grid. This isn't a packing problem. Each one is checked one at a time against this criteria.

Polyominoes are categorized by their area (i.e. the number of cells they contain). For example, you're probably familiar with Tetris. All of the tetris pieces have an area of 4. The number of pieces grows exponentially with the area.

Free polyominoes consider rotations, translations, and reflections identical. So the domino with two cells side by side is considered the same as the one with the cells arranged vertically. This is in contrast to fixed polyominoes where rotations, translations, and reflections are considered different. For my purposes, rotating a puzzle doesn't change it so it's duplicate data I don't want to save.

I mentioned the number of pieces grows exponentially and when the area is 64, the number of fixed polyominoes is 4,496,670,726,609,716,846,990,603,851,802,851,046. Burnside's Lemma will reduce this by about 8. That's still a lot but that set also contains a lot of aren't within the bounds of my constraints. This would include a lot that are mostly tall and thin or too wide and narrow. There's no way the upper bound could exceed the number of ways to arrange black and white pixels in an 8x8 image which would be 2^64 or 18,446,744,073,709,551,616.

I believe my program is correct because it produces free polyominoes from oeis A000105 up until n = 9. When n = 9, there's a single piece from that set that is too wide (the 1 by 9) and it gets removed. So I'm fairly confident that it works. The only problem is that it's slow and will take several 100 years to finish. I'm not sure what I can do to speed this up.

Lastly, it's saving the output as list of bitboards to a text file which I've commented out for you. I have another program that reads them.

import itertools
from collections import defaultdict
from timeit import default_timer as timer

# I'm only going to n = 63 because n = 64 is trivial.
n = 63
print("Number of Polyomino pieces up to size", n)
temp_output = []

smallest_piece = [(0, 0)]
pieces_of_size = {1: [smallest_piece]}


def possible_expansions(piece):
    neighbor_offsets = [(-1, 0), (1, 0), (0, -1), (0, 1)]
    positions = set()

    for (x, y) in piece:
        for dx, dy in neighbor_offsets:
            positions.add((x + dx, y + dy))

    expansions = []
    piece_set = set(piece)

    for p in positions:
        if p not in piece_set:
            expansions.append(piece + [p])
    return expansions

def rotate_90_cw(piece):
    return [(y, -x) for (x, y) in piece]

def flip(piece):
    return [(-x, y) for (x, y) in piece]

def canonical(piece):
    min_x = min(x for (x, y) in piece)
    min_y = min(y for (x, y) in piece)
    res = sorted((x - min_x, y - min_y) for (x, y) in piece)
    return res

def hash_piece(piece):
    return hash(tuple(piece))

def expand_pieces(pieces):
    expanded = []
    temp_output = []
    expanded_hashes = defaultdict(list)
    for piece in pieces:
        for e in possible_expansions(piece):
            exp = canonical(e)
            if (max(x for (x, y) in exp) >= 8) or (max(y for (x, y) in exp) >= 8):
                continue
            is_new = True
            if exp in expanded_hashes[hash_piece(exp)]:
                is_new = False
                continue
            for rotation in range(3):
                exp = canonical(rotate_90_cw(exp))
                if exp in expanded_hashes[hash_piece(exp)]:
                    is_new = False
                    continue
            exp = canonical(flip(exp))
            if exp in expanded_hashes[hash_piece(exp)]:
                is_new = False
                continue
            for rotation in range(3):
                exp = canonical(rotate_90_cw(exp))
                if exp in expanded_hashes[hash_piece(exp)]:
                    is_new = False
                    continue
            if is_new:
                temp_output.append(str(coordinates_list_to_bitboard(exp)) + "\n")
                expanded.append(exp)
                expanded_hashes[hash_piece(exp)].append(exp)
    # I've commented this out so it doesn't save anything to your hard drive.
    #output = "".join(temp_output)
    #f = open("C:\\Users\\rober\\Documents\\" + "size is " + str(len(exp)) +".txt", "w")
    #f.write(output)
    return expanded

def coordinates_list_to_bitboard(coordinates_list):
    bitboard = 0
    for x, y in coordinates_list:
        if 0 <= x < 8 and 0 <= y < 8:
            bit_position = y * 8 + x
            bitboard |= 1 << bit_position
        else:
            print(coordinates_list)
            raise ValueError(f"Coordinates ({x}, {y}) are out of bounds for an 8x8 board")
    return bitboard

for i in range(2, n + 1):
    start = timer()
    pieces_of_size[i] = expand_pieces(pieces_of_size[i - 1])
    print("Pieces with {} blocks: {}".format(i, len(pieces_of_size[i])))
    end = timer()
    print(end - start)

I think something about running it in Visual Studio Community is really slowing it down. I'm about 35x slower right out the gate compared with others.

Answer 1

Great first question!

I've been taking a crack at helping to reduce the runtime for this. So far, I seem to have been somewhat successful. (I'm averaging a ~21.5% increase in performance)

Here is what I suggest so far:

1.) Pre-sorting the expansions

By pre-sorting the expansions using a heap, we are making subsequent calls to canonical and it's use of sorted more efficient. This was the most noticeable performance increase. I chose a heap because it's a good option for sorting when appending to an empty list, which you are doing in possible_expansions with the expansions list. This is because heaps have an insertion time complexity of O(log n), which is very good.

Your code:

# Inside of `possible_expansions`
for p in positions:
    if p not in piece_set:
        expansions.append(piece + [p])

# Inside of `expand_pieces`
for piece in pieces:
    for e in possible_expansions(piece):
        exp = canonical(e)

My change:

# Inside of `possible_expansions`
for p in positions:
    if p not in piece_set:
        heappush(expansions, piece + [p])

# Inside of `expand_pieces`
exps = possible_expansions(piece)
for _ in range(len(exps)): # Use _ to show that we aren't using the loop var
    e = heappop(exps)
    exp = canonical(e)

2.) Reducing loop copies

Originally, canonical effectively re-uses the same list comprehension for getting min_x and min_y:

def canonical(piece):
    min_x = min(x for (x, y) in piece)
    min_y = min(y for (x, y) in piece)
    res = sorted((x - min_x, y - min_y) for (x, y) in piece)
    return res

By using a for loop instead of a list comprehension, we can more readably convert this into a single loop, reducing iterations by 50%:

def canonical(piece):
    min_x = float('inf')
    min_y = min_x
    for (x, y) in piece:
        if x < min_x:
            min_x = x
        if y < min_y:
            min_y = y
    res = sorted((x - min_x, y - min_y) for (x, y) in piece)
    return res

Since we got rid of the calls to min, we must keep track of the current minimums for x and y ourselves, and because this code could eventually reach extremely high numbers for x and y, we can set min_x and min_y to float('inf').

2.5) One more loop removal

On line 52 inside of expand_pieces you check if the current piece, exp, contains an x or y greater than or equal to8:

if (max(x for (x, y) in exp) >= 8) or (max(y for (x, y) in exp) >= 8):

We can remove one list comprehension, and replace max in the other with any, further reducing total iterations:

if any([(x > 7 or y > 7) for (x, y) in exp]):

I made another simplification by changing >= 8 to > 7, but this is more of a preference than anything.

We can replace the list comprehensions with a for loop, and break out as soon as we find a piece that is too big:

too_big = False
for (x, y) in exp:
    if (x > 7) or (y > 7) or (x < -7) or (y < -7):
        too_big = True
        break
if too_big:
    continue

2.75) Don't add bad permutations

In the first for loop of possible_expansions, we can reject pieces that are too large:

for (x, y) in piece:
    for dx, dy in neighbor_offsets:
        if abs(x + dx) > 7 or abs(y + dy) > 7:
            continue
        positions.add((x + dx, y + dy))

3.) Preventing re-initialization of constant variables

neighbor_offsets never changes but is initialized every time possible_expansions is called. You can move the declaration and initialization of neighbor_offsets to the global scope, just above the definition of possible_expansions.

Output sample comparison:

Number of Polyomino pieces up to size 63
    Output as-is                        Output after optimizations
Pieces with 10 blocks: 4645
0.36945559999730904         vs.     0.28241100000741426
Pieces with 11 blocks: 16929
1.572818599990569           vs.     1.208543800006737
Pieces with 12 blocks: 62149
6.77434189998894            vs.     5.203064699991955
Pieces with 13 blocks: 226701
28.675560999996378          vs.     22.347502399992663
Pieces with 14 blocks: 819508
116.19935209999676          vs.     91.290519400005

I will edit this if I find anything else!

Answer 2

There are 6,200,686,124,225,191 free polyominoes that touch all four edges of an 8x8 grid (A268311), plus however many are 8x7, 7x7, 8x6, etc.

That's too many for my disk. How about sampling some at random? Generate a random 64-bit integer and check that it represents a connected polyomino and that it is the minimal representation after translation, reflection and rotation. About 1 in 3000 integers should pass the tests, so you should have a large sample reasonably quickly.

Answer 3

I got nerdsniped in Discord with this one. No code yet, too distracted (/maybe later/someone else will need to translate to python), and the text is more suitable to a random post than a formal review. At a glance (and wholly ignoring the idea of saving anything to disk):

• At the surface-level, there is a lot of unnecessary generation of data. possible_expansions() generates 4 tuples to check for every pixel present in the piece. That gets collapsed by a set but the set check would be unnecessary with a bit of extra checking up front.
• Data structures are re-generated when needed instead of being cached somehow. Memory usage could balloon quickly if sets are stored with each piece though, so maybe something else (see below) would have to be used.
• I don't know python well enough to know how the in lookup works for sets (or list-as-sets?), but I'd imagine it has to be generic enough for the lack of types so it's probably a linear lookup as opposed to a hash or tree, which may in some cases be faster, but for this problem a custom solution would be better, especially if you stopped using tuples altogether.
• Keeping track of the bounds of pieces and maybe using multiple lookup tables for each boundary area could reduce the number of lookups and normalizations needed to check for existing pieces. Knowing bounds could also reduce the scope further by only allowing [square + wide] or [square + tall] bounds, as a piece in either set has an equivalent transformed piece in the other set. Normalizing bounds would only be needed when expanding left/up. Expansions could be reduced by only expanding the set of w*h pieces to a set of (w+1)*h or w*(h+1) pieces.
• The specific problem is to generate all possible polyominoes that fit within an 8*8 area. That means that up to n=64 polyominoes can be generated, but generating n=64 polyominoes in arbitrary space is massively inefficient as the problem describes. An 8x8 area has 64 possible [pixel] positions, all of those positions are constant, it should use a bit set instead of a set of tuples.
• Once you're in the land of bit sets, you can make optimizations on hash or tree lookups for pieces that already exist as difference checks become a single 64-bit integer comparison instead of an O(n^2) loop through all positions in two sets ("does any piece match this" as opposed to "does any piece have all positions matching all of the new positions").
• Bit sets make it easy to do transformations on the bits, you could use rotations on the i64 to shift individual bytes to imitate a vertical flip, and horizontal flips can be done easily enough with a [256] lookup table.
• Similarly, positions to expand to or functions used to expand can be placed in a lookup table for each direction, with the index taken by masking bytes for the horizontals and adding a few or-shifts for the verticals. Could maybe even do some clever hackery with the lookups by including two rows in the mask/shift and only doing vertical expansions where a horizontal would not have happened.