Сжатие координат

В базе GeoNames 45 227 городов, где живёт от десяти тысяч человек. Широта и долгота у каждого записаны с точностью до тысячной доли градуса, и если завести таблицу «сколько людей живёт в этой точке» с клеткой на каждую возможную пару координат, в ней будет 180 001 × 360 001 — шестьдесят пять миллиардов клеток, почти все пустые. Городов же всего 45 тысяч. Ниже — они все. Выберите прямоугольник, и сцена посчитает, сколько в нём городов и людей. Потом нажмите «Сжатые координаты».

Карта перекраивается. Каждый город переезжает с настоящей широты на её номер среди всех различных широт, а с долготы — на номер долготы. Различных широт у городов 32 407, долгот — 38 493, и каждая получает полосу одинаковой ширины. Поэтому раздуваются полосы, где городов много хоть где-нибудь на свете: долготы от Гибралтара до Волги, Индостана, востока Азии, широты от тропиков до пятьдесят пятой параллели. Тихий океан, где на тысячи километров нет ни одной долготы с городом, сжимается в узкие полосы у левого и правого краёв карты, север за полярным кругом и юг за сороковой параллелью — тоже в тонкие полосы, Антарктида — в линию у нижнего края. А ответ про прямоугольник не меняется: в прямоугольнике от Атлантики до Днепра на обеих картах 10 061 город и 470 893 364 человека. Числа другие, порядок тот же.

Идея

Многим алгоритмам важны не сами значения, а только какое из них больше. Тогда каждое значение можно заменить его местом в отсортированном списке различных значений: самое маленькое станет нулём, следующее — единицей и так далее. Номера идут подряд от 0 до m − 1, где m — число различных значений, и по ним можно заводить обычные массивы, деревья Фенвика и отрезков, сколь бы огромными ни были исходные числа.

Это и есть сжатие координат. Слово «координаты» здесь условное: сжимают что угодно, что можно сравнивать, — даты, суммы в тенге, номера узлов карты, строки. В статистике то же самое называют переходом к рангам, в SQL — функцией DENSE_RANK, в pandas — factorize.

Рецепт

Три шага. Собрать все значения, которые встретятся, — и в данных, и в запросах, если запросы известны заранее. Отсортировать и выбросить повторы. Номер значения — его индекс в получившемся списке; обратный перевод номера в значение — просто обращение к списку по индексу.

Номер можно искать двумя способами. Словарь «значение → номер» отвечает за O(1) в среднем, но знает только те значения, что были при сжатии. Двоичный поиск в отсортированном списке — bisect_left в Python, sort.SearchInts в Go — работает за O(log m) и годится для любого числа: он отвечает, сколько различных значений меньше данного. Для значения из списка это и есть его номер, а для чужого — номер места, куда оно встало бы. Последняя строка листинга проверяет, что оба способа согласны на всех 45 227 широтах.

Попутно листинг сообщает, что Алматы на широте 43,252 получает номер 25 008 из 32 407 и что южнее него лежат 33 996 городов мира — 75,2%. Номер в отсортированном списке — то же, что место в рейтинге, и сжатие сразу отвечает на вопросы о долях и местах.

В numpy все три шага делает одна функция: np.unique(a, return_inverse=True) возвращает и список различных значений, и номер каждого элемента. В pandas — pd.factorize(a, sort=True); без sort=True номера раздаются в порядке первого появления, а такой порядок со значениями не связан, и сравнивать номера между собой уже нельзя.

Почему ответы не меняются

Пусть x₀ < x₁ < … < x_m−1 — различные значения, а id(v) — номер значения v, то есть такое i, что x_i = v. Для значений, которых в списке нет, договоримся о двух функциях: lb(a) — число различных значений, меньших a, и ub(b) — число различных значений, не больших b. Это и есть bisect_left и bisect_right.

