Merge pull request #2360 from Ericgig/prepare-qutip-5.0.0

Prepare qutip 5.0.0

.github/workflows/tests.yml

@@ -69,6 +69,7 @@ jobs:
             python-version: "3.10"
             scipy-requirement: ">=1.10,<1.11"
             numpy-requirement: ">=1.24,<1.25"
+            oldcython: 1
             nocython: 1
 
           # Python 3.11 and recent numpy
@@ -128,15 +129,19 @@ jobs:
         # In the run, first we handle any special cases.  We do this in bash
         # rather than in the GitHub Actions file directly, because bash gives us
         # a proper programming language to use.
+        # We install without build isolation so qutip is compiled with the
+        # version of cython, scipy, numpy in the test matrix, not a temporary
+        # version use in the installation virtual environment.
         run: |
-          QUTIP_TARGET="tests,graphics,semidefinite,ipython,extras"
-          if [[ -z "${{ matrix.nocython }}" ]]; then
-            QUTIP_TARGET="$QUTIP_TARGET,runtime_compilation"
-          fi
-          if [[ "${{ matrix.oldcython }}" ]]; then
-            pip install cython==0.29.36
-          fi
-          export CI_QUTIP_WITH_OPENMP=${{ matrix.openmp }}
+          # Install the extra requirement
+          python -m pip install pytest>=5.2 pytest-rerunfailures  # tests
+          python -m pip install matplotlib>=1.2.1  # graphics
+          python -m pip install cvxpy>=1.0 cvxopt  # semidefinite
+          python -m pip install ipython  # ipython
+          python -m pip install loky tqdm  # extras
+          python -m pip install "coverage${{ matrix.coverage-requirement }}" chardet
+          python -m pip install pytest-cov coveralls pytest-fail-slow
+
           if [[ -z "${{ matrix.nomkl }}" ]]; then
             conda install blas=*=mkl "numpy${{ matrix.numpy-requirement }}" "scipy${{ matrix.scipy-requirement }}"
           elif [[ "${{ matrix.os }}" =~ ^windows.*$ ]]; then
@@ -153,9 +158,17 @@ jobs:
             # Use openmpi because mpich causes problems. Note, environment variable names change in v5
             conda install "openmpi<5" mpi4py
           fi
-          python -m pip install -e .[$QUTIP_TARGET]
-          python -m pip install "coverage${{ matrix.coverage-requirement }}"
-          python -m pip install pytest-cov coveralls pytest-fail-slow
+          if [[ "${{ matrix.oldcython }}" ]]; then
+            python -m pip install cython==0.29.36 filelock
+          else
+            python -m pip install cython filelock
+          fi
+
+          python -m pip install -e . -v --no-build-isolation
+
+          if [[ "${{ matrix.nocython }}" ]]; then
+            python -m pip uninstall cython -y
+          fi
 
       - name: Package information
         run: |

qutip/core/cy/qobjevo.pyx

@@ -613,11 +613,11 @@ cdef class QobjEvo:
 
     def __imatmul__(QobjEvo self, other):
         if isinstance(other, (Qobj, QobjEvo)):
-            if self.dims[1] != other.dims[0]:
+            if self._dims[1] != other._dims[0]:
                 raise TypeError("incompatible dimensions" +
                                 str(self.dims[1]) + ", " +
                                 str(other.dims[0]))
-            self._dims = Dimensions([self.dims[0], other.dims[1]])
+            self._dims = Dimensions([self._dims[0], other._dims[1]])
             self.shape = (self.shape[0], other.shape[1])
             if isinstance(other, Qobj):
                 other = _ConstantElement(other)

qutip/core/data/matmul.pyx

