Choudhury-Leung-Whitt normalization constant by numerical inversion of the generating function (JACM 42(5):935-970, 1995), and its limited load-dependent extension through the per-center transforms of Bertozzi and McKenna (SIAM Review 35(2):239-268, 1993). More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#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_clw.h:

Classes
struct	line::pfqn::ClwResult< T >
	Return value of pfqn_clw and pfqn_clw_lld, mirroring [G, lG]. More...
struct	line::pfqn::ClwOptions
	Optional lattice and aliasing parameters; empty means "use the CLW defaults". More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const std::vector< long > &m, const ClwOptions &opt)
	Choudhury-Leung-Whitt normalization constant by numerical inversion of the generating function (JACM 42(5):935-970, 1995), and its limited load-dependent extension through the per-center transforms of Bertozzi and McKenna (SIAM Review 35(2):239-268, 1993).
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z)
	Overload with unit multiplicities and the CLW default parameters.
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const std::vector< long > &m)
	Overload with the CLW default parameters.
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw_lld (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const Matrix< T > &mu, const ClwOptions &opt)
	Limited load-dependent form (matlab pfqn_clw_lld.m).
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw_lld (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const Matrix< T > &mu)
	Overload with the CLW default parameters.
template<class T>
ClwResult< T >	line::pfqn::pfqn_clw_lld (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z)
	Overload with all queues load independent.

Detailed Description

Choudhury-Leung-Whitt normalization constant by numerical inversion of the generating function (JACM 42(5):935-970, 1995), and its limited load-dependent extension through the per-center transforms of Bertozzi and McKenna (SIAM Review 35(2):239-268, 1993).

Templated port of matlab/src/api/pfqn/pfqn_clw.m and pfqn_clw_lld.m.

The generating function of g(K) is (CLW eq. 4.5)

G(z) = exp( sum_j rho_{j0} z_j ) / prod_i ( 1 - sum_j rho_{ji} z_j )^{m_i}

and, with limited load-dependent stations (Bertozzi-McKenna 2.17/2.23),

G(z) = exp( sum_j rho_{j0} z_j ) prod_i F_i( sum_j rho_{ji} z_j ), F_i(x) = [ c_i + sum_{n=1}^{l_i-1} (c_i - S_i(n)) / prod_{k<=n} S_i(k) x^n ] / (c_i - x),

analytic except for a simple pole at x = c_i. g(K) is the coefficient of prod_j z_j^{K_j}, recovered by p NESTED one-dimensional lattice-Poisson inversions (eq. 2.3) on contours of radius r_j = 10^{-gamma_j/(2 l_j K_j)}, with the restrictive static scaling of eqs. 5.41-5.46 and log-domain recovery (eq. 7.1). pfqn_clw applies both of the paper's speed-ups: dimension reduction by decomposition (Sec. 3, Sec. 5.4), which inverts the subset D minimizing |D| + max_i |S_i(D)| (eq. 3.3) and then each connected component of the remainder separately, and Euler summation of the inner sums (Sec. 2.4, eq. 2.22), which replaces 2 l_j K_j contour points by 2 l_j (n+m+1) wherever K_j > n+m, refining m until |E(m,n) - E(m,n+1)| settles. pfqn_clw_lld keeps the plain nested inversion of cost prod_j 2 l_j K_j.

COMPLEX ARITHMETIC WITHOUT std::complex. The contour integrand is genuinely complex, and std::complex is specified only for float, double and long double; instantiating it on a Boost.Multiprecision number is unspecified behavior. detail::Cx<T> below is a two-field complex with the six operations this routine needs, so the Real backends get real high precision on the inversion rather than silently falling back to double.

Arithmetic: TRANSCENDENTAL, double and Real only. exp, log, atan2, sqrt and a fractional power all appear; the contour radius alone is 10^{-gamma/(2 l K)}, which is not in the field of the inputs. This is the one member of the normalizing-constant family that CANNOT be instantiated exactly, and that is intrinsic to inverting a generating function numerically, not an artifact of the port. Accuracy against the exact pfqn_ca on the models tested is ~4e-10 for p = 2, ~7e-7 for p = 3, matching the reference's own claim of about 1e-9 and confirming that the port reproduces the reference's aliasing rather than adding error of its own.

REFERENCE DEFECTS

pfqn_clw returns NaN when SOME chain has zero population. The contour count is 2 l_j K_j, so a chain with K_j = 0 makes the final division acc / (2 l_j K_j r_j^{K_j}) a 0/0. pfqn_clw_lld guards against exactly this by dropping the zero-population chains up front ("the coefficient of z_j^0 equals the pgf restricted to z_j = 0, so chain j is removed exactly"); pfqn_clw never received that guard. Reproduce with pfqn_clw([0.1 0.2; 0.3 0.05], [2 0], [1.0 0.5]), which returns NaN where the answer is the p = 1 constant. THIS PORT APPLIES THE GUARD to both routines, so pfqn_clw here returns the finite value; that is the only input on which the two disagree.
Dead code in both routines: alpha0_j = exp(-alpha_j rho_{j0}) is computed in the scaling loop and never read. The recovery (eq. 7.1) is written in terms of sum_j alpha_j rho_{j0}, i.e. -sum_j log alpha0_j, so the array is redundant rather than wrong. It is not carried here.
denom(denom <= 0) = eps in the scaling loop silently substitutes 2.2e-16 for a nonpositive deflated denominator, which turns a chain whose predecessors have already saturated a queue into an enormous effective intensity rather than reporting the saturation. Reproduced, because the scaling only has to keep the contour inside the disc of analyticity and the recovery divides the choice back out, so the result is unaffected; but it is a silent branch, not a designed one.

Definition in file pfqn_clw.h.

Classes

Namespaces

Functions

Detailed Description