Commit ad427072 authored by Davis King's avatar Davis King

merged

parents 66d5a906 0694c005
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
#include "matrix.h"
#include "numerical_integration/integrate_function_adapt_simpson.h"
#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
......@@ -12,13 +12,16 @@ namespace dlib
// e -- Euler's Constant
const double e = 2.7182818284590452354;
// sqrt_2 -- the square root of 2
// sqrt_2 -- The square root of 2
const double sqrt_2 = 1.4142135623730950488;
// sqrt_3 -- the square root of 3
// sqrt_3 -- The square root of 3
const double sqrt_3 = 1.7320508075688772935;
// light_spd -- the speed of light in vacuum in meters per second
// log10_2 -- The logarithm base 10 of two
const double log10_2 = 0.30102999566398119521;
// light_spd -- The speed of light in vacuum in meters per second
const double light_spd = 2.99792458e8;
// newton_G -- Newton's gravitational constant (in metric units of m^3/(kg*s^2))
......@@ -27,6 +30,23 @@ namespace dlib
// planck_cst -- Planck's constant (in units of Joules * seconds)
const double planck_cst = 6.62606957e-34;
// golden_ratio -- The Golden Ratio
const double golden_ratio = 1.6180339887498948482;
// euler_gamma -- The Euler Mascheroni Constant
const double euler_gamma = 0.5772156649015328606065;
// catalan -- Catalan's Constant
const double catalan = 0.91596559417721901505;
// glaisher -- Glaisher Kinkelin constant
const double glaisher = 1.2824271291006226369;
// khinchin -- Khinchin's constant
const double khinchin = 2.6854520010653064453;
// apery -- Apery's constant
const double apery = 1.2020569031595942854;
}
#endif //DLIB_NUMERIC_CONSTANTs_H_
......
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
#include "matrix.h"
#include "numerical_integration/integrate_function_adapt_simpson.h"
#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for full license
#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__
#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__
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)
{
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::abs(is)) || (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 eps = std::numeric_limits<double>::epsilon();
if(tol < eps)
{
tol = eps;
}
const T ba = b-a;
const T fa = f(a);
const T fb = f(b);
const T fm = f((a+b)/2);
T is =ba/8*(fa+fb+fm+ f(a + 0.9501*ba) + f(a + 0.2311*ba) + f(a + 0.6068*ba)
+ f(a + 0.4860*ba) + f(a + 0.8913*ba));
if(is == 0)
{
is = b-a;
}
is = is*tol;
int cnt = 0;
T tstvl = impl_adapt_simp_stop(f, a, b, fa, fm, fb, is, cnt);
return tstvl;
}
#endif //DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON.h__
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__
#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__
template <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
ensures
- returns an approximation of the integral of f over the domain [a,b]
using the adaptive Simpson method outlined in
Gander, W. and W. Gautshi, "Adaptive Quadrature -- Revisited"
BIT, Vol. 40, (2000), pp.84-101
- tol is a tolerance parameter that typically determines
the overall accuracy of approximated integral. We suggest
a default value of 1e-10 for tol.
!*/
#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__
......@@ -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
......
......@@ -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
......
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for the full license.
// This function test battery is given in:
//
// Test functions taken from Pedro Gonnet's dissertation at ETH:
// Adaptive Quadrature Re-Revisited
// http://e-collection.library.ethz.ch/eserv/eth:65/eth-65-02.pdf
#include <math.h>
#include <dlib/matrix.h>
#include <dlib/numeric_constants.h>
#include <dlib/numerical_integration.h>
#include "tester.h"
namespace
{
using namespace test;
using namespace dlib;
using namespace std;
logger dlog("test.numerical_integration");
class numerical_integration_tester : public tester
{
public:
numerical_integration_tester (
) :
tester ("test_numerical_integration",
"Runs tests on the numerical integration function.",
0
)
{}
void perform_test()
{
dlog <<dlib::LINFO << "Testing integrate_function_adapt_simpson";
matrix<double,23,1> m;
double tol = 1e-10;
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) < 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) < 1e-11);
DLIB_TEST(abs(m(5) - 0.4000000000000000000) < 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);
DLIB_TEST(abs(m(10) - 0.3798854930417224753) < eps);
DLIB_TEST(abs(m(11) - 0.7775036341124982763) < eps);
DLIB_TEST(abs(m(12) - 0.5000000000000000000) < eps);
DLIB_TEST(abs(m(13) - 1.0000000000000000000) < eps);
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) < 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);
}
static double gg1(double x)
{
return pow(e,x);
}
static double gg2(double x)
{
if(x > 0.3)
{
return 1.0;
}
else
{
return 0;
}
}
static double gg3(double x)
{
return pow(x,0.5);
}
static double gg4(double x)
{
return 23.0/25.0*cosh(x)-cos(x);
}
static double gg5(double x)
{
return 1/(pow(x,4) + pow(x,2) + 0.9);
}
static double gg6(double x)
{
return pow(x,1.5);
}
static double gg7(double x)
{
return pow(x,-0.5);
}
static double gg8(double x)
{
return 1/(1 + pow(x,4));
}
static double gg9(double x)
{
return 2/(2 + sin(10*pi*x));
}
static double gg10(double x)
{
return 1/(1+x);
}
static double gg11(double x)
{
return 1.0/(1 + pow(e,x));
}
static double gg12(double x)
{
return x/(pow(e,x)-1.0);
}
static double gg13(double x)
{
return sqrt(50)*pow(e,-50.0*pi*x*x);
}
static double gg14(double x)
{
return 25.0*pow(e,-25.0*x);
}
static double gg15(double x)
{
return 50.0/(pi*(2500.0*x*x+1));
}
static double gg16(double x)
{
return 50.0*pow((sin(50.0*pi*x)/(50.0*pi*x)),2);
}
static double gg17(double x)
{
return cos(cos(x)+3*sin(x)+2*cos(2*x)+3*cos(3*x));
}
static double gg18(double x)
{
return log10(x);
}
static double gg19(double x)
{
return 1/(1.005+x*x);
}
static double gg20(double x)
{
return 1/cosh(20.0*(x-1.0/5.0)) + 1/cosh(400.0*(x-2.0/5.0))
+ 1/cosh(8000.0*(x-3.0/5.0));
}
static double gg21(double x)
{
return 1.0/(1.0+(230.0*x-30.0)*(230.0*x-30.0));
}
static double gg22(double x)
{
if(x < 1)
{
return (x + 1.0);
}
else if(x >= 1 && x <= 3)
{
return (3.0 - x);
}
else
{
return 2.0;
}
}
};
numerical_integration_tester a;
}
// Numerical Integration method based on the adaptive Simpson method in
// Gander, W. and W. Gautschi, "Adaptive Quadrature – Revisited,"
// BIT, Vol. 40, 2000, pp. 84-101
// Test functions taken from Pedro Gonnet's dissertation at ETH:
// Adaptive Quadrature Re-Revisited
// http://e-collection.library.ethz.ch/eserv/eth:65/eth-65-02.pdf
#include <iostream>
#include <iomanip>
#include <stdint.h>
#include <dlib/matrix.h>
using namespace std;
using namespace dlib;
//***************************************************************//
//*Begin definitions of test functions //
//Initial Test Function
double f(double x)
{
return pow(x,0.5);
}
// The Lyness - Kaganove test functions from page 167 of Gonnet's thesis.
// lambda in [0,1], alpha in [-0.5,0], x in [0,1]
double LK1(double x, double lambda, double alpha)
{
return pow(abs(x-lambda),alpha);
}
// lambda in [0,1], alpha in [0,1], x in [0,1]
double LK2(double x, double lambda, double alpha)
{
if(x > lambda)
{
return 0;
}
else
{
return pow(e, alpha*x);
}
}
// lambda in [0,1], alpha in [0,4], x in [0,1]
double LK3(double x, double lambda, double alpha)
{
return pow(e,-alpha*abs(x-lambda));
}
// lambda in [1,2], alpha in [-6,-3], x in [1,2]
double LK4(double x, double lambda, double alpha)
{
return pow(10,alpha)/((x-lambda)*(x-lambda)+pow(10,alpha));
}
// lambda_i in [1,2], alpha in [-5,-3], x in [1,2]
double LK5(double x, double lambda1, double lambda2, double lambda3, double lambda4, double alpha)
{
return pow(10,alpha)/((x-lambda1)*(x-lambda1)+pow(10,alpha))
+ pow(10,alpha)/((x-lambda2)*(x-lambda2)+pow(10,alpha))
+ pow(10,alpha)/((x-lambda3)*(x-lambda3)+pow(10,alpha))
+ pow(10,alpha)/((x-lambda4)*(x-lambda4)+pow(10,alpha));
}
// lambda in [0,1], alpha in [1.8,2], x in [0,1]
double LK6(double x, double lambda, double alpha)
{
double beta = pow(10,alpha)/max(lambda*lambda,(1-lambda)*(1-lambda));
return 2*beta*cos(beta*(x-lambda)*(x-lambda));
}
// Test Battery from reference [33] and p. 168 of Gonnet's thesis.
// x in [0,1]
double GG1(double x)
{
return pow(e,x);
}
// x in [0,1]
double GG2(double x)
{
if(x > 0.3)
{
return 1.0;
}
else
{
return 0;
}
}
// x in [0,1]
double GG3(double x)
{
return pow(x,0.5);
}
// x in [0,1]
double GG4(double x)
{
return 22/25*cosh(x)-cos(x);
}
// x in [-1,1]
double GG5(double x)
{
return 1/(pow(x,4) + pow(x,2) + 0.9);
}
// x in [0,1]
double GG6(double x)
{
return pow(x,1.5);
}
// x in [0,1]
double GG7(double x)
{
return pow(x,-0.5);
}
// x in [0,1]
double GG8(double x)
{
return 1/(1 + pow(x,4));
}
// x in [0,1]
double GG9(double x)
{
return 2/(2 + sin(10*pi*x));
}
// x in [0,1]
double GG10(double x)
{
return 1/(1+x);
}
// x in [0,1]
double GG11(double x)
{
1/(1 + pow(e,x));
}
// x in [0,1]
double GG12(double x)
{
return x/(pow(e,x)-1);
}
// x in [0.1, 1]
double GG13(double x)
{
return sin(100.0*pi*x)/(pi*x);
}
// x in [0, 10]
double GG14(double x)
{
return sqrt(50)*pow(e,-50.0*pi*x*x);
}
// x in [0, 10]
double GG15(double x)
{
return 25.0*pow(e,-25.0*x);
}
// x in [0, 10]
double GG16(double x)
{
return 50.0/(pi*(2500.0*x*x+1));
}
// x in [0.01, 1]
double GG17(double x)
{
return 50.0*pow((sin(50.0*pi*x)/(50.0*pi*x)),2);
}
// x in [0, pi]
double GG18(double x)
{
return cos(cos(x)+3*sin(x)+2*cos(2*x)+3*cos(3*x));
}
// x in [0,1]
double GG19(double x)
{
return log10(x);
}
// x in [-1,1]
double GG20(double x)
{
return 1/(1.005+x*x);
}
// x in [0,1]
double GG21(double x)
{
return 1/cosh(20.0*(x-1/5)) + 1/cosh(400.0*(x-2/5)) + 1/cosh(8000.0*(x-3/5));
}
// x in [0,1]
double GG22(double x)
{
return 4*pi*pi*x*sin(20.0*pi*x)*cos(2*pi*x);
}
// x in [0,1]
double GG23(double x)
{
return 1/(1+(230*x-30)*(230*x-30));
}
// x in [0,3]
double GG24(double x)
{
return floor(pow(e,x));
}
// x in [0,5]
double GG25(double x)
{
if(x < 1)
{
return (x + 1);
}
else if(x >= 1 && x <= 3)
{
return 3 - x;
}
else
{
return 2;
}
}
// Returns double machine precision
// Taken from Wikipedia en.wikipedia.org/wiki/Machine_epsilon
template<typename float_t, typename int_t>
float_t machine_eps()
{
union
{
float_t f;
int_t i;
} one, one_plus, little, last_little;
one.f = 1.0;
little.f = 1.0;
last_little.f = little.f;
while(true)
{
one_plus.f = one.f;
one_plus.f += little.f;
if( one.i != one_plus.i )
{
last_little.f = little.f;
little.f /= 2.0;
}
else
{
return last_little.f;
}
}
}
// Main Integration Function.
// Supporting Integration Function
template <typename T, typename funct>
T AdaptSimpstp(const funct& f, T a, T b, T fa, T fm, T fb, T is)
{
T m = (a + b)/2;
T h = (b - a)/4;
T fml = f(a + h);
T fmr = f(b - h);
T i1 = h/1.5*(fa+4*fm+fb);
T i2 = h/3.0*(fa+4*(fml+fmr)+2*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))
{
cout << "INT ERR" << endl;
}
Q = i1;
}
else
{
Q = AdaptSimpstp(f, a, m, fa, fml, fm, is) + AdaptSimpstp(f,m,b,fm,fmr,fb,is);
}
return Q;
}
// Main integration function.
// f -- function to integrate,
// a -- left end point
// b -- right end point
// tol -- error tolerance
template <typename T, typename funct>
T AdaptSimp(const funct& f, T a, T b, T tol)
{
T eps = machine_eps<T, uint64_t>();
if(tol < eps)
{
tol = eps;
}
const T ba = b-a;
const T fa = f(a);
const T fb = f(b);
const T fm = f((a+b)/2);
T is =ba/8*(fa+fb+fm+ f(a + 0.9501*ba) + f(a + 0.2311*ba) + f(a + 0.6068*ba)
+ f(a + 0.4860*ba) + f(a + 0.8913*ba));
if(is == 0)
{
is = b-a;
}
is = is*tol/eps;
T tstvl = AdaptSimpstp(f, a, b, fa, fm, fb, is);
return tstvl;
}
// Examples
int main()
{
typedef double T;
T tol = 1e-10;
T a = 0;
T b = 5;
T tstvl2 = AdaptSimp(&f, a, b, tol);
cout << "Integral Value is: " << std::setprecision(18) << tstvl2 << endl;
return 0;
}
// The contents of this file are in the public domain. See
// LICENSE_FOR_EXAMPLE_PROGRAMS.txt
/*
This example demonstrates the usage of the numerical quadrature function
integrate_function_adapt_simpson. This function takes as input a single variable
function, the endpoints of a domain over which the function will be integrated, and a
tolerance parameter. It outputs an approximation of the integral of this function
over the specified domain. The algorithm is based on the adaptive Simpson method outlined in:
Numerical Integration method based on the adaptive Simpson method in
Gander, W. and W. Gautschi, "Adaptive Quadrature – Revisited,"
BIT, Vol. 40, 2000, pp. 84-101
*/
#include <iostream>
#include <stdint.h>
#include <dlib/matrix.h>
#include <dlib/numeric_constants.h>
#include <dlib/numerical_integration.h>
using namespace std;
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.
// 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 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
// [1e-10, 1e-8] region.
double tol = 1e-10;
// Here we instantiate a class which contains the numerical quadrature method.
// adapt_simp ad;
// Here we compute the integrals of the five functions defined above using the same
// tolerance level for each.
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) << 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;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment