Localizando todas as combinações de poliomatos livres dentro de uma área específica com um SAT-solver (Python)

15

Eu sou novo no mundo dos solucionadores de SAT e precisaria de algumas orientações sobre o seguinte problema.

Considerando que:

❶ Eu tenho uma seleção de 14 células adjacentes em uma grade 4 * 4

Have Eu tenho 5 poliaminos (A, B, C, D, E) dos tamanhos 4, 2, 5, 2 e 1

Poly esses poliaminoinos são livres , ou seja, sua forma não é fixa e pode formar padrões diferentes

insira a descrição da imagem aqui

Como posso calcular todas as combinações possíveis desses 5 poliominos livres dentro da área selecionada (células em cinza) com um solucionador SAT?

Tomando emprestado da resposta perspicaz do @ spinkus e da documentação das ferramentas OR, eu poderia criar o seguinte código de exemplo (executado em um 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)

insira a descrição da imagem aqui

O problema é que eu codifiquei 5 poliaminoós únicos / fixos e não sei como definir as restrições para que cada padrão possível para cada polioino seja levado em consideração (desde que seja possível).

solub
fonte
Soube pela primeira vez sobre as ferramentas OR do Google. É possível utilizar as bibliotecas padrão do Python, como itertools, numpy, networkx?
mathfux 13/01
Eu preferiria usar um solucionador de sat, ou ferramentas, de preferência.
solub 13/01
@solub é muito fácil modelar / resolver esse tipo de problema usando a linguagem MiniZinc, já que existem restrições de alto nível para colocar objetos irregulares em uma superfície. Se você seguir o curso gratuito "Modelagem Avançada para Otimização Discreta" no Coursera , será realmente ensinado a fazê-lo e receberá alguns exemplos práticos (e mais complexos). O Or-Tools possui uma interface para o idioma MiniZinc, para que você ainda possa aproveitar seu poder para encontrar uma solução rápida.
Patrick Trentin
11
Parece interessante, obrigado pelo ponteiro. Não tenho certeza de que responderá ao problema específico que eu tenho (definindo restrições que envolvem poliaminoinos livres, não estáticos), mas definitivamente vou dar uma olhada nisso.
solub 14/01
11
Devo me desculpar, eu tinha esquecido completamente desta questão. Houve uma pergunta relacionada na minizinctag com uma resposta detalhada que aborda minha sugestão anterior sobre o uso minizinc.
Patrick Trentin

Respostas:

10

EDIT: Eu perdi a palavra "livre" na resposta original e dei resposta usando o OR-Tools para poliaminoinos fixos. Foi adicionada uma seção para responder a fim de incluir uma solução para poliaminoinos gratuitos - que o AFAICT se mostra bastante difícil de expressar com precisão na programação de restrições com o OR-Tools.

POLIOMINOS FIXOS COM FERRAMENTAS:

Sim, você pode fazer isso com programação de restrição no OR-Tools. O OR-Tools não sabe nada sobre geometria de grade 2D, portanto, é necessário codificar a geometria de cada forma que você possui em termos de restrições posicionais. Ou seja, uma forma é uma coleção de blocos / células que devem ter um certo relacionamento entre si, devem estar dentro dos limites da grade e não devem se sobrepor. Depois de ter seu modelo de restrição, basta solicitar ao CP-SAT Solver que o resolva, no seu caso, para todas as soluções possíveis.

Aqui está uma prova de conceito realmente simples, com duas formas de retângulo em uma grade 4x4 (você provavelmente também desejaria adicionar algum tipo de código de intérprete para passar das descrições de forma para um conjunto de variáveis ​​e restrições do OR-Tools em um problema de maior escala) já que inserir as restrições manualmente é um pouco 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()

Dá:

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

POLIOMINOS LIVRES:

Se considerarmos a grade de células como um gráfico, o problema pode ser reinterpretado ao encontrar uma partição k das células da grade em que cada partição tem um tamanho específico e, além disso, cada partição é um componente conectado . Ou seja, não há diferença entre um componente conectado e um poliomino, e o restante desta resposta faz essa suposição.

