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

Refactor so as to not store a reference member, make call threadsafe, compute entire vector in one go.

Browse Source
This commit is contained in:
Nick Thompson
2019-01-04 12:50:58 -07:00
parent 95f993c9bc
commit bfa7619954
3 changed files with 205 additions and 106 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
doc/differentiation/lanczos_smoothing.qbk
View File

@@ -14,18 +14,18 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
namespace boost::math::differentiation {
template <class RandomAccessContainer, size_t order=1>
template <class Real, size_t order=1>
class discrete_lanczos_derivative {
public:
using Real = typename RandomAccessContainer::value_type;
discrete_lanczos_derivative(RandomAccessContainer const & v,
Real spacing = 1,
discrete_lanczos_derivative(Real spacing,
size_t n = 18,
size_t approximation_order = 3);
Real operator[](size_t i) const;
template<class RandomAccessContainer>
Real operator()(RandomAccessContainer const & v, size_t i) const;
void reset_data(RandomAccessContainer const &v);
template<class RandomAccessContainer>
RandomAccessContainer operator()(RandomAccessContainer const & v) const;
void reset_spacing(Real spacing);
};
@@ -35,24 +35,29 @@ namespace boost::math::differentiation {
[heading Description]
The `discrete_lanczos_derivative` class calculates a finite-difference approximation to the derivative of a noisy sequence of equally-spaced values /v/ at an index /i/.
The `discrete_lanczos_derivative` class calculates a finite-difference approximation to the derivative of a noisy sequence of equally-spaced values /v/.
A basic usage is
std::vector<double> v(500);
// fill v with noisy data.
double spacing = 0.001;
using boost::math::differentiation::discrete_lanczos_derivative;
auto lanczos = discrete_lanczos_derivative(v, spacing);
double dvdt = lanczos[30];
auto lanczos = discrete_lanczos_derivative(spacing);
// Compute dvdt at index 30:
double dvdt30 = lanczos(v, 30);
// Compute derivative of entire vector:
std::vector<double> dvdt = lanczos(v);
Noise-suppressing second derivatives can also be computed.
The syntax is as follows:
std::vector<double> v(500);
// fill v with noisy data.
auto lanczos = lanczos_derivative<decltype(v), 2>(v, spacing);
// evaluate:
double dvdt = lanczos[25];
auto lanczos = lanczos_derivative<double, 2>(spacing);
// evaluate second derivative at a point:
double d2vdt2 = lanczos(v, 18);
// evaluate second derivative of entire vector:
std::vector<double> d2vdt2 = lanczos(v);
If the data has variance \u03C3[super 2],
then the variance of the computed derivative is roughly \u03C3[super 2]/p/[super 3] /n/[super -3] \u0394 /t/[super -2],
@@ -63,8 +68,8 @@ You can play around with the approximation order /p/ and the filter length /n/:
size_t n = 12;
size_t p = 2;
auto lanczos = lanczos_derivative(v, spacing, n, p);
double dvdt = lanczos[24];
auto lanczos = lanczos_derivative(spacing, n, p);
double dvdt = lanczos(v, i);
If /p=2n/, then the discrete Lanczos derivative is not smoothing:
It reduces to the standard /2n+1/-point finite-difference formula.
@@ -79,7 +84,7 @@ balanced against the danger of overfitting and averaging over non-stationarity.
The choice of approximation order /p/ for a given /n/ is more difficult.
If your signal is believed to be a polynomial,
it does not make sense to set /p/ to larger than the polynomial degree-
though it may be sensible to take /p/ less than the polynomial degree.
though it may be sensible to take /p/ less than this.
For a sinusoidal signal contaminated with AWGN, we ran a few tests showing that for SNR = 1,
p = n/8 gave the best results,
@@ -93,18 +98,17 @@ Since each filter has length /2n+1/ and there are /n/ filters, whose element eac
the complexity of the constructor call is O(/n/[super 2]/p/).
This is not cheap-though for most cases small /p/ and /n/ not too large (< 20) is desired.
However, for concreteness, on the author's 2.7GHz Intel laptop CPU, the /n=16/, /p=3/ filter takes 9 microseconds to compute.
This is far from negligible, and as such we provide API calls which allow the filters to be used with multiple data:
This is far from negligible, and as such the filters may be used with multiple data, and even shared between threads:
std::vector<double> v(500);
// fill v with noisy data.
auto lanczos = lanczos_derivative(v, spacing);
// use lanczos with v . . .
std::vector<double> w(500);
lanczos.reset_data(w);
// use lanczos with w . . .
std::vector<double> v1(500);
std::vector<double> v2(500);
// fill v1, v2 with noisy data.
auto lanczos = lanczos_derivative(spacing);
std::vector<double> dv1dt = lanczos(v1); // threadsafe
std::vector<double> dv2dt = lanczos(v2); // threadsafe
// need to use a different spacing?
lanczos.reset_spacing(0.02);
lanczos.reset_spacing(0.02); // not threadsafe
The implementation follows [@https://doi.org/10.1080/00207160.2012.666348 McDevitt, 2012],

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

@@ -258,20 +258,16 @@ std::vector<Real> acceleration_boundary_filter(size_t n, size_t p, int64_t s)
} // namespace detail
template <class RandomAccessContainer, size_t order = 1>
template <typename Real, size_t order = 1>
class discrete_lanczos_derivative {
public:
using Real = typename RandomAccessContainer::value_type;
discrete_lanczos_derivative(RandomAccessContainer const &v,
Real const & spacing = 1,
size_t n = 18,
size_t approximation_order = 3)
: m_v{v}, dt{spacing}
discrete_lanczos_derivative(Real const & spacing,
size_t n = 18,
size_t approximation_order = 3)
: m_dt{spacing}
{
static_assert(!std::is_integral_v<Real>, "Spacing must be a floating point type.");
BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0.");
using std::size;
BOOST_ASSERT_MSG(size(v) >= 2*n+1,
"Vector must be at least as long as the filter length");
if constexpr (order == 1)
{
@@ -281,12 +277,12 @@ public:
"The approximation order must be >= 2");
m_f = detail::interior_filter<Real>(n, approximation_order);
boundary_filters.resize(n);
m_boundary_filters.resize(n);
for (size_t i = 0; i < n; ++i)
{
// s = i - n;
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
boundary_filters[i] = detail::boundary_filter<Real>(n, approximation_order, s);
m_boundary_filters[i] = detail::boundary_filter<Real>(n, approximation_order, s);
}
}
else if constexpr (order == 2)
@@ -297,11 +293,11 @@ public:
{
m_f[i] = f[i+n];
}
boundary_filters.resize(n);
m_boundary_filters.resize(n);
for (size_t i = 0; i < n; ++i)
{
int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
boundary_filters[i] = detail::acceleration_boundary_filter<Real>(n, approximation_order, s);
m_boundary_filters[i] = detail::acceleration_boundary_filter<Real>(n, approximation_order, s);
}
}
else
@@ -310,96 +306,96 @@ public:
}
}
void reset_data(RandomAccessContainer const &v)
{
using std::size;
BOOST_ASSERT_MSG(size(v) >= m_f.size(), "Vector must be at least as long as the filter length");
m_v = v;
}
void reset_spacing(Real const & spacing)
{
BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0.");
dt = spacing;
m_dt = spacing;
}
Real spacing() const
{
return dt;
return m_dt;
}
Real operator[](size_t i) const
template<class RandomAccessContainer>
Real operator()(RandomAccessContainer const & v, size_t i) const
{
static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,
"The type of the values in the vector provided does not match the type in the filters.");
using std::size;
BOOST_ASSERT_MSG(size(v) >= m_boundary_filters[0].size(),
"Vector must be at least as long as the filter length");
if constexpr (order==1)
{
using std::size;
if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size())
if (i >= m_f.size() - 1 && i <= size(v) - m_f.size())
{
Real dv = 0;
for (size_t j = 1; j < m_f.size(); ++j)
{
dv += m_f[j] * (m_v[i + j] - m_v[i - j]);
dv += m_f[j] * (v[i + j] - v[i - j]);
}
return dv / dt;
return dv / m_dt;
}
// m_f.size() = N+1
if (i < m_f.size() - 1)
{
auto &f = boundary_filters[i];
auto &bf = m_boundary_filters[i];
Real dv = 0;
for (size_t j = 0; j < f.size(); ++j) {
dv += f[j] * m_v[j];
for (size_t j = 0; j < bf.size(); ++j)
{
dv += bf[j] * v[j];
}
return dv / dt;
return dv / m_dt;
}
if (i > size(m_v) - m_f.size() && i < size(m_v))
if (i > size(v) - m_f.size() && i < size(v))
{
int k = size(m_v) - 1 - i;
auto &f = boundary_filters[k];
int k = size(v) - 1 - i;
auto &bf = m_boundary_filters[k];
Real dv = 0;
for (size_t j = 0; j < f.size(); ++j)
for (size_t j = 0; j < bf.size(); ++j)
{
dv += f[j] * m_v[m_v.size() - 1 - j];
dv += bf[j] * v[size(v) - 1 - j];
}
return -dv / dt;
return -dv / m_dt;
}
}
else if constexpr (order==2)
{
using std::size;
if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size())
if (i >= m_f.size() - 1 && i <= size(v) - m_f.size())
{
Real d2v = m_f[0]*m_v[i];
Real d2v = m_f[0]*v[i];
for (size_t j = 1; j < m_f.size(); ++j)
{
d2v += m_f[j] * (m_v[i + j] + m_v[i - j]);
d2v += m_f[j] * (v[i + j] + v[i - j]);
}
return d2v / (dt*dt);
return d2v / (m_dt*m_dt);
}
// m_f.size() = N+1
if (i < m_f.size() - 1)
{
auto &f = boundary_filters[i];
auto &bf = m_boundary_filters[i];
Real d2v = 0;
for (size_t j = 0; j < f.size(); ++j) {
d2v += f[j] * m_v[j];
for (size_t j = 0; j < bf.size(); ++j)
{
d2v += bf[j] * v[j];
}
return d2v / (dt*dt);
return d2v / (m_dt*m_dt);
}
if (i > size(m_v) - m_f.size() && i < size(m_v))
if (i > size(v) - m_f.size() && i < size(v))
{
int k = size(m_v) - 1 - i;
auto &f = boundary_filters[k];
int k = size(v) - 1 - i;
auto &bf = m_boundary_filters[k];
Real d2v = 0;
for (size_t j = 0; j < f.size(); ++j)
for (size_t j = 0; j < bf.size(); ++j)
{
d2v += f[j] * m_v[m_v.size() - 1 - j];
d2v += bf[j] * v[size(v) - 1 - j];
}
return d2v / dt;
return d2v / (m_dt*m_dt);
}
}
@@ -408,11 +404,97 @@ public:
return std::numeric_limits<Real>::quiet_NaN();
}
template<class RandomAccessContainer>
RandomAccessContainer operator()(RandomAccessContainer const & v) const
{
static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,
"The type of the values in the vector provided does not match the type in the filters.");
using std::size;
BOOST_ASSERT_MSG(size(v) >= m_boundary_filters[0].size(),
"Vector must be at least as long as the filter length");
RandomAccessContainer w(size(v));
if constexpr (order==1)
{
for (size_t i = 0; i < m_f.size() - 1; ++i)
{
auto &bf = m_boundary_filters[i];
Real dv = 0;
for (size_t j = 0; j < bf.size(); ++j)
{
dv += bf[j] * v[j];
}
w[i] = dv / m_dt;
}
for(size_t i = m_f.size() - 1; i <= size(v) - m_f.size(); ++i)
{
Real dv = 0;
for (size_t j = 1; j < m_f.size(); ++j)
{
dv += m_f[j] * (v[i + j] - v[i - j]);
}
w[i] = dv / m_dt;
}
for(size_t i = size(v) - m_f.size() + 1; i < size(v); ++i)
{
int k = size(v) - 1 - i;
auto &f = m_boundary_filters[k];
Real dv = 0;
for (size_t j = 0; j < f.size(); ++j)
{
dv += f[j] * v[size(v) - 1 - j];
}
w[i] = -dv / m_dt;
}
}
else if constexpr (order==2)
{
// m_f.size() = N+1
for (size_t i = 0; i < m_f.size() - 1; ++i)
{
auto &bf = m_boundary_filters[i];
Real d2v = 0;
for (size_t j = 0; j < bf.size(); ++j)
{
d2v += bf[j] * v[j];
}
w[i] = d2v / (m_dt*m_dt);
}
for (size_t i = m_f.size() - 1; i <= size(v) - m_f.size(); ++i)
{
Real d2v = m_f[0]*v[i];
for (size_t j = 1; j < m_f.size(); ++j)
{
d2v += m_f[j] * (v[i + j] + v[i - j]);
}
w[i] = d2v / (m_dt*m_dt);
}
for (size_t i = size(v) - m_f.size() + 1; i < size(v); ++i)
{
int k = size(v) - 1 - i;
auto &bf = m_boundary_filters[k];
Real d2v = 0;
for (size_t j = 0; j < bf.size(); ++j)
{
d2v += bf[j] * v[size(v) - 1 - j];
}
w[i] = d2v / (m_dt*m_dt);
}
}
return w;
}
private:
const RandomAccessContainer &m_v;
std::vector<Real> m_f;
std::vector<std::vector<Real>> boundary_filters;
Real dt;
std::vector<std::vector<Real>> m_boundary_filters;
Real m_dt;
};
} // namespaces

