algorithm - peano - fractales
Asignación de valor N-dimensional a un punto en la curva de Hilbert (6)
Tengo un gran conjunto de puntos N-dimensionales (decenas de millones; N está cerca de 100).
Necesito asignar estos puntos a una sola dimensión mientras conserva la localidad espacial. Quiero usar la curva de relleno de espacio de Hilbert para hacerlo.
Para cada punto, quiero elegir el punto más cercano de la curva. El valor de Hilbert del punto (longitud de la curva desde el inicio de la curva hasta el punto elegido) es el valor de una sola dimensión que busco.
La computación no tiene que ser instantánea, pero espero que no sea más de varias horas en hardware decente para el hogar moderno.
¿Alguna sugerencia sobre la implementación? ¿Hay alguna biblioteca que me ayude? (El lenguaje no importa mucho)
Algoritmo para el mapeo de n-> 1 y 1-> n dado aquí "Cálculo de mapeos entre uno y valores n-dimensionales usando la curva de relleno de espacio de Hilbert" JK Lawder
Si busca "módulo SFC y superposición de Kademlia", encontrará un grupo que dice usarlo como parte de su sistema. Si ve la fuente, probablemente pueda extraer la función correspondiente.
Finalmente me derrumbé y desembolsé algo de dinero. El AIP (Instituto Americano de Física) tiene un artículo breve y bonito con el código fuente en C. "Programando la curva de Hilbert" de John Skilling (del AIP Conf. Proc. 707, 381 (2004)) tiene un apéndice con el código para mapeos en ambas direcciones. Funciona para cualquier número de dimensiones> 1, no es recursivo, no usa tablas de búsqueda de transición de estado que consumen grandes cantidades de memoria y, en su mayoría, utiliza operaciones de bits. Por lo tanto, es razonablemente rápido y tiene una buena huella de memoria.
Si elige comprar el artículo, descubrí un error en el código fuente.
La siguiente línea de código (que se encuentra en la función TransposetoAxes) tiene el error:
para (i = n-1; i> = 0; i--) X [i] ^ = X [i-1];
La corrección es cambiar el mayor que o igual (> =) a un mayor que (>). Sin esta corrección, se accede a la matriz X utilizando un índice negativo cuando la variable "i" se convierte en cero, lo que hace que el programa falle.
Recomiendo leer el artículo (que tiene siete páginas, incluido el código), ya que explica cómo funciona el algoritmo, que está lejos de ser obvio.
Traducí su código en C # para mi propio uso. El código sigue Skilling realiza la transformación en su lugar, sobrescribiendo el vector por el que pasas. Elegí hacer un clon del vector de entrada y devolver una nueva copia. Además, implementé los métodos como métodos de extensión.
El código de Skilling representa el índice de Hilbert como una transposición, almacenado como una matriz. Me parece más conveniente intercalar los bits y formar un único BigInteger (más útil en diccionarios, más fácil de iterar en bucles, etc.), pero optimicé esa operación y su inversa con números mágicos, operaciones de bits y similares, y el el código es extenso, así que lo he omitido.
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;
}
}
}
He publicado código de trabajo en C # a github.
No tengo claro cómo esto hará lo que quieras. Considere este caso trival en 3D:
001 ------ 101
|/ |/
| / | /
| 011 ------ 111
| | | |
| | | |
000 -|---- 100 |
/ | / |
/ | / |
010 ------ 110
que puede ser "Hilbertized" por la siguiente ruta:
001 -----> 101
/ /
/ /
011 111
^ |
| |
000 | 100 |
/ | / |
/ | / V
010 110
en el orden 1D:
000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100
Aquí está lo desagradable. Considere la lista de pares y distancias 1D a continuación:
000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1
En todos los casos, los valores de la mano izquierda y la derecha son la misma distancia 3D entre sí (+/- 1 en la primera posición), lo que parece implicar una "localidad espacial" similar. Pero linealizar mediante cualquier elección de ordenamiento dimensional (y, entonces z, luego z, en el ejemplo anterior) rompe esa localidad.
Otra forma de decir esto es que tomar un punto de partida y ordenar los puntos restantes por su distancia desde ese punto de partida proporcionará resultados significativamente diferentes. Tomando 000 como inicio, por ejemplo:
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
Este efecto crece exponencialmente con el número de dimensiones (suponiendo que cada dimensión tenga el mismo "tamaño").
Otra posibilidad sería construir un kd-tree en sus datos, y luego en un recorrido en orden del árbol para obtener el pedido. Construir el árbol kd solo requiere que tengas un buen algoritmo de búsqueda de mediana, del cual hay muchos.
Pasé un poco de tiempo traduciendo el código de Paul Chernoch a Java y limpiándolo. Es posible que haya un error en mi código, especialmente porque no tengo acceso al documento original. Sin embargo, pasa las pruebas unitarias que pude escribir. Está abajo.
Tenga en cuenta que he evaluado las curvas Z-Order y Hilbert para la indexación espacial en conjuntos de datos grandes. Debo decir que Z-Order ofrece mucha mejor calidad. Pero siéntete libre de probarlo por ti mismo.
/**
* 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;
}
No veo cómo puede usar una curva de Hilbert en una dimensión.
Si está interesado en asignar puntos a una dimensión inferior mientras conserva las distancias (con un error mínimo), entonces puede examinar los algoritmos de "Escalamiento multidimensional".
El recocido simulado es un enfoque.
Editar: Gracias por el comentario. Veo lo que querías decir con el enfoque de Hilbert Curve ahora. Sin embargo, este es un problema difícil, y dado N = 100 y 10 millones de puntos de datos, no creo que ningún enfoque preserve la localidad y se ejecute en un tiempo razonable. No creo que los kd funcionen aquí.
Si encontrar un pedido total no es importante para usted, entonces puede buscar el hashing basado en localidad y otros esquemas de vecinos más cercanos. El escalado jerárquico multidimensional con cubos de puntos para reducir el tamaño de entrada puede darle un buen orden, pero nuevamente es dudoso en una dimensión tan alta.