math/special_functions/logsumexp.hpp

0001 //  (C) Copyright Matt Borland 2022.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #include <cmath>
0007 #include <iterator>
0008 #include <utility>
0009 #include <algorithm>
0010 #include <type_traits>
0011 #include <initializer_list>
0012 #include <boost/math/special_functions/logaddexp.hpp>
0013
0014 namespace boost { namespace math {
0015
0016 // https://nhigham.com/2021/01/05/what-is-the-log-sum-exp-function/
0017 // See equation (#)
0018 template <typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>
0019 Real logsumexp(ForwardIterator first, ForwardIterator last)
0020 {
0021     using std::exp;
0022     using std::log1p;
0023
0024     const auto elem = std::max_element(first, last);
0025     const Real max_val = *elem;
0026
0027     Real arg = 0;
0028     while (first != last)
0029     {
0030         if (first != elem)
0031         {
0032             arg += exp(*first - max_val);
0033         }
0034
0035         ++first;
0036     }
0037
0038     return max_val + log1p(arg);
0039 }
0040
0041 template <typename Container, typename Real = typename Container::value_type>
0042 inline Real logsumexp(const Container& c)
0043 {
0044     return logsumexp(std::begin(c), std::end(c));
0045 }
0046
0047 template <typename... Args, typename Real = typename std::common_type<Args...>::type,
0048           typename std::enable_if<std::is_floating_point<Real>::value, bool>::type = true>
0049 inline Real logsumexp(Args&& ...args)
0050 {
0051     std::initializer_list<Real> list {std::forward<Args>(args)...};
0052
0053     if(list.size() == 2)
0054     {
0055         return logaddexp(*list.begin(), *std::next(list.begin()));
0056     }
0057     return logsumexp(list.begin(), list.end());
0058 }
0059
0060 }} // Namespace boost::math