Поиск всех комбинаций свободных полиомино в определенной области с помощью SAT-решателя (Python)

15

Я новичок в мире решателей SAT, и мне нужно некоторое руководство по следующей проблеме.

Учитывая, что:

❶ У меня есть выбор из 14 соседних ячеек в сетке 4 * 4

❷ У меня 5 полиомино (A, B, C, D, E) размеров 4, 2, 5, 2 и 1

Poly эти полиомино являются свободными , то есть их форма не является фиксированной и может формировать различные структуры

введите описание изображения здесь

Как я могу вычислить все возможные комбинации этих 5 свободных полиомино в выделенной области (ячейки серого цвета) с помощью SAT-решателя?

Заимствуя как проницательный ответ @ spinkus, так и документацию по OR-tools, я мог бы сделать следующий пример кода (работает в блокноте 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)

введите описание изображения здесь

Проблема в том, что у меня есть жестко запрограммированные 5 уникальных / фиксированных полиомино, и я не знаю, как определить ограничения, чтобы каждый возможный шаблон для каждого полиомино учитывался (при условии, что это возможно).

solub
источник
Я впервые слышу о Google OR-tools. Можно ли использовать стандартные библиотеки Python , такие как itertools, numpy, networkx?
mathfux
Я бы предпочел использовать спутниковый решатель или инструменты.
Solub
@solub довольно легко смоделировать / решить эту проблему, используя язык MiniZinc, поскольку существуют высокоуровневые ограничения для размещения неправильных объектов на поверхности. Если вы пройдете бесплатный курс «Расширенное моделирование для дискретной оптимизации» на Coursera , вас действительно научат делать это и приведут несколько практических (и более сложных) примеров. Or-Tools имеет интерфейс для языка MiniZinc, так что вы можете использовать его возможности для быстрого поиска решения.
Патрик Трентин
1
Кажется интересным, спасибо за указатель. Не уверен, что это ответит на конкретную проблему, которая у меня есть (определение ограничений, связанных со свободными полиамино, а не статическими), но я обязательно на это посмотрю.
Solub
1
Я должен извиниться, я полностью забыл об этом вопросе. Там был связанный вопрос в minizincтеге с подробным ответом, который охватывает мое предыдущее предложение об использовании minizinc.
Патрик Трентин

Ответы:

10

РЕДАКТИРОВАТЬ: я пропустил слово «свободный» в исходном ответе и дал ответ, используя OR-Tools для фиксированных polyominoes. В ответ добавили раздел, включающий решение для свободных полиомино, которое, как оказалось, довольно сложно выразить AFAICT точно в программировании с помощью OR-Tools.

ФИКСИРОВАННЫЕ ПОЛИОМИНО С ИЛИ ИНСТРУМЕНТАМИ:

Да, вы можете сделать это с помощью программирования ограничений в OR-Tools. OR-Tools ничего не знает о геометрии двумерной сетки, поэтому вам необходимо кодировать геометрию каждой фигуры в терминах позиционных ограничений. Т.е. форма - это набор блоков / ячеек, которые должны иметь определенное отношение друг к другу, должны находиться в пределах сетки и не должны пересекаться. Когда у вас есть модель ограничений, вы просто попросите, чтобы CP-SAT Solver решил ее, в вашем случае, для всех возможных решений.

Вот действительно простое доказательство концепции с двумя прямоугольными фигурами на сетке 4x4 (вы также, вероятно, захотите добавить некоторый код интерпретатора, чтобы перейти от описаний фигур к набору переменных и ограничений OR-Tools в задаче большего масштаба так как ввод ограничений вручную немного утомителен).

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

дает:

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

БЕСПЛАТНЫЕ ПОЛИОМИНО:

Если мы рассматриваем сетку ячеек как граф, проблему можно переосмыслить как нахождение k-разбиения ячеек сетки, где каждый раздел имеет определенный размер и, кроме того, каждый раздел является связанным компонентом . Т.е. AFAICT нет никакой разницы между связанным компонентом и polyomino, и остальная часть этого ответа делает это предположение.

