Квадродерево и сетка: когда простое лучше сложного

Если ты знаешь, когда сетка лучше дерева, от чего зависит глубина квадродерева и как с ним связаны кривая Мортона и геохеши, листай до «Теперь сам».

Как думаешь, что выбрать

Четыре тысячи точек и запросы «что в этом окне». После kd-деревьев и деревьев диапазонов хочется взять что-то из них. Но есть способ, который пишется за десять строк: разрезать поле на сетку 32 × 32 и разложить точки по корзинам.

Ответь себе до чтения: всегда ли сетка хуже дерева? И что случится, если девять десятых точек собраны в одном районе? Ответ: на равномерных данных сетка проверяет почти только то, что нашла. На скоплении запрос окном 100 × 100 находит 412 точек, а проверяет 3600 - всю кучу из одной корзины.

На что это похоже

Почтовые ящики в многоквартирном доме: по ящику на квартиру, и почтальону не нужно никакого дерева. Но если половина почты города приходит в один офис, ящик офиса переполнен, а остальные пустуют. Хорошо бы делить ящики там, где почты много, - как районы города: в центре улицы мелкие и частые, на окраинах кварталы огромные. Квадродерево - это такое деление по плотности.

Механизм: корзины, деление на четыре и кривая Мортона

Сетка. Корзина точки (x, y) - пара (x · m // S, y · m // S). Запрос перебирает корзины, которые задевает окно, и проверяет их точки. Цена - число корзин плюс число точек в них. Если плотность равномерна, точек в корзине в среднем n / m², и всё дёшево. Если нет - корзина с кучей стоит как перебор.

Квадродерево. Квадрат с точками: если точек больше одной, делим на четыре равных квадрата и продолжаем в каждом. Чем гуще точки, тем мельче квадраты.

Глубина. Квадрат на глубине d имеет сторону S / 2ᵈ. Квадрат делится, только если в нём две точки, а расстояние между ними не больше его диагонали. Значит, глубина не больше log₂(S / c) + 1,5, где c - расстояние между ближайшими точками. Две точки на расстоянии 1 в поле 2¹⁶ могут дать глубину 16 - если их разделяет только самый мелкий разрез, как (0, 0) и (1, 0), - даже когда всего точек две. Оценка - сверху: точки (32767, 0) и (32768, 0) разойдутся уже на первом разрезе. Число узлов - O((глубина + 1) · n).

Кривая Мортона. Чередуй биты x и y: x = 101, y = 011 дают ключ 1 0 0 1 1 1. Числа по такому ключу идут в порядке обхода квадродерева: каждый квадрат дерева - отрезок подряд идущих ключей. Поэтому квадродерево не обязательно хранить деревом - достаточно отсортировать точки по ключу Мортона и положить в обычный индекс базы данных. Так устроены геохеши: короткий префикс ключа - большая ячейка, длинный - маленькая. Ячейки у геохеша не всегда квадратные: долготу и широту он делит по очереди, и при нечётном числе бит ячейка вытянута.

Правило

Сетка - корзины одного размера: дёшево на равномерных данных, дорого на скоплениях. Квадродерево делит квадрат на четыре, пока в нём больше одной точки; глубина ≤ log₂(S / c) + 1,5, где c - наименьшее расстояние между точками. Обход квадродерева - кривая Мортона, и квадраты дерева - отрезки ключей в отсортированном порядке.

Разбор: глубина квадродерева

64 случайные точки в поле 2¹⁶, потом те же 64 плюс две точки на расстоянии 1.

import random

def kvadrodrevo(tochki, x, y, s, glubina, stat):
    # квадрат [x, x + s) × [y, y + s); делим на четыре, пока в нём
    # больше одной точки; stat = [узлов, наибольшая глубина]
    stat[0] = stat[0] + 1
    stat[1] = max(stat[1], glubina)
    if len(tochki) <= 1:
        return
    h = s // 2
    chetverti = [[], [], [], []]
    for p in tochki:
        k = 0
        # номер четверти - ключ Мортона: бит x старший
        if p[0] >= x + h:
            k = k + 2
        if p[1] >= y + h:
            k = k + 1
        chetverti[k].append(p)
    for k in range(4):
        kvadrodrevo(chetverti[k], x + h * (k // 2), y + h * (k % 2), h,
                    glubina + 1, stat)

random.seed(4)
S = 2 ** 16
ravnomerno = []
for i in range(64):
    x = random.randint(0, S - 1)
    ravnomerno.append([x, random.randint(0, S - 1)])
for nabor in [ravnomerno, ravnomerno + [[40000, 40000], [40001, 40000]]]:
    stat = [0, 0]
    kvadrodrevo(nabor, 0, 0, S, 0, stat)
    print(len(nabor), stat)

64 [181, 7]
66 [225, 16]

Шестьдесят четыре точки: 181 узел, глубина 7 - обычная глубина для случайных точек. Добавили две точки на расстоянии 1 - и глубина стала 16: чтобы разделить их, квадрат пришлось делить до стороны 1, это log₂ 65 536 = 16 уровней. Узлов стало 225 - лишние 44 узла: 11 лишних делений по четыре ребёнка, вдоль цепочки квадратов, где лежит только эта пара.

Здесь и граница квадродерева: оно чувствительно к близким парам. Для гарантии O(n) узлов используют сжатые квадродеревья (цепочки без ветвлений заменяют одним узлом) или kd-деревья, где глубина всегда log n.

Разбор: сетка на равномерных точках и на скоплении

Четыре тысячи точек, сетка 32 × 32 - корзины по 312,5 × 312,5. Два набора: равномерный и «скопление», где 90% точек в квадрате 300 × 300 у центра. Два запроса: маленькое окно в центре и большое в стороне.

import random

def setka(tochki, m, S):
    # m × m корзин по сторону S / m
    korziny = {}
    for p in tochki:
        k = (p[0] * m // S) * m + p[1] * m // S
        if k not in korziny:
            korziny[k] = []
        korziny[k].append(p)
    return korziny

def zapros(korziny, m, S, x1, x2, y1, y2):
    # [найдено, проверено точек]
    naideno = 0
    provereno = 0
    for i in range(x1 * m // S, x2 * m // S + 1):
        for j in range(y1 * m // S, y2 * m // S + 1):
            for p in korziny.get(i * m + j, []):
                provereno = provereno + 1
                if x1 <= p[0] <= x2 and y1 <= p[1] <= y2:
                    naideno = naideno + 1
    return [naideno, provereno]

random.seed(9)
S = 10000
ravnomerno = []
skoplenie = []
for i in range(4000):
    x = random.randint(0, S - 1)
    ravnomerno.append([x, random.randint(0, S - 1)])
    # 90% точек - в квадрате 300 × 300 у центра
    if i % 10 != 0:
        x = random.randint(4850, 5149)
        skoplenie.append([x, random.randint(4850, 5149)])
    else:
        x = random.randint(0, S - 1)
        skoplenie.append([x, random.randint(0, S - 1)])
for nabor in [ravnomerno, skoplenie]:
    k = setka(nabor, 32, S)
    print(zapros(k, 32, S, 4900, 5000, 4900, 5000),
          zapros(k, 32, S, 1000, 2000, 1000, 2000))

[1, 5] [36, 51]
[412, 3600] [6, 7]

Равномерные точки: окно 100 × 100 нашло одну точку и проверило 5, окно 1000 × 1000 - 36 и 51. Почти без лишней работы. Скопление: то же маленькое окно нашло 412 точек и проверило 3600 - скопление легло в четыре соседние корзины, и окно задело все четыре: оно платит за каждую точку скопления. Окно в стороне - 6 и 7: там пусто, и сетка снова хороша.

Вывод для выбора: если данные равномерны и размер окна известен, сетка с корзиной размером с окно почти идеальна: запрос проверяет немногим больше, чем находит. Если плотность неравномерна - дерево: квадродерево, kd-дерево или R-дерево.

Что ломается без этого

Сетка «на глазок». Корзина намного больше запроса - каждый запрос перебирает лишнее; намного меньше - запрос перебирает тысячи пустых корзин. Размер корзины подбирают под типичный запрос.

Квадродерево на данных с дубликатами. Две одинаковые точки нельзя разделить никаким квадратом - деление не кончится. Нужен предел глубины или правило «в листе может быть несколько точек».

Кривая Мортона как «расстояние». Соседи по ключу почти всегда соседи в пространстве, но не наоборот: две точки по разные стороны большой границы квадрата близко лежат, а ключи у них далеко. Запрос окна по ключам режут на несколько отрезков.

Зачем это в работе

Геохеши и квадратные тайлы карт - кривая Мортона и квадродерево: адрес тайла - путь от корня. H3 от Uber делит Землю на шестиугольники с той же иерархической идеей (на каждом уровне ровно 12 пятиугольников - сферу шестиугольниками не покрыть). Сбалансированные квадродеревья, где соседние квадраты различаются по размеру не больше чем вдвое, - заготовка для сеток в расчётах: мелко там, где сложная граница, крупно внутри. Игровые движки держат объекты сцены в квадродеревьях и октодеревьях - трёхмерном аналоге.

Теперь сам

1. Поле 1024 × 1024, ближайшие точки на расстоянии 4. Оценка глубины квадродерева?
2. Ключ Мортона для x = 2 (10) и y = 3 (11), если, как в уроке, в каждой паре бит x идёт первым?
3. Почему сетка - плохой выбор для точек ресторанов города?