Fix setting of sso.m_ops in heterodyne smesolver and passing through of sc_ops to photocurrent solver.

doc/changes/2081.bugfix

@@ -0,0 +1,2 @@
+Fix setting of sso.m_ops in heterodyne smesolver and passing through of sc_ops to photocurrent solver. (`#2081 <https://github.com/qutip/qutip/pull/2081>`_ by Bogdan Reznychenko and Simon Cross)
+Update calls to SciPy eigvalsh and eigsh to pass the range of eigenvalues to return using ``subset_by_index=`` instead of ``eigvals=``. ``eigvals=`` was deprecated in SciPy 1.11 and will be removed in SciPy 1.12. (`#2081 <https://github.com/qutip/qutip/pull/2081>`_ by Simon Cross)

qutip/sparse.py

@@ -28,12 +28,13 @@ def _orthogonalize(vec, other):
         norm = np.sum(np.conj(vec) * vec)**0.5
         vec /= norm
 
-    def eigh(mat, eigvals=[]):
+    def eigh(mat, subset_by_index=()):
         val, vec = la.eig(mat)
         val = np.real(val)
         idx = np.argsort(val)
         val = val[idx]
         vec = vec[:, idx]
+        eigvals = subset_by_index
         if eigvals:
             val = val[eigvals[0]:eigvals[1]+1]
             vec = vec[:, eigvals[0]:eigvals[1]+1]
@@ -47,9 +48,10 @@ def eigh(mat, eigvals=[]):
                 same_eigv = 0
         return val, vec
 
-    def eigvalsh(a, UPLO="L", eigvals=[]):
+    def eigvalsh(a, UPLO="L", subset_by_index=()):
         val = la.eigvals(a)
         val = np.sort(np.real(val))
+        eigvals = subset_by_index
         if eigvals:
             return val[eigvals[0]:eigvals[1]+1]
         return val
@@ -177,10 +179,10 @@ def _dense_eigs(data, isherm, vecs, N, eigvals, num_large, num_small):
             else:
                 if num_small > 0:
                     evals, evecs = eigh(
-                        data, eigvals=[0, num_small - 1])
+                        data, subset_by_index=[0, num_small - 1])
                 if num_large > 0:
                     evals, evecs = eigh(
-                        data, eigvals=[N - num_large, N - 1])
+                        data, subset_by_index=[N - num_large, N - 1])
         else:
             evals, evecs = la.eig(data)
     else:
@@ -189,9 +191,9 @@ def _dense_eigs(data, isherm, vecs, N, eigvals, num_large, num_small):
                 evals = eigvalsh(data)
             else:
                 if num_small > 0:
-                    evals = eigvalsh(data, eigvals=[0, num_small - 1])
+                    evals = eigvalsh(data, subset_by_index=[0, num_small - 1])
                 if num_large > 0:
-                    evals = eigvalsh(data, eigvals=[N - num_large, N - 1])
+                    evals = eigvalsh(data, subset_by_index=[N - num_large, N - 1])
         else:
             evals = la.eigvals(data)
 

qutip/stochastic.py

