amachine.am_solve

View Source

 1import numpy as np
 2import scipy as sp
 3import sympy
 4from fractions import Fraction
 5import warnings
 6
 7def solve_for_pi( T : np.ndarray, EPS : float = 1e-12 ) -> np.ndarray :
 8
 9    n = T.shape[0]
10
11    # helper ---------------------------------------------
12
13    def normalize_pi( pi ) :
14        if np.any( pi < -1e-8 ):
15            raise ValueError("Significant negative values in solve result")
16
17        pi = np.clip( pi, a_min=0.0, a_max=None )
18        s = pi.sum()
19
20        if abs(s) < EPS :
21            raise Exception("Degenerate stationary distribution (sum ≈ 0)")
22
23        return pi / s
24
25    # least squares method --------------------------------
26
27    try:
28        A = T.T - np.eye(n)
29        A = np.vstack([np.eye(n) - T.T, np.ones((1, n))])
30        b = np.zeros(n + 1)
31        b[-1] = 1.0
32        pi, *_ = np.linalg.lstsq(A, b, rcond=None)
33        return normalize_pi(pi)
34
35    except (np.linalg.LinAlgError, ValueError):
36        warnings.warn( "Falling back to linalg.solve" )
37
38    # Solve method ----------------------------------------
39
40    try:
41        A = T.T - np.eye(n)
42        A[-1, :] = 1.0
43        b = np.zeros(n)
44        b[-1] = 1.0
45        pi = np.linalg.solve(A, b)
46        return normalize_pi( pi )
47
48    except (np.linalg.LinAlgError, ValueError):
49        warnings.warn( "Falling back to nullspace method" )
50
51    # Nullspace fallback -------------------------------
52
53    try:
54        A = T.T - np.eye(n)
55        ns = sp.linalg.null_space( A )
56
57        if ns.shape[1] != 1:
58            raise ValueError("Unexpected nullspace dimension")
59
60        pi = ns[:, 0]
61        pi *= np.sign(pi[np.argmax(np.abs(pi))])
62        return normalize_pi( pi )
63
64    except (np.linalg.LinAlgError, ValueError):
65        warnings.warn( "Falling back to eigen method" )
66
67    # Eigen fallback ----------------------------------
68
69    eigvals, eigvecs = np.linalg.eig( T.T )
70    idx = np.argmin(np.abs(eigvals - 1))
71    pi = np.real_if_close(eigvecs[:, idx])
72    pi = np.real(pi)
73    pi *= np.sign(pi[np.argmax(np.abs(pi))])
74
75    return normalize_pi( pi )
76
77def solve_for_pi_fractional( T : sympy.Matrix ) -> tuple[Fraction, ...] :
78
79    n = T.shape[0]
80
81    # Setup the system (T.T - I)
82    A = T.transpose() - sympy.Matrix.eye(n)
83
84    # Replace the last row with 1s to enforce sum(pi) = 1
85    for j in range(n):
86        A[n-1, j] = 1
87
88    # Setup the target vector
89    b = sympy.Matrix([0] * (n - 1) + [1])
90
91    # Solve
92    pi_exact = A.LUsolve(b)
93    pi_exact = ( Fraction( int(x.p), int(x.q) ) for x in pi_exact )
94
95    return tuple(pi_exact)

def solve_for_pi(T: numpy.ndarray, EPS: float = 1e-12) -> numpy.ndarray: View Source

 8def solve_for_pi( T : np.ndarray, EPS : float = 1e-12 ) -> np.ndarray :
 9
10    n = T.shape[0]
11
12    # helper ---------------------------------------------
13
14    def normalize_pi( pi ) :
15        if np.any( pi < -1e-8 ):
16            raise ValueError("Significant negative values in solve result")
17
18        pi = np.clip( pi, a_min=0.0, a_max=None )
19        s = pi.sum()
20
21        if abs(s) < EPS :
22            raise Exception("Degenerate stationary distribution (sum ≈ 0)")
23
24        return pi / s
25
26    # least squares method --------------------------------
27
28    try:
29        A = T.T - np.eye(n)
30        A = np.vstack([np.eye(n) - T.T, np.ones((1, n))])
31        b = np.zeros(n + 1)
32        b[-1] = 1.0
33        pi, *_ = np.linalg.lstsq(A, b, rcond=None)
34        return normalize_pi(pi)
35
36    except (np.linalg.LinAlgError, ValueError):
37        warnings.warn( "Falling back to linalg.solve" )
38
39    # Solve method ----------------------------------------
40
41    try:
42        A = T.T - np.eye(n)
43        A[-1, :] = 1.0
44        b = np.zeros(n)
45        b[-1] = 1.0
46        pi = np.linalg.solve(A, b)
47        return normalize_pi( pi )
48
49    except (np.linalg.LinAlgError, ValueError):
50        warnings.warn( "Falling back to nullspace method" )
51
52    # Nullspace fallback -------------------------------
53
54    try:
55        A = T.T - np.eye(n)
56        ns = sp.linalg.null_space( A )
57
58        if ns.shape[1] != 1:
59            raise ValueError("Unexpected nullspace dimension")
60
61        pi = ns[:, 0]
62        pi *= np.sign(pi[np.argmax(np.abs(pi))])
63        return normalize_pi( pi )
64
65    except (np.linalg.LinAlgError, ValueError):
66        warnings.warn( "Falling back to eigen method" )
67
68    # Eigen fallback ----------------------------------
69
70    eigvals, eigvecs = np.linalg.eig( T.T )
71    idx = np.argmin(np.abs(eigvals - 1))
72    pi = np.real_if_close(eigvecs[:, idx])
73    pi = np.real(pi)
74    pi *= np.sign(pi[np.argmax(np.abs(pi))])
75
76    return normalize_pi( pi )

def solve_for_pi_fractional( T: sympy.matrices.dense.MutableDenseMatrix) -> tuple[fractions.Fraction, ...]: View Source

78def solve_for_pi_fractional( T : sympy.Matrix ) -> tuple[Fraction, ...] :
79
80    n = T.shape[0]
81
82    # Setup the system (T.T - I)
83    A = T.transpose() - sympy.Matrix.eye(n)
84
85    # Replace the last row with 1s to enforce sum(pi) = 1
86    for j in range(n):
87        A[n-1, j] = 1
88
89    # Setup the target vector
90    b = sympy.Matrix([0] * (n - 1) + [1])
91
92    # Solve
93    pi_exact = A.LUsolve(b)
94    pi_exact = ( Fraction( int(x.p), int(x.q) ) for x in pi_exact )
95
96    return tuple(pi_exact)