Индексируем небо

Потребность в пространственном поиске возникает довольно часто. Равно как и в пространственном JOIN’е — нахождении пересечения двух наборов пространственных объектов. Далеко не всегда хочется привлекать тяжелую артиллерию. Что ж, попробуем придумать способ решить проблему “малой кровью, могучим ударом”.

Введение

Известно, что индексирование более чем одномерных объектов представляет собой существенную проблему, хотя, потребность в этом достаточно велика. Возникающие трудности носят объективный характер т.к. что бы мы ни делали, файл индекса на физическом носителе последователен (одномерен) и мы обязаны предоставить некоторый порядок обхода плоскости (пространства). При этом есть интуитивное желание расположить близкие в пространстве точки не очень далеко и в индексном файле. К сожалению, без разрывов это невозможно.

Отношение к пространственному индексированию в GIS сообществе вообще двоякое. С одной стороны, изрядное количество членов данного сообщества уверено, что проблемы то никакой и нет, случись такая необходимость, они справились бы с ней на два счета. А с другой стороны, есть мнение, что “все уже украдено до нас”, большинство современных СУБД как-то умеют решать данную проблему, некоторые довольно давно и, в целом, успешно.

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

Здесь мы рассмотрим самый простой случай индексации пространственных данных — работу с точечными объектами. Будем использовать честное кодирование координат с помощью заметающей кривой. Конечно, есть и другие методы, pixel-block индексы, R-деревья, квадродеревья, на худой конец. Возможно, это тема для другой статьи, но сейчас мы описываем именно работу с точечными объектами именно с помощью заметающей кривой.
Не будем утверждать, что предложенный метод идеален, но у него как минимум два достоинства — он очень прост и весьма эффективен. С его помощью буквально «на коленке» можно решать вполне серьезные задачи.

В качестве данных возьмем звездный атлас и его имитации случайными числами. Каждый набор данных состоит из некоторого числа объектов, каждый из которых описывается:

• X – координата, 32-х битовое целочисленное значение
• Y – координата, 32-х битовое целочисленное значение
• VAL – 32-х битовое целочисленное значение

Тестовые данные

Наборов данных у нас несколько:

1. 55M — Данный набор состоит из 55,368,239 объектов, координаты которых являются случайными числами в интервале 0… 2**27 по X и 0… 2**26 по Y, VAL всегда 0.
2. 526M — Данный набор состоит из 526,280,881 объектов, координаты которых являются случайными числами в интервале 0… 2**27 по X и 0… 2**27 по Y, VAL всегда 0.
3. USNO-SA2.0 — Это каталог из 55.368.239 звезд, полученный из ftp://ftp.nofs.navy.mil ^[1].
   X- координата есть (RA (Right Ascension) in decimal hours)*15*3600*100
   Y – координата — ((DEC in decimal degrees)+90)*3600*100
   VAL – некоторое значение, в котором, помимо всего прочего, содержатся яркости в красном и синем диапазонах.
4. USNO-A2.0 — Это каталог из 526,280,881 звезды, который может быть получен из ftp://ftp.nofs.navy.mil ^[1]. Для этих данных приводятся лишь оценки. Значения для его объектов аналогичны USNO-SA2.0

Теория

Идея заметающей кривой заключается в том, чтобы пронумеровать дискретное пространство. Пространственный объект представляется точкой в N-мерном пространстве. Точка на плоскости представляет сама себя, прямоугольный экстент, описывающий плоскую фигуру — это 4-мерная точка etc. Такая точка – число фиксированной длины, с которым довольно просто иметь дело. В частности, такое число можно поместить в обычное B-дерево. Поисковый запрос разбивается на непрерывные отрезки заметающей кривой и для каждого отрезка выполняется подзапрос. Все точки, попавшие в интервалы подзапросов – искомые.

Таким образом, нам нужна кривая, которая посетит каждую точку нашего (двумерного в данном случае) дискретного пространства. Какие возможны варианты?

1. Igloo. 2-мерное пространство нумеруется строчка за строчкой подобно развертке в кинескопе. Для равномерно распределенных точечных объектов в двумерном пространстве и при небольших запросах этот метод идеален. Недостатками его являются:
   • Эффективно работает на 2-мерном пространстве, иначе порождает очень много подзапросов.
   • Для больших поисковых экстентов порождает большое количество подзапросов.
   • Неэффективен для данных, склонных к кластеризации

2. Самоподобные кривые. Заметающая кривая имеет фрактальную природу, элементарная ячейка (например, квадрат 2х2) нумеруется определенным образом, после чего подобным образом нумеруется каждая под-ячейка. Самой известной такой кривой является так называемый Z-order (или bit-interleaving), природа которого следует из названия. Нумерация плоскости с разрешением 64x64 выглядит как:
   
   Синий прямоугольник на картинке – типовой поисковый запрос, каждый раз, когда он пересекает кривую, порождается новый подзапрос. В сущности, при неосторожном использовании, данный метод порождает даже больше подзапросов, чем igloo. Этим и объясняется то, что данный метод практически не используется в СУБД. Хотя, само — подобие наделяет его чудесной особенностью – адаптацией к неравномерно распределенным данным. Кроме того, непрерывные отрезки кривой могут накрывать значительные области пространства, что полезно, если суметь этим воспользоваться.
