Courtois decomposition of a nearly completely decomposable (NCD) CTMC. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
#include "line/api/mc/ctmc_randomization.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/api/mc/ctmc_solve_reducible.h"
#include "line/api/mc/dtmc_solve.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for ctmc_courtois.h:

Classes
struct	line::mc::CourtoisResult< T >

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
CourtoisResult< T >	line::mc::ctmc_courtois (const Matrix< T > &Q, const std::vector< std::vector< std::size_t > > &MS, const T &q)
	Courtois decomposition of a nearly completely decomposable (NCD) CTMC.
template<class T>
CourtoisResult< T >	line::mc::ctmc_courtois (const Matrix< T > &Q, const std::vector< std::vector< std::size_t > > &MS)
	Overload deriving the rate as MATLAB does, q = (21/20) max\|Qperm\|.

Detailed Description

Courtois decomposition of a nearly completely decomposable (NCD) CTMC.

Templated port of matlab/src/api/mc/ctmc_courtois.m and jar/src/main/java/jline/api/mc/Ctmc_courtois.java. The states are permuted into macro-state order; the coupling between macro-states is deleted to give a block-diagonal generator Qdec whose blocks are solved independently for the conditional (micro) distributions; the macro-state chain G is assembled from the uniformized matrix weighted by those micro distributions; and the approximation is the product of the macro and micro probabilities, mapped back to the original state ordering.

The approximation error is governed by the degree of coupling eps, which is meaningful only when it is small against epsMAX = (1 - max_i mu_i)/2, mu_i being the subdominant eigenvalue modulus of the i-th diagonal block made stochastic. Both diagnostics are always computed here; MATLAB computes them only when more than three outputs are requested and returns eps = 0 otherwise, which is a trap for a caller who asks for two outputs and reads the second as a coupling measure.

WHICH SUM eps IS, and the divergence that used to live here. The NCD index of Courtois is the largest ROW sum of the coupling matrix, ||B||_inf, the largest probability of leaving a macro-state in one uniformized step. MATLAB once wrote it max(sum(B)), and sum(B) on a matrix is the vector of COLUMN sums, so its eps was the largest column sum instead; Ctmc_courtois.java computed the row sum (B.sumRows().elementMax()) and the two separated whenever any state carried coupling to more than one other macro-state. MATLAB is max(sum(B,2)) as of 2026-08-15 and the three codebases agree, so eps here is the ROW sum. epsColMax keeps the old column-sum value, because a golden or a comparison recorded before that date holds it: on the six-state fixture of test_mc_aggregation the row sum is 0.0076190 and the column sum 0.013333, a factor of 1.75, and the column sum is the OPTIMISTIC one – it understates the coupling, which is the wrong direction for a diagnostic answering "is this partition decomposable enough to trust".

GATED ON TRANSCENDENTAL ARITHMETIC. epsMAX is a subdominant eigenvalue modulus, obtained by the Francis QR iteration, which is iterative and stops on a tolerance; it has no closed form in the field of the rates (the eigenvalues of a rational matrix are algebraic, not rational). The rest of the method – permutation, block solves, aggregation, disaggregation – is finite and exact, so a caller who needs the approximate stationary vector at Rational and no diagnostics can assemble it from ctmc_solve_reducible and dtmc_solve directly.

Definition in file ctmc_courtois.h.

Classes

Namespaces

Functions

Detailed Description