Find2D-BL アルゴリズムの実装

find2d_bl.py

# An implementation of Find2D-BL.
#
# Shinji Imahori, Yuyao Chien, Yuma Tanaka, and Mutsunori Yagiura.
# Enumerating bottom-left stable positions for rectangle placements with overlap.
# Journal of the Operations Research Society of Japan. Vol.57, No.1, March 2014, pp.45--61.
# https://www.jstage.jst.go.jp/article/jorsj/57/1/57_KJ00009350944/_article/-char/ja/

from __future__ import annotations

from dataclasses import dataclass
from enum import IntEnum
from typing import Optional

import numpy as np
import pytest


@dataclass(frozen=True)
class Point2D:
    x: int
    y: int


@dataclass(frozen=True)
class Rect2D:
    width: int
    height: int


@dataclass(frozen=True)
class NFP2D:
    x1: int
    y1: int
    x2: int
    y2: int


class VEdgeType(IntEnum):  # right < left になるように
    RIGHT = 0
    LEFT = 1


class HEdgeType(IntEnum):  # top < bottom になるように
    TOP = 0
    BOTTOM = 1


@dataclass(frozen=True)
class VEdge:
    x: int
    y: int
    type_: VEdgeType
    placement_index: int


@dataclass(frozen=True)
class HEdge:
    x: int
    y: int
    type_: HEdgeType
    placement_index: int


## -----------------------------------------------------------------------------
##
## Find2D-BL で利用する完全二分木
##
## -----------------------------------------------------------------------------


@dataclass
class BLNode:
    """Find2D-BL で利用する完全二分木のノード"""

    p_self: int
    p_min: int
    p_max: int
    parent: Optional[BLNode]
    left: Optional[BLNode]
    right: Optional[BLNode]

    def add(self, lambda_: int):
        self.p_self += lambda_
        self.p_min += lambda_
        # self.p_max += lambda_  # 現論文には記載があるが、なくても BL 法は実装できる。

    def update_min_max(self):
        assert self.left is not None and self.right is not None
        self.p_min = self.p_self + min(self.left.p_min, self.right.p_min)
        # self.p_max = self.p_self + max(self.left.p_max, self.right.p_max)


class BLTree:
    root: BLNode
    leaves: list[BLNode]

    def __init__(self, num_intervals: int):
        """Find2D-BL で利用する完全二分木の初期状態を生成する。
        Parameters
        ----------
        num_intervals
            葉数
        Returns
        -------
            初期状態の完全二分木
        """

        depth = int(np.ceil(np.log2(num_intervals - 1)))
        num_leaves = 2**depth  # 木が持つすべての葉の数
        num_nodes = 2 ** (depth + 1) - 1  # 木が持つすべてのノード数
        num_internal_nodes = num_nodes - num_leaves  # 内部ノード数
        num_dummy_leaves = num_leaves - num_intervals + 1  # ダミーノードの数

        nodes = [BLNode(0, 0, 0, None, None, None) for _ in range(num_nodes)]

        # 親子の参照を初期化する
        for i in range(num_internal_nodes):
            left = nodes[2 * i + 1]
            right = nodes[2 * i + 2]
            nodes[i].left = left
            nodes[i].right = right
            left.parent = nodes[i]
            right.parent = nodes[i]

        # ダミー葉は p_self=1 で初期化する（常に長方形を配置したくないノード）
        for i in range(num_nodes - num_dummy_leaves, num_nodes):
            nodes[i].p_self = 1
            nodes[i].p_min = 1
            nodes[i].p_max = 1

        # 内部ノードは子の情報を元に値を決める
        for i in reversed(range(num_internal_nodes)):
            nodes[i].update_min_max()

        self.root = nodes[0]
        self.leaves = nodes[num_internal_nodes:]

    def update_values(
        self, sweep_line_type: HEdgeType, left_index: int, right_index: int
    ):
        """新たに遭遇した NFP の辺に応じて、葉ごとの NFP の重複数を更新する。
        Parameters
        ----------
        sweep_line_type
            走査線のタイプ（上辺 or 下辺）
        left_index
            新たに遭遇した NFP の左辺に対応する葉のインデックス
        right_index
            新たに遭遇した NFP の右辺に対応する葉のインデックス
        """

        assert left_index <= right_index

        lambda_ = -1 if sweep_line_type == HEdgeType.TOP else +1
        left = self.leaves[left_index]
        right = self.leaves[right_index]

        left.add(lambda_)
        right.add(lambda_)

        while left.parent is not right.parent:
            # Add lambda_ to sibling nodes if necessary.
            if left is not left.parent.right:
                left.parent.right.add(lambda_)

            if right is not right.parent.left:
                right.parent.left.add(lambda_)

            # Update parents.
            left.parent.update_min_max()
            right.parent.update_min_max()

            left = left.parent
            right = right.parent

        # Update all nodes on the path to the root.
        node = left.parent  # Note that left.parent == right.parent
        while node is not None:
            node.update_min_max()
            node = node.parent

    def find_no_overwrap(self, alpha: int) -> Optional[int]:
        """左の葉から探索を始め、最初に NFP の重複数がゼロになる葉を返す。
        Parameters
        ----------
        alpha
            探索範囲の左端の葉のインデックス
        Returns
        -------
            alpha より右側で最初に NFP 重複数がゼロになる葉のインデックス。
            そのような葉がない場合は、None を返す。"""

        def _find(node: BLNode, offset: int, num_leaves: int) -> Optional[int]:
            if alpha >= offset + num_leaves:  # 探索範囲を外れた
                return None

            if node.p_min > 0:  # node の配下に NFP の重複数が 0 の葉が存在しない
                return None

            if node.left is None:  # node は葉
                # p_min の定義より、葉について p_min == p_self なので node.p_self == 0 が保証される
                return offset
            else:  # node は中間ノード
                k = num_leaves // 2
                # 左のノードの後で右のノードを探索する
                result = _find(node.left, offset, k)
                if result is None:
                    result = _find(node.right, offset + k, k)

                return result

        return _find(self.root, 0, len(self.leaves))


