amachine.am_solve 

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