Encontrar todas las combinaciones de poliominoes libres dentro de un área específica con un solucionador SAT (Python)

Soy nuevo en el mundo de los solucionadores de SAT y necesitaría alguna orientación sobre el siguiente problema.

Teniendo en cuenta que:

❶ Tengo una selección de 14 celdas adyacentes en una cuadrícula de 4 * 4

❷ Tengo 5 poliominoes (A, B, C, D, E) de tamaños 4, 2, 5, 2 y 1

❸ estos poliominós son libres , es decir, su forma no es fija y pueden formar diferentes patrones

¿Cómo puedo calcular todas las combinaciones posibles de estos 5 poliominós libres dentro del área seleccionada (celdas en gris) con un solucionador SAT?

Tomando prestado tanto de la respuesta perspicaz de @ spinkus como de la documentación de las herramientas OR, podría hacer el siguiente código de ejemplo (se ejecuta en un 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)

El problema es que he codificado 5 polyominoes únicos / fijos y no sé cómo definir las restricciones para tener en cuenta cada posible patrón para cada poliomino (siempre que sea posible).

EDITAR: Me perdí la palabra "gratis" en la respuesta original y di respuesta usando las herramientas OR para poliominoes fijos. Se agregó una sección para responder para incluir una solución para poliominoes libres, que AFAICT resulta ser bastante difícil de expresar con precisión en la programación de restricciones con OR-Tools.

POLIOMINOES FIJOS CON HERRAMIENTAS OR:

Sí, puedes hacerlo con programación restringida en OR-Tools. OR-Tools no sabe nada sobre la geometría de cuadrícula 2D, por lo que debe codificar la geometría de cada forma que tenga en términos de restricciones de posición. Es decir, una forma es una colección de bloques / celdas que deben tener una cierta relación entre sí, deben estar dentro de los límites de la cuadrícula y no deben superponerse. Una vez que tenga su modelo de restricción, simplemente solicite al CP-SAT Solver que lo resuelva, en su caso, para todas las soluciones posibles.

Aquí hay una prueba de concepto realmente simple con dos formas rectangulares en una cuadrícula de 4x4 (probablemente también desee agregar algún tipo de código de intérprete para pasar de las descripciones de forma a un conjunto de variables y restricciones de herramientas OR en un problema de mayor escala ya que ingresar las restricciones a mano es un poco tedioso).

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()

Da:

...
------
|    |
| ###|
|  # |
|  # |
------
------
|    |
| ###|
|   #|
|   #|
------
Status = OPTIMAL
Number of solutions found: 60

POLIOMINOS GRATIS:

Si consideramos la cuadrícula de celdas como un gráfico, el problema se puede reinterpretar como encontrar una partición k de las celdas de la cuadrícula donde cada partición tiene un tamaño específico y, además, cada partición es un componente conectado . Es decir, AFAICT no hay diferencia entre un componente conectado y un poliomino y el resto de esta respuesta hace esa suposición.

Encontrar todas las posibles "particiones k de las celdas de la cuadrícula donde cada partición tiene un tamaño específico" es bastante trivial de expresar en la programación de restricciones de herramientas OR. Pero la parte de conectividad es difícil AFAICT (lo intenté y fallé durante bastante tiempo ...). Creo que la programación de restricciones de OR-Tools no es el enfoque correcto. Noté que la referencia de OR-Tools C ++ para las bibliotecas de optimización de red tiene algunas cosas sobre los componentes conectados que pueden valer la pena, pero no estoy familiarizado con eso. Por otro lado, la solución de búsqueda recursiva ingenua en Python es bastante factible.

Aquí hay una solución ingenua "a mano". Es bastante lento, pero es soportable para su estuche 4x4. Las direcciones se utilizan para identificar cada celda en la cuadrícula. (También tenga en cuenta que la página wiki alude a algo como este algoritmo como una solución ingenua y parece que sugiere algunos más eficientes para problemas de poliomino similares).

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()

Da:

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

Una forma relativamente sencilla de restringir una región simplemente conectada en OR-Tools es restringir su borde para que sea un circuito . Si todos sus poliominos tienen un tamaño inferior a 8, no debemos preocuparnos por los que no están simplemente conectados.

Este código encuentra todas las soluciones 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)

Para cada polyonomino, y cada posible celda superior izquierda, tiene una variable booleana que indica si esta celda es la parte superior izquierda del rectángulo delimitador.

Para cada celda y cada poliomino, tiene una variable booleana que indica si esta celda está ocupada por este poliomino.

Ahora, para cada celda y cada poliomino, tiene una serie de implicaciones: la celda superior izquierda seleccionada implica que cada celda está realmente ocupada por este poliomino.

Luego las restricciones: para cada celda, como máximo un poliomino lo ocupa para cada poliomino, hay exactamente una celda que es su parte superior izquierda.

Este es un problema booleano puro.