## -----------------------------------------------------------------------------
##
## Find2D-BL のメインルーチン
##
## -----------------------------------------------------------------------------


def _NFPs_to_edges(nfps: list[NFP2D]) -> tuple[list[VEdge], list[HEdge]]:
    """NFP を辺に分解する。"""

    left_right_edges = []
    top_bottom_edges = []

    for i, nfp in enumerate(nfps):
        left = VEdge(nfp.x1, nfp.y2, VEdgeType.LEFT, i)
        right = VEdge(nfp.x2, nfp.y2, VEdgeType.RIGHT, i)
        top = HEdge(nfp.x2, nfp.y2, HEdgeType.TOP, i)
        bottom = HEdge(nfp.x2, nfp.y1, HEdgeType.BOTTOM, i)

        left_right_edges.extend([left, right])
        top_bottom_edges.extend([top, bottom])

    return left_right_edges, top_bottom_edges


def _get_interval_indices(intervals: list[VEdge]) -> np.ndarray:
    """左辺・右辺の集合からインターバルを作る。"""

    # i 番目の既配置図形 placements[i] に対して、
    #   - 左辺の x 座標は left_right_edges[ interval_indices[i][0] ]
    #   - 右辺の x 座標は left_right_edges[ interval_indices[i][1] ]
    # になるように interval_indices を作る。

    num_placements = len(intervals) // 2
    placement_to_left_indices = [-1 for _ in range(num_placements)]
    placement_to_right_indices = [-1 for _ in range(num_placements)]

    for j, e in enumerate(intervals):
        if e.type_ == VEdgeType.LEFT:
            left_indices[e.placement_index] = j
        else:
            right_indices[e.placement_index] = j

    return placement_to_left_indices, placement_to_right_indices


def find2dbl(nfps: list[NFP2D], new_rect: Rect2D, box_height: int) -> Optional[Point2D]:
    # new_rect を配置可能な y 座標の上限値
    y_limit = box_height - new_rect.height
    # NFP から辺を抽出
    left_right_edges, top_bottom_edges = _NFPs_to_edges(nfps)
    # 辺をソートする
    intervals = sorted(left_right_edges, key=lambda e: (e.x, e.type_, e.y))
    sweep_lines = sorted(top_bottom_edges, key=lambda e: (e.y, e.type_, e.x))
    # sweep line の placement_index から左辺・右辺の index を引くための表を作る
    placement_to_left_indices, placement_to_right_indices = _get_interval_indices(
        intervals
    )
    # BLTree を作る
    tree = BLTree(len(intervals))
    # 下の辺から順に走査していく
    for sweep_line in sweep_lines:
        # 走査線に対応する既配置図形の index
        i = sweep_line.placement_index
        # 既配置図形の左辺・右辺に対応する葉ノードの index
        left_leaf = placement_to_left_indices[i]
        right_leaf = placement_to_right_indices[i] - 1
        # NFP の重複数を更新する
        tree.update_values(sweep_line.type_, left_leaf, right_leaf)
        # 箱の中に new_rect を配置できなくなった
        if sweep_line.y > y_limit:
            break
        # 走査線がまだ箱の内側に入っていない
        if sweep_line.y < 0:
            continue
        # 既配置図形の上辺ならば、new_rect を配置可能かもしれない
        if sweep_line.type_ == HEdgeType.TOP:
            # NFP 重複がゼロになる点を探す
            found_leaf = tree.find_no_overwrap(left_leaf)
            if found_leaf is not None:
                y = sweep_line.y
                x = intervals[found_leaf].x
                return Point2D(x, y)

    return None  # NFP 重複数がゼロになる点は見つからなかった


