Implemented a numerical quadrature method based on an adaptive

Simpson rule.  Added unit tests and supporting examples for this
function.

dlib/numerical_integration/integrate_function_adapt_simpson.h

@@ -4,14 +4,39 @@
 #ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__
 #define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__
 
-#include "../matrix.h"
-
-class adapt_simp
+template <typename T, typename funct>
+T impl_adapt_simp_stop(const funct& f, T a, T b, T fa, T fm, T fb, T is, int cnt)
 {
-public:
+    int MAXINT = 1000;
+   
+    T m   = (a + b)/2.0;
+    T h   = (b - a)/4.0;
+    T fml = f(a + h);
+    T fmr = f(b - h);
+    T i1 = h/1.5*(fa+4.0*fm+fb);
+    T i2 = h/3.0*(fa+4.0*(fml+fmr)+2.0*fm+fb);
+    i1 = (16.0*i2 - i1)/15.0;
+    T Q = 0;
+
+    if((std::abs(i1-i2) <= std::numeric_limits<double>::epsilon()) || (m <= a) || (b <= m))
+    {
+        Q = i1;
+    }
+    else 
+    {
+        if(cnt < MAXINT)
+        {cnt = cnt + 1;
+
+            Q = impl_adapt_simp_stop(f,a,m,fa,fml,fm,is,cnt) 
+              + impl_adapt_simp_stop(f,m,b,fm,fmr,fb,is,cnt); 
+        }
+    }
+
+    return Q;
+}
 
 template <typename T, typename funct>
-T integrate_function_adapt_simp(const funct& f, T a, T b, T tol)
+T integrate_function_adapt_simp(const funct& f, T a, T b, T tol = 1e-10)
 {
     T eps = std::numeric_limits<double>::epsilon();
 
@@ -37,52 +62,11 @@ T integrate_function_adapt_simp(const funct& f, T a, T b, T tol)
 
     int cnt = 0;
 
-    T tstvl = adapt_simp_stop(f, a, b, fa, fm, fb, is, cnt);
+    T tstvl = impl_adapt_simp_stop(f, a, b, fa, fm, fb, is, cnt);
 
     return tstvl;
 
 }
 
-private: 
-
-template <typename T, typename funct>
-T adapt_simp_stop(const funct& f, T a, T b, T fa, T fm, T fb, T is, int &cnt)
-{
-    int MAXINT = 1000;
-   
-    T m   = (a + b)/2.0;
-    T h   = (b - a)/4.0;
-    T fml = f(a + h);
-    T fmr = f(b - h);
-    T i1 = h/1.5*(fa+4.0*fm+fb);
-    T i2 = h/3.0*(fa+4.0*(fml+fmr)+2.0*fm+fb);
-    i1 = (16.0*i2 - i1)/15.0;
-    T Q = 0;
-
-    if((is+(i1-i2) == is) || (m <= a) || (b <= m))
-    {
-        if((m <= a) || (b <= m))
-        {
-            DLIB_ASSERT(1, "\tintegrate_function_adapt_simpson::adapt_simp_stop"
-                           << "\n\tmidpoint evaluation occurred at endpoint");
-        }
-    
-        Q = i1;
-    }
-    else 
-    {
-        if(cnt < MAXINT)
-        {cnt = cnt + 1;
-
-            Q = adapt_simp_stop(f,a,m,fa,fml,fm,is,cnt) 
-              + adapt_simp_stop(f,m,b,fm,fmr,fb,is,cnt); 
-        }
-    }
-
-    return Q;
-}
-
-};
-
 #endif //DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON.h__
 

dlib/numerical_integration/integrate_function_adapt_simpson_abstract.h

@@ -3,16 +3,17 @@
 #ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__
 #define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__
 
