GIntegral class interface definition. More...

#include <GIntegral.hpp>

Inheritance diagram for GIntegral:

Public Member Functions
	GIntegral (void)
	Void constructor. More...

	GIntegral (GFunction *kernel)
	Function kernel constructor. More...

	GIntegral (const GIntegral &integral)
	Copy constructor. More...

virtual	~GIntegral (void)
	Destructor. More...

GIntegral &	operator= (const GIntegral &integral)
	Assignment operator. More...

void	clear (void)
	Clear integral. More...

GIntegral *	clone (void) const
	Clone integral. More...

std::string	classname (void) const
	Return class name. More...

void	max_iter (const int &iter)
	Set maximum number of iterations. More...

const int &	max_iter (void) const
	Return maximum number of iterations. More...

void	fixed_iter (const int &iter)
	Set fixed number of iterations. More...

const int &	fixed_iter (void) const
	Return fixed number of iterations. More...

void	eps (const double &eps)
	Set relative precision. More...

const double &	eps (void) const
	Get relative precision. More...

void	silent (const bool &silent)
	Set silence flag. More...

const bool &	silent (void) const
	Get silence flag. More...

const int &	iter (void) const
	Return number of iterations. More...

const int &	calls (void) const
	Get number of function calls. More...

const bool &	is_valid (void) const
	Signal if integration result is valid. More...

const std::string &	message (void) const
	Return integration status message. More...

void	kernel (GFunction *kernel)
	Set function kernel. More...

const GFunction *	kernel (void) const
	Get function kernel. More...

double	romberg (std::vector< double > bounds, const int &order=5)
	Perform Romberg integration. More...

double	romberg (const double &a, const double &b, const int &order=5)
	Perform Romberg integration. More...

double	trapzd (const double &a, const double &b, const int &n=1, double result=0.0)
	Perform Trapezoidal integration. More...

double	adaptive_simpson (const double &a, const double &b) const
	Adaptive Simpson's integration. More...

double	gauss_kronrod (const double &a, const double &b) const
	Gauss-Kronrod integration. More...

std::string	print (const GChatter &chatter=NORMAL) const
	Print integral information. More...

Public Member Functions inherited from GBase
virtual	~GBase (void)
	Destructor. More...

Protected Member Functions
void	init_members (void)
	Initialise class members. More...

void	copy_members (const GIntegral &integral)
	Copy class members. More...

void	free_members (void)
	Delete class members. More...

double	polint (double xa, double ya, int n, double x, double *dy)
	Perform Polynomial interpolation. More...

double	adaptive_simpson_aux (const double &a, const double &b, const double &eps, const double &S, const double &fa, const double &fb, const double &fc, const int &bottom) const
	Auxiliary function for adaptive Simpson's method. More...

double	rescale_error (double err, const double &result_abs, const double &result_asc) const
	Rescale errors for Gauss-Kronrod integration. More...

Protected Attributes
GFunction *	m_kernel
	Pointer to function kernel. More...

double	m_eps
	Requested relative integration precision. More...

int	m_max_iter
	Maximum number of iterations. More...

int	m_fix_iter
	Fixed number of iterations. More...

bool	m_silent
	Suppress integration warnings in console. More...

int	m_iter
	Number of iterations used. More...

int	m_calls
	Number of function calls used. More...

bool	m_isvalid
	Integration result valid (true=yes) More...

bool	m_has_abserr
	Has absolute integration error. More...

bool	m_has_relerr
	Has relative integration error. More...

double	m_abserr
	Absolute integration error. More...

double	m_relerr
	Absolute integration error. More...

std::string	m_message
	Status message (if result is invalid) More...

Detailed Description

GIntegral class interface definition.

This class allows to perform integration using various methods. The integrand is implemented by a derived class of GFunction.

Definition at line 46 of file GIntegral.hpp.

Constructor & Destructor Documentation

GIntegral::GIntegral ( void )

explicit

Void constructor.

Definition at line 231 of file GIntegral.cpp.

References init_members().

Referenced by clone().

GIntegral::GIntegral ( GFunction * kernel )

explicit

Function kernel constructor.

Parameters

[in] kernel Pointer to function kernel.

The function kernel constructor assigns the function kernel pointer in constructing the object.

Definition at line 249 of file GIntegral.cpp.

References init_members(), kernel(), and m_kernel.