Из первой части следует главное: любой алгоритм, который смотрит на значения только через сравнения «меньше» и «равно», на номерах сделает ровно те же шаги и даст тот же ответ. Сортировка, поиск, подсчёт инверсий, наибольшая возрастающая подпоследовательность, дерево Фенвика по значениям — все такие. Вторая часть переводит запрос с границами, которых нет среди данных, на язык номеров: прямоугольник «широта от 42,252 до 44,252» превращается в «номера широт от lb(42 252) до ub(44 252) − 1». Сцена показывает этот перевод под картой.

Чего теорема не обещает, видно по формулировке: про разности и суммы значений в ней ничего нет. Номера 7 и 8 могут стоять за числами 10 и 11, а могут — за 10 и 10⁹. Алгоритм, которому нужны расстояния — длины, площади, «ближе ли на километр», — на номерах ошибётся. К этому вернёмся в разделе о ловушке.

Два измерения: одного сжатия мало

Вернёмся к прямоугольникам. Если сжать обе оси, таблица «сколько людей в клетке» уменьшится с 6,5·10¹⁰ клеток до 32 407 × 38 493 ≈ 1,25·10⁹. В пятьдесят раз меньше, но по восемь байт на клетку это всё ещё десять гигабайт, почти все — нули. В двух измерениях сжатие убирает пустые строки и столбцы, но не пустые клетки.

Поэтому сжимают одну ось, а вдоль другой идут обходом. Отсортируем города с запада на восток и будем по одному складывать их население в дерево Фенвика, где ячейка — номер широты. Когда в дереве лежат все города западнее какой-то долготы, сумма по номерам широт от lo до hi − 1 — это население полосы «западнее этой долготы, между этими широтами». Прямоугольник — разность двух таких полос: до его восточного края и до западного. Каждый прямоугольник задаёт два момента обхода, и если отсортировать эти моменты вместе с городами, хватит одного прохода на все прямоугольники сразу. Дерево при этом — 32 407 ячеек по восемь байт, четверть мегабайта.

rectangle                            latitude            longitude  towns      people
a degree around Almaty       42.252 ..  44.252     75.912 ..  77.912     10     2217955
the Nile delta               29.800 ..  31.600     29.800 ..  32.400    137    36672070
the island of Java           -8.800 ..  -5.800    105.100 .. 114.600    243    51941160
the whole world             -90.000 ..  90.000   -180.000 .. 180.000  45227  4210365101

В квадрате два на два градуса вокруг Алматы — 10 городов и 2,2 миллиона человек, в прямоугольнике вокруг острова Ява — 243 города и почти 52 миллиона. Во всех 45 227 городах, по данным GeoNames, живёт 4,2 миллиарда человек — около половины человечества. Числа сверены перебором: листинг в конце сравнивает каждый ответ с прямым подсчётом.

Обход с отложенными вопросами — частный случай общего приёма: если все запросы известны заранее, их можно задавать не в том порядке, в каком они пришли, а в том, в каком удобно отвечать. Сам обход — заметающая прямая; другие приёмы того же рода — чтение запросов с конца в статье о системе непересекающихся множеств и алгоритм Мо в корневой декомпозиции.

Код

Оба листинга выше написаны так, как пишут в жизни. В Python сжатие — sorted(set(values)) и словарь по enumerate, в Go — slices.Sort, затем slices.Compact, который выбрасывает соседние повторы, и map. Перевести границы запроса в номера — bisect_left и bisect_right в Python, sort.SearchInts(ys, a) и sort.SearchInts(ys, b + 1) в Go: для целых чисел «первое, что больше b» — то же, что «первое, что не меньше b + 1».

Две мелочи, на которых ошибаются чаще всего. Первая: сжимать нужно всё, что будет сравниваться, — и данные, и границы запросов, если по ним заводится массив. Вторая: номер нового значения, которого не было при сжатии, словарь не знает, а двоичный поиск даёт номер соседа. Если позже по этому номеру что-то записать, запись попадёт в чужую ячейку.

