Compute stationary distribution using power iteration to avoid numerical issue on Windows

Author: benjeffery

msprime/mutations.py

@@ -388,13 +388,29 @@ def __init__(
             raise ValueError("u must be between 0 and 1")
         alleles = [str(int(x)) for x in range(lo, hi + 1)]
         transition_matrix = general_microsat_rate_matrix(s, u, v, p, m, lo, hi)
+
+        # Print for debugging (keep existing print statements)
+        print("Transition matrix:")
+        print(transition_matrix, np.sum(transition_matrix))
+        print("Transition matrix T:")
+        print(transition_matrix.T, np.sum(transition_matrix.T))
+
         if root_distribution is None:
-            # solve for stationary distribution
-            S, U = np.linalg.eig(transition_matrix.T)
-            U = np.real_if_close(U, tol=1)
-            stationary = np.array(U[:, np.where(np.abs(S - 1.0) < 1e-8)[0][0]])
-            stationary = stationary / np.sum(stationary)
-            root_distribution = stationary
+            # Compute the stationary distribution using power iteration
+            n = transition_matrix.shape[0]
+            pi = np.ones(n) / n
+            max_iter = 1000
+            tol = 1e-10
+            for _ in range(max_iter):
+                pi_next = pi @ transition_matrix
+                if np.max(np.abs(pi_next - pi)) < tol:
+                    root_distribution = pi_next
+                    break
+                pi = pi_next
+            else:
+                # If the loop completes without breaking, use the last approximation
+                root_distribution = pi
+
         super().__init__(alleles, root_distribution, transition_matrix)