Qutip 4.4.1 (#1069)

* simdim debug

* removed mutable in mcsolve

* Merge pull request #1041 from Ericgig/CS_inter

QobjEvo do not need to start from 0 anymore

* Merge pull request #1049 from goerz/sesolvemsg

Don't print msg when mesolve falls back to sesolve

* Purity (#1045)

* updated readme installation version to latest

* fix tr() docstring

* add purity as Qobj method

* add tests for purity

* fix typo

* fix tests

* update purity docstring

* update test purity docs

* add PEP8 spaces

* all file limit to 80 characters

* fix typo err

* versions number and more lenient codeclimat

* picking bug on windows

* Compatibility of constant QobjEvo (#1064)

* Support cte Qobjevo addition

* strill allow Q_objec=None

* correct type for Q_object=[Qobj, Qobj]

* Time-dependent coefficients without cubic interpolation in QobjEvo (#1050)

* non-compiled case with interp1d

* step function coefficients in compiled code

* remove temporary functions

* Merge similar classes StepCoeff

* improve case identification

* allow tlist not from 0

* fix bug and add docs

* add tests for coeff call

* clearer error type

* delete redundant line and improve description

* Revised random ket generators for dims (#1058)

* Revised random ket generators for dims

Revised rand_ket, rand_ket_haar to allow specifying the ket dimensions.

Previously _check_ket_dims was broken, and definition of dims in the two functions was conflicting.

Changes here:
- removed _check_ket_dims, using the more general _check_dims instead.
- Changed the input method for the rand_ket functions to receive either int, as before, _or_ a dims specifications, rather than both. 
- If dims is specified, N is inferred from dims.

Since both functions were broken previously if dims was specified, and this is backwards compatible if only N is specified, it should be entirely backwards compatible.

* Update random_objects.py

* Revised rand_ket and rand_ket_haar

Restored functional shape to previous form.

* Update random_objects.py

* Update random_objects.py

* elif without condition

* Bugfix: permute and correlator (#1060)

* permute fix and shifted QobjEvo

* time shift in QobjEvo and permute bugfix

* style improvement

* more pep8

* Remove mutable default values in mesolve/sesolve (#1062)

This closes #1061

* some pep8

.codeclimate.yml

@@ -2,21 +2,21 @@ version: "2"
 checks:
   argument-count:
     config:
-      threshold: 6
+      threshold: 8
   complex-logic:
     config:
-      threshold: 5
+      threshold: 10
   file-lines:
     enabled: false
   method-complexity:
     config:
-      threshold: 5
+      threshold: 10
   method-count:
     config:
       threshold: 50
   method-lines:
     config:
-      threshold: 50
+      threshold: 70
   nested-control-flow:
     config:
       threshold: 4

qutip/correlation.py

@@ -51,6 +51,7 @@
 from qutip.mcsolve import mcsolve
 from qutip.operators import qeye
 from qutip.qobj import Qobj, isket, issuper
+from qutip.qobjevo import QobjEvo
 from qutip.rhs_generate import rhs_clear, _td_wrap_array_str
 from qutip.cy.utilities import _cython_build_cleanup
 from qutip.settings import debug
@@ -576,12 +577,12 @@ def spectrum_correlation_fft(tlist, y, inverse=False):
     dt = tlist[1] - tlist[0]
     if not np.allclose(np.diff(tlist), dt*np.ones(N-1,dtype=float)):
         raise Exception('tlist must be equally spaced for FFT.')
-    
+
     if inverse:
            F = N * scipy.fftpack.ifft(y)
     else:
            F = scipy.fftpack.fft(y)
-    
+
     # calculate the frequencies for the components in F
     f = scipy.fftpack.fftfreq(N, dt)
 
@@ -813,7 +814,7 @@ def correlation_4op_1t(H, state0, taulist, c_ops, a_op, b_op, c_op, d_op,
     References
     ----------
     See, Gardiner, Quantum Noise, Section 5.2.
-                       
+
     .. note:: Deprecated in QuTiP 3.1
               Use correlation_3op_1t() instead.
 
@@ -1100,7 +1101,7 @@ def _correlation_me_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
     rho_t = mesolve(H, rho0, tlist, c_ops, [],
                     args=args, options=options).states
     corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
-    H_shifted, c_ops_shifted, _args = _transform_L_t_shift(H, c_ops, args)
+    H_shifted, c_ops_shifted, _args = _transform_L_t_shift_new(H, c_ops, args)
     if config.tdname:
         _cython_build_cleanup(config.tdname)
     rhs_clear()
@@ -1229,7 +1230,7 @@ def _correlation_mc_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
     ).states
 
     corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
-    H_shifted, c_ops_shifted, _args = _transform_L_t_shift(H, c_ops, args)
+    H_shifted, c_ops_shifted, _args = _transform_L_t_shift_new(H, c_ops, args)
     if config.tdname:
         _cython_build_cleanup(config.tdname)
     rhs_clear()
@@ -1291,10 +1292,10 @@ def _correlation_mc_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
                     dtype=corr_mat.dtype
                 )
                 corr_mat[t_idx, :] += corr_mat_add
-                    
+
         if t_idx == 1:
             options.rhs_reuse = True
-    
+
     if config.tdname:
         _cython_build_cleanup(config.tdname)
     rhs_clear()
@@ -1303,7 +1304,6 @@ def _correlation_mc_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
 
 
 # pseudo-inverse solvers
-
 def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
     """
     Internal function for calculating the spectrum of the correlation function
@@ -1344,126 +1344,52 @@ def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
 
 
 # auxiliary
-
-def _transform_L_t_shift(H, c_ops, args={}):
-    """
-    Time shift the Hamiltonian with private time-shift variable _t0
-    """
-
-    # while this list doesn't seem exhaustive, mesolve has already been called
-    # and has hence called _td_format_check from qutip.rhs_generate. important
-    # to keep in mind is that mesolve already requires the types of
-    # time-dependence to be the same for the hamiltonian as for the collapse
-    # operators.
-
-
-    if isinstance(H, Qobj):
-        # constant hamiltonian
-        H_shifted = H  # not shifted!
-
-    c_ops_is_td = False
-    if isinstance(c_ops, list):
-        for i in range(len(c_ops)):
-            # test is collapse operators are time-dependent
-            if isinstance(c_ops[i], list):
-                c_ops_is_td = True
-
-    if not c_ops_is_td:
-        # constant collapse operators
-        c_ops_shifted = c_ops  # not shifted!
-
-    if isinstance(H, Qobj) and not c_ops_is_td:
-        # constant hamiltonian and collapse operators
-        _args = args  # not shifted!
-
-    if isinstance(H, types.FunctionType):
-        # function-callback based time-dependence
-        if isinstance(args, dict) or args is None:
-            if args is None:
-                _args = {"_t0": 0}
-            else:
-                _args = args.copy()
-                _args["_t0"] = 0
-            H_shifted = lambda t, args_i: H(t+args_i["_t0"], args_i)
+def _transform_shift_one_coeff(op, args):
+    if isinstance(op, types.FunctionType):
+        # function-list based time-dependence
+        if isinstance(args, dict):
+            def fn(t, args_i):
+                return op(t + args_i["_t0"], args_i)
+            fn = lambda t, args_i: \
+                op(t + args_i["_t0"], args_i)
         else:
-            raise TypeError("If using function-callback based Hamiltonian" +
-                            "time-dependence, args must be a dictionary")
-
-    if isinstance(H, list) or c_ops_is_td:
-        # string/function-list based time-dependence
-        if args is None:
-            _args = {"_t0": 0}
-        elif isinstance(args, dict):
-            _args = args.copy()
-            _args["_t0"] = 0
-        else:
-            _args = {"_user_args": args, "_t0": 0}
-
-        if isinstance(H, list):
-            # hamiltonian is time-dependent
-            H_shifted = []
-
-            for i in range(len(H)):
-                if isinstance(H[i], list):
-                    # modify Hamiltonian time dependence in accordance with the
-                    # quantum regression theorem
-                    if isinstance(args, dict) or args is None:
-                        if isinstance(H[i][1], types.FunctionType):
-                            # function-list based time-dependence
-                            fn = lambda t, args_i: \
-                                H[i][1](t + args_i["_t0"], args_i)
-                        else:
-                            # string-list based time-dependence
-                            # Again, note: _td_format_check already raises
-                            # errors formixed td formatting
-                            fn = sub("(?<=[^0-9a-zA-Z_])t(?=[^0-9a-zA-Z_])",
-                                     "(t+_t0)", H[i][1])
-                    else:
-                        if isinstance(H[i][1], types.FunctionType):
-                            # function-list based time-dependence
-                            fn = lambda t, args_i: \
-                                H[i][1](t + args_i["_t0"],
-                                        args_i["_user_args"])
-                        else:
-                            raise TypeError("If using string-list based" +
-                                            "Hamiltonian time-dependence, " +
-                                            "args must be a dictionary")
-                    H_shifted.append([H[i][0], fn])
-                else:
-                    H_shifted.append(H[i])
-
-        if c_ops_is_td:
-            # collapse operators are time-dependent
-            c_ops_shifted = []
-
-            for i in range(len(c_ops)):
-                if isinstance(c_ops[i], list):
-                    # modify collapse operators time dependence in accordance
-                    # with the quantum regression theorem
-                    if isinstance(args, dict) or args is None:
-                        if isinstance(c_ops[i][1], types.FunctionType):
-                            # function-list based time-dependence
-                            fn = lambda t, args_i: \
-                                c_ops[i][1](t + args_i["_t0"], args_i)
-                        else:
-                            # string-list based time-dependence
-                            # Again, note: _td_format_check already raises
-                            # errors formixed td formatting
-                            fn = sub("(?<=[^0-9a-zA-Z_])t(?=[^0-9a-zA-Z_])",
-                                     "(t+_t0)", c_ops[i][1])
-                    else:
-                        if isinstance(H[i][1], types.FunctionType):
-                            # function-list based time-dependence
-                            fn = lambda t, args_i: \
-                                c_ops[i][1](t + args_i["_t0"],
-                                            args_i["_user_args"])
-                        else:
-                            raise TypeError("If using string-list based" +
-                                            "collapse operator" +
-                                            "time-dependence, " +
-                                            "args must be a dictionary")
-                    c_ops_shifted.append([c_ops[i][0], fn])
-                else:
-                    c_ops_shifted.append(c_ops[i])
+            def fn(t, args_i):
+                return op(t + args_i["_t0"], args_i["_user_args"])
+    else:
+        fn = sub("(?<=[^0-9a-zA-Z_])t(?=[^0-9a-zA-Z_])",
+                 "(t+_t0)", " " + op + " ")
+    return fn
+
+
+def _transform_shift_one_op(op, args={}):
+    if isinstance(op, Qobj):
+        new_op = op
+    elif isinstance(op, QobjEvo):
+        new_op = op
+        new_op._shift
+    elif callable(op):
+        def new_op(t, args_i):
+            return op(t + args_i["_t0"], args_i)
+    elif isinstance(op, list):
+        new_op = []
+        for block in op:
+            if isinstance(block, list):
+                new_op.append([block[0],
+                               _transform_shift_one_coeff(block[1], args)])
+            else:
+                new_op.append(block)
+    return new_op
+
+
+def _transform_L_t_shift_new(H, c_ops, args={}):
+    H_shifted = _transform_shift_one_op(H, args)
+    c_ops_shifted = [_transform_shift_one_op(op, args) for op in c_ops]
+    if args is None:
+        _args = {"_t0": 0}
+    elif isinstance(args, dict):
+        _args = args.copy()
+        _args["_t0"] = 0
+    else:
+        _args = {"_user_args": args, "_t0": 0}
 
     return H_shifted, c_ops_shifted, _args

qutip/cy/cqobjevo_factor.pyx

@@ -39,7 +39,8 @@ cimport cython
 cimport libc.math
 from qutip.cy.inter import _prep_cubic_spline
 from qutip.cy.inter cimport (_spline_complex_cte_second,
-                             _spline_complex_t_second)
+                             _spline_complex_t_second,
+                             _step_complex_t, _step_complex_cte)
 from qutip.cy.interpolate cimport (interp, zinterp)
 from qutip.cy.cqobjevo cimport CQobjEvo
 
@@ -226,6 +227,50 @@ cdef class InterCoeffT(CoeffFunc):
         self.M = state[5]
 
 
+cdef class StepCoeff(CoeffFunc):
+    cdef int n_t
+    cdef double[::1] tlist
+    cdef complex[:,::1] y
+
+    def __init__(self, ops, args, tlist):
+        cdef int i, j
+        self._args = {}
+        self._num_ops = len(ops)
+        self.tlist = tlist
+        self.n_t = len(tlist)
+        self.y = np.zeros((self._num_ops, self.n_t), dtype=complex)
+        for i in range(self._num_ops):
+            for j in range(self.n_t):
+                self.y[i,j] = ops[i][2][j]
+
+    def set_arg(self, args):
+        pass
+
+    def __getstate__(self):
+        return (self._num_ops, self.n_t, None, np.array(self.tlist),
+                np.array(self.y))
+
+    def __setstate__(self, state):
+        self._num_ops = state[0]
+        self.n_t = state[1]
+        self.tlist = state[3]
+        self.y = state[4]
+
+
+cdef class StepCoeffT(StepCoeff):
+    cdef void _call_core(self, double t, complex* coeff):
+        cdef int i
+        for i in range(self._num_ops):
+            coeff[i] = _step_complex_t(t, self.tlist, self.y[i, :], self.n_t)
+
+
+cdef class StepCoeffCte(StepCoeff): 
+    cdef void _call_core(self, double t, complex* coeff):
+        cdef int i
+        for i in range(self._num_ops):
+            coeff[i] = _step_complex_cte(t, self.tlist, self.y[i, :], self.n_t)
+
+
 cdef class StrCoeff(CoeffFunc):
     def __init__(self, ops, args, tlist, dyn_args=[]):
         self._num_ops = len(ops)

qutip/cy/inter.pxd

@@ -48,3 +48,11 @@ cdef double _spline_float_t_second(double x, double[::1] t,
 cdef double _spline_float_cte_second(double x, double[::1] t,
                                      double[::1] y, double[::1] M,
                                      int N, double dt)
+
+cdef double _step_float_cte(double x, double[::1] t, double[::1] y, int n_t)
+
+cdef complex _step_complex_cte(double x, double[::1] t, complex[::1] y, int n_t)
+
+cdef double _step_float_t(double x, double[::1] t, double[::1] y, int n_t)
+
+cdef complex _step_complex_t(double x, double[::1] t, complex[::1] y, int n_t)