Сколько это стоит

Сжатие — сортировка, O(n log n), и проход для выбрасывания повторов, O(n). Номер для каждого из n значений — ещё O(n) словарём или O(n log m) двоичным поиском. Память — список различных значений и массив номеров, O(n). Окупается это там, где дальше нужна структура по значениям: дерево Фенвика по номерам весит O(m), а не O(размаха чисел), и каждая операция с ним стоит O(log m).

Но окупается не сразу. Один прямоугольник быстрее посчитать перебором: 45 227 проверок против сортировки всех городов. Замер ниже — q случайных прямоугольников на тех же городах, двумя способами. Оси логарифмические.

Перебор растёт ровно пропорционально q. Обход сначала стоит на месте: почти всё его время — сортировка и проход по городам, и прямоугольники ему почти ничего не добавляют. В Python кривые пересекаются около 35 прямоугольников, а в Go — около 250: перебор в Go тратит на город меньше наносекунды, и обогнать его труднее. Дальше разрыв растёт вместе с q: на тысяче прямоугольников обход в Python быстрее перебора в 23 раза, на четырёх тысячах в Go — в 13 раз. Это общее правило подготовки данных: она стоит как один-два прохода по ним и окупается, когда по подготовленным данным ходят много раз.

Ловушка: пропавшие промежутки

Классическая олимпиадная задача. Вдоль забора одну за другой клеят афиши, новая закрывает старые. Концы афиш — числа до миллиарда. Сколько афиш видно хотя бы краешком? Покрасить забор по метрам нельзя — он слишком длинный. Сжать концы афиш и красить массив номеров — можно. Но посмотрите, что из этого выходит.

Три афиши из листинга: [1, 10], поверх неё [1, 4] и [6, 10]. На заборе видны все три — от первой остался пятый метр. Сожмём концы: значения 1, 4, 6, 10 получают номера 0, 1, 2, 3, вторая афиша занимает клетки 0–1, третья — 2–3, и первая закрыта целиком. Пятого метра между номерами 1 и 2 просто нет: сжатие сохранило порядок концов, но не то, что лежало между ними. Это ровно то, чего теорема не обещала.

Лекарство — сжимать не точки, а элементарные отрезки. Афиша [l, r] занимает метры от l до r включительно, то есть полуинтервал [l, r + 1). Разрежем забор во всех точках l и r + 1. Каждый кусок между соседними разрезами либо целиком под афишей, либо целиком вне её — внутри куска ни одна афиша не начинается и не кончается. Значит, кусок можно считать одной клеткой, и раскраска клеток в точности повторяет раскраску забора. В примере разрезы 1, 5, 6, 11, клеток три, и средняя, [5, 6), — тот самый пятый метр.

На двух тысячах случайных афиш на заборе длиной в миллиард миллиметров наивное сжатие насчитало 642 видимые афиши, а верное — 708. Шестьдесят шесть афиш видны только в промежутках между чужими концами, и наивное сжатие их потеряло. Такой ответ не падает и выглядит правдоподобно, поэтому эта ошибка живёт долго.

Из того же рассуждения следует правило для задач, где важны длины: храните рядом с номерами сами значения. Клетка i — это отрезок от x_i до x_i+1, и его длина x_i+1 − x_i никуда не делась, просто лежит в другом массиве. Так считают длину объединения отрезков и площадь объединения прямоугольников: проходят по клеткам и складывают длины покрытых.

Ещё три ловушки помельче.

• Равные значения. Сжатие склеивает равные в один номер. Если алгоритму нужны строгие неравенства — «строго возрастающая подпоследовательность», «строго меньше», — это то, что нужно. Если нужно различать равные элементы, сжимают пары (значение, индекс).
• Дробные числа. Сжимать можно и их, но равенство дробей коварно: 0,1 + 0,2 не равно 0,3, и два «одинаковых» значения, посчитанных по-разному, получат разные номера. Перед сжатием дроби округляют до нужной точности или переводят в целые, как координаты городов здесь — в тысячные доли градуса.
• Значения, пришедшие позже. Номера выдаются один раз на весь список. Новое значение между x₅ и x₆ вставить некуда: пришлось бы перенумеровать всё, что больше.

