Finden aller Kombinationen von freien Polyominoes innerhalb eines bestimmten Bereichs mit einem SAT-Solver (Python)

15

Ich bin neu in der Welt der SAT-Löser und benötige eine Anleitung zum folgenden Problem.

Bedenkt, dass:

❶ Ich habe eine Auswahl von 14 benachbarten Zellen in einem 4 * 4-Raster

❷ Ich habe 5 Polyominoes (A, B, C, D, E) der Größen 4, 2, 5, 2 und 1

❸ Diese Polyominoes sind frei , dh ihre Form ist nicht festgelegt und kann unterschiedliche Muster bilden

Geben Sie hier die Bildbeschreibung ein

Wie kann ich mit einem SAT-Solver alle möglichen Kombinationen dieser 5 freien Polyominoes innerhalb des ausgewählten Bereichs (Zellen in Grau) berechnen ?

Aus der aufschlussreichen Antwort von @ spinkus und der Dokumentation zu den OP-Tools konnte ich den folgenden Beispielcode erstellen (läuft in einem Jupyter-Notizbuch):

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)

Geben Sie hier die Bildbeschreibung ein

Das Problem ist, dass ich 5 eindeutige / feste Polyominos fest codiert habe und nicht weiß, wie ich die Einschränkungen definieren soll, damit jedes mögliche Muster für jedes Polyomino berücksichtigt wird (sofern dies möglich ist).

solub
quelle
Ich höre zum ersten Mal von Google OR-Tools. Ist es möglich , Standard - Python - Bibliotheken zu verwenden , wie zum Beispiel itertools, numpy, networkx?
Mathfux
Ich würde es vorziehen, einen Satellitenlöser oder Werkzeuge zu verwenden.
Solub
@solub Es ist ziemlich einfach, diese Art von Problem mit der MiniZinc-Sprache zu modellieren / zu lösen, da es hohe Einschränkungen für das Platzieren unregelmäßiger Objekte auf einer Oberfläche gibt. Wenn Sie den kostenlosen Kurs "Erweiterte Modellierung für die diskrete Optimierung" auf Coursera durchgehen, wird Ihnen tatsächlich beigebracht, wie das geht, und Sie erhalten einige praktische (und komplexere) Beispiele. Or-Tools verfügt über eine Schnittstelle für die MiniZinc-Sprache, sodass Sie die Möglichkeiten nutzen können, um eine schnelle Lösung zu finden.
Patrick Trentin
1
Scheint interessant, danke für den Zeiger. Ich bin mir nicht sicher, ob es das spezifische Problem beantworten wird, das ich habe (Definieren von Einschränkungen mit freien Polyominoes, nicht statischen), aber ich werde es mir auf jeden Fall ansehen.
Solub
1
Ich muss mich entschuldigen, ich hatte diese Frage völlig vergessen. Das Tag enthält eine verwandte Frageminizinc mit einer detaillierten Antwort, die meinen vorherigen Vorschlag zur Verwendung abdeckt minizinc.
Patrick Trentin

Antworten:

10

BEARBEITEN: Ich habe das Wort "frei" in der ursprünglichen Antwort verpasst und mit OR-Tools für feste Polyominoes geantwortet. Ein zu beantwortender Abschnitt wurde hinzugefügt, um eine Lösung für freie Polyominoes aufzunehmen - was sich bei AFAICT als ziemlich schwierig herausstellt, um es in der Einschränkungsprogrammierung mit OR-Tools genau auszudrücken.

FESTE POLYOMINOES MIT OR-TOOLS:

Ja, Sie können dies mit der Einschränkungsprogrammierung in OR-Tools tun . OR-Tools weiß nichts über 2D-Gittergeometrie, daher müssen Sie die Geometrie jeder Form in Bezug auf Positionsbeschränkungen codieren. Das heißt, eine Form ist eine Sammlung von Blöcken / Zellen, die eine bestimmte Beziehung zueinander haben müssen, sich innerhalb der Grenzen des Gitters befinden müssen und sich nicht überlappen dürfen. Sobald Sie Ihr Constraint-Modell haben, bitten Sie den CP-SAT-Solver , es in Ihrem Fall für alle möglichen Lösungen zu lösen.

Hier ist ein wirklich einfacher Proof of Concept mit zwei Rechteckformen in einem 4x4-Raster (Sie möchten wahrscheinlich auch eine Art Interpreter-Code hinzufügen, um von Formbeschreibungen zu einer Reihe von OR-Tools-Variablen und -Beschränkungen in einem größeren Problem zu gelangen da die Eingabe der Einschränkungen von Hand etwas mühsam ist).

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

Gibt:

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

KOSTENLOSE POLYOMINOES:

Wenn wir das Gitter von Zellen als Diagramm betrachten, kann das Problem so interpretiert werden , dass eine k-Partition der Zellen des Gitters gefunden wird, wobei jede Partition eine bestimmte Größe hat und zusätzlich jede Partition eine verbundene Komponente ist . Dh AFAICT gibt es keinen Unterschied zwischen einer verbundenen Komponente und einem Polyomino, und der Rest dieser Antwort macht diese Annahme.

