From eddefa31b329b9b4eab3cb2de912ab25ed832c67 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 5 Dec 2024 13:16:26 -0700 Subject: [PATCH 1/3] add robust utils header --- doc/sphinx/src/using.rst | 59 ++++++++++++++++++++++++++ ports-of-call/robust_utils.hpp | 75 ++++++++++++++++++++++++++++++++++ 2 files changed, 134 insertions(+) create mode 100644 ports-of-call/robust_utils.hpp diff --git a/doc/sphinx/src/using.rst b/doc/sphinx/src/using.rst index 5abf037b..9781d052 100644 --- a/doc/sphinx/src/using.rst +++ b/doc/sphinx/src/using.rst @@ -128,6 +128,65 @@ Please note that none of these functions are thread or MPI aware. In a parallel as appropriate. +robust_utils.hpp +^^^^^^^^^^^^^^^^^^^ + +``robust_utils.hpp`` contains small utility functions for numerical +robustness, especially around floating point numbers. The available +functionality is contained in the namespace ``PortsOfCall::Robust`` and includes: + +* ``constexpr auto SMALL()`` returns a small number of type ``T``. +* ``constexpr auto EPS()`` returns a value of type ``T`` close to machine epsilon. +* ``constexpr auto min_exp_arg()`` returns the smallest of type ``T`` safe value to pass into an exponent. +* ``constexpr auto max_exp_exp_arg()`` returns the max value of type ``T`` to pass into an exponent. +* ``auto make_positive(const T val)`` makes the argument of type ``T`` positive. + +where here all functionality is templated on type ``T`` and marked +with ``PORTABLE_INLINE_FUNCTION``. The default type ``T`` is always +``Real``. + +The function + +.. code-block:: cpp + + template + PORTABLE_FORCEINLINE_FUNCTION + Real make_bounded(const T val, const T vmin, const T vmax); + +bounds ``val`` between ``vmin`` and ``vmax``. + +The function + +.. code-block:: cpp + + template + PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val); + +returns the sign of a quantity. + +The function + +.. code-block:: cpp + + template + PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) + +computes the ratio :math:`A/B` but in a way robust to divide by zero +errors. If both :math:`A` and :math:`B` are zero, this function will +return 0. + +The function + +.. code-block:: cpp + + template + PORTABLE_FORCEINLINE_FUNCTION T safe_arg_exp(const T &x) + +returns exponentiation in such a way that avoids floating point +exceptions. For very large negative inputs, it returns 0. For very +large positive ones, it returns +``std::numeric_limits::infinity()``. + macros_arrays.hpp ^^^^^^^^^^^^^^^^^^^ diff --git a/ports-of-call/robust_utils.hpp b/ports-of-call/robust_utils.hpp new file mode 100644 index 00000000..d5cd0272 --- /dev/null +++ b/ports-of-call/robust_utils.hpp @@ -0,0 +1,75 @@ +//------------------------------------------------------------------------------ +// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef PORTS_OF_CALL_ROBUST_UTILS_HPP_ +#define PORTS_OF_CALL_ROBUST_UTILS_HPP_ + +#include +#include +#include + +namespace PortsOfCall { +namespace Robust { + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() { + return 10 * std::numeric_limits::min(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() { + return 10 * std::numeric_limits::epsilon(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr T min_exp_arg() { + return (std::numeric_limits::min_exponent - 1) * M_LN2; +} + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr T max_exp_arg() { + return std::numeric_limits::max_exponent * M_LN2; +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) { + return std::max(val, EPS()); +} + +template +PORTABLE_FORCEINLINE_FUNCTION Real make_bounded(const T val, const T vmin, const T vmax) { + return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); +} + +template +PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { + return (T(0) <= val) - (val < T(0)); +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { + return a / (b + sgn(b) * SMALL()); +} + +template +PORTABLE_FORCEINLINE_FUNCTION T safe_arg_exp(const T &x) { + return x < min_exp_arg() ? 0.0 + : x > max_exp_arg() ? std::numeric_limits::infinity() + : std::exp(x); +} + +} // namespace Robust +} // namespace PortsOfCall + +#endif // PORTS_OF_CALL_ROBUST_UTILS_HPP_ From d23c5d931931bd6c01601c7edf5a1a74eb5231fc Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 5 Dec 2024 13:33:44 -0700 Subject: [PATCH 2/3] add note about sgn --- doc/sphinx/src/using.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/src/using.rst b/doc/sphinx/src/using.rst index 9781d052..648efd6b 100644 --- a/doc/sphinx/src/using.rst +++ b/doc/sphinx/src/using.rst @@ -162,7 +162,12 @@ The function template PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val); -returns the sign of a quantity. +returns the sign of a quantity ``val``. + +.. note:: + + Note this implementation **never** returns zero. It **always** + returns :math:`\pm 1`. The function From 08e3dfa3a38306d2430d0e04caa1d5c6324486ad Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 5 Dec 2024 20:55:51 -0700 Subject: [PATCH 3/3] Brendan's comments --- doc/sphinx/src/using.rst | 14 ++++++++------ ports-of-call/robust_utils.hpp | 2 ++ 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/sphinx/src/using.rst b/doc/sphinx/src/using.rst index 648efd6b..729830c3 100644 --- a/doc/sphinx/src/using.rst +++ b/doc/sphinx/src/using.rst @@ -137,8 +137,8 @@ functionality is contained in the namespace ``PortsOfCall::Robust`` and includes * ``constexpr auto SMALL()`` returns a small number of type ``T``. * ``constexpr auto EPS()`` returns a value of type ``T`` close to machine epsilon. -* ``constexpr auto min_exp_arg()`` returns the smallest of type ``T`` safe value to pass into an exponent. -* ``constexpr auto max_exp_exp_arg()`` returns the max value of type ``T`` to pass into an exponent. +* ``constexpr auto min_exp_arg()`` returns the smallest safe value of type ``T`` to pass into an exponent. +* ``constexpr auto max_exp_exp_arg()`` returns the max safe value of type ``T`` to pass into an exponent. * ``auto make_positive(const T val)`` makes the argument of type ``T`` positive. where here all functionality is templated on type ``T`` and marked @@ -153,7 +153,8 @@ The function PORTABLE_FORCEINLINE_FUNCTION Real make_bounded(const T val, const T vmin, const T vmax); -bounds ``val`` between ``vmin`` and ``vmax``. +bounds ``val`` between ``vmin`` and ``vmax``, exclusive. Note this is +slightly different than ``std::clamp``, which uses inclusive bounds. The function @@ -176,9 +177,10 @@ The function template PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) -computes the ratio :math:`A/B` but in a way robust to divide by zero -errors. If both :math:`A` and :math:`B` are zero, this function will -return 0. +computes the ratio :math:`A/B` but in a way robust to 0/0 errors. If +both :math:`A` and :math:`B` are zero, this function will return 0. If +:math:`|A| > 0` and :math:`B=0`, then it will return a very large, +possibly (but not guaranteed to be) infinite number. The function diff --git a/ports-of-call/robust_utils.hpp b/ports-of-call/robust_utils.hpp index d5cd0272..d95bc105 100644 --- a/ports-of-call/robust_utils.hpp +++ b/ports-of-call/robust_utils.hpp @@ -59,6 +59,8 @@ PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { template PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { + B mask = static_cast(std::abs(b) < SMALL()); + B denom = mask * sgn(b) * SMALL() + (1 - mask) * b; return a / (b + sgn(b) * SMALL()); }