Document wavelet transform. Submit (mildly) failing unit test for transforms. [CI SKIP]

doc/math.qbk

@@ -727,6 +727,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include quadrature/double_exponential.qbk]
 [include quadrature/ooura_fourier_integrals.qbk]
 [include quadrature/naive_monte_carlo.qbk]
+[include quadrature/wavelet_transforms.qbk]
 [include differentiation/numerical_differentiation.qbk]
 [include differentiation/autodiff.qbk]
 [include differentiation/lanczos_smoothing.qbk]

doc/quadrature/wavelet_transforms.qbk

@@ -0,0 +1,56 @@
+[/
+Copyright (c) 2019 Nick Thompson
+Copyright (c) 2019 Paul A. Bristow
+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)
+]
+
+[section:wavelet_transforms Wavelet Transforms]
+
+[heading Synopsis]
+
+```
+    #include <boost/math/quadrature/wavelet_transforms.hpp>
+    
+    namespace boost::math::quadrature {
+
+    template<class F, typename Real, int p>
+    class daubechies_wavelet_transform
+    {
+    public:
+        daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = 100*std::numeric_limits<Real>::epsilon(),
+        int max_refinements = 12) {}
+
+        daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet<Real, p> wavelet, Real tol = 100*std::numeric_limits<Real>::epsilon(),
+        int max_refinements = 12);
+
+        auto operator()(Real s, Real t)->decltype(std::declval<F>()(std::declval<Real>())) const;
+
+    };
+
+    } 
+```
+
+Wavelet transforms compute
+
+[$../graphs/wavelet_transform_definition.svg]
+
+For compactly supported Daubechies wavelets, the bounds can always be taken as finite, and we have 
+
+[$../graphs/daubechies_wavelet_transform_definition.svg]
+
+which gives definies the /s=0/ case.
+
+A basic usage is 
+
+    auto psi = daubechies_wavelet<double, 8>();
+    auto f = [](double x) {
+        return sin(2/x);
+    };
+    auto Wf = daubechies_wavelet_transform(f, psi);
+
+    auto z = Wf(0.8, 7.2);
+
+[endsect] [/section:wavelet_transforms]
+

include/boost/math/quadrature/wavelet_transforms.hpp

@@ -15,11 +15,11 @@ template<class F, typename Real, int p>
 class daubechies_wavelet_transform
 {
 public:
-    daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = boost::math::tools::root_epsilon<Real>(),
+    daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = 100*std::numeric_limits<Real>::epsilon(),
     int max_refinements = 12) : f_{f}, psi_(grid_refinements), tol_{tol}, max_refinements_{max_refinements}
     {}
 
-    daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet<Real, p> wavelet, Real tol = boost::math::tools::root_epsilon<Real>(),
+    daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet<Real, p> wavelet, Real tol = 100*std::numeric_limits<Real>::epsilon(),
     int max_refinements = 12) : f_{f}, psi_{wavelet}, tol_{tol}, max_refinements_{max_refinements}
     {}
 

test/wavelet_transform_test.cpp

@@ -27,6 +27,7 @@ using boost::multiprecision::float128;
 using boost::math::constants::pi;
 using boost::math::constants::root_two;
 using boost::math::quadrature::daubechies_wavelet_transform;
+using boost::math::quadrature::trapezoidal;
 
 template<typename Real, int p>
 void test_wavelet_transform()
@@ -38,21 +39,26 @@ void test_wavelet_transform()
         return abs(psi(x));
     };
     auto [a, b] = psi.support();
