天体ショーを整数計画問題で解いてみた。

ペンシルパズルの「天体ショー」を整数計画問題として定式化してソルバーを使って解いてみた。

これは数理最適化 Advent Calendar 2020の19日目の記事です。

天体ショー

ペンシルパズルだとピクロスや数独が有名だけど、天体ショーもその一つ。

いろいろなペンシルパズルがあるけど、自分はこの天体ショーが大好き。

問題の図が星を散りばめたようになっていて、まさに「天体ショー」
ルールが「点対称」を使ったものになっていて、「点対称」と「天体ショー」を掛けているのが素敵
解き終わったときにイラストが出てくるので、解く喜びが大きい
序盤はサクサクとテンポよく進み、中盤で考えどころがあって、大きなブロックを切り出せると興奮する
問題を作るのも比較的簡単（ドット絵を考えて、点対称な図形で区切っていけばいい）

そこで、天体ショーのソルバーを作ってみることを考えた。

人間が解くようにロジックを組んでいってもいいんだけど、整数計画問題として定式化さえできれば、汎用ソルバーで解くこともできる。なので、今回は整数計画問題として定式化することを考えてみた。

ただ、検索してみると、すでに先人が。

鉛筆パズルの整数計画法による解法定式化集
（52番目に天体ショーがある）

なので、これを実装するだけかなと思ったんだけど、ちゃんと読んでみると問題があることが分かった。

連結制約

基本的な方針は、あるマスがある星に属するかどうかを0-1変数で表現する先人の方法でOK。ただ、先人の定式化だと「ある星に属するマスがすべて連結している」という条件（以下、連結制約）をちゃんと表現できていないことに気がついた。

一応、「孤立したマスがない」「2マスだけが繋がったブロックがない」「L字型にだけ繋がったブロックがない」という条件は入っていて、これで大体OKだったらしい。

けど、次のような簡単な問題でも上記の制約だけではダメになる：

f:id:yamaimo0625:20201219224218p:plain

この問題の正しい答えは以下：

f:id:yamaimo0625:20201219225415p:plain

けど、先人の入れた制約だと、以下のような答えも正しいとなってしまう：

f:id:yamaimo0625:20201219230227p:plain

先人の制約はすべて守られているけど、明らかに上のブロックと下のブロックは真ん中のブロックと連結してないので、正しい答えになっていない。

じゃあ、どうすれば連結制約を正しく表現できるのか・・・？

これが想像以上に難しかった。ああでもない、こうでもないと悩み続け、昨日やっと答えが分かった。

全域木

答えは「ブロックが全域木になっていることを表現する」というもの。

任意のブロックについて、マスを頂点、マスとマスの境界を辺と考えると、これはグラフになっていて、「ブロックが連結である」ことと「グラフに全域木が存在する」ことは同値になっている。（※ここでは全域木の定義を「連結」で「閉路を持たない」グラフとする。単に「閉路を持たない」グラフは全域森）

そして、ブロックは星のあるマスから伸びていくと考えると、グラフに全域木が存在するなら、星のないマスはすべて、ただ1つの親を持つことになる。（逆に、星のないマスがすべて、ただ1つの親を持つなら、グラフに全域木が存在する）

そこで、 $D = \{\textrm{up}, \textrm{right}, \textrm{down}, \textrm{left}\}$ として、マス $n$ に対して、方向 $d \in D$ のマスが親かどうかを、変数 $p_{n, d} \in \{0, 1\}$ で表現することにする。さらに、星のないマス $n$ に対して、星のあるマスまでのパスの長さを変数 $l_n \in \mathbb{Z}_{\ge 0}$ で表現することにする。

このとき、次の3つが制約として必要になる：

星のない任意のマスについて、親はただ1つ存在する
$\sum_{d\in D} p_{n, d} = 1$
親のパスの長さは、子のパスの長さよりちょうど1小さい
$-M (1- p_{n, d}) \le l_{n} - l_{n'} - 1 \le M (1 - p_{n, d})$
（※ $M$ は十分に大きな正の整数で、子 $n$ と親 $n'$ は隣り合っていて、 $d$ は適切な方向とする）
親は子と同じ星に属する
$p_{n, d} + x_{n, m} - 1 \le x_{n', m}$
（※ 親 $n'$ は子 $n$ に隣接していて、 $d$ は適切な方向とし、変数 $x_{n, m} \in \{0, 1\}$ はマス $n$ が星 $m$ に属するかどうかを表すとする）

2つめの制約はちょっとテクニカルで、 $n'$ が親でないなら $-M \le l_{n} - l_{n'} - 1 \le M$ となるので実質意味のない制約になり、逆に $n'$ が親なら $l_{n} - l_{n'} - 1 = 0$ なので子のパスの長さが親よりちょうど1大きいという制約になる。

