网格上的2D箱包装

Question

我有一个n × m网格和一组polyominos.我想知道是否可以将它们打包到网格中:不允许重叠或旋转.

我希望像大多数包装问题一样,这个版本是NP难以近似的,所以我不期待任何疯狂,但是算法可以在25×25左右的网格上找到合理的包装并且相当全面,大约10×10会很好.(我的瓷砖大多是tetrominos - 四个块 - 但它们可能有5-9 +块.)

我会接受任何人提供的任何东西:算法,论文,可以调整的现有程序.

Answer 1

这是一个类似原型的SAT求解器方法,它解决了:

• 先验固定的多米诺骨牌模式(见Constants / Input代码)
  • 如果应该允许旋转,则必须将旋转的碎片添加到该组中
• 每个多聚体都可以放置0-inf次
• 除此之外没有得分机制:
  • 非覆盖瓷砖的数量最小化!

考虑到组合优化(SAT,CP,MIP)的经典现成方法,这个方法可能会扩展得最好(有根据的猜测).在设计定制启发式时,它也很难被击败!

如果需要,这些幻灯片在实践中为SAT解算器提供了一些实用的介绍.在这里,我们使用完整的基于CDCL的求解器(如果有的话,将始终在有限时间内找到解决方案;如果没有,将始终能够证明在有限时间内没有解决方案;内存当然也起作用!).

一般而言,更复杂(线性)的每瓦评分函数难以合并.这是(M)IP方法可以更好的地方.但就纯粹搜索而言,SAT解决方案总体上要快得多.

N=25我的polyomino-set 的问题需要大约1秒(并且一个可以在多个粒度级别上轻松地将其平行化 - > SAT-solver(threadings-param)与外部循环;后者将在后面解释).

当然以下是:

• 因为这是一个NP难题,所以会有简单且不容易的实例
• 我没有用许多不同的polyominos进行科学基准测试
  • 可以预期,有些集合比其他集合更容易解决
• 这是一个可能的SAT公式(不是最微不足道的!)无数的
  • 每种配方都有优点和缺点

理念

一般方法是创建一个决策问题并将其转换为CNF,然后由高效的SAT解算器解决(此处:cryptominisat; CNF将采用DIMCAS-CNF格式),它将用作黑盒解​​算器(没有参数调整!).

由于目标是优化填充瓦片的数量并且我们正在使用决策问题,我们需要一个外环,添加一个最小的瓦片使用约束并尝试解决它.如果不成功,请减少此数字.所以一般来说,我们多次调用SAT求解器(从头开始!).

CNF有许多不同的配方/转化可能.这里我们使用(二元)决策变量X来指示一个位置.甲放置就像是一个元组polyomino, x_index, y_index(这索引标记一些图案的左上字段).变量数量和所有多项式的可能位置数量之间存在一对一的映射关系.

核心思想是:在一个解决方案的所有可能的放置组合的空间中搜索,这不会使一些约束无效.

此外,我们还有决策变量Y,表示正在填充的图块.有M*N这样的变数.

当有权访问所有可能的展示位置时,很容易为每个tile-index(M*N)计算碰撞集.给定一些固定的磁贴,我们可以检查哪些展示位置可以填充这个展示位置,并将问题限制为仅选择<=1那些展示位置.这是活跃的X.在(M)IP世界中,这可能被称为凸壳,用于碰撞.

n<=k - 约束在SAT求解中无处不在,许多不同的配方都是可能的.朴素编码通常需要指数数量的子句,这很容易变得不可行.使用新变量,可以使用许多变量子句权衡(参见Tseitin编码).我正在重复使用一个(旧代码;只有我的代码只有python2的原因),这对我来说过去很有用.它基于将基于硬件的反逻辑描述为CNF,并提供良好的经验和理论性能(见论文).当然还有很多选择.

另外,我们需要强制SAT求解器不要使所有变量都为负.我们必须添加描述以下内容的约束(这是一种方法):

• 如果使用某个字段:必须至少有一个有效的放置位置(poly + x + y),这会导致覆盖此字段!
  • 这是一个基本的逻辑含义,很容易被表述为一个潜在的大逻辑或

然后只丢失核心循环,尝试填充N个字段,然后填写N-1直到成功.这再次使用n<=k前面提到的配方.

码

这是python2代码,它需要运行脚本的目录中的SAT-solver cryptominisat 5.

我也在使用python优秀科学堆栈的工具.

# 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!

输出控制台

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]]

输出图

此参数化不能涵盖一个字段!

其他一些具有更大模式的示例

Square M=N=61(prime - > intuition:harder)其中base-CNF有450.723个子句和185.462个变量.有一个最佳的包装!

非正方形 M,N =83,131(双素数),其中基数 - CNF具有1.346.511子句和553.748个变量.有一个最佳的包装!