We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
qutip.propagator and array with a single value
Hi,

I may have missed something but from the docstrings, I would expect the last two lines of this example script to actually return the same thing:

```python
>>> import numpy as np
>>> import qutip
>>> omega_0 = 1.0 * 2 * np.pi
>>> Omega = 0.1 * 2 * np.pi
>>> omega_d = 1.5 * 2 * np.pi
>>> T = 2 * np.pi / omega_d

>>> H0 = omega_0 / 2.0 * qutip.sigmaz()
>>> H1 = Omega * qutip.sigmax()
>>> H2 = Omega * qutip.sigmay()

>>> H = [H0, [H1, 'cos(w_d * t)'], [H2, 'sin(w_d * t)']]
>>> args = {'w_d': omega_d}

>>> (qutip.propagator(H, t, c_op_list=[], args=args) - qutip.propagator(H, [t], c_op_list=[], args=args)).norm()
0.04271636357907908

>>> (qutip.propagator(H, t, c_op_list=[], args=args) - qutip.propagator(H, [0, t], c_op_list=[], args=args)).norm()
0.0
```

That is, I would expect `qutip.propagator(H, t, c_op_list=[], args=args)` to give the same result as `qutip.propagator(H, [t], c_op_list=[], args=args)` (single `t` value or array of times with a single `t` element). This does not seem to be the case here, with a not negligible difference.

Is this an expected behavior (which would probably deserve more emphasis in the doc)?

