Matrix-Analytic Methods

QBD processes and matrix-analytic solutions.

The mam module implements matrix-analytic methods for analyzing queues with structured Markov chains, including QBD, MAP, RAP, and G/M/1 and M/G/1 type processes.

Key function categories:

  • MAP analysis: map_pie(), map_mean(), map_var(), map_scv(), map_skew()

  • MAP fitting: map2_fit(), mmpp2_fit(), aph_fit(), aph2_fit()

  • PH distributions: Phase-type analysis and fitting

  • Transformations: map_scale(), map_normalize(), map_timereverse()

  • QBD methods: qbd_R(), qbd_mapmap1(), qbd_raprap1()

  • Compression: compress_adaptive(), compress_spectral()

Matrix-Analytic Methods (MAM) for MAP/PH distributions.

Native Python implementations for analyzing Markovian Arrival Processes (MAPs), Phase-Type (PH) distributions, and related matrix-analytic methods.

Key algorithms:

map_piq: CTMC steady-state of MAP map_pie: Embedded DTMC steady-state map_lambda: Arrival rate computation map_mean, map_var, map_scv: Moment computations solver_mam_map_bmap_1: MAP/BMAP/1 queue solver using GI/M/1-type ETAQA solver_mam_bmap_map_1: BMAP/MAP/1 queue solver using M/G/1-type ETAQA

map_lambda(D0, D1=None)[source]

Compute the arrival rate (λ) of a MAP.

The arrival rate is λ = π * D1 * e where π is the steady-state and e is the column vector of ones.

Parameters:

  • D0 (ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (ndarray) – Visible transition matrix (optional)

Returns:

Arrival rate λ

Return type:

float

mmap_super_safe(mmap_list, maxorder=1000, method='default')[source]

Safely superpose multiple MMAPs into a single MMAP.

Combines multiple MMAPs while ensuring the order of the resulting MMAP does not exceed maxorder. If combining would exceed maxorder, alternative methods are used (fitting to exponential or simplified MAP).

Parameters:

  • mmap_list (List[Tuple[ndarray, List[ndarray]]]) – List of (D0, D_list) tuples, one per MMAP

  • maxorder (int) – Maximum allowed order (state space size) for result

  • method (str) – Combination method (“default” or “match”)

Returns:

Tuple of (D0_super, D_list_super) for the superposed MMAP

Return type:

Tuple[ndarray, List[ndarray]]

Algorithm:

  1. Sort MMAPs by squared coefficient of variation (SCV) of unmarked process

  2. Iteratively combine: start with simplest, add others in order

  3. If Kronecker product would exceed maxorder, use simpler approximations

mmap_compress(D0, D_list, method='default', target_order=None)[source]

Compress an MMAP using various approximation methods.

Reduces the state space of an MMAP while preserving key statistical properties (moments, inter-arrival time distribution, etc.).

The MMAP is [D0] + D_list with D_list = [D1, D11, …, D1C], i.e. the aggregate D1 followed by the C per-class marking matrices, so C = len(D_list) - 1. This is the layout of mmap_normalize and of every caller (npfqn_traffic_merge, mmap_fj), which pass mmap[1:].

Parameters:

  • D0 (ndarray) – Hidden transition matrix

  • D_list (List[ndarray]) – [D1, D11, …, D1C], the aggregate followed by the per-class marking matrices

  • method (str) –

    Compression method (same strings as mmap_compress.m). Options: - “default”, “mixture”, “mixture.order1”: order-1 mixture fitting - “mixture.order2”: order-2 mixture (mmap_mixture_fit_mmap) - “mamap2”: MAMAP(2,m) fit (mamap2m_fit_mmap) - “mamap2.fb”: MAMAP(2,m) forward/backward fit (mamap2m_fit_gamma_fb_mmap) - “m3pp.approx_cov”, “m3pp.approx_ag”, “m3pp.exact_delta”,

    ”m3pp.approx_delta”: M3PP(2,m) count fits (m3pp2m_fitc_theoretical)

    • ”exponential”: Fit to single-state Poisson (Python-only extra)

  • target_order (int | None) – If given, skip compression when the MMAP order is already at most this value. There is no implicit size guard: the order-1 mixture is a lossy approximation with a contract (see _compress_mixture_order1), and callers rely on it holding regardless of the input order.

Returns:

Tuple of (D0_compressed, D_list_compressed) in the same layout.

Return type:

Tuple[ndarray, List[ndarray]]

ldqbd(Q0, Q1, Q2, options=None)[source]

Solve a level-dependent QBD process.

Parameters:

  • Q0 (List[ndarray]) – List of upward transition matrices [Q0^(0), Q0^(1), …, Q0^(N-1)] Q0^(n) has shape (states_n, states_{n+1})

  • Q1 (List[ndarray]) – List of local transition matrices [Q1^(0), Q1^(1), …, Q1^(N)] Q1^(n) has shape (states_n, states_n)

  • Q2 (List[ndarray]) – List of downward transition matrices [Q2^(1), Q2^(2), …, Q2^(N)] Q2^(n) has shape (states_n, states_{n-1})

  • options (LdqbdOptions | None) – LdqbdOptions instance (uses defaults if None)

Returns:

LdqbdResult containing rate matrices R and stationary distribution pi

Return type:

LdqbdResult

Algorithm:

  1. Compute rate matrices backward using continued fraction recursion

  2. Compute stationary distribution forward using R matrices

class LdqbdOptions(epsilon=1e-10, max_iter=1000, verbose=False)[source]

Bases: object

Options for LDQBD solver.

Variables:

  • epsilon (float) – Convergence tolerance (default: 1e-10)

  • max_iter (int) – Maximum number of iterations (default: 1000)

  • verbose (bool) – Print debug information (default: False)

epsilon: float = 1e-10
max_iter: int = 1000
verbose: bool = False
ldqbd_mphc(D0, D1, alpha, c, arr_rate, sf=None, carry=False)[source]

Block-tridiagonal generator of an M/PH/c queue with level-dependent arrivals.

Parameters:

  • D0 – service sub-generator (p x p), phase changes without completion.

  • D1 – service completion block (p x p); D1 = (-D0 @ 1) @ alpha for a PH.

  • alpha – length-p vector a server starts each new job in.

  • c – number of identical servers (>= 1; capped at the top level).

  • arr_rate – length Nlev+1; arr_rate[n] is the arrival rate out of level n.

  • sf (Sequence[float] | None) – optional length-Nlev multiplier on the station’s TOTAL service rate at level n (load dependence). Each busy server then runs at sf[n-1]/min(n,c) of nominal, so sf[n-1] == min(n,c) reproduces the unscaled queue exactly.

  • carry (bool) – True for a correlated (non-renewal) MAP service, whose phase is carried from one service to the next in START order, idle periods included: each job’s service begins in the phase where the previous job’s service ended. That is how JMT and LDES sample a MAP service; alpha is then unused. See _ldqbd_mapc_carry().

Returns:

(Q0, Q1, Q2) with Q0[n] the upward block of level n (length Nlev), Q1[n] the local block of level n (length Nlev+1) and Q2[n] the downward block of level n+1 -> n (length Nlev, so Q2[n-1] is the block leaving level n, matching the argument order the native ldqbd() takes).

Return type:

Tuple[List[ndarray], List[ndarray], List[ndarray]]

Level sizes grow over the boundary levels 0..c and repeat above them, so the blocks joining differently sized neighbours are rectangular; ldqbd, ldqbd_R and ldqbd_pi all accept that heterogeneity.