mirror of
https://github.com/boostorg/math.git
synced 2025-05-11 21:33:52 +00:00
Implement logaddexp (#763)
* Implement logaddexp * Disable test for ASAN * Implement logsumexp * Add performance file and include results in the docs * Address review comments * Simplify overflow test and comply with min/max guidelines * Minor cleanup * FIxes to comments and docs [ci skip] * Return status code. Co-authored-by: Nick Thompson <nathompson7@protonmail.com>
This commit is contained in:
Matt Borland
GitHub
parent
1153427c3a
commit
e8c40e309c
124
doc/sf/logaddexp.qbk
Normal file
124
doc/sf/logaddexp.qbk
Normal file
@ -0,0 +1,124 @@
|
||||
[/
|
||||
Copyright Matt Borland 2022
|
||||
Distributed under the Boost Software License, Version 1.0.
|
||||
(See accompanying file LICENSE_1_0.txt or copy at
|
||||
http://www.boost.org/LICENSE_1_0.txt).
|
||||
]
|
||||
|
||||
[section:logaddexp logaddexp and logsumexp]
|
||||
|
||||
[h4 Synopsis]
|
||||
|
||||
#include <boost/math/special_functions/logaddexp.hpp>
|
||||
|
||||
namespace boost { namespace math {
|
||||
|
||||
template <typename Real>
|
||||
Real logaddexp (Real x1, Real x2);
|
||||
|
||||
}} // namespaces
|
||||
|
||||
#include <boost/math/special_functions/logsumexp.hpp>
|
||||
|
||||
namespace boost { namespace math {
|
||||
|
||||
template <typename ForwardIterator, typename Real = std::iterator_traits<ForwardIterator>::value_type>
|
||||
Real logsumexp (ForwardIterator first, ForwardIterator last);
|
||||
|
||||
template <typename ForwardContainer, typename Real = ForwardContainer::value_type>
|
||||
Real logsumexp (const ForwardContainer& c);
|
||||
|
||||
template <typename... Args, typename Real = typename std::common_type<Args...>::type>
|
||||
Real logsumexp(Args&& ...args)
|
||||
|
||||
}} // namespace
|
||||
|
||||
|
||||
The function `logaddexp(x1, x2)` computes log(exp(x1) + exp(x2)).
|
||||
The function `logsumexp(x1, x2, ...)` is generalized to compute log(exp(x1) + exp(x2) + ...).
|
||||
This function is useful in statistics where the calculated probabilities of events may be so small as to exceed the range of normal floating point numbers.
|
||||
In such cases the logarithm of the calculated probability is stored.
|
||||
|
||||
|
||||
The use is
|
||||
|
||||
using std::log;
|
||||
|
||||
double x1 = log(1e-50);
|
||||
double x2 = log(2.5e-50);
|
||||
double x3 = log(3e-50);
|
||||
std::vector<double> x = {x1, x2, x3};
|
||||
|
||||
double probability1 = boost::math::logaddexp(x1, x2);
|
||||
double probability2 = boost::math::logsumexp(x1, x2, x3);
|
||||
double probability3 = boost::math::logsumexp(x);
|
||||
double probability4 = boost::math::logsumexp(x.begin(), x.end());
|
||||
|
||||
Performance Data:
|
||||
|
||||
```
|
||||
Running ./logaddexp_performance
|
||||
Run on Apple M1 Pro
|
||||
CPU Caches:
|
||||
L1 Data 64 KiB (x10)
|
||||
L1 Instruction 128 KiB (x10)
|
||||
L2 Unified 4096 KiB (x5)
|
||||
Load Average: 2.23, 1.89, 1.88
|
||||
-----------------------------------------------------------------------------
|
||||
Benchmark Time CPU Iterations
|
||||
-----------------------------------------------------------------------------
|
||||
logaddexp_performance<float> 1.05 ns 1.05 ns 597983940
|
||||
logaddexp_performance<double> 1.03 ns 1.03 ns 672682369
|
||||
```
|
||||
|
||||
```
|
||||
Running ./logsumexp_performance
|
||||
Run on Apple M1 Pro
|
||||
CPU Caches:
|
||||
L1 Data 64 KiB (x10)
|
||||
L1 Instruction 128 KiB (x10)
|
||||
L2 Unified 4096 KiB (x5)
|
||||
Load Average: 1.56, 1.67, 1.81
|
||||
-----------------------------------------------------------------------------------------------
|
||||
Benchmark Time CPU Iterations
|
||||
-----------------------------------------------------------------------------------------------
|
||||
logsumexp_performance<float>/64/real_time 388 ns 388 ns 1797191
|
||||
logsumexp_performance<float>/128/real_time 761 ns 761 ns 890017
|
||||
logsumexp_performance<float>/256/real_time 1513 ns 1513 ns 460217
|
||||
logsumexp_performance<float>/512/real_time 3026 ns 3026 ns 231454
|
||||
logsumexp_performance<float>/1024/real_time 6043 ns 6043 ns 113901
|
||||
logsumexp_performance<float>/2048/real_time 12084 ns 12084 ns 57590
|
||||
logsumexp_performance<float>/4096/real_time 24240 ns 24240 ns 28835
|
||||
logsumexp_performance<float>/8192/real_time 48326 ns 48323 ns 14478
|
||||
logsumexp_performance<float>/16384/real_time 96645 ns 96642 ns 7181
|
||||
logsumexp_performance<float>/32768/real_time 193351 ns 193351 ns 3607
|
||||
logsumexp_performance<float>/65536/real_time 386537 ns 386538 ns 1810
|
||||
logsumexp_performance<float>/131072/real_time 774055 ns 774056 ns 894
|
||||
logsumexp_performance<float>/262144/real_time 1548913 ns 1548847 ns 451
|
||||
logsumexp_performance<float>/524288/real_time 3092771 ns 3092770 ns 226
|
||||
logsumexp_performance<float>/1048576/real_time 6188087 ns 6188089 ns 112
|
||||
logsumexp_performance<float>/real_time_BigO 5.90 N 5.90 N
|
||||
logsumexp_performance<float>/real_time_RMS 0 % 0 %
|
||||
logsumexp_performance<double>/64/real_time 388 ns 388 ns 1806669
|
||||
logsumexp_performance<double>/128/real_time 770 ns 770 ns 898340
|
||||
logsumexp_performance<double>/256/real_time 1534 ns 1534 ns 454768
|
||||
logsumexp_performance<double>/512/real_time 3063 ns 3063 ns 228057
|
||||
logsumexp_performance<double>/1024/real_time 6126 ns 6126 ns 112667
|
||||
logsumexp_performance<double>/2048/real_time 12243 ns 12243 ns 56963
|
||||
logsumexp_performance<double>/4096/real_time 24476 ns 24476 ns 28485
|
||||
logsumexp_performance<double>/8192/real_time 48979 ns 48978 ns 14215
|
||||
logsumexp_performance<double>/16384/real_time 97929 ns 97929 ns 7070
|
||||
logsumexp_performance<double>/32768/real_time 195826 ns 195826 ns 3560
|
||||
logsumexp_performance<double>/65536/real_time 391855 ns 391835 ns 1787
|
||||
logsumexp_performance<double>/131072/real_time 784119 ns 784110 ns 882
|
||||
logsumexp_performance<double>/262144/real_time 1566408 ns 1566386 ns 446
|
||||
logsumexp_performance<double>/524288/real_time 3151649 ns 3150955 ns 223
|
||||
logsumexp_performance<double>/1048576/real_time 6300578 ns 6299027 ns 110
|
||||
logsumexp_performance<double>/real_time_BigO 6.01 N 6.01 N
|
||||
```
|
||||
|
||||
[heading References]
|
||||
|
||||
* Higham, Nicholas J. [@https://nhigham.com/2021/01/05/what-is-the-log-sum-exp-function/ What is the Log-Sum-Exp Function?]
|
||||
|
||||
[endsect]
|
||||
41
include/boost/math/special_functions/logaddexp.hpp
Normal file
41
include/boost/math/special_functions/logaddexp.hpp
Normal file
@ -0,0 +1,41 @@
|
||||
// (C) Copyright Matt Borland 2022.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <boost/math/special_functions/fpclassify.hpp>
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
namespace boost { namespace math {
|
||||
|
||||
// Calculates log(exp(x1) + exp(x2))
|
||||
template <typename Real>
|
||||
Real logaddexp(Real x1, Real x2) noexcept
|
||||
{
|
||||
using std::log1p;
|
||||
using std::exp;
|
||||
using std::abs;
|
||||
|
||||
// Validate inputs first
|
||||
if (!(boost::math::isfinite)(x1))
|
||||
{
|
||||
return x1;
|
||||
}
|
||||
else if (!(boost::math::isfinite)(x2))
|
||||
{
|
||||
return x2;
|
||||
}
|
||||
|
||||
const Real temp = x1 - x2;
|
||||
|
||||
if (temp > 0)
|
||||
{
|
||||
return x1 + log1p(exp(-temp));
|
||||
}
|
||||
|
||||
return x2 + log1p(exp(temp));
|
||||
}
|
||||
|
||||
}} // Namespace boost::math
|
||||
60
include/boost/math/special_functions/logsumexp.hpp
Normal file
60
include/boost/math/special_functions/logsumexp.hpp
Normal file
@ -0,0 +1,60 @@
|
||||
// (C) Copyright Matt Borland 2022.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include <cmath>
|
||||
#include <iterator>
|
||||
#include <utility>
|
||||
#include <algorithm>
|
||||
#include <type_traits>
|
||||
#include <initializer_list>
|
||||
#include <boost/math/special_functions/logaddexp.hpp>
|
||||
|
||||
namespace boost { namespace math {
|
||||
|
||||
// https://nhigham.com/2021/01/05/what-is-the-log-sum-exp-function/
|
||||
// See equation (#)
|
||||
template <typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>
|
||||
Real logsumexp(ForwardIterator first, ForwardIterator last)
|
||||
{
|
||||
using std::exp;
|
||||
using std::log1p;
|
||||
|
||||
const auto elem = std::max_element(first, last);
|
||||
const Real max_val = *elem;
|
||||
|
||||
Real arg = 0;
|
||||
while (first != last)
|
||||
{
|
||||
if (first != elem)
|
||||
{
|
||||
arg += exp(*first - max_val);
|
||||
}
|
||||
|
||||
++first;
|
||||
}
|
||||
|
||||
return max_val + log1p(arg);
|
||||
}
|
||||
|
||||
template <typename Container, typename Real = typename Container::value_type>
|
||||
inline Real logsumexp(const Container& c)
|
||||
{
|
||||
return logsumexp(std::begin(c), std::end(c));
|
||||
}
|
||||
|
||||
template <typename... Args, typename Real = typename std::common_type<Args...>::type,
|
||||
typename std::enable_if<std::is_floating_point<Real>::value, bool>::type = true>
|
||||
inline Real logsumexp(Args&& ...args)
|
||||
{
|
||||
std::initializer_list<Real> list {std::forward<Args>(args)...};
|
||||
|
||||
if(list.size() == 2)
|
||||
{
|
||||
return logaddexp(*list.begin(), *std::next(list.begin()));
|
||||
}
|
||||
return logsumexp(list.begin(), list.end());
|
||||
}
|
||||
|
||||
}} // Namespace boost::math
|
||||
@ -35,6 +35,26 @@ std::vector<T> generate_random_vector(std::size_t size, std::size_t seed)
|
||||
return v;
|
||||
}
|
||||
|
||||
template<typename T, typename std::enable_if<std::is_floating_point<T>::value, bool>::type = true>
|
||||
std::vector<T> generate_random_vector(std::size_t size, std::size_t seed, T mean, T stddev)
|
||||
{
|
||||
if (seed == 0)
|
||||
{
|
||||
std::random_device rd;
|
||||
seed = rd();
|
||||
}
|
||||
std::vector<T> v(size);
|
||||
|
||||
std::mt19937 gen(seed);
|
||||
|
||||
std::normal_distribution<T> dis(mean, stddev);
|
||||
for (std::size_t i = 0; i < v.size(); ++i)
|
||||
{
|
||||
v[i] = dis(gen);
|
||||
}
|
||||
return v;
|
||||
}
|
||||
|
||||
template<typename T, typename std::enable_if<std::is_integral<T>::value, bool>::type = true>
|
||||
std::vector<T> generate_random_vector(std::size_t size, std::size_t seed)
|
||||
{
|
||||
|
||||
35
reporting/performance/logaddexp_performance.cpp
Normal file
35
reporting/performance/logaddexp_performance.cpp
Normal file
@ -0,0 +1,35 @@
|
||||
// (C) Copyright Matt Borland 2022.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include <random>
|
||||
#include <limits>
|
||||
#include <benchmark/benchmark.h>
|
||||
#include <boost/math/special_functions/logaddexp.hpp>
|
||||
|
||||
using boost::math::logaddexp;
|
||||
|
||||
template <typename Real>
|
||||
void logaddexp_performance(benchmark::State& state)
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937_64 mt {rd()};
|
||||
std::uniform_real_distribution<long double> dist(1e-50, 5e-50);
|
||||
|
||||
Real x = static_cast<Real>(dist(mt));
|
||||
Real y = static_cast<Real>(dist(mt));
|
||||
|
||||
for (auto _ : state)
|
||||
{
|
||||
benchmark::DoNotOptimize(logaddexp(x, y));
|
||||
x += std::numeric_limits<Real>::epsilon();
|
||||
y += std::numeric_limits<Real>::epsilon();
|
||||
}
|
||||
}
|
||||
|
||||
BENCHMARK_TEMPLATE(logaddexp_performance, float);
|
||||
BENCHMARK_TEMPLATE(logaddexp_performance, double);
|
||||
BENCHMARK_TEMPLATE(logaddexp_performance, long double);
|
||||
|
||||
BENCHMARK_MAIN();
|
||||
32
reporting/performance/logsumexp_performance.cpp
Normal file
32
reporting/performance/logsumexp_performance.cpp
Normal file
@ -0,0 +1,32 @@
|
||||
// (C) Copyright Matt Borland 2022.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include <vector>
|
||||
#include <benchmark/benchmark.h>
|
||||
#include <boost/math/special_functions/logsumexp.hpp>
|
||||
#include <boost/math/tools/random_vector.hpp>
|
||||
|
||||
using boost::math::logsumexp;
|
||||
using boost::math::generate_random_vector;
|
||||
|
||||
template <typename Real>
|
||||
void logsumexp_performance(benchmark::State& state)
|
||||
{
|
||||
constexpr std::size_t seed {};
|
||||
const std::size_t size = state.range(0);
|
||||
std::vector<Real> test_set = generate_random_vector<Real>(size, seed);
|
||||
|
||||
for(auto _ : state)
|
||||
{
|
||||
benchmark::DoNotOptimize(logsumexp(test_set));
|
||||
}
|
||||
state.SetComplexityN(state.range(0));
|
||||
}
|
||||
|
||||
BENCHMARK_TEMPLATE(logsumexp_performance, float)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity()->UseRealTime();
|
||||
BENCHMARK_TEMPLATE(logsumexp_performance, double)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity()->UseRealTime();
|
||||
BENCHMARK_TEMPLATE(logsumexp_performance, long double)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity()->UseRealTime();
|
||||
|
||||
BENCHMARK_MAIN();
|
||||
@ -121,6 +121,8 @@ test-suite special_fun :
|
||||
|
||||
[ run hypot_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ]
|
||||
[ run pow_test.cpp ../../test/build//boost_unit_test_framework ]
|
||||
[ run logaddexp_test.cpp ../../test/build//boost_unit_test_framework ]
|
||||
[ run logsumexp_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_hdr_initializer_list cxx11_variadic_templates ] ]
|
||||
[ run ccmath_sqrt_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ]
|
||||
[ run ccmath_isinf_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ]
|
||||
[ run ccmath_isnan_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ]
|
||||
|
||||
23
test/compile_test/sf_logaddexp_incl_test.cpp
Normal file
23
test/compile_test/sf_logaddexp_incl_test.cpp
Normal file
@ -0,0 +1,23 @@
|
||||
// Copyright Matt Borland 2022
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
//
|
||||
// Basic sanity check that header <boost/math/special_functions/log1p.hpp>
|
||||
// #includes all the files that it needs to.
|
||||
//
|
||||
#include <boost/math/special_functions/logaddexp.hpp>
|
||||
//
|
||||
// Note this header includes no other headers, this is
|
||||
// important if this test is to be meaningful:
|
||||
//
|
||||
#include "test_compile_result.hpp"
|
||||
|
||||
void compile_and_link_test()
|
||||
{
|
||||
check_result<float>(boost::math::logaddexp<float>(1.0f, 1.0f));
|
||||
check_result<double>(boost::math::logaddexp<double>(1.0, 1.0));
|
||||
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
|
||||
check_result<long double>(boost::math::logaddexp<long double>(1.0l, 1.0l));
|
||||
#endif
|
||||
}
|
||||
23
test/compile_test/sf_logsumexp_incl_test.cpp
Normal file
23
test/compile_test/sf_logsumexp_incl_test.cpp
Normal file
@ -0,0 +1,23 @@
|
||||
// Copyright Matt Borland 2022
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
//
|
||||
// Basic sanity check that header <boost/math/special_functions/log1p.hpp>
|
||||
// #includes all the files that it needs to.
|
||||
//
|
||||
#include <boost/math/special_functions/logsumexp.hpp>
|
||||
//
|
||||
// Note this header includes no other headers, this is
|
||||
// important if this test is to be meaningful:
|
||||
//
|
||||
#include "test_compile_result.hpp"
|
||||
|
||||
void compile_and_link_test()
|
||||
{
|
||||
check_result<float>(boost::math::logsumexp<float>(1.0f, 1.0f, 1.0f));
|
||||
check_result<double>(boost::math::logsumexp<double>(1.0, 1.0, 1.0));
|
||||
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
|
||||
check_result<long double>(boost::math::logsumexp<long double>(1.0l, 1.0l, 1.0l));
|
||||
#endif
|
||||
}
|
||||
54
test/logaddexp_test.cpp
Normal file
54
test/logaddexp_test.cpp
Normal file
@ -0,0 +1,54 @@
|
||||
// (C) Copyright Matt Borland 2022.
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include "math_unit_test.hpp"
|
||||
#include <limits>
|
||||
#include <boost/math/special_functions/logaddexp.hpp>
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
template <typename Real>
|
||||
void test()
|
||||
{
|
||||
using boost::math::logaddexp;
|
||||
using std::log;
|
||||
using std::exp;
|
||||
|
||||
constexpr Real nan_val = std::numeric_limits<Real>::quiet_NaN();
|
||||
constexpr Real inf_val = std::numeric_limits<Real>::infinity();
|
||||
|
||||
// NAN
|
||||
CHECK_NAN(logaddexp(nan_val, Real(1)));
|
||||
CHECK_NAN(logaddexp(Real(1), nan_val));
|
||||
CHECK_NAN(logaddexp(nan_val, nan_val));
|
||||
|
||||
// INF
|
||||
CHECK_EQUAL(logaddexp(inf_val, Real(1)), inf_val);
|
||||
CHECK_EQUAL(logaddexp(Real(1), inf_val), inf_val);
|
||||
CHECK_EQUAL(logaddexp(inf_val, inf_val), inf_val);
|
||||
|
||||
// Equal values
|
||||
constexpr Real ln2 = boost::math::constants::ln_two<Real>();
|
||||
CHECK_ULP_CLOSE(Real(2) + ln2, logaddexp(Real(2), Real(2)), 1);
|
||||
CHECK_ULP_CLOSE(Real(1e-50) + ln2, logaddexp(Real(1e-50), Real(1e-50)), 1);
|
||||
|
||||
// Spot check
|
||||
// https://numpy.org/doc/stable/reference/generated/numpy.logaddexp.html
|
||||
// Calculated at higher precision using wolfram alpha
|
||||
Real x1 = 1e-50l;
|
||||
Real x2 = 2.5e-50l;
|
||||
Real spot1 = static_cast<Real>(exp(x1));
|
||||
Real spot2 = static_cast<Real>(exp(x2));
|
||||
Real spot12 = logaddexp(x1, x2);
|
||||
|
||||
CHECK_ULP_CLOSE(log(spot1 + spot2), spot12, 1);
|
||||
}
|
||||
|
||||
int main (void)
|
||||
{
|
||||
test<float>();
|
||||
test<double>();
|
||||
test<long double>();
|
||||
return boost::math::test::report_errors();
|
||||
}
|
||||
113
test/logsumexp_test.cpp
Normal file
113
test/logsumexp_test.cpp
Normal file
@ -0,0 +1,113 @@
|
||||
// (C) Copyright Matt Borland and Nick Thompson 2022.
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#include "math_unit_test.hpp"
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
#include <boost/math/special_functions/logsumexp.hpp>
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
#include <boost/math/tools/random_vector.hpp>
|
||||
|
||||
template <typename Real>
|
||||
void test()
|
||||
{
|
||||
using boost::math::logsumexp;
|
||||
using std::log;
|
||||
using std::exp;
|
||||
|
||||
// Spot check 2 values
|
||||
// Also validate that 2 values does not attempt to instantiate the iterator version
|
||||
// https://numpy.org/doc/stable/reference/generated/numpy.logaddexp.html
|
||||
// Calculated at higher precision using wolfram alpha
|
||||
Real x1 = 1e-50l;
|
||||
Real x2 = 2.5e-50l;
|
||||
Real spot1 = static_cast<Real>(exp(x1));
|
||||
Real spot2 = static_cast<Real>(exp(x2));
|
||||
Real spot12 = logsumexp(x1, x2);
|
||||
CHECK_ULP_CLOSE(log(spot1 + spot2), spot12, 1);
|
||||
|
||||
// Spot check 3 values and compare result of each different interface
|
||||
Real x3 = 5e-50l;
|
||||
Real spot3 = static_cast<Real>(exp(x3));
|
||||
std::vector<Real> x_vals {x1, x2, x3};
|
||||
|
||||
Real spot123 = logsumexp(x1, x2, x3);
|
||||
Real spot123_container = logsumexp(x_vals);
|
||||
Real spot123_iter = logsumexp(x_vals.begin(), x_vals.end());
|
||||
|
||||
CHECK_EQUAL(spot123, spot123_container);
|
||||
CHECK_EQUAL(spot123_container, spot123_iter);
|
||||
CHECK_ULP_CLOSE(log(spot1 + spot2 + spot3), spot123, 1);
|
||||
|
||||
// Spot check 4 values with repeated largest value
|
||||
Real x4 = x3;
|
||||
Real spot4 = spot3;
|
||||
Real spot1234 = logsumexp(x1, x2, x3, x4);
|
||||
x_vals.emplace_back(x4);
|
||||
Real spot1234_container = logsumexp(x_vals);
|
||||
|
||||
CHECK_EQUAL(spot1234, spot1234_container);
|
||||
CHECK_ULP_CLOSE(log(spot1 + spot2 + spot3 + spot4), spot1234, 1);
|
||||
|
||||
// Check with a value of vastly different order of magnitude
|
||||
Real x5 = 1.0l;
|
||||
Real spot5 = static_cast<Real>(exp(x5));
|
||||
x_vals.emplace_back(x5);
|
||||
Real spot12345 = logsumexp(x_vals);
|
||||
CHECK_ULP_CLOSE(log(spot1 + spot2 + spot3 + spot4 + spot5), spot12345, 1);
|
||||
}
|
||||
|
||||
// The naive method of computation should overflow:
|
||||
template<typename Real>
|
||||
void test_overflow()
|
||||
{
|
||||
using boost::math::logsumexp;
|
||||
using std::exp;
|
||||
using std::log;
|
||||
|
||||
Real x = ((std::numeric_limits<Real>::max)()/2);
|
||||
|
||||
Real naive_result = log(exp(x) + exp(x));
|
||||
CHECK_EQUAL(std::isfinite(naive_result), false);
|
||||
|
||||
Real result = logsumexp(x, x);
|
||||
CHECK_EQUAL(std::isfinite(result), true);
|
||||
CHECK_ULP_CLOSE(result, x + boost::math::constants::ln_two<Real>(), 1);
|
||||
}
|
||||
|
||||
template <typename Real>
|
||||
void test_random()
|
||||
{
|
||||
using std::exp;
|
||||
using std::log;
|
||||
using boost::math::logsumexp;
|
||||
using boost::math::generate_random_vector;
|
||||
|
||||
std::vector<Real> test_values = generate_random_vector(128, 0, Real(1e-50l), Real(1e-40l));
|
||||
Real naive_exp_sum = 0;
|
||||
|
||||
for(const auto& val : test_values)
|
||||
{
|
||||
naive_exp_sum += exp(val);
|
||||
}
|
||||
|
||||
CHECK_ULP_CLOSE(log(naive_exp_sum), logsumexp(test_values), 1);
|
||||
}
|
||||
|
||||
int main (void)
|
||||
{
|
||||
test<float>();
|
||||
test<double>();
|
||||
test<long double>();
|
||||
|
||||
test_overflow<float>();
|
||||
test_overflow<double>();
|
||||
test_overflow<long double>();
|
||||
|
||||
test_random<float>();
|
||||
test_random<double>();
|
||||
test_random<long double>();
|
||||
return boost::math::test::report_errors();
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user