@@ -69,8 +69,10 @@ cdef int _check_shape(Data left, Data right, Data out=None) except -1 nogil:
         )
     if (
         out is not None
-        and out.shape[0] != left.shape[0]
-        and out.shape[1] != right.shape[1]
+        and (
+            out.shape[0] != left.shape[0]
+            or out.shape[1] != right.shape[1]
+        )
     ):
         raise ValueError(
             "incompatible output shape, got "

qutip/core/data/reshape.pyx

@@ -2,7 +2,7 @@
 #cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
 
 from libc.string cimport memcpy, memset
-
+from scipy.linalg cimport cython_blas as blas
 cimport cython
 
 import warnings
@@ -52,7 +52,7 @@ cpdef CSR reshape_csr(CSR matrix, idxint n_rows_out, idxint n_cols_out):
     return out
 
 
-cdef inline idxint _reshape_dense_reindex(idxint idx, idxint size):
+cdef inline size_t _reshape_dense_reindex(size_t idx, size_t size):
     return (idx // size) + (idx % size)
 
 
@@ -66,8 +66,9 @@ cpdef Dense reshape_dense(Dense matrix, idxint n_rows_out, idxint n_cols_out):
     out = dense.zeros(n_rows_out, n_cols_out)
     cdef size_t idx_in=0, idx_out=0
     cdef size_t size = n_rows_out * n_cols_out
+    cdef size_t tmp = (<size_t> matrix.shape[1]) * (<size_t> n_rows_out)
     # TODO: improve the algorithm here.
-    cdef size_t stride = _reshape_dense_reindex(matrix.shape[1]*n_rows_out, size)
+    cdef size_t stride = _reshape_dense_reindex(tmp, size)
     for idx_in in range(size):
         out.data[idx_out] = matrix.data[idx_in]
         idx_out = _reshape_dense_reindex(idx_out + stride, size)
@@ -99,7 +100,16 @@ cpdef Dense column_stack_dense(Dense matrix, bint inplace=False):
         return out
     if inplace:
         warnings.warn("cannot stack columns inplace for C-ordered matrix")
-    return reshape_dense(matrix.transpose(), matrix.shape[0]*matrix.shape[1], 1)
+    out = dense.zeros(matrix.shape[0] * matrix.shape[1], 1)
+    cdef idxint col
+    cdef int ONE=1
+    for col in range(matrix.shape[1]):
+        blas.zcopy(
+            &matrix.shape[0],
+            &matrix.data[col], &matrix.shape[1],
+            &out.data[col * matrix.shape[0]], &ONE
+        )
+    return out
 
 
 cpdef Dia column_stack_dia(Dia matrix):

qutip/core/qobj.py

@@ -312,7 +312,7 @@ def __init__(self, arg=None, dims=None,
     def copy(self):
         """Create identical copy"""
         return Qobj(arg=self._data,
-                    dims=self.dims,
+                    dims=self._dims,
                     isherm=self._isherm,
                     isunitary=self._isunitary,
                     copy=True)
@@ -533,6 +533,7 @@ def _str_header(self):
             "Quantum object: dims=" + str(self.dims),
             "shape=" + str(self._data.shape),
             "type=" + repr(self.type),
+            "dtype=" + self.dtype.__name__,
         ])
         if self.type in ('oper', 'super'):
             out += ", isherm=" + str(self.isherm)

qutip/core/states.py

@@ -571,7 +571,7 @@ def projection(dimensions, n, m, offset=None, *, dtype=None):
         Number of basis states in Hilbert space.  If a list, then the resultant
         object will be a tensor product over spaces with those dimensions.
 
-    n, m : float
+    n, m : int
         The number states in the projection.
 
     offset : int, default: 0

qutip/core/superoperator.py

@@ -12,6 +12,7 @@
 from . import data as _data
 from .dimensions import Compound, SuperSpace, Space
 
+
 def _map_over_compound_operators(f):
     """
     Convert a function which takes Qobj into one that can also take compound
@@ -49,7 +50,7 @@ def liouvillian(H=None, c_ops=None, data_only=False, chi=None):
         In some systems it is possible to determine the statistical moments
         (mean, variance, etc) of the probability distributions of occupation of
         various states by numerically evaluating the derivatives of the steady
-        state occupation probability as a function of artificial phase 
+        state occupation probability as a function of artificial phase
         parameters ``chi`` which are included in the
         :func:`lindblad_dissipator` for each collapse operator. See the
         documentation of :func:`lindblad_dissipator` for references and further
@@ -95,15 +96,10 @@ def liouvillian(H=None, c_ops=None, data_only=False, chi=None):
         L += sum(lindblad_dissipator(c_op, chi=chi_)
                  for c_op, chi_ in zip(c_ops, chi))
         return L
-
-    op_dims = H.dims
-    op_shape = H.shape
-    sop_dims = [[op_dims[0], op_dims[0]], [op_dims[1], op_dims[1]]]
-    sop_shape = [np.prod(op_dims), np.prod(op_dims)]
-    spI = _data.identity(op_shape[0], dtype=type(H.data))
-
+    spI = _data.identity_like(H.data)
     data = _data.mul(_data.kron(spI, H.data), -1j)
-    data = _data.add(data, _data.kron_transpose(H.data, spI), scale=1j)
+    data = _data.add(data, _data.kron_transpose(H.data, spI),
+                     scale=1j)
 
     for c_op, chi_ in zip(c_ops, chi):
         c = c_op.data
@@ -117,7 +113,7 @@ def liouvillian(H=None, c_ops=None, data_only=False, chi=None):
         return data
     else:
         return Qobj(data,
-                    dims=sop_dims,
+                    dims=[H._dims, H._dims],
                     superrep='super',
                     copy=False)
 
@@ -316,10 +312,9 @@ def spost(A):
     """
     if not A.isoper:
         raise TypeError('Input is not a quantum operator')
-    Id = _data.identity(A.shape[0], dtype=type(A.data))
     data = _data.kron_transpose(A.data, _data.identity_like(A.data))
     return Qobj(data,
-                dims=[A.dims, A.dims],
+                dims=[A._dims, A._dims],
                 superrep='super',
                 isherm=A._isherm,
                 copy=False)

qutip/simdiag.py

@@ -79,15 +79,15 @@ def simdiag(ops, evals: bool = True, *,
         A = ops[jj]
         shape = A.shape
         if shape[0] != shape[1]:
-            raise TypeError('Matricies must be square.')
+            raise TypeError('Matrices must be square.')
         if shape[0] != N:
             raise TypeError('All matrices. must be the same shape')
         if not A.isherm:
-            raise TypeError('Matricies must be Hermitian')
+            raise TypeError('Matrices must be Hermitian')
         for kk in range(jj):
             B = ops[kk]
             if (A * B - B * A).norm() / (A * B).norm() > tol:
-                raise TypeError('Matricies must commute.')
+                raise TypeError('Matrices must commute.')
 
     # TODO: rewrite using Data object
     eigvals, eigvecs = _data.eigs(ops[0].data, True, True)

qutip/solver/_feedback.py

@@ -133,18 +133,18 @@ def __repr__(self):
         return "CollapseFeedback"
 
 
-def _default_weiner(t):
+def _default_wiener(t):
     return np.zeros(1)
 
 
-class _WeinerFeedback(_Feedback):
-    code = "WeinerFeedback"
+class _WienerFeedback(_Feedback):
+    code = "WienerFeedback"
 
     def __init__(self, default=None):
-        self.default = default or _default_weiner
+        self.default = default or _default_wiener
 
     def check_consistency(self, dims):
         pass
 
     def __repr__(self):
-        return "WeinerFeedback"
+        return "WienerFeedback"

qutip/solver/mcsolve.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 from ..core import QobjEvo, spre, spost, Qobj, unstack_columns
-from .multitraj import MultiTrajSolver, _MTSystem
+from .multitraj import MultiTrajSolver, _MultiTrajRHS
 from .solver_base import Solver, Integrator, _solver_deprecation
 from .result import McResult, McTrajectoryResult, McResultImprovedSampling
 from .mesolve import mesolve, MESolver
@@ -167,27 +167,19 @@ def mcsolve(H, state, tlist, c_ops=(), e_ops=None, ntraj=500, *,
     return result
 
 
-class _MCSystem(_MTSystem):
+class _MCRHS(_MultiTrajRHS):
     """
     Container for the operators of the solver.
     """
 
-    def __init__(self, rhs, c_ops, n_ops):
-        self.rhs = rhs
+    def __init__(self, H, c_ops, n_ops):
+        self.rhs = H
         self.c_ops = c_ops
         self.n_ops = n_ops
-        self._collapse_key = ""
 
     def __call__(self):
         return self.rhs
 
-    def __getattr__(self, attr):
-        if attr == "rhs":
-            raise AttributeError
-        if hasattr(self.rhs, attr):
-            return getattr(self.rhs, attr)
-        raise AttributeError
-
     def arguments(self, args):
         self.rhs.arguments(args)
         for c_op in self.c_ops:
@@ -456,14 +448,16 @@ def __init__(self, H, c_ops, *, options=None):
         self._num_collapse = len(self._c_ops)
         self.options = options
 
-        system = _MCSystem(rhs, self._c_ops, self._n_ops)
+        system = _MCRHS(rhs, self._c_ops, self._n_ops)
         super().__init__(system, options=options)
 
     def _restore_state(self, data, *, copy=True):
         """
         Retore the Qobj state from its data.
         """
-        if self._state_metadata['dims'] == self.rhs.dims[1]:
+        # Duplicated from the Solver class, but removed the check for the
+        # normalize_output option, since MCSolver doesn't have that option.
+        if self._state_metadata['dims'] == self.rhs._dims[1]:
             state = Qobj(unstack_columns(data),
                          **self._state_metadata, copy=False)
         else:
@@ -480,37 +474,20 @@ def _initialize_stats(self):
         })
         return stats
 
-    def _initialize_run_one_traj(self, seed, state, tlist, e_ops,
-                                 no_jump=False, jump_prob_floor=0.0):
-        result = self._trajectory_resultclass(e_ops, self.options)
-        generator = self._get_generator(seed)
-        self._integrator.set_state(tlist[0], state, generator,
-                                   no_jump=no_jump,
-                                   jump_prob_floor=jump_prob_floor)
-        result.add(tlist[0], self._restore_state(state, copy=False))
-        return result
-
-    def _run_one_traj(self, seed, state, tlist, e_ops, no_jump=False,
-                      jump_prob_floor=0.0):
+    def _run_one_traj(self, seed, state, tlist, e_ops, **integrator_kwargs):
         """
         Run one trajectory and return the result.
         """
-        result = self._initialize_run_one_traj(seed, state, tlist, e_ops,
-                                               no_jump=no_jump,
-                                               jump_prob_floor=jump_prob_floor)
-        seed, result = self._integrate_one_traj(seed, tlist, result)
+        seed, result = super()._run_one_traj(seed, state, tlist, e_ops,
+                                             **integrator_kwargs)
         result.collapse = self._integrator.collapses
         return seed, result
 
     def run(self, state, tlist, ntraj=1, *,
             args=None, e_ops=(), timeout=None, target_tol=None, seeds=None):
-        """
-        Do the evolution of the Quantum system.
-        See the overridden method for further details. The modification
-        here is to sample the no-jump trajectory first. Then, the no-jump
-        probability is used as a lower-bound for random numbers in future
-        monte carlo runs
-        """
+        # Overridden to sample the no-jump trajectory first. Then, the no-jump
+        # probability is used as a lower-bound for random numbers in future
+        # monte carlo runs
         if not self.options.get("improved_sampling", False):
             return super().run(state, tlist, ntraj=ntraj, args=args,
                                e_ops=e_ops, timeout=timeout,
@@ -540,7 +517,8 @@ def run(self, state, tlist, ntraj=1, *,
         start_time = time()
         map_func(
             self._run_one_traj, seeds[1:],
-            (state0, tlist, e_ops, False, no_jump_prob),
+            task_args=(state0, tlist, e_ops),
+            task_kwargs={'no_jump': False, 'jump_prob_floor': no_jump_prob},
             reduce_func=result.add, map_kw=map_kw,
             progress_bar=self.options["progress_bar"],
             progress_bar_kwargs=self.options["progress_kwargs"]
@@ -557,9 +535,9 @@ def _get_integrator(self):
             integrator = method
         else:
             raise ValueError("Integrator method not supported.")
-        integrator_instance = integrator(self.system(), self.options)
+        integrator_instance = integrator(self.rhs(), self.options)
         mc_integrator = self._mc_integrator_class(
-            integrator_instance, self.system, self.options
+            integrator_instance, self.rhs, self.options
         )
         self._init_integrator_time = time() - _time_start
         return mc_integrator