## -----------------------------------------------------------------------------
##
## ユーティリティ
##
## -----------------------------------------------------------------------------


def create_placements_of_box(box: Rect2D) -> list[tuple[Point2D, Rect2D]]:
    c = 1  # 適当な定数（1 以上の任意の数）
    return [
        # (Point2D(0, box.height), Rect(box.width, c)),  # top
        (Point2D(0, -c), Rect2D(box.width, c)),  # bottom
        (Point2D(-c, 0), Rect2D(c, box.height)),  # left
        (Point2D(box.width, 0), Rect2D(c, box.height)),  # right
    ]


def create_NFP(point: Point2D, rect: Rect2D, new_rect: Rect2D) -> NFP2D:
    x1 = point.x - new_rect.width
    y1 = point.y - new_rect.height
    x2 = point.x + rect.width
    y2 = point.y + rect.height
    return NFP2D(x1, y1, x2, y2)


def find_BL_point(
    placements: list[tuple[Point2D, Rect2D]], new_rect: Rect2D, box: Rect2D
) -> Optional[Point2D]:
    placements = create_placements_of_box(box) + placements
    nfps = [create_NFP(point, rect, new_rect) for point, rect in placements]
    return find2dbl(nfps, new_rect, box.height)


## -----------------------------------------------------------------------------
##
## テスト
##
## -----------------------------------------------------------------------------


@pytest.mark.parametrize(
    "new_rect, expected",
    [
        (Rect2D(1, 1), Point2D(0, 0)),  # 箱より十分小さい図形
        (Rect2D(5, 5), Point2D(0, 0)),  # 箱と同じ大きさの図形
        (Rect2D(6, 1), None),  # 箱より幅の広い図形
        (Rect2D(1, 6), None),  # 箱より高い図形
    ],
)
def test_何もない空間に図形を配置する(new_rect: Rect2D, expected: Optional[Point2D]):
    box = Rect2D(5, 5)
    actual = find_BL_point([], new_rect, box)

    assert expected == actual


@pytest.mark.parametrize(
    "placements, new_rect, expected",
    [
        # 既配置の図形の横に配置
        (
            [
                (Point2D(0, 0), Rect2D(3, 1)),
                (Point2D(3, 0), Rect2D(1, 2)),
            ],
            Rect2D(1, 1),
            Point2D(4, 0),
        ),
        # 既配置の図形の直上に配置
        (
            [
                (Point2D(0, 0), Rect2D(3, 1)),
                (Point2D(3, 0), Rect2D(1, 2)),
            ],
            Rect2D(2, 2),
            Point2D(0, 1),
        ),
        # 既配置の図形の上に少し空間を開けて配置
        (
            [
                (Point2D(0, 0), Rect2D(3, 1)),
                (Point2D(3, 0), Rect2D(1, 3)),
                (Point2D(0, 1), Rect2D(1, 1)),
            ],
            Rect2D(4, 1),
            Point2D(0, 3),
        ),
        # BL 点に配置されていない図形がある場合
        (
            [
                (Point2D(1, 1), Rect2D(2, 1)),
                (Point2D(1, 2), Rect2D(2, 1)),
                (Point2D(3, 0), Rect2D(1, 2)),
            ],
            Rect2D(2, 2),
            Point2D(3, 2),
        ),
        # 配置されている図形に重なりがある場合
        (
            [
                (Point2D(0, 0), Rect2D(3, 1)),
                (Point2D(1, 0), Rect2D(3, 2)),
                (Point2D(2, 1), Rect2D(1, 2)),
            ],
            Rect2D(2, 2),
            Point2D(0, 2),
        ),
        # 入らない
        (
            [
                (Point2D(0, 0), Rect2D(3, 1)),
                (Point2D(3, 0), Rect2D(1, 2)),
            ],
            Rect2D(5, 4),
            None,
        ),
    ],
)
def test_既配置の図形に対しBL点を計算する(
    placements: list[tuple[Point2D, Rect2D]],
    new_rect: Rect2D,
    expected: Optional[Point2D],
):
    box = Rect2D(5, 5)
    actual = find_BL_point(placements, new_rect, box)

    assert expected == actual