@@ -597,13 +597,20 @@ def smesolve(H, rho0, times, c_ops=[], sc_ops=[], e_ops=[],
     elif sso.method == 'heterodyne':
         if sso.m_ops is None:
             m_ops = []
+        elif len(sso.m_ops) != 2 * len(sc_ops):
+            raise ValueError(
+                "When using the heterodyne method there should be two"
+                " measurement operators (m_ops) for each collapse operator"
+                " (sc_ops)."
+            )
         sso.sops = []
         for c in sso.sc_ops:
             if sso.m_ops is None:
                 m_ops += [c + c.dag(), -1j * (c - c.dag())]
             sso.sops += [(spre(c) + spost(c.dag())) / np.sqrt(2),
                          (spre(c) - spost(c.dag())) * -1j / np.sqrt(2)]
-        sso.m_ops = m_ops
+        if sso.m_ops is None:
+            sso.m_ops = m_ops
         if not isinstance(sso.dW_factors, list):
             sso.dW_factors = [np.sqrt(2)] * len(sso.sops)
         elif len(sso.dW_factors) == len(sso.m_ops):
@@ -689,7 +696,7 @@ def ssesolve(H, psi0, times, sc_ops=[], e_ops=[],
     if "method" in kwargs and kwargs["method"] == "photocurrent":
         print("stochastic solver with photocurrent method has been moved to "
               "it's own function: photocurrent_sesolve")
-        return photocurrent_sesolve(H, psi0, times, c_ops=c_ops,
+        return photocurrent_sesolve(H, psi0, times, sc_ops=sc_ops,
                                     e_ops=e_ops, _safe_mode=_safe_mode,
                                     args=args, **kwargs)
 

qutip/tests/test_stochastic_me.py

@@ -4,7 +4,8 @@
 
 from qutip import (smesolve, mesolve, photocurrent_mesolve, liouvillian,
                    QobjEvo, spre, spost, destroy, coherent, parallel_map,
-                   qeye, fock_dm, general_stochastic, ket2dm, num)
+                   qeye, fock_dm, general_stochastic, ket2dm, num,
+                   basis, sigmax, sigmay, sigmaz, sigmam)
 
 def f(t, args):
     return args["a"] * t
@@ -296,3 +297,66 @@ def test_smesolve_bad_e_ops():
         res = smesolve(H, psi0, times, sc_ops=sc_ops, e_ops=e_ops, noise=1,
                        ntraj=ntraj, nsubsteps=nsubsteps, method='homodyne',
                        map_func=parallel_map)
+
+
+def test_heterodyne_smesolve_custom_m_ops():
+    b = 1  # drive amplitude
+    gamma = 1  # spont.  emission rate
+    eta = 0.3  # coupling efficiency
+    n_steps = 1000
+    n_traj = 300
+    noise_seed = 0  # noise random seed
+
+    H = np.sqrt(eta * gamma) * b * sigmay()
+    c_ops = [np.sqrt(gamma) * sigmam()]
+    psi0 = basis(2, 0)
+    times = np.linspace(0, 2 * np.pi, n_steps)
+
+    sme_het = smesolve(
+        H,
+        psi0,
+        times,
+        [],
+        c_ops,
+        e_ops=[sigmax(), sigmay(), sigmaz()],
+        store_measurement=True,
+        dW_factors=[1e-5, 1e-5],
+        method="heterodyne",
+        m_ops=[np.sqrt(eta)*sigmax(), np.sqrt(eta)*sigmay()],
+        ntraj=n_traj,
+        noise=noise_seed,
+    )
+
+    assert np.array(sme_het.measurement).shape == (n_traj, n_steps, 1, 2)
+    np.testing.assert_allclose(
+        np.array(sme_het.measurement).mean(axis=0)[:, 0, 0].T,
+        np.sqrt(eta) * sme_het.expect[0],
+        atol=1e-2,
+    )
+
+
+def test_heterodyne_mesolve_incorrect_custom_m_ops():
+    eta = 0.3
+
+    with pytest.raises(ValueError) as err:
+        smesolve(
+            sigmax(),
+            basis(2),
+            np.linspace(0, 1, 10),
+            [],
+            [sigmam()],
+            e_ops=[],
+            store_measurement=True,
+            method="heterodyne",
+            # three m_ops, which is incorrect:
+            m_ops=[
+                np.sqrt(eta) * sigmax(),
+                np.sqrt(eta) * sigmay(),
+                np.sqrt(eta) * sigmaz(),
+            ],
+            ntraj=10,
+        )
+    assert str(err.value) == (
+        "When using the heterodyne method there should be two measurement"
+        " operators (m_ops) for each collapse operator (sc_ops)."
+    )