Menemukan semua kombinasi polyomino gratis dalam area tertentu dengan SAT-solver (Python)

15

Saya baru di dunia pemecah SAT dan perlu beberapa panduan tentang masalah berikut.

Mengingat bahwa:

❶ Saya memiliki pilihan 14 sel yang berdekatan dalam kisi 4 * 4

❷ Saya memiliki 5 poliomino (A, B, C, D, E) dengan ukuran 4, 2, 5, 2 dan 1

Poly poliomino ini gratis , artinya bentuknya tidak tetap dan dapat membentuk pola yang berbeda

masukkan deskripsi gambar di sini

Bagaimana saya bisa menghitung semua kombinasi yang mungkin dari 5 poliomino bebas ini di dalam area yang dipilih (sel abu-abu) dengan SAT-solver?

Meminjam baik dari jawaban berwawasan @ spinkus dan dokumentasi OR-tools saya bisa membuat kode contoh berikut (berjalan dalam 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)

masukkan deskripsi gambar di sini

Masalahnya adalah bahwa saya memiliki 5 polyomino unik / tetap yang dikodekan dan saya tidak tahu bagaimana mendefinisikan kendala sehingga setiap pola yang mungkin untuk setiap polyomino diperhitungkan (asalkan dimungkinkan).

solub
sumber
Saya mendengar tentang Google OR-tools pertama kali. Apakah mungkin untuk menggunakan pustaka Python standar seperti itertools,numpy , networkx?
mathfux
Saya lebih suka menggunakan sat-solver, atau-tools lebih disukai.
solub
@ Solub cukup mudah untuk memodelkan / memecahkan masalah semacam ini menggunakan bahasa MiniZinc, karena ada batasan tingkat tinggi untuk menempatkan objek tidak beraturan di permukaan. Jika Anda mengikuti kursus gratis "Pemodelan Tingkat Lanjut untuk Optimalisasi Diskrit" di Coursera , Anda sebenarnya akan diajari cara melakukannya dan diberikan beberapa contoh praktis (dan lebih rumit). Or-Tools memiliki antarmuka untuk bahasa MiniZinc, sehingga Anda masih dapat memanfaatkan kekuatannya untuk menemukan solusi cepat.
Patrick Trentin
1
Tampaknya menarik, terima kasih atas penunjuknya. Tidak yakin itu akan menjawab masalah spesifik yang saya miliki (mendefinisikan batasan yang melibatkan polyomino gratis, tidak statis) tapi saya pasti akan melihatnya.
solub
1
Saya harus minta maaf, saya benar-benar lupa tentang pertanyaan ini. Ada pertanyaan terkait dalam minizinctag dengan jawaban terinci yang mencakup saran saya sebelumnya tentang penggunaan minizinc.
Patrick Trentin

Jawaban:

10

EDIT: Saya melewatkan kata "gratis" dalam jawaban asli dan memberikan jawaban menggunakan OR-Tools untuk polyomino tetap. Menambahkan bagian untuk menjawab untuk menyertakan solusi untuk poliomino gratis - yang ternyata AFAICT cukup sulit untuk diungkapkan secara tepat dalam pemrograman kendala dengan OR-Tools.

POLYOMINO TETAP DENGAN ATAU ALAT:

Ya, Anda bisa melakukannya dengan pemrograman kendala di OR-Tools. OR-Tools tidak tahu apa-apa tentang geometri kisi 2D sehingga Anda harus menyandikan geometri dari setiap bentuk yang Anda miliki dalam hal batasan posisi. Yaitu bentuk adalah kumpulan blok / sel yang harus memiliki hubungan tertentu satu sama lain, harus dalam batas-batas grid dan tidak boleh tumpang tindih. Setelah Anda memiliki model kendala, Anda cukup bertanya Solver CP-SAT untuk menyelesaikannya, dalam kasus Anda, untuk semua solusi yang mungkin.

Ini adalah bukti konsep yang sangat sederhana dengan dua bentuk persegi panjang pada kisi 4x4 (Anda mungkin juga ingin menambahkan semacam kode juru bahasa untuk beralih dari deskripsi bentuk ke sekumpulan variabel OR-Tools dan kendala dalam masalah skala yang lebih besar karena memasukkan kendala dengan tangan agak membosankan).

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

Memberi:

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

POLYOMINO GRATIS:

Jika kita menganggap kisi-kisi sel sebagai grafik, masalahnya dapat ditafsirkan kembali sebagai menemukan k-partisi dari sel-sel kisi di mana setiap partisi memiliki ukuran tertentu dan di samping itu setiap partisi adalah komponen yang terhubung . Yaitu AFAICT tidak ada perbedaan antara komponen yang terhubung dan polyomino dan sisa jawaban ini membuat asumsi itu.

