Scipy update causes windows test fail in msprime.mutations.TPM() #2349
Description
The test tests/test_mutations.py::TestMatrixMutationModel::test_TPM[0.5-0.25-50-500] fails under Windows with recent numpy.
Line https://github.com/tskit-dev/msprime/blob/main/msprime/mutations.py#L395 finds the values in matrix S that are close to 1.
Under linux this array is:
[ 1.00000000e+00 9.99652651e-01 9.98613097e-01 9.96888767e-01 ....
Windows:
[ 0.99885139 0.98776449 0.98045235 0.99193321 0.96798442 0.99111088 ...
As seen on this CI run:
https://github.com/tskit-dev/msprime/actions/runs/13968809909/job/39105348422
This is likely due to a BLAS library used by np.linalg.eig accumulating floating-point error on Windows.
One way to solve this would be through power iteration:
if root_distribution is None:
# 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:
raise ValueError(
f"Failed to compute stationary distribution after {max_iter} iterations. "
"This may indicate either:\n"
"1. The range between lo and hi is too large (currently {hi-lo+1} states), creating numerical instability\n"
"2. The combination of parameters (particularly s={s}, p={p}, m={m}) creates extremely small transition probabilities\n"
"3. The matrix contains values below the precision threshold of floating-point arithmetic\n"
)See more discussion at
#2348