Nonogram puzzle solution space

Question

I've written some code to brute-force enumerate the solution space of nonogram puzzles up to a certain size. It does a good job up to the first ~4 billion puzzles, but I run out of disk space to store any bigger collection of the solution space set.

Because of some incredibly interesting patterns within the encoded set of solvable instances, I'd appreciate my code being reviewed for correctness. This is an NP-complete problem, but shows self-similarity and self-affinity in the solution space with the current implementation. Please, help me shoot down this idea as quickly as possible -- I'd like to regain a productive life :)

Here's how the code is supposed to work. It outputs an XYZ point cloud of the nonogram solution space given a puzzle's maximum width, great for visualization.

The code generates all the possible permutations of an arbitrarily sized boolean image. Each permutation is a puzzle solution to at least one, if not many, boolean image input pairs. The input pairs are "projections" of the solution from each lattice direction, while having the same axis-constrained continuous runs of set bits. The only difference between solution and input is very subtle: the unset padding between contiguous runs is flexible along an axis.

Here's an example pairing of inputs and solution. Note that the pictured top-right solution may not be unique for the given input images, and the given input images don't necessarily construct only that solution. It's merely an example of the "nonogram property."

You'll notice in the code I have a peculiar encoding of the inputs and solutions as integers, essentially a traversal of the image's cells converted to a bitstring. This encoding is chosen as a visualization convenience, and I've attained similar patterns with different traversal orders. The main goal with the encoding is to reduce dimensionality for plotting with a one-to-one correspondence of images to integers, avoiding the problem of colliding identifiers.

I'm including a montage of the first four iterations as subsets of the solution space shadow. The full shadow is essentially a look-up table for an NP oracle, so seeing patterns here could have remarkable consequences. Arranged from left to right are tables scaled to a common 512x512 resolution -- 4, 256, 262144, and 4294967296 puzzles accounted for, respectively. Each black pixel represents an input pair with no solution, white says a solution exists.

from sys import argv
from itertools import product, chain, groupby
from functools import partial
from multiprocessing import Pool

def indices(width):
    for a in range(width):
        for b in range(a+1):
            yield (b, a)
        for b in reversed(range(a)):
            yield (a, b)

def encode(matrix, width):
    return sum(1<<i for i, bit
               in enumerate(matrix[i][j] for i, j in indices(width))
               if bit)

def count_runs(row):
    return [sum(group) for key, group in groupby(row) if key]

def flex(solution, width):
    counts = list(map(count_runs, solution))
    for matrix in product((False, True), repeat=width**2):
        candidate = list(zip(*[iter(matrix)]*width))
        if list(map(count_runs, candidate)) == counts:
            yield candidate

def nonogram_solutions(solution, width):
    xy = solution
    yx = list(zip(*solution))
    enc_sol = encode(solution, width)
    return [(encode(xy, width), encode(yx, width), enc_sol)
            for xy, yx in product(flex(xy, width), flex(yx, width))]

def main(width):
    pool = Pool()
    sol_matrices = (list(zip(*[iter(matrix)]*width)) for matrix
                    in product((False, True), repeat=width**2))
    nonograms = partial(nonogram_solutions, width=width)
    solutions = pool.imap_unordered(nonograms, sol_matrices, 1)
    pool.close()
    for xy, yx, solution in chain.from_iterable(solutions):
        print(solution, xy, yx)
    pool.join()

if __name__ == "__main__":
    main(int(argv[1]))

Answer 1

1. There are no docstrings. What do these functions do? Lack of docstrings make the code hard to review, because we don't know what the functions are supposed to do.

2. The call indices(n) generates the Cartesian product of range(n) × range(n) in a particular order. A natural question is, does the order matter, or could we use itertools.product instead:

   itertools.product(range(width), repeat=2)
   

   The post explains why you've chosen the order. But how is someone reading the code supposed to know that? There needs to be a comment.

3. encode could be simplified to avoid the if:

   sum(bit<<i for i, bit in enumerate(...))
   

4. This code:

   for matrix in product((False, True), repeat=width**2):
       candidate = list(zip(*[iter(matrix)]*width))
   

   is effectively the same as:

   sol_matrices = (list(zip(*[iter(matrix)]*width)) for matrix
                   in product((False, True), repeat=width**2))
   

   and so would benefit from having its own function (which would have a name and docstring that would explain what it does).

5. The algorithm takes every nonogram puzzle, and then compares the run counts with the run counts for every nonogram puzzle of the same size. If width is \$ n \$, there are \$ 2^{n^2} \$ nonogram puzzles, and it takes \$ Ω(n^2) \$ to compute the run counts for a single puzzle, so the overall runtime is the ludicrous \$ Ω(n^24^{n^2}) \$.

   This could be brought down to \$ Ω(n^22^{n^2}) \$ if you prepared a dictionary mapping run counts to sets of encoded nonogram solutions, and down to \$ Ω(n2^{n^2}) \$ if you prepared the list of all rows, prepared a dictionary mapping each row to its run counts, and then generated the puzzles and their run counts simultaneously using the Cartesian product of the rows.