GIntegral::GIntegral ( const GIntegral & integral )

Copy constructor.

Parameters

[in] integral Integral.

Definition at line 267 of file GIntegral.cpp.

References copy_members(), and init_members().

GIntegral::~GIntegral ( void )

virtual

Destructor.

Definition at line 283 of file GIntegral.cpp.

References free_members().

Member Function Documentation

double GIntegral::adaptive_simpson	(	const double &	a,
		const double &	b
	)		const

Adaptive Simpson's integration.

Parameters

[in]	a	Left integration boundary.
[in]	b	Right integration boundary.

Integrates the function using an adaptive Simpson's rule. The initial interval [a,b] is split into two sub-intervals [a,c] and [c,b] for which the integral is computed using

\[ \frac{b-a}{6} f(a) + 4f(c) + f(b) \]

where \(c=(a+b)/2\) is the mid-point of interval [a,b]. Each sub-interval is then recursively divided into sub-interval and the process is repeated. Dividing of sub-intervals is stopped when the difference between subsequent intervals falls below the relative tolerance specified by eps(). The maximum recursion depth is set by the max_iter() method.

I almost do not dare to confess, but the code has been taken from http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method It's really pretty simple ...

< Mid-point of interval [a,b]

< Length of interval [a,b]

Definition at line 708 of file GIntegral.cpp.

References abs(), adaptive_simpson_aux(), GFunction::eval(), m_calls, m_eps, m_has_abserr, m_has_relerr, m_isvalid, m_iter, m_kernel, m_max_iter, m_message, m_silent, gammalib::str(), and gammalib::warning().

double GIntegral::adaptive_simpson_aux	(	const double &	a,
		const double &	b,
		const double &	eps,
		const double &	S,
		const double &	fa,
		const double &	fb,
		const double &	fc,
		const int &	bottom
	)		const

protected

Auxiliary function for adaptive Simpson's method.

Parameters

[in]	a	Left integration boundary.
[in]	b	Right integration boundary.
[in]	eps	Precision.
[in]	S	Integral of last computation.
[in]	fa	Function value at left integration boundary.
[in]	fb	Function value at right integration boundary.
[in]	fc	Function value at mid-point of interval [a,b]
[in]	bottom	Iteration counter (stop when 0)

Implements a recursive auxiliary method for the adative_simpson() integrator.

< Mid-point of interval [a,b]

< Length of interval [a,b]

< Mid-point of interval [a,c]

< Mid-point of interval [c,b]

Definition at line 1173 of file GIntegral.cpp.

References abs(), eps(), GFunction::eval(), m_calls, m_isvalid, m_iter, and m_kernel.

Referenced by adaptive_simpson().

const int & GIntegral::calls ( void ) const

inline

Get number of function calls.

Returns: Number of function calls.

Definition at line 233 of file GIntegral.hpp.

References m_calls.

Referenced by print(), and romberg().

std::string GIntegral::classname ( void ) const

inlinevirtual

Return class name.

Returns: String containing the class name ("GIntegral").

Implements GBase.

Definition at line 129 of file GIntegral.hpp.

void GIntegral::clear ( void )

virtual

Clear integral.

Implements GBase.

Definition at line 335 of file GIntegral.cpp.

References free_members(), and init_members().

GIntegral * GIntegral::clone ( void ) const

virtual

Clone integral.

Returns: Pointer to deep copy of integral.

Implements GBase.

Definition at line 353 of file GIntegral.cpp.

References GIntegral().

void GIntegral::copy_members ( const GIntegral & integral )

protected

Copy class members.

Parameters

[in] integral Integral.

Definition at line 1044 of file GIntegral.cpp.

References m_abserr, m_calls, m_eps, m_fix_iter, m_has_abserr, m_has_relerr, m_isvalid, m_iter, m_kernel, m_max_iter, m_message, m_relerr, and m_silent.

Referenced by GIntegral(), and operator=().

void GIntegral::eps ( const double & eps )

inline

Set relative precision.

Parameters

[in] eps Relative precision.

Definition at line 208 of file GIntegral.hpp.

References eps(), and m_eps.