3. Еще один вариант само — подобной кривой: П–образная кривая Гильберта:
   

   В отличие от Z-order’а, данная кривая обладает непрерывностью, что позволяет осуществлять поиск чуть эффективнее, значительным недостатком является относительно высокая стоимость вычисления координат из значения и наоборот.

   Третий вариант с самоподобием 2Х2: X-образная кривая пересекает сама себя и автор не знает прецедента ее использования для пространственной индексации.

   Известны также, кривая Мура ^[2]– аналог кривой Гильберта, кривая Серпинского ^[3], кривая Госпера ^[4], …

4. Мы же будем использовать т.н. каскадный Igloo:
   
   Для получения данной картинки мы разделили диапазон координат на куски по 3 бита и перемешали X & Y — т.е. значение заметающей кривой получается как:
   • биты с 0 по 3 — это биты с 0 по 3 X координаты
   • биты с 3 по 6 — это биты с 0 по 3 Y координаты
   • биты с 6 по 9 — это биты с 3 по 6 X координаты
   • биты с 9 по 12 — это биты с 3 по 6 Y координаты

   Таким образом получаем блоки 8Х8 со строчной разверткой внутри, которые сами образуют между собой строчную развертку.
   Но это для демонстрации принципа. Для индексации будем делать блоки по 17 бит. Если интересует, почему 17, об этом позже.

Результаты

Создание индекса

Где:

• Disk Size – размер получившегося индексного файла в килобайтах
• Build Time – время, потраченное на его создание в минутах/секундах
• MEM size – пиковый размер оперативной памяти в Mb, необходимый построителю индексов

Поиск в индексе

В нижеследующих таблицах приведены значения, характерные для запросов на построенном индексе. Приведены суммарные значения для 100,000 поисков по случайно выбранным (в заселенной области, конечно) квадратным экстентам.

• 55M
  

  Extent Size	Time (sec)	NObj	RDisk	σ(RDisk)
  2x2	58	23	118,435	0.39
  4x4	59	100	118,555	0.39
  10x10	59	600	118,918	0.4
  20x20	60	2,622	119,221	0.4
  120x120	62	89,032	123,575	0.45
  1200x1200	98	8,850,098	174,846	0.82
  7200x7200	408	318,664,063	567,137	1.3

  Где:

  • Extent Size – размер поискового экстента, единицей в случае USNO-SA(X).0 является 1 arc second
  • Time – время в секундах, потраченное на выполнение 100,000 запросов
  • NObj – число объектов, найденное всеми запросами
  • RDisk – суммарное число дисковых операций чтения (промахи мимо внутреннего кэша), данное значение гораздо показательнее времени т.к. даже файл размером в сотни мегабайт эффективно кэшируется операционной системой
  • σ(RDisk) – стандартное отклонение числа дисковых операций, отметим, что по результатам экспериментов худшее значение RDisk не превосходило среднего (RDisk ) + 12 * σ(RDisk)

• 526M
  

  Extent Size	Time (sec)	NObj	RDisk
  2x2	738	110	185,672
  4x4	741	445	185,734
  10x10	730	2,903	186,560
  20x20	774	11,615	187,190
  120x120	800	421,264	196,080
  1200x1200	1200	42,064,224	307,904
  7200x7200	3599	1,514,471,480	1,442,489

  Где колонки аналогичны 55M.
  Обращает на себя внимание, что, в то время, как число обращений к диску изменилось незначительно (на треть), время выполнения запросов выросло более чем в 10 раз – результат отсутствия эффективного кэширования файла индекса операционной системой.

• USNO-SA2.0
  

  Extent Size	Time (sec)	NObj	RDisk	σ(RDisk)
  2x2	48	28	143,582	0.5
  4x4	50	115	143,887	0.5
  10x10	45	657	144,085	0.5
  20x20	45	2,585	144,748	0.51
  120x120	47	94,963	151,223	0.56
  1200x1200	80	9,506,746	224,016	0.97
  7200x7200	387	345,165,845	842,853	2.97

  Где колонки аналогичны 55M

• USNO-A2.0 (оценка)
  

  Extent Size	Time (sec)	NObj	RDisk	σ(RDisk)
  2x2	~600	~130	~200,200	~0.4
  4x4	~600	~500	~200,200	~0.4
  10x10	~600	~3,000	~200,200	~0.4
  20x20	~600	~12,000	~200,200	~0.4
  120x120	~600	~450,000	~250,200	~0.5
  1200x1200	~1,000	~45,000,000	~300,200	~1.4
  7200x7200	~3500	~1,600,000,000	~1,500,200	~2.0

  Где колонки аналогичны 526M

