infer_gibbs.h File Reference

Gibbs sampling demand estimator for a closed delay-plus-queue model. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <vector>
#include "line/api/pfqn/pfqn_bs.h"
#include "line/api/pfqn/pfqn_mc_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for infer_gibbs.h:

Go to the source code of this file.

Classes
struct	line::infer::GibbsTrace< T >
	One class's trace. More...
struct	line::infer::GibbsStateProbs< T >
	The empirical state distribution built from the traces. More...
struct	line::infer::GibbsOptions
	MATLAB's hard-coded budgets, exposed with their MATLAB values as defaults. More...
struct	line::infer::GibbsSlice< T >
	The deterministic content of one coordinate update. More...

Namespaces
namespace	line
namespace	line::infer

Functions
template<class T>
GibbsStateProbs< T >	line::infer::gibbs_analyse_data (const std::vector< GibbsTrace< T > > &data, std::size_t data_needed)
	Empirical state distribution of the replayed traces (MATLAB's analyseData).
template<class T>
GibbsSlice< T >	line::infer::gibbs_slice (const std::vector< T > &think_time, const std::vector< T > &theta, const Matrix< long > &testset, std::size_t index, const std::vector< T > &N, const T &logG_init, double interval, const T &range_size)
	The deterministic half of one coordinate update: the log normalizing constant along the grid and the normalized slice it implies.
template<class T>
std::vector< T >	line::infer::infer_gibbs (const std::vector< GibbsTrace< T > > &data, const T &nbCores, const GibbsOptions &opts, pfqn::McRng &rng)
	Estimated per-class mean demands.
template<class T>
std::vector< T >	line::infer::infer_gibbs (const std::vector< GibbsTrace< T > > &data, const T &nbCores, double tol, pfqn::McRng &rng)
	MATLAB's three-argument form, with its hard-coded budgets.

Variables
static const std::size_t	line::infer::GIBBS_NNODES = 2
	The two nodes of the model the estimator assumes: delay, then queue.

Detailed Description

Gibbs sampling demand estimator for a closed delay-plus-queue model.

Templated port of matlab/src/api/infer/infer_gibbs.m. No JAR counterpart.

The observed data are per-class arrival and departure traces at one station. The estimator replays them into a continuous-time record of the population vector, turns that record into an empirical distribution over states, samples a test set from it, and then draws the per-class demands from their conditional posteriors one coordinate at a time. The posterior of a single demand is a slice: on a grid of candidate values,

log p(theta_j) = (sum over the test set of the class-j queue count) log theta_j - n_test log G(theta),

i.e. the product-form likelihood with the normalizing constant carried along. G is never evaluated directly. It is propagated along the grid with the exact derivative identity of the closed product-form constant, d log G / d theta_j = Q_j / theta_j, integrated one grid step at a time with Q from Bard-Schweitzer AMVA (the 'TE' method of the MATLAB). That is why every grid step costs one pfqn_bs solve, warm started from the previous point, and why the walk starts at the current theta with the running log G the caller carries between coordinate updates.

WHAT IS AND IS NOT PORTED. MATLAB's 'MCI' branch and its pdf_slice helper are unreachable: alg is assigned 'TE' as a literal with no way in, and pdf_slice is DECLARED with eleven parameters and CALLED with ten from the one dead call site, so the MCI branch would raise on its first use. Dead, broken code is not ported. Everything the 'TE' path executes is.

The four sample budgets that MATLAB hard-codes (data_needed, the test set size, the chain length, and the 50-sample convergence block) are options here, defaulting to MATLAB's values. That is a strict superset: the default construction reproduces the MATLAB exactly, and a caller that wants a short chain no longer has to edit the source.

RANDOMNESS. Two places draw: the test set is sampled from the empirical state distribution by inversion, and each coordinate update is drawn from its normalized slice. Both take an explicit McRng, so a run is reproducible from its seed. MATLAB's Mersenne Twister stream is deliberately NOT reproduced – same algorithm, different stream – so the two implementations agree in distribution and on every deterministic intermediate, not sample by sample.

MATLAB's eps in the derivative denominators is the DOUBLE machine epsilon, a fixed constant of the formula rather than a property of the working arithmetic, so the port carries the literal 2^-52 into T instead of asking the arithmetic for its own epsilon. Substituting Real50's epsilon would change the first grid step, where theta = 0 and the constant is the entire denominator.

ARITHMETIC: logarithms throughout, plus a fixed-point AMVA iteration per grid point, so the routine is gated on transcendental arithmetic and registered for Double only. Real is NOT registered: log(0) at the first grid point is a defined -infinity in IEEE double and the slice weight it produces underflows to exactly zero, which is the behaviour the algorithm relies on; that is a property of the double arithmetic MATLAB runs in, and the port does not claim a high-precision instantiation it has not validated.

Definition in file infer_gibbs.h.