diff --git a/backend/src/Utilities/NumericAt.hpp b/backend/src/Utilities/NumericAt.hpp new file mode 100644 index 00000000..68ff6b4c --- /dev/null +++ b/backend/src/Utilities/NumericAt.hpp @@ -0,0 +1,429 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include +#include +#include +#include +#include + +#include "Utilities/Invoke.hpp" +#include "Utilities/MakeArray.hpp" +#include "Utilities/NonCopyable.hpp" + +namespace elastica { + + namespace detail { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Base type for memoizing a compact set of integral indices + // \ingroup utils + */ + struct numeric_at_base_t : NonCopyable { + //**Constructors********************************************************** + /*!\name Constructors */ + //@{ + + //************************************************************************ + /*!\brief Default constructor. + */ + numeric_at_base_t() = default; + //************************************************************************ + + //************************************************************************ + /*!\brief Move constructor. + */ + numeric_at_base_t(numeric_at_base_t&&) = default; + //@} + //************************************************************************ + + //**Destructor************************************************************ + /*!\name Destructor */ + //@{ + virtual ~numeric_at_base_t() = default; + //@} + //************************************************************************ + + //**Utility*************************************************************** + /*!\name Utility methods */ + //@{ + + //************************************************************************ + /*!\brief Gets the indices buffer + */ + virtual std::size_t* data() noexcept = 0; + //************************************************************************ + + //************************************************************************ + /*!\brief Gets the number of indices + */ + virtual std::size_t size() const noexcept = 0; + //************************************************************************ + + //@} + //************************************************************************ + }; + /*! \endcond */ + //************************************************************************** + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief A compact set of integral indices known at compile-time + // \ingroup utils + */ + template + struct numeric_at_t final : numeric_at_base_t { + //**Type definitions****************************************************** + //! Type of array for integral indices + using type = std::array; + //************************************************************************ + + //**Constructors********************************************************** + /*!\name Constructors */ + //@{ + + //************************************************************************ + /*!\brief Array constructor. + // + // \param a array of integer indices + */ + explicit constexpr numeric_at_t(type a) noexcept : args{a} {}; + //************************************************************************ + + //************************************************************************ + /*!\brief Move constructor. + */ + numeric_at_t(numeric_at_t&& other) noexcept + : args{std::move(other.args)} {}; + //************************************************************************ + + //@} + //************************************************************************ + + //**Destructor************************************************************ + /*!\name Destructor */ + //@{ + ~numeric_at_t() override = default; + //@} + //********************************************************************** + + //**Utility*************************************************************** + /*!\name Utility methods */ + //@{ + + //************************************************************************ + /*!\brief Gets the indices buffer + */ + std::size_t* data() noexcept override { return args.data(); } + //************************************************************************ + + //************************************************************************ + /*!\brief Gets the number of indices + */ + std::size_t size() const noexcept override { return N; } + //************************************************************************ + + //@} + //************************************************************************ + + //**Member variables****************************************************** + //! The indices + type args; + //************************************************************************ + }; + /*! \endcond */ + //************************************************************************** + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief A compact set of integral indices only known at runtime + // \ingroup utils + */ + struct runtime_numeric_at_t final : numeric_at_base_t { + //**Type definitions****************************************************** + //! Type of storage for integral indices + using type = std::vector; + //************************************************************************ + + //**Constructors********************************************************** + /*!\name Constructors */ + //@{ + + //************************************************************************ + /*!\brief vector constructor. + // + // \param a vector of integer indices + */ + explicit runtime_numeric_at_t(type&& a) : args(std::move(a)){}; + //************************************************************************ + + //************************************************************************ + /*!\brief Move constructor. + */ + runtime_numeric_at_t(runtime_numeric_at_t&& other) noexcept + : args(std::move(other.args)){}; + //************************************************************************ + + //@} + //************************************************************************ + + //**Destructor************************************************************ + /*!\name Destructor */ + //@{ + ~runtime_numeric_at_t() override = default; + //@} + //********************************************************************** + + //**Utility*************************************************************** + /*!\name Utility methods */ + //@{ + + //************************************************************************ + /*!\brief Reserves space in the underlying storage + * + * \param size size to reserve + */ + void reserve(std::size_t size) { args.reserve(size); } + //************************************************************************ + + //************************************************************************ + /*!\brief Gets the indices buffer + */ + std::size_t* data() noexcept override { return args.data(); } + //************************************************************************ + + //************************************************************************ + /*!\brief Gets the number of indices + */ + std::size_t size() const noexcept override { return args.size(); } + //@} + //************************************************************************ + + //**Member variables****************************************************** + //! The runtime indices + type args{}; + //************************************************************************ + }; + /*! \endcond */ + //************************************************************************** + + } // namespace detail + + //============================================================================ + // + // UTILITY FUNCTIONS + // + //============================================================================ + + //**************************************************************************** + /*!\brief Constructs a compact set of integral indices from parameters + * \ingroup utils + * + * \example + * \snippet Test_NumericAt.cpp numeric_at_example + * + * + * \param x An integral index + * \param ...indices Other integral indices + */ + template + inline constexpr auto numeric_at(std::size_t x, Ints... ints) + -> detail::numeric_at_t { + return detail::numeric_at_t{make_array(x, ints...)}; + // sizeof...(Ints) > 0 ? make_array(x, ints...) : make_array<1UL, + // int>(x); return detail::at_t{ + // make_array(x, ints...)}; + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Sorts and checks duplicates within a compact set of integral indices + * \ingroup utils + * + * \details + * Sorts and checks for duplicates. If duplicates exist, an exception is + * raised + * + * \example + * \snippet Test_NumericAt.cpp scd_example + * + * \param ats Numeric indices to be sorted and checked + * \return ats + */ + template + inline auto sort_and_check_duplicates(detail::numeric_at_t&& ats) + -> detail::numeric_at_t { + std::sort(std::begin(ats.args), std::end(ats.args)); + if (std::adjacent_find(std::begin(ats.args), std::end(ats.args)) != + std::end(ats.args)) { + throw std::logic_error("Repeated elements in the set"); + } + return std::move(ats); + } + //**************************************************************************** + + //**************************************************************************** + /*!\brief Sorts and checks duplicates within a compact set of integral indices + * \ingroup utils + * + * \details + * Runtime equivalent of sort_and_check_duplicates() + * + * \param ats A unique pointer to numeric indices to be sorted and checked + * \return ats + */ + template + inline auto sort_and_check_duplicates( + std::unique_ptr&& ats) + -> std::unique_ptr { + std::sort(ats->data(), ats->data() + ats->size()); + if (std::adjacent_find(ats->data(), ats->data() + ats->size()) != + (ats->data() + ats->size())) { + throw std::logic_error("Repeated elements in the set"); + } + return std::move(ats); + } + //**************************************************************************** + + /////////// The following functions are no longer used /////////////////////// + //**************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /* + template + auto sanitize(detail::numeric_at_t const& index_set, + std::size_t n_dofs) noexcept -> detail::numeric_at_t { + // there may be overflow here, but its less likely + auto bound = [nd = static_cast(n_dofs)](int x) -> int { + return x < 0 ? (nd + x) : x; + }; + return index_apply([&bound, &args = index_set.args](auto... Is) { + return numeric_at(bound(std::get(args))...); + }); + } + + template + auto sanitize(detail::numeric_at_t&& index_set, + std::size_t n_dofs) noexcept -> detail::numeric_at_t { + auto bound = [n_dofs = static_cast(n_dofs)](int& x) -> void { + x = x < 0 ? (n_dofs + x) : x; + }; + index_apply([&bound, &args = index_set.args](auto... Is) { + return EXPAND_PACK_LEFT_TO_RIGHT(bound(std::get(args))); + }); + // force move construction, C++17 should fix this by automatic NRVO + return std::move(index_set); + } + + template + auto sanitize(std::unique_ptr&& index_set, + std::size_t n_dofs) noexcept + -> std::unique_ptr { + auto bound = [n_dofs = static_cast(n_dofs)](int& x) -> void { + x = x < 0 ? (n_dofs + x) : x; + }; + std::for_each(index_set->data(), index_set->data() + index_set->size(), + bound); + // force move construction, C++17 should fix this by automatic NRVO + return std::move(index_set); + } + + template + auto check_range(detail::numeric_at_t&& sanitized, std::size_t n_dofs) + -> detail::numeric_at_t { + // < 0 or == SIZE_MAX, then report error + // there may be overflow here, but its less likely + std::ostringstream s_stream; + auto checker = [&s_stream, + n_dofs = static_cast(n_dofs)](int x) -> bool { + bool failed = (x < 0) || (x >= n_dofs); + if (failed) + s_stream << x << ", "; + return failed; + }; + + bool failed = false; + /////// compile-time equivalent + for (const auto& x : sanitized.args) { + failed |= checker(x); + } + + if (failed) { + s_stream << "as the number of degrees of freedom is only " << n_dofs; + throw std::range_error(s_stream.str()); + } + + // force move construction, C++17 should fix this by automatic NRVO + return std::move(sanitized); + } + + // template + // auto check_range(detail::at_t&& sanitized, std::size_t n_dofs) + // -> detail::at_t { + template + auto check_range(std::unique_ptr&& sanitized, + std::size_t n_dofs) + -> std::unique_ptr { + // < 0 or == SIZE_MAX, then report error + // there may be overflow here, but its less likely + std::ostringstream s_stream; + auto checker = [&s_stream, + n_dofs = static_cast(n_dofs)](int const x) -> bool { + bool failed = (x < 0) || (x >= n_dofs); + if (failed) + s_stream << x << ", "; + return failed; + }; + + // runtime replacement + const bool failed = [&]() { + bool int_failed = false; + auto accumulator = [&int_failed, &checker](int const x) { + int_failed |= checker(x); + }; + std::for_each(sanitized->data(), sanitized->data() + sanitized->size(), + accumulator); + return int_failed; + }(); + + if (failed) { + s_stream << "as the number of degrees of freedom is only " << n_dofs; + throw std::range_error(s_stream.str()); + } + + // force move construction, C++17 should fix this by automatic NRVO + return std::move(sanitized); + } + + template + auto sanitize_and_check(detail::numeric_at_t const& index_set, + std::size_t n_dofs) -> detail::numeric_at_t { + return check_range(sanitize(index_set, n_dofs), n_dofs); + } + + template + auto sanitize_and_check(detail::numeric_at_t&& index_set, + std::size_t n_dofs) -> detail::numeric_at_t { + return check_range(sanitize(std::move(index_set), n_dofs), n_dofs); + } + + // runtime support + template + auto sanitize_and_check( + std::unique_ptr&& index_set, + std::size_t n_dofs) -> std::unique_ptr { + return check_range(sanitize(std::move(index_set), n_dofs), n_dofs); + } + */ + /*! \endcond */ + //**************************************************************************** + +} // namespace elastica diff --git a/backend/src/Utilities/NumericRange.hpp b/backend/src/Utilities/NumericRange.hpp new file mode 100644 index 00000000..7c86bdb1 --- /dev/null +++ b/backend/src/Utilities/NumericRange.hpp @@ -0,0 +1,202 @@ +#pragma once + +//****************************************************************************** +// Includes +//****************************************************************************** +#include + +#include "ErrorHandling/Assert.hpp" +#include "Utilities/NonCopyable.hpp" +#include "Utilities/Requires.hpp" +#include "Utilities/TypeTraits/Cpp17.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +namespace elastica { + + namespace detail { + + //========================================================================== + // + // CLASS DEFINITION + // + //========================================================================== + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief A numeric range of indices + // \ingroup utils + // + // \details + // Range models a range of indices on which one can iterate upon or evaluate + // in case of floating point ranges. + // + // \tparam F The arithmetic type of the range + */ + template + struct NumericRangeImpl : public NonCopyable { + //**Constructors********************************************************** + /*!\name Constructors */ + //@{ + + //************************************************************************ + /*!\brief Default constructor. + */ + constexpr NumericRangeImpl(F st, F so, Step se) noexcept( + std::is_nothrow_move_constructible::value) + : start_(std::move(st)), stop_(std::move(so)), step_(std::move(se)){}; + //************************************************************************ + + //************************************************************************ + /*!\brief Move constructor. + */ + constexpr NumericRangeImpl(NumericRangeImpl&& other) noexcept( + std::is_nothrow_move_constructible::value) + : start_(std::move(other.start_)), + stop_(std::move(other.stop_)), + step_(std::move(other.step_)) {} + //************************************************************************ + + //@} + //************************************************************************ + + //**Destructor************************************************************ + /*!\name Destructor */ + //@{ + ~NumericRangeImpl() = default; + //@} + //************************************************************************ + + //**Member variables****************************************************** + // makes little sense to modify a range after creation, hence declared + // const + //! The start index of range + F const start_; + //! The stop index of range + F const stop_; + //! The step of range + Step const step_; + //************************************************************************ + }; + /*! \endcond */ + //************************************************************************** + + //========================================================================== + // + // ALIAS DECLARATIONS + // + //========================================================================== + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Alias declaration for an integer NumericRangeImpl. + // \ingroup utils + */ + using NumericRange = NumericRangeImpl; + /*! \endcond */ + //************************************************************************** + + } // namespace detail + + namespace tt { + //************************************************************************** + /*!\brief Checks whether a type models \elastica NumericRange + // \ingroup utils type_traits + // + // \example + // \snippet Test_NumericRange.cpp is_numeric_range_example + */ + template + struct IsNumericRange : ::tt::is_a<::elastica::detail::NumericRangeImpl, T> {}; + //************************************************************************** + } // namespace tt + + namespace detail { + + //************************************************************************** + /*! \cond ELASTICA_INTERNAL */ + /*!\brief Idiomatic way to construct a half-open numeric range + // + // \details + // The range is closed from left and open from right and models + // \f$ [a, b) \f$, where \f$ a \f$ is the start and \f$ b \f$ is the stop. + // + // \example + // \snippet Test_NumericRange.cpp numeric_range_example + // + // \param start Start of the range + // \param stop Stop of the range + // \param step Step of the range + */ + template + inline constexpr auto numeric_range( + F start, F stop, + Step step = 1) noexcept(noexcept(NumericRangeImpl{ + start, stop, step})) -> NumericRangeImpl { + return detail::NumericRangeImpl{start, stop, step}; + } + /*! \endcond */ + //************************************************************************** + + } // namespace detail + + //============================================================================ + // + // UTILITY FUNCTIONS + // + //============================================================================ + + //**************************************************************************** + /*!\brief Idiomatic way to construct a half-open numeric range + // + // \details + // The range is closed from left and open from right and models + // \f$ [a, b) \f$, where \f$ a \f$ is the start and \f$ b \f$ is the stop. + // + // \example + // \snippet Test_NumericRange.cpp numeric_range_example + // + // \param start Start of the range + // \param stop Stop of the range + // \param step Step of the range + // + // \note + // - while the name is more verbose, it differentiates itself from interface:: + // range(). + // - The parameters are restricted to C++ integral but non-boolean types + */ + template + inline decltype(auto) numeric_range(F start, F stop, Step step = 1) { + static_assert( + std::is_integral::value && not std::is_same::value, + "Type must be an arithmetic, non-boolean type"); + static_assert( + std::is_integral::value && not std::is_same::value, + "Step must be an arithmetic, non-boolean type"); + + ELASTICA_ASSERT(start != stop, "Range start and end cannot be equal"); + ELASTICA_ASSERT(step != Step(0), "Step size cannot be zero"); + const bool has_positive_step(step > Step(0)); + ELASTICA_ASSERT(((has_positive_step && start <= stop) || + (!has_positive_step && start >= stop)), + "Step moves away from end"); + // switch arguments here so that we can work only with size_t in forcing + // etc. + return has_positive_step + ? detail::numeric_range(start, stop, static_cast(step)) + : detail::numeric_range(stop, start, static_cast(-step)); + } + //**************************************************************************** + + // template + // auto sanitize(detail::NumericRangeImpl const& index_range, std::size_t + // n_dofs) noexcept + // -> detail::NumericRangeImpl{ + // // there may be overflow here, but its less likely + // auto bound = [nd = static_cast(n_dofs)](int x) -> int { + // return x < 0 ? (nd + x) : x; + // }; + // return index_apply([&bound, &args = index_set.args](auto... Is) { + // return at(bound(std::get(args))...); + // }); + // } +} // namespace elastica