Merge branch 'master' into qutip-5.0.X

doc/changes/2352.feature

@@ -0,0 +1 @@
+Add dtype to printed ouput of qobj

doc/changes/2354.misc

@@ -0,0 +1 @@
+Allow scipy 1.12 to be used with qutip.

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/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

@@ -463,7 +463,7 @@ def _restore_state(self, data, *, copy=True):
         """
         Retore the Qobj state from its data.
         """
-        if self._state_metadata['dims'] == self.rhs.dims[1]:
+        if self._state_metadata['dims'] == self.rhs._dims[1]:
             state = Qobj(unstack_columns(data),
                          **self._state_metadata, copy=False)
         else:

qutip/solver/nonmarkov/transfertensor.py

@@ -71,7 +71,7 @@ def ttmsolve(dynmaps, state0, times, e_ops=(), num_learning=0, options=None):
     opt = {
         "store_final_state": False,
         "store_states": None,
-        "normalize_output": "ket",
+        "normalize_output": True,
         "threshold": 0.0,
         "num_learning": 0,
     }

qutip/solver/result.py

@@ -526,7 +526,7 @@ def __init__(
         self._post_init(**kw)
 
     @property
-    def _store_average_density_matricies(self) -> bool:
+    def _store_average_density_matrices(self) -> bool:
         return (
             self.options["store_states"]
             or (self.options["store_states"] is None and self._raw_ops == {})
@@ -536,7 +536,7 @@ def _store_average_density_matricies(self) -> bool:
     def _store_final_density_matrix(self) -> bool:
         return (
             self.options["store_final_state"]
-            and not self._store_average_density_matricies
+            and not self._store_average_density_matrices
             and not self.options["keep_runs_results"]
         )
 
@@ -552,7 +552,7 @@ def _add_first_traj(self, trajectory):
         """
         self.times = trajectory.times
 
-        if trajectory.states and self._store_average_density_matricies:
+        if trajectory.states and self._store_average_density_matrices:
             self._sum_states = [
                 qzero_like(self._to_dm(state)) for state in trajectory.states
             ]
@@ -668,7 +668,7 @@ def _post_init(self):
         store_trajectory = self.options["keep_runs_results"]
         if store_trajectory:
             self.add_processor(self._store_trajectory)
-        if self._store_average_density_matricies:
+        if self._store_average_density_matrices:
             self.add_processor(self._reduce_states)
         if self._store_final_density_matrix:
             self.add_processor(self._reduce_final_state)
@@ -1131,7 +1131,7 @@ def _average_computer(self):
 
     def _add_first_traj(self, trajectory):
         super()._add_first_traj(trajectory)
-        if trajectory.states and self._store_average_density_matricies:
+        if trajectory.states and self._store_average_density_matrices:
             del self._sum_states
             self._sum_states_no_jump = [
                 qzero_like(self._to_dm(state)) for state in trajectory.states

qutip/solver/solver_base.py

@@ -38,7 +38,7 @@ class Solver:
         "progress_kwargs": {"chunk_size": 10},
         "store_final_state": False,
         "store_states": None,
-        "normalize_output": "ket",
+        "normalize_output": True,
         "method": "adams",
     }
     _resultclass = Result
@@ -84,7 +84,7 @@ def _prepare_state(self, state):
                             f" and {state.dims}")
 
         self._state_metadata = {
-            'dims': state.dims,
+            'dims': state._dims,
             'isherm': state.isherm and not (self.rhs.dims == state.dims)
         }
         if self.rhs.dims[1] == state.dims:
@@ -95,7 +95,7 @@ def _restore_state(self, data, *, copy=True):
         """
         Retore the Qobj state from its data.
         """
-        if self._state_metadata['dims'] == self.rhs.dims[1]:
+        if self._state_metadata['dims'] == self.rhs._dims[1]:
             state = Qobj(unstack_columns(data),
                          **self._state_metadata, copy=False)
         else:

qutip/solver/stochastic.py

@@ -4,10 +4,11 @@
 from .result import MultiTrajResult, Result, ExpectOp
 from .multitraj import MultiTrajSolver
 from .. import Qobj, QobjEvo
+from ..core.dimensions import Dimensions
 import numpy as np
 from functools import partial
 from .solver_base import _solver_deprecation
-from ._feedback import _QobjFeedback, _DataFeedback, _WeinerFeedback
+from ._feedback import _QobjFeedback, _DataFeedback, _WienerFeedback
 
 
 class StochasticTrajResult(Result):
@@ -210,8 +211,10 @@ def __init__(self, issuper, H, sc_ops, c_ops, heterodyne):
 
         if self.issuper and not self.H.issuper:
             self.dims = [self.H.dims, self.H.dims]
+            self._dims = Dimensions([self.H._dims, self.H._dims])
         else:
             self.dims = self.H.dims
+            self._dims = self.H._dims
 
     def __call__(self, options):
         if self.issuper:
@@ -229,14 +232,14 @@ def arguments(self, args):
             sc_op.arguments(args)
 
     def _register_feedback(self, val):
-        self.H._register_feedback({"wiener_process": val}, "stochatic solver")
+        self.H._register_feedback({"wiener_process": val}, "stochastic solver")
         for c_op in self.c_ops:
             c_op._register_feedback(
-                {"WeinerFeedback": val}, "stochatic solver"
+                {"WienerFeedback": val}, "stochastic solver"
             )
         for sc_op in self.sc_ops:
             sc_op._register_feedback(
-                {"WeinerFeedback": val}, "stochatic solver"
+                {"WienerFeedback": val}, "stochastic solver"
             )
 
 
@@ -710,13 +713,13 @@ def options(self, new_options):
         MultiTrajSolver.options.fset(self, new_options)
 
     @classmethod
-    def WeinerFeedback(cls, default=None):
+    def WienerFeedback(cls, default=None):
         """
-        Weiner function of the trajectory argument for time dependent systems.
+        Wiener function of the trajectory argument for time dependent systems.
 
         When used as an args:
 
-            ``QobjEvo([op, func], args={"W": SMESolver.WeinerFeedback()})``
+            ``QobjEvo([op, func], args={"W": SMESolver.WienerFeedback()})``
 
         The ``func`` will receive a function as ``W`` that return an array of
         wiener processes values at ``t``. The wiener process for the i-th
@@ -726,7 +729,7 @@ def WeinerFeedback(cls, default=None):
 
         .. note::
 
-            WeinerFeedback can't be added to a running solver when updating
+            WienerFeedback can't be added to a running solver when updating
             arguments between steps: ``solver.step(..., args={})``.
 
         Parameters
@@ -736,7 +739,7 @@ def WeinerFeedback(cls, default=None):
             When not passed, a function returning ``np.array([0])`` is used.
 
         """
-        return _WeinerFeedback(default)
+        return _WienerFeedback(default)
 
     @classmethod
     def StateFeedback(cls, default=None, raw_data=False):

qutip/tests/core/data/test_dense.py

@@ -305,3 +305,12 @@ def test_like_keep_order(func, fortran):
     old = dense.zeros(3, 3, fortran=fortran)
     new = func(old)
     assert new.fortran == old.fortran
+
+
+@pytest.mark.parametrize("shape", [(4, 3), (3, 3), (2, 8)])
+def test_inplace_matmul_error(shape):
+    op = dense.zeros(4, 4)
+    out = dense.zeros(*shape)
+    with pytest.raises(ValueError) as err:
+        data.matmul_dense(op, op, 1., out=out)
+    assert str(err.value).startswith("incompatible output shape")

qutip/tests/core/test_qobj.py

@@ -1264,7 +1264,13 @@ def test_data_as():
     assert "dia_matrix" in str(err.value)
 
 
-@pytest.mark.parametrize('dtype', ["CSR", "Dense"])
+@pytest.mark.parametrize('dtype', ["CSR", "Dense", "Dia"])
 def test_qobj_dtype(dtype):
     obj = qutip.qeye(2, dtype=dtype)
-    assert obj.dtype == qutip.data.to.parse(dtype)
+    assert obj.dtype == qutip.data.to.parse(dtype)
+
+
+@pytest.mark.parametrize('dtype', ["CSR", "Dense", "Dia"])
+def test_dtype_in_info_string(dtype):
+    obj = qutip.qeye(2, dtype=dtype)
+    assert dtype.lower() in str(obj).lower()