2D bin packing on a grid

Asked 21/12, 2017 at 5:58 Answered 22/12, 2017 at 1:10

Solved algorithm mathematical-optimization discrete-mathematics bin-packing

I have an n × m grid and a collection of polyominos. I would like to know if it is possible to pack them into the grid: no overlapping or rotation is allowed.

I expect that like most packing problems this version is NP-hard and difficult to approximate, so I'm not expecting anything crazy, but an algorithm that could find reasonable packings on a grid around 25 × 25 and be fairly comprehensive around 10 × 10 would be great. (My tiles are mostly tetrominos -- four blocks -- but they could have 5–9+ blocks.)

I'll take whatever anyone has to offer: an algorithm, a paper, an existing program which can be adapted.

Hippocrene answered 21/12, 2017 at 5:58 Comment(14)

Im assuming you're discounting brute force approaches, similar to how one solves a knights tour, in favor of a more heuristic/algorithmic process? – Philipp 21/12, 2017 at 7:1

Could you share some sample input data and expected output? – Jubbah 21/12, 2017 at 7:11

You should search on bin packing, there is alot of articles on this. Maybe This is handy https://mcmap.net/q/944055/-block-layout-algorithm? – Roadhouse 21/12, 2017 at 12:19

Can you repeat polyominoes, or can you use each at most once? Is there any value to maximize (as in: each poliominoe may have a different value), or is this just a boolean "can it be done or not" question? – Ignatius 21/12, 2017 at 12:48

It would also be nice to know, how many polyominos you got. – Puga 21/12, 2017 at 14:52

@Ignatius I would be interested in any version you have information on. The most natural version for me is one in which I have infinitely many copies of each polyomino and an associated value with each, nonincreasing in each polyomino. The simplification where each shape has the same value is a good approximation. I've been working with the question "is a given arrangement possible" because it seems simpler, but if you can attack the harder optimization problem by all means go ahead. – Hippocrene 21/12, 2017 at 15:9

@Puga The packings are pretty tight, with just a few empty blocks. In a 12×12 I might have 32 4-ominos and 2 5-ominos, for example, leaving 6 blocks empty across the board. – Hippocrene 21/12, 2017 at 15:13

