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_transient: CTMC transient analysis ctmc_uniformization: Uniformization method for transient analysis ctmc_foxglynn: Fox-Glynn uniformization for transient analysis ctmc_gmres: Restarted GMRES with ILUT preconditioning 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:
Q (numpy.ndarray) – Infinitesimal generator matrix (row sums should be zero)
method (str | None) – ‘gmres’ or ‘direct’ to force a solution method; None or ‘default’ selects by size, GMRES above GMRES_MIN_STATES states
- 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.
- Parameters:
Q (numpy.ndarray) – Infinitesimal generator matrix (possibly reducible)
pin (numpy.ndarray | None) – Initial probability vector, or None when not available
- Returns:
Steady-state probability vector
- Return type:
- 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')[source]
Compute transient probabilities of a CTMC.
Calculates time-dependent state probabilities π(t) for specified time points using matrix exponential methods.
- Parameters:
Q (numpy.ndarray) – Infinitesimal generator matrix
initial_dist (numpy.ndarray) – Initial probability distribution π(0)
time_points (float | numpy.ndarray) – Array of time points to evaluate, or single time value
method (str) – ‘expm’ for matrix exponential, ‘ode’ for ODE solver
- 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=100)[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.
- Returns:
(piTimeAvg, piExit) as 1D arrays.
- ctmc_uniformization(Q, lambda_rate=None)[source]
Uniformize CTMC generator matrix.
Converts CTMC to an equivalent uniformized discrete-time chain for numerical analysis and simulation purposes.
The uniformized DTMC has transition matrix P = I + Q/λ where λ is the uniformization rate (max exit rate).
- Parameters:
Q (numpy.ndarray) – Infinitesimal generator matrix
lambda_rate (float | None) – Uniformization rate (optional, auto-computed if None)
- Returns:
‘P’: Uniformized transition matrix
’lambda’: Uniformization rate
- Return type:
dict containing
- ctmc_foxglynn(pi0, Q, t, tol=1e-12, maxiter=-1)[source]
Transient distribution of a CTMC by Fox-Glynn uniformization.
- Parameters:
pi0 (numpy.ndarray) – Initial probability distribution
Q (numpy.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)[source]
Fox-Glynn truncation window and normalized 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.
- Parameters:
- Returns:
Tuple of (left truncation point, right truncation point, weights)
- Return type:
Tuple[int, int, numpy.ndarray]
- 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:
Tuple[numpy.ndarray, int, float, int]
- 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_randomization(Q, initial_dist, time_points, precision=1e-10)[source]
Compute CTMC transient probabilities using randomization.
Uses Jensen’s randomization method (uniformization) to compute transient probabilities by converting the CTMC to a uniformized DTMC.
This method is numerically stable and avoids matrix exponentials.
- Parameters:
Q (numpy.ndarray) – Infinitesimal generator matrix
initial_dist (numpy.ndarray) – Initial probability distribution
time_points (numpy.ndarray) – Array of time points to evaluate
precision (float) – Numerical precision for truncation (Poisson tail)
- Returns:
Transient probabilities at each time point
- Return type:
- 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:
Q (numpy.ndarray) – Infinitesimal generator matrix
I (numpy.ndarray | None) – States to retain (array of indices). If None, defaults to 0..ceil(n/2)-1 (matching JAR/MATLAB).
- 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
- Parameters:
Q (numpy.ndarray) – Original infinitesimal generator matrix
pi (numpy.ndarray | None) – Steady-state distribution (optional, computed if None)
- Returns:
Time-reversed generator matrix
- Return type:
- ctmc_rand(n, density=0.3, max_rate=10.0)[source]
Generate random CTMC generator matrix.
- Parameters:
- Returns:
Random infinitesimal generator matrix
- Return type:
- 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:
Q (numpy.ndarray) – Infinitesimal generator matrix
initial_state (int) – Starting state (integer index)
max_time (float) – Maximum simulation time
max_events (int) – Maximum number of transitions (default: 10000)
seed (int | None) – Random seed for reproducibility (optional)
- 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
- Parameters:
Q (numpy.ndarray) – Candidate generator matrix
tolerance (float) – Numerical tolerance (default: 1e-10)
- Returns:
True if matrix is valid CTMC generator
- Return type:
- 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.
- state_space: numpy.ndarray
- state_space_aggr: numpy.ndarray
- state_space_hashed: numpy.ndarray
- node_state_space: Dict[int, numpy.ndarray]
- 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.
- Parameters:
P (numpy.ndarray) – Transition probability matrix (row stochastic)
- Returns:
Steady-state probability distribution (1D array)
- Return type:
- 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).
- Parameters:
P (numpy.ndarray) – Transition probability matrix (possibly reducible)
pin (numpy.ndarray) – Initial probability vector (optional)
- Returns:
Steady-state probability vector
- Return type:
- 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.
- Parameters:
A (numpy.ndarray) – Input matrix to normalize
- Returns:
Row-stochastic matrix
- Return type:
- 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
- Parameters:
P (numpy.ndarray) – Candidate transition matrix
tolerance (float) – Numerical tolerance (default: 1e-10)
- Returns:
True if matrix is valid DTMC transition matrix
- Return type:
- 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.
- Parameters:
P (numpy.ndarray) – Transition probability matrix
initial_state (int) – Starting state index
num_steps (int) – Number of simulation steps
- Returns:
Array of visited states (length num_steps + 1)
- Return type:
- dtmc_rand(n, density=0.5)[source]
Generate random DTMC transition matrix.
- Parameters:
- Returns:
Random transition probability matrix
- Return type:
- 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
- Parameters:
P (numpy.ndarray) – Original transition matrix
pi (numpy.ndarray | None) – Steady-state distribution (optional, computed if None)
- Returns:
Time-reversed transition matrix
- Return type:
- 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.
- Parameters:
P (numpy.ndarray) – Transition probability matrix
keep_states (numpy.ndarray) – States to retain in reduced model
eliminate_states (numpy.ndarray | None) – States to eliminate (optional, inferred if None)
- Returns:
Reduced transition matrix (stochastic complement)
- Return type:
- 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.
- Parameters:
P (numpy.ndarray) – Transition probability matrix
initial_dist (numpy.ndarray) – Initial probability distribution π(0)
steps (int) – Number of time steps
- Returns:
Array of shape (steps+1, n) with transient probabilities
- Return type:
- 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.
- Parameters:
P (numpy.ndarray) – Transition probability matrix
target_states (numpy.ndarray) – Array of target state indices
- Returns:
Array of mean hitting times from each state
- Return type:
- class CourtoisResult(p, Qperm, Qdec, eps, epsMAX, P, B, q)[source]
Bases:
objectResult of Courtois decomposition.
- Qperm: numpy.ndarray
- Qdec: numpy.ndarray
- class KMSResult(p, p_1, Qperm, eps, epsMAX, pcourt)[source]
Bases:
objectResult of KMS aggregation-disaggregation.
- p_1: numpy.ndarray
- Qperm: numpy.ndarray
- pcourt: numpy.ndarray
- class TakahashiResult(p, p_1, pcourt, Qperm, eps, epsMAX)[source]
Bases:
objectResult of Takahashi aggregation-disaggregation.
- p_1: numpy.ndarray
- pcourt: numpy.ndarray
- Qperm: numpy.ndarray
- 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:
Q (numpy.ndarray) – Infinitesimal generator matrix
MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i
q (float | None) – Randomization coefficient (optional)
- 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:
Q (numpy.ndarray) – Infinitesimal generator matrix
MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i
numSteps (int) – Number of iterative steps (default: 10)
- 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:
Q (numpy.ndarray) – Infinitesimal generator matrix
MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i
numSteps (int) – Number of iterative steps (default: 10)
- 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
Continuous-Time Markov Chains (line_solver.api.ctmc)
The ctmc module contains algorithms specifically for continuous-time Markov
chains (CTMCs), including steady-state solvers and transient analysis.
Discrete-Time Markov Chains (line_solver.api.dtmc)
The dtmc module provides algorithms for discrete-time Markov chains (DTMCs),
including steady-state and transient analysis.