Merge pull request #1478 from AGaliciaMartinez/issue#1460

Fix imaginary values in correlation functions being discarded.

qutip/correlation.py

@@ -116,7 +116,7 @@ def correlation_2op_1t(H, state0, taulist, c_ops, a_op, b_op,
     Returns
     -------
     corr_vec : ndarray
-        An array of correlation values for the times specified by `tlist`.
+        An array of correlation values for the times specified by `taulist`.
 
     References
     ----------

qutip/mesolve.py

@@ -510,22 +510,25 @@ def get_curr_state_data(r):
 
         if opt.store_states or expt_callback:
             cdata = get_curr_state_data(r)
+            fdata = dense2D_to_fastcsr_fmode(cdata, size, size)
 
-        if opt.store_states:
-            if issuper(rho0):
-                fdata = dense2D_to_fastcsr_fmode(cdata, size, size)
-                output.states.append(Qobj(fdata, dims=dims))
+            # Try to guess if there is a fast path for rho_t
+            if issuper(rho0) or not rho0.isherm:
+                rho_t = Qobj(fdata, dims=dims)
             else:
-                fdata = dense2D_to_fastcsr_fmode(cdata, size, size)
-                output.states.append(Qobj(fdata, dims=dims, fast="mc-dm"))
+                rho_t = Qobj(fdata, dims=dims, fast="mc-dm")
+
+        if opt.store_states:
+            output.states.append(rho_t)
 
         if expt_callback:
             # use callback method
-            output.expect.append(e_ops(t, Qobj(cdata, dims=dims)))
+            output.expect.append(e_ops(t, rho_t))
 
         for m in range(n_expt_op):
             output.expect[m][t_idx] = expect_rho_vec(e_ops_data[m], r.y,
-                                                     e_ops[m].isherm)
+                                                     e_ops[m].isherm
+                                                     and rho0.isherm)
 
         if t_idx < n_tsteps - 1:
             r.integrate(r.t + dt[t_idx])
@@ -534,6 +537,7 @@ def get_curr_state_data(r):
 
     if opt.store_final_state:
         cdata = get_curr_state_data(r)
-        output.final_state = Qobj(cdata, dims=dims, isherm=True)
+        output.final_state = Qobj(cdata, dims=dims,
+                                  isherm=rho0.isherm or None)
 
     return output

qutip/tests/test_correlation.py

@@ -33,6 +33,7 @@
 
 import pytest
 import functools
+from itertools import product
 import numpy as np
 import qutip
 
@@ -58,6 +59,7 @@ def test_correlation_solver_equivalence(solver, start, legacy):
     system.
     """
     a = qutip.destroy(_equivalence_dimension)
+    x = ( a +a.dag() )/np.sqrt(2)
     H = a.dag() * a
     G1 = 0.75
     n_th = 2
@@ -251,3 +253,56 @@ def test_hamiltonian_order_unimportant():
     backwards = qutip.correlation_2op_2t(H[::-1], start, times, times, [sp],
                                          sp.dag(), sp)
     np.testing.assert_allclose(forwards, backwards, atol=1e-6)
+
+
+@pytest.mark.parametrize(['solver', 'state'], [
+    pytest.param('me', _equivalence_fock, id="me-ket"),
+    pytest.param('me', _equivalence_coherent, id="me-dm"),
+    pytest.param('me', None, id="me-steady"),
+    pytest.param('es', _equivalence_fock, id="es-ket"),
+    pytest.param('es', _equivalence_coherent, id="es-dm"),
+    pytest.param('es', None, id="es-steady"),
+    pytest.param('mc', _equivalence_fock, id="mc-ket",
+                 marks=[pytest.mark.slow]),
+])
+@pytest.mark.parametrize("is_e_op_hermitian", [True, False],
+                         ids=["hermitian", "nonhermitian"])
+@pytest.mark.parametrize("w", [1, 2])
+@pytest.mark.parametrize("gamma", [1, 10])
+def test_correlation_2op_1t_known_cases(solver,
+                                        state,
+                                        is_e_op_hermitian,
+                                        w,
+                                        gamma,
+                                       ):
+    """This test compares the output correlation_2op_1 solution to an analytical
+    solution."""
+
+    a = qutip.destroy(_equivalence_dimension)
+    x = (a + a.dag())/np.sqrt(2)
+
+    H = w * a.dag() * a
+
+    a_op = x if is_e_op_hermitian else a
+    b_op = x if is_e_op_hermitian else a.dag()
+    c_ops = [np.sqrt(gamma) * a]
+
+    times = np.linspace(0, 1, 30)
+
+    # Handle the case state==None when computing expt values
+    rho0 = state if state else qutip.steadystate(H, c_ops)
+    if is_e_op_hermitian:
+        # Analitycal solution for x,x as operators.
+        base = 0
+        base += qutip.expect(a*x, rho0)*np.exp(-1j*w*times - gamma*times/2)
+        base += qutip.expect(a.dag()*x, rho0)*np.exp(1j*w*times - gamma*times/2)
+        base /= np.sqrt(2)
+    else:
+        # Analitycal solution for a,adag as operators.
+        base = qutip.expect(a*a.dag(), rho0)*np.exp(-1j*w*times - gamma*times/2)
+
+    cmp = qutip.correlation_2op_1t(H, state, times, c_ops, a_op, b_op, solver=solver)
+
+    np.testing.assert_allclose(base, cmp, atol=0.25 if solver == 'mc' else 2e-5)
+
+

qutip/tests/test_mesolve.py

@@ -920,5 +920,55 @@ def f(t, args):
             msg="evolution with feedback not proceding as expected")
 
 
+def test_non_hermitian_dm():
+    """Test that mesolve works correctly for density matrices that are
+    not Hermitian.
+    See Issue #1460
+    """
+    N = 2
+    a = destroy(N)
+    x = (a + a.dag())/np.sqrt(2)
+    H = a.dag() * a
+
+    # Create non-Hermitian initial state.
+    rho0 = x*fock_dm(N, 0)
+
+    tlist = np.linspace(0, 0.1, 2)
+
+    options = Options()
+    options.store_final_state = True
+    options.store_states = True
+
+    result = mesolve(H, rho0, tlist, e_ops=[x], options=options)
+
+    msg = ('Mesolve is not working properly with a non Hermitian density' +
+       ' matrix as input. Check computation of '
+      )
+
+    imag_part = np.abs(np.imag(result.expect[0][-1]))
+    # Since we used an initial state that is not Hermitian, the expectation of
+    # x must be imaginary for t>0.
+    assert_(imag_part > 0,
+            msg + "expectation values. They should be imaginary")
+
+    # Check that the output state is not hermitian since the input was not
+    # Hermitian either.
+    assert_(not result.final_state.isherm,
+            msg + " final density  matrix. It should not be hermitian")
+    assert_(not result.states[-1].isherm,
+            msg + " states. They should not be hermitian.")
+
+    # Check that when suing a callable we get imaginary expectation values.
+    def callable_x(t, rho):
+        "Dummy callable_x expectation operator."
+        return expect(rho, x)
+    result = mesolve(H, rho0, tlist, e_ops=callable_x)
+
+    imag_part = np.abs(np.imag(result.expect[-1]))
+    assert_(imag_part > 0,
+            msg + "expectation values when using callable operator." +
+            "They should be imaginary.")
+
+
 if __name__ == "__main__":
     run_module_suite()