backend/transformer/hilbert.hpp

0001 /*
0002  * SPDX-PackageName: "covfie, a part of the ACTS project"
0003  * SPDX-FileCopyrightText: 2022 CERN
0004  * SPDX-License-Identifier: MPL-2.0
0005  */
0006
0007 #pragma once
0008
0009 #include <algorithm>
0010 #include <climits>
0011 #include <iostream>
0012 #include <memory>
0013 #include <numeric>
0014 #include <type_traits>
0015
0016 #include <covfie/core/array.hpp>
0017 #include <covfie/core/backend/primitive/array.hpp>
0018 #include <covfie/core/concepts.hpp>
0019 #include <covfie/core/parameter_pack.hpp>
0020 #include <covfie/core/qualifiers.hpp>
0021 #include <covfie/core/utility/binary_io.hpp>
0022 #include <covfie/core/utility/nd_map.hpp>
0023 #include <covfie/core/utility/nd_size.hpp>
0024 #include <covfie/core/utility/numeric.hpp>
0025 #include <covfie/core/vector.hpp>
0026
0027 namespace covfie::backend {
0028 template <
0029     concepts::vector_descriptor _input_vector_t,
0030     concepts::field_backend _storage_t>
0031 struct hilbert {
0032     using this_t = hilbert<_input_vector_t, _storage_t>;
0033     static constexpr bool is_initial = false;
0034
0035     using backend_t = _storage_t;
0036
0037     using contravariant_input_t =
0038         covfie::vector::array_vector_d<_input_vector_t>;
0039     using contravariant_output_t = typename backend_t::contravariant_input_t;
0040     using covariant_input_t = typename backend_t::covariant_output_t;
0041     using covariant_output_t = covariant_input_t;
0042
0043     using coordinate_t = typename contravariant_input_t::vector_t;
0044     using array_t = backend_t;
0045
0046     using configuration_t = utility::nd_size<contravariant_input_t::dimensions>;
0047
0048     static constexpr uint32_t IO_MAGIC_HEADER = 0xAB020004;
0049
0050     static_assert(
0051         contravariant_input_t::dimensions == 2,
0052         "Number of dimensions for input must be exactly two."
0053     );
0054
0055     COVFIE_HOST_DEVICE static void
0056     rot(std::size_t n,
0057         std::size_t * x,
0058         std::size_t * y,
0059         std::size_t rx,
0060         std::size_t ry)
0061     {
0062         if (ry == 0) {
0063             if (rx == 1) {
0064                 *x = n - 1 - *x;
0065                 *y = n - 1 - *y;
0066             }
0067
0068             std::size_t t = *x;
0069             *x = *y;
0070             *y = t;
0071         }
0072     }
0073
0074     COVFIE_HOST_DEVICE static std::size_t calculate_index(
0075         coordinate_t c,
0076         utility::nd_size<contravariant_input_t::dimensions> sizes
0077     )
0078     {
0079         // Borrowed from https://en.wikipedia.org/wiki/Hilbert_curve
0080
0081         std::size_t rx, ry, s, d = 0;
0082
0083         std::size_t x = c[0];
0084         std::size_t y = c[1];
0085
0086         // TODO: `sizes[0]` has to equal `sizes[1]`.
0087         for (s = sizes[0] / 2; s > 0; s /= 2) {
0088             rx = (x & s) > 0;
0089             ry = (y & s) > 0;
0090             d += s * s * ((3 * rx) ^ ry);
0091             rot(sizes[0], &x, &y, rx, ry);
0092         }
0093
0094         return d;
0095     }
0096
0097     template <typename T>
0098     static std::unique_ptr<
0099         std::decay_t<typename backend_t::covariant_output_t::vector_t>[]>
0100     make_hilbert_copy(const T & other)
0101     {
0102         configuration_t sizes = other.get_configuration();
0103
0104         std::unique_ptr<
0105             std::decay_t<typename backend_t::covariant_output_t::vector_t>[]>
0106             res = std::make_unique<std::decay_t<
0107                 typename backend_t::covariant_output_t::vector_t>[]>(
0108                 utility::ipow(
0109                     utility::round_pow2(
0110                         *std::max_element(sizes.begin(), sizes.end())
0111                     ),
0112                     contravariant_input_t::dimensions
0113                 )
0114             );
0115         typename T::parent_t::non_owning_data_t nother(other);
0116
0117         utility::nd_map<decltype(sizes)>(
0118             [&nother, &res](decltype(sizes) t) {
0119                 coordinate_t c;
0120
0121                 for (std::size_t i = 0; i < contravariant_input_t::dimensions;
0122                      ++i) {
0123                     c[i] = t[i];
0124                 }
0125
0126                 std::size_t idx = calculate_index(c);
0127
0128                 for (std::size_t i = 0; i < covariant_output_t::dimensions; ++i)
0129                 {
0130                     res[idx][i] = nother.at(c)[i];
0131                 }
0132             },
0133             sizes
0134         );
0135
0136         return res;
0137     }
0138
0139     struct owning_data_t {
0140         using parent_t = this_t;
0141
0142         owning_data_t() = default;
0143         owning_data_t(const owning_data_t &) = default;
0144         owning_data_t(owning_data_t &&) = default;
0145         owning_data_t & operator=(const owning_data_t &) = default;
0146         owning_data_t & operator=(owning_data_t &&) = default;
0147
0148         template <typename T>
0149         requires(std::same_as<
0150                  typename T::parent_t::configuration_t,
0151                  configuration_t> &&
0152                      std::constructible_from<
0153                          typename backend_t::owning_data_t,
0154                          std::size_t,
0155                          std::add_rvalue_reference_t<std::unique_ptr<std::decay_t<
0156                              typename backend_t::covariant_output_t::
0157                                  vector_t>[]>>>) explicit owning_data_t(const T &
0158                                                                             o)
0159             : m_sizes(o.get_configuration())
0160             , m_storage(
0161                   utility::ipow(
0162                       utility::round_pow2(
0163                           *std::max_element(m_sizes.begin(), m_sizes.end())
0164                       ),
0165                       contravariant_input_t::dimensions
0166                   ),
0167                   make_hilbert_copy(o)
0168               )
0169         {
0170         }
0171
0172         template <typename... Args>
0173         requires(std::constructible_from<
0174                  typename backend_t::owning_data_t,
0175                  std::
0176                      size_t>) explicit owning_data_t(configuration_t conf, Args... args)
0177             : m_sizes(conf)
0178             , m_storage(std::forward<Args>(args)...)
0179         {
0180         }
0181
0182         template <typename... Args>
0183         explicit owning_data_t(parameter_pack<configuration_t, Args...> && args)
0184             : m_sizes(args.x)
0185             , m_storage(std::move(args.xs))
0186         {
0187         }
0188
0189         explicit owning_data_t(parameter_pack<owning_data_t> && conf)
0190             : owning_data_t(std::move(conf.x))
0191         {
0192         }
0193
0194         explicit owning_data_t(
0195             const configuration_t & c, typename backend_t::owning_data_t && b
0196         )
0197             : m_sizes(c)
0198             , m_storage(std::forward<typename backend_t::owning_data_t>(b))
0199         {
0200         }
0201
0202         typename backend_t::owning_data_t & get_backend(void)
0203         {
0204             return m_storage;
0205         }
0206
0207         const typename backend_t::owning_data_t & get_backend(void) const
0208         {
0209             return m_storage;
0210         }
0211
0212         configuration_t get_configuration(void) const
0213         {
0214             return m_sizes;
0215         }
0216
0217         static owning_data_t read_binary(std::istream & fs)
0218         {
0219             utility::read_io_header(fs, IO_MAGIC_HEADER);
0220
0221             auto sizes = utility::read_binary<decltype(m_sizes)>(fs);
0222             auto be = decltype(m_storage)::read_binary(fs);
0223
0224             utility::read_io_footer(fs, IO_MAGIC_HEADER);
0225
0226             return owning_data_t(sizes, std::move(be));
0227         }
0228
0229         static void write_binary(std::ostream & fs, const owning_data_t & o)
0230         {
0231             utility::write_io_header(fs, IO_MAGIC_HEADER);
0232
0233             fs.write(
0234                 reinterpret_cast<const char *>(&o.m_sizes),
0235                 sizeof(decltype(m_sizes))
0236             );
0237
0238             decltype(m_storage)::write_binary(fs, o.m_storage);
0239
0240             utility::write_io_footer(fs, IO_MAGIC_HEADER);
0241         }
0242
0243         configuration_t m_sizes;
0244         typename backend_t::owning_data_t m_storage;
0245     };
0246
0247     struct non_owning_data_t {
0248         using parent_t = this_t;
0249
0250         non_owning_data_t(const owning_data_t & o)
0251             : m_sizes(o.m_sizes)
0252             , m_storage(o.m_storage)
0253         {
0254         }
0255
0256         COVFIE_HOST_DEVICE typename covariant_output_t::vector_t
0257         at(coordinate_t c) const
0258         {
0259 #ifndef NDEBUG
0260             for (std::size_t i = 0; i < contravariant_input_t::dimensions; ++i)
0261             {
0262                 assert(c[i] < m_sizes[i]);
0263             }
0264 #endif
0265
0266             return m_storage.at(calculate_index(c));
0267         }
0268
0269         typename backend_t::non_owning_data_t & get_backend(void)
0270         {
0271             return m_storage;
0272         }
0273
0274         const typename backend_t::non_owning_data_t & get_backend(void) const
0275         {
0276             return m_storage;
0277         }
0278
0279         configuration_t m_sizes;
0280         typename backend_t::non_owning_data_t m_storage;
0281     };
0282 };
0283 }