apfel/ome/piecewise.h

0001 /******************************************************************************
0002  * This file is part of libome                                                *
0003  * Copyright (C) 2025 Arnd Behring, Kay Schoenwald                            *
0004  * SPDX-License-Identifier: GPL-3.0-or-later                                  *
0005  ******************************************************************************/
0006
0007 /**
0008  * \file
0009  * \brief Callable for modelling piecewise functions
0010  */
0011
0012 #ifndef LIBOME_PIECEWISE_H
0013 #define LIBOME_PIECEWISE_H
0014
0015 #include <stdexcept>
0016 #include <vector>
0017 #include <algorithm>
0018
0019 namespace apfel
0020 {
0021   namespace ome
0022   {
0023     /**
0024      * \brief Exception: definition of piecewise function is inconsistent
0025      */
0026     class inconsistent_piecewise : public std::logic_error
0027     {
0028     public:
0029       // Inherit constructors from base class
0030       using std::logic_error::logic_error;
0031     };
0032
0033     /**
0034      * \brief One-dimensional piecewise function
0035      *
0036      * \details
0037      * This class models a piecewise function defined on a subdivision of the
0038      * real axis into non-overlapping, contiguous itervals. Upon evaluation, the
0039      * first argument is used to look up on which interval the function is
0040      * evaluated and the corresponding function, stored as a callable is
0041      * evaluated.
0042      *
0043      * \tparam Tnum Numerical type for the first argument upon evaluation, the
0044      *         return type
0045      * \tparam Tfunc Callable type for the functions on each interval
0046      * \tparam Trest Pack of types for the remaining arguments of the wrapped
0047      *         callable
0048      */
0049     template<typename Tnum, typename Tfunc, typename... Trest>
0050     class piecewise
0051     {
0052     public:
0053       /// Type alias for the numerical type template parameter
0054       using numeric_type = Tnum;
0055       /// Type alias for the callable type template parameter
0056       using element_type = Tfunc;
0057
0058       /**
0059        * \brief Default constructor
0060        *
0061        * \details
0062        * Constructs a function with only one interval and default constructed
0063        * callable.
0064        */
0065       piecewise()
0066         : boundaries_(), pieces_{{element_type()}} {};
0067
0068       /**
0069        * \brief Construct piecewise function from interval boundaries and
0070        *        callables on those intervals.
0071        *
0072        * \details
0073        * The vector boundaries must be sorted in ascending order and have
0074        * exactly one element less that the vector pieces. The first interval
0075        * starts at negative infinity and extends up to, but excluding the first
0076        * boundary \f$b_0\f$. Each further interval is given by
0077        * \f$[b_i,b_{i+1})\f$ and the last interval extends to positive infinity.
0078        *
0079        * \param boundaries Vector of interval boundaries
0080        * \param pieces Vector of callable objects representing the function
0081        *        on each interval
0082        */
0083       piecewise(const std::vector<numeric_type> boundaries, const std::vector<element_type> pieces)
0084         : boundaries_(boundaries), pieces_(pieces)
0085       {
0086         if(boundaries_.size()+1 != pieces_.size())
0087           throw(inconsistent_piecewise("The must be one more function than boundaries in a piecewise function"));
0088
0089         if(!std::is_sorted(boundaries_.begin(), boundaries_.end()))
0090           throw(inconsistent_piecewise("Interval boundaries must be specified in ascending order"));
0091       };
0092
0093       /**
0094        * \brief Helper to allow construction directly from intializer lists
0095        *
0096        * \details
0097        * See \ref piecewise(const std::vector<numeric_type>,
0098        * const std::vector<element_type>) for details.
0099        */
0100       piecewise(std::initializer_list<numeric_type> boundaries, std::initializer_list<element_type> pieces)
0101         : piecewise(std::vector<numeric_type>(boundaries), std::vector<element_type>(pieces))
0102       {};
0103
0104       /**
0105        * \brief Evaluate piecewise function
0106        *
0107        * \param x First argument that decides on which interval we have to evaluate
0108        * \param rest Function parameter pack that is passed on to the wrapped
0109        *        callable unchanged
0110        *
0111        * \return The value of the piecewise function at the given point
0112        */
0113       numeric_type operator()(numeric_type x, Trest... rest) const
0114       {
0115         if(pieces_.empty())
0116           return(static_cast<numeric_type>(0));
0117
0118         auto piece = pieces_.cbegin();
0119         for(auto boundary = boundaries_.cbegin(); boundary != boundaries_.cend(); ++piece, ++boundary)
0120           {
0121             if(x < *boundary)
0122               {
0123                 return((*piece)(x, rest...));
0124               }
0125           }
0126         return((*piece)(x, rest...));
0127       };
0128
0129     private:
0130       std::vector<numeric_type> boundaries_;
0131       std::vector<element_type> pieces_;
0132     };
0133   };
0134 }
0135 #endif