Line data    Source code

       1             : // SPDX-License-Identifier: BSD-3-Clause
       2             : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
       3             : /// @file
       4             : /// @author Simon Heybrock
       5             : #pragma once
       6             : 
       7             : #include <cmath>
       8             : #include <type_traits>
       9             : 
      10             : #include "scipp/common/numeric.h"
      11             : #include "scipp/common/overloaded.h"
      12             : #include "scipp/core/eigen.h"
      13             : #include "scipp/core/element/arg_list.h"
      14             : #include "scipp/core/transform_common.h"
      15             : 
      16             : namespace scipp::core::element {
      17             : 
      18             : template <class... Extra>
      19             : constexpr auto special_value_args =
      20             :     arg_list<int32_t, int64_t, double, float, Extra...>;
      21             : 
      22             : // WARNING: When adding support for new spatial types (or other types containing
      23             : // floating-point elements) you must specialize numeric::isnan and friends.
      24             : constexpr auto isnan =
      25             :     overloaded{special_value_args<Eigen::Vector3d>,
      26          91 :                [](const auto x) {
      27             :                  using numeric::isnan;
      28          91 :                  return isnan(x);
      29             :                },
      30          39 :                [](const units::Unit &) { return units::none; }};
      31             : 
      32             : constexpr auto isinf =
      33             :     overloaded{special_value_args<Eigen::Vector3d>,
      34          21 :                [](const auto x) {
      35             :                  using numeric::isinf;
      36          21 :                  return isinf(x);
      37             :                },
      38          17 :                [](const units::Unit &) { return units::none; }};
      39             : 
      40             : constexpr auto isfinite =
      41             :     overloaded{special_value_args<Eigen::Vector3d>,
      42         870 :                [](const auto x) {
      43             :                  using numeric::isfinite;
      44         870 :                  return isfinite(x);
      45             :                },
      46         175 :                [](const units::Unit &) { return units::none; }};
      47             : 
      48             : namespace detail {
      49          19 : template <typename T> auto isposinf(T x) {
      50          19 :   return numeric::isinf(x) && !numeric::signbit(x);
      51             : }
      52             : 
      53          19 : template <typename T> auto isneginf(T x) {
      54          19 :   return numeric::isinf(x) && numeric::signbit(x);
      55             : }
      56             : } // namespace detail
      57             : 
      58             : constexpr auto isposinf =
      59             :     overloaded{special_value_args<>,
      60          21 :                [](const auto x) {
      61             :                  using detail::isposinf;
      62          21 :                  return isposinf(x);
      63             :                },
      64          17 :                [](const units::Unit &) { return units::none; }};
      65             : 
      66             : constexpr auto isneginf =
      67             :     overloaded{special_value_args<>,
      68          21 :                [](const auto x) {
      69             :                  using detail::isneginf;
      70          21 :                  return isneginf(x);
      71             :                },
      72          17 :                [](const units::Unit &) { return units::none; }};
      73             : 
      74             : constexpr auto replace_special = overloaded{
      75             :     arg_list<double, float>, transform_flags::expect_all_or_none_have_variance,
      76             :     transform_flags::force_variance_broadcast,
      77           5 :     [](const units::Unit &x, const units::Unit &repl) {
      78           5 :       expect::equals(x, repl);
      79           5 :       return x;
      80             :     }};
      81             : 
      82             : constexpr auto replace_special_out_arg = overloaded{
      83             :     arg_list<double, float>, transform_flags::expect_all_or_none_have_variance,
      84             :     transform_flags::force_variance_broadcast,
      85          16 :     [](units::Unit &a, const units::Unit &b, const units::Unit &repl) {
      86          16 :       expect::equals(b, repl);
      87          16 :       a = b;
      88          16 :     }};
      89             : 
      90             : constexpr auto nan_to_num =
      91           6 :     overloaded{replace_special, [](const auto x, const auto &repl) {
      92             :                  using std::isnan;
      93           6 :                  return isnan(x) ? repl : x;
      94             :                }};
      95             : 
      96             : constexpr auto nan_to_num_out_arg = overloaded{
      97          16 :     replace_special_out_arg, [](auto &x, const auto y, const auto &repl) {
      98             :       using numeric::isnan;
      99          16 :       x = isnan(y) ? repl : y;
     100          16 :     }};
     101             : 
     102             : constexpr auto positive_inf_to_num =
     103           3 :     overloaded{replace_special, [](const auto &x, const auto &repl) {
     104             :                  if constexpr (is_ValueAndVariance_v<std::decay_t<decltype(x)>>)
     105           0 :                    return isinf(x) && x.value > 0 ? repl : x;
     106             :                  else
     107           3 :                    return numeric::isinf(x) && x > 0 ? repl : x;
     108             :                }};
     109             : 
     110             : constexpr auto positive_inf_to_num_out_arg = overloaded{
     111          10 :     replace_special_out_arg, [](auto &x, const auto &y, const auto &repl) {
     112             :       if constexpr (is_ValueAndVariance_v<std::decay_t<decltype(y)>>)
     113           0 :         x = isinf(y) && y.value > 0 ? repl : y;
     114             :       else
     115          10 :         x = numeric::isinf(y) && y > 0 ? repl : y;
     116          10 :     }};
     117             : 
     118             : constexpr auto negative_inf_to_num =
     119           3 :     overloaded{replace_special, [](const auto &x, const auto &repl) {
     120             :                  if constexpr (is_ValueAndVariance_v<std::decay_t<decltype(x)>>)
     121           0 :                    return isinf(x) && x.value < 0 ? repl : x;
     122             :                  else
     123           3 :                    return numeric::isinf(x) && x < 0 ? repl : x;
     124             :                }};
     125             : 
     126             : constexpr auto negative_inf_to_num_out_arg = overloaded{
     127          10 :     replace_special_out_arg, [](auto &x, const auto &y, const auto &repl) {
     128             :       if constexpr (is_ValueAndVariance_v<std::decay_t<decltype(y)>>)
     129           0 :         x = isinf(y) && y.value < 0 ? repl : y;
     130             :       else
     131          10 :         x = numeric::isinf(y) && y < 0 ? repl : y;
     132          10 :     }};
     133             : 
     134             : } // namespace scipp::core::element