Das Auffinden aller möglichen "k-Partitionen der Zellen des Gitters, in denen jede Partition eine bestimmte Größe hat", ist in der OR-Tools-Einschränkungsprogrammierung ziemlich trivial auszudrücken. Aber der Teil der Verbundenheit ist schwer AFAICT (ich habe es eine ganze Weile versucht und bin gescheitert ...). Ich denke, OR-Tools Constraint-Programmierung ist nicht der richtige Ansatz. Ich habe festgestellt, dass die OR-Tools C ++ - Referenz für die Netzwerkoptimierungsbibliotheken einige Informationen zu verbundenen Komponenten enthält, die möglicherweise einen Blick wert sind, aber ich bin damit nicht vertraut. Andererseits ist eine naive rekursive Suchlösung in Python ziemlich machbar.

Hier ist eine "von Hand" naive Lösung. Es ist ziemlich langsam, aber für Ihren 4x4-Koffer erträglich. Adressen werden verwendet, um jede Zelle im Raster zu identifizieren. (Beachten Sie auch, dass die Wiki-Seite auf etwas wie diesen Algorithmus als naive Lösung anspielt und anscheinend einige effizientere für ähnliche Polyomino-Probleme vorschlägt).

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

Gibt:

...
-  -  -  -
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
quelle
Das ist sehr hilfreich, vielen Dank. Problematisch ist, dass Ihr Beispiel nur für Polyominoes mit festen Formen funktioniert. Die Frage bezieht sich auf freie Polyominoes (feste Anzahl von Zellen, aber mit unterschiedlichen Formen, die Frage wird aus Gründen der Klarheit bearbeitet). Nach Ihrem Beispiel müssten wir jede mögliche Form (+ Rotationen + Reflexionen) für jedes Polyomino der Größe S ... hart codieren, was nicht realisierbar ist. Es bleibt die Frage, ob es möglich ist, solche Einschränkungen mit OP-Tools zu implementieren.
Solub
Oh, ich habe den "freien" Teil verpasst. Hmmm, nun, das Problem kann gestellt werden: "Finde eine 5-Partition eines 25-omino, wobei das 25-omino auf ein WxH-Gitter beschränkt ist und jede der 5 Partitionen auch X-omino für X = (7,6,6 ist 4,2) .. ". Ich denke, es ist in OR-Tools möglich, aber es riecht so, als wäre es einfacher, nur die CSP-Back-Tracking-Tiefe zu implementieren. Suchen Sie zuerst direkt danach: Finden Sie mögliche 25-Ominos. Führen Sie für jedes mögliche 25-omino eine CSP-Suche durch, indem Sie ein X auswählen, das ein X-omino innerhalb des 25-Dominos erstellt, bis Sie eine vollständige Lösung finden oder zurückverfolgen müssen.
Spinkus
Der Vollständigkeit halber habe ich so etwas wie die naive, auf direkter Suche basierende Lösung hinzugefügt, auf die ich im vorherigen Kommentar hingewiesen habe.
Spinkus
5

Eine relativ einfache Möglichkeit, einen einfach verbundenen Bereich in OR-Tools einzuschränken, besteht darin, seine Grenze auf eine Schaltung zu beschränken . Wenn alle Ihre Polyominos eine Größe von weniger als 8 haben sollen, müssen wir uns nicht um nicht einfach verbundene kümmern.

Dieser Code findet alle 3884 Lösungen:

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
quelle
4

Für jedes Polyonomino und jede mögliche obere linke Zelle haben Sie eine boolesche Variable, die angibt, ob diese Zelle der obere linke Teil des umschließenden Rechtecks ​​ist.

Für jede Zelle und jedes Polyomino haben Sie eine boolesche Variable, die angibt, ob diese Zelle von diesem Polyomino belegt ist.

Nun haben Sie für jede Zelle und jedes Polyomino eine Reihe von Implikationen: Die Auswahl der oberen linken Zelle impliziert, dass jede Zelle tatsächlich von diesem Polyomino besetzt ist.

Dann die Einschränkungen: Für jede Zelle belegt höchstens ein Polyomino es für jedes Polyomino, es gibt genau eine Zelle, die ihr oberer linker Teil ist.

Dies ist ein reines boolesches Problem.

Laurent Perron
quelle
Vielen Dank für die Antwort! Ich habe ehrlich gesagt keine Ahnung, wie ich dies mit or-tools implementieren soll. Gibt es ein Beispiel (aus den verfügbaren Python-Beispielen), das Sie insbesondere vorschlagen würden, um mir den Einstieg zu erleichtern?
Solub
Es tut mir wirklich leid, da ich Ihre Antwort nicht wirklich verstehe. Nicht sicher, worauf sich "einschließendes Rechteck" bezieht oder wie "für jede Zelle und jedes Polyomino" in Code übersetzt wird (verschachtelte 'for'-Schleife?). Wie auch immer, würde es Ihnen etwas ausmachen, mir zu sagen, ob Ihre Erklärung den Fall der freien Polyominos anspricht (die Frage wurde aus Gründen der Klarheit bearbeitet).
Solub