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_infgen(D0, D1)[source]
Compute the infinitesimal generator of a MAP.
The generator Q = D0 + D1 represents the underlying CTMC.
- map_piq(D0, D1=None)[source]
Compute steady-state distribution of the underlying CTMC of a MAP.
Solves πQ = 0 where Q = D0 + D1 is the generator.
- map_prob(D0, D1=None)[source]
Stationary distribution of the underlying CTMC of a MAP (alias of
map_piq(), matching the MATLABmap_probname).
- map_pie(D0, D1=None)[source]
Compute equilibrium distribution of embedded DTMC.
The embedded DTMC has transition matrix P = (-D0)^{-1} * D1. Its steady-state is π_e = π * D1 / (π * D1 * e).
- 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.
- map_mean(D0, D1=None)[source]
Compute mean inter-arrival time of a MAP.
The mean is 1/λ where λ is the arrival rate.
- map_var(D0, D1=None)[source]
Compute variance of inter-arrival times of a MAP.
Var[X] = E[X²] - E[X]²
Uses map_moment for consistency with the JAR implementation.
- map_scv(D0, D1=None)[source]
Compute squared coefficient of variation (SCV) of a MAP.
SCV = Var[X] / E[X]² = (E[X²] - E[X]²) / E[X]²
- map_moment(D0, D1, k)[source]
Compute the k-th moment of inter-arrival time distribution.
E[X^k] = k! * π_e * (-D0)^{-k} * e
where π_e is the embedded DTMC steady-state.
- map_scale(D0, D1, new_mean)[source]
Rescale a MAP to a given MEAN inter-arrival time.
Port of map_scale.m: the third argument is the TARGET MEAN, not a multiplier. The rates are scaled by mean/new_mean, which leaves every normalized moment and every autocorrelation alone and moves only the first moment, and the result is passed through map_normalize (the feasibility repair) as the reference does.
THIS ARGUMENT USED TO BE A FACTOR here, and nowhere else: MATLAB, the JAR and the C++ port all take the new mean, and this module’s own private helper _map_scale in api/solvers/mam/mmap_fj.py already did too. The two conventions are each other’s reciprocal-ish (a factor c gives mean/c), so a call written for one and read by the other produces a MAP with the wrong rate and the right shape, which no moment check on the SCV would catch.
- map_normalize(D0, D1)[source]
Make a MAP feasible again: port of map_normalize.m.
Takes real parts, clips negative entries to zero and rebuilds D0’s diagonal so that (D0 + D1) e = 0. MATLAB, the JAR (Map_normalize.java) and the C++ port (map_transform.h) all do exactly this, and every caller here wants it: map_scale closes with it, kpcfit’s rescaling helper closes with it, and mmpp_rand uses it to turn two random matrices into a generator pair.
IT USED TO RESCALE THE MEAN TO ONE, which is a different operation altogether and left every one of those callers wrong in a way no moment check would show: mmpp_rand returned a D0 whose diagonal had never been repaired, so the pair was not a generator at all, and the two rescaling helpers had the mean they had just set pulled straight back to one. Use map_scale(D0, D1, 1.0) where unit mean is what is wanted.
- map_isfeasible(D0, D1, tolerance=1e-10)[source]
Check if (D0, D1) form a valid MAP.
A valid MAP requires: - D0 has non-positive diagonal and non-negative off-diagonal - D1 has non-negative elements - D0 + D1 is a valid generator (row sums = 0)
- map_exponential(mean)[source]
Create a MAP representation of an exponential distribution (MATLAB-style).
This is a MATLAB-compatible wrapper that takes mean instead of rate.
- Parameters:
mean (float) – Mean inter-arrival time (= 1/λ)
- Returns:
Tuple of (D0, D1) matrices representing Exp(1/mean)
- Return type:
Examples
>>> D0, D1 = map_exponential(2) # Poisson process with rate λ=0.5
- map_erlang(mean, k)[source]
Create a MAP representation of an Erlang-k distribution (MATLAB-style).
This is a MATLAB-compatible wrapper that takes mean as first argument.
- Parameters:
- Returns:
Tuple of (D0, D1) matrices representing Erlang-k with given mean
- Return type:
Examples
>>> D0, D1 = map_erlang(2, 3) # Erlang-3 with mean 2
- map_hyperexp(probs, means=None, p=0.99)[source]
Two-phase hyperexponential process as a MAP.
Two call forms, mirroring the two conventions in the codebase:
map_hyperexp(MEAN, SCV, p)with scalar arguments fits a two-phase hyperexponential to the given mean and squared coefficient of variation, selecting phase 1 with probabilityp(default 0.99). This is the MATLABmap_hyperexp.msignature. ReturnsNonewhen the request is outside the feasible set, as MATLAB returns{}: with a fixedpthe reachable SCV is bounded, e.g. SCV <= 3 at p = 0.5.map_hyperexp(probs, means)with array arguments builds the MAP of a hyperexponential with the given phase probabilities and phase means.
- Returns:
Tuple of (D0, D1) matrices, or None if the moment fit is infeasible
- map_gamma(D0, D1, limit=1000)[source]
Estimate the autocorrelation decay rate of a MAP.
Mirrors MATLAB map_gamma. For MAPs of order higher than 2 the ACF is not geometric, so the decay rate is obtained by fitting rho_k = RHO0*gamma^k in the least-squares sense, with RHO0 = (1 - 1/SCV)/2 held fixed.
This is NOT a Gamma-distribution constructor; use map_erlang for that. It is also distinct from map_gamma2, which returns the second largest eigenvalue of the embedded DTMC (the two agree for order 2 only).
- map_sumind(maps)[source]
Compute the sum of independent MAPs.
Creates a MAP representing the sum (concatenation) of independent random variables represented by the input MAPs.
- Parameters:
maps (list) – List of MAPs, each as (D0, D1) tuple
- Returns:
Tuple of (D0, D1) matrices representing the sum
- Return type:
Examples
>>> # Sum of exponential and Erlang-2 >>> MAP1 = map_exponential(1.0) >>> MAP2 = map_erlang(1.0, 2) >>> D0, D1 = map_sumind([MAP1, MAP2])
- map_cdf(D0, D1, points)[source]
Compute cumulative distribution function of inter-arrival times.
F(t) = 1 - π_e * exp(D0*t) * e
- Parameters:
- Returns:
CDF values at specified points
- Return type:
Examples
>>> map_cdf(D0, D1, 1.0) # Returns P(T <= 1) >>> map_cdf(D0, D1, [1.0, 5.0]) # Returns [P(T<=1), P(T<=5)]
- map_pdf(D0, D1, points)[source]
Compute probability density function of inter-arrival times.
f(t) = π_e * exp(D0*t) * (-D0) * e
- Parameters:
- Returns:
PDF values at specified points
- Return type:
Examples
>>> map_pdf(D0, D1, [0.5, 1.0, 2.0])
- map_acf(D0, D1, lags=1)[source]
Compute autocorrelation coefficients of inter-arrival times.
- Parameters:
- Returns:
Array of autocorrelation coefficients at specified lags
- Return type:
Examples
>>> map_acf(D0, D1) # lag-1 autocorrelation >>> map_acf(D0, D1, np.arange(1, 11)) # first 10 autocorrelations
- map_idc(D0, D1)[source]
Compute the asymptotic index of dispersion.
I = SCV * (1 + 2 * sum_{k=1}^{inf} rho_k)
where SCV is the squared coefficient of variation and rho_k is the lag-k autocorrelation coefficient.
- map_count_var(D0, D1, t)[source]
Compute variance of counting process at resolution t.
- Parameters:
- Returns:
Variance of arrivals in (0, t]
- Return type:
- Reference:
He and Neuts, “Markov chains with marked transitions”, 1998
- map_count_idc(D0, D1, t)[source]
Index of dispersion for counts (IDC) of a MAP at time point(s) t.
The IDC of the counting process A(t) associated to the MAP is I_a(t) = Var(A(t)) / E[A(t)], t > 0, i.e. the scaled variance-time curve. It interpolates between I_a(0+) = SCV of the interarrival time (renewal MAP) and the asymptotic value I_a(inf) = map_idc(MAP).
- Reference:
W. Whitt and W. You, “A Robust Queueing Network Analyzer Based on Indices of Dispersion”, eq. (1).
- map_varcount(D0, D1, tset)[source]
Compute variance of counting process (alternative implementation).
- map_count_moment(D0, D1, t, orders)[source]
Compute power moments of counts at resolution t.
Uses numerical differentiation of the moment generating function.
- map_mmpp2(mean, scv, skew=-1, acf1=-1)[source]
Fit an MMPP(2) as a MAP.
Matches the requested mean, SCV, skewness and lag-1 autocorrelation exactly. Raises ValueError when the request lies outside the MMPP(2) feasible set rather than returning a non-MAP.
- Parameters:
- Returns:
Tuple of (D0, D1) matrices
- Return type:
Examples
>>> D0, D1 = map_mmpp2(1, 2, -1, 0.2) # Minimal skewness, ACF=0.2 >>> D0, D1 = map_mmpp2(1, 2, -1, -1) # Minimal skewness, max ACF
- map_gamma2(D0, D1)[source]
Compute the second largest eigenvalue of embedded DTMC.
This is the autocorrelation decay rate.
- map_randn(k, mu=(1.0, 1.0), sigma=(0.5, 0.5))[source]
Generate a random MAP with normally distributed elements.
- map_max(D0_a, D1_a, D0_b, D1_b)[source]
Create MAP for max(X, Y) where X ~ MAP_A and Y ~ MAP_B.
The phase space is ordered as [(i,j) pairs, B-only phases, A-only phases]: in the first block both A and B are still running, in the second block A has already completed and B is awaited, in the third block B has completed and A is awaited. An arrival is recorded when the second of the two completes, i.e. only out of the last two blocks.
- map_pntiter(D0, D1, na, t, M=None)[source]
Compute probability of na arrivals in interval [0, t].
Uses iterative bisection method (Neuts and Li).
- map2_fit(e1, e2, e3=-1.0, g2=0.0)[source]
Fit a MAP(2) distribution to moments and autocorrelation.
Based on: A. Heindl, G. Horvath, K. Gross “Explicit inverse characterization of acyclic MAPs of second order”
- Parameters:
- Returns:
MAP: Tuple (D0, D1) if successful, None if failed
error_code: 0 for success, >0 for various errors
- Return type:
Tuple of (MAP, error_code) where
- Error codes:
0: Success 10: Mean out of bounds 20: Correlated exponential 30: h2 out of bounds 40: h3 out of bounds 51-54: g2 out of bounds
- map2_fit_idc(e1, e2, e3, idc)[source]
Fit a MAP(2) matching the first three moments and the index of dispersion.
A MAP(2) has a geometrically decaying autocorrelation, so its index of dispersion obeys I = SCV + (SCV-1)*g2/(1-g2), as reported in Section 5.2.2 of Casale, Mi, Cherkasova and Smirni, IEEE Trans. Soft. Eng. 37(5), 2011. The relation is inverted in closed form as g2 = (I-SCV)/(I-1) and the decay rate is passed to map2_fit. A third moment outside the feasible region is replaced by its lower limit (3/2)*e2^2/e1, the largest heavy-tail decay a MAP(2) admits.
The paper returns an exponential whenever SCV <= 1 or I < SCV. The rule does more than avoid an infeasible fit and must not be relaxed: a flow-equivalent server whose service is exponential and load dependent is exact for a product-form subnetwork by Norton’s theorem, whereas any MAP(2) fitted to the marginal inter-departure statistics is not, because the departure stream of the subnetwork is not independent of the rest of the model.
- Parameters:
- Returns:
Tuple (MAP, status) with status 0 all four descriptors matched, 1 exponential as burstiness is not representable, 2 third moment clamped, 3 third moment selected automatically, 4 fit failed and an exponential is returned
- Return type:
- map_joint(D0, D1, a, i)[source]
Compute joint moments of a MAP.
E[(X_{a1})^{i1} * (X_{a1+a2})^{i2} * …]
- map_issym(D0, D1=None)[source]
Check if MAP contains symbolic elements.
In Python, this always returns False as we use numpy arrays.
- map_feastol()[source]
Get the feasibility tolerance exponent for MAPs.
This is the exponent k of the toolbox feasibility tolerance 10^-k, so the tolerance itself is 10^-8. It is NOT the tolerance: map_feastol() == 8.
- Returns:
Tolerance exponent k, to be used as 10**(-map_feastol())
- Return type:
- map_largemap()[source]
Get threshold for “large” MAP where exact computation is expensive.
- Returns:
- Return type:
Order threshold (default
- mmap_normalize(D0, D_list)[source]
Fix MMAP feasibility by clamping negative values and closing the generator.
Port of matlab/src/api/mam/mmap_normalize.m (mirrored by the JAR jline.api.mam.Mmap_normalize). The MMAP is [D0] + D_list with D_list = [D1, D21, …, D2C], i.e. the aggregate D1 followed by the C per-class marking matrices, so C = len(D_list) - 1. This is the layout produced by _mmap_to_tuple and by the mmap_super / mmap_max callers, which pass result[1:].
The aggregate D1 is recomputed as the sum of the per-class matrices, and the diagonal of D0 closes each row of the generator D0 + D1:
D0[k,k] = -sum_{j!=k} D0[k,j] - sum_j D1[k,j]
- 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:
- Returns:
Tuple of (D0_super, D_list_super) for the superposed MMAP
- Return type:
- Algorithm:
Sort MMAPs by squared coefficient of variation (SCV) of unmarked process
Iteratively combine: start with simplest, add others in order
If Kronecker product would exceed maxorder, use simpler approximations
- mmap_mark(D0, D_list, prob)[source]
Reclassify arrivals in an MMAP according to a probability matrix.
Converts an MMAP with K classes to one with R classes based on a KxR probability matrix that describes how arrivals are reclassified.
- Parameters:
- Returns:
Tuple of (D0_new, D_list_new) where D_list_new has R matrices
- Return type:
Example
If prob = [[0.8, 0.2], [0.3, 0.7]], a 2-class MMAP becomes a 2-class MMAP where 80% of class-1 arrivals are now class-1, 20% become class-2, etc.
- mmap_scale(D0, D_list, M, max_iter=30)[source]
Scale the mean inter-arrival times of an MMAP.
Adjusts all matrices to achieve specified mean inter-arrival times.
- Parameters:
- Returns:
Tuple of (D0_scaled, D_list_scaled)
- Return type:
- Algorithm:
For single M: Scale all matrices uniformly by ratio = M_old / M For vector M: Use iterative coordinate descent to find scaling factors
- mmap_hide(D0, D_list, types)[source]
Hide (remove) specified arrival classes from an MMAP.
The hidden classes are set to zero matrices, effectively removing them from observation while maintaining the underlying stochastic process.
- Parameters:
- Returns:
Tuple of (D0_hidden, D_list_hidden) where specified classes are zero
- Return type:
Example
mmap_hide(D0, [D1, D2, D3], types=[1]) hides class 2 (index 1)
- aph2_assemble(l1, l2, p1)[source]
Build the APH(2) with the given parameters. Mirrors aph2_assemble.m.
- aph2_fitall(M1, M2, M3)[source]
Fit every APH(2) matching the given first three moments. Mirrors aph2_fitall.m.
Returns a list of 1 or 2 APH(2) fits. It is never empty: when no APH(2) solution is feasible the general APH fitter aph_fit(M1, M2, M3, 2) supplies a fallback, exactly as in the MATLAB and JAR sources.
- aph2_adjust(M1, M2, M3, method='simple')[source]
Find the APH(2)-feasible (M2, M3) closest to the given pair, holding M1 fixed. Mirrors the ‘simple’ branch of aph2_adjust.m [Telek and Heindl, 2002], which is the only branch aph2_fit uses.
- aph2_fit(M1, M2, M3)[source]
Fit an APH(2) to the first three moments, adjusting M2/M3 if they are infeasible. Mirrors aph2_fit.m.
M1 is always matched exactly: aph2_adjust never alters it, and the aph_fit fallback rescales its result to mean M1.
- aph2_fit_map(map_)[source]
Fit a MAP of arbitrary order by a canonical APH(2). Mirrors aph2_fit_map.m.
- aph2_fit_trace(T)[source]
Fit a trace of inter-arrival times by a canonical APH(2). Mirrors aph2_fit_trace.m: the first three raw sample moments drive aph2_fit.
- 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. Options: - “default”, “mixture”, “mixture.order1”: order-1 mixture fitting - “exponential”: Fit to single-state Poisson
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:
- mmap_exponential(lambda_rate, nclasses=None)[source]
Create an exponential MMAP (Poisson arrival process).
- mmap_isfeasible(mmap, tol=None)[source]
Check whether an MMAP is feasible up to the given tolerance.
Checks: - Elements are real (no imaginary parts) - D0 + D1 rows sum to zero (generator property) - Diagonal of D0 is negative - Off-diagonal of D0 is non-negative - All D1c matrices are non-negative - D1 = D11 + D12 + … + D1C
- mmap_count_lambda(mmap)[source]
Compute the arrival rate of the counting process for each class.
The rate for class k is λ_k = π * D_k * e where π is the steady-state of the underlying CTMC.
- mmap_pie(mmap)[source]
Compute the stationary probability of the DTMC embedded at restart instants after an arrival of each class.
- mmap_embedded(mmap)[source]
Compute the embedded DTMC transition matrices for each class.
For class k, the embedded matrix is Pc_k = (-D0)^{-1} * D_{2+k}
- mmap_sample(mmap, n_samples, pi=None, seed=None)[source]
Generate samples from a Marked MAP.
- Parameters:
- Returns:
T: Array of inter-arrival times
A: Array of class labels (1-indexed for MATLAB compatibility)
- Return type:
Tuple of (T, A) where
- mmap_sum(mmap, n)[source]
Create an MMAP representing the sum of n identical MMAPs (n-fold sum).
This creates an MMAP whose inter-arrival times are distributed as the sum of n inter-arrival times from the original MMAP.
- mmap_super(mmap_a, mmap_b=None, opt='default')[source]
Superpose two or more MMAPs.
- Parameters:
- Returns:
Superposed MMAP
- Return type:
- mmap_mixture(alpha, maps)[source]
Create a probabilistic mixture of MAPs.
Each MAP in the list is selected with probability alpha[i].
- mmap_max(mmap_a, mmap_b, k)[source]
Create MMAP for the maximum of arrivals from two synchronized MMAPs.
This models a synchronization queue where arrivals from both sources must be paired before release.
- mmap_maps(mmap)[source]
Extract K MAPs, one for each class of the MMAP[K] process.
For class k, the MAP is {D0 + D1 - D_{2+k}, D_{2+k}}
- mmap_count_moment(mmap, t, orders)[source]
Per-class power moments of the counting process at resolution t.
For class k the marginal counting process is that of the MAP
{D0 + sum_{j != k} D1_j, D1_k}: every arrival of another class is an unobserved transition of the marginal process. The moments of that MAP follow frommap_count_moment.
- mmap_count_idc(mmap, t)[source]
Compute the per-class Index of Dispersion of Counts at resolution t.
IDC = Var[N(t)] / E[N(t)]
- mmap_count_mcov(mmap, t)[source]
Compute the count covariance between each pair of classes at time scale t.
- mmap_idc(mmap)[source]
Compute the asymptotic Index of Dispersion of Counts for each class.
This is the limit of IDC(t) as t -> infinity.
- mmap_sigma(mmap)[source]
Compute one-step class transition probabilities.
p_{i,j} = P(C_k = j | C_{k-1} = i)
- mmap_sigma2(mmap)[source]
Compute two-step class transition probabilities.
p_{i,j,h} = P(C_k = h | C_{k-1} = j, C_{k-2} = i)
- mmap_forward_moment(mmap, orders, norm=True)[source]
Compute the theoretical forward moments of an MMAP.
Forward moments are E[X^k | previous arrival was class c].
- Parameters:
- Returns:
Matrix of shape (C, len(orders)) where element (c, k) is the order-k forward moment for class c
- Return type:
- mmap_backward_moment(mmap, orders, norm=True)[source]
Compute the theoretical backward moments of an MMAP.
Backward moments are E[X^k | next arrival will be class c].
- Parameters:
- Returns:
Matrix of shape (C, len(orders)) where element (c, k) is the order-k backward moment for class c
- Return type:
- mmap_cross_moment(mmap, k)[source]
Compute the k-th order moment of inter-arrival times between class pairs.
Computes E[X^k | previous = class i, next = class j] for all i, j.
- mmap_modulate(P, HT, MMAP)[source]
Modulate MMAPs in continuous time according to holding time distributions.
Creates an MMAP that switches between different MMAPs according to a continuous-time Markov modulated process.
- 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:
- Algorithm:
Compute rate matrices backward using continued fraction recursion
Compute stationary distribution forward using R matrices
- class LdqbdResult(R, pi, pi_cells=None)[source]
Bases:
objectResult of LDQBD solver.
- Variables:
R (List[numpy.ndarray]) – List of rate matrices R^(1), R^(2), …, R^(N)
pi (numpy.ndarray) – Stationary distribution vector [pi_0, pi_1, …, pi_N]
pi_cells (List[numpy.ndarray] | None) – Per-level stationary vectors, still resolved by phase; pi is their row sums. A caller that needs a phase-conditional quantity – the mean of something that depends on the phase and not only on the level – cannot recover it from the aggregated vector and must read these.
- class LdqbdOptions(epsilon=1e-10, max_iter=1000, verbose=False)[source]
Bases:
objectOptions for LDQBD solver.
- Variables:
- ldqbd_mphc(D0, D1, alpha, c, arr_rate, sf=None)[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) @ alphafor 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, sosf[n-1] == min(n,c)reproduces the unscaled queue exactly.
- Returns:
(Q0, Q1, Q2)withQ0[n]the upward block of level n (length Nlev),Q1[n]the local block of level n (length Nlev+1) andQ2[n]the downward block of level n+1 -> n (length Nlev, soQ2[n-1]is the block leaving level n, matching the argument order the nativeldqbd()takes).- Return type:
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.
- ph_multisets(p, k)[source]
Configurations of k identical servers over p service phases.
Rows are the compositions of k into p nonnegative parts:
M[r, i]is the number of the k busy servers sitting in phase i. There arecomb(k+p-1, p-1)of them, the multiset count of Asmussen and Moller (2001) – identical servers are exchangeable, so only the phase COUNTS carry information and the ordered space of sizep**kcollapses onto this one.The order is fixed and shared by every caller, so a configuration index means the same thing in each of them: the first part descends. k == 1 therefore yields the identity rows e_1 … e_p in phase order, which is what makes the c == 1 case of
ldqbd_mphc()coincide with plain phase indexing.
- class MAPBMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, R, mean_batch_size)[source]
Bases:
objectResult of MAP/BMAP/1 queue analysis.
- solver_mam_map_bmap_1(C0, C1, D)[source]
Solve a MAP/BMAP/1 queue using GI/M/1-type matrix-analytic methods.
The MAP is specified by matrices (C0, C1) where: - C0: transitions without arrivals - C1: transitions triggering arrivals
The BMAP for service is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without service completions - Dk: transitions triggering batch service of k customers (k >= 1)
The GI/M/1-type structure for MAP/BMAP/1 is:
B1 A0 0 0 ... B2 A1 A0 0 ... Q = B3 A2 A1 A0 ... ...
Where:
A0 = C1 otimes I_ms (MAP arrival, level +1) A1 = C0 otimes I_ms + I_ma otimes D0 (phase changes, level 0) A_{k+1} = I_ma otimes D_k (batch size k service, level -k)
- Parameters:
- Returns:
MAPBMAP1Result with performance metrics
- Return type:
- class BMAPMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, G, mean_batch_size)[source]
Bases:
objectResult of BMAP/MAP/1 queue analysis.
- solver_mam_bmap_map_1(D, S0, S1)[source]
Solve a BMAP/MAP/1 queue using M/G/1-type matrix-analytic methods.
The BMAP is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without arrivals - Dk: transitions triggering batch size k arrivals (k >= 1)
The MAP for service is specified by matrices (S0, S1) where: - S0: transitions without service completions - S1: transitions triggering service completions
The M/G/1-type structure for BMAP/MAP/1 is:
A0 = I_ma otimes S1 (service completion, level -1) A1 = D0 otimes I_ms + I_ma otimes S0 (phase changes, level 0) A_{k+1} = D_k otimes I_ms (batch arrival size k, level +k)
- Parameters:
- Returns:
BMAPMAP1Result with performance metrics
- Return type:
- class QbdRapResult(levelProb, QN, R, G, U, spr, pqueue, pi0)[source]
Bases:
objectResult of the equilibrium analysis of a QBD with RAP components.
- Variables:
levelProb (numpy.ndarray) – Marginal level probabilities, levels 0..numLevels
QN (float) – Mean queue length, computed exactly as pi0*R*inv(I-R)^2*e
R (numpy.ndarray) – Rate matrix R = A0*inv(-U)
G (numpy.ndarray) – Matrix G solving A0*G^2 + A1*G + A2 = 0
U (numpy.ndarray) – Matrix U = A1 + A0*G
spr (float) – Spectral radius Sp(R); positive recurrent iff Sp(R) < 1
pqueue (numpy.ndarray) – (numLevels+1) x m array whose n-th row is the level vector pi_n
pi0 (numpy.ndarray) – Level-0 vector pi_0, the boundary vector of Theorem 7
- qbd_R(B, L, F, iter_max=100000, tol=1e-12)[source]
Compute QBD rate matrix R using successive substitutions.
- Solves the matrix quadratic equation:
R^2 * A_{-1} + R * A_0 + A_1 = 0
where A_{-1} = B, A_0 = L, A_1 = F.
- Parameters:
- Returns:
Rate matrix R
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_R.m
- qbd_R_logred(B, L, F, iter_max=1000, tol=1e-14)[source]
Compute QBD rate matrix R using logarithmic reduction.
Uses the logarithmic reduction algorithm which has quadratic convergence compared to linear convergence of successive substitutions.
- Parameters:
- Returns:
Rate matrix R
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_R_logred.m Latouche & Ramaswami, Ch. 8
- qbd_rg(B, L, F, method='logred', iter_max=1000, tol=1e-14)[source]
Compute both R and G matrices for a QBD process.
- G is the minimal non-negative solution to:
A_1 * G^2 + A_0 * G + A_{-1} = 0
- R is the minimal non-negative solution to:
R^2 * A_{-1} + R * A_0 + A_1 = 0
- Parameters:
- Returns:
QBDResult with R, G, U, and eta (caudal characteristic)
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_rg.m
- qbd_blocks_mapmap1(D0_arr, D1_arr, D0_srv, D1_srv)[source]
Construct QBD blocks for a MAP/MAP/1 queue.
Builds the backward (B), local (L), and forward (F) transition blocks for the QBD representation of a MAP/MAP/1 queue.
- Parameters:
- Returns:
Tuple of (B, L, F) QBD blocks
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
- qbd_bmapbmap1(MAPa, pbatcha, MAPs)[source]
Compute QBD blocks for a BMAP/BMAP/1 queue.
Constructs the QBD (Quasi-Birth-Death) transition blocks for a BMAP/BMAP/1 queue with batch arrivals.
- Parameters:
- Returns:
A0: Local transition block A_1: Downward transition block A1_list: List of upward transition blocks for each batch size B0: Initial boundary local block B1_list: List of boundary upward blocks for each batch size
- Return type:
Tuple of (A0, A_1, A1_list, B0, B1_list) where
References
Original MATLAB: matlab/src/api/mam/qbd_bmapbmap1.m
- qbd_mapmap1(MAPa, MAPs, util=None)[source]
Analyze a MAP/MAP/1 queue using QBD methods.
Solves a MAP/MAP/1 queue using Quasi-Birth-Death process methods, computing throughput, queue length, utilization, and other metrics.
- Parameters:
- Returns:
Tuple of (XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs_scaled) where:
XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic (spectral radius of R) G: Rate matrix G A_1: Downward transition block A0: Local transition block A1: Upward transition block U: Matrix U MAPs_scaled: Scaled service process
- Return type:
Tuple[float, float, float, ndarray, ndarray, float | None, ndarray, ndarray, ndarray, ndarray, ndarray, Tuple[ndarray, ndarray]]
References
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
- qbd_rap(A0, A1, A2, B0=None, B1=None, numLevels=20)[source]
Equilibrium analysis of a Quasi-Birth-and-Death process with Rational Arrival Process (RAP) components.
The process is specified directly by its level-independent blocks (A0,A1,A2) and its boundary blocks (B0,B1), where A0 drives level increases, A2 drives level decreases and A1 the within-level evolution. Unlike a Markovian QBD the blocks need not be nonnegative: they are only required to be conservative, (A0+A1+A2)*e = 0, and to define a genuine RAP through the prediction-process interpretation. This makes qbd_rap strictly more general than qbd_raprap1, which builds a product-space QBD from two INDEPENDENT RAPs; here the arrival process and the sequence of service times may be driven from a shared phase space and therefore be cross-correlated.
This is the block-level core of the RAP QBD family. qbd_raprap1 is the thin wrapper over it that builds the product-space blocks of two independent RAPs; callers with a coupled model must call qbd_rap directly because no product form exists to factor out.
Algorithm (Theorem 7 of the reference):
1. Solve A0*G^2 + A1*G + A2 = 0 for G. 2. U = A1 + A0*G. 3. R = A0*inv(-U). 4. Find the row vector pihat0 with pihat0*(B1 + R*A2) = 0, pihat0*e = 1. 5. pi0 = K*pihat0 with K chosen so that pi0*inv(I-R)*e = 1. 6. pi_n = pi0*R^n, and the marginal level probability is pi_n*e.
The process is positive recurrent iff Sp(R) < 1 and step 4 has a solution.
Computation of G: the blocks are not nonnegative, so the probabilistic iterations used for Markovian QBDs (logarithmic reduction, cyclic reduction) carry no convergence guarantee here, and the paper explicitly leaves the general case open (“The issue of justifying algorithms for the evaluation of the matrix G for such processes has not been undertaken”, Section 6). See qbd_rap_g: the rank-one closed form is used when it applies, otherwise functional iteration followed by Newton’s method, and an unconverged G is never returned.
References
N. G. Bean and B. F. Nielsen, “Quasi-Birth-and-Death Processes with Rational Arrival Process Components”, Stochastic Models, 26(3), 2010, pp. 309-334 (DTU technical report IMM-2007-20). The argument rests on the prediction-process interpretation of a RAP due to Asmussen and Bladt, which is what allows a QBD argument to be carried over to matrices that are not nonnegative; the same prediction process underlies the conditional-vector RAP sampler in RAP.sample.
Original MATLAB: matlab/src/api/mam/qbd_rap.m
- Parameters:
- Returns:
QbdRapResult with levelProb, QN, R, G, U, spr, pqueue and pi0
- Raises:
ValueError – if the blocks are not conservative, if the process is not positive recurrent, or if G cannot be computed
- Return type:
- qbd_rap_g(A0, A1, A2, block_scale)[source]
Solve A0*G^2 + A1*G + A2 = 0 for the matrix G.
Uses the exact rank-one closed form when A2 has rank one, and otherwise natural functional iteration as a warm start followed by Newton’s method on the Sylvester-form Jacobian. Never returns an unconverged iterate.
- qbd_raprap1(RAPa, RAPs, util=None)[source]
Analyze a RAP/RAP/1 queue using QBD methods.
Solves a RAP/RAP/1 queue (Rational Arrival Process) using QBD methods, computing throughput, queue length, utilization, and other metrics.
References
N. G. Bean and B. F. Nielsen, “Quasi-Birth-and-Death Processes with Rational Arrival Process Components”, Stochastic Models, 26(3), 2010, pp. 309-334. The analysis rests on the prediction-process interpretation of a RAP due to Asmussen and Bladt, which is what allows a QBD argument to be carried over to matrices that are not nonnegative. The same prediction process underlies the conditional-vector RAP sampler in RAP.sample.
- Parameters:
- Returns:
XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic G: Rate matrix G B: Backward transition block L: Local transition block F: Forward transition block
- Return type:
Tuple of (XN, QN, UN, pqueue, R, eta, G, B, L, F) where
References
Original MATLAB: matlab/src/api/mam/qbd_raprap1.m
- qbd_setupdelayoff(lambda_val, mu, alpharate, alphascv, betarate, betascv)[source]
Analyze queue with setup delay and turn-off phases.
Performs queue-length analysis for a queueing system with setup delay (warm-up) and turn-off periods using QBD methods.
The system operates as follows: 1. When empty and job arrives, server enters setup phase 2. After setup, server becomes active and serves jobs 3. When queue empties, server enters turn-off phase 4. After turn-off, server becomes idle
- Parameters:
- Returns:
Average queue length QN
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_setupdelayoff.m
- qbd_depproc_etaqa(MAPa, MAPs, n)[source]
Construct MAP departure process for MAP/MAP/1-FCFS via ETAQA truncation.
Builds a Markovian Arrival Process (MAP) approximation of the departure process from a single-server FCFS queue with MAP arrival process and MAP service process. The QBD is truncated at level n using ETAQA boundary corrections.
- Parameters:
MAPa – Arrival MAP as list [D0, D1] of numpy arrays.
MAPs – Service MAP as list [D0, D1] of numpy arrays.
n – Truncation level (number of QBD levels beyond level 0).
- Returns:
List [D0, D1] representing the departure MAP, where D0 and D1 are numpy arrays of size (n+1)*lvlsz x (n+1)*lvlsz.
- qbd_depproc_etaqa_ps(MAPa, MAPs, n)[source]
Construct MAP departure process for MAP/MAP/1-PS via ETAQA truncation.
Same as qbd_depproc_etaqa but for Processor Sharing (PS) discipline. Under PS, when j customers are present, each is served at rate mu/j. This leads to rate-dependent service completion splitting in D0 and D1.
- Parameters:
MAPa – Arrival MAP as list [D0, D1] of numpy arrays.
MAPs – Service MAP as list [D0, D1] of numpy arrays.
n – Truncation level (number of QBD levels beyond level 0).
- Returns:
List [D0, D1] representing the departure MAP, where D0 and D1 are numpy arrays of size (n+1)*lvlsz x (n+1)*lvlsz.
- qbd_depproc_jointmom(MAPa, MAPs, iset)[source]
Compute joint moments E[X_0^i * X_1^j] of consecutive inter-departure times.
Uses matrix-analytic methods on the QBD representation of a MAP/MAP/1 queue to compute joint factorial moments of pairs of consecutive inter-departure times.
- Parameters:
MAPa – Arrival MAP as list [D0, D1] of numpy arrays.
MAPs – Service MAP as list [D0, D1] of numpy arrays.
iset – Array of shape (K, 2) where each row [i, j] specifies the moment orders for E[X_0^i * X_1^j].
- Returns:
1D numpy array of length K with the computed joint moments.
- map_ccdf_derivative(MAP, i)[source]
Compute derivative at zero of a MAP’s complementary CDF.
Calculates the i-th derivative at t=0 of the complementary cumulative distribution function (CCDF) of the inter-arrival time distribution.
Formula: ν_i = π_e * D_0^i * e
where π_e is the embedded stationary vector and e is the column vector of ones.
- Parameters:
- Returns:
Value of the i-th derivative at zero
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_ccdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”
- map_jointpdf_derivative(MAP, iset)[source]
Compute partial derivative at zero of a MAP’s joint PDF.
Calculates the partial derivative at t=0 of the joint probability density function of consecutive inter-arrival times.
For index set {i_1, i_2, …, i_k}: γ = π_e * D_0^{i_1} * D_1 * D_0^{i_2} * D_1 * … * D_0^{i_k} * D_1 * e
- Parameters:
- Returns:
Value of the partial derivative at zero
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_jointpdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”
- map_factorial_moment(MAP, k)[source]
Compute the k-th factorial moment of a MAP.
The k-th factorial moment is computed using derivatives of the inter-arrival time distribution.
- Parameters:
- Returns:
k-th factorial moment
- Return type:
References
Based on MAP moment formulas from matrix-analytic methods
- map_joint_moment(MAP, k, l)[source]
Compute the (k,l)-th joint moment of consecutive inter-arrival times.
E[X_n^k * X_{n+1}^l] for a MAP with inter-arrival times X_n.
Uses the embedded transition kernel P = (-D0)^{-1} D1 between the two resolvents. Sanity check: for a Poisson process the interarrivals are independent, so map_joint_moment(MAP, 1, 1) == map_mean(MAP) ** 2.
- Parameters:
- Returns:
Joint moment E[X_n^k * X_{n+1}^l]
- Return type:
References
Based on joint moment formulas for MAPs
- map_m1ps_cdf_respt(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]
Compute complementary sojourn time CDF for MAP/M/1-PS queue.
- Based on:
Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.
- Parameters:
- Returns:
Complementary CDF values W_bar(x) = Pr[W > x] at each point in x.
- Return type:
Example
>>> # Exponential arrivals (Poisson process with rate lambda=0.5) >>> C = np.array([[-0.5]]) >>> D = np.array([[0.5]]) >>> mu = 1.0 >>> x = np.array([0.0, 0.5, 1.0, 2.0, 5.0]) >>> cdf = map_m1ps_cdf_respt(C, D, mu, x)
- map_compute_R(C, D, mu)[source]
Compute rate matrix R for MAP/M/1 queue.
- Computes the minimal nonnegative solution of the matrix equation:
D + R(C - mu*I) + mu*R^2 = 0
This matrix is used in the analysis of MAP/M/1 queues based on quasi-birth-death processes.
- Parameters:
- Returns:
M x M rate matrix (minimal nonnegative solution)
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_compute_R.m Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.
- map_m1ps_h_recursive(C, D, mu, N, K)[source]
Recursive computation of h_{n,k} coefficients for MAP/M/1-PS.
The h_{n,k} vectors satisfy the recursion (Theorem 1 in paper):
h_{n,0} = e (vector of ones), for n = 0, 1, ... h_{n,k+1} = 1/(theta+mu) * [n*mu/(n+1) * h_{n-1,k} + (theta*I + C) * h_{n,k} + D * h_{n+1,k}] where h_{-1,k} = 0 for all k
- Parameters:
C (ndarray | list) – M x M matrix governing MAP transitions without arrivals
D (ndarray | list) – M x M matrix governing MAP transitions with arrivals
mu (float) – Service rate (scalar)
N (int) – Maximum value of n to compute (determines rows)
K (int) – Maximum value of k to compute (determines columns)
- Returns:
2D list h[n][k] of (M,1) arrays containing h_{n,k}
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_m1ps_h_recursive.m Masuyama, H., & Takine, T. (2003).
- map_m1ps_sojourn(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]
Compute sojourn time distribution in MAP/M/1-PS queue.
Alias for map_m1ps_cdf_respt.
- Parameters:
- Returns:
Complementary CDF values W_bar(x) = Pr[W > x]
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_m1ps_sojourn.m
- mmdp_isfeasible(Q, R, tol=1e-10)[source]
Check if (Q, R) defines a valid MMDP.
Requirements: - Q must be a valid generator (square, row sums = 0, proper signs) - R must be diagonal with non-negative entries - Q and R must have compatible dimensions
- Parameters:
- Returns:
True if (Q, R) defines a valid MMDP, False otherwise
- Return type:
Examples
>>> Q = np.array([[-0.5, 0.5], [0.3, -0.3]]) >>> R = np.array([[2.0, 0], [0, 5.0]]) >>> mmdp_isfeasible(Q, R) True
- hyperexp_fit_longtail(ccdf, k=None, c1=None, b=1.5, decade=4.0, points=None)[source]
Fit a hyperexponential to a long-tail distribution, recursively over time scales.
WHY MOMENTS ARE THE WRONG HANDLE. A Pareto law with tail index below 2 has infinite variance, so no two- or three-moment fit exists at all; and even when the moments are finite, matching them says nothing about the several ORDERS OF MAGNITUDE of time scale over which a long-tail distribution actually acts. This procedure matches the CCDF ITSELF at points spread across those decades.
THE RECURSION. Order the components so that
lambda_1 < ... < lambda_k. In the far tail only the slowest component survives, so(p_1, lambda_1)can be fitted there alone, from the ccdf atc_1andb c_1:lambda_1 = ln(F^c(c_1)/F^c(b c_1)) / ((b-1)c_1), p_1 = F^c(c_1) exp(lambda_1 c_1).
Subtract that component from the ccdf and repeat one decade lower, and so on (eqs. 4.6-4.11). The last component takes whatever probability is left,
p_k = 1 - sum_{j<k} p_j, and its rate follows from the ccdf atc_k(eqs. 4.12-4.14). This is Prony’s method applied to a ccdf.- Parameters:
ccdf (Callable[[float], float]) – F^c(t) = P(X > t) of the distribution to approximate
k (int | None) – number of exponential components;
Nonetakes one per decade between the 0.9 quantile and the 1e-6 quantile, which is the range the spacingdecadecan actually coverc1 (float) – the largest fitting argument; defaults to the point where the ccdf falls below 1e-6, which puts the slowest component in the real tail
b (float) – the within-scale spacing, 1 < b < c_i/c_{i+1}. The default pair (b, decade) = (1.5, 4) is not the paper’s illustrative (2, 10): the algorithm is exact AT the fitting arguments and free between them, and measured on a Weibull(0.3) the tighter grid cuts the worst between-point error from about 54% to 12%, at the cost of more components. Pass (2, 10) for the paper’s own figures
decade (float) – the ratio between successive fitting arguments,
c_i = c_1 decade^-(i-1); recomputed automatically when c1 is None so that the k arguments span from the 0.9 quantile to the 1e-6 onepoints (Sequence[float] | None) – explicit decreasing fitting arguments, overriding c1 and decade
- Returns:
Dict with
p(the mixing probabilities),lambda(the rates),points(the c_i used),meanof the fitted law,targetMeanof the original one (integrated over the covered range),coverage(the interval the fit is constrained on) andmaxRelErrorof the fitted ccdf at the fitting arguments. The last component matches only atc_k, its weight being fixed by the total probability, so the error atb c_kis not zero by construction.- Return type:
References
A. Feldmann, W. Whitt (1998). Fitting mixtures of exponentials to long-tail distributions to analyze network performance models. Performance Evaluation 31, 245-279.