639 lines · plain
1// (C) Copyright Nick Thompson 2018.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_NORMS_HPP7#define BOOST_MATH_TOOLS_NORMS_HPP8#include <algorithm>9#include <iterator>10#include <complex>11#include <cmath>12#include <boost/math/tools/assert.hpp>13#include <boost/math/tools/complex.hpp>14 15#include <boost/math/tools/is_standalone.hpp>16#ifndef BOOST_MATH_STANDALONE17#include <boost/config.hpp>18#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR19#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."20#endif21#endif22 23 24namespace boost::math::tools {25 26// Mallat, "A Wavelet Tour of Signal Processing", equation 2.60:27template<class ForwardIterator>28auto total_variation(ForwardIterator first, ForwardIterator last)29{30 using T = typename std::iterator_traits<ForwardIterator>::value_type;31 using std::abs;32 BOOST_MATH_ASSERT_MSG(first != last && std::next(first) != last, "At least two samples are required to compute the total variation.");33 auto it = first;34 if constexpr (std::is_unsigned<T>::value)35 {36 T tmp = *it;37 double tv = 0;38 while (++it != last)39 {40 if (*it > tmp)41 {42 tv += *it - tmp;43 }44 else45 {46 tv += tmp - *it;47 }48 tmp = *it;49 }50 return tv;51 }52 else if constexpr (std::is_integral<T>::value)53 {54 double tv = 0;55 double tmp = *it;56 while(++it != last)57 {58 double tmp2 = *it;59 tv += abs(tmp2 - tmp);60 tmp = *it;61 }62 return tv;63 }64 else65 {66 T tmp = *it;67 T tv = 0;68 while (++it != last)69 {70 tv += abs(*it - tmp);71 tmp = *it;72 }73 return tv;74 }75}76 77template<class Container>78inline auto total_variation(Container const & v)79{80 return total_variation(v.cbegin(), v.cend());81}82 83 84template<class ForwardIterator>85auto sup_norm(ForwardIterator first, ForwardIterator last)86{87 BOOST_MATH_ASSERT_MSG(first != last, "At least one value is required to compute the sup norm.");88 using T = typename std::iterator_traits<ForwardIterator>::value_type;89 using std::abs;90 if constexpr (boost::math::tools::is_complex_type<T>::value)91 {92 auto it = std::max_element(first, last, [](T a, T b) { return abs(b) > abs(a); });93 return abs(*it);94 }95 else if constexpr (std::is_unsigned<T>::value)96 {97 return *std::max_element(first, last);98 }99 else100 {101 auto pair = std::minmax_element(first, last);102 if (abs(*pair.first) > abs(*pair.second))103 {104 return abs(*pair.first);105 }106 else107 {108 return abs(*pair.second);109 }110 }111}112 113template<class Container>114inline auto sup_norm(Container const & v)115{116 return sup_norm(v.cbegin(), v.cend());117}118 119template<class ForwardIterator>120auto l1_norm(ForwardIterator first, ForwardIterator last)121{122 using T = typename std::iterator_traits<ForwardIterator>::value_type;123 using std::abs;124 if constexpr (std::is_unsigned<T>::value)125 {126 double l1 = 0;127 for (auto it = first; it != last; ++it)128 {129 l1 += *it;130 }131 return l1;132 }133 else if constexpr (std::is_integral<T>::value)134 {135 double l1 = 0;136 for (auto it = first; it != last; ++it)137 {138 double tmp = *it;139 l1 += abs(tmp);140 }141 return l1;142 }143 else144 {145 decltype(abs(*first)) l1 = 0;146 for (auto it = first; it != last; ++it)147 {148 l1 += abs(*it);149 }150 return l1;151 }152 153}154 155template<class Container>156inline auto l1_norm(Container const & v)157{158 return l1_norm(v.cbegin(), v.cend());159}160 161 162template<class ForwardIterator>163auto l2_norm(ForwardIterator first, ForwardIterator last)164{165 using T = typename std::iterator_traits<ForwardIterator>::value_type;166 using std::abs;167 using std::norm;168 using std::sqrt;169 using std::is_floating_point;170 using std::isfinite;171 if constexpr (boost::math::tools::is_complex_type<T>::value)172 {173 typedef typename T::value_type Real;174 Real l2 = 0;175 for (auto it = first; it != last; ++it)176 {177 l2 += norm(*it);178 }179 Real result = sqrt(l2);180 if (!isfinite(result))181 {182 Real a = sup_norm(first, last);183 l2 = 0;184 for (auto it = first; it != last; ++it)185 {186 l2 += norm(*it/a);187 }188 return a*sqrt(l2);189 }190 return result;191 }192 else if constexpr (is_floating_point<T>::value ||193 std::numeric_limits<T>::max_exponent)194 {195 T l2 = 0;196 for (auto it = first; it != last; ++it)197 {198 l2 += (*it)*(*it);199 }200 T result = sqrt(l2);201 // Higham, Accuracy and Stability of Numerical Algorithms,202 // Problem 27.5 presents a different algorithm to deal with overflow.203 // The algorithm used here takes 3 passes *if* there is overflow.204 // Higham's algorithm is 1 pass, but more requires operations than the no overflow case.205 // I'm operating under the assumption that overflow is rare since the dynamic range of floating point numbers is huge.206 if (!isfinite(result))207 {208 T a = sup_norm(first, last);209 l2 = 0;210 for (auto it = first; it != last; ++it)211 {212 T tmp = *it/a;213 l2 += tmp*tmp;214 }215 return a*sqrt(l2);216 }217 return result;218 }219 else220 {221 double l2 = 0;222 for (auto it = first; it != last; ++it)223 {224 double tmp = *it;225 l2 += tmp*tmp;226 }227 return sqrt(l2);228 }229}230 231template<class Container>232inline auto l2_norm(Container const & v)233{234 return l2_norm(v.cbegin(), v.cend());235}236 237template<class ForwardIterator>238size_t l0_pseudo_norm(ForwardIterator first, ForwardIterator last)239{240 using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;241 size_t count = 0;242 for (auto it = first; it != last; ++it)243 {244 if (*it != RealOrComplex(0))245 {246 ++count;247 }248 }249 return count;250}251 252template<class Container>253inline size_t l0_pseudo_norm(Container const & v)254{255 return l0_pseudo_norm(v.cbegin(), v.cend());256}257 258template<class ForwardIterator>259size_t hamming_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)260{261 size_t count = 0;262 auto it1 = first1;263 auto it2 = first2;264 while (it1 != last1)265 {266 if (*it1++ != *it2++)267 {268 ++count;269 }270 }271 return count;272}273 274template<class Container>275inline size_t hamming_distance(Container const & v, Container const & w)276{277 return hamming_distance(v.cbegin(), v.cend(), w.cbegin());278}279 280template<class ForwardIterator>281auto lp_norm(ForwardIterator first, ForwardIterator last, unsigned p)282{283 using std::abs;284 using std::pow;285 using std::is_floating_point;286 using std::isfinite;287 using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;288 if constexpr (boost::math::tools::is_complex_type<RealOrComplex>::value)289 {290 using std::norm;291 using Real = typename RealOrComplex::value_type;292 Real lp = 0;293 for (auto it = first; it != last; ++it)294 {295 lp += pow(abs(*it), p);296 }297 298 auto result = pow(lp, Real(1)/Real(p));299 if (!isfinite(result))300 {301 auto a = boost::math::tools::sup_norm(first, last);302 Real lp = 0;303 for (auto it = first; it != last; ++it)304 {305 lp += pow(abs(*it)/a, p);306 }307 result = a*pow(lp, Real(1)/Real(p));308 }309 return result;310 }311 else if constexpr (is_floating_point<RealOrComplex>::value || std::numeric_limits<RealOrComplex>::max_exponent)312 {313 BOOST_MATH_ASSERT_MSG(p >= 0, "For p < 0, the lp norm is not a norm");314 RealOrComplex lp = 0;315 316 for (auto it = first; it != last; ++it)317 {318 lp += pow(abs(*it), p);319 }320 321 RealOrComplex result = pow(lp, RealOrComplex(1)/RealOrComplex(p));322 if (!isfinite(result))323 {324 RealOrComplex a = boost::math::tools::sup_norm(first, last);325 lp = 0;326 for (auto it = first; it != last; ++it)327 {328 lp += pow(abs(*it)/a, p);329 }330 result = a*pow(lp, RealOrComplex(1)/RealOrComplex(p));331 }332 return result;333 }334 else335 {336 double lp = 0;337 338 for (auto it = first; it != last; ++it)339 {340 double tmp = *it;341 lp += pow(abs(tmp), p);342 }343 double result = pow(lp, 1.0/static_cast<double>(p));344 if (!isfinite(result))345 {346 double a = boost::math::tools::sup_norm(first, last);347 lp = 0;348 for (auto it = first; it != last; ++it)349 {350 double tmp = *it;351 lp += pow(abs(tmp)/a, p);352 }353 result = a*pow(lp, static_cast<double>(1)/static_cast<double>(p));354 }355 return result;356 }357}358 359template<class Container>360inline auto lp_norm(Container const & v, unsigned p)361{362 return lp_norm(v.cbegin(), v.cend(), p);363}364 365 366template<class ForwardIterator>367auto lp_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2, unsigned p)368{369 using std::pow;370 using std::abs;371 using std::is_floating_point;372 using std::isfinite;373 using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;374 auto it1 = first1;375 auto it2 = first2;376 377 if constexpr (boost::math::tools::is_complex_type<RealOrComplex>::value)378 {379 using Real = typename RealOrComplex::value_type;380 using std::norm;381 Real dist = 0;382 while(it1 != last1)383 {384 auto tmp = *it1++ - *it2++;385 dist += pow(abs(tmp), p);386 }387 return pow(dist, Real(1)/Real(p));388 }389 else if constexpr (is_floating_point<RealOrComplex>::value || std::numeric_limits<RealOrComplex>::max_exponent)390 {391 RealOrComplex dist = 0;392 while(it1 != last1)393 {394 auto tmp = *it1++ - *it2++;395 dist += pow(abs(tmp), p);396 }397 return pow(dist, RealOrComplex(1)/RealOrComplex(p));398 }399 else400 {401 double dist = 0;402 while(it1 != last1)403 {404 double tmp1 = *it1++;405 double tmp2 = *it2++;406 // Naively you'd expect the integer subtraction to be faster,407 // but this can overflow or wraparound:408 //double tmp = *it1++ - *it2++;409 dist += pow(abs(tmp1 - tmp2), p);410 }411 return pow(dist, 1.0/static_cast<double>(p));412 }413}414 415template<class Container>416inline auto lp_distance(Container const & v, Container const & w, unsigned p)417{418 return lp_distance(v.cbegin(), v.cend(), w.cbegin(), p);419}420 421 422template<class ForwardIterator>423auto l1_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)424{425 using std::abs;426 using std::is_floating_point;427 using std::isfinite;428 using T = typename std::iterator_traits<ForwardIterator>::value_type;429 auto it1 = first1;430 auto it2 = first2;431 if constexpr (boost::math::tools::is_complex_type<T>::value)432 {433 using Real = typename T::value_type;434 Real sum = 0;435 while (it1 != last1) {436 sum += abs(*it1++ - *it2++);437 }438 return sum;439 }440 else if constexpr (is_floating_point<T>::value || std::numeric_limits<T>::max_exponent)441 {442 T sum = 0;443 while (it1 != last1)444 {445 sum += abs(*it1++ - *it2++);446 }447 return sum;448 }449 else if constexpr (std::is_unsigned<T>::value)450 {451 double sum = 0;452 while(it1 != last1)453 {454 T x1 = *it1++;455 T x2 = *it2++;456 if (x1 > x2)457 {458 sum += (x1 - x2);459 }460 else461 {462 sum += (x2 - x1);463 }464 }465 return sum;466 }467 else if constexpr (std::is_integral<T>::value)468 {469 double sum = 0;470 while(it1 != last1)471 {472 double x1 = *it1++;473 double x2 = *it2++;474 sum += abs(x1-x2);475 }476 return sum;477 }478 else479 {480 BOOST_MATH_ASSERT_MSG(false, "Could not recognize type.");481 }482 483}484 485template<class Container>486auto l1_distance(Container const & v, Container const & w)487{488 using std::size;489 BOOST_MATH_ASSERT_MSG(size(v) == size(w),490 "L1 distance requires both containers to have the same number of elements");491 return l1_distance(v.cbegin(), v.cend(), w.begin());492}493 494template<class ForwardIterator>495auto l2_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)496{497 using std::abs;498 using std::norm;499 using std::sqrt;500 using std::is_floating_point;501 using std::isfinite;502 using T = typename std::iterator_traits<ForwardIterator>::value_type;503 auto it1 = first1;504 auto it2 = first2;505 if constexpr (boost::math::tools::is_complex_type<T>::value)506 {507 using Real = typename T::value_type;508 Real sum = 0;509 while (it1 != last1) {510 sum += norm(*it1++ - *it2++);511 }512 return sqrt(sum);513 }514 else if constexpr (is_floating_point<T>::value || std::numeric_limits<T>::max_exponent)515 {516 T sum = 0;517 while (it1 != last1)518 {519 T tmp = *it1++ - *it2++;520 sum += tmp*tmp;521 }522 return sqrt(sum);523 }524 else if constexpr (std::is_unsigned<T>::value)525 {526 double sum = 0;527 while(it1 != last1)528 {529 T x1 = *it1++;530 T x2 = *it2++;531 if (x1 > x2)532 {533 double tmp = x1-x2;534 sum += tmp*tmp;535 }536 else537 {538 double tmp = x2 - x1;539 sum += tmp*tmp;540 }541 }542 return sqrt(sum);543 }544 else545 {546 double sum = 0;547 while(it1 != last1)548 {549 double x1 = *it1++;550 double x2 = *it2++;551 double tmp = x1-x2;552 sum += tmp*tmp;553 }554 return sqrt(sum);555 }556}557 558template<class Container>559auto l2_distance(Container const & v, Container const & w)560{561 using std::size;562 BOOST_MATH_ASSERT_MSG(size(v) == size(w),563 "L2 distance requires both containers to have the same number of elements");564 return l2_distance(v.cbegin(), v.cend(), w.begin());565}566 567template<class ForwardIterator>568auto sup_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)569{570 using std::abs;571 using std::norm;572 using std::sqrt;573 using std::is_floating_point;574 using std::isfinite;575 using T = typename std::iterator_traits<ForwardIterator>::value_type;576 auto it1 = first1;577 auto it2 = first2;578 if constexpr (boost::math::tools::is_complex_type<T>::value)579 {580 using Real = typename T::value_type;581 Real sup_sq = 0;582 while (it1 != last1) {583 Real tmp = norm(*it1++ - *it2++);584 if (tmp > sup_sq) {585 sup_sq = tmp;586 }587 }588 return sqrt(sup_sq);589 }590 else if constexpr (is_floating_point<T>::value || std::numeric_limits<T>::max_exponent)591 {592 T sup = 0;593 while (it1 != last1)594 {595 T tmp = *it1++ - *it2++;596 if (sup < abs(tmp))597 {598 sup = abs(tmp);599 }600 }601 return sup;602 }603 else // integral values:604 {605 double sup = 0;606 while(it1 != last1)607 {608 T x1 = *it1++;609 T x2 = *it2++;610 double tmp;611 if (x1 > x2)612 {613 tmp = x1-x2;614 }615 else616 {617 tmp = x2 - x1;618 }619 if (sup < tmp) {620 sup = tmp;621 }622 }623 return sup;624 }625}626 627template<class Container>628auto sup_distance(Container const & v, Container const & w)629{630 using std::size;631 BOOST_MATH_ASSERT_MSG(size(v) == size(w),632 "sup distance requires both containers to have the same number of elements");633 return sup_distance(v.cbegin(), v.cend(), w.begin());634}635 636 637}638#endif639