A descoberta de todas as partições-k possíveis das células da grade em que cada partição possui um tamanho específico é bastante trivial para expressar na programação de restrições do OR-Tools. Mas a parte da conexão é difícil AFAICT (tentei e falhei por um bom tempo ...). Acho que a programação de restrições do OR-Tools não é a abordagem correta. Notei que a referência do OR-Tools C ++ para as bibliotecas de otimização de rede tem algumas coisas nos componentes conectados que podem valer uma olhada, mas não estou familiarizado com isso. Por outro lado, a solução de pesquisa recursiva ingênua em Python é bastante factível.

Aqui está uma solução ingênua "manual". É muito lento, mas é suportável para o seu caso 4x4. Os endereços são usados ​​para identificar cada célula na grade. (Observe também que a página da wiki faz alusão a algo como esse algoritmo como uma solução ingênua e parece sugerir alguns mais eficientes para problemas similares de polioma).

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

Dá:

...
-  -  -  -
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
spinkus
fonte
Isso é muito útil, muito obrigado. Uma coisa que é problemática é que o seu exemplo funciona apenas para poliaminos de formas fixas; a questão é sobre poliaminos livres (número fixo de células, mas com formas diferentes, a pergunta será editada para maior clareza). Seguindo o seu exemplo, teríamos que codificar com firmeza todas as formas possíveis (+ rotações + reflexões) para cada poliomino de tamanho S ... que não seja viável. As questões permanecem: é possível implementar essas restrições com as ferramentas de OR?
solub 14/01
Oh, perdi a parte "grátis". Hmmm, bem, o problema pode ser colocado "encontre uma partição 5 de um 25-omino em que o 25-omino esteja restrito a uma grade WxH, e cada uma das 5 partições também seja X-omino para X = (7,6,6 , 4,2) .. ". Eu acho que é possível fazer no OR-Tools, mas parece que seria mais fácil implementar a profundidade de rastreamento de retorno do CSP primeiro, procure diretamente isso: Encontre 25 ominos possíveis. Para cada 25-omino possível, faça uma busca de CSP de retorno, escolhendo um X que constrói um X-omino dentro do dominó de 25, até encontrar uma solução completa ou ter que voltar atrás.
spinkus 15/01
Adicionado algo como a ingênua solução baseada em pesquisa direta a que me referi no comentário anterior por ser completo.
spinkus 20/01
5

Uma maneira relativamente direta de restringir uma região simplesmente conectada no OR-Tools é restringir sua borda a ser um circuito . Se todos os seus polyominos tiverem tamanho menor que 8, não precisamos nos preocupar com os não conectados simplesmente.

Este código encontra todas as soluções 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)
Anders Kaseorg
fonte
4

Para cada polyonomino e cada célula superior esquerda possível, você tem uma variável booleana que indica se essa célula é a parte superior esquerda do retângulo anexo.

Para cada célula e cada poliomino, você tem uma variável booleana que indica se essa célula está ocupada por esse poliomino.

Agora, para cada célula e cada poliomino, você tem uma série de implicações: a célula superior esquerda é selecionada implica que cada célula é realmente ocupada por esse poliomino.

Em seguida, as restrições: para cada célula, no máximo um polioino ocupa-o para cada polioino, existe exatamente uma célula que é sua parte superior esquerda.

este é um problema booleano puro.

Laurent Perron
fonte
Muito obrigado pela resposta ! Sinceramente, não tenho idéia de como implementar isso com o or-tools. Existe algum exemplo (dos exemplos disponíveis em python fornecidos) que você sugeriria em particular para me ajudar a começar?
solub
Sinto muito, pois não estou realmente entendendo sua resposta. Não sei ao que "retângulo delimitador" está se referindo ou como "para cada célula e cada poliomino" seria traduzido no código (aninhado 'for' loop?). De qualquer forma, você se importaria de me dizer se a sua explicação aborda o caso de poliaminoinos livres (a pergunta foi editada para maior clareza).
solub