61 lines · plain
1// (C) Copyright Matt Borland 2022.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#include <cmath>7#include <iterator>8#include <utility>9#include <algorithm>10#include <type_traits>11#include <initializer_list>12#include <boost/math/special_functions/logaddexp.hpp>13 14namespace boost { namespace math {15 16// https://nhigham.com/2021/01/05/what-is-the-log-sum-exp-function/17// See equation (#)18template <typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>19Real logsumexp(ForwardIterator first, ForwardIterator last)20{21 using std::exp;22 using std::log1p;23 24 const auto elem = std::max_element(first, last);25 const Real max_val = *elem;26 27 Real arg = 0;28 while (first != last)29 {30 if (first != elem) 31 {32 arg += exp(*first - max_val);33 }34 35 ++first;36 }37 38 return max_val + log1p(arg);39}40 41template <typename Container, typename Real = typename Container::value_type>42inline Real logsumexp(const Container& c)43{44 return logsumexp(std::begin(c), std::end(c));45}46 47template <typename... Args, typename Real = typename std::common_type<Args...>::type, 48 typename std::enable_if<std::is_floating_point<Real>::value, bool>::type = true>49inline Real logsumexp(Args&& ...args)50{51 std::initializer_list<Real> list {std::forward<Args>(args)...};52 53 if(list.size() == 2)54 {55 return logaddexp(*list.begin(), *std::next(list.begin()));56 }57 return logsumexp(list.begin(), list.end());58}59 60}} // Namespace boost::math61