Merge branch 'master' of https://github.com/qutip/qutip into 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/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/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,13 +448,15 @@ 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.
         """
+        # 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)
@@ -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

qutip/solver/multitraj.py

@@ -8,23 +8,22 @@
 __all__ = ["MultiTrajSolver"]
 
 
-class _MTSystem:
+class _MultiTrajRHS:
     """
     Container for the operators of the solver.
     """
     def __init__(self, rhs):
         self.rhs = rhs
 
-    def __call__(self):
-        return self.rhs
-
     def arguments(self, args):
         self.rhs.arguments(args)
 
     def _register_feedback(self, type, val):
         pass
 
     def __getattr__(self, attr):
+        if attr == "rhs":
+            raise AttributeError
         if hasattr(self.rhs, attr):
             return getattr(self.rhs, attr)
         raise AttributeError
@@ -71,12 +70,11 @@ class MultiTrajSolver(Solver):
 
     def __init__(self, rhs, *, options=None):
         if isinstance(rhs, QobjEvo):
-            self.system = _MTSystem(rhs)
-        elif isinstance(rhs, _MTSystem):
-            self.system = rhs
+            self.rhs = _MultiTrajRHS(rhs)
+        elif isinstance(rhs, _MultiTrajRHS):
+            self.rhs = rhs
         else:
             raise TypeError("The system should be a QobjEvo")
-        self.rhs = self.system()
         self.options = options
         self.seed_sequence = np.random.SeedSequence()
         self._integrator = self._get_integrator()
@@ -233,23 +231,25 @@ def run(self, state, tlist, ntraj=1, *,
         result.stats['run time'] = time() - start_time
         return result
 
-    def _initialize_run_one_traj(self, seed, state, tlist, e_ops):
+    def _initialize_run_one_traj(self, seed, state, tlist, e_ops,
+                                 **integrator_kwargs):
         result = self._trajectory_resultclass(e_ops, self.options)
         generator = self._get_generator(seed)
-        self._integrator.set_state(tlist[0], state, generator)
+        self._integrator.set_state(tlist[0], state, generator,
+                                   **integrator_kwargs)
         result.add(tlist[0], self._restore_state(state, copy=False))
         return result
 
-    def _run_one_traj(self, seed, state, tlist, e_ops):
+    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)
+        result = self._initialize_run_one_traj(seed, state, tlist, e_ops,
+                                               **integrator_kwargs)
         return self._integrate_one_traj(seed, tlist, result)
 
     def _integrate_one_traj(self, seed, tlist, result):
-        for t in tlist[1:]:
-            t, state = self._integrator.integrate(t, copy=False)
+        for t, state in self._integrator.run(tlist):
             result.add(t, self._restore_state(state, copy=False))
         return seed, result
 
@@ -278,7 +278,8 @@ def _read_seed(self, seed, ntraj):
     def _argument(self, args):
         """Update the args, for the `rhs` and `c_ops` and other operators."""
         if args:
-            self.system.arguments(args)
+            self.rhs.arguments(args)
+            self._integrator.arguments(args)
 
     def _get_generator(self, seed):
         """

qutip/solver/stochastic.py

@@ -2,7 +2,7 @@
 
 from .sode.ssystem import StochasticOpenSystem, StochasticClosedSystem
 from .result import MultiTrajResult, Result, ExpectOp
-from .multitraj import MultiTrajSolver
+from .multitraj import _MultiTrajRHS, MultiTrajSolver
 from .. import Qobj, QobjEvo
 from ..core.dimensions import Dimensions
 import numpy as np
@@ -26,7 +26,7 @@ def _post_init(self, m_ops=(), dw_factor=(), heterodyne=False):
             self.m_ops.append(ExpectOp(op, f, self.m_expect[-1].append))
             self.add_processor(self.m_ops[-1]._store)
 
-    def add(self, t, state, noise):
+    def add(self, t, state, noise=None):
         super().add(t, state)
         if noise is not None and self.options["store_measurement"]:
             for i, dW in enumerate(noise):
@@ -166,7 +166,7 @@ def wiener_process(self):
         return self._trajectories_attr("wiener_process")
 
 
-class _StochasticRHS:
+class _StochasticRHS(_MultiTrajRHS):
     """
     In between object to store the stochastic system.
 
@@ -181,7 +181,7 @@ class _StochasticRHS:
     def __init__(self, issuper, H, sc_ops, c_ops, heterodyne):
 
         if not isinstance(H, (Qobj, QobjEvo)) or not H.isoper:
-            raise TypeError("The Hamiltonian must be am operator")
+            raise TypeError("The Hamiltonian must be an operator")
         self.H = QobjEvo(H)
 
         if isinstance(sc_ops, (Qobj, QobjEvo)):
@@ -500,8 +500,8 @@ class StochasticSolver(MultiTrajSolver):
     name = "StochasticSolver"
     _resultclass = StochasticResult
     _avail_integrators = {}
-    system = None
     _open = None
+
     solver_options = {
         "progress_bar": "text",
         "progress_kwargs": {"chunk_size": 10},
@@ -517,20 +517,22 @@ class StochasticSolver(MultiTrajSolver):
         "store_measurement": False,
     }
 
+    def _trajectory_resultclass(self, e_ops, options):
+        return StochasticTrajResult(
+            e_ops,
+            options,
+            m_ops=self.m_ops,
+            dw_factor=self.dW_factors,
+            heterodyne=self.heterodyne,
+        )
+
     def __init__(self, H, sc_ops, heterodyne, *, c_ops=(), options=None):
-        self.options = options
         self._heterodyne = heterodyne
         if self.name == "ssesolve" and c_ops:
             raise ValueError("c_ops are not supported by ssesolve.")
 
         rhs = _StochasticRHS(self._open, H, sc_ops, c_ops, heterodyne)
-        self.rhs = rhs
-        self.system = rhs
-        self.options = options
-        self.seed_sequence = np.random.SeedSequence()
-        self._integrator = self._get_integrator()
-        self._state_metadata = {}
-        self.stats = self._initialize_stats()
+        super().__init__(rhs, options=options)
 
         if heterodyne:
             self._m_ops = []
@@ -619,25 +621,9 @@ def dW_factors(self, new_dW_factors):
                 )
         self._dW_factors = new_dW_factors
 
-    def _run_one_traj(self, seed, state, tlist, e_ops):
-        """
-        Run one trajectory and return the result.
-        """
-        result = StochasticTrajResult(
-            e_ops,
-            self.options,
-            m_ops=self.m_ops,
-            dw_factor=self.dW_factors,
-            heterodyne=self.heterodyne,
-        )
-        generator = self._get_generator(seed)
-        self._integrator.set_state(tlist[0], state, generator)
-        state_t = self._restore_state(state, copy=False)
-        result.add(tlist[0], state_t, None)
-        for t in tlist[1:]:
-            t, state, noise = self._integrator.integrate(t, copy=False)
-            state_t = self._restore_state(state, copy=False)
-            result.add(t, state_t, noise)
+    def _integrate_one_traj(self, seed, tlist, result):
+        for t, state, noise in self._integrator.run(tlist):
+            result.add(t, self._restore_state(state, copy=False), noise)
         return seed, result
 
     @classmethod