また、最後の制約は、 $p_{n, d} = 1$ （つまり $n'$ が $n$ の親）かつ $x_{n, m} = 1$ （つまり $n$ が $m$ に属する）のときに限り、左辺は1となって（それ以外は0以下）、そのときに限り $x_{n', m} = 1$ （つまり $n'$ も $m$ に属する）となる（それ以外は右辺は0でも1でもいい）ので、ちゃんと条件を表現できていることが分かると思う。

（2020-12-26追記）
2.は元の制約だと閉路になる場合があったので、修正した。

PICOS

あとは他の制約なども数式で表現してソルバーに食わせればいいだけで、PythonからPICOS（＋GLPK）を使って解いてみた。

ちなみに、PICOSはPuLPやOR-Toolsのようなライブラリなんだけど、類似したライブラリの中で一番Pythonっぽいコードが書けるので、自分は気に入ってる。ソルバーもデフォルトでCVXOPTがついてきて、他にもGLPKやSCIP、あと持ってればCPLEXやGurobiも使える。さらにはLPやMIPだけじゃなくてSOCPやSDPも解ける優れもの。（知名度だけは低い・・・）

まず、問題はテキストで以下のように書くことにした：

-----------
| | |o|*|o|
--o--------
| | | * | |
---------o-
|o|*|o| | |
-----------
| * | |*| |
-----------
| |o| | |o|
-----------

oが白い星で、*が黒い星。

これを読み込んで情報を保持するために、MarkerクラスとProblemクラスを定義：

from enum import Enum

class Marker:
    class Color(Enum):
        WHITE = 0
        BLACK = 1
    
    def __init__(self, color, row, col):
        self.color = color
        self.row = row
        self.col = col

class Problem:
    @classmethod
    def read_str(cls, contents):
        lines = contents.splitlines()
        row_size = len(lines) // 2
        col_size = len(lines[0]) // 2
        markers = []
        for r_idx, line in enumerate(lines[1:-1]):  # 最初と最後の行は冗長なので無視
            for c_idx, char in enumerate(line[1:-1]): # 最初と最後の列は冗長なので無視
                row = r_idx / 2
                col = c_idx / 2
                if char == 'o':
                    marker = Marker(Marker.Color.WHITE, row, col)
                    markers.append(marker)
                elif char == '*':
                    marker = Marker(Marker.Color.BLACK, row, col)
                    markers.append(marker)
        return cls(row_size, col_size, markers)
    
    @classmethod
    def read_file(cls, filepath):
        with open(filepath) as f:
            contents = f.read()
        return cls.read_str(contents)    
    
    def __init__(self, row_size, col_size, markers):
        self.row_size = row_size
        self.col_size = col_size
        self.markers = markers

そして、解く関数を以下のように定義：
※めっちゃ長い

import numpy as np
import math
import picos

# 方向
directions = ['up', 'right', 'down', 'left']