Menemukan semua kemungkinan "k-partisi dari sel-sel grid di mana setiap partisi memiliki ukuran tertentu" cukup sepele untuk diekspresikan dalam pemrograman kendala OR-Tools. Tetapi bagian keterhubungannya adalah AFAICT yang sulit (saya sudah mencoba dan gagal cukup lama ...). Saya pikir pemrograman kendala OR-Tools bukan pendekatan yang tepat. Saya perhatikan referensi OR-Tools C ++ untuk pustaka optimasi jaringan memiliki beberapa hal pada komponen yang terhubung yang mungkin patut dilihat, tapi saya tidak terbiasa dengannya. Di sisi lain, solusi pencarian rekursif naif dengan Python cukup bisa dilakukan.

Inilah solusi naif "dengan tangan". Ini cukup lambat tetapi dapat ditoleransi untuk kasing 4x4 Anda. Alamat digunakan untuk mengidentifikasi setiap sel dalam kisi. (Juga perhatikan jenis halaman wiki menyinggung sesuatu seperti algoritma ini sebagai solusi naif dan sepertinya itu menyarankan beberapa yang lebih efisien untuk masalah polyomino serupa).

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

Memberi:

...
-  -  -  -
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
sumber
Ini sangat membantu, terima kasih banyak. Satu hal yang bermasalah adalah bahwa contoh Anda hanya berfungsi untuk poliomino bentuk tetap, pertanyaannya adalah tentang poliomino bebas (jumlah sel tetap tetapi dengan bentuk berbeda, pertanyaan akan diedit untuk kejelasan). Mengikuti contoh Anda, kami harus membuat hard-code setiap bentuk yang memungkinkan (+ rotasi + refleksi) untuk setiap polyomino ukuran S ... yang tidak layak. Pertanyaannya tetap, apakah mungkin untuk mengimplementasikan kendala tersebut dengan alat-OR?
solub
Oh ketinggalan bagian "bebas". Hmmm, well masalahnya bisa diletakkan "temukan 5-partisi dari 25-omino di mana 25-omino dibatasi ke kisi WxH, dan masing-masing 5 partisi juga X-omino untuk X = (7,6,6 , 4,2) .. ". Saya kira itu mungkin dilakukan di OR-Tools tapi baunya akan lebih mudah hanya dengan menerapkan CSP back tracking depth untuk mencari ini secara langsung: Temukan 25-ominos yang mungkin. Untuk setiap 25-omino yang memungkinkan, lakukan penelusuran CSP mundur dengan memilih X yang membangun X-omino dalam 25 domino, hingga Anda menemukan solusi lengkap atau harus mundur.
spinkus
Menambahkan sesuatu seperti solusi berbasis pencarian langsung yang naif yang saya singgung di komentar sebelumnya untuk kelengkapan.
spinkus
5

Salah satu cara yang relatif mudah untuk membatasi wilayah yang terhubung sederhana di OR-Tools adalah membatasi perbatasannya menjadi sirkuit . Jika semua polyominos Anda berukuran kurang dari 8, kami tidak perlu khawatir tentang yang tidak terhubung secara mudah.

Kode ini menemukan semua 3884 solusi:

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

Untuk setiap polyonomino, dan setiap sel kiri atas yang mungkin, Anda memiliki variabel boolean yang menunjukkan apakah sel ini adalah bagian kiri atas dari persegi panjang terlampir.

Untuk setiap sel dan setiap polyomino, Anda memiliki variabel boolean yang menunjukkan apakah sel ini ditempati oleh polyomino ini.

Sekarang, untuk setiap sel dan setiap polyomino, Anda memiliki serangkaian implikasi: sel kiri atas dipilih menyiratkan setiap sel sebenarnya ditempati oleh polyomino ini.

Kemudian kendala: untuk setiap sel, paling banyak satu polyomino menempatinya untuk setiap polyomino, ada persis satu sel yang merupakan bagian kiri atas.

ini adalah masalah boolean murni.

Laurent Perron
sumber
Terima kasih banyak atas jawabannya ! Jujur saya tidak tahu bagaimana menerapkan ini dengan atau-alat, apakah ada contoh (dari contoh python yang tersedia disediakan) yang akan Anda sarankan secara khusus untuk membantu saya memulai?
solub
Saya benar-benar minta maaf karena saya tidak benar-benar memahami jawaban Anda. Tidak yakin apa yang dimaksud dengan "melampirkan persegi panjang" atau bagaimana "untuk setiap sel dan setiap poliomino" akan diterjemahkan dalam kode (nested 'for' loop?). Bagaimanapun, maukah Anda memberi tahu saya jika penjelasan Anda membahas kasus poliomino gratis (pertanyaan telah diedit untuk kejelasan).
solub