Markov Chain Utilities
CTMC and DTMC analysis tools.
The mc module provides general tools for analyzing Markov chains, including
steady-state and transient analysis for both continuous-time and discrete-time
Markov chains.
Key function categories:
Steady-state analysis:
ctmc_solve_reducible(),dtmc_solve_reducible(),ctmc_stochcomp()Transient analysis:
ctmc_transient(),ctmc_uniformization()Simulation:
ctmc_simulate(),ctmc_rand()Aggregation methods:
ctmc_courtois(),ctmc_takahashi(),ctmc_kms()State-space generation:
ctmc_ssg(),ctmc_ssg_reachability()Generator matrices:
ctmc_makeinfgen(),ctmc_multi()
Markov Chain Analysis (line_solver.api.mc)
Markov Chain analysis algorithms.
Native Python implementations for continuous-time and discrete-time Markov chain analysis.
- Key algorithms:
ctmc_solve: CTMC steady-state distribution ctmc_sens: Steady-state sensitivity to a scalar parameter ctmc_transient_sens: Transient sensitivity to a scalar parameter ctmc_transient: CTMC transient analysis ctmc_uniformization: Transient distribution by Jensen uniformization ctmc_randomization: Randomized (uniformized) DTMC P = I + Q/q ctmc_foxglynn: Fox-Glynn uniformization for transient analysis ctmc_fau: Fast adaptive uniformization for transient analysis ctmc_gmres: Restarted GMRES with ILUT preconditioning ctmc_bicgstab: BiCGSTAB with the same ILUT preconditioner ctmc_saddlepoint: Pr{N(t)=k} of a MAP counting process by saddlepoint ctmc_stochcomp: Stochastic complementation dtmc_solve: DTMC steady-state distribution
- ctmc_solve(Q, method=None)[source]
Solve for steady-state probabilities of a CTMC.
Computes the stationary distribution π by solving πQ = 0 with normalization constraint Σπ = 1.
Handles reducible CTMCs by decomposing into strongly connected components and solving each separately.
- Parameters:
- Returns:
Steady-state probability distribution (1D array)
- Return type:
- ctmc_solve_reducible(Q, pin=None)[source]
Solve reducible CTMCs by converting to DTMC via uniformization.
Port of MATLAB ctmc_solve_reducible.m, which is a thin delegate to dtmc_solve_reducible(ctmc_randomization(Q), pin, tol=1e-12). The whole convention for a reducible chain therefore lives in dtmc_solve_reducible: the limiting vector of a chain with several closed communicating classes is NOT unique, and it is resolved by starting uniformly over the SCCs and propagating through the lumped limiting matrix (so, when every SCC is closed, each class carries equal weight).
This formerly delegated to ctmc_solve, which just returns whichever null vector the linear solver happens to land on – for two closed classes that is all the mass on the first, disagreeing with MATLAB and the JAR.
- ctmc_solve_reducible_blkdecomp(Q, pin=None)[source]
Solve reducible CTMCs via direct block decomposition on the generator.
- Algorithm:
Decompose states into transient and recurrent classes via SCC
For transient states: solve sojourn * Q_tt = -p0_t for expected sojourn
Compute hitting probabilities: hit = sojourn * Q_ta + p0_r
For each recurrent class: solve pi_c * Q_cc = 0, scale by hitting prob
This avoids the randomization to DTMC used in ctmc_solve_reducible.
- ctmc_sens(Q, dQ, pi=None)[source]
Sensitivity of the steady-state distribution of a CTMC to a scalar parameter theta, given the generator Q, its derivative dQ = dQ/dtheta, and the steady-state vector pi.
Differentiating the balance equations pi*Q = 0 and pi*e = 1 with respect to theta gives the linear system
(dpi/dtheta) * Q = -pi * (dQ/dtheta), sum_i dpi_i/dtheta = 0,
i.e. Trivedi and Bobbio (2017), Eq. (9.81). The system has the same coefficient matrix as the steady-state solve itself, so obtaining a sensitivity costs one extra solve against a matrix that is already assembled. The normalization replaces one row of the singular Q’, exactly as in the steady-state solve.
- ctmc_transient_sens(Q, dQ, pi0=None, t0=None, t1=None)[source]
Sensitivity of the transient distribution of a CTMC to a scalar parameter theta, given the generator Q and its derivative dQ = dQ/dtheta.
Differentiating the forward equations d pi(t)/dt = pi(t) Q with respect to theta, and assuming the initial vector does not depend on theta, gives
d/dt (dpi(t)/dtheta) = (dpi(t)/dtheta) Q + pi(t) (dQ/dtheta), dpi(0)/dtheta = 0,
i.e. Trivedi and Bobbio (2017), Eq. (9.82). The state and its sensitivity are integrated as one augmented system of size 2n, since the sensitivity equation is driven by pi(t) and the two cannot be advanced separately.
- Parameters:
- Returns:
(dpi, pi, t) with dpi the sensitivity of the distribution at each time point (len(t), n), pi the distribution at each time point (len(t), n) and t the vector of time points.
- Return type:
References
Original MATLAB: matlab/src/api/mc/ctmc_transient_sens.m
- ctmc_makeinfgen(Q)[source]
Convert a matrix into a valid infinitesimal generator for a CTMC.
An infinitesimal generator has: - Row sums equal to zero - Non-positive diagonal elements - Non-negative off-diagonal elements
- Parameters:
Q – Candidate infinitesimal generator matrix
- Returns:
Valid infinitesimal generator matrix with corrected diagonal
- ctmc_transient(Q, initial_dist, time_points, method='expm', epsilon=1e-6, delta=1e-12)[source]
Compute transient probabilities of a CTMC.
Calculates time-dependent state probabilities π(t) for specified time points using matrix exponential methods.
- Parameters:
Q (ndarray) – Infinitesimal generator matrix
initial_dist (ndarray) – Initial probability distribution π(0)
time_points (float | ndarray) – Array of time points to evaluate, or single time value
method (str) – ‘expm’ for matrix exponential, ‘ode’ for ODE solver, ‘fau’ for fast adaptive uniformization (ctmc_fau) MARCHED over the grid
epsilon (float) – ‘fau’ only, total probability mass the whole grid may discard; it is divided by the number of steps, each step removing mass and none putting any back, so the accumulated defect stays below it
delta (float) – ‘fau’ only, occupancy below which a state is dropped
- Returns:
Transient probabilities at each time point. Shape: (len(time_points), n) if multiple times, (n,) if single time
- Return type:
- ctmc_timeaverage(pi0, Q, t, tol=1e-12, maxiter=-1)[source]
Time-averaged transient distribution of a CTMC over [0, t] via uniformization.
Companion of the endpoint pi0*exp(Q*t); additionally returns the time average
piTimeAvg = pi0 * (1/t) * int_0^t exp(Q*tau) d(tau)
as well as the endpoint piExit = pi0*exp(Q*t), both from the same Jensen uniformization series. Used by the SolverENV state-vector analyzer (deterministic-sojourn option). Mirrors matlab ctmc_timeaverage.m, including the horizon splitting that keeps exp(-q*t) from underflowing and the adaptive truncation depth (maxiter <= 0).
- Returns:
(piTimeAvg, piExit) as 1D arrays.
- ctmc_uniformization(pi0, Q, time_points, precision=1e-10)[source]
Transient distribution of a CTMC by uniformization (Jensen’s method).
Numerically stable and free of matrix exponentials: the CTMC is randomized into a DTMC and the distribution is the Poisson-weighted sum of its powers. Argument order follows MATLAB ctmc_uniformization(pi0, Q, t) and the kpctoolbox twin; for the randomized matrix itself see ctmc_randomization.
- Parameters:
- Returns:
shape (len(time_points), n), or (n,) when a single time point is given.
- Return type:
Transient probabilities
- ctmc_foxglynn(pi0, Q, t, tol=1e-12, maxiter=-1)[source]
Transient distribution of a CTMC by Fox-Glynn uniformization.
- Parameters:
pi0 (ndarray) – Initial probability distribution
Q (ndarray) – Infinitesimal generator matrix
t (float) – Transient analysis period boundary [0,t]
tol (float) – Poisson tail-mass truncation tolerance
maxiter (int) – Maximum truncation depth; pass a nonpositive value to let the Fox-Glynn right truncation point size it
- Returns:
Transient probability vector at time t
- Return type:
- ctmc_foxglynn_weights(lam, tol=1e-12, maxiter=-1, normalize=True)[source]
Fox-Glynn truncation window and Poisson weights.
Following Fox-Glynn, the weights are built by the two-sided recursion w[k-1] = w[k]*k/lam and w[k+1] = w[k]*lam/(k+1) anchored at the mode, so neither exp(-lam) nor lam^k/k! is ever evaluated and the overflow and underflow that limit the direct series cannot occur. Anchoring at w[mode] = 1 keeps the extreme weights near tol, far above the denormal threshold, making Fox and Glynn’s rescaling of the mode weight unnecessary here. The normalizing sum is accumulated in increasing order of magnitude.
With normalize left at its default the window is rescaled to sum to one, as Fox and Glynn prescribe, so the truncated tails are redistributed over the window. Passing normalize=False scales the anchor by the true mode probability through a log-gamma instead, returning the Poisson probabilities themselves so that 1 - sum(w) is the discarded tail rather than being absorbed; ctmc_fau needs that, its error being reported as missing mass rather than as a bound.
- Parameters:
lam (float) – Poisson rate, that is the uniformization constant times the horizon
tol (float) – Poisson tail-mass truncation tolerance
maxiter (int) – Cap on the right truncation point; nonpositive leaves it uncapped
normalize (bool) – Rescale the window to sum to one, rather than returning the true Poisson probabilities
- Returns:
Tuple of (left truncation point, right truncation point, weights)
- Return type:
- ctmc_fau(pi0, Q, t, epsilon=1e-6, delta=1e-12, maxsteps=-1)[source]
Transient distribution of a CTMC at time t by fast adaptive uniformization.
Ordinary uniformization fixes one rate q >= max_i |q_ii| over the whole state space and mixes the powers of P = I + Q/q against a Poisson(q*t) law, so its cost is set by the fastest state anywhere, including states that carry no probability at time t. Adaptive uniformization instead picks a rate per step from the states the iterate occupies,
Lambda_n >= max{|q_ii| : i in supp(u^(n))}, u^(n+1) = u^(n)(I + Q/Lambda_n),
which keeps every entry of u^(n+1) nonnegative. The subordinating process is then the pure birth process N(t) with rates Lambda_0, Lambda_1, … and
pi(t) = sum_{n>=0} P{N(t) = n} u^(n).
The fast variant drops an entry of u^(n) below delta rather than propagating it, so the support tracks the states of non-negligible occupancy instead of the reachable set. Nothing is renormalized anywhere, so the error is not estimated but measured: the birth index truncated at K, the Poisson window of the weight computation and the delta threshold each remove mass and none puts any back, whence
0 <= pi(t) - pit componentwise, and |pi(t) - pit|_1 = sum(pi0) - sum(pit) = info.error_bound.
The birth weights are computed exactly rather than quadratured, by uniformizing the bidiagonal birth generator; see _fau_weights. The sweep runs twice because b_n(t) needs the rates up to n, which are not known before the sweep ends, while u^(n) is needed after them, and storing every iterate would cost K times the support. Stopping is certified by stochastic domination of the birth epochs by an Erlang, so this method never takes more steps than uniformization at the largest rate it visited.
This is a transient method: it produces no stationary distribution.
- Parameters:
pi0 (ndarray) – Initial probability distribution
Q – Infinitesimal generator matrix, dense or scipy sparse (CSR is used as given, other sparse formats are converted)
t (float) – Time horizon, t >= 0
epsilon (float) – Birth-process truncation tolerance
delta (float) – Occupancy threshold below which a state is dropped
maxsteps (int) – Cap on birth steps; nonpositive for the default cap
- Returns:
Tuple of (defective distribution at time t, diagnostics)
- Return type:
- class CtmcFauInfo(steps, lambda_min, lambda_max, uniform_rate, weight_tail, weight_window, dropped_mass, error_bound, support_max, support_final, truncated, absorbed)[source]
Bases:
objectDiagnostics of a fast adaptive uniformization sweep.
- Variables:
steps (int) – number of birth steps K+1 actually taken
lambda_min (float) – smallest adaptive rate used
lambda_max (float) – largest adaptive rate used, the Lstar of the weights
uniform_rate (float) – max_i |q_ii|, the rate ordinary uniformization would use
weight_tail (float) – mass reaching the overflow index, that is P{N(t) > K}
weight_window (float) – Poisson mass outside the Fox-Glynn window of the weights
dropped_mass (float) – probability removed by the occupancy threshold
error_bound (float) – sum(pi0) - sum(pit), which IS the L1 error
support_max (int) – largest occupied support over the sweep
support_final (int) – support at the last step
truncated (bool) – True if maxsteps stopped the sweep
absorbed (bool) – True if the support emptied or became absorbing
- ctmc_gmres(A, b, tol=None, restart=None, maxit=None, x0=None)[source]
Solve the sparse nonsymmetric system A*x = b by restarted GMRES.
- Parameters:
A – Coefficient matrix, dense or sparse. Converted to CSC internally.
b – Right-hand side
tol (float | None) – Relative residual tolerance (default 1e-12)
restart (int | None) – Krylov subspace dimension between restarts (default min(n, 50))
maxit (int | None) – Maximum number of restart cycles (default ceil(n/restart))
x0 – Initial guess (default uniform 1/n)
- Returns:
(x, flag, relres, iter), where flag follows the MATLAB gmres convention: 0 converged, 1 iteration limit reached, 2 preconditioner ill-conditioned, 3 stagnation or breakdown. Callers must check flag and fall back to the direct solve when it is nonzero.
- Return type:
- ctmc_gmres_multi(A, B, tol=None, restart=None, maxit=None)[source]
Solve A*X = B for every column of B, reusing one ILUT factorization across all of them and starting each column from the previous solution.
This is the shape of the stochastic complement, whose right-hand side is a whole block of the generator: refactorizing per column would cost more than the direct solve it replaces.
- Parameters:
- Returns:
(X, flag). flag is 0 only if every column converged; on any other value X is None and the caller must fall back to the direct solve. Returning a partial block would leave that fallback ambiguous.
- ctmc_bicgstab(A, b, tol=None, maxit=None, x0=None)[source]
Solve the sparse nonsymmetric system A*x = b by preconditioned BiCGSTAB.
- Parameters:
- Returns:
(x, flag, relres, iter), where flag follows the MATLAB bicgstab convention: 0 converged, 1 iteration limit reached, 2 preconditioner ill-conditioned, 3 stagnation, 4 a scalar quantity became too small or too large to continue. Callers must check flag and fall back to another solve when it is nonzero. iter counts matrix-vector products with A: two per complete iteration, which is what makes it comparable with the iter of ctmc_gmres and across the four codebases. scipy reports only complete iterations, so a solve that converges at a half step is counted here as the full pair.
- Return type:
- ctmc_bicgstab_multi(A, B, tol=None, maxit=None)[source]
Solve A*X = B for every column of B, reusing one ILUT factorization across all of them and starting each column from the previous solution.
This is the shape of the stochastic complement, whose right-hand side is a whole block of the generator: refactorizing per column would cost more than the direct solve it replaces.
- Parameters:
- Returns:
(X, flag). flag is 0 only if every column converged; on any other value X is None and the caller must fall back to another solve. Returning a partial block would leave that fallback ambiguous.
- ctmc_randomization(Q, q=None)[source]
Randomize (uniformize) a CTMC generator into a DTMC.
The randomized DTMC has transition matrix P = I + Q/q, where q is the randomization rate; it carries the same stationary vector as Q. Same name and meaning as MATLAB ctmc_randomization and as the kpctoolbox twin line_solver.lib.kpctoolbox.mc.ctmc_randomization. For the transient distribution by Jensen’s method see ctmc_uniformization.
- ctmc_stochcomp(Q, I=None)[source]
Compute stochastic complement of CTMC.
Reduces the CTMC by eliminating states not in I while preserving the steady-state distribution restricted to the kept states.
- Parameters:
- Returns:
‘S’: Stochastic complement (reduced generator)
’Q11’: Submatrix for kept states
’Q12’: Transitions from kept to eliminated
’Q21’: Transitions from eliminated to kept
’Q22’: Submatrix for eliminated states
’T’: Transient contribution matrix
- Return type:
dict containing
- ctmc_timereverse(Q, pi=None)[source]
Compute time-reversed CTMC generator.
The time-reversed generator Q* has elements: Q*_{ij} = π_j * Q_{ji} / π_i
- ctmc_simulate(Q, initial_state, max_time, max_events=10000, seed=None)[source]
Simulate CTMC sample path using Gillespie algorithm.
Generates a realization of the continuous-time Markov chain using the next-reaction method.
- Parameters:
- Returns:
‘states’: Array of visited states
’times’: Array of transition times
’sojourn_times’: Time spent in each state
- Return type:
dict with
- ctmc_isfeasible(Q, tolerance=1e-10)[source]
Check if matrix is a valid CTMC infinitesimal generator.
Validates: - Off-diagonal elements are non-negative - Row sums are zero - Diagonal elements are non-positive
- ctmc_ssg(sn, options=None)[source]
Generate complete CTMC state space for a queueing network.
Creates all possible network states including those not reachable from the initial state. For open classes, a cutoff parameter limits the maximum population to keep state space finite.
The state space is aggregated to show per-station-class job counts.
- Parameters:
- Returns:
state_space: Complete state space matrix (rows=states, cols=state components)
state_space_aggr: Aggregated state space (rows=states, cols=stations*classes)
state_space_hashed: Hashed state indices for lookup
node_state_space: Dictionary of per-node state spaces
sn: Updated network structure with space field populated
- Return type:
CtmcSsgResult containing
References
MATLAB: matlab/src/api/mc/ctmc_ssg.m
- ctmc_ssg_reachability(sn, options=None)[source]
Generate reachable CTMC state space for a queueing network.
Creates only the states reachable from the initial state through valid transitions. This is more efficient than ctmc_ssg for networks with constrained reachability.
- Parameters:
- Returns:
state_space: Reachable state space matrix
state_space_aggr: Aggregated state space (per station-class)
state_space_hashed: Hashed state indices
node_state_space: Dictionary of per-node state spaces
sn: Updated network structure
- Return type:
CtmcSsgResult containing
References
MATLAB: matlab/src/api/mc/ctmc_ssg_reachability.m
- ctmc_memory_gate(log_nstates, force=False, verbose=False, safety_fraction=0.6)[source]
Hardware-aware, profiling-calibrated CTMC memory pre-gate.
Decides whether a CTMC steady-state solve of a state space of worst-case size exp(log_nstates) is safe on the current host. The budget is a fraction of available memory; the per-state cost is calibrated by profiling sparse LU factorization and cached per machine.
- Parameters:
- Returns:
ok (bool): True if solve is safe, False if memory exceeded (and force=False)
msg (str): Status or warning message
- Return type:
Tuple (ok, msg) where
- class CtmcSsgResult(state_space, state_space_aggr, state_space_hashed, node_state_space, sn)[source]
Bases:
objectResult from CTMC state space generation.
- dtmc_solve(P)[source]
Solve for steady-state probabilities of a DTMC.
Computes the stationary distribution π by solving π(P - I) = 0 with normalization constraint Σπ = 1.
This leverages the CTMC solver by treating (P - I) as an infinitesimal generator.
- dtmc_solve_reducible(P, pin=None)[source]
Solve reducible DTMCs with transient states.
Handles DTMCs with multiple recurrent classes and transient states by: 1. Decomposing into strongly connected components (SCCs) 2. Identifying recurrent vs transient SCCs 3. Computing limiting distribution considering absorption from transient states
For a reducible DTMC with a single transient SCC, this computes the limiting distribution when starting from the transient states (e.g., class switching networks where jobs start in a transient class).
- dtmc_makestochastic(A)[source]
Convert matrix to row-stochastic transition matrix.
Normalizes each row to sum to 1. Rows with zero sum are replaced with uniform distribution.
- dtmc_isfeasible(P, tolerance=1e-10)[source]
Check if matrix is a valid DTMC transition matrix.
Validates: - All elements are non-negative - All row sums equal 1
- dtmc_simulate(P, initial_state, num_steps, seed=None)[source]
Simulate DTMC sample path.
Generates a realization of the discrete-time Markov chain for a specified number of steps.
- dtmc_timereverse(P, pi=None)[source]
Compute time-reversed DTMC transition matrix.
The time-reversed chain has transition probabilities: P*_{ij} = π_j * P_{ji} / π_i
- dtmc_stochcomp(P, keep_states, eliminate_states=None)[source]
Compute stochastic complement of DTMC.
Reduces the DTMC by eliminating specified states while preserving the steady-state distribution restricted to the kept states.
- dtmc_stochcomp_full(P, keep_states, eliminate_states=None)[source]
Stochastic complement of a DTMC together with its blocks.
Same computation as dtmc_stochcomp, additionally returning the four blocks of the transition matrix partitioned by the kept and the eliminated states. Twin of the MATLAB [S,P11,P12,P21,P22] = dtmc_stochcomp(P,I) and of the JAR dtmc_stochcomp_full.
- dtmc_transient(P, initial_dist, steps)[source]
Compute transient probabilities of a DTMC.
Calculates π(n) = π(0) * P^n for each step from 0 to steps.
- dtmc_hitting_time(P, target_states)[source]
Compute mean hitting times to target states.
Calculates the expected number of steps to reach any target state from each starting state.
- ctmc_passage_ph(Q, pi0, target)[source]
Phase-type representation (alpha, S, s0, keep, atom) of the first passage time from pi0 into the target state set.
With A the complement of the target, S = Q(A,A), s0 = -S*1 and alpha = pi0(A), so L(s) = alpha (sI-S)^{-1} s0 + atom and F(t) = 1 - alpha exp(St)1. That is Eqs. 1-2 written as a phase-type law rather than as n scalar equations.
ALPHA IS DELIBERATELY NOT NORMALIZED. Its mass is 1 - atom; the missing mass is the ATOM AT ZERO carried by initial states already inside the target. A caller that normalizes alpha and forgets the atom reports F(0) = 0 for a passage that has already completed with probability atom.
- ctmc_passage_lst(Q, pi0, target, s)[source]
Laplace-Stieltjes transform of the first passage time, L(s) = alpha (sI-S)^{-1} s0 + atom, at the (possibly complex) points s. Eqs. 1-2: one linear system per value of s.
ONE SOLVE PER s, NOT PER (s,t) PAIR. The saving over a dense matrix exponential is that the solves are sparse, so this route reaches chains a dense expm cannot hold. It is NOT a saving in the number of time points: every Abate-Whitt inverter places its nodes at s = beta/t, so a grid of T points costs T*|beta| solves. On a small chain ctmc_passage_time’s default ‘expm’ route is faster.
- ctmc_passage_moments(Q, pi0, target, nmax=1)[source]
Moments of order 1..nmax of the first passage time into the target set.
Returns (mall, m): mall is (nstates, nmax), row i for a passage started in state i, zero on target states and inf where the target cannot be reached; m is the pi0-weighted moment vector.
This is Eq. 3, -q_ii M_i(n) = sum_{k not in B} q_ik M_k(n) + n M_i(n-1), i.e. (-S) M(n) = n M(n-1) with M(0) = 1: nmax linear solves and no transform inversion at all. The equivalent closed form n! alpha (-S)^{-n} 1 is NOT how it is evaluated – forming the inverse of the sub-generator destroys the sparsity the recursion preserves.
- ctmc_passage_time(Q, pi0, target, tset, method='expm', lti_method='euler')[source]
CDF and density of the first passage time into the target state set: F(t) = 1 - alpha exp(St) 1 and f(t) = alpha exp(St) s0.
method=’expm’ (default) is exact and reuses one matrix exponential along a uniform grid. method=’lt’ inverts the transform of Eqs. 1-2 through api/lti; it exists for chains whose non-target block is too large for a dense exp(St), not because it needs fewer time points (see ctmc_passage_lst). On a small chain ‘expm’ is both faster and more accurate, hence the default.
The returned dict carries ‘atom’, the mass of pi0 already inside the target, which is F(0).
- ctmc_hitting_time(Q, target_states)[source]
Mean time to reach any state in target_states from each state of a CTMC.
Continuous-time twin of dtmc_hitting_time and the first-moment special case of ctmc_passage_moments: (-S) h = 1 on the non-target block, where dtmc_hitting_time solves (I - P_NT) h = 1. Unreachable states give inf.
- smp_passage_lst(P, hlst, pi0, target, s)[source]
Laplace-Stieltjes transform of the semi-Markov first passage time, Eqs. 4-5:
L_i(s) = sum_{k not in B} r*_ik(s) L_k(s) + sum_{k in B} r*_ik(s)
so (I - R*_AA(s)) L_A(s) = R*_AB(s) 1, one linear system per value of s.
hlst is either a length-nstates sequence of handles h*_i(s), in which case r*_ik(s) = P[i,k] h*_i(s) and the complex numbers stay on the DIAGONAL of the system (Eq. 5), or an (nstates, nstates) nested sequence of handles r*_ik(s) for the full Markov-renewal kernel, the harder case the paper flags.
Distribution objects supply their own transform: Markovian.evalLST gives the closed form pie (sI-D0)^{-1} (-D0) e for the phase-type family, so hlst[i] = dist.evalLST is the intended way to build these.
- smp_passage_moments(P, hmom, pi0, target, nmax=1)[source]
Moments of order 1..nmax of the first passage time into the target state set for a semi-Markov chain with embedded transition matrix P.
hmom selects which of the paper’s two recursions runs, and they are NOT the same computation:
- (nstates, nmax) array m_i(r), the holding time in i depends only on i.
Eq. 7 with the u_i(r) recurrence of Eq. 8,
u_i(r) = -sum_{j=1..r} C(r,j) m_i(j) u_i(r-j), u_i(0) = 1,
which are the derivatives at the origin of 1/h*_i(s). Cheaper: no per-pair moments.
- (nstates, nstates) list-of-lists hmom[i][k] = [m_ik(1) … m_ik(nmax)],
the r-th moment of the holding time in i WHEN THE NEXT STATE IS k. Eq. 6, the full Markov-renewal kernel. m_ik(0) = P[i,k] is implied and must not be supplied.
Unlike the Markov case the n-th moment needs every moment from 1 to n, so nmax cannot be raised for free.
- smp_passage_time(P, hlst, pi0, target, tset, lti_method='euler')[source]
CDF and density of the semi-Markov first passage time, by inverting smp_passage_lst through api/lti.
There is no matrix-exponential route here: a semi-Markov chain has no generator to exponentiate, which is exactly the case uniformization does not reach and the transform does.
lti_method defaults to ‘euler’ RATHER THAN ‘weeks’. Semi-Markov passage densities are the case Sec. 4.2 singles out as slow-converging for a Laguerre series: a kernel with a deterministic or discontinuous holding time gives a density whose derivatives jump, and laplace_weeks_scaling then refuses by name rather than returning noise.
- class CourtoisResult(p, Qperm, Qdec, eps, epsMAX, P, B, q)[source]
Bases:
objectResult of Courtois decomposition.
- class KMSResult(p, p_1, Qperm, eps, epsMAX, pcourt)[source]
Bases:
objectResult of KMS aggregation-disaggregation.
- class TakahashiResult(p, p_1, pcourt, Qperm, eps, epsMAX)[source]
Bases:
objectResult of Takahashi aggregation-disaggregation.
- ctmc_courtois(Q, MS, q=None)[source]
Courtois decomposition for near-completely decomposable CTMCs.
Decomposes a large CTMC into macrostates and computes approximate steady-state probabilities using hierarchical aggregation.
- Parameters:
- Returns:
CourtoisResult with approximate solution and diagnostics
- Return type:
References
Original MATLAB: matlab/src/api/mc/ctmc_courtois.m Courtois, “Decomposability: Queueing and Computer System Applications”, 1977
- ctmc_kms(Q, MS, numSteps=10)[source]
Koury-McAllister-Stewart aggregation-disaggregation method.
Iteratively refines the Courtois decomposition solution using aggregation and disaggregation steps.
- Parameters:
- Returns:
KMSResult with refined solution
- Return type:
References
Original MATLAB: matlab/src/api/mc/ctmc_kms.m Koury, McAllister, Stewart, “Iterative Methods for Computing Stationary Distributions of Nearly Completely Decomposable Markov Chains”, 1984
- ctmc_takahashi(Q, MS, numSteps=10)[source]
Takahashi’s aggregation-disaggregation method.
Iteratively refines the Courtois decomposition solution using a different aggregation-disaggregation scheme.
- Parameters:
- Returns:
TakahashiResult with refined solution
- Return type:
References
Original MATLAB: matlab/src/api/mc/ctmc_takahashi.m Takahashi, “A Lumping Method for Numerical Calculations of Stationary Distributions of Markov Chains”, 1975
- ctmc_multi(Q, MS, MSS)[source]
Multigrid aggregation-disaggregation method.
Two-level hierarchical decomposition using nested macrostates.
- Parameters:
- Returns:
TakahashiResult with multigrid solution
- Return type:
References
Original MATLAB: matlab/src/api/mc/ctmc_multi.m
- ctmc_saddlepoint(D0, D1=None, t=None, k=None, method=None, pi0=None)[source]
Saddlepoint approximation of Pr{N(t)=k} for the MAP counting process.
- Parameters:
D0 (
array (K,K), orthe MAP pair (D0,D1) / [D0,D1],in which case the) – remaining arguments shift left by one.D1 (
array (K,K),nonnegative. D0+D1 must be an irreducible generator.)t (
floatorarray. Time horizon,broadcast against k.)k (
intorarray. Event count,a nonnegative integer,broadcast against t.)method (
'daniels2'(default) second-order saddlepoint,error O(1/K2**2);) – ‘daniels’ first order with the Perron amplitude, error O(1/K2); ‘plain’ the bare first-order form with the amplitude set to 1.pi0 (
array (K,). Initial phase distribution; the stationary distribution) – of D0+D1 if None.
- Returns:
p (
ndarray. ApproximationofPr```{N(t)=k}`)logp (
ndarray. Its natural logarithm,evaluated without forming p,so it) – stays accurate below the smallest positive double.theta (
ndarray. The saddle theta*,-inf where k=0.)info (
dict with per-point arrays eta,deta,d2eta,d3eta,d4eta,ampl,) – corr, k2 (the expansion parameter), iter, exact, and the scalar lambda.
The ctmc_* and dtmc_* families (steady-state solvers, transient analysis,
uniformization) live in line_solver.api.mc and are listed above.