Thanks!
Best,

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `qutip/solver/propagator.py`
Content:
```
1 __all__ = ['Propagator', 'propagator', 'propagator_steadystate']
2 
3 import numbers
4 import numpy as np
5 
6 from .. import Qobj, qeye, qeye_like, unstack_columns, QobjEvo, liouvillian
7 from ..core import data as _data
8 from .mesolve import mesolve, MESolver
9 from .sesolve import sesolve, SESolver
10 from .heom.bofin_solvers import HEOMSolver
11 from .solver_base import Solver
12 from .multitraj import MultiTrajSolver
13 
14 
15 def propagator(H, t, c_ops=(), args=None, options=None, **kwargs):
16     r"""
17     Calculate the propagator U(t) for the density matrix or wave function such
18     that :math:`\psi(t) = U(t)\psi(0)` or
19     :math:`\rho_{\mathrm vec}(t) = U(t) \rho_{\mathrm vec}(0)`
20     where :math:`\rho_{\mathrm vec}` is the vector representation of the
21     density matrix.
22 
23     Parameters
24     ----------
25     H : :obj:`.Qobj`, :obj:`.QobjEvo`, :obj:`.QobjEvo` compatible format
26         Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or
27         QobjEvo. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`] or callable
28         that can be made into :obj:`.QobjEvo` are also accepted.
29 
30     t : float or array-like
31         Time or list of times for which to evaluate the propagator.
32 
33     c_ops : list, optional
34         List of Qobj or QobjEvo collapse operators.
35 
36     args : dictionary, optional
37         Parameters to callback functions for time-dependent Hamiltonians and
38         collapse operators.
39 
40     options : dict, optional
41         Options for the solver.
42 
43     **kwargs :
44         Extra parameters to use when creating the
45         :obj:`.QobjEvo` from a list format ``H``.
46 
47     Returns
48     -------
49     U : :obj:`.Qobj`, list
50         Instance representing the propagator(s) :math:`U(t)`. Return a single
51         Qobj when ``t`` is a number or a list when ``t`` is a list.
52 
53     """
54     if isinstance(t, numbers.Real):
55         tlist = [0, t]
56         list_output = False
57     else:
58         tlist = t
59         list_output = True
60 
61     if not isinstance(H, (Qobj, QobjEvo)):
62         H = QobjEvo(H, args=args, **kwargs)
63 
64     if c_ops:
65         H = liouvillian(H, c_ops)
66 
67     U0 = qeye_like(H)
68 
69     if H.issuper:
70         out = mesolve(H, U0, tlist, args=args, options=options).states
71     else:
72         out = sesolve(H, U0, tlist, args=args, options=options).states
73 
74     if list_output:
75         return out
76     else:
77         return out[-1]
78 
79 
80 def propagator_steadystate(U):
81     r"""Find the steady state for successive applications of the propagator
82     :math:`U`.
83 
84     Parameters
85     ----------
86     U : :obj:`.Qobj`
87         Operator representing the propagator.
88 
89     Returns
90     -------
91     a : :obj:`.Qobj`
92         Instance representing the steady-state density matrix.
93     """
94     evals, estates = U.eigenstates()
95     shifted_vals = np.abs(evals - 1.0)
96     ev_idx = np.argmin(shifted_vals)
97     rho_data = unstack_columns(estates[ev_idx].data)
98     rho_data = _data.mul(rho_data, 0.5 / _data.trace(rho_data))
99     return Qobj(_data.add(rho_data, _data.adjoint(rho_data)),
100                 dims=U.dims[0],
101                 isherm=True,
102                 copy=False)
103 
104 
105 class Propagator:
106     """
107     A generator of propagator for a system.
108 
109     Usage:
110 
111         U = Propagator(H, c_ops)
112 
113         psi_t = U(t) @ psi_0
114 
115     Save some previously computed propagator are stored to speed up subsequent
116     computation. Changing ``args`` will erase these stored probagator.
117 
118     Parameters
119     ----------
120     system : :obj:`.Qobj`, :obj:`.QobjEvo`, :class:`.Solver`
121         Possibly time-dependent system driving the evolution, either already
122         packaged in a solver, such as :class:`.SESolver` or :class:`.BRSolver`,
123         or the Liouvillian or Hamiltonian as a :obj:`.Qobj`,
124         :obj:`.QobjEvo`. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`]
125         or callable that can be made into :obj:`.QobjEvo` are also accepted.
126 
127         Solvers that run non-deterministacilly, such as :class:`.MCSolver`, are
128         not supported.
129 
130     c_ops : list, optional
131         List of :obj:`.Qobj` or :obj:`.QobjEvo` collapse operators.
132 
133     args : dictionary, optional
134         Parameters to callback functions for time-dependent Hamiltonians and
135         collapse operators.
136 
137     options : dict, optional
138         Options for the solver.
139 
140     memoize : int, default: 10
141         Max number of propagator to save.
142 
143     tol : float, default: 1e-14
144         Absolute tolerance for the time. If a previous propagator was computed
145         at a time within tolerance, that propagator will be returned.
146 
147     Notes
148     -----
149     The :class:`Propagator` is not a :obj:`.QobjEvo` so
150     it cannot be used for operations with :obj:`.Qobj` or
151     :obj:`.QobjEvo`. It can be made into a
152     :obj:`.QobjEvo` with ::
153 
154         U = QobjEvo(Propagator(H))
155 
156     """
157     def __init__(self, system, *, c_ops=(), args=None, options=None,
158                  memoize=10, tol=1e-14):
159         if isinstance(system, MultiTrajSolver):
160             raise TypeError("Non-deterministic solvers cannot be used "
161                             "as a propagator system")
162         elif isinstance(system, HEOMSolver):
163             raise NotImplementedError(
164                 "HEOM is not supported by Propagator. "
165                 "Please, tell us on GitHub issues if you need it!"
166             )
167         elif isinstance(system, Solver):
168             self.solver = system
169         else:
170             Hevo = QobjEvo(system, args=args)
171             c_ops = [QobjEvo(op, args=args) for op in c_ops]
172             if Hevo.issuper or c_ops:
173                 self.solver = MESolver(Hevo, c_ops=c_ops, options=options)
174             else:
175                 self.solver = SESolver(Hevo, options=options)
176 
177         self.times = [0]
178         self.invs = [None]
179         self.props = [qeye(self.solver.sys_dims)]
180         self.solver.start(self.props[0], self.times[0])
181         self.cte = self.solver.rhs.isconstant
182         H_0 = self.solver.rhs(0)
183         self.unitary = not H_0.issuper and H_0.isherm
184         self.args = args
185         self.memoize = max(3, int(memoize))
186         self.tol = tol
187 
188     def _lookup_or_compute(self, t):
189         """
190         Get U(t) from cache or compute it.
191         """
192         idx = np.searchsorted(self.times, t)
193         if idx < len(self.times) and abs(t-self.times[idx]) <= self.tol:
194             U = self.props[idx]
195         elif idx > 0 and abs(t-self.times[idx-1]) <= self.tol:
196             U = self.props[idx-1]
197         else:
198             U = self._compute(t, idx)
199             self._insert(t, U, idx)
200         return U
201 
202     def __call__(self, t, t_start=0, **args):
203         """
204         Get the propagator from ``t_start`` to ``t``.
205 
206         Parameters
207         ----------
208         t : float
209             Time at which to compute the propagator.
210         t_start: float [0]
211             Time at which the propagator start such that:
212                 ``psi[t] = U.prop(t, t_start) @ psi[t_start]``
213         args : dict
214             Argument to pass to a time dependent Hamiltonian.
215             Updating ``args`` take effect since ``t=0`` and the new ``args``
216             will be used in future call.
217         """
218         # We could improve it when the system is constant using U(2t) = U(t)**2
219         if not self.cte and args and args != self.args:
220             self.args = args
221             self.solver._argument(args)
222             self.times = [0]
223             self.props = [qeye_like(self.props[0])]
224             self.solver.start(self.props[0], self.times[0])
225 
226         if t_start:
227             if t == t_start:
228                 U = self._lookup_or_compute(0)
229             if self.cte:
230                 U = self._lookup_or_compute(t - t_start)
231             else:
232                 Uinv = self._inv(self._lookup_or_compute(t_start))
233                 U = self._lookup_or_compute(t) @ Uinv
234         else:
235             U = self._lookup_or_compute(t)
236         return U
237 
238     def inv(self, t, **args):
239         """
240         Get the inverse of the propagator at ``t``, such that
241             ``psi_0 = U.inv(t) @ psi_t``
242 
243         Parameters
244         ----------
245         t : float
246             Time at which to compute the propagator.
247         args : dict
248             Argument to pass to a time dependent Hamiltonian.
249             Updating ``args`` take effect since ``t=0`` and the new ``args``
250             will be used in future call.
251         """
252         return self._inv(self(t, **args))
253 
254     def _compute(self, t, idx):
255         """
256         Compute the propagator at ``t``, ``idx`` point to a pair of
257         (time, propagator) close to the desired time.
258         """
259         t_last = self.solver._integrator.get_state(copy=False)[0]
260         if self.times[idx-1] <= t_last <= t:
261             U = self.solver.step(t)
262         elif idx > 0:
263             self.solver.start(self.props[idx-1], self.times[idx-1])
264             U = self.solver.step(t)
265         else:
266             # Evolving backward in time is not supported by all integrator.
267             self.solver.start(qeye_like(self.props[0]), t)
268             Uinv = self.solver.step(self.times[idx])
269             U = self._inv(Uinv)
270         return U
271 
272     def _inv(self, U):
273         return U.dag() if self.unitary else U.inv()
274 
275     def _insert(self, t, U, idx):
276         """
277         Insert a new pair of (time, propagator) to the memorized states.
278         """
279         while len(self.times) >= self.memoize:
280             rm_idx = self.memoize // 2
281             if self.times[rm_idx] < t:
282                 idx -= 1
283             del self.times[rm_idx]
284             del self.props[rm_idx]
285         self.times.insert(idx, t)
286         self.props.insert(idx, U)
287
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/qutip/solver/propagator.py b/qutip/solver/propagator.py
--- a/qutip/solver/propagator.py
+++ b/qutip/solver/propagator.py
@@ -28,7 +28,11 @@
         that can be made into :obj:`.QobjEvo` are also accepted.
 
     t : float or array-like
-        Time or list of times for which to evaluate the propagator.
+        Time or list of times for which to evaluate the propagator. If a single
+        time ``t`` is passed, the propagator from ``0`` to ``t`` is computed.
+        When ``t`` is a list, the propagators from the first time in the list
+        to each elements in ``t`` is returned. In that case, the first output
+        will always be the identity matrix.
 
     c_ops : list, optional
         List of Qobj or QobjEvo collapse operators.

{"golden_diff": "diff --git a/qutip/solver/propagator.py b/qutip/solver/propagator.py\n--- a/qutip/solver/propagator.py\n+++ b/qutip/solver/propagator.py\n@@ -28,7 +28,11 @@\n         that can be made into :obj:`.QobjEvo` are also accepted.\n \n     t : float or array-like\n-        Time or list of times for which to evaluate the propagator.\n+        Time or list of times for which to evaluate the propagator. If a single\n+        time ``t`` is passed, the propagator from ``0`` to ``t`` is computed.\n+        When ``t`` is a list, the propagators from the first time in the list\n+        to each elements in ``t`` is returned. In that case, the first output\n+        will always be the identity matrix.\n \n     c_ops : list, optional\n         List of Qobj or QobjEvo collapse operators.\n", "issue": "qutip.propagator and array with a single value\nHi,\r\n\r\nI may have missed something but from the docstrings, I would expect the last two lines of this example script to actually return the same thing:\r\n\r\n```python\r\n>>> import numpy as np\r\n>>> import qutip\r\n>>> omega_0 = 1.0 * 2 * np.pi\r\n>>> Omega = 0.1 * 2 * np.pi\r\n>>> omega_d = 1.5 * 2 * np.pi\r\n>>> T = 2 * np.pi / omega_d\r\n\r\n>>> H0 = omega_0 / 2.0 * qutip.sigmaz()\r\n>>> H1 = Omega * qutip.sigmax()\r\n>>> H2 = Omega * qutip.sigmay()\r\n\r\n>>> H = [H0, [H1, 'cos(w_d * t)'], [H2, 'sin(w_d * t)']]\r\n>>> args = {'w_d': omega_d}\r\n\r\n>>> (qutip.propagator(H, t, c_op_list=[], args=args) - qutip.propagator(H, [t], c_op_list=[], args=args)).norm()\r\n0.04271636357907908\r\n\r\n>>> (qutip.propagator(H, t, c_op_list=[], args=args) - qutip.propagator(H, [0, t], c_op_list=[], args=args)).norm()\r\n0.0\r\n```\r\n\r\nThat is, I would expect `qutip.propagator(H, t, c_op_list=[], args=args)` to give the same result as `qutip.propagator(H, [t], c_op_list=[], args=args)` (single `t` value or array of times with a single `t` element). This does not seem to be the case here, with a not negligible difference.\r\n\r\nIs this an expected behavior (which would probably deserve more emphasis in the doc)?\r\n\r\nThanks!\r\nBest,\n", "before_files": [{"content": "__all__ = ['Propagator', 'propagator', 'propagator_steadystate']\n\nimport numbers\nimport numpy as np\n\nfrom .. import Qobj, qeye, qeye_like, unstack_columns, QobjEvo, liouvillian\nfrom ..core import data as _data\nfrom .mesolve import mesolve, MESolver\nfrom .sesolve import sesolve, SESolver\nfrom .heom.bofin_solvers import HEOMSolver\nfrom .solver_base import Solver\nfrom .multitraj import MultiTrajSolver\n\n\ndef propagator(H, t, c_ops=(), args=None, options=None, **kwargs):\n    r\"\"\"\n    Calculate the propagator U(t) for the density matrix or wave function such\n    that :math:`\\psi(t) = U(t)\\psi(0)` or\n    :math:`\\rho_{\\mathrm vec}(t) = U(t) \\rho_{\\mathrm vec}(0)`\n    where :math:`\\rho_{\\mathrm vec}` is the vector representation of the\n    density matrix.\n\n    Parameters\n    ----------\n    H : :obj:`.Qobj`, :obj:`.QobjEvo`, :obj:`.QobjEvo` compatible format\n        Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or\n        QobjEvo. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`] or callable\n        that can be made into :obj:`.QobjEvo` are also accepted.\n\n    t : float or array-like\n        Time or list of times for which to evaluate the propagator.\n\n    c_ops : list, optional\n        List of Qobj or QobjEvo collapse operators.\n\n    args : dictionary, optional\n        Parameters to callback functions for time-dependent Hamiltonians and\n        collapse operators.\n\n    options : dict, optional\n        Options for the solver.\n\n    **kwargs :\n        Extra parameters to use when creating the\n        :obj:`.QobjEvo` from a list format ``H``.\n\n    Returns\n    -------\n    U : :obj:`.Qobj`, list\n        Instance representing the propagator(s) :math:`U(t)`. Return a single\n        Qobj when ``t`` is a number or a list when ``t`` is a list.\n\n    \"\"\"\n    if isinstance(t, numbers.Real):\n        tlist = [0, t]\n        list_output = False\n    else:\n        tlist = t\n        list_output = True\n\n    if not isinstance(H, (Qobj, QobjEvo)):\n        H = QobjEvo(H, args=args, **kwargs)\n\n    if c_ops:\n        H = liouvillian(H, c_ops)\n\n    U0 = qeye_like(H)\n\n    if H.issuper:\n        out = mesolve(H, U0, tlist, args=args, options=options).states\n    else:\n        out = sesolve(H, U0, tlist, args=args, options=options).states\n\n    if list_output:\n        return out\n    else:\n        return out[-1]\n\n\ndef propagator_steadystate(U):\n    r\"\"\"Find the steady state for successive applications of the propagator\n    :math:`U`.\n\n    Parameters\n    ----------\n    U : :obj:`.Qobj`\n        Operator representing the propagator.\n\n    Returns\n    -------\n    a : :obj:`.Qobj`\n        Instance representing the steady-state density matrix.\n    \"\"\"\n    evals, estates = U.eigenstates()\n    shifted_vals = np.abs(evals - 1.0)\n    ev_idx = np.argmin(shifted_vals)\n    rho_data = unstack_columns(estates[ev_idx].data)\n    rho_data = _data.mul(rho_data, 0.5 / _data.trace(rho_data))\n    return Qobj(_data.add(rho_data, _data.adjoint(rho_data)),\n                dims=U.dims[0],\n                isherm=True,\n                copy=False)\n\n\nclass Propagator:\n    \"\"\"\n    A generator of propagator for a system.\n\n    Usage:\n\n        U = Propagator(H, c_ops)\n\n        psi_t = U(t) @ psi_0\n\n    Save some previously computed propagator are stored to speed up subsequent\n    computation. Changing ``args`` will erase these stored probagator.\n\n    Parameters\n    ----------\n    system : :obj:`.Qobj`, :obj:`.QobjEvo`, :class:`.Solver`\n        Possibly time-dependent system driving the evolution, either already\n        packaged in a solver, such as :class:`.SESolver` or :class:`.BRSolver`,\n        or the Liouvillian or Hamiltonian as a :obj:`.Qobj`,\n        :obj:`.QobjEvo`. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`]\n        or callable that can be made into :obj:`.QobjEvo` are also accepted.\n\n        Solvers that run non-deterministacilly, such as :class:`.MCSolver`, are\n        not supported.\n\n    c_ops : list, optional\n        List of :obj:`.Qobj` or :obj:`.QobjEvo` collapse operators.\n\n    args : dictionary, optional\n        Parameters to callback functions for time-dependent Hamiltonians and\n        collapse operators.\n\n    options : dict, optional\n        Options for the solver.\n\n    memoize : int, default: 10\n        Max number of propagator to save.\n\n    tol : float, default: 1e-14\n        Absolute tolerance for the time. If a previous propagator was computed\n        at a time within tolerance, that propagator will be returned.\n\n    Notes\n    -----\n    The :class:`Propagator` is not a :obj:`.QobjEvo` so\n    it cannot be used for operations with :obj:`.Qobj` or\n    :obj:`.QobjEvo`. It can be made into a\n    :obj:`.QobjEvo` with ::\n\n        U = QobjEvo(Propagator(H))\n\n    \"\"\"\n    def __init__(self, system, *, c_ops=(), args=None, options=None,\n                 memoize=10, tol=1e-14):\n        if isinstance(system, MultiTrajSolver):\n            raise TypeError(\"Non-deterministic solvers cannot be used \"\n                            \"as a propagator system\")\n        elif isinstance(system, HEOMSolver):\n            raise NotImplementedError(\n                \"HEOM is not supported by Propagator. \"\n                \"Please, tell us on GitHub issues if you need it!\"\n            )\n        elif isinstance(system, Solver):\n            self.solver = system\n        else:\n            Hevo = QobjEvo(system, args=args)\n            c_ops = [QobjEvo(op, args=args) for op in c_ops]\n            if Hevo.issuper or c_ops:\n                self.solver = MESolver(Hevo, c_ops=c_ops, options=options)\n            else:\n                self.solver = SESolver(Hevo, options=options)\n\n        self.times = [0]\n        self.invs = [None]\n        self.props = [qeye(self.solver.sys_dims)]\n        self.solver.start(self.props[0], self.times[0])\n        self.cte = self.solver.rhs.isconstant\n        H_0 = self.solver.rhs(0)\n        self.unitary = not H_0.issuper and H_0.isherm\n        self.args = args\n        self.memoize = max(3, int(memoize))\n        self.tol = tol\n\n    def _lookup_or_compute(self, t):\n        \"\"\"\n        Get U(t) from cache or compute it.\n        \"\"\"\n        idx = np.searchsorted(self.times, t)\n        if idx < len(self.times) and abs(t-self.times[idx]) <= self.tol:\n            U = self.props[idx]\n        elif idx > 0 and abs(t-self.times[idx-1]) <= self.tol:\n            U = self.props[idx-1]\n        else:\n            U = self._compute(t, idx)\n            self._insert(t, U, idx)\n        return U\n\n    def __call__(self, t, t_start=0, **args):\n        \"\"\"\n        Get the propagator from ``t_start`` to ``t``.\n\n        Parameters\n        ----------\n        t : float\n            Time at which to compute the propagator.\n        t_start: float [0]\n            Time at which the propagator start such that:\n                ``psi[t] = U.prop(t, t_start) @ psi[t_start]``\n        args : dict\n            Argument to pass to a time dependent Hamiltonian.\n            Updating ``args`` take effect since ``t=0`` and the new ``args``\n            will be used in future call.\n        \"\"\"\n        # We could improve it when the system is constant using U(2t) = U(t)**2\n        if not self.cte and args and args != self.args:\n            self.args = args\n            self.solver._argument(args)\n            self.times = [0]\n            self.props = [qeye_like(self.props[0])]\n            self.solver.start(self.props[0], self.times[0])\n\n        if t_start:\n            if t == t_start:\n                U = self._lookup_or_compute(0)\n            if self.cte:\n                U = self._lookup_or_compute(t - t_start)\n            else:\n                Uinv = self._inv(self._lookup_or_compute(t_start))\n                U = self._lookup_or_compute(t) @ Uinv\n        else:\n            U = self._lookup_or_compute(t)\n        return U\n\n    def inv(self, t, **args):\n        \"\"\"\n        Get the inverse of the propagator at ``t``, such that\n            ``psi_0 = U.inv(t) @ psi_t``\n\n        Parameters\n        ----------\n        t : float\n            Time at which to compute the propagator.\n        args : dict\n            Argument to pass to a time dependent Hamiltonian.\n            Updating ``args`` take effect since ``t=0`` and the new ``args``\n            will be used in future call.\n        \"\"\"\n        return self._inv(self(t, **args))\n\n    def _compute(self, t, idx):\n        \"\"\"\n        Compute the propagator at ``t``, ``idx`` point to a pair of\n        (time, propagator) close to the desired time.\n        \"\"\"\n        t_last = self.solver._integrator.get_state(copy=False)[0]\n        if self.times[idx-1] <= t_last <= t:\n            U = self.solver.step(t)\n        elif idx > 0:\n            self.solver.start(self.props[idx-1], self.times[idx-1])\n            U = self.solver.step(t)\n        else:\n            # Evolving backward in time is not supported by all integrator.\n            self.solver.start(qeye_like(self.props[0]), t)\n            Uinv = self.solver.step(self.times[idx])\n            U = self._inv(Uinv)\n        return U\n\n    def _inv(self, U):\n        return U.dag() if self.unitary else U.inv()\n\n    def _insert(self, t, U, idx):\n        \"\"\"\n        Insert a new pair of (time, propagator) to the memorized states.\n        \"\"\"\n        while len(self.times) >= self.memoize:\n            rm_idx = self.memoize // 2\n            if self.times[rm_idx] < t:\n                idx -= 1\n            del self.times[rm_idx]\n            del self.props[rm_idx]\n        self.times.insert(idx, t)\n        self.props.insert(idx, U)\n", "path": "qutip/solver/propagator.py"}], "after_files": [{"content": "__all__ = ['Propagator', 'propagator', 'propagator_steadystate']\n\nimport numbers\nimport numpy as np\n\nfrom .. import Qobj, qeye, qeye_like, unstack_columns, QobjEvo, liouvillian\nfrom ..core import data as _data\nfrom .mesolve import mesolve, MESolver\nfrom .sesolve import sesolve, SESolver\nfrom .heom.bofin_solvers import HEOMSolver\nfrom .solver_base import Solver\nfrom .multitraj import MultiTrajSolver\n\n\ndef propagator(H, t, c_ops=(), args=None, options=None, **kwargs):\n    r\"\"\"\n    Calculate the propagator U(t) for the density matrix or wave function such\n    that :math:`\\psi(t) = U(t)\\psi(0)` or\n    :math:`\\rho_{\\mathrm vec}(t) = U(t) \\rho_{\\mathrm vec}(0)`\n    where :math:`\\rho_{\\mathrm vec}` is the vector representation of the\n    density matrix.\n\n    Parameters\n    ----------\n    H : :obj:`.Qobj`, :obj:`.QobjEvo`, :obj:`.QobjEvo` compatible format\n        Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or\n        QobjEvo. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`] or callable\n        that can be made into :obj:`.QobjEvo` are also accepted.\n\n    t : float or array-like\n        Time or list of times for which to evaluate the propagator. If a single\n        time ``t`` is passed, the propagator from ``0`` to ``t`` is computed.\n        When ``t`` is a list, the propagators from the first time in the list\n        to each elements in ``t`` is returned. In that case, the first output\n        will always be the identity matrix.\n\n    c_ops : list, optional\n        List of Qobj or QobjEvo collapse operators.\n\n    args : dictionary, optional\n        Parameters to callback functions for time-dependent Hamiltonians and\n        collapse operators.\n\n    options : dict, optional\n        Options for the solver.\n\n    **kwargs :\n        Extra parameters to use when creating the\n        :obj:`.QobjEvo` from a list format ``H``.\n\n    Returns\n    -------\n    U : :obj:`.Qobj`, list\n        Instance representing the propagator(s) :math:`U(t)`. Return a single\n        Qobj when ``t`` is a number or a list when ``t`` is a list.\n\n    \"\"\"\n    if isinstance(t, numbers.Real):\n        tlist = [0, t]\n        list_output = False\n    else:\n        tlist = t\n        list_output = True\n\n    if not isinstance(H, (Qobj, QobjEvo)):\n        H = QobjEvo(H, args=args, **kwargs)\n\n    if c_ops:\n        H = liouvillian(H, c_ops)\n\n    U0 = qeye_like(H)\n\n    if H.issuper:\n        out = mesolve(H, U0, tlist, args=args, options=options).states\n    else:\n        out = sesolve(H, U0, tlist, args=args, options=options).states\n\n    if list_output:\n        return out\n    else:\n        return out[-1]\n\n\ndef propagator_steadystate(U):\n    r\"\"\"Find the steady state for successive applications of the propagator\n    :math:`U`.\n\n    Parameters\n    ----------\n    U : :obj:`.Qobj`\n        Operator representing the propagator.\n\n    Returns\n    -------\n    a : :obj:`.Qobj`\n        Instance representing the steady-state density matrix.\n    \"\"\"\n    evals, estates = U.eigenstates()\n    shifted_vals = np.abs(evals - 1.0)\n    ev_idx = np.argmin(shifted_vals)\n    rho_data = unstack_columns(estates[ev_idx].data)\n    rho_data = _data.mul(rho_data, 0.5 / _data.trace(rho_data))\n    return Qobj(_data.add(rho_data, _data.adjoint(rho_data)),\n                dims=U.dims[0],\n                isherm=True,\n                copy=False)\n\n\nclass Propagator:\n    \"\"\"\n    A generator of propagator for a system.\n\n    Usage:\n\n        U = Propagator(H, c_ops)\n\n        psi_t = U(t) @ psi_0\n\n    Save some previously computed propagator are stored to speed up subsequent\n    computation. Changing ``args`` will erase these stored probagator.\n\n    Parameters\n    ----------\n    system : :obj:`.Qobj`, :obj:`.QobjEvo`, :class:`.Solver`\n        Possibly time-dependent system driving the evolution, either already\n        packaged in a solver, such as :class:`.SESolver` or :class:`.BRSolver`,\n        or the Liouvillian or Hamiltonian as a :obj:`.Qobj`,\n        :obj:`.QobjEvo`. ``list`` of [:obj:`.Qobj`, :obj:`.Coefficient`]\n        or callable that can be made into :obj:`.QobjEvo` are also accepted.\n\n        Solvers that run non-deterministacilly, such as :class:`.MCSolver`, are\n        not supported.\n\n    c_ops : list, optional\n        List of :obj:`.Qobj` or :obj:`.QobjEvo` collapse operators.\n\n    args : dictionary, optional\n        Parameters to callback functions for time-dependent Hamiltonians and\n        collapse operators.\n\n    options : dict, optional\n        Options for the solver.\n\n    memoize : int, default: 10\n        Max number of propagator to save.\n\n    tol : float, default: 1e-14\n        Absolute tolerance for the time. If a previous propagator was computed\n        at a time within tolerance, that propagator will be returned.\n\n    Notes\n    -----\n    The :class:`Propagator` is not a :obj:`.QobjEvo` so\n    it cannot be used for operations with :obj:`.Qobj` or\n    :obj:`.QobjEvo`. It can be made into a\n    :obj:`.QobjEvo` with ::\n\n        U = QobjEvo(Propagator(H))\n\n    \"\"\"\n    def __init__(self, system, *, c_ops=(), args=None, options=None,\n                 memoize=10, tol=1e-14):\n        if isinstance(system, MultiTrajSolver):\n            raise TypeError(\"Non-deterministic solvers cannot be used \"\n                            \"as a propagator system\")\n        elif isinstance(system, HEOMSolver):\n            raise NotImplementedError(\n                \"HEOM is not supported by Propagator. \"\n                \"Please, tell us on GitHub issues if you need it!\"\n            )\n        elif isinstance(system, Solver):\n            self.solver = system\n        else:\n            Hevo = QobjEvo(system, args=args)\n            c_ops = [QobjEvo(op, args=args) for op in c_ops]\n            if Hevo.issuper or c_ops:\n                self.solver = MESolver(Hevo, c_ops=c_ops, options=options)\n            else:\n                self.solver = SESolver(Hevo, options=options)\n\n        self.times = [0]\n        self.invs = [None]\n        self.props = [qeye(self.solver.sys_dims)]\n        self.solver.start(self.props[0], self.times[0])\n        self.cte = self.solver.rhs.isconstant\n        H_0 = self.solver.rhs(0)\n        self.unitary = not H_0.issuper and H_0.isherm\n        self.args = args\n        self.memoize = max(3, int(memoize))\n        self.tol = tol\n\n    def _lookup_or_compute(self, t):\n        \"\"\"\n        Get U(t) from cache or compute it.\n        \"\"\"\n        idx = np.searchsorted(self.times, t)\n        if idx < len(self.times) and abs(t-self.times[idx]) <= self.tol:\n            U = self.props[idx]\n        elif idx > 0 and abs(t-self.times[idx-1]) <= self.tol:\n            U = self.props[idx-1]\n        else:\n            U = self._compute(t, idx)\n            self._insert(t, U, idx)\n        return U\n\n    def __call__(self, t, t_start=0, **args):\n        \"\"\"\n        Get the propagator from ``t_start`` to ``t``.\n\n        Parameters\n        ----------\n        t : float\n            Time at which to compute the propagator.\n        t_start: float [0]\n            Time at which the propagator start such that:\n                ``psi[t] = U.prop(t, t_start) @ psi[t_start]``\n        args : dict\n            Argument to pass to a time dependent Hamiltonian.\n            Updating ``args`` take effect since ``t=0`` and the new ``args``\n            will be used in future call.\n        \"\"\"\n        # We could improve it when the system is constant using U(2t) = U(t)**2\n        if not self.cte and args and args != self.args:\n            self.args = args\n            self.solver._argument(args)\n            self.times = [0]\n            self.props = [qeye_like(self.props[0])]\n            self.solver.start(self.props[0], self.times[0])\n\n        if t_start:\n            if t == t_start:\n                U = self._lookup_or_compute(0)\n            if self.cte:\n                U = self._lookup_or_compute(t - t_start)\n            else:\n                Uinv = self._inv(self._lookup_or_compute(t_start))\n                U = self._lookup_or_compute(t) @ Uinv\n        else:\n            U = self._lookup_or_compute(t)\n        return U\n\n    def inv(self, t, **args):\n        \"\"\"\n        Get the inverse of the propagator at ``t``, such that\n            ``psi_0 = U.inv(t) @ psi_t``\n\n        Parameters\n        ----------\n        t : float\n            Time at which to compute the propagator.\n        args : dict\n            Argument to pass to a time dependent Hamiltonian.\n            Updating ``args`` take effect since ``t=0`` and the new ``args``\n            will be used in future call.\n        \"\"\"\n        return self._inv(self(t, **args))\n\n    def _compute(self, t, idx):\n        \"\"\"\n        Compute the propagator at ``t``, ``idx`` point to a pair of\n        (time, propagator) close to the desired time.\n        \"\"\"\n        t_last = self.solver._integrator.get_state(copy=False)[0]\n        if self.times[idx-1] <= t_last <= t:\n            U = self.solver.step(t)\n        elif idx > 0:\n            self.solver.start(self.props[idx-1], self.times[idx-1])\n            U = self.solver.step(t)\n        else:\n            # Evolving backward in time is not supported by all integrator.\n            self.solver.start(qeye_like(self.props[0]), t)\n            Uinv = self.solver.step(self.times[idx])\n            U = self._inv(Uinv)\n        return U\n\n    def _inv(self, U):\n        return U.dag() if self.unitary else U.inv()\n\n    def _insert(self, t, U, idx):\n        \"\"\"\n        Insert a new pair of (time, propagator) to the memorized states.\n        \"\"\"\n        while len(self.times) >= self.memoize:\n            rm_idx = self.memoize // 2\n            if self.times[rm_idx] < t:\n                idx -= 1\n            del self.times[rm_idx]\n            del self.props[rm_idx]\n        self.times.insert(idx, t)\n        self.props.insert(idx, U)\n", "path": "qutip/solver/propagator.py"}]}