boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Compute filters in higher precision and cast back to low precision so that high-p filters are accurate. [CI SKIP]

Browse Source
This commit is contained in:
Nick Thompson
2019-01-22 15:46:19 -07:00
parent 9aee84a1ec
commit b2f3054e2f
3 changed files with 74 additions and 130 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/differentiation/lanczos_smoothing.qbk
View File

@@ -27,7 +27,9 @@ namespace boost::math::differentiation {
template<class RandomAccessContainer>
RandomAccessContainer operator()(RandomAccessContainer const & v) const;
void reset_spacing(Real spacing);
void set_spacing(Real spacing);
Real get_spacing() const;
};
} // namespaces

56
include/boost/math/differentiation/lanczos_smoothing.hpp
View File

@@ -1,4 +1,4 @@
// (C) Copyright Nick Thompson 2018.
// (C) Copyright Nick Thompson 2019.
// 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)
@@ -7,7 +7,6 @@
#define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP
#include <vector>
#include <boost/assert.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
namespace boost::math::differentiation {
@@ -276,23 +275,53 @@ public:
"The approximation order must be <= 2n");
BOOST_ASSERT_MSG(approximation_order >= 2,
"The approximation order must be >= 2");
m_f = detail::interior_filter<Real>(n, approximation_order);
if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)
{
auto interior = detail::interior_filter<long double>(n, approximation_order);
m_f.resize(interior.size());
for (size_t j = 0; j < interior.size(); ++j)
{
m_f[j] = static_cast<Real>(interior[j]);
}
}
else
{
m_f = detail::interior_filter<Real>(n, approximation_order);
}
m_boundary_filters.resize(n);
// This for loop is a natural candidate for parallelization.
// But does it matter? Probably not.
for (size_t i = 0; i < n; ++i)
{
// s = i - n;
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
m_boundary_filters[i] = detail::boundary_filter<Real>(n, approximation_order, s);
if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)
{
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
auto bf = detail::boundary_filter<long double>(n, approximation_order, s);
m_boundary_filters[i].resize(bf.size());
for (size_t j = 0; j < bf.size(); ++j)
{
m_boundary_filters[i][j] = static_cast<Real>(bf[j]);
}
}
else
{
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
m_boundary_filters[i] = detail::boundary_filter<Real>(n, approximation_order, s);
}
}
}
else if constexpr (order == 2)
{
// High precision isn't warranted for small p; only for large p.
// (The computation appears stable for large n.)
// But given that the filters are reusable for many vectors,
// it's better to do a high precision computation and then cast back,
// since the resulting cost is a factor of 2, and the cost of the filters not working is hours of debugging.
if constexpr (std::is_same_v<Real, double> || std::is_same_v<Real, float>)
{
std::cout << "We're using high precision!\n";
using boost::multiprecision::cpp_bin_float_50;
auto f = detail::acceleration_boundary_filter<cpp_bin_float_50>(n, approximation_order, 0);
auto f = detail::acceleration_boundary_filter<long double>(n, approximation_order, 0);
m_f.resize(n+1);
for (size_t i = 0; i < m_f.size(); ++i)
{
@@ -302,7 +331,7 @@ public:
for (size_t i = 0; i < n; ++i)
{
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
auto bf = detail::acceleration_boundary_filter<cpp_bin_float_50>(n, approximation_order, s);
auto bf = detail::acceleration_boundary_filter<long double>(n, approximation_order, s);
m_boundary_filters[i].resize(bf.size());
for (size_t j = 0; j < bf.size(); ++j)
{
@@ -312,7 +341,8 @@ public:
}
else
{
std::cout << "No high precision.\n";
// Given that the purpose is denoising, for higher precision calculations,
// the default precision should be fine.
auto f = detail::acceleration_boundary_filter<Real>(n, approximation_order, 0);
m_f.resize(n+1);
for (size_t i = 0; i < m_f.size(); ++i)
@@ -333,13 +363,13 @@ public:
}
}
void reset_spacing(Real const & spacing)
void set_spacing(Real const & spacing)
{
BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0.");
m_dt = spacing;
}
Real spacing() const
Real get_spacing() const
{
return m_dt;
}

144
test/lanczos_smoothing_test.cpp
View File

@@ -1,5 +1,5 @@
/*
* Copyright Nick Thompson, 2017
* Copyright Nick Thompson, 2019
* 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)
@@ -27,54 +27,6 @@ using boost::math::differentiation::detail::discrete_legendre;
using boost::math::differentiation::detail::interior_filter;
using boost::math::differentiation::detail::boundary_filter;
template<class RandomAccessContainer1, class RandomAccessContainer2>
size_t l1_ulp_error(RandomAccessContainer1 const & v1, RandomAccessContainer2 const & v2)
{
using std::abs;
using Real1 = typename RandomAccessContainer1::value_type;
using Real2 = typename RandomAccessContainer2::value_type;
BOOST_ASSERT_MSG(sizeof(Real1) <= sizeof(Real2),
"Second container must contain more accurate value than the first container.");
BOOST_ASSERT_MSG(v1.size() == v2.size(),
"Containers must have same size.");
size_t error = 0;
for (size_t i = 0; i < v1.size(); ++i)
{
auto err = abs(boost::math::float_distance<Real1>(v1[i], static_cast<Real1>(v2[i])));
if (abs(v2[i]) > std::numeric_limits<Real1>::epsilon())
{
error += static_cast<size_t>(err);
}
}
return error;
}
template<class RandomAccessContainer1, class RandomAccessContainer2>
std::tuple<size_t, typename RandomAccessContainer1::value_type, typename RandomAccessContainer2::value_type> sup_ulp_error(RandomAccessContainer1 const & v1, RandomAccessContainer2 const & v2)
{
using std::abs;
using Real1 = typename RandomAccessContainer1::value_type;
using Real2 = typename RandomAccessContainer2::value_type;
BOOST_ASSERT_MSG(sizeof(Real1) <= sizeof(Real2),
"Second container must contain more accurate value than the first container.");
BOOST_ASSERT_MSG(v1.size() == v2.size(),
"Containers must have same size.");
size_t error = 0;
Real1 worst_val1 = std::numeric_limits<Real1>::quiet_NaN();
Real2 worst_val2 = std::numeric_limits<Real2>::quiet_NaN();
for (size_t i = 0; i < v1.size(); ++i)
{
auto err = abs(boost::math::float_distance<Real1>(v1[i], static_cast<Real1>(v2[i])));
if (err > error)
{
error = err;
worst_val1 = v1[i];
worst_val2 = v2[i];
}
}
return {error, worst_val1, worst_val2};
}
template<class Real>
void test_dlp_norms()
{
@@ -526,40 +478,18 @@ template<class Real>
void test_acceleration_filters()
{
Real eps = std::numeric_limits<Real>::epsilon();
std::cout << std::setprecision(std::numeric_limits<double>::digits10);
std::cout << std::scientific;
for (size_t n = 1; n < 100; ++n)
{
for(size_t p = 3; p <= 2*n; ++p)
{
for(int64_t s = -int64_t(n); s <= 0 /*int64_t(n)*/; ++s)
{
auto v1 = boost::math::differentiation::detail::acceleration_boundary_filter<double>(n,p,s);
auto v2 = boost::math::differentiation::detail::acceleration_boundary_filter<long double>(n,p,s);
size_t dist1 = l1_ulp_error(v1, v2);
auto [dist2, worst_val1, worst_val2] = sup_ulp_error(v1, v2);
std::cout << "(n,p,s) = (" << n << ", " << p << ", " << s << ") = "
<< dist1 << ", sup = " << dist2 << ", worst = " << worst_val1 << ", actual = " << worst_val2 << "\n";
}
}
}
/*Real eps = std::numeric_limits<Real>::epsilon();
for (size_t n = 1; n < 15; ++n)
for (size_t n = 1; n < 5; ++n)
{
for(size_t p = 3; p <= 2*n; ++p)
{
for(int64_t s = -int64_t(n); s <= 0; ++s)
{
auto f = boost::math::differentiation::detail::acceleration_boundary_filter<Real>(n,p,s);
if (std::is_same_v<Real, double>)
auto g = boost::math::differentiation::detail::acceleration_boundary_filter<long double>(n,p,s);
std::vector<Real> f(g.size());
for (size_t i = 0; i < g.size(); ++i)
{
//std::cout << "We're using high precision!\n";
auto f1 = boost::math::differentiation::detail::acceleration_boundary_filter<cpp_bin_float_100>(n,p,s);
for (size_t i = 0; i < f1.size(); ++i)
{
f[i] = static_cast<Real>(f1[i]);
}
f[i] = static_cast<Real>(g[i]);
}
Real sum = 0;
@@ -590,6 +520,9 @@ void test_acceleration_filters()
BOOST_CHECK_CLOSE_FRACTION(sum, 2, 4*sqrt(eps));
// See unlabelled equation in McDevitt, 2012, just after equation 26:
// It appears that there is an off-by-one error in that equation, since p + 1 moments don't vanish, only p.
// This test is itself suspect; the condition number of the moment sum is infinite.
// So the *slightest* error in the filter gets amplified by the test; in terms of the
// behavior of the actual filter, it's not a big deal.
for (size_t l = 3; l <= p; ++l)
{
sum = 0;
@@ -603,11 +536,11 @@ void test_acceleration_filters()
c = (t - sum) - y;
sum = t;
}
BOOST_CHECK_SMALL(abs(sum), 500*l*sqrt(eps));
BOOST_CHECK_SMALL(abs(sum), 500*sqrt(eps));
}
}
}
}*/
}
}
template<class Real>
@@ -636,56 +569,32 @@ void test_lanczos_acceleration()
}
for (size_t i = 0; i < v.size(); ++i)
{
BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1000*eps);
BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1500*eps);
}
v.resize(2048);
Real step = two_pi<Real>()/v.size();
for(size_t i = 0; i < v.size(); ++i)
{
Real x = i*step;
v[i] = sin(x);
}
std::random_device rd{};
auto seed = rd();
std::cout << "seed = " << seed << "\n";
// Now add noise, and kick up the smoothing of the Lanzcos derivative (increase n):
//std::random_device rd{};
//auto seed = rd();
//std::cout << "seed = " << seed << "\n";
size_t seed = 2507134629;
std::mt19937 gen(seed);
std::normal_distribution<> dis{0, 0.00001};
Real std_dev = 0.1;
std::normal_distribution<Real> dis{0, std_dev};
for (size_t i = 0; i < v.size(); ++i)
{
v[i] += dis(gen);
}
size_t n = 100;
lanczos = discrete_lanczos_derivative<Real, 2>(step, n, 3);
lanczos = discrete_lanczos_derivative<Real, 2>(Real(1), 18, 3);
auto w = lanczos(v);
BOOST_TEST(w.size() == v.size());
BOOST_CHECK_SMALL(w[0], 0.01);
for(size_t i = 1; i < n; ++i)
for (size_t i = 0; i < v.size(); ++i)
{
BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01);
BOOST_CHECK_CLOSE_FRACTION(w[i], 14, std_dev/200);
}
for(size_t i = n; i < v.size() -n; ++i)
{
BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01);
}
/*for(size_t i = v.size() - n; i < v.size(); ++i)
{
BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01);
}*/
}
BOOST_AUTO_TEST_CASE(lanczos_smoothing_test)
{
//test_acceleration_filters<double>();
//test_acceleration_filters<cpp_bin_float_50>();
/*test_dlp_second_derivative<double>();
test_dlp_second_derivative<double>();
test_dlp_norms<double>();
test_dlp_evaluation<double>();
test_dlp_derivatives<double>();
@@ -701,7 +610,10 @@ BOOST_AUTO_TEST_CASE(lanczos_smoothing_test)
test_interior_filter<double>();
test_interior_filter<long double>();
test_interior_lanczos<double>();*/
test_interior_lanczos<double>();
test_acceleration_filters<double>();
test_lanczos_acceleration<float>();
test_lanczos_acceleration<double>();
}