65
test/lanczos_smoothing_test.cpp
View File

@@ -208,12 +208,17 @@ void test_interior_lanczos()
{
for (size_t p = 2; p < 2*n; p += 2)
{
auto lsd = discrete_lanczos_derivative(v, 0.1, n, p);
auto dld = discrete_lanczos_derivative(Real(0.1), n, p);
for (size_t m = n; m < v.size() - n; ++m)
{
Real dvdt = lsd[m];
Real dvdt = dld(v, m);
BOOST_CHECK_SMALL(dvdt, tol);
}
auto dvdt = dld(v);
for (size_t m = n; m < v.size() - n; ++m)
{
BOOST_CHECK_SMALL(dvdt[m], tol);
}
}
}
@@ -227,12 +232,17 @@ void test_interior_lanczos()
{
for (size_t p = 2; p < 2*n; p += 2)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto dld = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = n; m < v.size() - n; ++m)
{
Real dvdt = lsd[m];
Real dvdt = dld(v, m);
BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 2000*tol);
}
auto dvdt = dld(v);
for (size_t m = n; m < v.size() - n; ++m)
{
BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 7, 2000*tol);
}
}
}
@@ -241,7 +251,6 @@ void test_interior_lanczos()
//std::cout << "Seed = " << seed << "\n";
std::mt19937 gen(4172378669);
std::normal_distribution<> dis{0, 0.01};
std::cout << std::fixed;
for (size_t i = 0; i < v.size(); ++i)
{
v[i] = 7*i+8 + dis(gen);
@@ -251,10 +260,10 @@ void test_interior_lanczos()
{
for (size_t p = 2; p < 2*n; p += 2)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto dld = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = n; m < v.size() - n; ++m)
{
BOOST_CHECK_CLOSE_FRACTION(lsd[m], Real(7), Real(0.0042));
BOOST_CHECK_CLOSE_FRACTION(dld(v, m), Real(7), Real(0.0042));
}
}
}
@@ -269,10 +278,10 @@ void test_interior_lanczos()
{
for (size_t p = 2; p < 2*n; p += 2)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto dld = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = n; m < v.size() - n; ++m)
{
BOOST_CHECK_CLOSE_FRACTION(lsd[m], Real(30*m + 7), Real(0.00008));
BOOST_CHECK_CLOSE_FRACTION(dld(v,m), Real(30*m + 7), Real(0.00008));
}
}
}
@@ -288,11 +297,11 @@ void test_interior_lanczos()
{
for (size_t p = 3; p < 100 && p < n/2; p += 2)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto dld = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = n; m < v.size() - n && m < n + 10; ++m)
{
BOOST_CHECK_CLOSE_FRACTION(lsd[m], omega*cos(omega*m), Real(0.03));
BOOST_CHECK_CLOSE_FRACTION(dld(v,m), omega*cos(omega*m), Real(0.03));
}
}
}
@@ -377,15 +386,15 @@ void test_boundary_lanczos()
{
for (size_t p = 2; p < 2*n; ++p)
{
auto lsd = discrete_lanczos_derivative(v, 0.0125, n, p);
auto lsd = discrete_lanczos_derivative(Real(0.0125), n, p);
for (size_t m = 0; m < n; ++m)
{
Real dvdt = lsd[m];
Real dvdt = lsd(v,m);
BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol));
}
for (size_t m = v.size() - n; m < v.size(); ++m)
{
Real dvdt = lsd[m];
Real dvdt = lsd(v,m);
BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol));
}
}
@@ -400,16 +409,16 @@ void test_boundary_lanczos()
{
for (size_t p = 2; p < 2*n; ++p)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto lsd = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = 0; m < n; ++m)
{
Real dvdt = lsd[m];
Real dvdt = lsd(v,m);
BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, sqrt(tol));
}
for (size_t m = v.size() - n; m < v.size(); ++m)
{
Real dvdt = lsd[m];
Real dvdt = lsd(v,m);
BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 4*sqrt(tol));
}
}
@@ -424,10 +433,10 @@ void test_boundary_lanczos()
{
for (size_t p = 2; p < 2*n; ++p)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto lsd = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = 0; m < v.size(); ++m)
{
BOOST_CHECK_CLOSE_FRACTION(lsd[m], 30*m+7, 30*sqrt(tol));
BOOST_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 30*sqrt(tol));
}
}
}
@@ -447,14 +456,18 @@ void test_boundary_lanczos()
{
for (size_t p = 2; p < n; ++p)
{
auto lsd = discrete_lanczos_derivative(v, Real(1), n, p);
auto lsd = discrete_lanczos_derivative(Real(1), n, p);
for (size_t m = 0; m < v.size(); ++m)
{
BOOST_CHECK_CLOSE_FRACTION(lsd[m], 30*m+7, 0.005);
BOOST_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 0.005);
}
auto dvdt = lsd(v);
for (size_t m = 0; m < v.size(); ++m)
{
BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 30*m+7, 0.005);
}
}
}
}
template<class Real>
@@ -502,10 +515,10 @@ void test_lanczos_acceleration()
{
Real eps = std::numeric_limits<Real>::epsilon();
std::vector<Real> v(100, 7);
auto lanczos = discrete_lanczos_derivative<decltype(v), 2>(v, Real(1), 4, 3);
auto lanczos = discrete_lanczos_derivative<Real, 2>(Real(1), 4, 3);
for (size_t i = 0; i < v.size(); ++i)
{
BOOST_CHECK_SMALL(lanczos[i], eps);
BOOST_CHECK_SMALL(lanczos(v, i), eps);
}
for(size_t i = 0; i < v.size(); ++i)
@@ -514,7 +527,7 @@ void test_lanczos_acceleration()
}
for (size_t i = 0; i < v.size(); ++i)
{
BOOST_CHECK_SMALL(lanczos[i], 200*eps);
BOOST_CHECK_SMALL(lanczos(v,i), 200*eps);
}
for(size_t i = 0; i < v.size(); ++i)
@@ -523,7 +536,7 @@ void test_lanczos_acceleration()
}
for (size_t i = 0; i < v.size(); ++i)
{
BOOST_CHECK_CLOSE_FRACTION(lanczos[i], 14, 1000*eps);
BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1000*eps);
}