Нахождение всех возможных «k-разделов ячеек сетки, где каждый раздел имеет определенный размер» довольно тривиально, чтобы выразить это в программировании ограничений OR-Tools. Но связность - это тяжело (я пытался и терпел неудачу довольно долго ...). Я думаю, что программирование ограничений OR-Tools не является правильным подходом. Я заметил, что в справочнике OR-Tools C ++ для библиотек оптимизации сети есть кое-что о подключенных компонентах, на которые стоит обратить внимание, но я не знаком с этим. С другой стороны, наивное рекурсивное решение поиска в Python вполне выполнимо.

Вот это «от руки» наивное решение. Это довольно медленно, но терпимо для вашего случая 4х4. Адреса используются для идентификации каждой ячейки в сетке. (Также обратите внимание, что вики-страница ссылается на что-то вроде этого алгоритма как наивное решение и выглядит так, как будто предлагает несколько более эффективных для подобных проблем с полиомино).

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

дает:

...
-  -  -  -
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
источник
Это очень полезно, большое спасибо. Одна проблема, которая проблематична, это то, что ваш пример работает только для полиомино фиксированной формы, вопрос о свободных полиомино (фиксированное число ячеек, но с разными формами, вопрос будет отредактирован для ясности). Следуя вашему примеру, мы должны были бы жестко закодировать каждую возможную форму (+ повороты + отражения) для каждого polyomino размера S ... который не является жизнеспособным. Остается вопрос, возможно ли реализовать такие ограничения с помощью OR-инструментов?
Solub
Ой пропустил "свободную" часть. Хммм, ну проблема может быть поставлена ​​"найти 5-секцию 25-омино, где 25-омино ограничено сеткой WxH, и каждый из 5 секций также является X-омино для X = (7,6,6 , 4,2) .. ". Я полагаю, что это можно сделать в OR-Tools, но пахнет так, как будто бы проще реализовать глубину обратного отслеживания CSP, прежде чем искать это напрямую: Найти возможные 25-юноши. Для каждого возможного 25-омино выполните поиск CSP с обратным отслеживанием, выбрав X, строящий X-омино в пределах 25 домино, пока вы не найдете полное решение или вам не придется возвращаться.
Спинкус
Добавлено что-то наподобие простого решения, основанного на прямом поиске, которое я упоминал в предыдущем комментарии для полноты.
Спинкус
5

Одним из относительно простых способов ограничения односвязной области в OR-Tools является ограничение ее границы как схемы . Если все ваши polyominos должны иметь размер меньше 8, нам не нужно беспокоиться о не просто связанных.

Этот код находит все 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)
Андерс Касеорг
источник
4

Для каждого polyonomino и каждой возможной верхней левой ячейки у вас есть логическая переменная, которая указывает, является ли эта ячейка верхней левой частью окружающего прямоугольника.

Для каждой ячейки и каждого polyomino у вас есть логическая переменная, которая указывает, занята ли эта ячейка этим polyomino.

Теперь для каждой ячейки и каждого полиомино у вас есть ряд значений: верхняя левая ячейка выбрана означает, что каждая клетка фактически занята этим полиомино.

Тогда ограничения: для каждой ячейки, самое большее одно polyomino занимает это для каждого polyomino, есть точно одна ячейка, которая является ее верхней левой частью.

это чисто логическая проблема.

Лоран Перрон
источник
Большое спасибо за ответ! Я, честно говоря, понятия не имею, как реализовать это с помощью or-tools, есть ли какой-нибудь пример (из доступных приведенных примеров Python), который вы бы предложили, в частности, чтобы помочь мне начать?
Solub
Мне очень жаль, потому что я не совсем понимаю ваш ответ. Не уверен, что подразумевается под «вмещающим прямоугольником» или как «для каждой ячейки и каждого polyomino» переводиться в код (вложенный цикл for?). В любом случае, не могли бы вы сказать мне, если ваше объяснение касается случая свободных полиомино (вопрос был отредактирован для ясности).
Solub