Znajdowanie wszystkich kombinacji wolnych poliominoów w określonym obszarze za pomocą SAT-solvera (Python)

15

Jestem nowy w świecie solverów SAT i potrzebuję wskazówek dotyczących następującego problemu.

Biorąc pod uwagę, że:

❶ Mam wybór 14 sąsiadujących komórek w siatce 4 * 4

❷ Mam 5 poliominoów (A, B, C, D, E) o rozmiarach 4, 2, 5, 2 i 1

Poly te poliominoes są wolne , tzn. Ich kształt nie jest ustalony i mogą tworzyć różne wzory

wprowadź opis zdjęcia tutaj

Jak obliczyć wszystkie możliwe kombinacje tych 5 wolnych poliominoów w wybranym obszarze (komórki w kolorze szarym) za pomocą solvera SAT?

Pożyczając zarówno z wnikliwej odpowiedzi @ spinkus, jak i dokumentacji narzędzi OR, mogłem utworzyć następujący przykładowy kod (działa w Notatniku Jupyter):

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)

wprowadź opis zdjęcia tutaj

Problem polega na tym, że mam zakodowane na stałe 5 unikalnych / stałych poliomino i nie wiem, jak zdefiniować ograniczenia, aby uwzględnić każdy możliwy wzór dla każdego poliomino (pod warunkiem, że jest to możliwe).

solub
źródło
Po raz pierwszy słyszę o narzędziach Google OR. Czy możliwe jest wykorzystanie standardowych bibliotek Pythona, takich jak itertools, numpy, networkx?
mathfux
Wolałbym użyć solvera lub narzędzi.
solub
@solub dość łatwo jest modelować / rozwiązywać tego rodzaju problemy za pomocą języka MiniZinc, ponieważ istnieją wysokie ograniczenia dotyczące umieszczania nieregularnych obiektów na powierzchni. Jeśli przejdziesz przez bezpłatny kurs „Zaawansowane modelowanie dla dyskretnej optymalizacji” na Coursera , tak naprawdę nauczysz się , jak to robić i podasz kilka praktycznych (i bardziej złożonych) przykładów. Or-Tools ma interfejs dla języka MiniZinc, więc nadal możesz wykorzystać jego moc, aby znaleźć szybkie rozwiązanie.
Patrick Trentin
1
Wydaje się interesujące, dzięki za wskaźnik. Nie jestem pewien, czy rozwiąże konkretny problem, który mam (definiowanie ograniczeń związanych z wolnymi poliominoami, a nie statycznymi), ale na pewno się temu przyjrzę.
solub
1
Muszę przeprosić, zupełnie zapomniałem o tym pytaniu. Nastąpił pokrewne pytanie w minizincznaczniku o szczegółową odpowiedź, która obejmuje moje wcześniejsze sugestie na temat korzystania minizinc.
Patrick Trentin

Odpowiedzi:

10

EDYCJA: Brakowało mi słowa „darmowy” w oryginalnej odpowiedzi i dałem odpowiedź za pomocą OR-Tools dla stałych poliomino. Dodano sekcję z odpowiedzią zawierającą rozwiązanie dla darmowych poliominoesów - które AFAICT okazuje się dość trudne do wyrażenia w programowaniu ograniczeń za pomocą OR-Tools.

STAŁE POLIOMINO Z OR-NARZĘDZIAMI:

Tak, możesz to zrobić z programowaniem ograniczeń w OR-Tools. OR-Tools nic nie wie o geometrii siatki 2D, więc musisz zakodować geometrię każdego kształtu pod względem ograniczeń pozycyjnych. Tj. Kształt to zbiór bloków / komórek, które muszą mieć ze sobą określoną relację, muszą znajdować się w granicach siatki i nie mogą się pokrywać. Gdy masz już model ograniczenia, po prostu poproś CP-SAT Solver o jego rozwiązanie we wszystkich możliwych rozwiązaniach.

Oto naprawdę prosty dowód koncepcji z dwoma prostokątnymi kształtami na siatce 4x4 (prawdopodobnie prawdopodobnie chciałbyś także dodać kod interpretera, aby przejść od opisów kształtów do zestawu zmiennych i ograniczeń OR-Tools w problemie na większą skalę ponieważ ręczne wprowadzanie ograniczeń jest nieco uciążliwe).

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

Daje:

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

DARMOWE POLIOMINO:

Jeśli weźmiemy pod uwagę siatkę komórek jako wykres, problem można ponownie zinterpretować jako znalezienie k-partycji komórek siatki, gdzie każda partycja ma określony rozmiar, a ponadto każda partycja jest połączonym komponentem . Tzn. AFAICT nie ma różnicy między połączonym komponentem a poliomino, a reszta tej odpowiedzi przyjmuje to założenie.