def solve_by_picos(problem, solver='glpk'):
    # 整数計画問題
    ip = picos.Problem()
    ip.set_objective('find')
    
    # 変数 ==========

    # (row,col)のマスがmarkerに属するかどうか
    belong = {}
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            for i, marker in enumerate(problem.markers):
                belong[row,col,marker] = picos.BinaryVariable(f'belong_{row}_{col}_{i}')

    # (row,col)のマスについて、directionのマスが全域木で親かどうか
    parent = {}
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            for direction in directions:
                parent[row,col,direction] = picos.BinaryVariable(f'parent_{row}_{col}_{direction}')
    # 上端、右端、下端、左端については、各方向が親になることはないので、変数を0に置き換える
    for row in range(problem.row_size):
        # 左端は左が親にならない
        parent[row,0,'left'] = 0
        # 右端は右が親にならない
        parent[row,problem.col_size-1,'right'] = 0
    for col in range(problem.col_size):
        # 上端は上が親にならない
        parent[0,col,'up'] = 0
        # 下端は下が親にならない
        parent[problem.row_size-1,col,'down'] = 0

    # (row,col)のマスの、全域木でのルートまでの距離（明確に属するマスをルートとする）
    distance = {}
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            distance[row,col] = picos.IntegerVariable(f'distance_{row}_{col}', lower=0)

    # 制約 ==========
    
    # 各マーカーについて、明確に属するマスを制約として表現する
    # （距離が0になるという制約も入れる）
    # このとき、明確に属するマスは親が不要なので、変数を0で置き換える
    for marker in problem.markers:
        row_range = range(math.floor(marker.row), math.ceil(marker.row)+1)
        col_range = range(math.floor(marker.col), math.ceil(marker.col)+1)
        for row in row_range:
            for col in col_range:
                ip.add_constraint(belong[row,col,marker] == 1)
                ip.add_constraint(distance[row,col] == 0)
                for direction in directions:
                    parent[row,col,direction] = 0

    # 各マスはいずれかのマーカーに属する
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            belong_sum = sum(belong[row,col,marker] for marker in problem.markers)
            ip.add_constraint(belong_sum == 1)

    # 点対称の位置にあるマスは同じマーカーに属する
    # NOTE:
    # 点対称なマスに関して、以下が成り立つ：
    #   (row + opposite_row) / 2 = marker.row, (col + opposite_col) / 2 = marker.col
    # なので、点対称なマスの位置は、以下で計算できる：
    #   opposite_row = 2*marker.row - row, opposite_col = 2*marker_col - col
    # まず、点対称なマスが枠外になる場合、そのマスがマーカーに属することはない。
    # そして、
    #   opposite_row < row
    # であれば、行についてマーカーを跨いで反対側に来ているので、すでに登録済み（再度登録はしない）であり、
    #   opposite_row == row かつ opposite_col <= col
    # のときも、マーカーのある行の列についてマーカーを跨いで反対側に来ているので、すでに登録済み（再度登録はしない）。
    for marker in problem.markers:
        for row in range(problem.row_size):
            opposite_row = int(2*marker.row - row)
            for col in range(problem.col_size):
                opposite_col = int(2*marker.col - col)
                
                opposite_is_out = (
                    (opposite_row < 0) or (problem.row_size <= opposite_row)
                    or (opposite_col < 0) or (problem.col_size <= opposite_col)
                )
                if opposite_is_out:
                    ip.add_constraint(belong[row,col,marker] == 0)
                else:
                    already_added = (
                        (opposite_row < row)
                        or ((opposite_row == row) and (opposite_col <= col))
                    )
                    if not already_added:
                        ip.add_constraint(belong[row,col,marker] == belong[opposite_row,opposite_col,marker])

    # 親が不要でないマスは、いずれか1つの方向が親になる
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            parent_sum = sum(parent[row,col,direction] for direction in directions)
            if type(parent_sum) != int:
                ip.add_constraint(parent_sum == 1)

    # 親の距離は子よりちょうど1小さい
    # NOTE:
    # distance[row,col]とdistance[parent_row,parent_col]について、
    # a. 親になっているなら 0 <= distance[row,col] - distance[parent_row,parent_col] - 1 <= 0
    # b. 親になっていないなら -∞ <= distance[row,col] - distance[parent_row,parent_col] - 1 <= ∞
    # が成り立つので、十分大きなMに対して
    # -M*(1-parent[row,col,direction]) <= distance[row,col] - distance[parent_row,parent_col] - 1 <= M*(1-parent[row,col,direction])
    # という制約で表現できる。
    M = problem.row_size * problem.col_size
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            for direction in directions:
                if type(parent[row,col,direction]) != int:
                    if direction == 'up':
                        parent_row = row - 1
                        parent_col = col
                    elif direction == 'right':
                        parent_row = row
                        parent_col = col + 1
                    elif direction == 'down':
                        parent_row = row + 1
                        parent_col = col
                    else:
                        parent_row = row
                        parent_col = col - 1
                    dist_expr = distance[row,col] - distance[parent_row,parent_col] - 1
                    ip.add_constraint(dist_expr >= -M * (1-parent[row,col,direction]))
                    ip.add_constraint(dist_expr <= M * (1-parent[row,col,direction]))

    # 親は子と同じマーカーに属する
    # NOTE:
    # directionの方向のマスが親（parent[row,col,direction] == 1）で、
    # 子のマスがmarkerに属する（belong[row,col,marker] == 1）ならば、
    # 親のマスもmarkerに属さなければならない（belong[parent_row,parent_col,marker] == 1）。
    # この論理は以下の制約で表現できる：
    #    parent[row,col,direction] + belong[row,col,marker] - 1 <= belong[parent_row,parent_col,marker]
    # （両方の条件が満たされたときだけ左辺は1になる（それ以外は0以下）ことに注意）
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            for direction in directions:
                if type(parent[row,col,direction]) != int:
                    if direction == 'up':
                        parent_row = row - 1
                        parent_col = col
                    elif direction == 'right':
                        parent_row = row
                        parent_col = col + 1
                    elif direction == 'down':
                        parent_row = row + 1
                        parent_col = col
                    else:
                        parent_row = row
                        parent_col = col - 1
                    for marker in problem.markers:
                        ip.add_constraint(parent[row,col,direction] + belong[row,col,marker] - 1 <= belong[parent_row,parent_col,marker])

    ip.solve(solver=solver)
    
    # 答えの構築
    # 各マスの属するマーカーが入ったndarrayを返す
    answer = np.zeros((problem.row_size, problem.col_size), dtype=np.object)
    for row in range(problem.row_size):
        for col in range(problem.col_size):
            for i, marker in enumerate(problem.markers):
                if belong[row,col,marker].value == 1:
                    answer[row,col] = marker
                    break
    
    return answer