diff options
Diffstat (limited to 'boost/math/tools')
| -rw-r--r-- | boost/math/tools/bivariate_statistics.hpp | 97 | ||||
| -rw-r--r-- | boost/math/tools/condition_numbers.hpp | 139 | ||||
| -rw-r--r-- | boost/math/tools/config.hpp | 2 | ||||
| -rw-r--r-- | boost/math/tools/detail/is_const_iterable.hpp | 40 | ||||
| -rw-r--r-- | boost/math/tools/norms.hpp | 640 | ||||
| -rw-r--r-- | boost/math/tools/numerical_differentiation.hpp | 261 | ||||
| -rw-r--r-- | boost/math/tools/polynomial.hpp | 73 | ||||
| -rw-r--r-- | boost/math/tools/roots.hpp | 283 | ||||
| -rw-r--r-- | boost/math/tools/signal_statistics.hpp | 346 | ||||
| -rw-r--r-- | boost/math/tools/univariate_statistics.hpp | 393 |
10 files changed, 1991 insertions, 283 deletions
diff --git a/boost/math/tools/bivariate_statistics.hpp b/boost/math/tools/bivariate_statistics.hpp new file mode 100644 index 0000000000..20b7500ed3 --- /dev/null +++ b/boost/math/tools/bivariate_statistics.hpp @@ -0,0 +1,97 @@ +// (C) Copyright Nick Thompson 2018. +// 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) + +#ifndef BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP +#define BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP + +#include <iterator> +#include <tuple> +#include <boost/assert.hpp> +#include <boost/multiprecision/detail/number_base.hpp> + + +namespace boost{ namespace math{ namespace tools { + +template<class Container> +auto means_and_covariance(Container const & u, Container const & v) +{ + using Real = typename Container::value_type; + using std::size; + BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance."); + BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample."); + + // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al. + Real cov = 0; + Real mu_u = u[0]; + Real mu_v = v[0]; + + for(size_t i = 1; i < size(u); ++i) + { + Real u_tmp = (u[i] - mu_u)/(i+1); + Real v_tmp = v[i] - mu_v; + cov += i*u_tmp*v_tmp; + mu_u = mu_u + u_tmp; + mu_v = mu_v + v_tmp/(i+1); + } + + return std::make_tuple(mu_u, mu_v, cov/size(u)); +} + +template<class Container> +auto covariance(Container const & u, Container const & v) +{ + auto [mu_u, mu_v, cov] = boost::math::tools::means_and_covariance(u, v); + return cov; +} + +template<class Container> +auto correlation_coefficient(Container const & u, Container const & v) +{ + using Real = typename Container::value_type; + using std::size; + BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance."); + BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples."); + + Real cov = 0; + Real mu_u = u[0]; + Real mu_v = v[0]; + Real Qu = 0; + Real Qv = 0; + + for(size_t i = 1; i < size(u); ++i) + { + Real u_tmp = u[i] - mu_u; + Real v_tmp = v[i] - mu_v; + Qu = Qu + (i*u_tmp*u_tmp)/(i+1); + Qv = Qv + (i*v_tmp*v_tmp)/(i+1); + cov += i*u_tmp*v_tmp/(i+1); + mu_u = mu_u + u_tmp/(i+1); + mu_v = mu_v + v_tmp/(i+1); + } + + // If both datasets are constant, then they are perfectly correlated. + if (Qu == 0 && Qv == 0) + { + return Real(1); + } + // If one dataset is constant and the other isn't, then they have no correlation: + if (Qu == 0 || Qv == 0) + { + return Real(0); + } + + // Make sure rho in [-1, 1], even in the presence of numerical noise. + Real rho = cov/sqrt(Qu*Qv); + if (rho > 1) { + rho = 1; + } + if (rho < -1) { + rho = -1; + } + return rho; +} + +}}} +#endif diff --git a/boost/math/tools/condition_numbers.hpp b/boost/math/tools/condition_numbers.hpp new file mode 100644 index 0000000000..afbea75cb5 --- /dev/null +++ b/boost/math/tools/condition_numbers.hpp @@ -0,0 +1,139 @@ +// (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) + +#ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP +#define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP +#include <cmath> +#include <boost/math/differentiation/finite_difference.hpp> + +namespace boost::math::tools { + +template<class Real, bool kahan=true> +class summation_condition_number { +public: + summation_condition_number(Real const x = 0) + { + using std::abs; + m_l1 = abs(x); + m_sum = x; + m_c = 0; + } + + void operator+=(Real const & x) + { + using std::abs; + // No need to Kahan the l1 calc; it's well conditioned: + m_l1 += abs(x); + if constexpr(kahan) + { + Real y = x - m_c; + Real t = m_sum + y; + m_c = (t-m_sum) -y; + m_sum = t; + } + else + { + m_sum += x; + } + } + + inline void operator-=(Real const & x) + { + this->operator+=(-x); + } + + // Is operator*= relevant? Presumably everything gets rescaled, + // (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c), + // but is this sensible? More important is it useful? + // In addition, it might change the condition number. + + [[nodiscard]] Real operator()() const + { + using std::abs; + if (m_sum == Real(0) && m_l1 != Real(0)) + { + return std::numeric_limits<Real>::infinity(); + } + return m_l1/abs(m_sum); + } + + [[nodiscard]] Real sum() const + { + // Higham, 1993, "The Accuracy of Floating Point Summation": + // "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected + // thus s=s+e is appended to the algorithm above)." + return m_sum + m_c; + } + + [[nodiscard]] Real l1_norm() const + { + return m_l1; + } + +private: + Real m_l1; + Real m_sum; + Real m_c; +}; + +template<class F, class Real> +auto evaluation_condition_number(F const & f, Real const & x) +{ + using std::abs; + using std::isnan; + using std::sqrt; + using boost::math::differentiation::finite_difference_derivative; + + Real fx = f(x); + if (isnan(fx)) + { + return std::numeric_limits<Real>::quiet_NaN(); + } + bool caught_exception = false; + Real fp; + try + { + fp = finite_difference_derivative(f, x); + } + catch(...) + { + caught_exception = true; + } + + if (isnan(fp) || caught_exception) + { + // Check if the right derivative exists: + fp = finite_difference_derivative<decltype(f), Real, 1>(f, x); + if (isnan(fp)) + { + // Check if a left derivative exists: + const Real eps = (std::numeric_limits<Real>::epsilon)(); + Real h = - 2 * sqrt(eps); + h = boost::math::differentiation::detail::make_xph_representable(x, h); + Real yh = f(x + h); + Real y0 = f(x); + Real diff = yh - y0; + fp = diff / h; + if (isnan(fp)) + { + return std::numeric_limits<Real>::quiet_NaN(); + } + } + } + + if (fx == 0) + { + if (x==0 || fp==0) + { + return std::numeric_limits<Real>::quiet_NaN(); + } + return std::numeric_limits<Real>::infinity(); + } + + return abs(x*fp/fx); +} + +} +#endif diff --git a/boost/math/tools/config.hpp b/boost/math/tools/config.hpp index 17bfec16fe..47fdea65c8 100644 --- a/boost/math/tools/config.hpp +++ b/boost/math/tools/config.hpp @@ -11,7 +11,7 @@ #endif #include <boost/config.hpp> -#include <boost/predef.h> +#include <boost/predef/architecture/x86.h> #include <boost/cstdint.hpp> // for boost::uintmax_t #include <boost/detail/workaround.hpp> #include <boost/type_traits/is_integral.hpp> diff --git a/boost/math/tools/detail/is_const_iterable.hpp b/boost/math/tools/detail/is_const_iterable.hpp new file mode 100644 index 0000000000..805b66e22a --- /dev/null +++ b/boost/math/tools/detail/is_const_iterable.hpp @@ -0,0 +1,40 @@ +// (C) Copyright John Maddock 2018. +// 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) + +#ifndef BOOST_MATH_TOOLS_IS_CONST_ITERABLE_HPP +#define BOOST_MATH_TOOLS_IS_CONST_ITERABLE_HPP + +#if !defined(BOOST_NO_CXX14_VARIABLE_TEMPLATES) && !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_SFINAE_EXPR) + +#define BOOST_MATH_HAS_IS_CONST_ITERABLE + +#include <boost/type_traits/is_detected.hpp> +#include <utility> + +namespace boost { + namespace math { + namespace tools { + namespace detail { + + template<class T> + using begin_t = decltype(std::declval<const T&>().begin()); + template<class T> + using end_t = decltype(std::declval<const T&>().end()); + template<class T> + using const_iterator_t = typename T::const_iterator; + + template <class T> + struct is_const_iterable + : public boost::integral_constant<bool, + boost::is_detected<begin_t, T>::value + && boost::is_detected<end_t, T>::value + && boost::is_detected<const_iterator_t, T>::value + > {}; + +} } } } + +#endif + +#endif // BOOST_MATH_TOOLS_IS_CONST_ITERABLE_HPP diff --git a/boost/math/tools/norms.hpp b/boost/math/tools/norms.hpp new file mode 100644 index 0000000000..478fe04db2 --- /dev/null +++ b/boost/math/tools/norms.hpp @@ -0,0 +1,640 @@ +// (C) Copyright Nick Thompson 2018. +// 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) + +#ifndef BOOST_MATH_TOOLS_NORMS_HPP +#define BOOST_MATH_TOOLS_NORMS_HPP +#include <algorithm> +#include <iterator> +#include <boost/type_traits/is_complex.hpp> +#include <boost/assert.hpp> +#include <boost/multiprecision/detail/number_base.hpp> + + +namespace boost::math::tools { + +// Mallat, "A Wavelet Tour of Signal Processing", equation 2.60: +template<class ForwardIterator> +auto total_variation(ForwardIterator first, ForwardIterator last) +{ + using T = typename std::iterator_traits<ForwardIterator>::value_type; + using std::abs; + BOOST_ASSERT_MSG(first != last && std::next(first) != last, "At least two samples are required to compute the total variation."); + auto it = first; + if constexpr (std::is_unsigned<T>::value) + { + T tmp = *it; + double tv = 0; + while (++it != last) + { + if (*it > tmp) + { + tv += *it - tmp; + } + else + { + tv += tmp - *it; + } + tmp = *it; + } + return tv; + } + else if constexpr (std::is_integral<T>::value) + { + double tv = 0; + double tmp = *it; + while(++it != last) + { + double tmp2 = *it; + tv += abs(tmp2 - tmp); + tmp = *it; + } + return tv; + } + else + { + T tmp = *it; + T tv = 0; + while (++it != last) + { + tv += abs(*it - tmp); + tmp = *it; + } + return tv; + } +} + +template<class Container> +inline auto total_variation(Container const & v) +{ + return total_variation(v.cbegin(), v.cend()); +} + + +template<class ForwardIterator> +auto sup_norm(ForwardIterator first, ForwardIterator last) +{ + BOOST_ASSERT_MSG(first != last, "At least one value is required to compute the sup norm."); + using T = typename std::iterator_traits<ForwardIterator>::value_type; + using std::abs; + if constexpr (boost::is_complex<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + auto it = std::max_element(first, last, [](T a, T b) { return abs(b) > abs(a); }); + return abs(*it); + } + else if constexpr (std::is_unsigned<T>::value) + { + return *std::max_element(first, last); + } + else + { + auto pair = std::minmax_element(first, last); + if (abs(*pair.first) > abs(*pair.second)) + { + return abs(*pair.first); + } + else + { + return abs(*pair.second); + } + } +} + +template<class Container> +inline auto sup_norm(Container const & v) +{ + return sup_norm(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +auto l1_norm(ForwardIterator first, ForwardIterator last) +{ + using T = typename std::iterator_traits<ForwardIterator>::value_type; + using std::abs; + if constexpr (std::is_unsigned<T>::value) + { + double l1 = 0; + for (auto it = first; it != last; ++it) + { + l1 += *it; + } + return l1; + } + else if constexpr (std::is_integral<T>::value) + { + double l1 = 0; + for (auto it = first; it != last; ++it) + { + double tmp = *it; + l1 += abs(tmp); + } + return l1; + } + else + { + decltype(abs(*first)) l1 = 0; + for (auto it = first; it != last; ++it) + { + l1 += abs(*it); + } + return l1; + } + +} + +template<class Container> +inline auto l1_norm(Container const & v) +{ + return l1_norm(v.cbegin(), v.cend()); +} + + +template<class ForwardIterator> +auto l2_norm(ForwardIterator first, ForwardIterator last) +{ + using T = typename std::iterator_traits<ForwardIterator>::value_type; + using std::abs; + using std::norm; + using std::sqrt; + using std::is_floating_point; + if constexpr (boost::is_complex<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + typedef typename T::value_type Real; + Real l2 = 0; + for (auto it = first; it != last; ++it) + { + l2 += norm(*it); + } + Real result = sqrt(l2); + if (!isfinite(result)) + { + Real a = sup_norm(first, last); + l2 = 0; + for (auto it = first; it != last; ++it) + { + l2 += norm(*it/a); + } + return a*sqrt(l2); + } + return result; + } + else if constexpr (is_floating_point<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point) + { + T l2 = 0; + for (auto it = first; it != last; ++it) + { + l2 += (*it)*(*it); + } + T result = sqrt(l2); + // Higham, Accuracy and Stability of Numerical Algorithms, + // Problem 27.5 presents a different algorithm to deal with overflow. + // The algorithm used here takes 3 passes *if* there is overflow. + // Higham's algorithm is 1 pass, but more requires operations than the no oveflow case. + // I'm operating under the assumption that overflow is rare since the dynamic range of floating point numbers is huge. + if (!isfinite(result)) + { + T a = sup_norm(first, last); + l2 = 0; + for (auto it = first; it != last; ++it) + { + T tmp = *it/a; + l2 += tmp*tmp; + } + return a*sqrt(l2); + } + return result; + } + else + { + double l2 = 0; + for (auto it = first; it != last; ++it) + { + double tmp = *it; + l2 += tmp*tmp; + } + return sqrt(l2); + } +} + +template<class Container> +inline auto l2_norm(Container const & v) +{ + return l2_norm(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +size_t l0_pseudo_norm(ForwardIterator first, ForwardIterator last) +{ + using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type; + size_t count = 0; + for (auto it = first; it != last; ++it) + { + if (*it != RealOrComplex(0)) + { + ++count; + } + } + return count; +} + +template<class Container> +inline size_t l0_pseudo_norm(Container const & v) +{ + return l0_pseudo_norm(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +size_t hamming_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2) +{ + size_t count = 0; + auto it1 = first1; + auto it2 = first2; + while (it1 != last1) + { + if (*it1++ != *it2++) + { + ++count; + } + } + return count; +} + +template<class Container> +inline size_t hamming_distance(Container const & v, Container const & w) +{ + return hamming_distance(v.cbegin(), v.cend(), w.cbegin()); +} + +template<class ForwardIterator> +auto lp_norm(ForwardIterator first, ForwardIterator last, unsigned p) +{ + using std::abs; + using std::pow; + using std::is_floating_point; + using std::isfinite; + using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type; + if constexpr (boost::is_complex<RealOrComplex>::value || + boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex) + { + using std::norm; + using Real = typename RealOrComplex::value_type; + Real lp = 0; + for (auto it = first; it != last; ++it) + { + lp += pow(abs(*it), p); + } + + auto result = pow(lp, Real(1)/Real(p)); + if (!isfinite(result)) + { + auto a = boost::math::tools::sup_norm(first, last); + Real lp = 0; + for (auto it = first; it != last; ++it) + { + lp += pow(abs(*it)/a, p); + } + result = a*pow(lp, Real(1)/Real(p)); + } + return result; + } + else if constexpr (is_floating_point<RealOrComplex>::value || + boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point) + { + BOOST_ASSERT_MSG(p >= 0, "For p < 0, the lp norm is not a norm"); + RealOrComplex lp = 0; + + for (auto it = first; it != last; ++it) + { + lp += pow(abs(*it), p); + } + + RealOrComplex result = pow(lp, RealOrComplex(1)/RealOrComplex(p)); + if (!isfinite(result)) + { + RealOrComplex a = boost::math::tools::sup_norm(first, last); + lp = 0; + for (auto it = first; it != last; ++it) + { + lp += pow(abs(*it)/a, p); + } + result = a*pow(lp, RealOrComplex(1)/RealOrComplex(p)); + } + return result; + } + else + { + double lp = 0; + + for (auto it = first; it != last; ++it) + { + double tmp = *it; + lp += pow(abs(tmp), p); + } + double result = pow(lp, 1.0/double(p)); + if (!isfinite(result)) + { + double a = boost::math::tools::sup_norm(first, last); + lp = 0; + for (auto it = first; it != last; ++it) + { + double tmp = *it; + lp += pow(abs(tmp)/a, p); + } + result = a*pow(lp, double(1)/double(p)); + } + return result; + } +} + +template<class Container> +inline auto lp_norm(Container const & v, unsigned p) +{ + return lp_norm(v.cbegin(), v.cend(), p); +} + + +template<class ForwardIterator> +auto lp_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2, unsigned p) +{ + using std::pow; + using std::abs; + using std::is_floating_point; + using std::isfinite; + using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type; + auto it1 = first1; + auto it2 = first2; + + if constexpr (boost::is_complex<RealOrComplex>::value || + boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex) + { + using Real = typename RealOrComplex::value_type; + using std::norm; + Real dist = 0; + while(it1 != last1) + { + auto tmp = *it1++ - *it2++; + dist += pow(abs(tmp), p); + } + return pow(dist, Real(1)/Real(p)); + } + else if constexpr (is_floating_point<RealOrComplex>::value || + boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point) + { + RealOrComplex dist = 0; + while(it1 != last1) + { + auto tmp = *it1++ - *it2++; + dist += pow(abs(tmp), p); + } + return pow(dist, RealOrComplex(1)/RealOrComplex(p)); + } + else + { + double dist = 0; + while(it1 != last1) + { + double tmp1 = *it1++; + double tmp2 = *it2++; + // Naively you'd expect the integer subtraction to be faster, + // but this can overflow or wraparound: + //double tmp = *it1++ - *it2++; + dist += pow(abs(tmp1 - tmp2), p); + } + return pow(dist, 1.0/double(p)); + } +} + +template<class Container> +inline auto lp_distance(Container const & v, Container const & w, unsigned p) +{ + return lp_distance(v.cbegin(), v.cend(), w.cbegin(), p); +} + + +template<class ForwardIterator> +auto l1_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2) +{ + using std::abs; + using std::is_floating_point; + using std::isfinite; + using T = typename std::iterator_traits<ForwardIterator>::value_type; + auto it1 = first1; + auto it2 = first2; + if constexpr (boost::is_complex<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + using Real = typename T::value_type; + Real sum = 0; + while (it1 != last1) { + sum += abs(*it1++ - *it2++); + } + return sum; + } + else if constexpr (is_floating_point<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point) + { + T sum = 0; + while (it1 != last1) + { + sum += abs(*it1++ - *it2++); + } + return sum; + } + else if constexpr (std::is_unsigned<T>::value) + { + double sum = 0; + while(it1 != last1) + { + T x1 = *it1++; + T x2 = *it2++; + if (x1 > x2) + { + sum += (x1 - x2); + } + else + { + sum += (x2 - x1); + } + } + return sum; + } + else if constexpr (std::is_integral<T>::value) + { + double sum = 0; + while(it1 != last1) + { + double x1 = *it1++; + double x2 = *it2++; + sum += abs(x1-x2); + } + return sum; + } + else + { + BOOST_ASSERT_MSG(false, "Could not recognize type."); + } + +} + +template<class Container> +auto l1_distance(Container const & v, Container const & w) +{ + using std::size; + BOOST_ASSERT_MSG(size(v) == size(w), + "L1 distance requires both containers to have the same number of elements"); + return l1_distance(v.cbegin(), v.cend(), w.begin()); +} + +template<class ForwardIterator> +auto l2_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2) +{ + using std::abs; + using std::norm; + using std::sqrt; + using std::is_floating_point; + using std::isfinite; + using T = typename std::iterator_traits<ForwardIterator>::value_type; + auto it1 = first1; + auto it2 = first2; + if constexpr (boost::is_complex<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + using Real = typename T::value_type; + Real sum = 0; + while (it1 != last1) { + sum += norm(*it1++ - *it2++); + } + return sqrt(sum); + } + else if constexpr (is_floating_point<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point) + { + T sum = 0; + while (it1 != last1) + { + T tmp = *it1++ - *it2++; + sum += tmp*tmp; + } + return sqrt(sum); + } + else if constexpr (std::is_unsigned<T>::value) + { + double sum = 0; + while(it1 != last1) + { + T x1 = *it1++; + T x2 = *it2++; + if (x1 > x2) + { + double tmp = x1-x2; + sum += tmp*tmp; + } + else + { + double tmp = x2 - x1; + sum += tmp*tmp; + } + } + return sqrt(sum); + } + else + { + double sum = 0; + while(it1 != last1) + { + double x1 = *it1++; + double x2 = *it2++; + double tmp = x1-x2; + sum += tmp*tmp; + } + return sqrt(sum); + } +} + +template<class Container> +auto l2_distance(Container const & v, Container const & w) +{ + using std::size; + BOOST_ASSERT_MSG(size(v) == size(w), + "L2 distance requires both containers to have the same number of elements"); + return l2_distance(v.cbegin(), v.cend(), w.begin()); +} + +template<class ForwardIterator> +auto sup_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2) +{ + using std::abs; + using std::norm; + using std::sqrt; + using std::is_floating_point; + using std::isfinite; + using T = typename std::iterator_traits<ForwardIterator>::value_type; + auto it1 = first1; + auto it2 = first2; + if constexpr (boost::is_complex<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + using Real = typename T::value_type; + Real sup_sq = 0; + while (it1 != last1) { + Real tmp = norm(*it1++ - *it2++); + if (tmp > sup_sq) { + sup_sq = tmp; + } + } + return sqrt(sup_sq); + } + else if constexpr (is_floating_point<T>::value || + boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point) + { + T sup = 0; + while (it1 != last1) + { + T tmp = *it1++ - *it2++; + if (sup < abs(tmp)) + { + sup = abs(tmp); + } + } + return sup; + } + else // integral values: + { + double sup = 0; + while(it1 != last1) + { + T x1 = *it1++; + T x2 = *it2++; + double tmp; + if (x1 > x2) + { + tmp = x1-x2; + } + else + { + tmp = x2 - x1; + } + if (sup < tmp) { + sup = tmp; + } + } + return sup; + } +} + +template<class Container> +auto sup_distance(Container const & v, Container const & w) +{ + using std::size; + BOOST_ASSERT_MSG(size(v) == size(w), + "sup distance requires both containers to have the same number of elements"); + return sup_distance(v.cbegin(), v.cend(), w.begin()); +} + + +} +#endif diff --git a/boost/math/tools/numerical_differentiation.hpp b/boost/math/tools/numerical_differentiation.hpp index 34fef0db87..52d2788895 100644 --- a/boost/math/tools/numerical_differentiation.hpp +++ b/boost/math/tools/numerical_differentiation.hpp @@ -2,266 +2,11 @@ // 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) - #ifndef BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP #define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP +#include <boost/math/differentiation/finite_difference.hpp> +#include <boost/config/header_deprecated.hpp> -/* - * Performs numerical differentiation by finite-differences. - * - * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception. - * A simple argument demonstrates that the error is unbounded as h->0. - * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h. - * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving - * |f'(x) - (f(x+h) -f(x))/h| < h|f''(x)|/2 + (|f(x)|+|f(x+h)|)eps/h =: g(h), where eps is the unit roundoff for the type. - * It is reasonable to choose h in a way that minimizes the maximum error bound g(h). - * The value of h that minimizes g is h = sqrt(2eps(|f(x)| + |f(x+h)|)/|f''(x)|), and for this value of h the error bound is - * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)). - * In fact it is not necessary to compute the ratio (|f(x+h)| + |f(x)|)/|f''(x)|; the error bound of ~\sqrt{\epsilon} still holds if we set it to one. - * - * - * For more details on this method of analysis, see - * - * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf - * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf - * - * - * It can be shown on general grounds that when choosing the optimal h, the maximum error in f'(x) is ~(|f(x)|eps)^k/k+1|f^(k-1)(x)|^1/k+1. - * From this we can see that full precision can be recovered in the limit k->infinity. - * - * References: - * - * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706. - * - * - * The second algorithm, the complex step derivative, is not ill-conditioned. - * However, it requires that your function can be evaluated at complex arguments. - * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h. - * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that. - * - * References: - * - * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112. - */ - -#include <complex> -#include <boost/math/special_functions/next.hpp> - -// Make sure everyone is informed that C++17 is required: -namespace boost{ namespace math{ namespace tools { - -namespace detail { - template<class Real> - Real make_xph_representable(Real x, Real h) - { - using std::numeric_limits; - // Redefine h so that x + h is representable. Not using this trick leads to large error. - // The compiler flag -ffast-math evaporates these operations . . . - Real temp = x + h; - h = temp - x; - // Handle the case x + h == x: - if (h == 0) - { - h = boost::math::nextafter(x, (numeric_limits<Real>::max)()) - x; - } - return h; - } -} - -template<class F, class Real> -Real complex_step_derivative(const F f, Real x) -{ - // Is it really this easy? Yes. - // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min(). - // This idea was tested over a few billion test cases and found the make the error *much* worse. - // Even 2eps and eps/2 made the error worse, which was surprising. - using std::complex; - using std::numeric_limits; - constexpr const Real step = (numeric_limits<Real>::epsilon)(); - constexpr const Real inv_step = 1/(numeric_limits<Real>::epsilon)(); - return f(complex<Real>(x, step)).imag()*inv_step; -} - -namespace detail { - - template <unsigned> - struct fd_tag {}; - - template<class F, class Real> - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits<Real>::epsilon)(); - // Error bound ~eps^1/2 - // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|). - // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1. - // This approximation will get better as we move to higher orders of accuracy. - Real h = 2 * sqrt(eps); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real y0 = f(x); - Real diff = yh - y0; - if (error) - { - Real ym = f(x - h); - Real ypph = abs(yh - 2 * y0 + ym) / h; - // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h - *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h; - } - return diff / h; - } - - template<class F, class Real> - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits<Real>::epsilon)(); - // Error bound ~eps^2/3 - // See the previous discussion to understand determination of h and the error bound. - // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}] - Real h = pow(3 * eps, static_cast<Real>(1) / static_cast<Real>(3)); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real diff = yh - ymh; - if (error) - { - Real yth = f(x + 2 * h); - Real ymth = f(x - 2 * h); - *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h); - } - - return diff / (2 * h); - } - - template<class F, class Real> - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits<Real>::epsilon)(); - // Error bound ~eps^4/5 - Real h = pow(11.25*eps, (Real)1 / (Real)5); - h = detail::make_xph_representable(x, h); - Real ymth = f(x - 2 * h); - Real yth = f(x + 2 * h); - Real yh = f(x + h); - Real ymh = f(x - h); - Real y2 = ymth - yth; - Real y1 = yh - ymh; - if (error) - { - // Mathematica code to extract the remainder: - // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}] - Real y_three_h = f(x + 3 * h); - Real y_m_three_h = f(x - 3 * h); - // Error from fifth derivative: - *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h); - // Error from function evaluation: - *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h); - } - return (y2 + 8 * y1) / (12 * h); - } - - template<class F, class Real> - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits<Real>::epsilon)(); - // Error bound ~eps^6/7 - // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h - Real h = pow(eps / 168, (Real)1 / (Real)7); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real y1 = yh - ymh; - Real y2 = f(x - 2 * h) - f(x + 2 * h); - Real y3 = f(x + 3 * h) - f(x - 3 * h); - - if (error) - { - // Mathematica code to generate fd scheme for 7th derivative: - // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}] - // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: - // Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}] - Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2; - *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h; - } - return (y3 + 9 * y2 + 45 * y1) / (60 * h); - } - - template<class F, class Real> - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits<Real>::epsilon)(); - // Error bound ~eps^8/9. - // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations. - // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions. - // Mathematica code to get the error: - // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}] - // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h. - Real h = pow(551.25*eps, (Real)1 / (Real)9); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real y1 = yh - ymh; - Real y2 = f(x - 2 * h) - f(x + 2 * h); - Real y3 = f(x + 3 * h) - f(x - 3 * h); - Real y4 = f(x - 4 * h) - f(x + 4 * h); - - Real tmp1 = 3 * y4 / 8 + 4 * y3; - Real tmp2 = 21 * y2 + 84 * y1; - - if (error) - { - // Mathematica code to generate fd scheme for 7th derivative: - // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}] - // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: - // Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h] - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}] - Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1; - *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h; - } - return (tmp1 + tmp2) / (105 * h); - } - - template<class F, class Real, class tag> - Real finite_difference_derivative(const F f, Real x, Real* error, const tag&) - { - // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated. - BOOST_STATIC_ASSERT_MSG(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8"); - } - -} - -template<class F, class Real, size_t order=6> -inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr) -{ - return detail::finite_difference_derivative(f, x, error, detail::fd_tag<order>()); -} +BOOST_HEADER_DEPRECATED("<boost/math/differentiation/finite_difference.hpp>"); -}}} // namespaces #endif diff --git a/boost/math/tools/polynomial.hpp b/boost/math/tools/polynomial.hpp index f6c7ef11ee..357f1313fd 100644 --- a/boost/math/tools/polynomial.hpp +++ b/boost/math/tools/polynomial.hpp @@ -23,6 +23,8 @@ #include <boost/math/policies/error_handling.hpp> #include <boost/math/special_functions/binomial.hpp> #include <boost/core/enable_if.hpp> +#include <boost/type_traits/is_convertible.hpp> +#include <boost/math/tools/detail/is_const_iterable.hpp> #include <vector> #include <ostream> @@ -197,7 +199,7 @@ division(polynomial<T> u, const polynomial<T>& v) BOOST_ASSERT(u); typedef typename polynomial<T>::size_type N; - + N const m = u.size() - 1, n = v.size() - 1; N k = m - n; polynomial<T> q; @@ -301,8 +303,15 @@ public: normalize(); } +#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES + polynomial(std::vector<T>&& p) : m_data(std::move(p)) + { + normalize(); + } +#endif + template <class U> - explicit polynomial(const U& point) + explicit polynomial(const U& point, typename boost::enable_if<boost::is_convertible<U, T> >::type* = 0) { if (point != U(0)) m_data.push_back(point); @@ -321,17 +330,24 @@ public: template <class U> polynomial(const polynomial<U>& p) { + m_data.resize(p.size()); for(unsigned i = 0; i < p.size(); ++i) { - m_data.push_back(boost::math::tools::real_cast<T>(p[i])); + m_data[i] = boost::math::tools::real_cast<T>(p[i]); } } - +#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE + template <class Range> + explicit polynomial(const Range& r, typename boost::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = 0) + : polynomial(r.begin(), r.end()) + { + } +#endif #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500) polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l)) { } - + polynomial& operator=(std::initializer_list<T> l) { @@ -358,9 +374,15 @@ public: { return m_data[i]; } + T evaluate(T z) const { - return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : 0; + return this->operator()(z); + } + + T operator()(T z) const + { + return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : T(0); } std::vector<T> chebyshev() const { @@ -377,8 +399,34 @@ public: return m_data; } - // operators: #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES + polynomial<T> prime() const + { + if (m_data.size() == 0) + { + return polynomial<T>({}); + } + + std::vector<T> p_data(m_data.size() - 1); + for (size_t i = 0; i < p_data.size(); ++i) { + p_data[i] = m_data[i+1]*static_cast<T>(i+1); + } + return polynomial<T>(std::move(p_data)); + } + + polynomial<T> integrate() const + { + std::vector<T> i_data(m_data.size() + 1); + // Choose integration constant such that P(0) = 0. + i_data[0] = T(0); + for (size_t i = 1; i < i_data.size(); ++i) + { + i_data[i] = m_data[i-1]/static_cast<T>(i); + } + return polynomial<T>(std::move(i_data)); + } + + // operators: polynomial& operator =(polynomial&& p) BOOST_NOEXCEPT { m_data = std::move(p.m_data); @@ -497,13 +545,13 @@ public: normalize(); return *this; } - + // Convenient and efficient query for zero. bool is_zero() const { return m_data.empty(); } - + // Conversion to bool. #ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS typedef bool (polynomial::*unmentionable_type)() const; @@ -524,13 +572,13 @@ public: { m_data.clear(); } - + /** Remove zero coefficients 'from the top', that is for which there are no * non-zero coefficients of higher degree. */ void normalize() { #ifndef BOOST_NO_CXX11_LAMBDAS - m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != 0; }).base(), m_data.end()); + m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end()); #else using namespace boost::lambda; m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end()); @@ -834,6 +882,3 @@ inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, #include <boost/math/tools/polynomial_gcd.hpp> #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP - - - diff --git a/boost/math/tools/roots.hpp b/boost/math/tools/roots.hpp index 81e8466138..e16294dc15 100644 --- a/boost/math/tools/roots.hpp +++ b/boost/math/tools/roots.hpp @@ -9,7 +9,10 @@ #ifdef _MSC_VER #pragma once #endif +#include <boost/multiprecision/detail/number_base.hpp> // test for multiprecision types. +#include <boost/type_traits/is_complex.hpp> // test for complex types +#include <iostream> #include <utility> #include <boost/config/no_tr1/cmath.hpp> #include <stdexcept> @@ -65,7 +68,7 @@ inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T) { using dummy::get; // Rely on ADL to find the correct overload of get: - val = get<0>(t); + val = get<0>(t); } template <class T, class U, class V> @@ -268,11 +271,15 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ if(fabs(delta * 2) > fabs(delta2)) { // last two steps haven't converged. - delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; - if (fabs(delta) > fabs(result)) + T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2; + if ((result != 0) && (fabs(shift) > fabs(result))) + { delta = sign(delta) * result; // protect against huge jumps! - // reset delta2 so we don't take this branch next time round: - delta1 = delta; + } + else + delta = shift; + // reset delta1/2 so we don't take this branch next time round: + delta1 = 3 * delta; delta2 = 3 * delta; } guess = result; @@ -421,12 +428,12 @@ namespace detail{ { // last two steps haven't converged. delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; - if (fabs(delta) > result) + if ((result != 0) && (fabs(delta) > result)) delta = sign(delta) * result; // protect against huge jumps! // reset delta2 so that this branch will *not* be taken on the // next iteration: delta2 = delta * 3; - delta1 = delta; + delta1 = delta * 3; BOOST_MATH_INSTRUMENT_VARIABLE(delta); } guess = result; @@ -436,7 +443,10 @@ namespace detail{ // check for out of bounds step: if(result < min) { - T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min); + T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) + ? T(1000) + : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result)) + ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min); if(fabs(diff) < 1) diff = 1 / diff; if(!out_of_bounds_sentry && (diff > 0) && (diff < 3)) @@ -514,9 +524,10 @@ namespace detail{ template <class T> static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T)) { + using std::fabs; T ratio = f0 / f1; T delta; - if(ratio / x < 0.1) + if((x != 0) && (fabs(ratio / x) < 0.1)) { delta = ratio + (f2 / (2 * f1)) * ratio * ratio; // check second derivative doesn't over compensate: @@ -559,10 +570,262 @@ inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEP return schroder_iterate(f, guess, min, max, digits, m); } +#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS +/* + * Why do we set the default maximum number of iterations to the number of digits in the type? + * Because for double roots, the number of digits increases linearly with the number of iterations, + * so this default should recover full precision even in this somewhat pathological case. + * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all. + */ +template<class Complex, class F> +Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits<typename Complex::value_type>::digits) +{ + typedef typename Complex::value_type Real; + using std::norm; + using std::abs; + using std::max; + // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method: + Complex z0 = guess + Complex(1,0); + Complex z1 = guess + Complex(0,1); + Complex z2 = guess; + + do { + auto pair = g(z2); + if (norm(pair.second) == 0) + { + // Muller's method. Notation follows Numerical Recipes, 9.5.2: + Complex q = (z2 - z1)/(z1 - z0); + auto P0 = g(z0); + auto P1 = g(z1); + Complex qp1 = static_cast<Complex>(1)+q; + Complex A = q*(pair.first - qp1*P1.first + q*P0.first); + + Complex B = (static_cast<Complex>(2)*q+static_cast<Complex>(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first; + Complex C = qp1*pair.first; + Complex rad = sqrt(B*B - static_cast<Complex>(4)*A*C); + Complex denom1 = B + rad; + Complex denom2 = B - rad; + Complex correction = (z1-z2)*static_cast<Complex>(2)*C; + if (norm(denom1) > norm(denom2)) + { + correction /= denom1; + } + else + { + correction /= denom2; + } + + z0 = z1; + z1 = z2; + z2 = z2 + correction; + } + else + { + z0 = z1; + z1 = z2; + z2 = z2 - (pair.first/pair.second); + } + + // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root + // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0. + // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered. + Real tol = max(abs(z2)*std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon()); + bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol; + bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol; + if (real_close && imag_close) + { + return z2; + } + + } while(max_iterations--); + + // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations + // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps + // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code, + // I found this condition generates correct roots, whereas the scale invariant condition discussed here: + // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method + // allows nonroots to be passed off as roots. + auto pair = g(z2); + if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon())) + { + return z2; + } + + return {std::numeric_limits<Real>::quiet_NaN(), + std::numeric_limits<Real>::quiet_NaN()}; +} +#endif + + +#if !defined(BOOST_NO_CXX17_IF_CONSTEXPR) +// https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711 +namespace detail +{ + template<class T> + inline T discriminant(T const & a, T const & b, T const & c) + { + T w = 4*a*c; + T e = std::fma(-c, 4*a, w); + T f = std::fma(b, b, -w); + return f + e; + } +} + +template<class T> +auto quadratic_roots(T const& a, T const& b, T const& c) +{ + using std::copysign; + using std::sqrt; + if constexpr (std::is_integral<T>::value) + { + // What I want is to write: + // return quadratic_roots(double(a), double(b), double(c)); + // but that doesn't compile. + double nan = std::numeric_limits<double>::quiet_NaN(); + if(a==0) + { + if (b==0 && c != 0) + { + return std::pair<double, double>(nan, nan); + } + else if (b==0 && c==0) + { + return std::pair<double, double>(0,0); + } + return std::pair<double, double>(-c/b, -c/b); + } + if (b==0) + { + double x0_sq = -double(c)/double(a); + if (x0_sq < 0) { + return std::pair<double, double>(nan, nan); + } + double x0 = sqrt(x0_sq); + return std::pair<double, double>(-x0,x0); + } + double discriminant = detail::discriminant(double(a), double(b), double(c)); + if (discriminant < 0) + { + return std::pair<double, double>(nan, nan); + } + double q = -(b + copysign(sqrt(discriminant), double(b)))/T(2); + double x0 = q/a; + double x1 = c/q; + if (x0 < x1) { + return std::pair<double, double>(x0, x1); + } + return std::pair<double, double>(x1, x0); + } + else if constexpr (std::is_floating_point<T>::value) + { + T nan = std::numeric_limits<T>::quiet_NaN(); + if(a==0) + { + if (b==0 && c != 0) + { + return std::pair<T, T>(nan, nan); + } + else if (b==0 && c==0) + { + return std::pair<T, T>(0,0); + } + return std::pair<T, T>(-c/b, -c/b); + } + if (b==0) + { + T x0_sq = -c/a; + if (x0_sq < 0) { + return std::pair<T, T>(nan, nan); + } + T x0 = sqrt(x0_sq); + return std::pair<T, T>(-x0,x0); + } + T discriminant = detail::discriminant(a, b, c); + // Is there a sane way to flush very small negative values to zero? + // If there is I don't know of it. + if (discriminant < 0) + { + return std::pair<T, T>(nan, nan); + } + T q = -(b + copysign(sqrt(discriminant), b))/T(2); + T x0 = q/a; + T x1 = c/q; + if (x0 < x1) + { + return std::pair<T, T>(x0, x1); + } + return std::pair<T, T>(x1, x0); + } + else if constexpr (boost::is_complex<T>::value || boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex) + { + typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN(); + if(a.real()==0 && a.imag() ==0) + { + using std::norm; + if (b.real()==0 && b.imag() && norm(c) != 0) + { + return std::pair<T, T>({nan, nan}, {nan, nan}); + } + else if (b.real()==0 && b.imag() && c.real() ==0 && c.imag() == 0) + { + return std::pair<T, T>({0,0},{0,0}); + } + return std::pair<T, T>(-c/b, -c/b); + } + if (b.real()==0 && b.imag() == 0) + { + T x0_sq = -c/a; + T x0 = sqrt(x0_sq); + return std::pair<T, T>(-x0, x0); + } + // There's no fma for complex types: + T discriminant = b*b - T(4)*a*c; + T q = -(b + sqrt(discriminant))/T(2); + return std::pair<T, T>(q/a, c/q); + } + else // Most likely the type is a boost.multiprecision. + { //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation. + T nan = std::numeric_limits<T>::quiet_NaN(); + if(a==0) + { + if (b==0 && c != 0) + { + return std::pair<T, T>(nan, nan); + } + else if (b==0 && c==0) + { + return std::pair<T, T>(0,0); + } + return std::pair<T, T>(-c/b, -c/b); + } + if (b==0) + { + T x0_sq = -c/a; + if (x0_sq < 0) { + return std::pair<T, T>(nan, nan); + } + T x0 = sqrt(x0_sq); + return std::pair<T, T>(-x0,x0); + } + T discriminant = b*b - 4*a*c; + if (discriminant < 0) + { + return std::pair<T, T>(nan, nan); + } + T q = -(b + copysign(sqrt(discriminant), b))/T(2); + T x0 = q/a; + T x1 = c/q; + if (x0 < x1) + { + return std::pair<T, T>(x0, x1); + } + return std::pair<T, T>(x1, x0); + } +} +#endif } // namespace tools } // namespace math } // namespace boost #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP - diff --git a/boost/math/tools/signal_statistics.hpp b/boost/math/tools/signal_statistics.hpp new file mode 100644 index 0000000000..74f9dfd031 --- /dev/null +++ b/boost/math/tools/signal_statistics.hpp @@ -0,0 +1,346 @@ +// (C) Copyright Nick Thompson 2018. +// 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) + +#ifndef BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP +#define BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP + +#include <algorithm> +#include <iterator> +#include <boost/type_traits/is_complex.hpp> +#include <boost/assert.hpp> +#include <boost/multiprecision/detail/number_base.hpp> +#include <boost/math/tools/roots.hpp> +#include <boost/math/tools/univariate_statistics.hpp> + + +namespace boost::math::tools { + +template<class ForwardIterator> +auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last) +{ + using std::abs; + using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type; + BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples."); + + std::sort(first, last, [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); }); + + + decltype(abs(*first)) i = 1; + decltype(abs(*first)) num = 0; + decltype(abs(*first)) denom = 0; + for (auto it = first; it != last; ++it) + { + decltype(abs(*first)) tmp = abs(*it); + num += tmp*i; + denom += tmp; + ++i; + } + + // If the l1 norm is zero, all elements are zero, so every element is the same. + if (denom == 0) + { + decltype(abs(*first)) zero = 0; + return zero; + } + return ((2*num)/denom - i)/(i-1); +} + +template<class RandomAccessContainer> +inline auto absolute_gini_coefficient(RandomAccessContainer & v) +{ + return boost::math::tools::absolute_gini_coefficient(v.begin(), v.end()); +} + +template<class ForwardIterator> +auto sample_absolute_gini_coefficient(ForwardIterator first, ForwardIterator last) +{ + size_t n = std::distance(first, last); + return n*boost::math::tools::absolute_gini_coefficient(first, last)/(n-1); +} + +template<class RandomAccessContainer> +inline auto sample_absolute_gini_coefficient(RandomAccessContainer & v) +{ + return boost::math::tools::sample_absolute_gini_coefficient(v.begin(), v.end()); +} + + +// The Hoyer sparsity measure is defined in: +// https://arxiv.org/pdf/0811.4706.pdf +template<class ForwardIterator> +auto hoyer_sparsity(const ForwardIterator first, const ForwardIterator last) +{ + using T = typename std::iterator_traits<ForwardIterator>::value_type; + using std::abs; + using std::sqrt; + BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Hoyer sparsity requires at least two samples."); + + if constexpr (std::is_unsigned<T>::value) + { + T l1 = 0; + T l2 = 0; + size_t n = 0; + for (auto it = first; it != last; ++it) + { + l1 += *it; + l2 += (*it)*(*it); + n += 1; + } + + double rootn = sqrt(n); + return (rootn - l1/sqrt(l2) )/ (rootn - 1); + } + else { + decltype(abs(*first)) l1 = 0; + decltype(abs(*first)) l2 = 0; + // We wouldn't need to count the elements if it was a random access iterator, + // but our only constraint is that it's a forward iterator. + size_t n = 0; + for (auto it = first; it != last; ++it) + { + decltype(abs(*first)) tmp = abs(*it); + l1 += tmp; + l2 += tmp*tmp; + n += 1; + } + if constexpr (std::is_integral<T>::value) + { + double rootn = sqrt(n); + return (rootn - l1/sqrt(l2) )/ (rootn - 1); + } + else + { + decltype(abs(*first)) rootn = sqrt(static_cast<decltype(abs(*first))>(n)); + return (rootn - l1/sqrt(l2) )/ (rootn - 1); + } + } +} + +template<class Container> +inline auto hoyer_sparsity(Container const & v) +{ + return boost::math::tools::hoyer_sparsity(v.cbegin(), v.cend()); +} + + +template<class Container> +auto oracle_snr(Container const & signal, Container const & noisy_signal) +{ + using Real = typename Container::value_type; + BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(), + "Signal and noisy_signal must be have the same number of elements."); + if constexpr (std::is_integral<Real>::value) + { + double numerator = 0; + double denominator = 0; + for (size_t i = 0; i < signal.size(); ++i) + { + numerator += signal[i]*signal[i]; + denominator += (noisy_signal[i] - signal[i])*(noisy_signal[i] - signal[i]); + } + if (numerator == 0 && denominator == 0) + { + return std::numeric_limits<double>::quiet_NaN(); + } + if (denominator == 0) + { + return std::numeric_limits<double>::infinity(); + } + return numerator/denominator; + } + else if constexpr (boost::is_complex<Real>::value || + boost::multiprecision::number_category<Real>::value == boost::multiprecision::number_kind_complex) + + { + using std::norm; + typename Real::value_type numerator = 0; + typename Real::value_type denominator = 0; + for (size_t i = 0; i < signal.size(); ++i) + { + numerator += norm(signal[i]); + denominator += norm(noisy_signal[i] - signal[i]); + } + if (numerator == 0 && denominator == 0) + { + return std::numeric_limits<typename Real::value_type>::quiet_NaN(); + } + if (denominator == 0) + { + return std::numeric_limits<typename Real::value_type>::infinity(); + } + + return numerator/denominator; + } + else + { + Real numerator = 0; + Real denominator = 0; + for (size_t i = 0; i < signal.size(); ++i) + { + numerator += signal[i]*signal[i]; + denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]); + } + if (numerator == 0 && denominator == 0) + { + return std::numeric_limits<Real>::quiet_NaN(); + } + if (denominator == 0) + { + return std::numeric_limits<Real>::infinity(); + } + + return numerator/denominator; + } +} + +template<class Container> +auto mean_invariant_oracle_snr(Container const & signal, Container const & noisy_signal) +{ + using Real = typename Container::value_type; + BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(), "Signal and noisy signal must be have the same number of elements."); + + Real mu = boost::math::tools::mean(signal); + Real numerator = 0; + Real denominator = 0; + for (size_t i = 0; i < signal.size(); ++i) + { + Real tmp = signal[i] - mu; + numerator += tmp*tmp; + denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]); + } + if (numerator == 0 && denominator == 0) + { + return std::numeric_limits<Real>::quiet_NaN(); + } + if (denominator == 0) + { + return std::numeric_limits<Real>::infinity(); + } + + return numerator/denominator; + +} + +template<class Container> +auto mean_invariant_oracle_snr_db(Container const & signal, Container const & noisy_signal) +{ + using std::log10; + return 10*log10(boost::math::tools::mean_invariant_oracle_snr(signal, noisy_signal)); +} + + +// Follows the definition of SNR given in Mallat, A Wavelet Tour of Signal Processing, equation 11.16. +template<class Container> +auto oracle_snr_db(Container const & signal, Container const & noisy_signal) +{ + using std::log10; + return 10*log10(boost::math::tools::oracle_snr(signal, noisy_signal)); +} + +// A good reference on the M2M4 estimator: +// D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR estimation techniques for the AWGN channel," IEEE Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000. +// A nice python implementation: +// https://github.com/gnuradio/gnuradio/blob/master/gr-digital/examples/snr_estimators.py +template<class ForwardIterator> +auto m2m4_snr_estimator(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3) +{ + BOOST_ASSERT_MSG(estimated_signal_kurtosis > 0, "The estimated signal kurtosis must be positive"); + BOOST_ASSERT_MSG(estimated_noise_kurtosis > 0, "The estimated noise kurtosis must be positive."); + using Real = typename std::iterator_traits<ForwardIterator>::value_type; + using std::sqrt; + if constexpr (std::is_floating_point<Real>::value || + boost::multiprecision::number_category<Real>::value == boost::multiprecision::number_kind_floating_point) + { + // If we first eliminate N, we obtain the quadratic equation: + // (ka+kw-6)S^2 + 2M2(3-kw)S + kw*M2^2 - M4 = 0 =: a*S^2 + bs*N + cs = 0 + // If we first eliminate S, we obtain the quadratic equation: + // (ka+kw-6)N^2 + 2M2(3-ka)N + ka*M2^2 - M4 = 0 =: a*N^2 + bn*N + cn = 0 + // I believe these equations are totally independent quadratics; + // if one has a complex solution it is not necessarily the case that the other must also. + // However, I can't prove that, so there is a chance that this does unnecessary work. + // Future improvements: There are algorithms which can solve quadratics much more effectively than the naive implementation found here. + // See: https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711#50065711 + auto [M1, M2, M3, M4] = boost::math::tools::first_four_moments(first, last); + if (M4 == 0) + { + // The signal is constant. There is no noise: + return std::numeric_limits<Real>::infinity(); + } + // Change to notation in Pauluzzi, equation 41: + auto kw = estimated_noise_kurtosis; + auto ka = estimated_signal_kurtosis; + // A common case, since it's the default: + Real a = (ka+kw-6); + Real bs = 2*M2*(3-kw); + Real cs = kw*M2*M2 - M4; + Real bn = 2*M2*(3-ka); + Real cn = ka*M2*M2 - M4; + auto [S0, S1] = boost::math::tools::quadratic_roots(a, bs, cs); + if (S1 > 0) + { + auto N = M2 - S1; + if (N > 0) + { + return S1/N; + } + if (S0 > 0) + { + N = M2 - S0; + if (N > 0) + { + return S0/N; + } + } + } + auto [N0, N1] = boost::math::tools::quadratic_roots(a, bn, cn); + if (N1 > 0) + { + auto S = M2 - N1; + if (S > 0) + { + return S/N1; + } + if (N0 > 0) + { + S = M2 - N0; + if (S > 0) + { + return S/N0; + } + } + } + // This happens distressingly often. It's a limitation of the method. + return std::numeric_limits<Real>::quiet_NaN(); + } + else + { + BOOST_ASSERT_MSG(false, "The M2M4 estimator has not been implemented for this type."); + return std::numeric_limits<Real>::quiet_NaN(); + } +} + +template<class Container> +inline auto m2m4_snr_estimator(Container const & noisy_signal, typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3) +{ + return m2m4_snr_estimator(noisy_signal.cbegin(), noisy_signal.cend(), estimated_signal_kurtosis, estimated_noise_kurtosis); +} + +template<class ForwardIterator> +inline auto m2m4_snr_estimator_db(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3) +{ + using std::log10; + return 10*log10(m2m4_snr_estimator(first, last, estimated_signal_kurtosis, estimated_noise_kurtosis)); +} + + +template<class Container> +inline auto m2m4_snr_estimator_db(Container const & noisy_signal, typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3) +{ + using std::log10; + return 10*log10(m2m4_snr_estimator(noisy_signal, estimated_signal_kurtosis, estimated_noise_kurtosis)); +} + +} +#endif diff --git a/boost/math/tools/univariate_statistics.hpp b/boost/math/tools/univariate_statistics.hpp new file mode 100644 index 0000000000..226fdf46d2 --- /dev/null +++ b/boost/math/tools/univariate_statistics.hpp @@ -0,0 +1,393 @@ +// (C) Copyright Nick Thompson 2018. +// 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) + +#ifndef BOOST_MATH_TOOLS_UNIVARIATE_STATISTICS_HPP +#define BOOST_MATH_TOOLS_UNIVARIATE_STATISTICS_HPP + +#include <algorithm> +#include <iterator> +#include <boost/type_traits/is_complex.hpp> +#include <boost/assert.hpp> +#include <boost/multiprecision/detail/number_base.hpp> + + +namespace boost::math::tools { + +template<class ForwardIterator> +auto mean(ForwardIterator first, ForwardIterator last) +{ + using Real = typename std::iterator_traits<ForwardIterator>::value_type; + BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean."); + if constexpr (std::is_integral<Real>::value) + { + double mu = 0; + double i = 1; + for(auto it = first; it != last; ++it) { + mu = mu + (*it - mu)/i; + i += 1; + } + return mu; + } + else + { + Real mu = 0; + Real i = 1; + for(auto it = first; it != last; ++it) { + mu = mu + (*it - mu)/i; + i += 1; + } + return mu; + } +} + +template<class Container> +inline auto mean(Container const & v) +{ + return mean(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +auto variance(ForwardIterator first, ForwardIterator last) +{ + using Real = typename std::iterator_traits<ForwardIterator>::value_type; + BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance."); + // Higham, Accuracy and Stability, equation 1.6a and 1.6b: + if constexpr (std::is_integral<Real>::value) + { + double M = *first; + double Q = 0; + double k = 2; + for (auto it = std::next(first); it != last; ++it) + { + double tmp = *it - M; + Q = Q + ((k-1)*tmp*tmp)/k; + M = M + tmp/k; + k += 1; + } + return Q/(k-1); + } + else + { + Real M = *first; + Real Q = 0; + Real k = 2; + for (auto it = std::next(first); it != last; ++it) + { + Real tmp = (*it - M)/k; + Q += k*(k-1)*tmp*tmp; + M += tmp; + k += 1; + } + return Q/(k-1); + } +} + +template<class Container> +inline auto variance(Container const & v) +{ + return variance(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +auto sample_variance(ForwardIterator first, ForwardIterator last) +{ + size_t n = std::distance(first, last); + BOOST_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance."); + return n*variance(first, last)/(n-1); +} + +template<class Container> +inline auto sample_variance(Container const & v) +{ + return sample_variance(v.cbegin(), v.cend()); +} + + +// Follows equation 1.5 of: +// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf +template<class ForwardIterator> +auto skewness(ForwardIterator first, ForwardIterator last) +{ + using Real = typename std::iterator_traits<ForwardIterator>::value_type; + BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute skewness."); + if constexpr (std::is_integral<Real>::value) + { + double M1 = *first; + double M2 = 0; + double M3 = 0; + double n = 2; + for (auto it = std::next(first); it != last; ++it) + { + double delta21 = *it - M1; + double tmp = delta21/n; + M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2); + M2 = M2 + tmp*(n-1)*delta21; + M1 = M1 + tmp; + n += 1; + } + + double var = M2/(n-1); + if (var == 0) + { + // The limit is technically undefined, but the interpretation here is clear: + // A constant dataset has no skewness. + return double(0); + } + double skew = M3/(M2*sqrt(var)); + return skew; + } + else + { + Real M1 = *first; + Real M2 = 0; + Real M3 = 0; + Real n = 2; + for (auto it = std::next(first); it != last; ++it) + { + Real delta21 = *it - M1; + Real tmp = delta21/n; + M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2); + M2 += tmp*(n-1)*delta21; + M1 += tmp; + n += 1; + } + + Real var = M2/(n-1); + if (var == 0) + { + // The limit is technically undefined, but the interpretation here is clear: + // A constant dataset has no skewness. + return Real(0); + } + Real skew = M3/(M2*sqrt(var)); + return skew; + } +} + +template<class Container> +inline auto skewness(Container const & v) +{ + return skewness(v.cbegin(), v.cend()); +} + +// Follows equation 1.5/1.6 of: +// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf +template<class ForwardIterator> +auto first_four_moments(ForwardIterator first, ForwardIterator last) +{ + using Real = typename std::iterator_traits<ForwardIterator>::value_type; + BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the first four moments."); + if constexpr (std::is_integral<Real>::value) + { + double M1 = *first; + double M2 = 0; + double M3 = 0; + double M4 = 0; + double n = 2; + for (auto it = std::next(first); it != last; ++it) + { + double delta21 = *it - M1; + double tmp = delta21/n; + M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3); + M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2); + M2 = M2 + tmp*(n-1)*delta21; + M1 = M1 + tmp; + n += 1; + } + + return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1)); + } + else + { + Real M1 = *first; + Real M2 = 0; + Real M3 = 0; + Real M4 = 0; + Real n = 2; + for (auto it = std::next(first); it != last; ++it) + { + Real delta21 = *it - M1; + Real tmp = delta21/n; + M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3); + M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2); + M2 = M2 + tmp*(n-1)*delta21; + M1 = M1 + tmp; + n += 1; + } + + return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1)); + } +} + +template<class Container> +inline auto first_four_moments(Container const & v) +{ + return first_four_moments(v.cbegin(), v.cend()); +} + + +// Follows equation 1.6 of: +// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf +template<class ForwardIterator> +auto kurtosis(ForwardIterator first, ForwardIterator last) +{ + auto [M1, M2, M3, M4] = first_four_moments(first, last); + if (M2 == 0) + { + return M2; + } + return M4/(M2*M2); +} + +template<class Container> +inline auto kurtosis(Container const & v) +{ + return kurtosis(v.cbegin(), v.cend()); +} + +template<class ForwardIterator> +auto excess_kurtosis(ForwardIterator first, ForwardIterator last) +{ + return kurtosis(first, last) - 3; +} + +template<class Container> +inline auto excess_kurtosis(Container const & v) +{ + return excess_kurtosis(v.cbegin(), v.cend()); +} + + +template<class RandomAccessIterator> +auto median(RandomAccessIterator first, RandomAccessIterator last) +{ + size_t num_elems = std::distance(first, last); + BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined."); + if (num_elems & 1) + { + auto middle = first + (num_elems - 1)/2; + std::nth_element(first, middle, last); + return *middle; + } + else + { + auto middle = first + num_elems/2 - 1; + std::nth_element(first, middle, last); + std::nth_element(middle, middle+1, last); + return (*middle + *(middle+1))/2; + } +} + + +template<class RandomAccessContainer> +inline auto median(RandomAccessContainer & v) +{ + return median(v.begin(), v.end()); +} + +template<class RandomAccessIterator> +auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) +{ + using Real = typename std::iterator_traits<RandomAccessIterator>::value_type; + BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples."); + + std::sort(first, last); + if constexpr (std::is_integral<Real>::value) + { + double i = 1; + double num = 0; + double denom = 0; + for (auto it = first; it != last; ++it) + { + num += *it*i; + denom += *it; + ++i; + } + + // If the l1 norm is zero, all elements are zero, so every element is the same. + if (denom == 0) + { + return double(0); + } + + return ((2*num)/denom - i)/(i-1); + } + else + { + Real i = 1; + Real num = 0; + Real denom = 0; + for (auto it = first; it != last; ++it) + { + num += *it*i; + denom += *it; + ++i; + } + + // If the l1 norm is zero, all elements are zero, so every element is the same. + if (denom == 0) + { + return Real(0); + } + + return ((2*num)/denom - i)/(i-1); + } +} + +template<class RandomAccessContainer> +inline auto gini_coefficient(RandomAccessContainer & v) +{ + return gini_coefficient(v.begin(), v.end()); +} + +template<class RandomAccessIterator> +inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) +{ + size_t n = std::distance(first, last); + return n*gini_coefficient(first, last)/(n-1); +} + +template<class RandomAccessContainer> +inline auto sample_gini_coefficient(RandomAccessContainer & v) +{ + return sample_gini_coefficient(v.begin(), v.end()); +} + +template<class RandomAccessIterator> +auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN()) +{ + using std::abs; + using Real = typename std::iterator_traits<RandomAccessIterator>::value_type; + using std::isnan; + if (isnan(center)) + { + center = boost::math::tools::median(first, last); + } + size_t num_elems = std::distance(first, last); + BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined."); + auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);}; + if (num_elems & 1) + { + auto middle = first + (num_elems - 1)/2; + std::nth_element(first, middle, last, comparator); + return abs(*middle); + } + else + { + auto middle = first + num_elems/2 - 1; + std::nth_element(first, middle, last, comparator); + std::nth_element(middle, middle+1, last, comparator); + return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2)); + } +} + +template<class RandomAccessContainer> +inline auto median_absolute_deviation(RandomAccessContainer & v, typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN()) +{ + return median_absolute_deviation(v.begin(), v.end(), center); +} + +} +#endif |