Znalezienie wszystkich możliwych „k-partycji komórek siatki, gdzie każda partycja ma określony rozmiar” jest dość trywialne do wyrażenia w programowaniu ograniczeń OR-Tools. Ale część związana z połączeniem jest trudna AFAICT (próbowałem i przez jakiś czas nie udawało mi się ...). Myślę, że programowanie ograniczeń OR-Tools nie jest właściwym podejściem. Zauważyłem, że odwołanie OR-Tools C ++ do bibliotek optymalizacji sieci zawiera pewne elementy na podłączonych komponentach, które mogą być warte obejrzenia, ale nie jestem z tym zaznajomiony. Z drugiej strony naiwne rozwiązanie wyszukiwania rekurencyjnego w Pythonie jest całkiem wykonalne.

Oto naiwne rozwiązanie „ręcznie”. Jest dość wolny, ale można go znieść w przypadku 4x4. Adresy służą do identyfikacji każdej komórki w siatce. (Zauważ też, że strona wiki nawiązuje do czegoś takiego jak ten algorytm jako naiwne rozwiązanie i wygląda na to, że sugeruje kilka bardziej wydajnych dla podobnych problemów z poliomino).

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

Daje:

...
-  -  -  -
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
źródło
To bardzo pomocne, dziękuję bardzo. Jedną z problematycznych rzeczy jest to, że twój przykład działa tylko w przypadku poliominoów o ustalonych kształtach, pytanie dotyczy wolnych poliominoesów (stała liczba komórek, ale o różnych kształtach, pytanie zostanie zredagowane dla jasności). Idąc za twoim przykładem, musielibyśmy zakodować na stałe każdy możliwy kształt (+ obroty + odbicia) dla każdego poliomino o rozmiarze S ... co nie jest wykonalne. Pozostaje pytanie, czy możliwe jest wdrożenie takich ograniczeń za pomocą narzędzi OR?
solub
Och, tęskniłem za „darmową” częścią. Hmmm, no cóż, problem można postawić "znaleźć 5-partycję 25-omino, gdzie 25-omino jest ograniczone do siatki WxH, a każda z 5 partycji jest również X-omino dla X = (7,6,6 , 4,2) .. ”. Wydaje mi się, że można to zrobić w OR-Tools, ale pachnie tak, jakby łatwiej byłoby po prostu wdrożyć głębokość śledzenia wstecz CSP, najpierw wyszukaj to bezpośrednio: Znajdź możliwe 25-ominos. Dla każdego możliwego 25-omino przeprowadź wyszukiwanie CSP z powrotem, wybierając X budujący X-omino w dominie 25, aż znajdziesz kompletne rozwiązanie lub musisz cofnąć się.
spinkus
Dodano coś w rodzaju rozwiązania opartego na naiwnym wyszukiwaniu bezpośrednim, o którym wspomniałem w poprzednim komentarzu, aby uzyskać kompletność.
spinkus
5

Jednym stosunkowo prostym sposobem ograniczenia prostego połączenia regionu w OR-Tools jest ograniczenie jego granicy do obwodu . Jeśli wszystkie twoje poliaminy mają mieć rozmiar mniejszy niż 8, nie musimy się martwić o nie tylko połączone.

Ten kod znajduje wszystkie 3884 rozwiązania:

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
źródło
4

Dla każdej polionomino i każdej możliwej lewej górnej komórki masz zmienną logiczną, która wskazuje, czy ta komórka jest lewą górną częścią otaczającego prostokąta.

Dla każdej komórki i każdego poliomina masz zmienną boolowską, która wskazuje, czy komórka jest zajęta przez to poliomino.

Teraz, dla każdej komórki i każdego poliomino, masz szereg implikacji: wybrana jest lewa górna komórka oznacza, że ​​każda komórka jest faktycznie zajęta przez to poliomino.

Następnie ograniczenia: dla każdej komórki, co najwyżej jedno poliomino zajmuje ją dla każdego poliomina, jest dokładnie jedna komórka, która jest jego górną lewą częścią.

jest to czysty problem boolowski.

Laurent Perron
źródło
Dziękuję bardzo za odpowiedź! Naprawdę nie mam pojęcia, jak zaimplementować to za pomocą or-tools. Czy istnieje jakiś przykład (z dostępnych przykładów w Pythonie), który sugerowałbyś w szczególności, aby pomóc mi zacząć?
solub
Jest mi naprawdę przykro, ponieważ tak naprawdę nie rozumiem twojej odpowiedzi. Nie jestem pewien, do czego odnosi się „otaczający prostokąt” lub w jaki sposób „dla każdej komórki i każdego poliomino” zostałby przetłumaczony w kodzie (zagnieżdżony „dla” pętli?). W każdym razie, czy mógłbyś mi powiedzieć, czy twoje wyjaśnienie dotyczy przypadku wolnych poliominoesów (pytanie zostało zredagowane dla jasności).
solub