-    auto psil1 = boost::math::quadrature::trapezoidal(abs_psi, a, b);
-    std::cout << "L1 norm of psi = " << psil1 << "\n";
-    auto [x_min, psi_min] = boost::math::tools::brent_find_minima(psi, a, b, std::numeric_limits<Real>::digits10);
-    std::cout << "Minimum value is " << psi_min <<  "  which occurs at x = " << x_min << "\n";
-    auto neg_psi = [&](Real x) { return -psi(x); };
-    auto [x_max, neg_psi_max] = boost::math::tools::brent_find_minima(neg_psi, a, b, std::numeric_limits<Real>::digits10);
-    std::cout << "Maximum value of psi is " << -neg_psi_max << "\n";
-    Real psi_sup_norm = std::max(abs(psi_min), std::abs(neg_psi_max));
+    auto psil1 = trapezoidal(abs_psi, a, b);
+    Real psi_sup_norm = 0;
+    for (double x = a; x < b; x += 0.0001)
+    {
+        Real y = psi(x);
+        if (std::abs(y) > psi_sup_norm)
+        {
+            psi_sup_norm = std::abs(y);
+        }
+    }
     std::cout << "psi sup norm = " << psi_sup_norm << "\n";
     // An even function:
     auto f = [](Real x) {
-        return std::exp(-abs(x))*std::cos(10*x);
+        return std::exp(-abs(x));
     };
     Real fmax = 1;
-
+    Real fl2 = 1;
+    Real fl1 = 2;
+    std::cout << "||f||_1 = " << fl1 << "\n";
+    std::cout << "||f||_2 = " << fl2 << "\n";
 
     auto Wf = daubechies_wavelet_transform(f, psi);
     for (double s = 0; s < 10; s += 0.01)
@@ -61,22 +67,58 @@ void test_wavelet_transform()
         Real w2 = Wf(-s, 0.0);
         // Since f is an even function, we get w1 = w2:
         CHECK_ULP_CLOSE(w1, w2, 12);
-        // Integral inequality:
-        Real r1 = sqrt(abs(s))*fmax*psil1;
-        //std::cout << "|w| = " << abs(w1) << ", ||f||_infty||psi||_1 = " << r1 << "\n";
-        if(abs(w1) > r1)
+    }
+
+    // The wavelet transform with respect to Daubechies wavelets 
+    for (double s = -10; s < 10; s += 0.1)
+    {
+        for (double t = -10; t < 10; t+= 0.1)
         {
-            std::cerr << " Integral inequality | int fg| <= ||f||_infty ||g||_infty is violated.\n";
+            Real w = Wf(s, t);
+            // Integral inequality:
+            Real r1 = sqrt(abs(s))*fmax*psil1;
+            if(abs(w) > r1)
+            {
+                std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1 is violated.\n";
+            }
+            if (s != 0)
+            {
+                Real r2 = fl1*psi_sup_norm/sqrt(abs(s));
+                if(abs(w) > r2)
+                {
+                    std::cerr << " Integral inequality | int fg| <= ||f||_1 ||psi||_infty/sqrt(|s|) is violated.\n";
+                    std::cerr << " Violation: " << abs(w) << " !<= " << r2 << "\n";
+                }
+                Real r3 = fmax*psil1/sqrt(abs(s));
+                if(abs(w) > r3)
+                {
+                    std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1/sqrt(|s|) is violated.\n";
+                    std::cerr << " Computed = " << abs(w) << ", expected " << r3 << "\n";
+                }
+            }
+            if (abs(w) > fl2)
+            {
+                std::cerr << "  Integral inequality |f psi_s,t| <= ||f||_2 ||psi||2 violated.\n";
+            }
+            Real r4 = sqrt(abs(s))*fl1*psi_sup_norm;
+            if (abs(w) > r4)
+            {
+                std::cerr << "  Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_1 ||psi||_infty is violated.\n";
+            }
+            Real r5 = sqrt(abs(s))*fmax*psil1;
+            if (abs(w) > r5)
+            {
+                std::cerr << "  Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_infty ||psi||_1 is violated.\n";
+            }
+
         }
     }
-    
 }
 
 int main()
 {
-    test_wavelet_transform<double, 8>();
-    /*boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
+    boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
         test_wavelet_transform<double, i+3>();
-    });*/
+    });
     return boost::math::test::report_errors();
 }