r"""Routines that can be passed as `propagator` to :func:`.optimize_pulses`
The numerical effort involved in the optimization is almost entirely within the
simulation of the system dynamics. In every iteration and for every objective,
the system must be "propagated" once forwards in time and once backwards in
time, see also :mod:`krotov.parallelization`.
The implementation of this time propagation must be inside the user-supplied
routine `propagator` that is passed to :func:`.optimize_pulses` and must
calculate the propagation over a single time step. In particular,
:func:`qutip.mesolve.mesolve` is not automatically used for simulating any
dynamics within the optimization. The signature for any `propagator`
must be the same as the "reference" :func:`expm` propagator:
>>> str(inspect.signature(krotov.propagators.expm))
'(H, state, dt, c_ops=None, backwards=False, initialize=False)'
The arguments are as follows (cf. :class:`Propagator`):
* `H` is the system Hamiltonian or Liouvillian, in a nested-list format similar
to that used by :func:`qutip.mesolve.mesolve`, e.g., for a Hamiltonian
$\Op{H} = \Op{H}_0 + c \Op{H}_1$, where $c$ is the value of a control field
at a particular point in time, `propagator` would receive a list ``[H0, [H1,
c]]`` where ``H0`` and ``H1`` are :class:`qutip.Qobj` operators.
The nested-list for `H` used here, with scalar values for the controls, is
obtained internally from the format used by :func:`~qutip.mesolve.mesolve`,
with time-dependent controls over the entire time grid, via
:func:`krotov.conversions.plug_in_pulse_values`.
* `state` is the :class:`qutip.Qobj` state that should be propagated, either a
Hilbert space state, or a density matrix.
* `dt` is the time step (a float). It is always positive, even for
``backwards=True``.
* `c_ops` is None, or a list of collapse (Lindblad) operators, where each list
element is a :class:`qutip.Qobj` instance (or possibly a nested list, for
time-dependent Lindblad operators. Note that is generally preferred for `H`
to be a Liouvillian, for dissipative dynamics.
* `backwards` (:class:`bool`): If passed as `True`, the `propagator` should
propagate backwards in time. In Hilbert space, this means using -`dt` instead
of `dt`. In Liouville space, there is no difference between forward and
backward propagation. In the context of Krotov's method, the backward
propagation uses the conjugate Hamiltonian, respectively Liouvillian.
However, the `propagator` routine does not need to be aware of this fact: it
will receive the appropriate `H` and `c_ops`.
* `initialize` (:class:`bool`): A flag to indicate the beginning of a
propagation over a time grid. If False in subsequent calls, the `propagator`
may assume that the input `state` is the result of the previous call to
`propagator`.
.. warning::
The routines in this module are provided with no guarantee to be either
general or efficient. The :func:`expm` propagator is exact to machine
precision, but generally extremely slow. For "production use", it is
recommended to supply a problem-specific `propagator` that is highly
optimized for speed. You might consider the use of Cython_. This is key to
minimize the runtime of the optimization.
The `initialize` flag enables "stateful" propagators that cache data between
calls. This can significantly improve numerical efficiency.
:class:`DensityMatrixODEPropagator` is an example for such a propagator. In
general, any stateful `propagator` should be an instance of
:class:`Propagator`.
.. _Cython: https://cython.org
"""
from abc import ABC, abstractmethod
import numpy as np
import qutip
import scipy
import threadpoolctl
from qutip.cy.spconvert import dense2D_to_fastcsr_fmode
from qutip.cy.spmatfuncs import spmvpy_csr
from qutip.superoperator import mat2vec, vec2mat
__all__ = ['expm', 'Propagator', 'DensityMatrixODEPropagator']
[docs]def expm(H, state, dt, c_ops=None, backwards=False, initialize=False):
"""Propagate using matrix exponentiation.
This supports `H` being a Hamiltonian (for a Hilbert space `state`) or a
Liouvillian (for `state` being a density matrix) in nested-list format.
Collapse operators `c_ops` are not supported. The propagator is not
stateful, thus `initialize` is ignored. The matrix exponentiation is
evaluated in single-threaded mode, to prevent accidental nested
parallelization.
"""
if c_ops is None:
c_ops = []
if len(c_ops) > 0:
raise NotImplementedError("Liouville exponentiation not implemented")
assert isinstance(H, list) and len(H) > 0
eqm_factor = -1j # factor in front of H on rhs of the equation of motion
if isinstance(H[0], list):
if H[0][1].type == 'super':
eqm_factor = 1
if backwards:
eqm_factor = eqm_factor.conjugate()
A = (eqm_factor * H[0][1]) * H[0][0]
else:
if H[0].type == 'super':
eqm_factor = 1
if backwards:
eqm_factor = eqm_factor.conjugate()
A = eqm_factor * H[0]
for part in H[1:]:
if isinstance(part, list):
A += (eqm_factor * part[1]) * part[0]
else:
A += eqm_factor * part
ok_types = (state.type == 'oper' and A.type == 'super') or (
state.type in ['ket', 'bra'] and A.type == 'oper'
)
if ok_types:
with threadpoolctl.threadpool_limits(limits=1):
return ((A * dt).expm())(state)
else:
raise NotImplementedError(
"Cannot handle argument types A:%s, state:%s"
% (A.type, state.type)
)
[docs]class Propagator(ABC):
"""Abstract base class for stateful propagators"""
[docs] @abstractmethod
def __call__(
self, H, state, dt, c_ops=None, backwards=False, initialize=False
):
"""Evaluation of a single propagation step
Args:
H (list): A Hamiltonian or Liouvillian in qutip's nested-list
format, with a scalar value in the place of a time-dependency.
For example, ``[H0, [H1, u]]`` for a drift Hamiltonian ``H0``,
a control Hamiltonian ``H1``, and a scalar value ``u`` that is
a time-dependent control evaluated for a particular point in
time.
state (qutip.Qobj): The state to propagate
dt (float): The time step over which to propagate
c_ops (list or None): A list of Lindblad operators. Using explicit
Lindblad operators should be avoided: it is usually more
efficient to convert them into a Lindbladian, passed as `H`
backwards (bool): Whether the propagation is forward in time or
backward in time
initialize (bool): Whether the propagator should (re-)initialize
for a new propagation, when the propagator is used to advance
on a time grid, `initialize` should be passed as True for the
initial time step (0 to `dt` in a forward propagation, or T to
T-dt for a backward propagation), and False otherwise.
Note:
A propagator may assume the propagation to be "sequential"
when `initialize` is False. That is, the state to propagate is the
result of the previous call to the propagator.
"""
pass
[docs]class DensityMatrixODEPropagator(Propagator):
"""Propagator for density matrix evolution under a Lindbladian
See :class:`qutip.solver.Options` for all arguments except `reentrant`.
Passing True for the `reentrant` re-initializes the propagator in every
time step.
Warning:
By default, the propagator is not "re-entrant". That is, you cannot use
more than one instance of :class:`DensityMatrixODEPropagator` in the
same process at the same time. This limitation is due to
:class:`scipy.integrate.ode` with the "zvode" integrator not being
re-entrant. Passing ``reentrant=True`` side-steps this problem by
re-initializating :class:`scipy.integrate.ode` in every time step. This
makes it possible to use :class:`DensityMatrixODEPropagator` in the
optimization of multiple objectives, but creates a significant
overhead.
"""
def __init__(
self,
method='adams',
order=12,
atol=1e-8,
rtol=1e-6,
nsteps=1000,
first_step=0,
min_step=0,
max_step=0,
reentrant=False,
):
self.method = method
self.order = order
self.atol = atol
self.rtol = rtol
self.nsteps = nsteps
self.first_step = first_step
self.min_step = min_step
self.max_step = max_step
self._L_list = None # cached Liouvillian data
self._control_indices = None # which indices in `L` have a control val
self._r = None # the integrator
self._t = 0.0 # time up to which we've integrated
self._y = None # current vectorized state
self.reentrant = reentrant
[docs] def __call__(
self, H, state, dt, c_ops=None, backwards=False, initialize=False
):
"""Evaluation of a single propagation step
Args:
H (list): A Liouvillian superoperator in qutip's nested-list
format, with a scalar value in the place of a time-dependency.
For example, ``[L0, [L1, u]]`` for a drift Liouvillian ``L0``,
a control Liouvillian ``H1``, and a scalar value ``u`` that is
a time-dependent control evaluated for a particular point in
time. If `initialize` is False, only the control values are
taken into account; any operators are assumed to be identical
to the internally cached values of `H` during initialization.
state (qutip.Qobj): The density matrix to propagate. The passed
value is ignored unless `initialize` is given as True.
Otherwise, it is assumed that `state` matches the (internally
stored) state that was the result from the previous propagation
step.
dt (float): The time step over which to propagate
c_ops (list or None): An empty list, or None. Since this propagator
assumes a full Liouvillian, it cannot be combined with Lindblad
operators.
backwards (bool): Whether the propagation is forward in time or
backward in time. Since the equation of motion for a
Liouvillian and conjugate Liouvillian is the same, this
parameter has no effect. Instead, for the backward propagation,
the conjugate Liouvillian must be passed for `L`.
initialize (bool): Whether to (re-)initialize for a new
propagation. This caches `H` (except for the control values)
and `state` internally.
"""
# H is really an L, but it's a very bad idea for a subclass not to
# follow the interface of the parent (including argument names).
# Internally, however, we'll use L instead of H
if initialize or self.reentrant:
self._initialize(H, state, dt, c_ops, backwards)
else:
if self.reentrant:
self._initialize_integrator(self._y)
# only update the control values
for i in self._control_indices:
self._L_list[i][1] = H[i][1]
self._t += dt
self._r.integrate(self._t)
self._y = self._r.y
return qutip.Qobj(
dense2D_to_fastcsr_fmode(
vec2mat(self._y), state.shape[0], state.shape[1]
),
dims=state.dims,
isherm=True,
)
@staticmethod
def _rhs(t, rho, L_list):
# _rhs being a staticmethod enables the propagator to be pickled (for
# parallelization)
out = np.zeros(rho.shape[0], dtype=complex)
L = L_list[0][0]
L_coeff = L_list[0][1]
spmvpy_csr(L.data, L.indices, L.indptr, rho, L_coeff, out)
for n in range(1, len(L_list)):
L = L_list[n][0]
L_coeff = L_list[n][1]
spmvpy_csr(L.data, L.indices, L.indptr, rho, L_coeff, out)
return out
def _initialize(self, L, rho, dt, c_ops, backwards):
self._initialize_data(L, rho, dt, c_ops, backwards)
self._initialize_integrator(self._y)
def _initialize_data(self, L, rho, dt, c_ops, backwards):
L_list = []
control_indices = []
if not (c_ops is None or len(c_ops) == 0):
# in principle, we could convert c_ops to a Lindbladian, here
raise NotImplementedError("c_ops not implemented")
for (i, spec) in enumerate(L):
if isinstance(spec, qutip.Qobj):
l_op = spec
l_coeff = 1
elif isinstance(spec, list) and isinstance(spec[0], qutip.Qobj):
l_op = spec[0]
l_coeff = spec[1]
control_indices.append(i)
else:
raise ValueError(
"Incorrect specification of time-dependent Liouvillian"
)
if l_op.type == 'super':
L_list.append([l_op.data, l_coeff, False])
else:
raise ValueError(
"Incorrect specification of time-dependent Liouvillian"
)
self._L_list = L_list
self._control_indices = control_indices
if rho.type == 'oper':
self._y = mat2vec(rho.full()).ravel('F') # initial state
else:
raise ValueError("rho must be a density matrix")
def _initialize_integrator(self, initial_vector):
r = scipy.integrate.ode(self._rhs)
r.set_integrator(
'zvode',
method=self.method,
order=self.order,
atol=self.atol,
rtol=self.rtol,
nsteps=self.nsteps,
first_step=self.first_step,
min_step=self.min_step,
max_step=self.max_step,
)
r.set_initial_value(initial_vector)
r.set_f_params(self._L_list)
self._r = r
self._t = 0.0