From 954cb4a1aba596d4e8d3c31ec7ee331fa0b1b676 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 5 Oct 2018 11:53:11 +0200 Subject: [PATCH] add inverse_Hilbert_space_filling_curve --- doc/news/changes/minor/20181005DenisDavydov | 4 + include/deal.II/base/utilities.h | 31 +++ source/base/utilities.cc | 206 +++++++++++++++++++- tests/base/utilities_13.cc | 145 ++++++++++++++ tests/base/utilities_13.output | 16 ++ 5 files changed, 400 insertions(+), 2 deletions(-) create mode 100644 doc/news/changes/minor/20181005DenisDavydov create mode 100644 tests/base/utilities_13.cc create mode 100644 tests/base/utilities_13.output diff --git a/doc/news/changes/minor/20181005DenisDavydov b/doc/news/changes/minor/20181005DenisDavydov new file mode 100644 index 0000000000..d8e1aa30c9 --- /dev/null +++ b/doc/news/changes/minor/20181005DenisDavydov @@ -0,0 +1,4 @@ +New: Add inverse_Hilbert_space_filling_curve() to map points in dim-dimensional +space to line index of the Hilbert space filling curve. +
+(Denis Davydov 2018/10/05) diff --git a/include/deal.II/base/utilities.h b/include/deal.II/base/utilities.h index 44d83cd612..1174563fa8 100644 --- a/include/deal.II/base/utilities.h +++ b/include/deal.II/base/utilities.h @@ -52,6 +52,9 @@ DEAL_II_NAMESPACE_OPEN +// forward declare Point +template +class Point; /** * A namespace for utility functions that are not particularly specific to @@ -72,6 +75,34 @@ namespace Utilities std::string dealii_version_string(); + /** + * Assign to each point in @p points an index using the Hilbert space filling curve. + * To that end, a bounding box for @p points will be determined, based on which their + * integer coordinates are calculated. + * The linear index is given as a dim-collection of bits, from high to low. + * This is done in order to keep the maximum resolution in terms of bit depth + * along each axis. Note that this dim-integer index can still be easily used + * for sorting and ordering, for example using the lexicographic ordering of + * tuples of integers. + * + * The depth of the Hilbert curve (i.e. the number of bits per dimension) by + * default is equal to 64. + */ + template + std::vector> + inverse_Hilbert_space_filling_curve( + const std::vector> &points, + const int bits_per_dim = 64); + + /** + * Same as above, but for points in integer coordinates. + */ + template + std::vector> + inverse_Hilbert_space_filling_curve( + const std::vector> &points, + const int bits_per_dim = 64); + /** * Convert a number @p value to a string, with as many digits as given to * fill with leading zeros. diff --git a/source/base/utilities.cc b/source/base/utilities.cc index d3bae9270b..baee0da002 100644 --- a/source/base/utilities.cc +++ b/source/base/utilities.cc @@ -24,6 +24,7 @@ #include #include +#include #include #include @@ -64,8 +65,6 @@ # include #endif - - DEAL_II_NAMESPACE_OPEN @@ -91,6 +90,184 @@ namespace Utilities + template + std::vector> + inverse_Hilbert_space_filling_curve( + const std::vector> &points, + const int bits_per_dim) + { + using Integer = std::uint64_t; + // take floating point number hopefully with mantissa >= 64bit + using LongDouble = long double; + + // return if there is nothing to do + if (points.size() == 0) + return std::vector>(); + + // get bounding box: + Point bl = points[0], tr = points[0]; + for (const auto &p : points) + for (unsigned int d = 0; d < dim; ++d) + { + const double cid = p[d]; + bl[d] = std::min(cid, bl[d]); + tr[d] = std::max(cid, tr[d]); + } + + std::array extents; + for (unsigned int i = 0; i < dim; ++i) + { + extents[i] = + static_cast(tr[i]) - static_cast(bl[i]); + Assert(extents[i] > 0., ExcMessage("Bounding box is degenerate.")); + } + + // make sure our conversion from fractional coordinates to + // Integers work as expected, namely our cast (LongDouble)max_int + const int min_bits = + std::min(bits_per_dim, + std::min(std::numeric_limits::digits, + std::numeric_limits::digits)); + + // based on that get the maximum integer: + const Integer max_int = (min_bits == std::numeric_limits::digits ? + std::numeric_limits::max() : + ((Integer)1 << min_bits) - 1); + + std::vector> int_points(points.size()); + + for (unsigned int i = 0; i < points.size(); ++i) + { + // convert into integers: + for (unsigned int d = 0; d < dim; ++d) + { + const LongDouble v = (static_cast(points[i][d]) - + static_cast(bl[d])) / + extents[d]; + Assert(v >= 0. && v <= 1., ExcInternalError()); + int_points[i][d] = + static_cast(v * static_cast(max_int)); + } + } + + // note that we call this with "min_bits" + return inverse_Hilbert_space_filling_curve(int_points, min_bits); + } + + + + template + std::vector> + inverse_Hilbert_space_filling_curve( + const std::vector> &points, + const int bits_per_dim) + { + using Integer = std::uint64_t; + + std::vector> int_points(points); + + std::vector> res(int_points.size()); + + // follow + // J. Skilling, Programming the Hilbert curve, AIP Conf. Proc. 707, 381 + // (2004); http://dx.doi.org/10.1063/1.1751381 also see + // https://stackoverflow.com/questions/499166/mapping-n-dimensional-value-to-a-point-on-hilbert-curve + // https://gitlab.com/octopus-code/octopus/blob/develop/src/grid/hilbert.c + // https://github.com/trilinos/Trilinos/blob/master/packages/zoltan/src/hsfc/hsfc_hilbert.c + // (Zoltan_HSFC_InvHilbertXd) + // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp + + // now we can map to 1D coordinate stored in Transpose format + // adopt AxestoTranspose function from the paper, that + // transforms in-place between geometrical axes and Hilbert transpose. + // 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] + + // Depth of the Hilbert curve + Assert(bits_per_dim <= std::numeric_limits::digits, + ExcMessage("This integer type can not hold " + + std::to_string(bits_per_dim) + " bits.")); + + const Integer M = 1 << (bits_per_dim - 1); // largest bit + + for (unsigned int index = 0; index < int_points.size(); ++index) + { + auto &X = int_points[index]; + auto &L = res[index]; + + // Inverse undo + for (Integer q = M; q > 1; q >>= 1) + { + const Integer p = q - 1; + for (unsigned int i = 0; i < dim; i++) + { + // invert + if (X[i] & q) + { + X[0] ^= p; + } + // exchange + else + { + const Integer t = (X[0] ^ X[i]) & p; + X[0] ^= t; + X[i] ^= t; + } + } + } + + // Gray encode (inverse of decode) + for (unsigned int i = 1; i < dim; i++) + X[i] ^= X[i - 1]; + + Integer t = 0; + for (Integer q = M; q > 1; q >>= 1) + if (X[dim - 1] & q) + t ^= q - 1; + for (unsigned int i = 0; i < dim; i++) + X[i] ^= t; + + // now we need to go from index stored in transpose format to + // consecutive format, which is better suited for comparators. + // we could interleave into some big unsigned int... + // https://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ + // https://stackoverflow.com/questions/4431522/given-2-16-bit-ints-can-i-interleave-those-bits-to-form-a-single-32-bit-int + // ...but we would loose spatial resolution! + + // interleave using brute force, follow TransposetoLine from + // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp + { + Integer p = M; + unsigned int j = 0; + for (unsigned int i = 0; i < dim; ++i) + { + L[i] = 0; + // go through bits using a mask q + for (Integer q = M; q > 0; q >>= 1) + { + if (X[j] & p) + L[i] |= q; + if (++j == dim) + { + j = 0; + p >>= 1; + } + } + } + } + + } // end of the loop over points + + return res; + } + + + std::string int_to_string(const unsigned int value, const unsigned int digits) { @@ -900,6 +1077,31 @@ namespace Utilities template std::string to_string(long double, unsigned int); + template std::vector> + inverse_Hilbert_space_filling_curve<1, double>( + const std::vector> &, + const int); + template std::vector> + inverse_Hilbert_space_filling_curve<1>( + const std::vector> &, + const int); + template std::vector> + inverse_Hilbert_space_filling_curve<2, double>( + const std::vector> &, + const int); + template std::vector> + inverse_Hilbert_space_filling_curve<2>( + const std::vector> &, + const int); + template std::vector> + inverse_Hilbert_space_filling_curve<3, double>( + const std::vector> &, + const int); + template std::vector> + inverse_Hilbert_space_filling_curve<3>( + const std::vector> &, + const int); + } // namespace Utilities DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/utilities_13.cc b/tests/base/utilities_13.cc new file mode 100644 index 0000000000..10fa8b65b9 --- /dev/null +++ b/tests/base/utilities_13.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test Utilities::inverse_Hilbert_space_filling_curve in 2D +// on the following nine points with bits=4: + +// | +// 15 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | | | | | | | | | | | | | +// | @ @---@ @ @ @---@ @ @ @---@ @ @ @---@ @ +// | | | | | | | | | +// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | | | | | +// | @---@ @---@---@---@ @---@ @---@ @---@---@---@ @---@ +// | | | | | +// | @ @---@---@ @---@---@ @ @ @---@---@ @---@---@ @ +// | | | | | | | | | | | | | +// Axes[1]| @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | +// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | | | | | | | | | +// | @ @---@---@ @---@---@ @---@ @---@---@ @---@---@ @ +// | | | +// | @---@ @---@---@ @---@---@ @---@---@ @---@---@ @---@ +// | | | | | | | | | | | +// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | | | +// | @ @---@ @ @---@ @---@ @---@ @---@ @ @---@ @ +// | | | | | | | | | | | | | | | +// | @---@ @---@ @ @---@---@ @---@---@ @ @---@ @---@ +// | | | +// 3 | 5---6 9---@ @ @---@---@ @---@---@ @ @---@ @---@ +// | | | | | | | | | | | | | | | +// 2 | 4 7---8 @ @---@ @---@ @---@ @---@ @ @---@ @ +// | | | | | | | +// 1 | 3---2 @---@ @---@ @---@ @---@ @---@ @---@ @---@ +// | | | | | | | | | | | +// 0 | 0---1 @---@---@ @---@---@ @---@---@ @---@---@ @--255 +// | +// ------------------------------------------------------------------- +// 0 1 2 3 ---> Axes[0] 15 + + +#include + +#include "../tests.h" + + +void +test1() +{ + deallog << "test1:" << std::endl; + const std::vector> points = {Point<2>(0, 0), + Point<2>(1, 0), + Point<2>(1, 1), + Point<2>(0, 1), + Point<2>(0, 2), + Point<2>(0, 3), + Point<2>(1, 3), + Point<2>(1, 2), + Point<2>(2, 2), + Point<2>(2, 3), + Point<2>(15, 15), + Point<2>(15, 0)}; + + const auto res = Utilities::inverse_Hilbert_space_filling_curve(points, 4); + + for (const auto &p : res) + deallog << p[0] << " " << p[1] << std::endl; +} + +// https://github.com/aditi137/Hilbert/blob/master/Hilbert/test.py +// test int version in 3D +void +test2() +{ + deallog << "test2:" << std::endl; + std::vector> points = {{5, 10, 20}}; + + const auto res = Utilities::inverse_Hilbert_space_filling_curve<3>(points, 5); + + for (const auto &p : res) + deallog << p[0] << " " << p[1] << " " << p[2] << std::endl; +} + + +// finally, a quick test of fractional coordinates -> int conversion inside +// inverse_Hilbert_space_filling_curve(): +template +void +test3() +{ + const int bits_int = std::numeric_limits::digits; + const int bits_num = std::numeric_limits::digits; + + const int min_bits = std::min(bits_int, bits_num); + + const Integer max = + (min_bits == bits_int ? std::numeric_limits::max() : + ((Integer)1 << min_bits) - 1); + + // make sure conversion is within bounds: + const std::vector points = {0.2, 0.7}; + for (const auto &v : points) + { + const Integer i = (Integer)(v * (Number)max); + Assert(i > 0 && i < max, ExcInternalError()); + } + + { + const Number v = 0.; + const Integer i = (Integer)(v * (Number)max); + Assert(i == 0, ExcInternalError()); + } + + { + const Number v = 1.; + const Integer i = (Integer)(v * (Number)max); + Assert(i == max, ExcInternalError()); + } +} + +int +main() +{ + initlog(); + + test1(); + test2(); + + test3(); +} diff --git a/tests/base/utilities_13.output b/tests/base/utilities_13.output new file mode 100644 index 0000000000..a5c97ad1b5 --- /dev/null +++ b/tests/base/utilities_13.output @@ -0,0 +1,16 @@ + +DEAL::test1: +DEAL::0 0 +DEAL::0 1 +DEAL::0 2 +DEAL::0 3 +DEAL::0 4 +DEAL::0 5 +DEAL::0 6 +DEAL::0 7 +DEAL::0 8 +DEAL::0 9 +DEAL::10 10 +DEAL::15 15 +DEAL::test2: +DEAL::7 21 25 -- 2.39.5