Referenced by GModelSpectralMultiplicative::eflux(), GModelSpectralExponential::eflux(), GModelSpectralLogParabola::eflux(), GModelSpectralSmoothBrokenPlaw::eflux(), GModelSpectralMultiplicative::flux(), GModelSpectralExponential::flux(), GModelSpectralLogParabola::flux(), GModelSpectralSmoothBrokenPlaw::flux(), GCOMIaq::location_smearing(), GCTAEdisp2D::normalize_table(), and GCOMD2Response::update_response_vector().

const double & GIntegral::eps ( void ) const

inline

Get relative precision.

Returns: Relative precision.

Definition at line 221 of file GIntegral.hpp.

References m_eps.

Referenced by adaptive_simpson_aux(), eps(), and print().

void GIntegral::fixed_iter ( const int & iter )

inline

Set fixed number of iterations.

Parameters

[in] iter Fixed number of iterations.

If the fixed number of iterations is set, the integration algorithm will always performed the given number of iterations, irrespectively of the precision that is reached. This feature is relevant for computing numerical derivates from numerically integrated functions.

Definition at line 183 of file GIntegral.hpp.

References iter(), and m_fix_iter.

const int & GIntegral::fixed_iter ( void ) const

inline

Return fixed number of iterations.

Returns: Fixed number of iterations.

Definition at line 196 of file GIntegral.hpp.

References m_fix_iter.

Referenced by print().

void GIntegral::free_members ( void )

protected

Delete class members.

Definition at line 1071 of file GIntegral.cpp.

Referenced by clear(), operator=(), and ~GIntegral().

double GIntegral::gauss_kronrod	(	const double &	a,
		const double &	b
	)		const

Gauss-Kronrod integration.

Parameters

[in]	a	Left integration boundary.
[in]	b	Right integration boundary.

Definition at line 766 of file GIntegral.cpp.

References abs(), GFunction::eval(), gammalib::gkw10, gammalib::gkw21a, gammalib::gkw21b, gammalib::gkw43a, gammalib::gkw43b, gammalib::gkw87a, gammalib::gkw87b, gammalib::gkx1, gammalib::gkx2, gammalib::gkx3, gammalib::gkx4, m_abserr, m_calls, m_eps, m_has_abserr, m_has_relerr, m_isvalid, m_iter, m_kernel, m_message, m_relerr, m_silent, rescale_error(), gammalib::str(), and gammalib::warning().

void GIntegral::init_members ( void )

protected

Initialise class members.

Definition at line 1015 of file GIntegral.cpp.

References m_abserr, m_calls, m_eps, m_fix_iter, m_has_abserr, m_has_relerr, m_isvalid, m_iter, m_kernel, m_max_iter, m_message, m_relerr, and m_silent.

Referenced by clear(), GIntegral(), and operator=().

const bool & GIntegral::is_valid ( void ) const

inline

Signal if integration result is valid.

Returns: True is integration result is valid.

Definition at line 297 of file GIntegral.hpp.

References m_isvalid.

Referenced by print().

const int & GIntegral::iter ( void ) const

inline

Return number of iterations.

Returns: Number of iterations.

Definition at line 141 of file GIntegral.hpp.

References m_iter.

Referenced by fixed_iter(), max_iter(), and print().

void GIntegral::kernel ( GFunction * kernel )

inline

Set function kernel.

Parameters

[in] kernel Function kernel.

Sets the function kernel for which the integral should be determined.

Definition at line 272 of file GIntegral.hpp.

References kernel(), and m_kernel.

const GFunction * GIntegral::kernel ( void ) const

inline

Get function kernel.

Returns: Function kernel.

Definition at line 285 of file GIntegral.hpp.

References m_kernel.

Referenced by GIntegral(), and kernel().

void GIntegral::max_iter ( const int & iter )

inline

Set maximum number of iterations.

Parameters

[in] iter Maximum number of iterations.

Definition at line 153 of file GIntegral.hpp.

References iter(), and m_max_iter.

const int & GIntegral::max_iter ( void ) const

inline

Return maximum number of iterations.

Returns: Maximum number of iterations.

Definition at line 166 of file GIntegral.hpp.

References m_max_iter.

Referenced by print(), and romberg().

const std::string & GIntegral::message ( void ) const

inline

Return integration status message.

Returns: Integration status message.

Definition at line 309 of file GIntegral.hpp.

References m_message.

Referenced by print().

GIntegral & GIntegral::operator= ( const GIntegral & integral )

Assignment operator.

Parameters

[in] integral Integral.

Returns: Integral.