-#include "matrix.h"
-
 templte <typename T, typename funct>
 T integrate_function_adapt_simp(const funct& f, T a, T b, T tol);
 /*!
     requires 
         - b > a
-        - tol > 0
-        - f to be real valued single variable function
+        - tol > 0, tol is a tolerance parameter that typically determines 
+                   the overall accuracy of approximated integral.  We suggest 
+                   a default value of 1e-10 for tol. 
 
+        - f to be real valued single variable function
+    
     ensures
         - returns an approximation of the integral of f over the domain [a,b] 
           using the adaptive Simpson method outlined in 

dlib/test/CMakeLists.txt

@@ -81,6 +81,7 @@ set (tests
    member_function_pointer.cpp
    metaprogramming.cpp
    multithreaded_object.cpp
+   numerical_integration.cpp
    object_detector.cpp
    oca.cpp
    one_vs_all_trainer.cpp

dlib/test/makefile

@@ -96,6 +96,7 @@ SRC += md5.cpp
 SRC += member_function_pointer.cpp
 SRC += metaprogramming.cpp
 SRC += multithreaded_object.cpp
+SRC += numerical_integration.cpp
 SRC += object_detector.cpp
 SRC += oca.cpp
 SRC += one_vs_all_trainer.cpp

dlib/test/numerical_integration.cpp

@@ -40,45 +40,43 @@ namespace
 
             matrix<double,23,1> m;
             double tol = 1e-10;
-            double eps = 1e-4;
-
-            adapt_simp ad;
-
-            m(0) = ad.integrate_function_adapt_simp(&gg1, 0.0, 1.0, tol);
-            m(1) = ad.integrate_function_adapt_simp(&gg2, 0.0, 1.0, tol);
-            m(2) = ad.integrate_function_adapt_simp(&gg3, 0.0, 1.0, tol);
-            m(3) = ad.integrate_function_adapt_simp(&gg4, 0.0, 1.0, tol);
-            m(4) = ad.integrate_function_adapt_simp(&gg5, -1.0, 1.0, tol);
-            m(5) = ad.integrate_function_adapt_simp(&gg6, 0.0, 1.0, tol);
-            m(6) = ad.integrate_function_adapt_simp(&gg7, 0.0, 1.0, tol);
-            m(7) = ad.integrate_function_adapt_simp(&gg8, 0.0, 1.0, tol);
-            m(8) = ad.integrate_function_adapt_simp(&gg9, 0.0, 1.0, tol);
-            m(9) = ad.integrate_function_adapt_simp(&gg10, 0.0, 1.0, tol);
-            m(10) = ad.integrate_function_adapt_simp(&gg11, 0.0, 1.0, tol);
-            m(11) = ad.integrate_function_adapt_simp(&gg12, 1e-6, 1.0, tol);
-            m(12) = ad.integrate_function_adapt_simp(&gg13, 0.0, 10.0, tol);
-            m(13) = ad.integrate_function_adapt_simp(&gg14, 0.0, 10.0, tol);
-            m(14) = ad.integrate_function_adapt_simp(&gg15, 0.0, 10.0, tol);
-            m(15) = ad.integrate_function_adapt_simp(&gg16, 0.01, 1.0, tol);
-            m(16) = ad.integrate_function_adapt_simp(&gg17, 0.0, pi, tol);
-            m(17) = ad.integrate_function_adapt_simp(&gg18, 0.0, 1.0, tol);
-            m(18) = ad.integrate_function_adapt_simp(&gg19, -1.0, 1.0, tol);
-            m(19) = ad.integrate_function_adapt_simp(&gg20, 0.0, 1.0, tol);
-            m(20) = ad.integrate_function_adapt_simp(&gg21, 0.0, 1.0, tol);
-            m(21) = ad.integrate_function_adapt_simp(&gg22, 0.0, 5.0, tol);
+            double eps = 1e-8;
+
+            m(0) = integrate_function_adapt_simp(&gg1, 0.0, 1.0, tol);
+            m(1) = integrate_function_adapt_simp(&gg2, 0.0, 1.0, tol);
+            m(2) = integrate_function_adapt_simp(&gg3, 0.0, 1.0, tol);
+            m(3) = integrate_function_adapt_simp(&gg4, 0.0, 1.0, tol);
+            m(4) = integrate_function_adapt_simp(&gg5, -1.0, 1.0, tol);
+            m(5) = integrate_function_adapt_simp(&gg6, 0.0, 1.0, tol);
+            m(6) = integrate_function_adapt_simp(&gg7, 0.0, 1.0, tol);
+            m(7) = integrate_function_adapt_simp(&gg8, 0.0, 1.0, tol);
+            m(8) = integrate_function_adapt_simp(&gg9, 0.0, 1.0, tol);
+            m(9) = integrate_function_adapt_simp(&gg10, 0.0, 1.0, tol);
+            m(10) = integrate_function_adapt_simp(&gg11, 0.0, 1.0, tol);
+            m(11) = integrate_function_adapt_simp(&gg12, 1e-6, 1.0, tol);
+            m(12) = integrate_function_adapt_simp(&gg13, 0.0, 10.0, tol);
+            m(13) = integrate_function_adapt_simp(&gg14, 0.0, 10.0, tol);
+            m(14) = integrate_function_adapt_simp(&gg15, 0.0, 10.0, tol);
+            m(15) = integrate_function_adapt_simp(&gg16, 0.01, 1.0, tol);
+            m(16) = integrate_function_adapt_simp(&gg17, 0.0, pi, tol);
+            m(17) = integrate_function_adapt_simp(&gg18, 0.0, 1.0, tol);
+            m(18) = integrate_function_adapt_simp(&gg19, -1.0, 1.0, tol);
+            m(19) = integrate_function_adapt_simp(&gg20, 0.0, 1.0, tol);
+            m(20) = integrate_function_adapt_simp(&gg21, 0.0, 1.0, tol);
+            m(21) = integrate_function_adapt_simp(&gg22, 0.0, 5.0, tol);
 
             // Here we compare the approximated integrals against 
             // highly accurate approximations generated either from 
             // the exact integral values or Mathematica's NIntegrate 
             // function using a working precision of 20. 
 
-            DLIB_TEST(abs(m(0) - 1.7182818284590452354) < eps);
+            DLIB_TEST(abs(m(0) - 1.7182818284590452354) < 1e-11);
             DLIB_TEST(abs(m(1) - 0.7000000000000000000) < eps);
             DLIB_TEST(abs(m(2) - 0.6666666666666666667) < eps);
             DLIB_TEST(abs(m(3) - 0.2397141133444008336) < eps);
-            DLIB_TEST(abs(m(4) - 1.5822329637296729331) < eps);
+            DLIB_TEST(abs(m(4) - 1.5822329637296729331) < 1e-12);
             DLIB_TEST(abs(m(5) - 0.4000000000000000000) < eps);
-            DLIB_TEST(abs(m(6) - 2.0000000000000000000) < eps);
+            DLIB_TEST(abs(m(6) - 2.0000000000000000000) < 1e-4);
             DLIB_TEST(abs(m(7) - 0.8669729873399110375) < eps);
             DLIB_TEST(abs(m(8) - 1.1547005383792515290) < eps);
             DLIB_TEST(abs(m(9) - 0.6931471805599453094) < eps);
@@ -89,7 +87,7 @@ namespace
             DLIB_TEST(abs(m(14) - 0.4993633810764567446) < eps);
             DLIB_TEST(abs(m(15) - 0.1121393035410217   ) < eps);
             DLIB_TEST(abs(m(16) - 0.2910187828600526985) < eps);
-            DLIB_TEST(abs(m(17) + 0.4342944819032518276) < eps);
+            DLIB_TEST(abs(m(17) + 0.4342944819032518276) < 1e-5);
             DLIB_TEST(abs(m(18) - 1.56439644406905     ) < eps);
             DLIB_TEST(abs(m(19) - 0.1634949430186372261) < eps);
             DLIB_TEST(abs(m(20) - 0.0134924856494677726) < eps);

examples/integrate_function_adapt_simp_ex.cpp

-// Copyright (C) 2013 Steve Taylor (steve98654@gmai.com)
 // The contents of this file are in the public domain.  See
 // LICENSE_FOR_EXAMPLE_PROGRAMS.txt
 
@@ -28,50 +27,40 @@ using namespace dlib;
 // Here we define a class that consists of the set of functions that we 
 // wish to integrate and comment in the domain of integration.
 
-class int_fcts {
-    
-public:
-
-    // x in [0,1]
-    static double gg1(double x)
-    {
-        return pow(e,x);
-    }   
-
-    // x in [0,1]
-    static double gg2(double x)
-    {
-        return x*x;
-    }
-
-    // x in [0, pi]
-    static double gg3(double x)
-    {
-        return 1/(x*x + cos(x)*cos(x));
-    }
-
-    // x in [-pi, pi]
-    static double gg4(double x)
-    {
-        return sin(x);
-    }
-
-    // x in [0,2]
-    static double gg5(double x)
-    {
-        return 1/(1 + x*x);
-    }
-};
+// x in [0,1]
+static double gg1(double x)
+{
+    return pow(e,x);
+}   
+
+// x in [0,1]
+static double gg2(double x)
+{
+    return x*x;
+}
+
+// x in [0, pi]
+static double gg3(double x)
+{
+    return 1/(x*x + cos(x)*cos(x));
+}
+
+// x in [-pi, pi]
+static double gg4(double x)
+{
+    return sin(x);
+}
+
+// x in [0,2]
+static double gg5(double x)
+{
+    return 1/(1 + x*x);
+}
 
 // Examples 
 int main()
 {
-    // We first define a matrix m to which we will store the approximated values of the
-    // integrals of our int_fcts functions.
-    
-    matrix<double,5,1> m;
-    
-    // Next we define a tolerance parameter.  Roughly speaking, a lower tolerance will
+    // We first define a tolerance parameter.  Roughly speaking, a lower tolerance will
     // result in a more accurate approximation of the true integral.  However, there are 
     // instances where too small of a tolerance may yield a less accurate approximation
     // than a larger tolerance.  We recommend taking the tolerance to be in the
@@ -80,25 +69,25 @@ int main()
     double tol = 1e-10;
 
     // Here we instantiate a class which contains the numerical quadrature method.
-    adapt_simp ad;
+    //    adapt_simp ad;
 
     // Here we compute the integrals of the five functions defined above using the same 
     // tolerance level for each.
 
-    m(0) = ad.integrate_function_adapt_simp(&int_fcts::gg1, 0.0, 1.0, tol);
-    m(1) = ad.integrate_function_adapt_simp(&int_fcts::gg2, 0.0, 1.0, tol);
-    m(2) = ad.integrate_function_adapt_simp(&int_fcts::gg3, 0.0, pi, tol);
-    m(3) = ad.integrate_function_adapt_simp(&int_fcts::gg4, -pi, pi, tol);
-    m(4) = ad.integrate_function_adapt_simp(&int_fcts::gg5, 0.0, 2.0, tol);
+    double m1 = integrate_function_adapt_simp(&gg1, 0.0, 1.0, tol);
+    double m2 = integrate_function_adapt_simp(&gg2, 0.0, 1.0, tol);
+    double m3 = integrate_function_adapt_simp(&gg3, 0.0, pi, tol);
+    double m4 = integrate_function_adapt_simp(&gg4, -pi, pi, tol);
+    double m5 = integrate_function_adapt_simp(&gg5, 0.0, 2.0, tol);
 
     // We finally print out the values of each of the approximated integrals to ten
     // significant digits.
 
-    cout << "\nThe integral of exp(x) for x in [0,1] is "          << std::setprecision(10) <<  m(0)  << endl; 
-    cout << "The integral of x^2 for in [0,1] is "                 << std::setprecision(10) <<  m(1)  << endl; 
-    cout << "The integral of 1/(x^2 + cos(x)^2) for in [0,pi] is " << std::setprecision(10) <<  m(2)  << endl;
-    cout << "The integral of sin(x) for in [-pi,pi] is "           << std::setprecision(10) <<  m(3)  << endl;
-    cout << "The integral of 1/(1+x^2) for in [0,2] is "           << std::setprecision(10) <<  m(4)  << endl;
+    cout << "\nThe integral of exp(x) for x in [0,1] is "          << std::setprecision(10) <<  m1  << endl; 
+    cout << "The integral of x^2 for in [0,1] is "                 << std::setprecision(10) <<  m2  << endl; 
+    cout << "The integral of 1/(x^2 + cos(x)^2) for in [0,pi] is " << std::setprecision(10) <<  m3  << endl;
+    cout << "The integral of sin(x) for in [-pi,pi] is "           << std::setprecision(10) <<  m4  << endl;
+    cout << "The integral of 1/(1+x^2) for in [0,2] is "           << std::setprecision(10) <<  m5  << endl;
     cout << "" << endl;
 
     return 0;