使用 SAT 求解器 (Python) 查找特定区域内的所有自由多联骨牌组合
Finding all the combinations of free polyominoes within a specific area with a SAT-solver (Python)
我是 SAT 求解器领域的新手，需要有关以下问题的一些指导。
考虑到：
❶ 我在一个 4*4 的网格中选择了 14 个相邻的单元格
❷ 我有 5 polyominoes (A, B, C, D, E) 尺码 4, 2, 5, 2 和 1
❸ 这些多联骨牌是自由的，即它们的形状不固定，可以形成不同的图案
我如何使用 SAT 求解器计算所选区域（灰色单元格）内这 5 个自由多联骨牌的所有可能组合？
借用@spinkus 富有洞察力的回答和 OR-tools 文档，我可以制作以下示例代码（在 Jupyter Notebook 中运行）：
from ortools.sat.python import cp_model
import numpy as np
import more_itertools as mit
import matplotlib.pyplot as plt
%matplotlib inline
W, H = 4, 4 #Dimensions of grid
sizes = (4, 2, 5, 2, 1) #Size of each polyomino
labels = np.arange(len(sizes)) #Label of each polyomino
colors = ('#FA5454', '#21D3B6', '#3384FA', '#FFD256', '#62ECFA')
cdict = dict(zip(labels, colors)) #Color dictionary for plotting
inactiveCells = (0, 1) #Indices of disabled cells (in 1D)
activeCells = set(np.arange(W*H)).difference(inactiveCells) #Cells where polyominoes can be fitted
ranges = [(next(g), list(g)[-1]) for g in mit.consecutive_groups(activeCells)] #All intervals in the stack of active cells
def main():
model = cp_model.CpModel()
#Create an Int var for each cell of each polyomino constrained to be within Width and Height of grid.
pminos = [[] for s in sizes]
for idx, s in enumerate(sizes):
for i in range(s):
pminos[idx].append([model.NewIntVar(0, W-1, 'p%i'%idx + 'c%i'%i + 'x'), model.NewIntVar(0, H-1, 'p%i'%idx + 'c%i'%i + 'y')])
#Define the shapes by constraining the cells relative to each other
## 1st polyomino -> tetromino ##
# #
# #
# # #
# ### #
# #
################################
p0 = pminos[0]
model.Add(p0[1][0] == p0[0][0] + 1) #'x' of 2nd cell == 'x' of 1st cell + 1
model.Add(p0[2][0] == p0[1][0] + 1) #'x' of 3rd cell == 'x' of 2nd cell + 1
model.Add(p0[3][0] == p0[0][0] + 1) #'x' of 4th cell == 'x' of 1st cell + 1
model.Add(p0[1][1] == p0[0][1]) #'y' of 2nd cell = 'y' of 1st cell
model.Add(p0[2][1] == p0[1][1]) #'y' of 3rd cell = 'y' of 2nd cell
model.Add(p0[3][1] == p0[1][1] - 1) #'y' of 3rd cell = 'y' of 2nd cell - 1
## 2nd polyomino -> domino ##
# #
# #
# # #
# # #
# #
#############################
p1 = pminos[1]
model.Add(p1[1][0] == p1[0][0])
model.Add(p1[1][1] == p1[0][1] + 1)
## 3rd polyomino -> pentomino ##
# #
# ## #
# ## #
# # #
# #
################################
p2 = pminos[2]
model.Add(p2[1][0] == p2[0][0] + 1)
model.Add(p2[2][0] == p2[0][0])
model.Add(p2[3][0] == p2[0][0] + 1)
model.Add(p2[4][0] == p2[0][0])
model.Add(p2[1][1] == p2[0][1])
model.Add(p2[2][1] == p2[0][1] + 1)
model.Add(p2[3][1] == p2[0][1] + 1)
model.Add(p2[4][1] == p2[0][1] + 2)
## 4th polyomino -> domino ##
# #
# #
# # #
# # #
# #
#############################
p3 = pminos[3]
model.Add(p3[1][0] == p3[0][0])
model.Add(p3[1][1] == p3[0][1] + 1)
## 5th polyomino -> monomino ##
# #
# #
# # #
# #
# #
###############################
#No constraints because 1 cell only
#No blocks can overlap:
block_addresses = []
n = 0
for p in pminos:
for c in p:
n += 1
block_address = model.NewIntVarFromDomain(cp_model.Domain.FromIntervals(ranges),'%i' % n)
model.Add(c[0] + c[1] * W == block_address)
block_addresses.append(block_address)
model.AddAllDifferent(block_addresses)
#Solve and print solutions as we find them
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(pminos)
status = solver.SearchForAllSolutions(model, solution_printer)
print('Status = %s' % solver.StatusName(status))
print('Number of solutions found: %i' % solution_printer.count)
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
''' Print a solution. '''
def __init__(self, variables):
cp_model.CpSolverSolutionCallback.__init__(self)
self.variables = variables
self.count = 0
def on_solution_callback(self):
self.count += 1
plt.figure(figsize = (2, 2))
plt.grid(True)
plt.axis([0,W,H,0])
plt.yticks(np.arange(0, H, 1.0))
plt.xticks(np.arange(0, W, 1.0))
for i, p in enumerate(self.variables):
for c in p:
x = self.Value(c[0])
y = self.Value(c[1])
rect = plt.Rectangle((x, y), 1, 1, fc = cdict[i])
plt.gca().add_patch(rect)
for i in inactiveCells:
x = i%W
y = i//W
rect = plt.Rectangle((x, y), 1, 1, fc = 'None', hatch = '///')
plt.gca().add_patch(rect)
问题是我硬编码了 5 个 unique/fixed 多联骨牌，我不知道 如何定义约束，以便考虑到每个多联骨牌的每个可能模式（前提是可以）。
对于每个 polyonomino 和每个可能的左上角单元格，您都有一个布尔变量，用于指示此单元格是否是封闭矩形的左上角部分。
对于每个单元格和每个 polyomino，您都有一个布尔变量，用于指示该单元格是否被该 polyomino 占用。
现在，对于每个单元格和每个 polyomino，你有一系列的含义：选择左上角的单元格意味着每个单元格实际上都被这个 polyomino 占用。
然后是约束条件：
对于每个单元格，最多有一个多联骨牌占据它
对于每个 polyomino，它的左上角恰好有一个单元格。
这是一道纯布尔题。
编辑： 我在原始答案中漏掉了 "free" 这个词，并使用 OR-Tools 给出了固定多联骨牌的答案。添加了一个部分来回答以包含免费多联骨牌的解决方案 - 结果证明 AFAICT 很难在使用 OR-Tools 的约束编程中精确表达。
使用 OR-TOOLS 修复的多联骨牌：
是的，你可以用 constraint programming in OR-Tools. OR-Tools knows nothing about 2D grid geometry so you have to encode the geometry of each shape you have in terms of positional constraints. I.e. a shape is a collection of blocks / cells that must have a certain relationship to each other, must be within the bounds of the grid and must not overlap. Once you have your constraint model you just ask the CP-SAT Solver 解决它，在你的情况下，所有可能的解决方案。
这是一个非常简单的概念证明，在 4x4 网格上有两个矩形（您可能还想添加某种解释器代码以从形状描述到一组 OR-Tools 变量和约束在更大规模的问题中，因为手动输入约束有点乏味）。
from ortools.sat.python import cp_model
(W, H) = (3, 3) # Width and height of our grid.
(X, Y) = (0, 1) # Convenience constants.
def main():
model = cp_model.CpModel()
# Create an Int var for each block of each shape constrained to be within width and height of grid.
shapes = [
[
[ model.NewIntVar(0, W, 's1b1_x'), model.NewIntVar(0, H, 's1b1_y') ],
[ model.NewIntVar(0, W, 's1b2_x'), model.NewIntVar(0, H, 's1b2_y') ],
[ model.NewIntVar(0, W, 's1b3_x'), model.NewIntVar(0, H, 's1b3_y') ],
],
[
[ model.NewIntVar(0, W, 's2b1_x'), model.NewIntVar(0, H, 's2b1_y') ],
[ model.NewIntVar(0, W, 's2b2_x'), model.NewIntVar(0, H, 's2b2_y') ],
]
]
# Define the shapes by constraining the blocks relative to each other.
# 3x1 rectangle:
s0 = shapes[0]
model.Add(s0[0][Y] == s0[1][Y])
model.Add(s0[0][Y] == s0[2][Y])
model.Add(s0[0][X] == s0[1][X] - 1)
model.Add(s0[0][X] == s0[2][X] - 2)
# 1x2 rectangle:
s1 = shapes[1]
model.Add(s1[0][X] == s1[1][X])
model.Add(s1[0][Y] == s1[1][Y] - 1)
# No blocks can overlap:
block_addresses = []
for i, block in enumerate(blocks(shapes)):
block_address = model.NewIntVar(0, (W+1)*(H+1), 'b%d' % (i,))
model.Add(block[X] + (H+1)*block[Y] == block_address)
block_addresses.append(block_address)
model.AddAllDifferent(block_addresses)
# Solve and print solutions as we find them
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(shapes)
status = solver.SearchForAllSolutions(model, solution_printer)
print('Status = %s' % solver.StatusName(status))
print('Number of solutions found: %i' % solution_printer.count)
def blocks(shapes):
''' Helper to enumerate all blocks. '''
for shape in shapes:
for block in shape:
yield block
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
''' Print a solution. '''
def __init__(self, variables):
cp_model.CpSolverSolutionCallback.__init__(self)
self.variables = variables
self.count = 0
def on_solution_callback(self):
self.count += 1
solution = [(self.Value(block[X]), self.Value(block[Y])) for shape in self.variables for block in shape]
print((W+3)*'-')
for y in range(0, H+1):
print('|' + ''.join(['#' if (x,y) in solution else ' ' for x in range(0, W+1)]) + '|')
print((W+3)*'-')
if __name__ == '__main__':
main()
给出：
...
------
| |
| ###|
| # |
| # |
------
------
| |
| ###|
| #|
| #|
------
Status = OPTIMAL
Number of solutions found: 60
免费多项式：
如果我们将单元格网格视为图形，则问题可以重新解释为 找到网格单元格的 k 分区，其中每个分区都有特定的大小，另外每个分区是 connected component. I.e. AFAICT there is no difference between a connected component and a polyomino 并且此答案的其余部分做出了该假设。
发现所有可能 "k-partitions of the cells of the grid where each partition has a specific size" 在 OR-Tools 约束规划中表达起来非常简单。但是 connectedness 部分是 hard AFAICT（我尝试了很长时间但失败了......）。我认为 OR-Tools 约束编程不是正确的方法。我注意到网络优化库的 OR-Tools C++ 参考在 connected components 上有一些内容可能值得一看，但我不熟悉它。另一方面，Python 中的朴素递归搜索解决方案非常可行。
这是一个 "by hand" 天真的解决方案。它非常慢，但对于您的 4x4 情况来说可以忍受。地址用于标识网格中的每个单元格。 （另请注意，wiki page 暗示像这种算法是一种天真的解决方案，看起来它建议一些更有效的解决类似的 polyomino 问题）。
import numpy as np
from copy import copy
from tabulate import tabulate
D = 4 # Dimension of square grid.
KCC = [5,4,2,2] # List of the sizes of the required k connected components (KCCs).
assert(sum(KCC) <= D*D)
VALID_CELLS = range(2,D*D)
def search():
solutions = set() # Stash of unique solutions.
for start in VALID_CELLS: # Try starting search from each possible starting point and expand out.
marked = np.zeros(D*D).tolist()
_search(start, marked, set(), solutions, 0, 0)
for solution in solutions: # Print results.
print(tabulate(np.array(solution).reshape(D, D)))
print('Number of solutions found:', len(solutions))
def _search(i, marked, fringe, solutions, curr_count, curr_part):
''' Recursively find each possible KCC in the remaining available cells the find the next, until none left '''
marked[i] = curr_part+1
curr_count += 1
if curr_count == KCC[curr_part]: # If marked K cells for the current CC move onto the next one.
curr_part += 1
if curr_part == len(KCC): # If marked K cells and there's no more CCs left we have a solution - not necessarily unique.
solutions.add(tuple(marked))
else:
for start in VALID_CELLS:
if marked[start] == 0:
_search(start, copy(marked), set(), solutions, 0, curr_part)
else:
fringe.update(neighbours(i, D))
while(len(fringe)):
j = fringe.pop()
if marked[j] == 0:
_search(j, copy(marked), copy(fringe), solutions, curr_count, curr_part)
def neighbours(i, D):
''' Find the address of all cells neighbouring the i-th cell in a DxD grid. '''
row = int(i/D)
n = []
n += [i-1] if int((i-1)/D) == row and (i-1) >= 0 else []
n += [i+1] if int((i+1)/D) == row and (i+1) < D**2 else []
n += [i-D] if (i-D) >=0 else []
n += [i+D] if (i+D) < D**2 else []
return filter(lambda x: x in VALID_CELLS, n)
if __name__ == '__main__':
search()
给出：
...
- - - -
0 0 1 1
2 2 1 1
4 2 3 1
4 2 3 0
- - - -
- - - -
0 0 4 3
1 1 4 3
1 2 2 2
1 1 0 2
- - - -
Number of solutions found: 3884
在 OR-Tools 中限制简单连接区域的一种相对直接的方法是将其边界限制为 circuit。如果你所有的多联骨牌的大小都小于 8，我们就不需要担心非单连的。
此代码找到所有 3884 个解决方案：
from ortools.sat.python import cp_model
cells = {(x, y) for x in range(4) for y in range(4) if x > 1 or y > 0}
sizes = [4, 2, 5, 2, 1]
num_polyominos = len(sizes)
model = cp_model.CpModel()
# Each cell is a member of one polyomino
member = {
(cell, p): model.NewBoolVar(f"member{cell, p}")
for cell in cells
for p in range(num_polyominos)
}
for cell in cells:
model.Add(sum(member[cell, p] for p in range(num_polyominos)) == 1)
# Each polyomino contains the given number of cells
for p, size in enumerate(sizes):
model.Add(sum(member[cell, p] for cell in cells) == size)
# Find the border of each polyomino
vertices = {
v: i
for i, v in enumerate(
{(x + i, y + j) for x, y in cells for i in [0, 1] for j in [0, 1]}
)
}
edges = [
edge
for x, y in cells
for edge in [
((x, y), (x + 1, y)),
((x + 1, y), (x + 1, y + 1)),
((x + 1, y + 1), (x, y + 1)),
((x, y + 1), (x, y)),
]
]
border = {
(edge, p): model.NewBoolVar(f"border{edge, p}")
for edge in edges
for p in range(num_polyominos)
}
for (((x0, y0), (x1, y1)), p), border_var in border.items():
left_cell = ((x0 + x1 + y0 - y1) // 2, (y0 + y1 - x0 + x1) // 2)
right_cell = ((x0 + x1 - y0 + y1) // 2, (y0 + y1 + x0 - x1) // 2)
left_var = member[left_cell, p]
model.AddBoolOr([border_var.Not(), left_var])
if (right_cell, p) in member:
right_var = member[right_cell, p]
model.AddBoolOr([border_var.Not(), right_var.Not()])
model.AddBoolOr([border_var, left_var.Not(), right_var])
else:
model.AddBoolOr([border_var, left_var.Not()])
# Each border is a circuit
for p in range(num_polyominos):
model.AddCircuit(
[(vertices[v0], vertices[v1], border[(v0, v1), p]) for v0, v1 in edges]
+ [(i, i, model.NewBoolVar(f"vertex_loop{v, p}")) for v, i in vertices.items()]
)
# Print all solutions
x_range = range(min(x for x, y in cells), max(x for x, y in cells) + 1)
y_range = range(min(y for x, y in cells), max(y for x, y in cells) + 1)
solutions = 0
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
def OnSolutionCallback(self):
global solutions
solutions += 1
for y in y_range:
print(
*(
next(
p
for p in range(num_polyominos)
if self.Value(member[(x, y), p])
)
if (x, y) in cells
else "-"
for x in x_range
)
)
print()
solver = cp_model.CpSolver()
solver.SearchForAllSolutions(model, SolutionPrinter())
print("Number of solutions found:", solutions)
我是 SAT 求解器领域的新手，需要有关以下问题的一些指导。
考虑到：
❶ 我在一个 4*4 的网格中选择了 14 个相邻的单元格
❷ 我有 5 polyominoes (A, B, C, D, E) 尺码 4, 2, 5, 2 和 1
❸ 这些多联骨牌是自由的，即它们的形状不固定，可以形成不同的图案
我如何使用 SAT 求解器计算所选区域（灰色单元格）内这 5 个自由多联骨牌的所有可能组合？
借用@spinkus 富有洞察力的回答和 OR-tools 文档，我可以制作以下示例代码（在 Jupyter Notebook 中运行）：
from ortools.sat.python import cp_model
import numpy as np
import more_itertools as mit
import matplotlib.pyplot as plt
%matplotlib inline
W, H = 4, 4 #Dimensions of grid
sizes = (4, 2, 5, 2, 1) #Size of each polyomino
labels = np.arange(len(sizes)) #Label of each polyomino
colors = ('#FA5454', '#21D3B6', '#3384FA', '#FFD256', '#62ECFA')
cdict = dict(zip(labels, colors)) #Color dictionary for plotting
inactiveCells = (0, 1) #Indices of disabled cells (in 1D)
activeCells = set(np.arange(W*H)).difference(inactiveCells) #Cells where polyominoes can be fitted
ranges = [(next(g), list(g)[-1]) for g in mit.consecutive_groups(activeCells)] #All intervals in the stack of active cells
def main():
model = cp_model.CpModel()
#Create an Int var for each cell of each polyomino constrained to be within Width and Height of grid.
pminos = [[] for s in sizes]
for idx, s in enumerate(sizes):
for i in range(s):
pminos[idx].append([model.NewIntVar(0, W-1, 'p%i'%idx + 'c%i'%i + 'x'), model.NewIntVar(0, H-1, 'p%i'%idx + 'c%i'%i + 'y')])
#Define the shapes by constraining the cells relative to each other
## 1st polyomino -> tetromino ##
# #
# #
# # #
# ### #
# #
################################
p0 = pminos[0]
model.Add(p0[1][0] == p0[0][0] + 1) #'x' of 2nd cell == 'x' of 1st cell + 1
model.Add(p0[2][0] == p0[1][0] + 1) #'x' of 3rd cell == 'x' of 2nd cell + 1
model.Add(p0[3][0] == p0[0][0] + 1) #'x' of 4th cell == 'x' of 1st cell + 1
model.Add(p0[1][1] == p0[0][1]) #'y' of 2nd cell = 'y' of 1st cell
model.Add(p0[2][1] == p0[1][1]) #'y' of 3rd cell = 'y' of 2nd cell
model.Add(p0[3][1] == p0[1][1] - 1) #'y' of 3rd cell = 'y' of 2nd cell - 1
## 2nd polyomino -> domino ##
# #
# #
# # #
# # #
# #
#############################
p1 = pminos[1]
model.Add(p1[1][0] == p1[0][0])
model.Add(p1[1][1] == p1[0][1] + 1)
## 3rd polyomino -> pentomino ##
# #
# ## #
# ## #
# # #
# #
################################
p2 = pminos[2]
model.Add(p2[1][0] == p2[0][0] + 1)
model.Add(p2[2][0] == p2[0][0])
model.Add(p2[3][0] == p2[0][0] + 1)
model.Add(p2[4][0] == p2[0][0])
model.Add(p2[1][1] == p2[0][1])
model.Add(p2[2][1] == p2[0][1] + 1)
model.Add(p2[3][1] == p2[0][1] + 1)
model.Add(p2[4][1] == p2[0][1] + 2)
## 4th polyomino -> domino ##
# #
# #
# # #
# # #
# #
#############################
p3 = pminos[3]
model.Add(p3[1][0] == p3[0][0])
model.Add(p3[1][1] == p3[0][1] + 1)
## 5th polyomino -> monomino ##
# #
# #
# # #
# #
# #
###############################
#No constraints because 1 cell only
#No blocks can overlap:
block_addresses = []
n = 0
for p in pminos:
for c in p:
n += 1
block_address = model.NewIntVarFromDomain(cp_model.Domain.FromIntervals(ranges),'%i' % n)
model.Add(c[0] + c[1] * W == block_address)
block_addresses.append(block_address)
model.AddAllDifferent(block_addresses)
#Solve and print solutions as we find them
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(pminos)
status = solver.SearchForAllSolutions(model, solution_printer)
print('Status = %s' % solver.StatusName(status))
print('Number of solutions found: %i' % solution_printer.count)
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
''' Print a solution. '''
def __init__(self, variables):
cp_model.CpSolverSolutionCallback.__init__(self)
self.variables = variables
self.count = 0
def on_solution_callback(self):
self.count += 1
plt.figure(figsize = (2, 2))
plt.grid(True)
plt.axis([0,W,H,0])
plt.yticks(np.arange(0, H, 1.0))
plt.xticks(np.arange(0, W, 1.0))
for i, p in enumerate(self.variables):
for c in p:
x = self.Value(c[0])
y = self.Value(c[1])
rect = plt.Rectangle((x, y), 1, 1, fc = cdict[i])
plt.gca().add_patch(rect)
for i in inactiveCells:
x = i%W
y = i//W
rect = plt.Rectangle((x, y), 1, 1, fc = 'None', hatch = '///')
plt.gca().add_patch(rect)
问题是我硬编码了 5 个 unique/fixed 多联骨牌，我不知道 如何定义约束，以便考虑到每个多联骨牌的每个可能模式（前提是可以）。
对于每个 polyonomino 和每个可能的左上角单元格，您都有一个布尔变量，用于指示此单元格是否是封闭矩形的左上角部分。
对于每个单元格和每个 polyomino，您都有一个布尔变量，用于指示该单元格是否被该 polyomino 占用。
现在，对于每个单元格和每个 polyomino，你有一系列的含义：选择左上角的单元格意味着每个单元格实际上都被这个 polyomino 占用。
然后是约束条件： 对于每个单元格，最多有一个多联骨牌占据它 对于每个 polyomino，它的左上角恰好有一个单元格。
这是一道纯布尔题。
编辑： 我在原始答案中漏掉了 "free" 这个词，并使用 OR-Tools 给出了固定多联骨牌的答案。添加了一个部分来回答以包含免费多联骨牌的解决方案 - 结果证明 AFAICT 很难在使用 OR-Tools 的约束编程中精确表达。
使用 OR-TOOLS 修复的多联骨牌：
是的，你可以用 constraint programming in OR-Tools. OR-Tools knows nothing about 2D grid geometry so you have to encode the geometry of each shape you have in terms of positional constraints. I.e. a shape is a collection of blocks / cells that must have a certain relationship to each other, must be within the bounds of the grid and must not overlap. Once you have your constraint model you just ask the CP-SAT Solver 解决它，在你的情况下，所有可能的解决方案。
这是一个非常简单的概念证明，在 4x4 网格上有两个矩形（您可能还想添加某种解释器代码以从形状描述到一组 OR-Tools 变量和约束在更大规模的问题中，因为手动输入约束有点乏味）。
from ortools.sat.python import cp_model
(W, H) = (3, 3) # Width and height of our grid.
(X, Y) = (0, 1) # Convenience constants.
def main():
model = cp_model.CpModel()
# Create an Int var for each block of each shape constrained to be within width and height of grid.
shapes = [
[
[ model.NewIntVar(0, W, 's1b1_x'), model.NewIntVar(0, H, 's1b1_y') ],
[ model.NewIntVar(0, W, 's1b2_x'), model.NewIntVar(0, H, 's1b2_y') ],
[ model.NewIntVar(0, W, 's1b3_x'), model.NewIntVar(0, H, 's1b3_y') ],
],
[
[ model.NewIntVar(0, W, 's2b1_x'), model.NewIntVar(0, H, 's2b1_y') ],
[ model.NewIntVar(0, W, 's2b2_x'), model.NewIntVar(0, H, 's2b2_y') ],
]
]
# Define the shapes by constraining the blocks relative to each other.
# 3x1 rectangle:
s0 = shapes[0]
model.Add(s0[0][Y] == s0[1][Y])
model.Add(s0[0][Y] == s0[2][Y])
model.Add(s0[0][X] == s0[1][X] - 1)
model.Add(s0[0][X] == s0[2][X] - 2)
# 1x2 rectangle:
s1 = shapes[1]
model.Add(s1[0][X] == s1[1][X])
model.Add(s1[0][Y] == s1[1][Y] - 1)
# No blocks can overlap:
block_addresses = []
for i, block in enumerate(blocks(shapes)):
block_address = model.NewIntVar(0, (W+1)*(H+1), 'b%d' % (i,))
model.Add(block[X] + (H+1)*block[Y] == block_address)
block_addresses.append(block_address)
model.AddAllDifferent(block_addresses)
# Solve and print solutions as we find them
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(shapes)
status = solver.SearchForAllSolutions(model, solution_printer)
print('Status = %s' % solver.StatusName(status))
print('Number of solutions found: %i' % solution_printer.count)
def blocks(shapes):
''' Helper to enumerate all blocks. '''
for shape in shapes:
for block in shape:
yield block
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
''' Print a solution. '''
def __init__(self, variables):
cp_model.CpSolverSolutionCallback.__init__(self)
self.variables = variables
self.count = 0
def on_solution_callback(self):
self.count += 1
solution = [(self.Value(block[X]), self.Value(block[Y])) for shape in self.variables for block in shape]
print((W+3)*'-')
for y in range(0, H+1):
print('|' + ''.join(['#' if (x,y) in solution else ' ' for x in range(0, W+1)]) + '|')
print((W+3)*'-')
if __name__ == '__main__':
main()
给出：
...
------
| |
| ###|
| # |
| # |
------
------
| |
| ###|
| #|
| #|
------
Status = OPTIMAL
Number of solutions found: 60
免费多项式：
如果我们将单元格网格视为图形，则问题可以重新解释为 找到网格单元格的 k 分区，其中每个分区都有特定的大小，另外每个分区是 connected component. I.e. AFAICT there is no difference between a connected component and a polyomino 并且此答案的其余部分做出了该假设。
发现所有可能 "k-partitions of the cells of the grid where each partition has a specific size" 在 OR-Tools 约束规划中表达起来非常简单。但是 connectedness 部分是 hard AFAICT（我尝试了很长时间但失败了......）。我认为 OR-Tools 约束编程不是正确的方法。我注意到网络优化库的 OR-Tools C++ 参考在 connected components 上有一些内容可能值得一看，但我不熟悉它。另一方面，Python 中的朴素递归搜索解决方案非常可行。
这是一个 "by hand" 天真的解决方案。它非常慢，但对于您的 4x4 情况来说可以忍受。地址用于标识网格中的每个单元格。 （另请注意，wiki page 暗示像这种算法是一种天真的解决方案，看起来它建议一些更有效的解决类似的 polyomino 问题）。
import numpy as np
from copy import copy
from tabulate import tabulate
D = 4 # Dimension of square grid.
KCC = [5,4,2,2] # List of the sizes of the required k connected components (KCCs).
assert(sum(KCC) <= D*D)
VALID_CELLS = range(2,D*D)
def search():
solutions = set() # Stash of unique solutions.
for start in VALID_CELLS: # Try starting search from each possible starting point and expand out.
marked = np.zeros(D*D).tolist()
_search(start, marked, set(), solutions, 0, 0)
for solution in solutions: # Print results.
print(tabulate(np.array(solution).reshape(D, D)))
print('Number of solutions found:', len(solutions))
def _search(i, marked, fringe, solutions, curr_count, curr_part):
''' Recursively find each possible KCC in the remaining available cells the find the next, until none left '''
marked[i] = curr_part+1
curr_count += 1
if curr_count == KCC[curr_part]: # If marked K cells for the current CC move onto the next one.
curr_part += 1
if curr_part == len(KCC): # If marked K cells and there's no more CCs left we have a solution - not necessarily unique.
solutions.add(tuple(marked))
else:
for start in VALID_CELLS:
if marked[start] == 0:
_search(start, copy(marked), set(), solutions, 0, curr_part)
else:
fringe.update(neighbours(i, D))
while(len(fringe)):
j = fringe.pop()
if marked[j] == 0:
_search(j, copy(marked), copy(fringe), solutions, curr_count, curr_part)
def neighbours(i, D):
''' Find the address of all cells neighbouring the i-th cell in a DxD grid. '''
row = int(i/D)
n = []
n += [i-1] if int((i-1)/D) == row and (i-1) >= 0 else []
n += [i+1] if int((i+1)/D) == row and (i+1) < D**2 else []
n += [i-D] if (i-D) >=0 else []
n += [i+D] if (i+D) < D**2 else []
return filter(lambda x: x in VALID_CELLS, n)
if __name__ == '__main__':
search()
给出：
...
- - - -
0 0 1 1
2 2 1 1
4 2 3 1
4 2 3 0
- - - -
- - - -
0 0 4 3
1 1 4 3
1 2 2 2
1 1 0 2
- - - -
Number of solutions found: 3884
在 OR-Tools 中限制简单连接区域的一种相对直接的方法是将其边界限制为 circuit。如果你所有的多联骨牌的大小都小于 8，我们就不需要担心非单连的。
此代码找到所有 3884 个解决方案：
from ortools.sat.python import cp_model
cells = {(x, y) for x in range(4) for y in range(4) if x > 1 or y > 0}
sizes = [4, 2, 5, 2, 1]
num_polyominos = len(sizes)
model = cp_model.CpModel()
# Each cell is a member of one polyomino
member = {
(cell, p): model.NewBoolVar(f"member{cell, p}")
for cell in cells
for p in range(num_polyominos)
}
for cell in cells:
model.Add(sum(member[cell, p] for p in range(num_polyominos)) == 1)
# Each polyomino contains the given number of cells
for p, size in enumerate(sizes):
model.Add(sum(member[cell, p] for cell in cells) == size)
# Find the border of each polyomino
vertices = {
v: i
for i, v in enumerate(
{(x + i, y + j) for x, y in cells for i in [0, 1] for j in [0, 1]}
)
}
edges = [
edge
for x, y in cells
for edge in [
((x, y), (x + 1, y)),
((x + 1, y), (x + 1, y + 1)),
((x + 1, y + 1), (x, y + 1)),
((x, y + 1), (x, y)),
]
]
border = {
(edge, p): model.NewBoolVar(f"border{edge, p}")
for edge in edges
for p in range(num_polyominos)
}
for (((x0, y0), (x1, y1)), p), border_var in border.items():
left_cell = ((x0 + x1 + y0 - y1) // 2, (y0 + y1 - x0 + x1) // 2)
right_cell = ((x0 + x1 - y0 + y1) // 2, (y0 + y1 + x0 - x1) // 2)
left_var = member[left_cell, p]
model.AddBoolOr([border_var.Not(), left_var])
if (right_cell, p) in member:
right_var = member[right_cell, p]
model.AddBoolOr([border_var.Not(), right_var.Not()])
model.AddBoolOr([border_var, left_var.Not(), right_var])
else:
model.AddBoolOr([border_var, left_var.Not()])
# Each border is a circuit
for p in range(num_polyominos):
model.AddCircuit(
[(vertices[v0], vertices[v1], border[(v0, v1), p]) for v0, v1 in edges]
+ [(i, i, model.NewBoolVar(f"vertex_loop{v, p}")) for v, i in vertices.items()]
)
# Print all solutions
x_range = range(min(x for x, y in cells), max(x for x, y in cells) + 1)
y_range = range(min(y for x, y in cells), max(y for x, y in cells) + 1)
solutions = 0
class SolutionPrinter(cp_model.CpSolverSolutionCallback):
def OnSolutionCallback(self):
global solutions
solutions += 1
for y in y_range:
print(
*(
next(
p
for p in range(num_polyominos)
if self.Value(member[(x, y), p])
)
if (x, y) in cells
else "-"
for x in x_range
)
)
print()
solver = cp_model.CpSolver()
solver.SearchForAllSolutions(model, SolutionPrinter())
print("Number of solutions found:", solutions)