17227928141

Committed 26 Aug 2025 04:28AM UTC coverage: 92.626%. Remained the same

Build # 17227928141

Build Type

Pull #787

github

Committed by

web-flow

Commit Message

Merge 2d53230c0 into 4a889ad1b

Pull Request Pull Request #787: Migrate np.dot and .dot() method calls to @ operator in library code only

Run Details

113 of 158 new or added lines in 14 files covered. (71.52%)

1 existing line in 1 file now uncovered.

7512 of 8110 relevant lines covered (92.63%)

2.78 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

52.0

/quantecon/markov/gth_solve.py

"""
Routine to compute the stationary distribution of an irreducible Markov
chain by the Grassmann-Taksar-Heyman (GTH) algorithm.

"""
import numpy as np
from numba import jit

from ..util.compat import copy_if_needed


def gth_solve(A, overwrite=False, use_jit=True):
    r"""
    This routine computes the stationary distribution of an irreducible
    Markov transition matrix (stochastic matrix) or transition rate
    matrix (generator matrix) `A`.

    More generally, given a Metzler matrix (square matrix whose
    off-diagonal entries are all nonnegative) `A`, this routine solves
    for a nonzero solution `x` to `x (A - D) = 0`, where `D` is the
    diagonal matrix for which the rows of `A - D` sum to zero (i.e.,
    :math:`D_{ii} = \sum_j A_{ij}` for all :math:`i`). One (and only
    one, up to normalization) nonzero solution exists corresponding to
    each reccurent class of `A`, and in particular, if `A` is
    irreducible, there is a unique solution; when there are more than
    one solution, the routine returns the solution that contains in its
    support the first index `i` such that no path connects `i` to any
    index larger than `i`. The solution is normalized so that its 1-norm
    equals one. This routine implements the Grassmann-Taksar-Heyman
    (GTH) algorithm [1]_, a numerically stable variant of Gaussian
    elimination, where only the off-diagonal entries of `A` are used as
    the input data. For a nice exposition of the algorithm, see Stewart
    [2]_, Chapter 10.

    Parameters
    ----------
    A : array_like(float, ndim=2)
        Stochastic matrix or generator matrix. Must be of shape n x n.
    
    overwrite : bool, optional(default=False)
        Whether to overwrite `A`.
    
    use_jit : bool, optional(default=True)
        Whether to use JIT compilation for performance.

    Returns
    -------
    x : numpy.ndarray(float, ndim=1)
        Stationary distribution of `A`.

    References
    ----------
    .. [1] W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative
       Analysis and Steady State Distributions for Markov Chains,"
       Operations Research (1985), 1107-1116.

    .. [2] W. J. Stewart, Probability, Markov Chains, Queues, and
       Simulation, Princeton University Press, 2009.

    """
    copy = copy_if_needed if overwrite else True

    A1 = np.array(A, dtype=float, copy=copy, order='C')
    # `order='C'` is for use with Numba <= 0.18.2
    # See issue github.com/numba/numba/issues/1103

    if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:
        raise ValueError('matrix must be square')

    n = A1.shape[0]
    x = np.zeros(n)

    if use_jit:
        _gth_solve_jit(A1, x)
        return x

    # if not using jit
    # === Reduction === #
    for k in range(n-1):
        scale = np.sum(A1[k, k+1:n])
        if scale <= 0:
            # There is one (and only one) recurrent class contained in
            # {0, ..., k};
            # compute the solution associated with that recurrent class.
            n = k+1
            break
        A1[k+1:n, k] /= scale

        A1[k+1:n, k+1:n] += A1[k+1:n, k:k+1] @ A1[k:k+1, k+1:n]

    # === Backward substitution === #
    x[n-1] = 1
    for k in range(n-2, -1, -1):
        x[k] = x[k+1:n] @ A1[k+1:n, k]

    # === Normalization === #
    x /= np.sum(x)

    return x


@jit(nopython=True)
def _gth_solve_jit(A, out):
    """
    JIT complied version of the main routine of gth_solve.

    Parameters
    ----------
    A : numpy.ndarray(float, ndim=2)
        Stochastic matrix or generator matrix. Must be of shape n x n.
        Data will be overwritten.

    out : numpy.ndarray(float, ndim=1)
        Output array in which to place the stationary distribution of A.

    """
    n = A.shape[0]

    # === Reduction === #
    for k in range(n-1):
        scale = np.sum(A[k, k+1:n])
        if scale <= 0:
            # There is one (and only one) recurrent class contained in
            # {0, ..., k};
            # compute the solution associated with that recurrent class.
            n = k+1
            break
        for i in range(k+1, n):
            A[i, k] /= scale

            for j in range(k+1, n):
                A[i, j] += A[i, k] * A[k, j]

    # === Backward substitution === #
    out[n-1] = 1
    for k in range(n-2, -1, -1):
        for i in range(k+1, n):
            out[k] += out[i] * A[i, k]

    # === Normalization === #
    norm = np.sum(out)
    for k in range(n):
        out[k] /= norm