Когда значения приходят по одному

Последняя ловушка — главное ограничение приёма: сжатие требует знать все значения заранее. В олимпиадных задачах так почти всегда и есть — вход читается целиком, и все запросы тоже. В живой системе, где события идут потоком, обходятся иначе. Можно хранить структуру не по номерам, а по самим значениям, но создавать узлы дерева отрезков только там, где что-то лежит: такое «неявное» дерево на отрезке [0, 10⁹] тратит около 30 узлов на каждое добавление. Можно взять сбалансированное дерево поиска, которое хранит сами значения и умеет считать, сколько из них меньше данного. Оба способа медленнее массива, зато не требуют знать будущее.

Ту же трудность знают колоночные базы данных. Они хранят повторяющиеся значения столбца — например, названия городов — номерами из словаря, и если словарь сохраняет порядок, то сортировать и искать можно прямо по номерам. Карстен Биннинг, Штефан Хильденбранд и Франц Фербер в статье 2009 года о таких словарях отмечают ту же беду: коды фиксированной длины, сохраняющие порядок, трудно расширять новыми значениями.

Где работает

• Карты и маршруты. Номера узлов OpenStreetMap давно больше 2³²; в исходниках маршрутизатора OSRM так и написано в комментарии, и номер узла OSM там занимает 64 бита. А внутри OSRM нумерует узлы подряд 32-битными числами, и графы для поиска путей хранятся в массивах по этим номерам. Массив с ячейкой на каждый исходный номер — а их больше четырёх миллиардов — занял бы десятки гигабайт.
• Данные и базы. np.unique(…, return_inverse=True), pd.factorize, категориальные столбцы pandas, оконная функция DENSE_RANK в SQL, словарное кодирование столбцов в колоночных хранилищах — всё это сжатие значений в номера.
• Статистика. Ранговая корреляция Спирмена — обычная корреляция, посчитанная по номерам значений вместо самих значений. Поэтому её не сбивают выбросы: миллиардер в выборке получает номер на единицу больше соседа, а не в тысячу раз больший доход. Выравнивание гистограммы фотографии устроено похоже: яркость пикселя заменяется долей пикселей, которые не ярче его.
• Олимпиады. Ограничение «координаты до 10⁹, точек до 2·10⁵» — почти всегда подсказка: дерево Фенвика или отрезков строится по номерам. Так считают инверсии (дерево Фенвика), наибольшую возрастающую подпоследовательность за O(n log n), площадь объединения прямоугольников (заметающая прямая).

История

У приёма нет автора: им пользовались с тех пор, как массивы стали индексировать чем-то, кроме счётчика цикла. Но самая известная ранняя его запись связана с задачей Клее. В 1977 году Виктор Клее спросил в журнальной заметке, можно ли найти длину объединения n отрезков быстрее, чем за O(n log n). В том же году Джон Бентли в неопубликованных заметках «Solutions to Klee's rectangle problems» решил двумерный вариант — площадь объединения прямоугольников — заметающей прямой. Для неё он придумал дерево отрезков, а дерево отрезков в классическом виде строится не по координатам, а по элементарным отрезкам между отсортированными концами — тем самым кускам, которые спасли афиши в разделе о ловушке. Майкл Фредман и Брюс Вейде в 1978 году показали, что для отрезков в естественной модели вычислений быстрее O(n log n) нельзя.

Задачи

Первая — прямоугольники на карте, как в листинге, но запросов десятки тысяч, и перебор не уложится во время. Во второй — афиши на заборе; самое естественное решение теряет промежутки, и тесты это заметят.