Definition at line 305 of file GIntegral.cpp.

References copy_members(), free_members(), and init_members().

double GIntegral::polint	(	double *	xa,
		double *	ya,
		int	n,
		double	x,
		double *	dy
	)

protected

Perform Polynomial interpolation.

Parameters

[in]	xa	Pointer to array of X values.
[in]	ya	Pointer to array of Y values.
[in]	n	Number of elements in arrays.
[in]	x	X value at which interpolations should be performed.
[out]	dy	Error estimate for interpolated values.

Given arrays xa[1,..,n] and ya[1,..,n], and given a value x, this method returns a value y, and an error estimate dy. If P(x) is the polynomial of degree n-1, then the returned value y=P(x).

Todo:

Implement exceptions instead of screen dump.

Use std::vector for xa and ya and start at 0

Definition at line 1094 of file GIntegral.cpp.

References abs(), G_POLINT, m_isvalid, m_message, m_silent, gammalib::str(), and gammalib::warning().

Referenced by romberg().

std::string GIntegral::print ( const GChatter & chatter = NORMAL ) const

virtual

Print integral information.

Parameters

[in] chatter Chattiness.

Returns: String containing integral information.

Implements GBase.

Definition at line 946 of file GIntegral.cpp.

References calls(), eps(), fixed_iter(), is_valid(), iter(), m_abserr, m_fix_iter, m_has_abserr, m_has_relerr, m_relerr, max_iter(), message(), gammalib::parformat(), SILENT, silent(), and gammalib::str().

Referenced by cta_irf_diffuse_kern_theta::eval(), and cta_nroi_diffuse_kern_theta::eval().

double GIntegral::rescale_error	(	double	err,
		const double &	result_abs,
		const double &	result_asc
	)		const

protected

Rescale errors for Gauss-Kronrod integration.

Parameters

[in]	err	Error estimate.
[in]	result_abs	???.
[in]	result_asc	???.

Returns: Rescaled error estimate.

Definition at line 1236 of file GIntegral.cpp.

References abs(), and pow().

Referenced by gauss_kronrod().

double GIntegral::romberg	(	std::vector< double >	bounds,
		const int &	order = `5`
	)

Perform Romberg integration.

Parameters

[in]	bounds	Integration boundaries.
[in]	order	Integration order (default: 5)

Exceptions

GException::invalid_argument Integration order incompatible with number of iterations.

Returns the integral of the integrand, computed over a number of intervals [a0,a1], [a1,a2], ... that are given as an unordered vector by the bounds argument.

Integration is performed by Romberg's method of order 2*order, where

order=1 is equivalent to the trapezoidal rule,
order=2 is equivalent to Simpson's rule, and
order=3 is equivalent to Boole's rule.

The number of iterations is limited by m_max_iter. m_eps specifies the requested fractional accuracy. By default it is set to 1e-6.

Definition at line 381 of file GIntegral.cpp.

References calls(), and m_calls.

Referenced by GCTAModelAeffBackground::aeff_integral(), GModelSpectralMultiplicative::eflux(), GModelSpectralExponential::eflux(), GModelSpectralLogParabola::eflux(), GModelSpectralSmoothBrokenPlaw::eflux(), GModelSky::eflux(), cta_psf_radial_kerns_delta::eval(), GCTAModelIrfBackground::npred_roi_kern_theta::eval(), GCTAModelAeffBackground::npred_roi_kern_theta::eval(), GCTAModelSpatial::npred_roi_kern_theta::eval(), GModelSpatial::circle_int_kern_rho::eval(), cta_irf_radial_kern_rho::eval(), GResponse::irf_radial_kern_theta::eval(), GResponse::irf_elliptical_kern_theta::eval(), cta_nroi_radial_kern_rho::eval(), cta_irf_elliptical_kern_rho::eval(), cta_nroi_elliptical_kern_rho::eval(), cta_irf_diffuse_kern_theta::eval(), cta_nroi_diffuse_kern_theta::eval(), cta_psf_radial_kern_rho::eval(), cta_psf_radial_kern_delta::eval(), cta_psf_elliptical_kern_rho::eval(), cta_psf_diffuse_kern_delta::eval(), GModelSpectralMultiplicative::flux(), GModelSpectralExponential::flux(), GModelSpectralLogParabola::flux(), GModelSpatial::flux(), GModelSpectralSmoothBrokenPlaw::flux(), GModelSky::flux(), GLATPsfV3::integrate_psf(), GCTAResponseIrf::irf_diffuse(), GResponse::irf_elliptical(), GCTAResponseIrf::irf_elliptical(), GResponse::irf_radial(), GCTAResponseIrf::irf_radial(), GCOMIaq::location_smearing(), GLATPsfV3::normalize_psf(), GCTAEdisp2D::normalize_table(), GCTAModelIrfBackground::npred(), GCTAModelRadialAcceptance::npred(), GCTAModelSpatial::npred(), GObservation::npred(), GObservation::npred_spec(), GCTAResponseIrf::npsf(), GCTAResponseIrf::nroi(), GCTAResponseIrf::nroi_diffuse(), GCTAResponseIrf::nroi_elliptical(), GCTAResponseIrf::nroi_radial(), GCTAModelRadialGauss::omega(), GCTAModelRadialPolynom::omega(), GCTAModelRadialProfile::omega(), GLATPsfV1::psf(), GCTAResponseCube::psf_diffuse(), GCTAResponseCube::psf_elliptical(), GCTAResponseCube::psf_radial(), GCTAAeffArf::remove_thetacut(), GCTABackgroundPerfTable::solidangle(), and GCOMD2Response::update_response_vector().

