pfqn_mmint2.h File Reference

McKenna-Mitra integral form of the normalizing constant of a repairman model (one queueing station, R classes, per-class think time), in its three MATLAB quadratures. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/api/qsys/qsys_quadrature.h"
#include "line/num/number.h"
#include "line/util/error.h"

Include dependency graph for pfqn_mmint2.h:

Go to the source code of this file.

Classes
struct	line::pfqn::MmintResult< T >
	Return value of the McKenna-Mitra quadratures, mirroring [G, lG]. More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
MmintResult< T >	line::pfqn::pfqn_mmint2 (const std::vector< T > &L, const std::vector< T > &N, const std::vector< T > &Z)
	Adaptive form (MATLAB pfqn_mmint2): Gauss-Kronrod on [0, 27.63] with absolute tolerance 1e-12.
template<class T>
MmintResult< T >	line::pfqn::pfqn_mmint2_gausslegendre (const std::vector< T > &L, const std::vector< T > &N, const std::vector< T > &Z, int m, std::size_t nodecap)
	Gauss-Legendre form on [0, 1e6] (MATLAB pfqn_mmint2_gausslegendre).
template<class T>
MmintResult< T >	line::pfqn::pfqn_mmint2_gausslegendre (const std::vector< T > &L, const std::vector< T > &N, const std::vector< T > &Z)
template<class T>
MmintResult< T >	line::pfqn::pfqn_mmint2_gausslaguerre (const std::vector< T > &L, const std::vector< T > &N, const std::vector< T > &Z, int m, std::size_t npts)
	Gauss-Laguerre form (MATLAB pfqn_mmint2_gausslaguerre).
template<class T>
MmintResult< T >	line::pfqn::pfqn_mmint2_gausslaguerre (const std::vector< T > &L, const std::vector< T > &N, const std::vector< T > &Z)

Detailed Description

McKenna-Mitra integral form of the normalizing constant of a repairman model (one queueing station, R classes, per-class think time), in its three MATLAB quadratures.

Templated port of matlab/src/api/pfqn/pfqn_mmint2.m, pfqn_mmint2_gausslegendre.m and pfqn_mmint2_gausslaguerre.m. All three evaluate

G = 1/prod_r N_r! * int_0^inf u^{m-1} e^{-u} prod_r (Z_r + L_r u)^{N_r} du

and differ only in the rule: adaptive Gauss-Kronrod on a truncated interval, a fixed Gauss-Legendre rule on [0, 1e6], or Gauss-Laguerre on [0, inf). The integrand is a polynomial times e^{-u}, so Gauss-Laguerre with enough nodes is exact up to rounding, which makes it the natural cross-check on the other two and on pfqn_ca.

NODES. MATLAB loads a Julia-generated table (gausslegendre-data.mat, gausslaguerre-data.mat). A table cannot be carried across arithmetics – it would pin every instantiation to the precision it was generated at – so the rules are regenerated in T by pfqn_asympt_common.h.

The Legendre form needs care. MATLAB's node count is n = max(300, min(tablesize, 2(sum N + m - 1) - 1)) and it takes the FIRST n entries of a 20000-point rule on [0, 1e6], which is not the same thing as a fresh n-point rule: the 20000-point prefix spans [0.0036, 557.8] and resolves the e^{-u} factor, whereas a genuine 300-point rule on [0, 1e6] has its first node at u = 13.7 and misses the mass entirely (it returns log G = -3.44 where the answer is 1.63). The port therefore generates the 20000-point rule and takes the same prefix, computing only the prefix since each Newton iteration is independent of the other nodes. nodecap is that table length and defaults to 20000, the length of MATLAB's gausslegendre-nodes.txt.

TRUNCATION of the adaptive form. MATLAB integrates over [0, -log(1 - (1 - 1e-12))] = [0, 27.63...], i.e. the 1 - 10^-order quantile of the unit exponential, and asks for AbsTol 1e-12. The port keeps both constants. That truncation is the dominant error for large populations, where the polynomial factor pushes mass well beyond the cutoff; the tests record where it starts to bite.

ARITHMETIC. Quadrature, so all three are gated on num_traits<T>::has_transcendental.

Definition in file pfqn_mmint2.h.