1	"""
2	Routine to compute the stationary distribution of an irreducible Markov
3	chain by the Grassmann-Taksar-Heyman (GTH) algorithm.
4
5	"""
6	import numpy as np	3✔
7	from numba import jit	3✔
8
9	from ..util.compat import copy_if_needed	3✔
10
11
12	def gth_solve(A, overwrite=False, use_jit=True):	3✔
13	r"""
14	This routine computes the stationary distribution of an irreducible
15	Markov transition matrix (stochastic matrix) or transition rate
16	matrix (generator matrix) `A`.
17
18	More generally, given a Metzler matrix (square matrix whose
19	off-diagonal entries are all nonnegative) `A`, this routine solves
20	for a nonzero solution `x` to `x (A - D) = 0`, where `D` is the
21	diagonal matrix for which the rows of `A - D` sum to zero (i.e.,
22	:math:`D_{ii} = \sum_j A_{ij}` for all :math:`i`). One (and only
23	one, up to normalization) nonzero solution exists corresponding to
24	each reccurent class of `A`, and in particular, if `A` is
25	irreducible, there is a unique solution; when there are more than
26	one solution, the routine returns the solution that contains in its
27	support the first index `i` such that no path connects `i` to any
28	index larger than `i`. The solution is normalized so that its 1-norm
29	equals one. This routine implements the Grassmann-Taksar-Heyman
30	(GTH) algorithm [1]_, a numerically stable variant of Gaussian
31	elimination, where only the off-diagonal entries of `A` are used as
32	the input data. For a nice exposition of the algorithm, see Stewart
33	[2]_, Chapter 10.
34
35	Parameters
36	----------
37	A : array_like(float, ndim=2)
38	Stochastic matrix or generator matrix. Must be of shape n x n.
39
40	overwrite : bool, optional(default=False)
41	Whether to overwrite `A`.
42
43	use_jit : bool, optional(default=True)
44	Whether to use JIT compilation for performance.
45
46	Returns
47	-------
48	x : numpy.ndarray(float, ndim=1)
49	Stationary distribution of `A`.
50
51	References
52	----------
53	.. [1] W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative
54	Analysis and Steady State Distributions for Markov Chains,"
55	Operations Research (1985), 1107-1116.
56
57	.. [2] W. J. Stewart, Probability, Markov Chains, Queues, and
58	Simulation, Princeton University Press, 2009.
59
60	"""
61	copy = copy_if_needed if overwrite else True	3✔
62
63	A1 = np.array(A, dtype=float, copy=copy, order='C')	3✔
64	# `order='C'` is for use with Numba <= 0.18.2
65	# See issue github.com/numba/numba/issues/1103
66
67	if len(A1.shape) != 2 or A1.shape[0] != A1.shape[1]:	3✔
68	raise ValueError('matrix must be square')	3✔
69
70	n = A1.shape[0]	3✔
71	x = np.zeros(n)	3✔
72
73	if use_jit:	3✔
74	_gth_solve_jit(A1, x)	3✔
75	return x	3✔
76
77	# if not using jit
78	# === Reduction === #
79	for k in range(n-1):	×
80	scale = np.sum(A1[k, k+1:n])	×
81	if scale <= 0:	×
82	# There is one (and only one) recurrent class contained in
83	# {0, ..., k};
84	# compute the solution associated with that recurrent class.
85	n = k+1	×
86	break	×
87	A1[k+1:n, k] /= scale	×
88
NEW 89	A1[k+1:n, k+1:n] += A1[k+1:n, k:k+1] @ A1[k:k+1, k+1:n]	×
90
91	# === Backward substitution === #
92	x[n-1] = 1	×
93	for k in range(n-2, -1, -1):	×
NEW 94	x[k] = x[k+1:n] @ A1[k+1:n, k]	×
95
96	# === Normalization === #
97	x /= np.sum(x)	×
98
99	return x	×
100
101
102	@jit(nopython=True)
103	def _gth_solve_jit(A, out):
104	"""
105	JIT complied version of the main routine of gth_solve.
106
107	Parameters
108	----------
109	A : numpy.ndarray(float, ndim=2)
110	Stochastic matrix or generator matrix. Must be of shape n x n.
111	Data will be overwritten.
112
113	out : numpy.ndarray(float, ndim=1)
114	Output array in which to place the stationary distribution of A.
115
116	"""
117	n = A.shape[0]
118
119	# === Reduction === #
120	for k in range(n-1):
121	scale = np.sum(A[k, k+1:n])
122	if scale <= 0:
123	# There is one (and only one) recurrent class contained in
124	# {0, ..., k};
125	# compute the solution associated with that recurrent class.
126	n = k+1
127	break
128	for i in range(k+1, n):
129	A[i, k] /= scale
130
131	for j in range(k+1, n):
132	A[i, j] += A[i, k] * A[k, j]
133
134	# === Backward substitution === #
135	out[n-1] = 1
136	for k in range(n-2, -1, -1):
137	for i in range(k+1, n):
138	out[k] += out[i] * A[i, k]
139
140	# === Normalization === #
141	norm = np.sum(out)
142	for k in range(n):
143	out[k] /= norm

QuantEcon / QuantEcon.py / 17227928141

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous