From 17d39896b07979ecba40eecfc252f8f5e56f4a61 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Thu, 25 Sep 2025 16:36:40 +0900
Subject: [PATCH 1/6] Hopefully fix loky MPI deadlock

---
 qutip/solver/parallel.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/qutip/solver/parallel.py b/qutip/solver/parallel.py
index 425a317f21..dd9e84d41d 100644
--- a/qutip/solver/parallel.py
+++ b/qutip/solver/parallel.py
@@ -483,6 +483,26 @@ def mpi_pmap(task, values, task_args=None, task_kwargs=None,
 
     from mpi4py.futures import MPIPoolExecutor
 
+    # If loky is used (see loky_pmap), it starts two "resource tracker" child
+    # processes that are normally never shut down. However, using mpi4py
+    # when these child processes are still alive may lead to a deadlock.
+    # -- WE RECOMMEND NOT TO MIX loky_pmap AND mpi_pmap IN THE SAME PROGRAM --
+    # However, it happens in our test runs. For this reason, we here try to
+    # shut down the resource tracker processes manually
+    try:
+        import loky
+        loky_resource_tracker = loky.backend.resource_tracker._resource_tracker
+        if loky_resource_tracker:
+            loky_resource_tracker._stop()
+
+        if hasattr(multiprocessing, 'resource_tracker'):
+            mp_resource_tracker =\
+                multiprocessing.resource_tracker._resource_tracker
+            if mp_resource_tracker:
+                mp_resource_tracker._stop()
+    except ImportError:
+        pass
+
     # If the provided num_cpus is None, we use the default value instead.
     # We thus intentionally make it impossible to call
     #   MPIPoolExecutor(max_workers=None, ...)

From 3a936e4bbab51c729fe6004887fa82895cfb1d98 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Thu, 25 Sep 2025 16:37:05 +0900
Subject: [PATCH 2/6] Correct calculation of root mean square error

---
 qutip/utilities.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qutip/utilities.py b/qutip/utilities.py
index 52ceaa55fd..79f18a9afb 100644
--- a/qutip/utilities.py
+++ b/qutip/utilities.py
@@ -490,7 +490,7 @@ def _rmse(fun, xdata, ydata, params):
     if (yhat == ydata).all():
         return 0
     return (
-        np.sqrt(np.mean((yhat - ydata) ** 2) / len(ydata))
+        np.sqrt(np.mean(np.abs(yhat - ydata) ** 2))
         / (np.max(ydata) - np.min(ydata))
     )
 

From 9c7318453c71212067373abd9ed0f576dff2ca8e Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Thu, 25 Sep 2025 16:37:43 +0900
Subject: [PATCH 3/6] Avoid calling np.random.seed in tests, update fitting
 tolerances for fixed rmse calculation

---
 qutip/tests/core/test_brtools.py     | 11 +++++-----
 qutip/tests/core/test_environment.py | 32 ++++++++++++++--------------
 qutip/tests/solver/test_results.py   |  1 +
 qutip/tests/test_mkl.py              |  4 ++--
 qutip/tests/test_utilities.py        | 23 ++++++++++++--------
 5 files changed, 38 insertions(+), 33 deletions(-)

diff --git a/qutip/tests/core/test_brtools.py b/qutip/tests/core/test_brtools.py
index 095e40644a..77822366b2 100644
--- a/qutip/tests/core/test_brtools.py
+++ b/qutip/tests/core/test_brtools.py
@@ -10,9 +10,8 @@
 )
 
 
-def _make_rand_data(shape):
-    np.random.seed(11)
-    array = np.random.rand(*shape) + 1j*np.random.rand(*shape)
+def _make_rand_data(shape, rng):
+    array = rng.random(shape) + 1j*rng.random(shape)
     return qutip.data.Dense(array)
 
 
@@ -31,9 +30,9 @@ def _make_rand_data(shape):
                          ids=['', 'transpose', 'conj', 'dag'])
 def test_matmul_var(datatype, transleft, transright):
     shape = (5, 5)
-    np.random.seed(11)
-    left = qutip.data.to(datatype, _make_rand_data(shape))
-    right = qutip.data.to(datatype, _make_rand_data(shape))
+    rng = np.random.default_rng(seed=11)
+    left = qutip.data.to(datatype, _make_rand_data(shape, rng))
+    right = qutip.data.to(datatype, _make_rand_data(shape, rng))
 
     expected = qutip.data.matmul(
         transform[transleft](left),
diff --git a/qutip/tests/core/test_environment.py b/qutip/tests/core/test_environment.py
index 70f8253fa7..22399a9a43 100644
--- a/qutip/tests/core/test_environment.py
+++ b/qutip/tests/core/test_environment.py
@@ -532,8 +532,8 @@ def test_fixed_cf_fit(self, reference, tMax, full_ansatz, tol):
 
         assert info["Nr"] == 2
         assert info["Ni"] == 2
-        assert info["rmse_real"] < 5e-3
-        assert info["rmse_imag"] < 5e-3
+        assert info["rmse_real"] < 5e-2
+        assert info["rmse_imag"] < 5e-2
         for key in ["fit_time_real", "fit_time_imag",
                     "params_real", "params_imag", "summary"]:
             assert key in info
@@ -551,7 +551,7 @@ def test_dynamic_cf_fit(self, reference, tMax, tol, full_ansatz):
         )
         tlist = np.linspace(0, tMax, 100)[1:]  # exclude t=0
         fit, info = env.approximate(
-            "cf", tlist, target_rmse=0.01, Nr_max=3, Ni_max=3,
+            "cf", tlist, target_rmse=0.1, Nr_max=3, Ni_max=3,
             full_ansatz=full_ansatz
         )
 
@@ -564,8 +564,8 @@ def test_dynamic_cf_fit(self, reference, tMax, tol, full_ansatz):
 
         assert info["Nr"] == 1
         assert info["Ni"] == 1
-        assert info["rmse_real"] <= 0.01
-        assert info["rmse_imag"] <= 0.01
+        assert info["rmse_real"] <= 0.1
+        assert info["rmse_imag"] <= 0.1
         for key in ["fit_time_real", "fit_time_imag",
                     "params_real", "params_imag", "summary"]:
             assert key in info
@@ -594,7 +594,7 @@ def test_fixed_sd_fit(self, reference, wMax, tol):
 
         assert info["N"] == 4
         assert info["Nk"] == 1
-        assert info["rmse"] < 5e-3
+        assert info["rmse"] < 5e-2
         for key in ["fit_time", "params", "summary"]:
             assert key in info
 
@@ -611,7 +611,7 @@ def test_dynamic_sd_fit(self, reference, wMax, tol, params):
         )
         wlist = np.linspace(0, wMax, 100)
         fit, info = env.approximate(
-            "sd", wlist, Nk=1, target_rmse=0.01, Nmax=5, **params
+            "sd", wlist, Nk=1, target_rmse=0.1, Nmax=5, **params
         )
 
         assert isinstance(fit, ExponentialBosonicEnvironment)
@@ -623,14 +623,14 @@ def test_dynamic_sd_fit(self, reference, wMax, tol, params):
 
         assert info["N"] < 5
         assert info["Nk"] == 1
-        assert info["rmse"] < 0.01
+        assert info["rmse"] < 0.1
         for key in ["fit_time", "params", "summary"]:
             assert key in info
 
     @pytest.mark.parametrize(["reference", "tMax", "N", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1),
-                     15, 8, 1e-3, id="Ohmic Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 4, 1e-3, id='UD Example'),
+                     15, 8, 1e-2, id="Ohmic Example"),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 4, 1e-2, id='UD Example'),
     ])
     @pytest.mark.parametrize("separate", [True, False])
     def test_fixed_prony_fit(self, reference, tMax, N, tol, separate):
@@ -665,8 +665,8 @@ def test_fixed_prony_fit(self, reference, tMax, N, tol, separate):
 
     @pytest.mark.parametrize(["reference", "tMax", "N", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1),
-                     15, 8, 1e-3, id="Ohmic Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 4, 1e-3, id='UD Example'),
+                     15, 8, 1e-2, id="Ohmic Example"),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 4, 1e-2, id='UD Example'),
     ])
     @pytest.mark.parametrize("separate", [True, False])
     def test_fixed_esprit_fit(self, reference, tMax, N, tol, separate):
@@ -701,8 +701,8 @@ def test_fixed_esprit_fit(self, reference, tMax, N, tol, separate):
 
     @pytest.mark.parametrize(["reference", "tMax", "N", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1),
-                     15, 8, 1e-3, id="Ohmic Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 2, 1e-3, id='UD Example'),
+                     15, 8, 1e-2, id="Ohmic Example"),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 2, 1e-2, id='UD Example'),
     ])
     @pytest.mark.parametrize("separate", [True, False])
     def test_fixed_espira1_fit(self, reference, tMax, N, tol, separate):
@@ -737,8 +737,8 @@ def test_fixed_espira1_fit(self, reference, tMax, N, tol, separate):
 
     @pytest.mark.parametrize(["reference", "tMax", "N", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1),
-                     15, 8, 1e-3, id="Ohmic Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 2, 1e-3, id='UD Example'),
+                     15, 8, 1e-2, id="Ohmic Example"),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 2, 1e-2, id='UD Example'),
     ])
     @pytest.mark.parametrize("separate", [True, False])
     def test_fixed_espira2_fit(self, reference, tMax, N, tol, separate):
diff --git a/qutip/tests/solver/test_results.py b/qutip/tests/solver/test_results.py
index fd3752dce1..2835857d6e 100644
--- a/qutip/tests/solver/test_results.py
+++ b/qutip/tests/solver/test_results.py
@@ -400,6 +400,7 @@ def test_repr(self, keep_runs_results):
         if keep_runs_results:
             assert "Trajectories saved." in repr
 
+    @pytest.mark.flaky(reruns=2)
     @pytest.mark.parametrize('keep_runs_results1', [True, False])
     @pytest.mark.parametrize('keep_runs_results2', [True, False])
     def test_merge_result(self, keep_runs_results1, keep_runs_results2):
diff --git a/qutip/tests/test_mkl.py b/qutip/tests/test_mkl.py
index 734f33e25e..647d3225c5 100644
--- a/qutip/tests/test_mkl.py
+++ b/qutip/tests/test_mkl.py
@@ -19,8 +19,8 @@ def test_single_rhs_vector_real(self):
                            [1, 0, 1],
                            [0, 0, 1]])
         As = scipy.sparse.csr_matrix(Adense)
-        np.random.seed(1234)
-        x = np.random.randn(3)
+        rng = np.random.default_rng(seed=1234)
+        x = rng.standard_normal(3)
         b = As * x
         x2 = mkl_spsolve(As, b, verbose=True)
         np.testing.assert_allclose(x, x2)
diff --git a/qutip/tests/test_utilities.py b/qutip/tests/test_utilities.py
index dab481ef30..8fd8755dc9 100644
--- a/qutip/tests/test_utilities.py
+++ b/qutip/tests/test_utilities.py
@@ -124,6 +124,8 @@ def test_cpu_count(monkeypatch):
 
 
 class TestFitting:
+    rng = np.random.default_rng(seed=42)
+
     def model(self, x, a, b, c):
         return np.real(a * np.exp(-(b + 1j * c) * x))
 
@@ -145,8 +147,7 @@ def generate_data(self, request):
         fparams2 = [3, 2, .5]
         y = self.model(x, *fparams1) + self.model(x, *fparams2)
         if noisy:
-            np.random.seed(42)  
-            noise = np.random.normal(0, 0.01, len(x))
+            noise = self.rng.normal(0, 0.01, len(x))
             y += noise
         return x, y, fparams1, fparams2, noisy
 
@@ -156,11 +157,15 @@ def test_fit(self, generate_data):
             self.model, num_params=3, xdata=x, ydata=y,
             lower=[-np.inf, -np.inf, 0], target_rmse=1e-8, Nmax=2
         )
+        fit_y = self.model(x, *params[0]) + self.model(x, *params[1])
+        assert rmse == pytest.approx(
+            np.sqrt(np.mean((y - fit_y)**2)) / (np.max(y) - np.min(y))
+        )
+
         if noisy:
