Propagator interface

[goerz]

At the moment, QuantumPropagators is entirely focused on piecewise-constant control fields. In order to serve as the simulation backend for optimal control methods beyond Krotov and GRAPE (e.g., CRAB, GROUP, GOAT), this will have to be expanded.

Furthermore, the current interface of translating time-dependent Hamiltonians into static operators is a bit idiosyncratic; this has hindered adoption by users that already have their own time propagation methods, e.g. OpenQuantumTools.

The package can be generalized by introducing a "Propagator" interface somewhat similar to the Integrator interface in DifferentialEquations.jl – in fact, a Propagator could very directly wrap an ODE Integrator when doing time propagation via a generic ODE solver for time-continous controls.

As the main user interface of QuantumPropagtors, we would have the following:

• propagator = init_prop(state, generator, tlist; method=:auto, backward=false, inplace=true, kwargs...) creates a method-dependent Propagator object that encapsulates the dynamics of the initial state on a time grid tlist under a generator (i.e., a time-dependent Hamiltonian/Liouvillian). The default structure for generator is a nested tuple, e.g. (H₀, (H₁, ϵ(t)) inspired by QuTiP, but it will be easy to have arbitrary custom types as the generator

• Every propagator has at least the properties state (read-only: the current state), t (the current time), and parameters. For piecewise-constant methods, parameters will be the values of the control field (ϵ(t)) discretized to the midpoints of tlist, but in general, this might depend both on the propagation method and the type of generator/controls.

• Each Propagator type is a subtype of AbstractPropagator, and may be a subtype of PiecewisePropagator (a propagator that takes has one parameter per time grid interval, but isn't necessarily constant on the interval – e.g. allowing for linear/spline interpolation), or PWCPropagator for propagators that are piecewise constant in addition to that. We will have ChebyPropagator, NewtonPropagator, and ExpPropagator for the current (PWC) functionality, and further propagators for methods to be added in the future. It should be relatively straightforward to implement a custom Propagator type. Beyond implementing a custom Propagator, a user will not normally have to worry about the internals of the Propagator object, but will interact with it only via the minimal API defined below.

• Propagation happens by calling propstep!(propagator), which advances the propagator to the next point in tlist, either forward or backward. For convenience, it returns the resulting propagator.state or nothing if tlist has been exhausted.

• A new propagation can be initialized with reinitprop!(propagator, state; kwargs...), which resets the propagator to the initial or final time and resets the initial state.

• The propagator.state can be mutated by calling set_state!(propagator, state). This is useful e.g. for removing population from forbidden levels (in lieu of a non-Hermitian Hamiltonian), for implementing quantum jumps in MCWF, or for the propagation of the extended gradient-vector (GradVector)

• The parameters can be mutated freely, allowing e.g. sequential updates in Krotov's method without overhead

The management of piecewise constant generators (i.e., converting time-dependent generators into static operators in each time interval) will be handled inside the propagator object. Right now, this is done inside the GRAPE/Krotov implementations, so moving the functionality into QuantumPropagators will reduce code duplication and simplify the optimal control code.

@seba-car @alastair-marshall Any thoughts on this? Does this seem like it will enable anything one might want to do in the future (like implementing CRAB/GROUP/GOAT?). I might still have to work through a concrete example for how parameters will be used in the non-piecewise case. Presumably, the CRAB coefficients would be stored in parameters, but then what happens if I still want to use a PWC-propagator like :cheby for a CRAB pulse?

@naezzell Does this seem like something you could work with?

[alastair-marshall]

I like the idea of moving closer to the DiffEq.jl interface, as they've put a lot of thought into making everything extensible.

A few questions spring to mind:

but it will be easy to have arbitrary custom types as the generator

will we define that interface somewhere/somehow try to enforce it? so someone can check if their custom type will work by just running some tests on it? Will we also define the interface for the custom Propagator type somewhere? Maybe we can have an unexported function inside this package test_proapagtor_interface that a user can call on their custom type to check that they've implemented everything properly.

Do the propagators not also need at least to store the generators/Hamiltonians too?

Presumably, the CRAB coefficients would be stored in parameters, but then what happens if I still want to use a PWC-propagator like :cheby for a CRAB pulse?

if we made parameters a custom type (maybe overkill) or a dict, we can dispatch if its for CRAB to an intermediate function that does the discretisation before the evolution?

[goerz]

but it will be easy to have arbitrary custom types as the generator

will we define that interface somewhere/somehow try to enforce it? so someone can check if their custom type will work by just running some tests on it?

Yes. The interface for the generator is actually quite straightforward. A custom type will have to specify the methods in the new QuantumPropagators.Controls modules, specifically:

• getcontrols
• evalcontrols / evalcontrols!
• getcontrolderiv / getcontrolderivs
• substitute_controls

A generator may also contain arbitrary custom types as controls (the objects returned by getcontrols). For example, we'll definitely want a CRABcontrol with the spectral coefficients as the natural parameters.

These at the very least have to define the other routines from QuantumPropagators.Controls:

• discretize
• discretize_on_midpoints
• get_control_parameters
• (to be added when we implement the ODE Propagator): eval_control_parameters — the inverse to get_control_parameters

What's more tricky is that e.g. evalcontrols converts a time-dependent generator into a static operator which can also be of an arbitrary custom type. That one is less well-defined. It mostly has to implement various linear algebra operations, e.g. "matrix-vector multiplication" with the state (which may also be of a custom type). It very much depends on how exactly a given propagator is implemented.

For now, the deal with custom types is "Set up your Objective with the custom types, run the optimization, and fix any MethodError that gets thrown".

That being said, I would like to add some checkers along the lines of what you describe to make things a little more user-friendly. Those would be run when the optimization initializes, or manually. I'll keep the implementation of this for a bit later

Will we also define the interface for the custom Propagator type somewhere? Maybe we can have an unexported function inside this package test_proapagtor_interface that a user can call on their custom type to check that they've implemented everything properly.

Yes, same thing. The interface is described in full detail in the documentation now, but having an automatic checker to make it more user-friendly would be something I'd consider in the future.

Do the propagators not also need at least to store the generators/Hamiltonians too?

Internally, probably, but not as a public property. In fact, I've taken active measures to prevent anyone from accessing the generator. All hell would break loose if someone were to mutate the generator after the propagator has been initialized. The mutable part of it (the pulse parameters or amplitudes) is exposed through the mutable parameters property only.

Presumably, the CRAB coefficients would be stored in parameters, but then what happens if I still want to use a PWC-propagator like :cheby for a CRAB pulse?

if we made parameters a custom type (maybe overkill) or a dict, we can dispatch if its for CRAB to an intermediate function that does the discretisation before the evolution?

To be clear, the parameters are always going to be a dict control => Vector{Float64}. There's a lot of room for dispatching on custom types (in the function in QuantumPropagators.Controls that I list above).

I've also given a good thinking to what happens with a pulse with custom parameters like a CRABcontrol in the various optimization/propagator combinations, and I'm pretty sure everything works out (or, can be dealt with in the future):

• Piecewise(PW)-optimization (GRAPE/Krotov) with an ODE Propagator — no problem: GRAPE/Krotov force the "parametrization" of discretizing to the time grid: get_control_parameters and eval_control_parameters aren't even called. The optimized pulses are just going to be amplitudes, of course. The ODE Propagator might be less efficient than Cheby/Newton, but it will work fine
• PW-optimization with PW-propagator (Cheby/Newton) — no problem: same as above
• gradient-free optimization (CRAB) with ODE Propagator — no problem, this is the standard case
• gradient-free optimization with PW-propagator — no problem if using propagate internally, a small problem if using a propagator/propstep! for more efficiency. Can be solved with WrappedPiecewisePropagator (see below)
• Full-AD optimization with ODE Propagator — no problem, standard case
• Full-AD with PW-propagator — no problem if using propagate, use a WrappedPiecewisePropagator otherwise
• GOAT/GROUP ("Schirmer Gradient") with ODE Propagator — no problem, standard case
• GOAT/GROUP with PW-propagator — no problem if using propagate, use a WrappedPiecewisePropagator otherwise

So the situation where there's a slight problem is whenever you want to use a PW-propagator (which requires that the parameters in the propagator are the amplitudes of the controls at the midpoints of the time grid) with an optimization method that acts on the natural parameters of the control, e.g. the spectral coefficients for a CRABcontrol. All these methods require only forward propagation, so if you use the propagate function internally, everything is fine. However, that means re-initializing the propagation again for every iteration of the control algorithm, so it's a bit inefficient. If you want to reuse an existing propagator object, you have to somehow translate between the "natural parameters" and the "pulse discretization parameters". In most cases, maybe you could do that as part of reinitprop!, but for automatic differentiation, you'd actually need the conversion to be traceable. So I think the most elegant way to do this would be to define a WrappedPiecewisePropagator that wraps around an existing PiecewisePropagator and uses a custom AbstractVector in the parameters for the PiecewisePropagator. Any access to the values in that vector dynamically refers back to the original "natural" parameters, e.g. the spectral coefficients for CRAB. I'm pretty sure this will work quite well, but we can revisit it when we get there.

[goerz]

So I think I'm pretty happy with #27 at this point. I'd like to merge and release this relatively quickly (maybe tomorrow), to make it easier for people to work on top of the breaking change.

The next step would be to add the ODE Propagator that wraps DifferentialEquations, but let's push the basic interface out first.

Have a look at the zipped documentation for QuantumPropagators and QuantumControl, in terms of if that explains the interface sufficiently (the "Overview" and "API')

• QuantumPropagators.zip
• QuantumControl.zip