pfqn_cub.h File Reference

Normalizing constant by Grundmann-Moeller cubature over the simplex. More...

#include <cmath>
#include <cstddef>
#include <functional>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for pfqn_cub.h:

Go to the source code of this file.

Classes
struct	line::pfqn::CubResult< T >
	Return value of pfqn_cub, mirroring [Gn, lGn]. More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
std::vector< T >	line::pfqn::grnmol (const std::function< T(const std::vector< T > &)> &f, std::size_t n, int s, const T &tol)
	Grundmann-Moeller rule of degrees 1, 3, ..., 2s+1 over the n-simplex with vertices the columns of the identity (MATLAB's grnmol on V = eye(n,n+1)).
template<class T>
CubResult< T >	line::pfqn::pfqn_cub (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, int order, const T &atol)
	Normalizing constant by Grundmann-Moeller cubature over the simplex.
template<class T>
CubResult< T >	line::pfqn::pfqn_cub (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z)

Detailed Description

Normalizing constant by Grundmann-Moeller cubature over the simplex.

Templated port of matlab/src/api/pfqn/pfqn_cub.m together with its two local functions simplexquad and grnmol (Grundmann and Moller, SIAM J. Numer. Anal. 15 (1978) 282-290). With Z = 0 the constant is a single integral of prod_r (u' L(:,r))^{N_r} over the (M-1)-simplex; the degree-(2s+1) rule is EXACT once s >= ceil((sum N - 1)/2), because the integrand is then a polynomial of degree sum(N) <= 2s+1. With Z > 0 an outer integral over the McKenna-Mitra scale variable v is added on a uniform grid, which is where the method stops being exact.

The grnmol rule is exported, because it is the only working Grundmann-Moeller implementation in the reference tree; see the note on pfqn_grnmol in the report accompanying this port.

ARITHMETIC. The rule itself is a weighted sum of integrand values at rational barycentric points, so at Z = 0 and full degree the whole computation would be exact in a field – were it not for the exp(gammaln(1+sum N+M-1) - sum gammaln(1+N)) prefactor, which MATLAB forms in logarithms. Since that prefactor is a ratio of factorials it could be formed exactly, but the reference does not, and reproducing the reference is the contract; the routine is therefore gated on num_traits<T>::has_transcendental, and the exactness claim above is about the cubature, not about the returned scalar.

Definition in file pfqn_cub.h.