-            assert rmse < 1e-3
-            # The atol is the std of noise times the maximum of the signal
-            assert (np.all(np.isclose(params, [fparams1, fparams2], atol=1e-2*np.max(y))) or
-                    np.all(np.isclose(params, [fparams2, fparams1], atol=1e-2*np.max(y))))
+            assert rmse < 1e-2
+            assert (np.all(np.isclose(params, [fparams1, fparams2], atol=.2)) or
+                    np.all(np.isclose(params, [fparams2, fparams1], atol=.2)))
         else:
             assert rmse < 1e-8
             assert (np.all(np.isclose(params, [fparams1, fparams2], atol=1e-3)) or
@@ -181,7 +186,7 @@ def test_espira_I(self, generate_data):
         x, y, _, _, noisy = generate_data
         rmse, params = utils.espira1(y, 4, tol=1e-16)
         if noisy:
-            assert rmse < 1e-3
+            assert rmse < 1e-2
             np.testing.assert_allclose(self.eval_prony(len(x), params), y, atol=1e-2*np.max(y))
         else:
             assert rmse < 1e-8
@@ -191,7 +196,7 @@ def test_espira_II(self, generate_data):
         x, y, _, _, noisy = generate_data
         rmse, params = utils.espira2(y, 4, tol=1e-16)
         if noisy:
-            assert rmse < 1e-3
+            assert rmse < 1e-2
             np.testing.assert_allclose(self.eval_prony(len(x), params), y, atol=1e-2*np.max(y))
         else:
             assert rmse < 1e-8
@@ -202,7 +207,7 @@ def test_prony_methods(self, generate_data, method):
         x, y, _, _, noisy = generate_data
         rmse, params = utils.prony_methods(method, y, 4)
         if noisy:
-            assert rmse < 1e-3
+            assert rmse < 1e-2
             np.testing.assert_allclose(self.eval_prony(len(x), params), y, atol=2e-2*np.max(y))
         else:
             assert rmse < 1e-8

From 4d1d03bfc7faba5b8334659f6ed92f3dd3fdb8f9 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Fri, 26 Sep 2025 11:01:17 +0900
Subject: [PATCH 4/6] Add another flaky mark

---
 qutip/tests/solver/test_results.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/qutip/tests/solver/test_results.py b/qutip/tests/solver/test_results.py
index 2835857d6e..7a33316919 100644
--- a/qutip/tests/solver/test_results.py
+++ b/qutip/tests/solver/test_results.py
@@ -283,6 +283,7 @@ def test_NmmcResult(self, include_no_jump, martingale, result_trace):
         for i, (s1, s2) in enumerate(zip(m_res.average_states, result_trace)):
             assert s1 == s2 * qutip.fock_dm(10, i)
 
+    @pytest.mark.flaky(reruns=2)
     @pytest.mark.parametrize('keep_runs_results', [True, False])
     @pytest.mark.parametrize('include_no_jump', [True, False])
     @pytest.mark.parametrize(["e_ops", "results"], [

From d201a03e12a652ce4cf67fa5384394ce4290f863 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Fri, 26 Sep 2025 13:11:06 +0900
Subject: [PATCH 5/6] Update tolerances

---
 qutip/tests/core/test_environment.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/qutip/tests/core/test_environment.py b/qutip/tests/core/test_environment.py
index 22399a9a43..2c60a3e45f 100644
--- a/qutip/tests/core/test_environment.py
+++ b/qutip/tests/core/test_environment.py
@@ -775,7 +775,7 @@ def test_fixed_espira2_fit(self, reference, tMax, N, tol, separate):
 
     @pytest.mark.parametrize(["reference", "wMax", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1), 15, .2, id="Ohmic Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 1e-4, id='UD Example'),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 1e-3, id='UD Example'),
     ])
     def test_fixed_aaa_fit(self, reference, wMax, tol):
         env = BosonicEnvironment.from_spectral_density(
@@ -797,7 +797,7 @@ def test_fixed_aaa_fit(self, reference, wMax, tol):
 
     @pytest.mark.parametrize(["reference", "wMax", "tol"], [
         pytest.param(OhmicReference(3, .75, 10, 1), 15, .2, id="DL Example"),
-        pytest.param(UDReference(1, .5, .1, 1), 2, 1e-4, id='UD Example'),
+        pytest.param(UDReference(1, .5, .1, 1), 2, 1e-3, id='UD Example'),
     ])
     def test_fixed_ps_fit(self, reference, wMax, tol):
         env = BosonicEnvironment.from_spectral_density(

From 2295f5a4657f46e4ca70dc42cc14323bbdbab412 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Fri, 26 Sep 2025 15:08:37 +0900
Subject: [PATCH 6/6] Update power spectrum fitting test again since it failed
 on MacOS

---
 qutip/tests/core/test_environment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qutip/tests/core/test_environment.py b/qutip/tests/core/test_environment.py
index 2c60a3e45f..34b3d16959 100644
--- a/qutip/tests/core/test_environment.py
+++ b/qutip/tests/core/test_environment.py
@@ -804,7 +804,7 @@ def test_fixed_ps_fit(self, reference, wMax, tol):
             reference.spectral_density, T=reference.T, tag="test"
         )
         wlist = np.linspace(-wMax, wMax, 200)
-        fit, info = env.approximate("ps", wlist,Nmax=6)
+        fit, info = env.approximate("ps", wlist, Nmax=6, maxfev=1e5)
 
         assert isinstance(fit, ExponentialBosonicEnvironment)
         assert fit.T == env.T