Question: I am new to the world of SAT solvers and would need some guidance regarding the following problem. Considering that: ❶ I have a selection of 14 adjacent cells in a 4*4 grid  ❷ I have 5 polyominoes (A, B, C, D, E) of sizes 4, 2, 5, 2 and 1  ❸ these polyominoes are free, i.e. their shape is not fixed and can form different patterns  How can I compute all the possible combinations of these 5 free polyominoes inside the selected area (cells in grey) with a SAT-solver ? Borrowing both from @spinkus's insightful answer and the OR-tools documentation I could make the following example code (runs in a 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)   The problem is that I have hard-coded 5 unique/fixed polyominoes and I don't know to how define the constraints so as each possible pattern for each polyomino is taken into account (provided it is possible).

Answer: One relatively straightforward way to constrain a simply connected region in OR-Tools is to constrain its border to be a circuit. If all your polyominos are to have size less than 8, we don’t need to worry about non-simply connected ones. This code finds all 3884 solutions: 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)