backend/transformer/linear.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 <cmath>
0010 #include <cstddef>
0011 #include <iostream>
0012 #include <type_traits>
0013 #include <variant>
0014
0015 #include <covfie/core/concepts.hpp>
0016 #include <covfie/core/definitions.hpp>
0017 #include <covfie/core/parameter_pack.hpp>
0018 #include <covfie/core/qualifiers.hpp>
0019 #include <covfie/core/utility/binary_io.hpp>
0020 #include <covfie/core/vector.hpp>
0021
0022 namespace covfie::backend {
0023 template <
0024     concepts::field_backend _backend_t,
0025     concepts::vector_descriptor _input_vector_d = covfie::vector::
0026         vector_d<float, _backend_t::contravariant_input_t::dimensions>>
0027 struct linear {
0028     using this_t = linear<_backend_t, _input_vector_d>;
0029     static constexpr bool is_initial = false;
0030
0031     using input_scalar_type = typename _input_vector_d::type;
0032     using backend_t = _backend_t;
0033
0034     using contravariant_input_t =
0035         covfie::vector::array_vector_d<_input_vector_d>;
0036     using contravariant_output_t = typename backend_t::contravariant_input_t;
0037     using covariant_input_t = typename backend_t::covariant_output_t;
0038     using covariant_output_t =
0039         covfie::vector::array_vector_d<typename covariant_input_t::vector_d>;
0040
0041     static_assert(
0042         std::is_floating_point_v<typename _input_vector_d::type>,
0043         "Linear interpolation contravariant input must have a "
0044         "floating point scalar type."
0045     );
0046     static_assert(
0047         std::is_floating_point_v<typename covariant_input_t::scalar_t>,
0048         "Linear interpolation covariant input must have a "
0049         "floating point scalar type."
0050     );
0051     static_assert(
0052         _input_vector_d::size == backend_t::contravariant_input_t::dimensions,
0053         "Linear interpolation contravariant input must have the "
0054         "same size as the backend contravariant input."
0055     );
0056     static_assert(
0057         std::is_object_v<typename covariant_output_t::vector_t>,
0058         "Covariant input type of linear interpolator must be an object type."
0059     );
0060
0061     using configuration_t = std::monostate;
0062
0063     static constexpr uint32_t IO_MAGIC_HEADER = 0xAB020005;
0064
0065     struct owning_data_t {
0066         using parent_t = this_t;
0067
0068         owning_data_t() = default;
0069
0070         template <typename... Args>
0071         explicit owning_data_t(configuration_t, Args... args)
0072             : m_backend(std::forward<Args>(args)...)
0073         {
0074         }
0075
0076         template <typename T>
0077         requires(std::same_as<
0078                  typename T::parent_t::configuration_t,
0079                  configuration_t>) explicit owning_data_t(const T & o)
0080             : m_backend(o.m_backend)
0081         {
0082         }
0083
0084         explicit owning_data_t(const typename backend_t::owning_data_t & o)
0085             : m_backend(o)
0086         {
0087         }
0088
0089         explicit owning_data_t(
0090             const configuration_t &, typename backend_t::owning_data_t && b
0091         )
0092             : m_backend(std::forward<typename backend_t::owning_data_t>(b))
0093         {
0094         }
0095
0096         template <typename... Args>
0097         explicit owning_data_t(parameter_pack<configuration_t, Args...> && args)
0098             : m_backend(std::move(args.xs))
0099         {
0100         }
0101
0102         explicit owning_data_t(parameter_pack<owning_data_t> && conf)
0103             : owning_data_t(std::move(conf.x))
0104         {
0105         }
0106
0107         typename backend_t::owning_data_t & get_backend(void)
0108         {
0109             return m_backend;
0110         }
0111
0112         const typename backend_t::owning_data_t & get_backend(void) const
0113         {
0114             return m_backend;
0115         }
0116
0117         configuration_t get_configuration(void) const
0118         {
0119             return {};
0120         }
0121
0122         static owning_data_t read_binary(std::istream & fs)
0123         {
0124             auto be = decltype(m_backend)::read_binary(fs);
0125
0126             return owning_data_t(configuration_t{}, std::move(be));
0127         }
0128
0129         static void write_binary(std::ostream & fs, const owning_data_t & o)
0130         {
0131             decltype(m_backend)::write_binary(fs, o.m_backend);
0132         }
0133
0134         typename backend_t::owning_data_t m_backend;
0135     };
0136
0137     struct non_owning_data_t {
0138         using parent_t = this_t;
0139
0140         non_owning_data_t(const owning_data_t & src)
0141             : m_backend(src.m_backend)
0142         {
0143         }
0144
0145         template <std::size_t... Is>
0146         COVFIE_HOST_DEVICE typename contravariant_output_t::vector_t
0147         _backend_index_helper(typename contravariant_output_t::vector_t coord, std::size_t n, std::index_sequence<Is...>)
0148             const
0149         {
0150             return {static_cast<typename decltype(m_backend
0151             )::parent_t::contravariant_input_t::scalar_t>(
0152                 coord[Is] + ((n & (std::size_t(1) << Is)) ? 1 : 0)
0153             )...};
0154         }
0155
0156         COVFIE_HOST_DEVICE typename covariant_output_t::vector_t
0157         at(typename contravariant_input_t::vector_t coord) const
0158         {
0159             if constexpr (covariant_output_t::dimensions == 1) {
0160                 typename contravariant_output_t::scalar_t i =
0161                     static_cast<typename contravariant_output_t::scalar_t>(
0162                         coord[0]
0163                     );
0164
0165                 input_scalar_type a = coord[0] - std::trunc(coord[0]);
0166
0167                 input_scalar_type ra = static_cast<input_scalar_type>(1.) - a;
0168
0169                 std::remove_reference_t<typename covariant_output_t::vector_t>
0170                     pc[2];
0171
0172                 for (std::size_t n = 0; n < 2; ++n) {
0173                     pc[n] =
0174                         m_backend.at({static_cast<typename decltype(m_backend
0175                         )::parent_t::contravariant_input_t::scalar_t>(
0176                             i + ((n & 1) ? 1 : 0)
0177                         )});
0178                 }
0179
0180                 typename covariant_output_t::vector_t rv;
0181
0182                 for (std::size_t q = 0; q < covariant_output_t::dimensions; ++q)
0183                 {
0184                     rv[q] = ra * static_cast<input_scalar_type>(pc[0][q]) +
0185                             a * static_cast<input_scalar_type>(pc[1][q]);
0186                 }
0187
0188                 return rv;
0189             } else if constexpr (covariant_output_t::dimensions == 2) {
0190                 typename contravariant_output_t::scalar_t i =
0191                     static_cast<typename contravariant_output_t::scalar_t>(
0192                         coord[0]
0193                     );
0194                 typename contravariant_output_t::scalar_t j =
0195                     static_cast<typename contravariant_output_t::scalar_t>(
0196                         coord[1]
0197                     );
0198
0199                 input_scalar_type a = coord[0] - std::trunc(coord[0]);
0200                 input_scalar_type b = coord[1] - std::trunc(coord[1]);
0201
0202                 input_scalar_type ra = static_cast<input_scalar_type>(1.) - a;
0203                 input_scalar_type rb = static_cast<input_scalar_type>(1.) - b;
0204
0205                 std::remove_reference_t<typename covariant_output_t::vector_t>
0206                     pc[4];
0207
0208                 for (std::size_t n = 0; n < 4; ++n) {
0209                     pc[n] = m_backend.at(
0210                         {static_cast<typename decltype(m_backend
0211                          )::parent_t::contravariant_input_t::scalar_t>(
0212                              i + ((n & 2) ? 1 : 0)
0213                          ),
0214                          static_cast<typename decltype(m_backend
0215                          )::parent_t::contravariant_input_t::scalar_t>(
0216                              j + ((n & 1) ? 1 : 0)
0217                          )}
0218                     );
0219                 }
0220
0221                 typename covariant_output_t::vector_t rv;
0222
0223                 for (std::size_t q = 0; q < covariant_output_t::dimensions; ++q)
0224                 {
0225                     rv[q] = ra * rb * static_cast<input_scalar_type>(pc[0][q]) +
0226                             ra * b * static_cast<input_scalar_type>(pc[1][q]) +
0227                             a * rb * static_cast<input_scalar_type>(pc[2][q]) +
0228                             a * b * static_cast<input_scalar_type>(pc[3][q]);
0229                 }
0230
0231                 return rv;
0232             } else if constexpr (covariant_output_t::dimensions == 3) {
0233                 typename contravariant_output_t::scalar_t i =
0234                     static_cast<typename contravariant_output_t::scalar_t>(
0235                         coord[0]
0236                     );
0237                 typename contravariant_output_t::scalar_t j =
0238                     static_cast<typename contravariant_output_t::scalar_t>(
0239                         coord[1]
0240                     );
0241                 typename contravariant_output_t::scalar_t k =
0242                     static_cast<typename contravariant_output_t::scalar_t>(
0243                         coord[2]
0244                     );
0245
0246                 input_scalar_type a = coord[0] - std::trunc(coord[0]);
0247                 input_scalar_type b = coord[1] - std::trunc(coord[1]);
0248                 input_scalar_type c = coord[2] - std::trunc(coord[2]);
0249
0250                 input_scalar_type ra = static_cast<input_scalar_type>(1.) - a;
0251                 input_scalar_type rb = static_cast<input_scalar_type>(1.) - b;
0252                 input_scalar_type rc = static_cast<input_scalar_type>(1.) - c;
0253
0254                 std::remove_reference_t<typename covariant_output_t::vector_t>
0255                     pc[8];
0256
0257                 for (std::size_t n = 0; n < 8; ++n) {
0258                     pc[n] = m_backend.at(
0259                         {static_cast<typename decltype(m_backend
0260                          )::parent_t::contravariant_input_t::scalar_t>(
0261                              i + ((n & 4) ? 1 : 0)
0262                          ),
0263                          static_cast<typename decltype(m_backend
0264                          )::parent_t::contravariant_input_t::scalar_t>(
0265                              j + ((n & 2) ? 1 : 0)
0266                          ),
0267                          static_cast<typename decltype(m_backend
0268                          )::parent_t::contravariant_input_t::scalar_t>(
0269                              k + ((n & 1) ? 1 : 0)
0270                          )}
0271                     );
0272                 }
0273
0274                 typename covariant_output_t::vector_t rv;
0275
0276                 for (std::size_t q = 0; q < covariant_output_t::dimensions; ++q)
0277                 {
0278                     rv[q] =
0279                         ra * rb * rc *
0280                             static_cast<input_scalar_type>(pc[0][q]) +
0281                         ra * rb * c * static_cast<input_scalar_type>(pc[1][q]) +
0282                         ra * b * rc * static_cast<input_scalar_type>(pc[2][q]) +
0283                         ra * b * c * static_cast<input_scalar_type>(pc[3][q]) +
0284                         a * rb * rc * static_cast<input_scalar_type>(pc[4][q]) +
0285                         a * rb * c * static_cast<input_scalar_type>(pc[5][q]) +
0286                         a * b * rc * static_cast<input_scalar_type>(pc[6][q]) +
0287                         a * b * c * static_cast<input_scalar_type>(pc[7][q]);
0288                 }
0289
0290                 return rv;
0291             } else {
0292                 typename contravariant_output_t::vector_t is;
0293
0294                 for (std::size_t n = 0; n < contravariant_output_t::dimensions;
0295                      ++n)
0296                 {
0297                     is[n] =
0298                         static_cast<contravariant_output_t::scalar_t>(coord[n]);
0299                 }
0300
0301                 input_scalar_type vs[contravariant_output_t::dimensions];
0302
0303                 for (std::size_t n = 0; n < contravariant_output_t::dimensions;
0304                      ++n)
0305                 {
0306                     vs[n] = coord[n] - std::trunc(coord[n]);
0307                 }
0308
0309                 input_scalar_type rs[contravariant_output_t::dimensions];
0310
0311                 for (std::size_t n = 0; n < contravariant_output_t::dimensions;
0312                      ++n)
0313                 {
0314                     rs[n] = static_cast<input_scalar_type>(1.) - vs[n];
0315                 }
0316
0317                 std::remove_reference_t<typename covariant_output_t::vector_t>
0318                     pc[std::size_t(1) << covariant_output_t::dimensions];
0319
0320                 for (std::size_t n = 0;
0321                      n < std::size_t(1) << covariant_output_t::dimensions;
0322                      ++n)
0323                 {
0324                     pc[n] = m_backend.at(_backend_index_helper(
0325                         is,
0326                         n,
0327                         std::make_index_sequence<
0328                             contravariant_input_t::dimensions>{}
0329                     ));
0330                 }
0331
0332                 typename covariant_output_t::vector_t rv;
0333
0334                 for (std::size_t q = 0; q < covariant_output_t::dimensions; ++q)
0335                 {
0336                     rv[q] = 0.f;
0337
0338                     for (std::size_t n = 0;
0339                          n < std::size_t(1) << covariant_output_t::dimensions;
0340                          ++n)
0341                     {
0342                         input_scalar_type f{1.};
0343
0344                         for (std::size_t m = 0;
0345                              m < covariant_output_t::dimensions;
0346                              ++m)
0347                         {
0348                             if (n & (std::size_t(1) << m)) {
0349                                 f *= vs[m];
0350                             } else {
0351                                 f *= rs[m];
0352                             }
0353                         }
0354
0355                         rv[q] += f * static_cast<input_scalar_type>(pc[n][q]);
0356                     }
0357                 }
0358
0359                 return rv;
0360             }
0361         }
0362
0363         typename backend_t::non_owning_data_t & get_backend(void)
0364         {
0365             return m_backend;
0366         }
0367
0368         const typename backend_t::non_owning_data_t & get_backend(void) const
0369         {
0370             return m_backend;
0371         }
0372
0373         typename backend_t::non_owning_data_t m_backend;
0374     };
0375 };
0376 }