Отображение N-мерного значения в точку на кривой Гильберта
У меня огромный набор N-мерных точек (десятки миллионов; N близко к 100).
Мне нужно отобразить эти точки в одном измерении, сохраняя пространственную локализацию. Я хочу использовать кривую заполнения пространства Гильберта, чтобы сделать это.
Для каждой точки я хочу выбрать ближайшую точку на кривой. Значение Гильберта точки (длина кривой от начала кривой до выбранной точки) - это значение одного измерения, которое я ищу.
Вычисления не должны быть мгновенными, но я ожидаю, что они будут не более нескольких часов на приличном оборудовании современного домашнего ПК.
Есть предложения по внедрению? Есть ли библиотеки, которые бы мне помогли? (Язык не имеет большого значения.)
8 ответов
Я наконец сломался и выложил немного денег. В AIP (Американский институт физики) есть хорошая короткая статья с исходным кодом на языке C. "Программирование кривой Гильберта" Джона Скиллинга (из AIP Conf. Proc. 707, 381 (2004)), есть приложение с кодом для отображения в обоих направлениях. Он работает для любого числа измерений> 1, не является рекурсивным, не использует таблицы поиска с изменением состояния, которые поглощают огромные объемы памяти, и в основном использует битовые операции. Таким образом, это достаточно быстро и имеет хороший объем памяти.
Если вы решили приобрести статью, я обнаружил ошибку в исходном коде.
Следующая строка кода (находится в функции TransposetoAxes) содержит ошибку:
для ( i = n-1; i >= 0; i--) X[i] ^= X[i-1];
Коррекция заключается в изменении значения больше или равно (>=) на значение больше (>). Без этого исправления доступ к массиву X осуществляется с использованием отрицательного индекса, когда переменная "i" становится равной нулю, что приводит к сбою программы.
Я рекомендую прочитать статью (которая состоит из семи страниц, включая код), так как она объясняет, как работает алгоритм, что далеко не очевидно.
Я перевел его код на C# для собственного использования. Код следует. Skilling выполняет преобразование на месте, перезаписывая вектор, который вы передаете. Я решил сделать клон входного вектора и вернуть новую копию. Также я реализовал методы как методы расширения.
Код Скиллинга представляет индекс Гильберта в виде транспонирования, хранящегося в виде массива. Мне удобнее чередовать биты и формировать один BigInteger (более полезный в словарях, легче перебирать циклы и т. Д.), Но я оптимизировал эту операцию и ее инверсию с помощью магических чисел, битовых операций и тому подобного, а также код длинный, поэтому я его опустил.
namespace HilbertExtensions
{
/// <summary>
/// Convert between Hilbert index and N-dimensional points.
///
/// The Hilbert index is expressed as an array of transposed bits.
///
/// Example: 5 bits for each of n=3 coordinates.
/// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
/// as its Transpose ^
/// X[0] = A D G J M X[2]| 7
/// X[1] = B E H K N <-------> | /X[1]
/// X[2] = C F I L O axes |/
/// high low 0------> X[0]
///
/// NOTE: This algorithm is derived from work done by John Skilling and published in "Programming the Hilbert curve".
/// (c) 2004 American Institute of Physics.
///
/// </summary>
public static class HilbertCurveTransform
{
/// <summary>
/// Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
///
/// Note: In Skilling's paper, this function is named TransposetoAxes.
/// </summary>
/// <param name="transposedIndex">The Hilbert index stored in transposed form.</param>
/// <param name="bits">Number of bits per coordinate.</param>
/// <returns>Coordinate vector.</returns>
public static uint[] HilbertAxes(this uint[] transposedIndex, int bits)
{
var X = (uint[])transposedIndex.Clone();
int n = X.Length; // n: Number of dimensions
uint N = 2U << (bits - 1), P, Q, t;
int i;
// Gray decode by H ^ (H/2)
t = X[n - 1] >> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (i = n - 1; i > 0; i--)
X[i] ^= X[i - 1];
X[0] ^= t;
// Undo excess work
for (Q = 2; Q != N; Q <<= 1)
{
P = Q - 1;
for (i = n - 1; i >= 0; i--)
if ((X[i] & Q) != 0U)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
return X;
}
/// <summary>
/// Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
/// That distance will be transposed; broken into pieces and distributed into an array.
///
/// The number of dimensions is the length of the hilbertAxes array.
///
/// Note: In Skilling's paper, this function is called AxestoTranspose.
/// </summary>
/// <param name="hilbertAxes">Point in N-space.</param>
/// <param name="bits">Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.</param>
/// <returns>The Hilbert distance (or index) as a transposed Hilbert index.</returns>
public static uint[] HilbertIndexTransposed(this uint[] hilbertAxes, int bits)
{
var X = (uint[])hilbertAxes.Clone();
var n = hilbertAxes.Length; // n: Number of dimensions
uint M = 1U << (bits - 1), P, Q, t;
int i;
// Inverse undo
for (Q = M; Q > 1; Q >>= 1)
{
P = Q - 1;
for (i = 0; i < n; i++)
if ((X[i] & Q) != 0)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
// Gray encode
for (i = 1; i < n; i++)
X[i] ^= X[i - 1];
t = 0;
for (Q = M; Q > 1; Q >>= 1)
if ((X[n - 1] & Q)!=0)
t ^= Q - 1;
for (i = 0; i < n; i++)
X[i] ^= t;
return X;
}
}
}
Я разместил рабочий код в C# на github.
Смотрите https://github.com/paulchernoch/HilbertTransformation
Алгоритм для отображения из n->1 и 1->n приведен здесь "Расчет отображений между одномерными и n-мерными значениями по кривой заполнения Гильберта" Дж. К. Лоудер
Если вы используете Google для "SFC module and Kademlia overlay", вы найдете группу, которая утверждает, что использует его как часть своей системы. Если вы просматриваете источник, вы можете извлечь соответствующую функцию.
Я потратил немного времени на перевод кода Пола Черноха на Java и его очистку. Возможно, в моем коде есть ошибка, особенно потому, что у меня нет доступа к статье, из которой она изначально. Тем не менее, он проходит, какие юнит-тесты мне удалось написать. Это ниже.
Обратите внимание, что я оценил кривые Z-порядка и Гильберта для пространственной индексации на больших наборах данных. Я должен сказать, что Z-Order обеспечивает гораздо лучшее качество. Но не стесняйтесь попробовать для себя.
/**
* Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
*
* Note: In Skilling's paper, this function is named TransposetoAxes.
* @param transposedIndex The Hilbert index stored in transposed form.
* @param bits Number of bits per coordinate.
* @return Point in N-space.
*/
static long[] HilbertAxes(final long[] transposedIndex, final int bits) {
final long[] result = transposedIndex.clone();
final int dims = result.length;
grayDecode(result, dims);
undoExcessWork(result, dims, bits);
return result;
}
static void grayDecode(final long[] result, final int dims) {
final long swap = result[dims - 1] >>> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (int i = dims - 1; i > 0; i--)
result[i] ^= result[i - 1];
result[0] ^= swap;
}
static void undoExcessWork(final long[] result, final int dims, final int bits) {
for (long bit = 2, n = 1; n != bits; bit <<= 1, ++n) {
final long mask = bit - 1;
for (int i = dims - 1; i >= 0; i--)
if ((result[i] & bit) != 0)
result[0] ^= mask; // invert
else
swapBits(result, mask, i);
}
}
/**
* Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
* That distance will be transposed; broken into pieces and distributed into an array.
*
* The number of dimensions is the length of the hilbertAxes array.
*
* Note: In Skilling's paper, this function is called AxestoTranspose.
* @param hilbertAxes Point in N-space.
* @param bits Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.
* @return The Hilbert distance (or index) as a transposed Hilbert index.
*/
static long[] HilbertIndexTransposed(final long[] hilbertAxes, final int bits) {
final long[] result = hilbertAxes.clone();
final int dims = hilbertAxes.length;
final long maxBit = 1L << (bits - 1);
inverseUndo(result, dims, maxBit);
grayEncode(result, dims, maxBit);
return result;
}
static void inverseUndo(final long[] result, final int dims, final long maxBit) {
for (long bit = maxBit; bit != 0; bit >>>= 1) {
final long mask = bit - 1;
for (int i = 0; i < dims; i++)
if ((result[i] & bit) != 0)
result[0] ^= mask; // invert
else
swapBits(result, mask, i);
} // exchange
}
static void grayEncode(final long[] result, final int dims, final long maxBit) {
for (int i = 1; i < dims; i++)
result[i] ^= result[i - 1];
long mask = 0;
for (long bit = maxBit; bit != 0; bit >>>= 1)
if ((result[dims - 1] & bit) != 0)
mask ^= bit - 1;
for (int i = 0; i < dims; i++)
result[i] ^= mask;
}
static void swapBits(final long[] array, final long mask, final int index) {
final long swap = (array[0] ^ array[index]) & mask;
array[0] ^= swap;
array[index] ^= swap;
}
Мне не ясно, как это будет делать то, что вы хотите. Рассмотрим этот случай с трехмерной графикой:
001 ------ 101
|\ |\
| \ | \
| 011 ------ 111
| | | |
| | | |
000 -|---- 100 |
\ | \ |
\ | \ |
010 ------ 110
который может быть "Hilbertized" по следующему пути:
001 -----> 101
\ \
\ \
011 111
^ |
| |
000 | 100 |
\ | \ |
\ | \ V
010 110
в 1D порядке:
000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100
Вот немного неприятно. Рассмотрим список пар и 1D расстояния ниже:
000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1
Во всех случаях левые и правые значения находятся на одном и том же трехмерном расстоянии друг от друга (+/- 1 в первой позиции), что, по-видимому, подразумевает сходную "пространственную локализацию". Но линеаризация любого выбора размерного порядка (y, затем z, затем z, в приведенном выше примере) нарушает эту локальность.
Еще один способ сказать, что если взять начальную точку и упорядочить оставшиеся точки по их расстоянию от этой начальной точки, то результаты будут существенно отличаться. принятие 000 в качестве начала, например:
1D ordering : distance 3D ordering : distance
---------------------- ----------------------
010 : 1 001,010,100 : 1
011,101,110 : sqrt(2)
111 : sqrt(3)
011 : 2
001 : 3
101 : 4
111 : 5
110 : 6
100 : 7
Этот эффект растет экспоненциально с увеличением количества измерений (при условии, что каждое измерение имеет одинаковый "размер").
Вот исходный код C Джона Скиллинга для кодирования/декодирования координат Гильберта в произвольных измерениях. Это из статьи, процитированной Полом Чернохом выше, «Программирование кривой Гильберта» Джона Скиллинга (из AIP Conf. Proc. 707, 381 (2004)).
В этом коде применено исправление ошибки.
Я также расширил main(), чтобы показать как кодирование, так и декодирование. Я также добавил функции interleaveBits() и uninterleaveBits(), которые преобразуют гильбертовские координаты в единый код и обратно, что интересует большинство людей.
Код Скиллинга работает для произвольных измерений, но мои функции interleaveBits() специфичны для трех измерений. Зато легко продлевается.
//+++++++++++++++++++++++++++ PUBLIC-DOMAIN SOFTWARE ++++++++++++++++++++++++++
// Functions: TransposetoAxes AxestoTranspose
// Purpose: Transform in-place between Hilbert transpose and geometrical axes
// Example: b=5 bits for each of n=3 coordinates.
// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
// as its Transpose
// X[0] = A D G J M X[2]|
// X[1] = B E H K N <-------> | /X[1]
// X[2] = C F I L O axes |/
// high low 0------ X[0]
// Axes are stored conventially as b-bit integers.
// Author: John Skilling 20 Apr 2001 to 11 Oct 2003
//-----------------------------------------------------------------------------
#include <cstdio>
typedef unsigned int coord_t; // char,short,int for up to 8,16,32 bits per word
void TransposetoAxes(coord_t* X, int b, int n) // Position, #bits, dimension
{
coord_t N = 2 << (b - 1), P, Q, t;
// Gray decode by H ^ (H/2)
t = X[n - 1] >> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (int i = n - 1; i > 0; i--) X[i] ^= X[i - 1];
X[0] ^= t;
// Undo excess work
for (Q = 2; Q != N; Q <<= 1) {
P = Q - 1;
for (int i = n - 1; i >= 0; i--)
if (X[i] & Q) // Invert
X[0] ^= P;
else { // Exchange
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
}
void AxestoTranspose(coord_t* X, int b, int n) // Position, #bits, dimension
{
coord_t M = 1 << (b - 1), P, Q, t;
// Inverse undo
for (Q = M; Q > 1; Q >>= 1) {
P = Q - 1;
for (int i = 0; i < n; i++)
if (X[i] & Q) // Invert
X[0] ^= P;
else { // Exchange
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
// Gray encode
for (int i = 1; i < n; i++) X[i] ^= X[i - 1];
t = 0;
for (Q = M; Q > 1; Q >>= 1)
if (X[n - 1] & Q) t ^= Q - 1;
for (int i = 0; i < n; i++) X[i] ^= t;
}
int interleaveBits(coord_t* X, int b, int n) // Position, #bits, dimension
{
unsigned int codex = 0, codey = 0, codez = 0;
const int nbits2 = 2 * b;
for (int i = 0, andbit = 1; i < nbits2; i += 2, andbit <<= 1) {
codex |= (unsigned int)(X[0] & andbit) << i;
codey |= (unsigned int)(X[1] & andbit) << i;
codez |= (unsigned int)(X[2] & andbit) << i;
}
return (codex << 2) | (codey << 1) | codez;
}
// From https://github.com/Forceflow/libmorton/blob/main/include/libmorton/morton3D.h
void uninterleaveBits(coord_t* X, int b, int n, unsigned int code) // Position, #bits, dimension
{
X[0] = X[1] = X[2] = 0;
for (unsigned int i = 0; i <= b; ++i) {
unsigned int selector = 1;
unsigned int shift_selector = 3 * i;
unsigned int shiftback = 2 * i;
X[2] |= (code & (selector << shift_selector)) >> (shiftback);
X[1] |= (code & (selector << (shift_selector + 1))) >> (shiftback + 1);
X[0] |= (code & (selector << (shift_selector + 2))) >> (shiftback + 2);
}
}
int main()
{
coord_t X[3] = {5, 10, 20}; // Any position in 32x32x32 cube
printf("Input coords = %d,%d,%d\n", X[0], X[1], X[2]);
AxestoTranspose(X, 5, 3); // Hilbert transpose for 5 bits and 3 dimensions
printf("Hilbert coords = %d,%d,%d\n", X[0], X[1], X[2]);
unsigned int code = interleaveBits(X, 5, 3);
printf("Hilbert integer = %d = %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d = 7865 check\n", code,
X[0] >> 4 & 1, X[1] >> 4 & 1, X[2] >> 4 & 1,
X[0] >> 3 & 1, X[1] >> 3 & 1, X[2] >> 3 & 1,
X[0] >> 2 & 1, X[1] >> 2 & 1, X[2] >> 2 & 1,
X[0] >> 1 & 1, X[1] >> 1 & 1, X[2] >> 1 & 1,
X[0] >> 0 & 1, X[1] >> 0 & 1, X[2] >> 0 & 1);
uninterleaveBits(X, 5, 3, code);
printf("Reconstructed Hilbert coords = %d,%d,%d\n", X[0], X[1], X[2]);
TransposetoAxes(X, 5, 3);
printf("Orig coords = %d,%d,%d\n", X[0], X[1], X[2]);
return 0;
}
EDIT: я взял этот код и соединил его с аналогичным кодом для других заполняющих пространство кривых (Morton, Raster, Boustrophedonic и Tiled) и сделал их доступными на GitHub. Я включил как прямое, так и обратное преобразование для всех кривых, а также некоторый код, который сравнивает их друг с другом для различных показателей качества. См. https://github.com/davemc0/DMcTools/blob/main/Math/SpaceFillCurve.h и код тестирования https://github.com/davemc0/DMcTools/blob/main/Test/SpaceFillCurveTest.cpp .
Другой возможностью было бы построить дерево kd на ваших данных, а затем выполнить упорядоченный обход дерева, чтобы получить порядок. Для построения дерева kd требуется только хороший алгоритм поиска медианы, которых много.
Я не понимаю, как вы можете использовать кривую Гильберта в одном измерении.
Если вы заинтересованы в отображении точек в более низкое измерение при сохранении расстояний (с минимальной ошибкой), вы можете обратиться к алгоритмам "Многомерное масштабирование".
Имитация отжига является одним из подходов.
Изменить: Спасибо за комментарий. Теперь я понимаю, что вы имели в виду под подходом к кривой Гильберта. Тем не менее, это сложная проблема, и, учитывая N=100 и 10 миллионов точек данных, я не думаю, что какой-либо подход хорошо сохранит локальность и будет работать в разумные сроки. Я не думаю, что kd-деревья будут работать здесь.
Если поиск общего порядка не важен для вас, вы можете изучить хэширование на основе локальных условий и другие приблизительные схемы ближайших соседей. Иерархическое многомерное масштабирование с набором точек для уменьшения входного размера может дать вам хороший порядок, но опять же сомнительно в таком высоком измерении.
На самом деле это не «ответ», но я поигрался с приведенным выше кодом и просмотрел документ, и я не уверен, что он работает.
Помимо ошибки, обнаруженной Paul Chernoch, следующие
forпетля segfaults тоже.После исправления он все еще рисует вторую половину 2-битной 2-мерной кривой Гильберта в обратном направлении.
Никто не говорит о том, что в документе приводится не отображение между единичным интервалом и единичными кубами, а отображение целых чисел в большие кубы, и
Я не вижу здесь упоминания о странном представлении, которое в статье используется для беззнаковых целых чисел «Транспонировать».
Если вы нарисуете двумерные координаты в последнем столбце, вы увидите вторую половину кривой (позиция 8 и далее) в обратном направлении.
NB: HTranspose is the representation of H they use in the paper
H, HTranspose, mappedCoordinates
------------------------------------
0: (00, 00), (0, 0)
1: (00, 01), (1, 0)
2: (01, 00), (1, 1)
3: (01, 01), (0, 1)
4: (00, 10), (0, 2)
5: (00, 11), (0, 3)
6: (01, 10), (1, 3)
7: (01, 11), (1, 2)
8: (10, 00), (3, 2)
9: (10, 01), (3, 3)
10: (11, 00), (2, 3)
11: (11, 01), (2, 2)
12: (10, 10), (2, 1)
13: (10, 11), (3, 1)
14: (11, 10), (3, 0)
15: (11, 11), (2, 0)
Вот код (на C++). Я не сильно изменил код C в статье. Однако C ++ немного строже к типам, поэтому мне нужно было использовать
std::bitset. Единственная моя работа - это функция
makeHTranspose, который принимает обычное представление целого числа, а затем преобразует его в «транспонирование», которое все здесь используют.
#include <array>
#include <bitset>
#include <iostream>
#include <cmath>
namespace hilbert {
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> TransposeToAxes(std::array<std::bitset<num_bits>,num_dims> X)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
coord_t N = 2 << (num_bits-1);
// Gray decode by H ^ (H/2)
coord_t t = X[num_dims-1] >> 1;
for(size_t i = num_dims-1; i > 0; i-- ) // https://stackoverflow.com/a/10384110
X[i] ^= X[i-1];
X[0] ^= t;
// Undo excess work
for( coord_t Q = 2; Q != N; Q <<= 1 ) {
coord_t P = Q.to_ulong() - 1;
for( size_t i = num_dims-1; i > 0 ; i-- ){ // did the same stackoverflow thing
if( (X[i] & Q).any() )
X[0] ^= P;
else{
t = (X[0]^X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
}
return X;
}
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> AxesToTranspose(std::array<std::bitset<num_bits>, num_dims> X)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
coord_t M = 1 << (num_bits-1);
// Inverse undo
for( coord_t Q = M; Q.to_ulong() > 1; Q >>= 1 ) {
coord_t P = Q.to_ulong() - 1;
for(size_t i = 0; i < num_bits; i++ ){
if( (X[i] & Q).any() )
X[0] ^= P;
else{
coord_t t = (X[0]^X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
} // exchange
// Gray encode
for( size_t i = 1; i < num_bits; i++ )
X[i] ^= X[i-1];
coord_t t = 0;
for( coord_t Q = M; Q.to_ulong() > 1; Q >>= 1 ){
if( (X[num_dims-1] & Q).any() )
t ^= Q.to_ulong()-1;
}
for( size_t i = 0; i < num_bits; i++ )
X[i] ^= t;
return X;
}
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> makeHTranspose(unsigned int H)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
using big_coord_t = std::bitset<num_bits*num_dims>;
big_coord_t Hb = H;
coords_t X;
for(size_t dim = 0; dim < num_dims; ++dim){
coord_t tmp;
unsigned c = num_dims - 1;
for(size_t rbit = dim; rbit < num_bits*num_dims; rbit += num_dims){
tmp[c] =Hb[num_bits*num_dims - 1 - rbit];
c--;
}
X[dim] = tmp;
}
return X;
}
} //namespace hilbert
int main()
{
constexpr unsigned nb = 2;
constexpr unsigned nd = 2;
using coord_t = std::bitset<nb>;
using coords_t = std::array<coord_t,nd>;
std::cout << "NB: HTranspose is the representation of H they use in the paper\n";
std::cout << "H, HTranspose, mappedCoordinates \n";
std::cout << "------------------------------------\n";
for(unsigned H = 0; H < pow(2,nb*nd); ++H){
// H with the representation they use in the paper
coords_t weirdH = hilbert::makeHTranspose<nb,nd>(H);
std::cout << H << ": ("
<< weirdH[0] << ", "
<< weirdH[1] << "), ("
<< hilbert::TransposeToAxes<nb,nd>(weirdH)[0].to_ulong() << ", "
<< hilbert::TransposeToAxes<nb,nd>(weirdH)[1].to_ulong() << ")\n";
}
}