algebra/utils/approximately_equal.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s)
0012 #include "detray/algebra/common/math.hpp"
0013 #include "detray/algebra/concepts.hpp"
0014 #include "detray/definitions/detail/qualifiers.hpp"
0015
0016 // System include(s)
0017 #include <concepts>
0018 #include <limits>
0019
0020 namespace detray::algebra {
0021
0022 /// Compare two scalars according to a max relative error tolerance
0023 /// @see
0024 /// https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
0025 ///
0026 /// @note This is by no means safe for all comparisons. Use with caution!
0027 ///
0028 /// @param a first value
0029 /// @param b second value
0030 /// @param rel_error maximal relative error
0031 ///
0032 /// @returns true if the two values are approximaterly equal
0033 template <concepts::scalar scalar_t, concepts::value value_t>
0034   requires std::convertible_to<scalar_t, value_t>
0035 DETRAY_HOST_DEVICE constexpr auto approx_equal(
0036     const scalar_t a, const scalar_t b,
0037     const value_t rel_error = 16.f * std::numeric_limits<value_t>::epsilon(),
0038     const value_t max_error = std::numeric_limits<value_t>::epsilon()) {
0039   if constexpr (std::integral<scalar_t>) {
0040     return a == b;
0041   } else {
0042     // Calculate the difference.
0043     const scalar_t diff{math::fabs(a - b)};
0044
0045     // If the numbers are close to zero
0046     if (diff <= max_error) {
0047       return true;
0048     }
0049
0050     // Find the largest value to scale the epsilon
0051     const scalar_t largest = math::max(math::fabs(a), math::fabs(b));
0052
0053     return (diff <= largest * rel_error);
0054   }
0055 }
0056
0057 /// Elementwise compare two vectors according to a max relative error tolerance
0058 ///
0059 /// @note This is by no means safe for all comparisons. Use with caution!
0060 ///
0061 /// @param v1 first vector
0062 /// @param v2 second vector
0063 /// @param rel_error maximal relative error
0064 ///
0065 /// @returns true if the two vectors are approximaterly equal
0066 template <typename vector1_t, typename vector2_t>
0067   requires((concepts::vector<vector1_t> || concepts::point<vector1_t>) &&
0068            (concepts::vector<vector2_t> || concepts::point<vector2_t>) &&
0069            !concepts::simd_scalar<vector1_t> &&
0070            !concepts::simd_scalar<vector2_t>)
0071 DETRAY_HOST_DEVICE constexpr auto approx_equal(
0072     const vector1_t& v1, const vector2_t& v2,
0073     const detray::traits::value_t<vector1_t> rel_error =
0074         16.f *
0075         std::numeric_limits<detray::traits::value_t<vector1_t>>::epsilon(),
0076     const detray::traits::value_t<vector1_t> max_error =
0077         std::numeric_limits<detray::traits::value_t<vector1_t>>::epsilon()) {
0078   static_assert(std::same_as<detray::traits::value_t<vector1_t>,
0079                              detray::traits::value_t<vector2_t>>);
0080
0081   using index_t = detray::traits::index_t<vector1_t>;
0082
0083   constexpr index_t size{detray::traits::size<vector1_t>};
0084
0085   bool ret{true};
0086   for (index_t i = 0; i < size; ++i) {
0087     ret = ret &&
0088           ::detray::algebra::approx_equal(v1[i], v2[i], rel_error, max_error);
0089   }
0090
0091   return ret;
0092 }
0093
0094 /// Elementwise compare two column matrices according to a max relative error
0095 /// tolerance
0096 ///
0097 /// @note This is by no means safe for all comparisons. Use with caution!
0098 ///
0099 /// @param v1 first column matrix
0100 /// @param v2 second column matrix
0101 /// @param rel_error maximal relative error
0102 ///
0103 /// @returns true if the two column matrices are approximaterly equal
0104 template <concepts::column_matrix vector1_t, concepts::column_matrix vector2_t>
0105 DETRAY_HOST_DEVICE constexpr auto approx_equal(
0106     const vector1_t& v1, const vector2_t& v2,
0107     const detray::traits::value_t<vector1_t> rel_error =
0108         16.f *
0109         std::numeric_limits<detray::traits::value_t<vector1_t>>::epsilon(),
0110     const detray::traits::value_t<vector1_t> max_error =
0111         std::numeric_limits<detray::traits::value_t<vector1_t>>::epsilon()) {
0112   static_assert(std::same_as<detray::traits::value_t<vector1_t>,
0113                              detray::traits::value_t<vector2_t>>);
0114
0115   using index_t = detray::traits::index_t<vector1_t>;
0116   using element_getter_t = detray::traits::element_getter_t<vector1_t>;
0117
0118   constexpr index_t rows{detray::traits::rows<vector1_t>};
0119
0120   bool ret{true};
0121   for (index_t i = 0; i < rows; ++i) {
0122     ret = ret && ::detray::algebra::approx_equal(element_getter_t{}(v1, i, 0),
0123                                                  element_getter_t{}(v2, i, 0),
0124                                                  rel_error, max_error);
0125   }
0126
0127   return ret;
0128 }
0129
0130 /// Elementwise compare two matrices according to a max relative error tolerance
0131 ///
0132 /// @note This is by no means safe for all comparisons. Use with caution!
0133 ///
0134 /// @param m1 first matrix
0135 /// @param m2 second matrix
0136 /// @param rel_error maximal relative error
0137 ///
0138 /// @returns true if the two matrices are approximaterly equal
0139 template <concepts::matrix matrix1_t, concepts::matrix matrix2_t>
0140 DETRAY_HOST_DEVICE constexpr auto approx_equal(
0141     const matrix1_t& m1, const matrix2_t& m2,
0142     const detray::traits::value_t<matrix1_t> rel_error =
0143         16.f *
0144         std::numeric_limits<detray::traits::value_t<matrix1_t>>::epsilon(),
0145     const detray::traits::value_t<matrix1_t> max_error =
0146         std::numeric_limits<detray::traits::value_t<matrix1_t>>::epsilon()) {
0147   static_assert(std::same_as<detray::traits::value_t<matrix1_t>,
0148                              detray::traits::value_t<matrix2_t>>);
0149
0150   using index_t = detray::traits::index_t<matrix1_t>;
0151   using element_getter_t = detray::traits::element_getter_t<matrix1_t>;
0152
0153   constexpr index_t rows{detray::traits::rows<matrix1_t>};
0154   constexpr index_t columns{detray::traits::columns<matrix1_t>};
0155
0156   bool ret{true};
0157   for (index_t j = 0; j < columns; ++j) {
0158     for (index_t i = 0; i < rows; ++i) {
0159       ret = ret && ::detray::algebra::approx_equal(element_getter_t{}(m1, i, j),
0160                                                    element_getter_t{}(m2, i, j),
0161                                                    rel_error, max_error);
0162     }
0163   }
0164
0165   return ret;
0166 }
0167
0168 /// Elementwise compare two transforms according to a max relative error
0169 /// tolerance
0170 ///
0171 /// @note This is by no means safe for all comparisons. Use with caution!
0172 ///
0173 /// @param trf1 first transform
0174 /// @param trf2 second transform
0175 /// @param rel_error maximal relative error
0176 ///
0177 /// @returns true if the two transforms are approximaterly equal
0178 template <concepts::transform3D transform_t>
0179 DETRAY_HOST_DEVICE constexpr auto approx_equal(
0180     const transform_t& trf1, const transform_t& trf2,
0181     const detray::traits::value_t<typename transform_t::scalar_type> rel_error =
0182         16.f * std::numeric_limits<detray::traits::value_t<
0183                    typename transform_t::scalar_type>>::epsilon(),
0184     const detray::traits::value_t<typename transform_t::scalar_type> max_error =
0185         std::numeric_limits<detray::traits::value_t<
0186             typename transform_t::scalar_type>>::epsilon()) {
0187   return (::detray::algebra::approx_equal(trf1.matrix(), trf2.matrix(),
0188                                           rel_error, max_error) &&
0189           ::detray::algebra::approx_equal(trf1.matrix_inverse(),
0190                                           trf2.matrix_inverse(), rel_error,
0191                                           max_error));
0192 }
0193
0194 }  // namespace detray::algebra