Well the candidates are the usual discrete-opt tools like SAT, CP, MIP, (Metaheuristics). But without a clear mathematical-specification it feels way too broad (go tackle whatever you want ... i don't like it ... imagine the spread of potential answers; i might also note, that the problem changed a lot from the original question to the last comment, which can be annoying too if someone already started working given incomplete assumptions). – Puga 21/12, 2017 at 15:13

@Puga I asked the narrower question (given a multiset of polyominos, can they fit on the grid?) because I didn't think the question I wanted could reasonably be addressed. But tucuxi asked for more details so I gave the actual underlying problem of interest. I don't want to let that distract from the question I asked though. I think your last comment is getting close to an answer though. – Hippocrene 21/12, 2017 at 15:19

@GijsDenHollander I'm familiar with bin packing, and I've even found some tools which could be bludgeoned for this purpose, but they're not a very good fit for what I'm doing. I skimmed a few academic papers on 2D bin packing but most focus on the special case of rectangles and my polyominos are more general. So I thought I'd tap the community wisdom here and see what was known rather than continue to flounder around. – Hippocrene 21/12, 2017 at 15:23

@JakeHeidt I don't think raw brute force would work because there's a combinatorial explosion in terms of number of pieces and placements. But I'd be happy for an approach that was a clever application of lots of cycles; I could see spending 10^13 or so on a given instance. that would be fine. But heuristics seem like a likely approach given what I know. – Hippocrene 21/12, 2017 at 15:26

I recommend SAT (for non-weight/score based feasibility; will be pretty pretty hard to beat imho) first. MIP for the more general problem. Don't start with heuristics which probably need weeks of tuning. At least that's my opinion. – Puga 21/12, 2017 at 15:28

Charles, it would be pretty helpful if you could specify a set of tetrominoes and a board size that you think would be interesting. Then we can test solutions and determine if we can solve in reasonable runtimes for the cases that you find interesting. – Jubbah 21/12, 2017 at 18:33

I've added a follow-up question on whether it's possible to better in the special case where all polyominos are in fact rectangles. Perhaps introduce non-binary integer variables for the location of each polyomino, then add constraints representing non-overlap for each of them? – Piccalilli 17/4, 2022 at 11:15

Here is a prototype-like SAT-solver approach, which tackles:

a-priori fixed polyomino patterns (see Constants / Input in code)
- if rotations should be allowed, rotated pieces have to be added to the set
every polyomino can be placed 0-inf times
there is no scoring-mechanic besides:
- the number of non-covered tiles is minimized!

Considering classic off-the-shelf methods for combinatorial-optimization (SAT, CP, MIP), this one will probably scale best (educated guess). It will also be very hard to beat when designing customized heuristics!

If needed, these slides provide some practical introduction to SAT-solvers in practice. Here we are using CDCL-based solvers which are complete (will always find a solution in finite time if there is one; will always be able to prove there is no solution in finite time if there is none; memory of course also plays a role!).

More complex (linear) per-tile scoring-functions are hard to incorporate in general. This is where a (M)IP-approach can be better. But in terms of pure search SAT-solving is much faster in general.

The N=25 problem with my polyomino-set takes ~ 1 second (and one could easily parallize this on multiple granularity-levels -> SAT-solver (threadings-param) vs. outer-loop; the latter will be explained later).

Of course the following holds:

as this is an NP-hard problem, there will be easy and non-easy instances
i did not do scientific benchmarks with many different sets of polyominos
- it's to be expected that some sets are easier to solve than others
this is one possible SAT-formulation (not the most trivial!) of infinite many
- each formulation has advantages and disadvantages

Idea

The general approach is creating a decision-problem and transforming it into CNF, which is then solved by highly efficient SAT-solvers (here: cryptominisat; CNF will be in DIMCAS-CNF format), which will be used as black-box solvers (no parameter-tuning!).

As the goal is to optimize the number of filled tiles and we are using a decision-problem, we need an outer-loop, adding a minimum tile-used constraint and try to solve it. If not successful, decrease this number. So in general we are calling the SAT-solver multiple times (from scratch!).

There are many different formulations / transformations to CNF possible. Here we use (binary) decision-variables X which indicate a placement. A placement is a tuple like polyomino, x_index, y_index (this index marks the top-left field of some pattern). There is a one-to-one mapping between the number of variables and the number of possible placements of all polyominos.

The core idea is: search in the space of all possible placement-combinations for one solution, which is not invalidating some constraints.

Additionally, we have decision-variables Y, which indicate a tile being filled. There are M*N such variables.

When having access to all possible placements, it's easy to calculate a collision-set for each tile-index (M*N). Given some fixed tile, we can check which placements can fill this one and constrain the problem to only select <=1 of those. This is active on X. In the (M)IP world this probably would be called convex-hull for the collisions.

n<=k-constraints are ubiquitous in SAT-solving and many different formulations are possible. Naive-encoding would need an exponential number of clauses in general which easily becomes infeasibly. Using new variables, there are many variable-clause trade-offs (see Tseitin-encoding) possible. I'm reusing one (old code; only reason why my code is python2-only) which worked good for me in the past. It's based on describing hardware-based counter-logic into CNF and provides good empirical- and theoretical performance (see paper). Of course there are many alternatives.

Additionally, we need to force the SAT-solver not to make all variables negative. We have to add constraints describing the following (that's one approach):

if some field is used: there has to be at least one placement active (poly + x + y), which results in covering this field!
- this is a basic logical implication easily formulated as one potentially big logical or

Then only the core-loop is missing, trying to fill N fields, then N-1 until successful. This is again using the n<=k formulation mentioned earlier.

Code

This is python2-code, which needs the SAT-solver cryptominisat 5 in the directory the script is run from.

I'm also using tools from python's excellent scientific-stack.

# PYTHON 2!
import math
import copy
import subprocess
import numpy as np
import matplotlib.pyplot as plt      # plotting-only
import seaborn as sns                # plotting-only
np.set_printoptions(linewidth=120)   # more nice console-output

""" Constants / Input
        Example: 5 tetrominoes; no rotation """
M, N = 25, 25
polyominos = [np.array([[1,1,1,1]]),
              np.array([[1,1],[1,1]]),
              np.array([[1,0],[1,0], [1,1]]),
              np.array([[1,0],[1,1],[0,1]]),
              np.array([[1,1,1],[0,1,0]])]

""" Preprocessing
        Calculate:
        A: possible placements
        B: covered positions
        C: collisions between placements
"""
placements = []
covered = []
for p_ind, p in enumerate(polyominos):
    mP, nP = p.shape
    for x in range(M):
        for y in range(N):
            if x + mP <= M:          # assumption: no zero rows / cols in each p
                if y + nP <= N:      # could be more efficient
                    placements.append((p_ind, x, y))
                    cover = np.zeros((M,N), dtype=bool)
                    cover[x:x+mP, y:y+nP] = p
                    covered.append(cover)                           
covered = np.array(covered)

collisions = []
for m in range(M):
    for n in range(N):
        collision_set = np.flatnonzero(covered[:, m, n])
        collisions.append(collision_set)

""" Helper-function: Cardinality constraints """
# K-ARY CONSTRAINT GENERATION
# ###########################
# SINZ, Carsten. Towards an optimal CNF encoding of boolean cardinality constraints.
# CP, 2005, 3709. Jg., S. 827-831.

def next_var_index(start):
    next_var = start
    while(True):
        yield next_var
        next_var += 1

class s_index():
    def __init__(self, start_index):
        self.firstEnvVar = start_index

    def next(self,i,j,k):
        return self.firstEnvVar + i*k +j

def gen_seq_circuit(k, input_indices, next_var_index_gen):
    cnf_string = ''
    s_index_gen = s_index(next_var_index_gen.next())

    # write clauses of first partial sum (i.e. i=0)
    cnf_string += (str(-input_indices[0]) + ' ' + str(s_index_gen.next(0,0,k)) + ' 0\n')
    for i in range(1, k):
        cnf_string += (str(-s_index_gen.next(0, i, k)) + ' 0\n')

    # write clauses for general case (i.e. 0 < i < n-1)
    for i in range(1, len(input_indices)-1):
        cnf_string += (str(-input_indices[i]) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
        cnf_string += (str(-s_index_gen.next(i-1, 0, k)) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
        for u in range(1, k):
            cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, u-1, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
            cnf_string += (str(-s_index_gen.next(i-1, u, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
        cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, k-1, k)) + ' 0\n')

    # last clause for last variable
    cnf_string += (str(-input_indices[-1]) + ' ' + str(-s_index_gen.next(len(input_indices)-2, k-1, k)) + ' 0\n')

    return (cnf_string, (len(input_indices)-1)*k, 2*len(input_indices)*k + len(input_indices) - 3*k - 1)

def gen_at_most_n_constraints(vars, start_var, n):
    constraint_string = ''
    used_clauses = 0
    used_vars = 0
    index_gen = next_var_index(start_var)
    circuit = gen_seq_circuit(n, vars, index_gen)
    constraint_string += circuit[0]
    used_clauses += circuit[2]
    used_vars += circuit[1]
    start_var += circuit[1]

    return [constraint_string, used_clauses, used_vars, start_var]

def parse_solution(output):
    # assumes there is one
    vars = []
    for line in output.split("\n"):
        if line:
            if line[0] == 'v':
                line_vars = list(map(lambda x: int(x), line.split()[1:]))
                vars.extend(line_vars)
    return vars

def solve(CNF):
    p = subprocess.Popen(["cryptominisat5.exe"], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
    result = p.communicate(input=CNF)[0]
    sat_line = result.find('s SATISFIABLE')
    if sat_line != -1:
        # solution found!
        vars = parse_solution(result)
        return True, vars
    else:
        return False, None

""" SAT-CNF: BASE """
X = np.arange(1, len(placements)+1)                                     # decision-vars
                                                                        # 1-index for CNF
Y = np.arange(len(placements)+1, len(placements)+1 + M*N).reshape(M,N)
next_var = len(placements)+1 + M*N                                      # aux-var gen
n_clauses = 0

cnf = ''                                                                # slow string appends
                                                                        # int-based would be better
# <= 1 for each collision-set
for cset in collisions:
    constraint_string, used_clauses, used_vars, next_var = \
        gen_at_most_n_constraints(X[cset].tolist(), next_var, 1)
    n_clauses += used_clauses
    cnf += constraint_string

# if field marked: one of covering placements active
for x in range(M):
    for y in range(N):
        covering_placements = X[np.flatnonzero(covered[:, x, y])]  # could reuse collisions
        clause = str(-Y[x,y])
        for i in covering_placements:
            clause += ' ' + str(i)
        clause += ' 0\n'
        cnf += clause
        n_clauses += 1

print('BASE CNF size')
print('clauses: ', n_clauses)
print('vars: ', next_var - 1)

""" SOLVE in loop -> decrease number of placed-fields until SAT """
print('CORE LOOP')
N_FIELD_HIT = M*N
while True:
    print(' N_FIELDS >= ', N_FIELD_HIT)
    # sum(y) >= N_FIELD_HIT
    # == sum(not y) <= M*N - N_FIELD_HIT
    cnf_final = copy.copy(cnf)
    n_clauses_final = n_clauses

    if N_FIELD_HIT == M*N:  # awkward special case
        constraint_string = ''.join([str(y) + ' 0\n' for y in Y.ravel()])
        n_clauses_final += N_FIELD_HIT
    else:
        constraint_string, used_clauses, used_vars, next_var = \
            gen_at_most_n_constraints((-Y).ravel().tolist(), next_var, M*N - N_FIELD_HIT)
        n_clauses_final += used_clauses

    n_vars_final = next_var - 1
    cnf_final += constraint_string
    cnf_final = 'p cnf ' + str(n_vars_final) + ' ' + str(n_clauses) + \
        ' \n' + cnf_final  # header

    status, sol = solve(cnf_final)
    if status:
        print(' SOL found: ', N_FIELD_HIT)

        """ Print sol """
        res = np.zeros((M, N), dtype=int)
        counter = 1
        for v in sol[:X.shape[0]]:
            if v>0:
                p, x, y = placements[v-1]
                pM, pN = polyominos[p].shape
                poly_nnz = np.where(polyominos[p] != 0)
                x_inds, y_inds = x+poly_nnz[0], y+poly_nnz[1]
                res[x_inds, y_inds] = p+1
                counter += 1
        print(res)

        """ Plot """
        # very very ugly code; too lazy
        ax1 = plt.subplot2grid((5, 12), (0, 0), colspan=11, rowspan=5)
        ax_p0 = plt.subplot2grid((5, 12), (0, 11))
        ax_p1 = plt.subplot2grid((5, 12), (1, 11))
        ax_p2 = plt.subplot2grid((5, 12), (2, 11))
        ax_p3 = plt.subplot2grid((5, 12), (3, 11))
        ax_p4 = plt.subplot2grid((5, 12), (4, 11))
        ax_p0.imshow(polyominos[0] * 1, vmin=0, vmax=5)
        ax_p1.imshow(polyominos[1] * 2, vmin=0, vmax=5)
        ax_p2.imshow(polyominos[2] * 3, vmin=0, vmax=5)
        ax_p3.imshow(polyominos[3] * 4, vmin=0, vmax=5)
        ax_p4.imshow(polyominos[4] * 5, vmin=0, vmax=5)
        ax_p0.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p1.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p2.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p3.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p4.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p0.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p1.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p2.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p3.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p4.yaxis.set_major_formatter(plt.NullFormatter())

        mask = (res==0)
        sns.heatmap(res, cmap='viridis', mask=mask, cbar=False, square=True, linewidths=.1, ax=ax1)
        plt.tight_layout()
        plt.show()
        break

    N_FIELD_HIT -= 1  # binary-search could be viable in some cases
                      # but beware the empirical asymmetry in SAT-solvers:
                      #    finding solution vs. proving there is none!

Output console

BASE CNF size
('clauses: ', 31509)
('vars: ', 13910)
CORE LOOP
(' N_FIELDS >= ', 625)
(' N_FIELDS >= ', 624)
(' SOL found: ', 624)
[[3 2 2 2 2 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 2 2]
 [3 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 2 2]
 [3 3 3 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 1 1 1 1 2 2]
 [2 2 3 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [2 2 3 3 3 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [1 1 1 1 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2]
 [1 1 1 1 3 3 3 2 2 1 1 1 1 2 2 2 2 2 2 2 2 1 1 1 1]
 [2 2 1 1 1 1 3 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [2 2 2 2 2 2 3 3 3 2 2 2 2 1 1 1 1 2 2 2 2 2 2 2 2]
 [2 2 2 2 2 2 2 2 3 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2]
 [2 2 1 1 1 1 2 2 3 3 3 2 2 2 2 2 2 1 1 1 1 2 2 2 2]
 [1 1 1 1 1 1 1 1 2 2 3 2 2 1 1 1 1 1 1 1 1 1 1 1 1]
 [2 2 3 1 1 1 1 3 2 2 3 3 4 1 1 1 1 2 2 1 1 1 1 2 2]
 [2 2 3 1 1 1 1 3 1 1 1 1 4 4 3 2 2 2 2 1 1 1 1 2 2]
 [2 2 3 3 5 5 5 3 3 1 1 1 1 4 3 2 2 1 1 1 1 1 1 1 1]
 [2 2 2 2 4 5 1 1 1 1 1 1 1 1 3 3 3 2 2 1 1 1 1 2 2]
 [2 2 2 2 4 4 2 2 1 1 1 1 1 1 1 1 3 2 2 1 1 1 1 2 2]
 [2 2 2 2 3 4 2 2 2 2 2 2 1 1 1 1 3 3 3 2 2 2 2 2 2]
 [3 4 2 2 3 5 5 5 2 2 2 2 1 1 1 1 2 2 3 2 2 2 2 2 2]
 [3 4 4 3 3 3 5 5 5 5 1 1 1 1 2 2 2 2 3 3 3 2 2 2 2]
 [3 3 4 3 1 1 1 1 5 1 1 1 1 4 2 2 2 2 2 2 3 2 2 2 2]
 [2 2 3 3 3 1 1 1 1 1 1 1 1 4 4 4 2 2 2 2 3 3 0 2 2]
 [2 2 3 1 1 1 1 1 1 1 1 5 5 5 4 4 4 1 1 1 1 2 2 2 2]
 [2 2 3 3 1 1 1 1 1 1 1 1 5 5 5 5 4 1 1 1 1 2 2 2 2]
 [2 2 1 1 1 1 1 1 1 1 1 1 1 1 5 1 1 1 1 1 1 1 1 2 2]]

Output plot

One field cannot be covered in this parameterization!

Some other examples with a bigger set of patterns

Square M=N=61 (prime -> intuition: harder) where the base-CNF has 450.723 clauses and 185.462 variables. There is an optimal packing!

Non-square M,N =83,131 (double prime) where the base-CNF has 1.346.511 clauses and 553.748 variables. There is an optimal packing!

Puga answered 22/12, 2017 at 1:10 Comment(6)

This is quite impressive, since you seem so knowledgeable about using a SAT solver, could you give me a hand in recommending how to learn the basics of them and when they should be used. I always have a hard time getting started with algo systems like a sat solver, i dont know where to even begin to wrap my head around it, and i give up on it. – Philipp 22/12, 2017 at 13:32

@JakeHeidt In general i would say: SAT is the most powerful on purely combinatoric tasks (not much arithmetic). It's hard to draw line between SAT/CP although those are very different. If one does not have good use of some powerful global-constraint i don't think CP will be competetive. On the plus side: CP allows easier modelling / more lang-stuctures like all-diff. (M)IP only feels right if there are enough linear-terms to make SAT-struggle while IP can take off. Pseudo-boolean tasks are often well-suited for MIP. – Puga 22/12, 2017 at 13:53

When one decided to use any of those 3 techniques, it's all about modelling. There are trade-offs and modelling is not perfectly understood in all of those domains. It's often based on intuition and empirical-tests (the same i would say for deep learning which has even less theory). For SAT for example: redundant constraint often help, although it feels like wasted clauses.In SAT-world, k-ary constraints are probably the most important and there are survey-papers with many approaches. But of course some models need k-ary, some do not.Besides that, solvers more or less are used as black-box – Puga 22/12, 2017 at 13:55

Theoretically in SAT/CP propagation efficiency is analyzed. Simple example: the paper of my k-ary constraints (Sinz): there is a nk* sequential-counter-approach and a n parallel-counter approach. Only the former, while being bigger can use unit-propagation (most important component of SAT-solver) efficiently. In IP it's probably easier to understand some models for experts (integrality-gap), but still hard. In IP there is one more thing entirely different: the best tools are commercial and closed-source. – Puga 22/12, 2017 at 14:0

I did use SAT-solvers for a lot of private tasks, often with success. But i'm surely not an expert. When people with a more deep understand tackle problems, things can get interesting. I always had a hard time to follow this one for example. I formulated the same problem before, many different encodings, but it was completely different. – Puga 22/12, 2017 at 14:6

Thank you so much for such a detailed reply! As far as not being an expert, posting answers to SO questions utilizing the technology in question sounds like an expert to me ;) Hopefully I'll be able to surmount my apparent inability to learn this stuff, with your help. – Philipp 22/12, 2017 at 14:44

One approach could be using integer programming. I'll implement this using the python pulp package, though packages are available for pretty much any programming language.

The basic idea is to define a decision variable for every possible placement location for every tile. If a decision variable takes value 1, then its associated tile is placed there. If it takes value 0, then it is not placed there. The objective is therefore to maximize the sum of the decision variables times the number of squares in the variable's tile --- this corresponds to placing the maximum number of squares possible on the board.

My code implements two constraints:

Each tile can only be placed once (below we will relax this constraint)
Each square can have at most one tile on it

Here's the output for a set of five fixed tetrominoes on a 4x5 grid:

import itertools
import pulp
import string

def covered(tile, base):
    return {(base[0] + t[0], base[1] + t[1]): True for t in tile}

tiles = [[(0,0), (1,0), (0,1), (0,2)],
         [(0,0), (1,0), (2,0), (3,0)],
         [(1,0), (0,1), (1,1), (2,0)],
         [(0,0), (1,0), (0,1), (1,1)],
         [(1,0), (0,1), (1,1), (2,1)]]
rows = 25
cols = 25
squares = {x: True for x in itertools.product(range(rows), range(cols))}
vars = list(itertools.product(range(rows), range(cols), range(len(tiles))))
vars = [x for x in vars if all([y in squares for y in covered(tiles[x[2]], (x[0], x[1])).keys()])]
x = pulp.LpVariable.dicts('tiles', vars, lowBound=0, upBound=1, cat=pulp.LpInteger)
mod = pulp.LpProblem('polyominoes', pulp.LpMaximize)
# Objective value is number of squares in tile
mod += sum([len(tiles[p[2]]) * x[p] for p in vars])
# Don't use any shape more than once
for tnum in range(len(tiles)):
    mod += sum([x[p] for p in vars if p[2] == tnum]) <= 1
# Each square can be covered by at most one shape
for s in squares:
    mod += sum([x[p] for p in vars if s in covered(tiles[p[2]], (p[0], p[1]))]) <= 1
# Solve and output
mod.solve()
out = [['-'] * cols for rep in range(rows)]
chars = string.ascii_uppercase + string.ascii_lowercase
numset = 0
for p in vars:
    if x[p].value() == 1.0:
        for off in tiles[p[2]]:
            out[p[0] + off[0]][p[1] + off[1]] = chars[numset]
        numset += 1
for row in out:
    print(''.join(row))

It obtains the following optimal solution:

AAAB-
A-BBC
DDBCC
DD--C

If we allow repeats (comment out the constraint limiting to one copy of each shape), then we can completely tile the grid:

ABCDD
ABCDD
ABCEE
ABCEE

It worked near-instantaneously for a 10x10 grid:

ABCCDDEEFF
ABCCDDEEFF
ABGHHIJJKK
ABGHHIJJKK
LLGMMINOPP
LLGMMINOPP
QQRRSTNOUV
QQRRSTNOUV
WWXXSTYYUV
WWXXSTYYUV

The code obtains an optimal solution for the 25x25 grid in 100 seconds of runtime, though unfortunately there aren't enough letter and numbers for my output code to print the solution.

Jubbah answered 21/12, 2017 at 16:48 Comment(2)

Ok. i cleaned up my comments and will play a bit with your code. It's rather surprising to me, that a default CBC will be efficient here (not even tuning cuts). – Puga 21/12, 2017 at 17:23

@Puga Yeah -- I suspect that the problem is just super easy (solveable by hand in the case of repeats allowed for a lot of grid shapes). I suspect if I had a more diverse set of tetrominoes and didn't allow repeats then things would get a lot more challenging for the solver. – Jubbah 21/12, 2017 at 17:30

I don't know if its of any use to you but I coded up a small sketchy frame in Python. It doesn't place polyminos yet but the functions are there - checking for dead empty spaces is primitive, though and needs a better approach. Then again, maybe it is all rubbish...

import functools
import itertools

M = 4 # x
N = 5 # y

field = [[9999]*(N+1)]+[[9999]+[0]*N+[9999] for _ in range(M)]+[[9999]*(N+1)]

def field_rd(p2d):
    return field[p2d[0]+1][p2d[1]+1]

def field_add(p2d,val):
    field[p2d[0]+1][p2d[1]+1] += val

def add2d(p,k):
    return p[0]+k[0],p[1]+k[1]

def norm(polymino_2d):
    x0,y0 = min(x for x,y in polymino_2d),min(y for x,y in polymino_2d)
    return tuple(sorted(map(lambda p: add2d(p,(-x0,-y0)), polymino_2d)))

def create_cutoff(occupied):
    """Receive a polymino and create the outer area of squares which could be cut off by a placement of this polymino"""
    cutoff = set(itertools.chain.from_iterable(map(lambda p: add2d(p,(x,y)),occupied) for (x,y) in [(-1,0),(1,0),(0,-1),(0,1)])) #(-1,-1),(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1)]))
    return tuple(cutoff.difference(occupied))

def is_occupied(p2d):
    return field_rd(p2d) == 0

def is_cutoff(p2d):
    return not is_occupied(p2d) and all(map(is_occupied,map(lambda p: add2d(p,p2d),[(-1,0),(1,0),(0,-1),(0,1)])))

def polym_colliding(p2d,occupied):
    return any(map(is_occupied,map(lambda p: add2d(p,p2d),occupied)))

def polym_cutoff(p2d,cutoff):
    return any(map(is_cutoff,map(lambda p: add2d(p,p2d),cutoff)))

def put(p2d,occupied,polym_nr):
    for p in occupied:
        field_add(add2d(p2d,p),polym_nr)

def remove(p2d,occupied,polym_nr):
    for p in polym:
        field_add(add2d(p2d,p),-polym_nr)

def place(p2d,polym_nr):
    """Try to place a polymino at point p2d. If it fits without cutting off unreachable single cells return True else False"""
    occupied = polym[polym_nr][0]
    if polym_colliding(p2d,occupied):
        return False
    put(p2d,occupied,polym_nr)
    cutoff = polym[polym_nr][1]
    if polym_cutoff(p2d,cutoff):
        remove(p2d,occupied,polym_nr)
        return False
    return True

def NxM_array(N,M):
    return [[0]*N for _ in range(M)]


def generate_all_polyminos(n):
    """Create all polyminos with size n"""
    def gen_recur(polymino,i,result):
        if i > 1:
            new_pts = set(itertools.starmap(add2d,itertools.product(polymino,[(-1,0),(1,0),(0,-1),(0,1)])))
            new_pts = new_pts.difference(polymino)
            for p in new_pts:
                gen_recur(polymino.union({p}),i-1,result)
        else:
            result.add(norm(polymino))
    #---------------------------------------
    all_polyminos = set()
    gen_recur({(0,0)},n,all_polyminos)
    return all_polyminos

print("All possible Tetris blocks (all orientations): ",generate_all_polyminos(4))

Phratry answered 21/12, 2017 at 16:15 Comment(0)

Hot tags

Godot Unity Godot Help Programming Godot 4.X GUI GDScript 3D 2D Physics CSharp Godot 3.X VR XR Projects C++