double GIntegral::romberg	(	const double &	a,
		const double &	b,
		const int &	order = `5`
	)

Perform Romberg integration.

Parameters

[in]	a	Left integration boundary.
[in]	b	Right integration boundary.
[in]	order	Integration order (default: 5)

Exceptions

GException::invalid_argument Integration order incompatible with number of iterations.

Returns the integral of the integrand from a to b. Integration is performed by Romberg's method of order 2*order, where

order=1 is equivalent to the trapezoidal rule,
order=2 is equivalent to Simpson's rule, and
order=3 is equivalent to Boole's rule.

The number of iterations is limited by m_max_iter. m_eps specifies the requested fractional accuracy. By default it is set to 1e-6.

Definition at line 426 of file GIntegral.cpp.

References abs(), G_ROMBERG, gammalib::is_infinite(), gammalib::is_notanumber(), m_abserr, m_calls, m_eps, m_fix_iter, m_has_abserr, m_has_relerr, m_isvalid, m_iter, m_max_iter, m_message, m_relerr, m_silent, max_iter(), polint(), gammalib::str(), trapzd(), and gammalib::warning().

void GIntegral::silent ( const bool & silent )

inline

Set silence flag.

Parameters

[in] silent Silence flag.

Definition at line 245 of file GIntegral.hpp.

References m_silent, and silent().

Referenced by GModelSpatial::flux(), GCOMIaq::location_smearing(), and GCOMD2Response::update_response_vector().

const bool & GIntegral::silent ( void ) const

inline

Get silence flag.

Returns: True is class is silent, false otherwise.

Definition at line 258 of file GIntegral.hpp.

References m_silent.

Referenced by print(), and silent().

double GIntegral::trapzd	(	const double &	a,
		const double &	b,
		const int &	n = `1`,
		double	result = `0.0`
	)

Perform Trapezoidal integration.

Parameters

[in]	a	Left integration boundary.
[in]	b	Right integration boundary.
[in]	n	Number of steps.
[in]	result	Result from a previous trapezoidal integration step.

Returns: Integration results.

Exceptions

GException::invalid_value Function kernel not set.

The method performs a trapezoidal integration of the function kernel for the integration interval [a,b].

If n = 1 the integral is computed using

\[ \int_a^b f(x) \, dx = \frac{1}{2} (b-a) (f(a) + f(b)) \]

For n > 1 the integral is computed using

\[ \int_a^b f(x) \, dx = \frac{1}{2} \left[{\tt result} + \frac{b-a}{2^{n-1}} \sum_{i=0}^{2^{n-1}-1} f \left( a + (0.5+i) \frac{b-a}{2^{n-1}} \right) \right] \]

where \({\tt result}\) is the integration result from a previous call to the method with n = n - 1.

Definition at line 585 of file GIntegral.cpp.

References GFunction::eval(), G_TRAPZD, m_calls, m_isvalid, m_kernel, m_message, m_silent, gammalib::str(), sum(), and gammalib::warning().

Referenced by GCTAResponseIrf::npsf(), and romberg().