image

Индекс на основе Z-order кривой в сравнении с R-деревом имеет массу преимуществ, он:

  • реализован как обычное B-дерево, а мы знаем что
  • страницы B-дерева имеют лучшую заполняемость, кроме того,
  • Z-ключи сами по себе более компактны
  • B-дерево имеет естественный порядок обхода, в отличие от R-дерева
  • B-дерево быстрее строится
  • B-дерево лучше сбалансировано
  • B-дерево понятнее, не зависит от эвристики расщепления/слияния страниц
  • B-дерево не деградирует при постоянных изменениях
  • ...

Впрочем, у индексов на основе Z-order есть и недостаток — сравнительно низкая производительность :). Под катом мы попробуем разобраться с чем связан этот недостаток и можно ли что-то с этим сделать.

Грубо, смысл Z-кривой таков — мы чередуем разряды x & y координат в одном значении, в результате. Пример — (x = 33 (010 0001), y=22 (01 0110), Z=1577 (0110 0010 1001). Если у двух точек близкие координаты, то скорее всего и Z- значения у них будут близки.
Трехмерный (или более) вариант устроен так же, мы чередуем разряды трех (или более) координат.

А для поиска в некотором экстенте мы должны “растеризовать” все значения Z-кривой для этого экстента, найти непрерывные интервалы и в каждом осуществить поиск.

Вот иллюстрация для экстента [111000...111000][111400...111400] (675 интервалов), правый верхний угол (каждая непрерывная ломаная — один интервал):

image

А, например, для экстента [111000...111000][112000...112000] мы получаем 1688 непрерывных интервалов, очевидно, их число в основном зависит от длины периметра. Забегая вперед, на тестовых данных в этот экстент попало 15 точек в 6 интервалах.

Да, большая часть этих интервалов мала, вплоть до вырожденного — из одного значения. И тем не менее, даже если во всём этом экстенте всего одно значение, оно может быть в любом из интервалов. И, нравится нам это или нет, придется делать все 1688 подзапросов чтобы узнать сколько там точек в действительности.

Просмотреть всё от значения левого нижнего угла до правого верхнего — не вариант, разница между углами в данном случае — 3 144 768, нам придется просмотреть в три с лишним раза больше данных, и это далеко не самый худший случай. Например, экстент [499500...499500][500500...500500] даст диапазон в 135 263 808 значений, в ~135 раз больше площади экстента.

И тут мы можем задать традиционный вопрос —

А что, если ...


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

Теперь предположим, что экстент достаточно мал, а данные разрежены и шансов найти что-то мало. Выполним тот же запрос от угла до угла. Если ничего не нашлось, значит ничего и нет. В противном случае есть шанс. Но как мы видели, площадь заметаемая запросом от угла до угла может многократно превышать поисковый экстент и у нас нет никакого желания вычитывать заведомо ненужные данные. Поэтому мы не будем просматривать весь курсор, а возьмем из него лишь минимальное Z-значение. Для этого запрос выполняется с (order by) и (top 1).

Итак, у нас есть некоторое значение. Допустим, это значение не из нашего экстента, что это может дать? Самое время вспомнить, что у нас есть сортированный массив [1...N] диапазонов подзапросов. Мы выполним двоичный поиск и узнаем, между какими подзапросами втиснулось это значение, пусть, между m и m+1. Замечательно, значит запросы от 1 до m можно не выполнять, там заведомо ничего нет.

Если значение принадлежит нашему экстенту, значит оно попадает в один из наших диапазонов и мы также можем найти в какой, пусть тоже m. Как и раньше, запросы с номерами 1 … m-1 можно не выполнять. А вот интервал с номером m заслуживает отдельного запроса, который выдаст нам всё, что в нем расположено.

Что же, возможно, еще остались данные, продолжим. Опять выполним запрос, но уже не от угла до угла, а от начала интервала m+1 до правого верхнего угла. И будем так поступать, пока не дойдем до конца списка интервалов.

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

Самое время проверить идею на практике. В качестве тестовой площадки станем использовать PostgeSQL 9.6.1, GiST. Замеры проводились на скромной виртуальной машине с двумя ядрами и 4 Гб ОЗУ, поэтому времена не имеют абсолютной ценности, а вот числам прочитанных страниц можно доверять.

Исходные данные


В качестве данных использованы 100 млн случайных точек в экстенте [0...1 000 000][0...1 000 000].

Заведем таблицу для 2-мерных точечных данных:

create table test_points (x integer,y integer);

Данные создадим:

gawk скриптом
BEGIN{
for (i = 0; i < 100000000; i++)
{
x = int(1000000 * rand());
y = int(1000000 * rand());
print x"\t«z;
}
exit(0);
}

Отсортируем полученный файл (объяснение далее) и оператором COPY зальем его в таблицу:

COPY test_points from '/home/.../data.csv';

Заливка таблицы занимает несколько минут. Размер данных (\dt+) — 4358 Mb

R-tree


Соответствующий индекс создается командой:

create index test_gist_idx on test_points using gist ((point(x,y)));

Но есть нюанс. На случайных данных индекс строится очень долго (во всяком случае у автора он не успел построиться за ночь). Построение на заранее отсортированных данных заняло около часа.

Размер индекса (\di+) — 9031 Mb

В сущности, для нас порядок данных в таблице не важен, но он должен быть общим для разных методов, поэтому пришлось использовать отсортированную таблицу.

Тестовый запрос выглядит так:

select count(1) from test_points where 
    point(x,y) <@ box(point("xmin","ymin"),point("xmax","ymax"));

Проверка на обычных индексах


Для проверки работоспособности будем выполнять пространственные запросы и на отдельных индексах по x & y. Изготавливаются они так:

create index x_test_points on test_points (x);
create index y_test_points on test_points (y);

Занимает это несколько минут.

Тестовый запрос:

select count(1) from test_points where 
      x >= "xmin" and x <= "xmax" and y >= "ymin" and y <= "ymax";

Z-индекс


Теперь нам потребуется функция, умеющая преобразовывать x,y координаты в Z-значение.

Для начала создадим расширение, в нём функцию:

CREATE FUNCTION zcurve_val_from_xy(bigint, bigint) RETURNS bigint
AS 'MODULE_PATHNAME'
LANGUAGE C IMMUTABLE STRICT;

Её тело:
static uint32 stoBits[8] = {0x0001, 0x0002, 0x0004, 0x0008, 
                  0x0010, 0x0020, 0x0040, 0x0080};
uint64
zcurve_fromXY (uint32 ix, uint32 iy)
{
  uint64 val = 0;
  int curmask = 0xf;
  unsigned char *ptr = (unsigned char *)&val;
  int i;
  for (i = 0; i < 8; i++)
    {
      int xp = (ix & curmask) >> (i<<2);
      int yp = (iy & curmask) >> (i<<2);
      int tmp = (xp & stoBits[0]) | ((yp & stoBits[0])<<1) |
            ((xp & stoBits[1])<<1) | ((yp & stoBits[1])<<2) |
            ((xp & stoBits[2])<<2) | ((yp & stoBits[2])<<3) |
            ((xp & stoBits[3])<<3) | ((yp & stoBits[3])<<4);
      curmask <<= 4;
      ptr[i] = (unsigned char)tmp;
    }
  return val;
}
Datum
zcurve_val_from_xy(PG_FUNCTION_ARGS)
{
   uint64 v1  = PG_GETARG_INT64(0);
   uint64 v2  = PG_GETARG_INT64(1);
   PG_RETURN_INT64(zcurve_fromXY(v1, v2));
}


Теперь (после CREATE EXTENSION, конечно) Z-индекс строится следующим образом:

create index zcurve_test_points on test_points(zcurve_val_from_xy(x, y));

Это занимает несколько минут (и не требует сортировки данных).

Размер индекса (\di+) — 2142 Mb (в ~4 раза меньше чем у R-tree)

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


Итак, в нашем первом (назовем его “наивным”) варианте мы будем поступать так:

  1. Для экстента размером dx*dy заводим массив идентификаторов соответствующего размера
  2. Для каждой точки в экстенте вычисляес Z-значение
  3. Сортируем массив идентификаторов
  4. Находим непрерывные интервалы
  5. Для каждого интервала выполняем подзапрос вида:

    select * from test_points where zcurve_val_from_xy(x, y) between $1 and $2

  6. Получаем результат

Для поиска с помощью этого варианта будем использовать функцию (тело ниже):

CREATE FUNCTION zcurve_oids_by_extent(bigint, bigint, bigint, bigint)
RETURNS bigint
AS 'MODULE_PATHNAME'
LANGUAGE C IMMUTABLE STRICT;

Пространственный запрос с использованием этой функции выглядит так:

select zcurve_oids_by_extent("xmin","ymin","xmax","ymax") as count;

Функция возвращает только число попаданий, сами данные при необходимости могут быть выведены с помощью “elog(INFO…)”.

Второй, улучшенный, (назовём его “с пробами”) вариант выглядит следующим образом:

  1. для экстента размером dx*dy заводим массив идентификаторов соответствующего размера
  2. для каждой точки в экстенте вычисляес Z-значение
  3. сортируем массив идентификаторов
  4. находим непрерывные интервалы
  5. стартуем с первого найденного интервала:

    1. Выполняем “пробный” запрос типа (параметры — границы интервала):

      select * from test_points where 
          zcurve_val_from_xy(x, y) between $1 and $2 
          order by zcurve_val_from_xy(x, y) limit 1
      

    2. Этот запрос выдает нам строку таблицы с наименьшим Z-значением от начала текущего тестируемого интервала до конца поискового экстента.
    3. Если запрос ничего не нашел, значит и данных в поисковом экстенте не осталось, выходим.
    4. Теперь мы можем проанализировать найденное Z-значение:

      1. Смотрим, попало ли оно в какой — либо из наших интервалов.
      2. Если нет, находим номер следующего оставшегося интервала и переходим на пункт 5.
      3. Если попало, выполняем запрос для этого интервала типа:

        select * from test_points 
            where zcurve_val_from_xy(x, y) between $1 and $2
        

      4. Берем следующий интервал и переходим на пункт 5.

Для поиска с помощью этого варианта будем использовать функцию:

CREATE FUNCTION zcurve_oids_by_extent_ii(bigint, bigint, bigint, bigint)
RETURNS bigint
AS 'MODULE_PATHNAME'
LANGUAGE C IMMUTABLE STRICT;

Пространственный запрос с использованием этой функции выглядит так:

select zcurve_oids_by_extent_ii("xmin","ymin","xmax","ymax") as count;

Функция также возвращает только число попаданий.

Растеризация


В описанных алгоритмах использован весьма простой и неэффективный алгоритм “растеризации” — получения списка интервалов.

С другой стороны, нетрудно выполнить замеры средних времен его работы на случайных экстентах нужного размера. Вот они:
Экстент dx * dy Ожидаемое среднее число точек в выдаче Время, мсек
100X100 1 .96 (.37 + .59)
316X316 10 11 (3.9 + 7.1)
1000X1000 100 119.7 (35 + 84.7)
3162X3162 1000 1298 (388 + 910)
10000X10000 10 000 14696 (3883 + 10813)
В скобках отдельно указаны две фазы — вычисление Z-значений и сортировка

Результаты


Вот сводная таблица с данными.
Npoints Type Time(ms) Reads Shared Hits
1 X&Y
rtree
Z-value
Z-value-ii
43.6
.5
8.3(9.4)
1.1(2.2)
59.0173
4.2314
4.0988
4.1984 (12.623)
6.1596
2.6565
803.171
20.1893 (57.775)
10 X&Y
rtree
Z-value
Z-value-ii
83.5
.6
15(26)
4(15)
182.592
13.7341
14.834
14.832 (31.439)
9.24363
2.72466
2527.56
61.814 (186.314)
100 X&Y
rtree
Z-value
Z-value-ii
220
2.1
80(200)
10(130)
704.208
95.8179
95.215
96.198 (160.443)
16.528
5.3754
8007.3
208.214 (600.049)
1000 X&Y
rtree
Z-value
Z-value-ii
740
12
500(1800)
200(1500)
3176.06
746.617
739.32
739.58 (912.631)
55.135
25.439
25816
842.88 (2028.81)
10 000 X&Y
rtree
Z-value
Z-value-ii
2 500
70 … 1 200
4700(19000)
1300(16000)
12393.2
4385.64
4284.45
4305.78 (4669)
101.43
121.56
86274.9
5785.06 (9188)
Npoints — среднее число точек в выдаче.

Type -
  • ‘X&Y’ — использование отдельных индексов по x & y
  • ‘rtree’ — запрос через R-tree
  • Z-value — честный поиск через интервалы
  • ‘Z-value-ii’ — поиск через интервалы с пробами

Time(ms) — среднее время выполнения запроса. В данных условиях величина весьма нестабильная, зависит и от кэша СУБД и от дискового кэша виртуальной машины и от дискового кэша хост-системы. Здесь приведена скорее для справки. Для Z-value и Z-value-ii приведены 2 числа. В скобках — фактическое время. Без скобок — время за вычетом затрат на “растеризацию”.

Reads — среднее число чтений на запрос (получено через EXPLAIN (ANALYZE,BUFFERS))

Shared hits — число обращений к буферам (...)
Для Z-value-ii в колонках Reads & Shared hits приведены 2 числа. В скобках — полное число чтений. Без скобок — за вычетом зондирующих запросов с order by и limit 1. Сделано это из за мнения, что такой запрос вычитывает все данные в интервале, сортирует и отдает минимальное, вместо того, чтобы просто отдать 1 значение из уже отсортированного индекса. Поэтому статистика по таким запросам была сочтена излишней, но приведена для справки.

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

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

Выводы


  1. R-дерево всё же очень хорошо в статике, КПД чтений страниц весьма высок.
  2. А вот Z-order индексу иногда приходится читать страницы, на которых нет нужных данных. Это происходит, когда тестирующий курсор попадает между интервалов, есть вероятность, что в этом промежутке окажется много чужих точек и конкретная страница не содержит никаких итоговых данных.
  3. И тем не менее, за счет более плотной упаковки Z-order близок к R-дереву по числу реально прочитанных страниц. Это говорит о том, что потенциально Z-order способен выдавать близкую производительность.
  4. Z-order индекс читает огромное количество страниц из кэша т.к. многократно делаются запросы по одному и тому же месту. С другой стороны, эти чтения стоят относительно недорого.
  5. На больших запросах Z-order сильно проигрывает в скорости R-дереву. Объясняется это тем, что для выполнения подзапросов мы пользуемся SPI — высокоуровневым и не слишком быстрым механизмом. И с “растеризацией”, конечно, надо что-то делать.
  6. На первый взгляд, использование проб интервалов не сильно то и ускорило работу, а формально даже ухудшило статистику чтений страниц. Но надо понимать, что это издержки высокоуровневых средств, которыми пришлось воспользоваться. Потенциально индекс на основе Z-order не хуже R-дерева по производительности и сильно лучше по остальным тактико-техническим характеристикам.

Перспективы


Для создания на основе Z-order кривой полноценного пространственного индекса, который был бы способен на равных соревноваться с R-деревом, необходимо решить следующие проблемы:

  • придумать недорогой алгоритм получения списка подинтервалов по экстенту
  • перейти на низкоуровневый доступ к дереву индекса

К счастью, и то и другое не выглядит чем-то невозможным.

Исходники
#include "postgres.h"
#include "catalog/pg_type.h"
#include "fmgr.h"
#include <string.h>
#include "executor/spi.h"

PG_MODULE_MAGIC;

uint64 zcurve_fromXY (uint32 ix, uint32 iy);
void zcurve_toXY (uint64 al, uint32 *px, uint32 *py);

static uint32 stoBits[8] = {0x0001, 0x0002, 0x0004, 0x0008, 
                  0x0010, 0x0020, 0x0040, 0x0080};
uint64
zcurve_fromXY (uint32 ix, uint32 iy)
{
  uint64 val = 0;
  int curmask = 0xf;
  unsigned char *ptr = (unsigned char *)&val;
  int i;
  for (i = 0; i < 8; i++)
    {
      int xp = (ix & curmask) >> (i<<2);
      int yp = (iy & curmask) >> (i<<2);
      int tmp = (xp & stoBits[0]) | ((yp & stoBits[0])<<1) |
            ((xp & stoBits[1])<<1) | ((yp & stoBits[1])<<2) |
            ((xp & stoBits[2])<<2) | ((yp & stoBits[2])<<3) |
            ((xp & stoBits[3])<<3) | ((yp & stoBits[3])<<4);
      curmask <<= 4;
      ptr[i] = (unsigned char)tmp;
    }
  return val;
}

void 
zcurve_toXY (uint64 al, uint32 *px, uint32 *py)
{
  unsigned char *ptr = (unsigned char *)&al;
  int ix = 0;
  int iy = 0;
  int i;

  if (!px || !py)
    return;

  for (i = 0; i < 8; i++)
    {
      int tmp = ptr[i];
      int tmpx = (tmp & stoBits[0]) + ((tmp & stoBits[2])>>1) + ((tmp & stoBits[4])>>2) + ((tmp & stoBits[6])>>3);
      int tmpy = ((tmp & stoBits[1])>>1) + ((tmp & stoBits[3])>>2) + ((tmp & stoBits[5])>>3) + ((tmp & stoBits[7])>>4);
      ix |= tmpx << (i << 2);
      iy |= tmpy << (i << 2);
    }        
  *px = ix;
  *py = iy;
}

PG_FUNCTION_INFO_V1(zcurve_val_from_xy);

Datum
zcurve_val_from_xy(PG_FUNCTION_ARGS)
{
   uint64 v1  = PG_GETARG_INT64(0);
   uint64 v2  = PG_GETARG_INT64(1);

   PG_RETURN_INT64(zcurve_fromXY(v1, v2));
}


static const int s_maxx = 1000000;
static const int s_maxy = 1000000;

#ifndef MIN
#define MIN(a,b) ((a)<(b)?(a):(b))
#endif


static int compare_uint64( const void *arg1, const void *arg2 )
{
  const uint64 *a = (const uint64 *)arg1;
  const uint64 *b = (const uint64 *)arg2;
  if (*a == *b)
    return 0;
  return *a > *b ? 1: -1;
}

SPIPlanPtr prep_interval_request();
int fin_interval_request(SPIPlanPtr pplan);
int run_interval_request(SPIPlanPtr pplan, uint64 v0, uint64 v1);

SPIPlanPtr prep_interval_request()
{
	SPIPlanPtr	pplan;
	char		sql[8192];
	int 		nkeys = 2;
	Oid		argtypes[2] = {INT8OID, INT8OID};	/* key types to prepare execution plan */
	int ret =0;

   	if ((ret = SPI_connect()) < 0)
   	/* internal error */
     	  elog(ERROR, "check_primary_key: SPI_connect returned %d", ret);

	snprintf(sql, sizeof(sql), "select * from test_points where zcurve_val_from_xy(x, y) between $1 and $2");

	/* Prepare plan for query */
	pplan = SPI_prepare(sql, nkeys, argtypes);
	if (pplan == NULL)
	/* internal error */
		elog(ERROR, "check_primary_key: SPI_prepare returned %d", SPI_result);
	return pplan;
}

int fin_interval_request(SPIPlanPtr pplan)
{
   SPI_finish();
   return 0;
}


int 
run_interval_request(SPIPlanPtr pplan, uint64 v0, uint64 v1)
{
	Datum		values[2];	/* key types to prepare execution plan */
	Portal 		portal;
	int cnt = 0, i;

	values[0] = Int64GetDatum(v0);
	values[1] = Int64GetDatum(v1);

	portal = SPI_cursor_open(NULL, pplan, values, NULL, true);
	if (NULL == portal)
		/* internal error */
		elog(ERROR, "check_primary_key: SPI_cursor_open");
	for (;;)
	{
		SPI_cursor_fetch(portal, true, 8);
		if (0 == SPI_processed || NULL == SPI_tuptable)
			break;
		{
			TupleDesc tupdesc = SPI_tuptable->tupdesc;
			for (i = 0; i < SPI_processed; i++)
			{
				HeapTuple tuple = SPI_tuptable->vals[i];
				//elog(INFO, "%s, %s", SPI_getvalue(tuple, tupdesc, 1), SPI_getvalue(tuple, tupdesc, 2));
				cnt++;
			}
		}
	}
	SPI_cursor_close(portal);

	return cnt;
}

PG_FUNCTION_INFO_V1(zcurve_oids_by_extent);
Datum
zcurve_oids_by_extent(PG_FUNCTION_ARGS)
{
   SPIPlanPtr	pplan;

   uint64 x0  = PG_GETARG_INT64(0);
   uint64 y0  = PG_GETARG_INT64(1);
   uint64 x1  = PG_GETARG_INT64(2);
   uint64 y1  = PG_GETARG_INT64(3);
   uint64 *ids = NULL;
   int cnt = 0;
   int sz = 0, ix, iy;

   x0 = MIN(x0, s_maxx);
   y0 = MIN(y0, s_maxy);
   x1 = MIN(x1, s_maxx);
   y1 = MIN(y1, s_maxy);

   if (x0 > x1)
     elog(ERROR, "xmin > xmax");
   if (y0 > y1)
     elog(ERROR, "ymin > ymax");

   sz = (x1 - x0 + 1) * (y1 - y0 + 1);
   ids = (uint64*)palloc(sz * sizeof(uint64));
   if (NULL == ids)
   /* internal error */
     elog(ERROR, "cant alloc %d bytes in zcurve_oids_by_extent", sz);
   for (ix = x0; ix <= x1; ix++)
    for (iy = y0; iy <= y1; iy++)
    {
     	ids[cnt++] = zcurve_fromXY(ix, iy);
    }
   qsort (ids, sz, sizeof(*ids), compare_uint64);

   cnt = 0;

   pplan = prep_interval_request();
   {
//	FILE *fl = fopen("/tmp/ttt.sql", "w");
	int cur_start = 0; 
	int ix;
	for (ix = cur_start + 1; ix < sz; ix++)
	{
		if (ids[ix] != ids[ix - 1] + 1)
		{
			cnt += run_interval_request(pplan, ids[cur_start], ids[ix - 1]);
//			fprintf(fl, "EXPLAIN (ANALYZE,BUFFERS) select * from test_points where zcurve_val_from_xy(x, y) between %ld and %ld;\n", ids[cur_start], ids[ix - 1]);
//			elog(INFO, "%d -> %d (%ld -> %ld)", cur_start, ix - 1, ids[cur_start], ids[ix - 1]);
//			cnt++;
			cur_start = ix;
		}
	}
    	if (cur_start != ix)
	{
			cnt += run_interval_request(pplan, ids[cur_start], ids[ix - 1]);
//			fprintf(fl, "EXPLAIN (ANALYZE,BUFFERS) select * from test_points where zcurve_val_from_xy(x, y) between %ld and %ld;\n", ids[cur_start], ids[ix - 1]);
//			elog(INFO, "%d -> %d (%ld -> %ld)", cur_start, ix - 1, ids[cur_start], ids[ix - 1]);
	}
//	fclose(fl);
   }
   fin_interval_request(pplan);
   pfree(ids);

   PG_RETURN_INT64(cnt);
}

//------------------------------------------------------------------------------------------------

struct interval_ctx_s {
	SPIPlanPtr	cr_;
	SPIPlanPtr	probe_cr_;
	uint64 		cur_val_;
	uint64 		top_val_;
	FILE * 		fl_;
};
typedef struct interval_ctx_s interval_ctx_t;

int prep_interval_request_ii(interval_ctx_t *ctx);
int run_interval_request_ii(interval_ctx_t *ctx, uint64 v0, uint64 v1);
int probe_interval_request_ii(interval_ctx_t *ctx, uint64 v0);
int fin_interval_request_ii(interval_ctx_t *ctx);


int prep_interval_request_ii(interval_ctx_t *ctx)
{
	char		sql[8192];
	int 		nkeys = 2;
	Oid		argtypes[2] = {INT8OID, INT8OID};	/* key types to prepare execution plan */
	int ret =0;

   	if ((ret = SPI_connect()) < 0)
   	/* internal error */
     	  elog(ERROR, "check_primary_key: SPI_connect returned %d", ret);

	snprintf(sql, sizeof(sql), "select * from test_points where zcurve_val_from_xy(x, y) between $1 and $2");
	ctx->cr_ = SPI_prepare(sql, nkeys, argtypes);
	if (ctx->cr_ == NULL)
	/* internal error */
		elog(ERROR, "check_primary_key: SPI_prepare returned %d", SPI_result);

	snprintf(sql, sizeof(sql), "select * from test_points where zcurve_val_from_xy(x, y) between $1 and %ld order by zcurve_val_from_xy(x::int4, y::int4) limit 1", ctx->top_val_);
	ctx->probe_cr_ = SPI_prepare(sql, 1, argtypes);
	if (ctx->probe_cr_ == NULL)
	/* internal error */
		elog(ERROR, "check_primary_key: SPI_prepare returned %d", SPI_result);

	return 1;
}

int 
probe_interval_request_ii(interval_ctx_t *ctx, uint64 v0)
{
	Datum		values[1];	/* key types to prepare execution plan */
	Portal 		portal;

	values[0] = Int64GetDatum(v0);
	{
//		uint32 lx, ly;
//		zcurve_toXY (v0, &lx, &ly);
//
//		elog(INFO, "probe(%ld:%d,%d)", v0, lx, ly);
	}

	if (ctx->fl_)
		fprintf(ctx->fl_, "EXPLAIN (ANALYZE,BUFFERS) select * from test_points where zcurve_val_from_xy(x, y) between %ld and %ld order by zcurve_val_from_xy(x::int4, y::int4) limit 1;\n", v0, ctx->top_val_);
	portal = SPI_cursor_open(NULL, ctx->probe_cr_, values, NULL, true);
	if (NULL == portal)
		/* internal error */
		elog(ERROR, "check_primary_key: SPI_cursor_open");
	{
		SPI_cursor_fetch(portal, true, 1);
		if (0 != SPI_processed && NULL != SPI_tuptable)
		{
			TupleDesc tupdesc = SPI_tuptable->tupdesc;

			bool isnull;
			HeapTuple tuple = SPI_tuptable->vals[0];
			Datum dx, dy;
			uint64 zv = 0;
			dx = SPI_getbinval(tuple, tupdesc, 1, &isnull);
			dy = SPI_getbinval(tuple, tupdesc, 2, &isnull);
			zv = zcurve_fromXY(DatumGetInt64(dx), DatumGetInt64(dy));
//			elog(INFO, "%ld %ld -> %ld", DatumGetInt64(dx), DatumGetInt64(dy), zv);

			ctx->cur_val_ = zv;
			SPI_cursor_close(portal);
			return 1;
		}
		SPI_cursor_close(portal);
	}
	return 0;
}


int 
run_interval_request_ii(interval_ctx_t *ctx, uint64 v0, uint64 v1)
{
	Datum		values[2];	/* key types to prepare execution plan */
	Portal 		portal;
	int cnt = 0, i;

	values[0] = Int64GetDatum(v0);
	values[1] = Int64GetDatum(v1);
//	elog(INFO, "[%ld %ld]", v0, v1);

	if (ctx->fl_)
		fprintf(ctx->fl_, "EXPLAIN (ANALYZE,BUFFERS) select * from test_points where zcurve_val_from_xy(x, y) between %ld and %ld;\n", v0, v1);
	portal = SPI_cursor_open(NULL, ctx->cr_, values, NULL, true);
	if (NULL == portal)
		/* internal error */
		elog(ERROR, "check_primary_key: SPI_cursor_open");
	for (;;)
	{
		SPI_cursor_fetch(portal, true, 8);
		if (0 == SPI_processed || NULL == SPI_tuptable)
			break;
		{
			TupleDesc tupdesc = SPI_tuptable->tupdesc;
			for (i = 0; i < SPI_processed; i++)
			{
				HeapTuple tuple = SPI_tuptable->vals[i];
//				elog(INFO, "%s, %s", SPI_getvalue(tuple, tupdesc, 1), SPI_getvalue(tuple, tupdesc, 2));
				cnt++;
			}
		}
	}
	SPI_cursor_close(portal);

	return cnt;
}


PG_FUNCTION_INFO_V1(zcurve_oids_by_extent_ii);
Datum
zcurve_oids_by_extent_ii(PG_FUNCTION_ARGS)
{
   uint64 x0  = PG_GETARG_INT64(0);
   uint64 y0  = PG_GETARG_INT64(1);
   uint64 x1  = PG_GETARG_INT64(2);
   uint64 y1  = PG_GETARG_INT64(3);
   uint64 *ids = NULL;
   int cnt = 0;
   int sz = 0, ix, iy;
   interval_ctx_t ctx;

   x0 = MIN(x0, s_maxx);
   y0 = MIN(y0, s_maxy);
   x1 = MIN(x1, s_maxx);
   y1 = MIN(y1, s_maxy);

   if (x0 > x1)
     elog(ERROR, "xmin > xmax");
   if (y0 > y1)
     elog(ERROR, "ymin > ymax");

   sz = (x1 - x0 + 1) * (y1 - y0 + 1);
   ids = (uint64*)palloc(sz * sizeof(uint64));
   if (NULL == ids)
   /* internal error */
     elog(ERROR, "can't alloc %d bytes in zcurve_oids_by_extent_ii", sz);
   for (ix = x0; ix <= x1; ix++)
     for (iy = y0; iy <= y1; iy++)
     {
     	ids[cnt++] = zcurve_fromXY(ix, iy);
     }
   qsort (ids, sz, sizeof(*ids), compare_uint64);

   ctx.top_val_ = ids[sz - 1];
   ctx.cur_val_ = 0;
   ctx.cr_ = NULL;
   ctx.probe_cr_ = NULL;
   ctx.fl_ = NULL;//fopen("/tmp/ttt.sql", "w");
   
   cnt = 0;

   prep_interval_request_ii(&ctx);
   {
	int cur_start = 0; 
	int ix;
	for (ix = cur_start + 1; ix < sz; ix++)
	{
		if (0 == probe_interval_request_ii(&ctx, ids[cur_start]))
			break;
		for (; cur_start < sz && ids[cur_start] < ctx.cur_val_; cur_start++);
//		if (ctx.cur_val_ != ids[cur_start])
//		{
//			cur_start++;
//			continue;
//		}
		ix = cur_start + 1;
                if (ix >= sz)
			break;
		for (; ix < sz && ids[ix] == ids[ix - 1] + 1; ix++);
//elog(INFO, "%d %d %d", ix, cur_start, sz);
		cnt += run_interval_request_ii(&ctx, ids[cur_start], ids[ix - 1]);
		cur_start = ix;
	}
   }
   if (ctx.fl_)
   	fclose(ctx.fl_);
   fin_interval_request(NULL);
   pfree(ids);

   PG_RETURN_INT64(cnt);
}

Выкладываю здесь, а не выложены на github т.к. код сугубо экспериментальный и не имеет практической ценности.

PPS: Огромное спасибо ребятам из PostgresPro за то что сподвигли меня на эту работу.
Поделиться с друзьями
-->

Комментарии (7)


  1. neolink
    10.01.2017 02:32
    +1

    Для полноты картины таки стоит проверить решение в лоб с составным индексом

    create index test_points_x_y on test_points (x, y);
    


    1. zzeng
      10.01.2017 10:49
      +1

      глянул одним глазком,
      составной индекс занимает 2 Гб,
      запрос select count(1)… отрабатывается как index only scan, т.е. без обращения к таблице, читает страниц от 2 до 5 раз меньше по сравнению с X&Y.
      Но это скорее хороший пример качественной работы оптимизатора, чем заслуга собственно индекса.


  1. kashey
    10.01.2017 12:47
    +2

    Очень хорошо, но есть пара проблем:
    1. Z-code это, фактически linear quad-tree. Не R. И совсем не B. Все, что вы обьявили бонусами в начале статьи — не работает. Никакой балансировки нет. И дерева никакого нет. Это что-то среднее между фракталом и кривой – spatial filling curve, И их много бывает разных. Смысл один — понижают размерность.
    2. Попробуйте Hilbert. Считать его сложнее, но более «близкие» значения в нем расположены «ближе»
    3. «Секрет» работы Z-value-ii хорошо описан в интернетах как BIGMIN/LITMAX (https://en.wikipedia.org/wiki/Z-order_curve#Use_with_one-dimensional_data_structures_for_range_searching)

    В общем крайне рекомендую активно дружить с этой математикой. Она реально хорошо работает, но надо понимать ограничения.
    И не забывать — Z это точка. В то время как R-tree продолжает жрать кактус.

    PS: Еще в интернетах можно найти множество более быстрых функций построения Z-кодов
    PPS: Некоторые модификации R-tree сами по Z строятся
    PPPS: Вэлкам ту зе клаб!


    1. zzeng
      10.01.2017 13:31
      +1

      по П.1 Z-code — это число, а числа в данном случае мы укладываем в B-дерево.
      То, что вы (по видимому) имели ввиду, называется k-d-деревья и это совсем другая история.
      И не забывать — Z это точка, прямоугольник — четырехмерная Z точка.


      1. kashey
        10.01.2017 13:59
        +2

        1. В принципе да, на стороне БД это уложиться в B-дерево. Получается достаточно честно. Смотрите на это без БД.
        2. Кстати – Z-code, который иногда еще называют «bit interleaving», в принципе KD. Каждый бит режет пространство своей плоскостью.
        3. Только вот поиск в окрестности такой четырех-мерной точки будет немного затруднен :)

        В общем случае Z(и компания) просто позволяет в начале как-то уложить в некий интервал данные, а потом в нем искать.
        При это самый «быстрый» поиск по кодам, когда вам нужна конкретная «ячейка» виртуального quadtree, лучше всего (технически) работает через маску. Вы говорите с чего код должен начинаться – фиксируете первые шаги погружения в дерево, а после некого момента «отпускаете» код, позволяя ему свободно менятся.
        Если весь код 4 бита, те 2 уровня, то фиксация 0b1000 с маской 0b0111 матчится на все элементы в левом(или правом) поддереве, те между 0b1000 и 0b1111.
        В принципе — это — возможность использовать эти «числа» как NestedSet — одна из основых фичей этих кривых. За прошедшие годы для разных spatial filling придумали очень много наркомании. Текстуры в памяти видеокарты по ним разложены. Gray Code, на котором 99% электроники работает — тоже spatial filling.

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


        1. zzeng
          10.01.2017 14:09
          +2

          Спасибо за ссылку, в принципе, я в курсе некоторой части этой магии.
          И эту статью рассматривал как вводную для описания алгоритма поиска именно по B-дереву, как самому честному способу хранить одномерные данные в (по большому счету) одномерном дисковом пространстве.
          Через недельку — другую, добро пожаловать.
          А одна старая статейка найдётся не только у вас :)


          1. kashey
            10.01.2017 14:21
            +1

            О! Месье знает толк колбасных обрезках. Жму вашу руку, коллега.