Пространственный JOIN

В нижеследующих таблицах приведены значения, характерные для join’ов индексов. Под join’ом мы понимаем поиск всех объектов, расположенных в пределах некоторого экстента. Например, если параметром запроса является 0.25 arc second, то для каждого элемента из одного индекса будет выполнен поиск во втором индексе с экстентом +- 0.25, т.е. квадрат со сторонами 0.5x0.5 arc second и центром в референтной точке.

• USNO-SA2.0 vs USNO-SA2.0
  

  Extent Size	Time (sec)	NObj
  0.5x0.5	175	2
  2x2	191	410
  6x6	212	4,412

  Где:

  • Extent Size – поисковый экстент в arc seconds
  • Time (sec) – время выполнения запроса в секундах
  • NObj – число найденных объектов (не включая самих себя), в данном случае надо разделить пополам т.к. это auto-join

• 55M vs USNO-SA2.0
  

  Extent Size	Time (sec)	NObj
  0.5x0.5	150	925
  2x2	165	13,815
  6x6	181	122,295

  Где колонки аналогичны USNO-SA2.0 vs USNO-SA2.0.

• 526M vs 526M
  

  Extent Size	Time (sec)	NObj
  0.5x0.5	1601	0
  2x2	1813	916
  6x6	2180	9943

  Где колонки аналогичны USNO-SA2.0 vs USNO-SA2.0.

Как и почему это работает ?

Для начала отметим важные моменты:

1. Для больших файлов дисковое кэширование не работает, не стоит на него рассчитывать.
2. Время обработки запроса состоит преимущественно из времён дисковых чтений, всё остальное относительно мало.
3. Для случайно распределенных запросов не работает никакое кэширование, если индекс устроен как B-дерево, не стоит рассчитывать на кэш страниц и тратить на него излишнюю память.
4. Для пространственно близких запросов достаточно небольшого (несколько десятков страниц) кэша.

Итак, попробуем понять как устроены данные.

• Глобальный экстент — 2**26 Х 2**27 или 180Х360 градусов.
• 180 градусов это 180*60*60 = 648000 угловых секунд.
• 1<<26 = 67 108 864. Значит, одна угловая секунда ~ 100 отсчетов.
• Один объект приходится в среднем на 1600 = (2*648000*648000/500000000) квадратных угловых секунд ~(40X40) .
• Или на 40*100 = 4000 отсчетов.
• На страницу B-дерева помещается ~1000 объектов.
• Если бы мы хотели сделать страницу дерева соизмеримой с блоком каскадного Igloo, то сделали бы этот блок размером 4000 * sqrt(1000) => ~1<< 17
• Вот и магическое число 17.
• Итак, на среднестатистическую страницу приходится экстент (40*32) ~1200Х1200 угловых секунд.
• В данных такая страница физически скорее всего расположена на двух смежных блоках каскадного Igloo
• Поэтому, поисковые запросы размером до 600Х600 угловых секунд скорее всего потребуют только одного чтения листовой страницы.
• Учитывая, что дерево индекса 4-х уровневое и верхние два уровня эффективно кэшируются, средний запрос такого размера требует двух физических чтений с диска.
• Что и требовалось объяснить.

Так как же всё-таки осуществляется поиск?

1. мы разрезаем поисковый экстент на куски в соответствии с блочным устройством индекса.
2. каждый кусок экстента соответствует одному блоку индекса и обрабатывается независимо.
3. внутри блока мы находим границы диапазона поиска
4. вычитываем весь диапазон, отфильтровывая все, что попадает за пределы поискового экстента.

Даже немного обидно, до чего всё просто.

Итоги

1. Даже самый крупный из известных нам атласов пространственных объектов (~3,000,000,000 штук) может быть проиндексирован за время, измеряемое часами
2. Объем дисковой памяти, необходимый для хранения такого индекса, измеряется десятками (~30) гигабайт и почти линеен числу объектов
3. Объем оперативной памяти, необходимой для поиска, весьма невелик (единицы Mb), индекс может быть использован, в том числе и во встроенных системах
4. Для пространственного поиска по малому экстенту (размеры которого сравнимы со средним расстоянием между объектами) время поиска при этом приближается ко времени 2-3 дисковых операций, т.е. практически к физическому пределу
5. Для относительно больших запросов время поиска больше зависит от числа попавших в выборку объектов
6. В любом случае, время остается предсказуемым, что позволяет использовать данный индекс, в том числе и в системах реального времени

PS

Описанная работа проводилась в 2004 году в соавторстве с Александром Артюшиным из замечательной компании DataEast ^[5]. На тот момент автор и сам был сотрудником этой компании и, да, публикация осуществляется с их согласия, спасибо Евгению Моисееву emoiseev ^[6].
Железо в те времена было не чета нынешнему, но принципиальные моменты не изменились.

Автор: zzeng

Источник ^[7]