QuantumControlBase

A full control problem with multiple trajectories.

ControlProblem(
   trajectories,
   tlist;
   kwargs...
)

The trajectories are a list of Trajectory instances, each defining an initial state and a dynamical generator for the evolution of that state. Usually, the trajectory will also include a target state (see Trajectory) and possibly a weight. The trajectories may also be given together with tlist as a mandatory keyword argument.

The tlist is the time grid on which the time evolution of the initial states of each trajectory should be propagated. It may also be given as a (mandatory) keyword argument.

The remaining kwargs are keyword arguments that are passed directly to the optimal control method. These typically include e.g. the optimization functional.

The control problem is solved by finding a set of controls that minimize an optimization functional over all trajectories.

Description of a state's time evolution.

Trajectory(
    initial_state,
    generator;
    target_state=nothing,
    weight=1.0,
    kwargs...
)

describes the time evolution of the initial_state under a time-dependent dynamical generator (e.g., a Hamiltonian or Liouvillian).

Trajectories are central to quantum control problems: an optimization functional depends on the result of propagating one or more trajectories. For example, when optimizing for a quantum gate, the optimization considers the trajectories of all logical basis states.

In addition to the initial_state and generator, a Trajectory may include data relevant to the propagation and to evaluating a particular optimization functional. Most functionals have the notion of a "target state" that the initial_state should evolve towards, which can be given as the target_state keyword argument. In some functionals, different trajectories enter with different weights [1], which can be given as a weight keyword argument. Any other keyword arguments are also available to a functional as properties of the Trajectory .

A Trajectory can also be instantiated using all keyword arguments.

Properties

All keyword arguments used in the instantiation are available as properties of the Trajectory. At a minimum, this includes initial_state, generator, target_state, and weight.

By convention, properties with a prop_ prefix, e.g., prop_method, will be taken into account when propagating the trajectory. See propagate_trajectory for details.

Construct the adjoint of a Trajectory.

adj_trajectory = adjoint(trajectory)

The adjoint trajectory contains the adjoint of the dynamical generator traj.generator. All other fields contain a copy of the original field value.

The primary purpose of this adjoint is to facilitate the backward propagation under the adjoint generator that is central to gradient-based optimization methods such as GRAPE and Krotov's method.

Combine multiple info_hook functions.

chain_infohooks(funcs...)

combines funcs into a single Function that can be passes as info_hook to ControlProblem or any optimize-function.

Each function in func must be a suitable info_hook by itself. This means that it should receive the optimization workspace object as its first positional parameter, then positional parameters specific to the optimization method, and then an arbitrary number of data parameters. It must return either nothing or a tuple of "info" objects (which will end up in the records field of the optimization result).

When chaining infohooks, the funcs will be called in series, and the "info" objects will be accumulated into a single result tuple. The combined results from previous funcs will be given to the subsequent funcs as data parameters. This allows for the infohooks in the chain to communicate.

The chain will return the final combined result tuple, or nothing if all funcs return nothing.

Check an amplitude in a Generator in the context of optimal control.

@test check_amplitude(
    ampl; tlist, for_gradient_optimization=true, quiet=false
)

verifies that the given ampl is a valid element in the list of amplitudes of a Generator object. This checks all the conditions of QuantumPropagators.Interfaces.check_amplitude. In addition, the following conditions must be met.

If for_gradient_optimization:

• The function get_control_deriv(ampl, control) must be defined
• If ampl does not depend on control, get_control_deriv(ampl, control) must return 0.0
• If ampl depends on control, get_control_deriv(ampl, control) must return an object u so that evaluate(u, tlist, n) returns a Number. In most cases, u itself will be a Number.

The function returns true for a valid amplitude and false for an invalid amplitude. Unless quiet=true, it will log an error to indicate which of the conditions failed.

Check the dynamical generator in the context of optimal control.

@test check_generator(
    generator; state, tlist,
    for_mutable_operator=true, for_immutable_operator=true,
    for_mutable_state=true, for_immutable_state=true,
    for_expval=true, for_pwc=true, for_time_continuous=false,
    for_parameterization=false, for_gradient_optimization=true,
    atol=1e-15, quiet=false
)

verifies the given generator. This checks all the conditions of QuantumPropagators.Interfaces.check_generator. In addition, the following conditions must be met.

If for_gradient_optimization:

• get_control_derivs(generator, controls) must be defined and return a vector containing the result of get_control_deriv(generator, control) for every control in controls.
• get_control_deriv(generator, control) must return an object that passes the less restrictive QuantumPropagators.Interfaces.check_generator if control is in get_controls(generator).
• get_control_deriv(generator, control) must return nothing if control is not in get_controls(generator)
• If generator is a Generator instance, every ampl in generator.amplitudes must pass check_amplitude(ampl; tlist).

The function returns true for a valid generator and false for an invalid generator. Unless quiet=true, it will log an error to indicate which of the conditions failed.

a = get_control_deriv(ampl, control)

returns the derivative $∂a_l(t)/∂ϵ_{l'}(t)$ of the given amplitude $a_l(\{ϵ_{l''}(t)\}, t)$ with respect to the given control $ϵ_{l'}(t)$. For "trivial" amplitudes, where $a_l(t) ≡ ϵ_l(t)$, the result with be either 1.0 or 0.0 (depending on whether ampl ≡ control). For non-trivial amplitudes, the result may be another amplitude that depends on the controls and potentially on time, but can be evaluated to a constant with evaluate.

Get the derivative of the generator $G$ w.r.t. the control $ϵ(t)$.

μ  = get_control_deriv(generator, control)

returns nothing if the generator (Hamiltonian or Liouvillian) does not depend on control, or a generator

\[μ = \frac{∂G}{∂ϵ(t)}\]

otherwise. For linear control terms, μ will be a static operator, e.g. an AbstractMatrix or an Operator. For non-linear controls, μ will be time-dependent, e.g. a Generator. In either case, evaluate should be used to evaluate μ into a constant operator for particular values of the controls and a particular point in time.

For constant generators, e.g. an Operator, the result is always nothing.

Get a vector of the derivatives of generator w.r.t. each control.

get_control_derivs(generator, controls)

return as vector containing the derivative of generator with respect to each control in controls. The elements of the vector are either nothing if generator does not depend on that particular control, or a function μ(α) that evaluates the derivative for a particular value of the control, see get_control_deriv.

Initialize a propagator for a given Trajectory.

propagator = init_prop_trajectory(
    traj,
    tlist;
    initial_state=traj.initial_state,
    kwargs...
)

initializes a Propagator for the propagation of the initial_state under the dynamics described by traj.generator.

All keyword arguments are forwarded to QuantumPropagators.init_prop, with default values from any property of traj with a prop_ prefix. That is, the keyword arguments for the underlying QuantumPropagators.init_prop are determined as follows:

• For any property of traj whose name starts with the prefix prop_, strip the prefix and use that property as a keyword argument for init_prop. For example, if traj.prop_method is defined, method=traj.prop_method will be passed to init_prop. Similarly, traj.prop_inplace would be passed as inplace=traj.prop_inplace, etc.
• Any explicitly keyword argument to init_prop_trajectory overrides the values from the properties of traj.

Note that the propagation method in particular must be specified, as it is a mandatory keyword argument in QuantumPropagators.propagate). Thus, either traj must have a property prop_method of the trajectory, or method must be given as an explicit keyword argument.

Return a function that evaluates $|χ_k⟩ = -∂J_T/∂⟨ϕ_k|$.

chi! = make_chi(
    J_T,
    trajectories;
    mode=:any,
    automatic=:default,
    via=(any(isnothing(t.target_state) for t in trajectories) ? :phi : :tau),
)

creates a function chi!(χ, ϕ, trajectories; τ) that sets the k'th element of χ to $|χ_k⟩ = -∂J_T/∂⟨ϕ_k|$, where $|ϕ_k⟩$ is the k'th element of ϕ. These are the states used as the boundary condition for the backward propagation propagation in Krotov's method and GRAPE. Each $|χₖ⟩$ is defined as a matrix calculus Wirtinger derivative,

\[|χ_k(T)⟩ = -\frac{∂J_T}{∂⟨ϕ_k|} = -\frac{1}{2} ∇_{ϕ_k} J_T\,;\qquad ∇_{ϕ_k} J_T ≡ \frac{∂J_T}{\Re[ϕ_k]} + i \frac{∂J_T}{\Im[ϕ_k]}\,.\]

The function J_T must take a vector of states ϕ and a vector of trajectories as positional parameters, and a vector τ as a keyword argument, see e.g. J_T_sm). If all trajectories define a target_state, then τ will be the overlap of the states ϕ with those target states. The functional J_T may or may not use those overlaps. Likewise, the resulting chi! may or may not use the keyword parameter τ.

The derivative can be calculated analytically of automatically (via automatic differentiation) depending on the value of mode. For mode=:any, an analytic derivative is returned if available, with a fallback to an automatic derivative.

If mode=:analytic, return an analytically known $-∂J_T/∂⟨ϕ_k|$, e.g.,

• J_T_sm → chi_sm!,
• J_T_re → chi_re!,
• J_T_ss → chi_ss!.

and throw an error if no analytic derivative is known.

If mode=:automatic, return an automatic derivative (even if an analytic derivative is known). The calculation of an automatic derivative (whether via mode=:any or mode=:automatic) requires that a suitable framework (e.g., Zygote or FiniteDifferences) has been loaded. The loaded module must be passed as automatic keyword argument. Alternatively, it can be registered as a default value for automatic by calling QuantumControl.set_default_ad_framework.

When evaluating $|χ_k⟩$ automatically, if via=:phi is given , $|χ_k(T)⟩$ is calculated directly as defined above from the gradient with respect to the states $\{|ϕ_k(T)⟩\}$. The resulting function chi! ignores any passed τ keyword argument.

If via=:tau is given instead, the functional $J_T$ is considered a function of overlaps $τ_k = ⟨ϕ_k^\tgt|ϕ_k(T)⟩$. This requires that all trajectories define a target_state and that J_T calculates the value of the functional solely based on the values of τ passed as a keyword argument. With only the complex conjugate $τ̄_k = ⟨ϕ_k(T)|ϕ_k^\tgt⟩$ having an explicit dependency on $⟨ϕ_k(T)|$, the chain rule in this case is

\[|χ_k(T)⟩ = -\frac{∂J_T}{∂⟨ϕ_k|} = -\left( \frac{∂J_T}{∂τ̄_k} \frac{∂τ̄_k}{∂⟨ϕ_k|} \right) = - \frac{1}{2} (∇_{τ_k} J_T) |ϕ_k^\tgt⟩\,.\]

Again, we have used the definition of the Wirtinger derivatives,

\[\begin{align*} \frac{∂J_T}{∂τ_k} &≡ \frac{1}{2}\left( \frac{∂ J_T}{∂ \Re[τ_k]} - i \frac{∂ J_T}{∂ \Im[τ_k]} \right)\,,\\ \frac{∂J_T}{∂τ̄_k} &≡ \frac{1}{2}\left( \frac{∂ J_T}{∂ \Re[τ_k]} + i \frac{∂ J_T}{∂ \Im[τ_k]} \right)\,, \end{align*}\]

and the definition of the Zygote gradient with respect to a complex scalar,

\[∇_{τ_k} J_T = \left( \frac{∂ J_T}{∂ \Re[τ_k]} + i \frac{∂ J_T}{∂ \Im[τ_k]} \right)\,.\]

QuantumControlBase.make_analytic_chi(::typeof(J_T_sm), trajectories) = chi_sm!

Return a function to evaluate $∂J_a/∂ϵ_{ln}$ for a pulse value running cost.

grad_J_a! = make_grad_J_a(
    J_a,
    tlist;
    mode=:any,
    automatic=:default,
)

returns a function so that grad_J_a!(∇J_a, pulsevals, tlist) sets $∂J_a/∂ϵ_{ln}$ as the elements of the (vectorized) ∇J_a. The function J_a must have the interface J_a(pulsevals, tlist), see, e.g., J_a_fluence.

The parameters mode and automatic are handled as in make_chi, where mode is one of :any, :analytic, :automatic, and automatic is he loaded module of an automatic differentiation framework, where :default refers to the framework set with QuantumControl.set_default_ad_framework.

make_analytic_grad_J_a(::typeof(J_a_fluence), tlist) = grad_J_a_fluence!

Optimize a quantum control problem.

result = optimize(problem; method, check=true, kwargs...)

optimizes towards a solution of given problem with the given method, which should be a Module implementing the method, e.g.,

using Krotov
result = optimize(problem; method=Krotov)

Note that method is a mandatory keyword argument.

If check is true (default), the initial_state and generator of each trajectory is checked with check_state and check_generator. Any other keyword argument temporarily overrides the corresponding keyword argument in problem. These arguments are available to the optimizer, see each optimization package's documentation for details.

To obtain the documentation for which options a particular method uses, run, e.g.,

? optimize(problem, ::Val{:Krotov})

where :Krotov is the name of the module implementing the method. The above is also the method signature that a Module wishing to implement a control method must define.

The returned result object is specific to the optimization method.

Propagate multiple trajectories in parallel.

result = propagate_trajectories(
    trajectories, tlist; use_threads=true, kwargs...
)

runs propagate_trajectory for every trajectory in trajectories, collects and returns a vector of results. The propagation happens in parallel if use_threads=true (default). All keyword parameters are passed to propagate_trajectory, except that if initial_state is given, it must be a vector of initial states, one for each trajectory. Likewise, to pass pre-allocated storage arrays to storage, a vector of storage arrays must be passed. A simple storage=true will still work to return a vector of storage results.

Propagate a Trajectory.

propagate_trajectory(
    traj,
    tlist;
    initial_state=traj.initial_state,
    kwargs...
)

propagates initial_state under the dynamics described by traj.generator. It takes the same keyword arguments as QuantumPropagators.propagate, with default values from any property of traj with a prop_ prefix (prop_method, prop_inplace, prop_callback, …). See init_prop_trajectory for details.

Note that method (a mandatory keyword argument in QuantumPropagators.propagate) must be specified, either as a property prop_method of the trajectory, or by passing a method keyword argument explicitly.

Register a callback to dump a running optimization to disk on unexpected exit.

A long-running optimization routine may use

if !isnothing(atexit_filename)
    set_atexit_save_optimization(
        atexit_filename, result; msg_property=:message, msg="Abort: ATEXIT"
    )
    # ...
    popfirst!(Base.atexit_hooks)  # remove callback
end

to register a callback that writes the given result object to the given filename in JLD2 format in the event that the program terminates unexpectedly. The idea is to avoid data loss if the user presses CTRL-C in a non-interactive program (SIGINT), or if the process receives a SIGTERM from an HPC scheduler because the process has reached its allocated runtime limit. Note that the callback cannot protect against data loss in all possible scenarios, e.g., a SIGKILL will terminate the program without giving the callback a chance to run (as will yanking the power cord).

As in the above example, the optimization routine should make set_atexit_save_optimization conditional on an atexit_filename keyword argument, which is what QuantumControl.@optimize_or_load will pass to the optimization routine. The optimization routine must remove the callback from Base.atexit_hooks when it exits normally. Note that in an interactive context, CTRL-C will throw an InterruptException, but not cause a shutdown. Optimization routines that want to prevent data loss in this situation should handle the InterruptException and return result, in addition to using set_atexit_save_optimization.

If msg_property is not nothing, the given msg string will be stored in the corresponding property of the (mutable) result object before it is written out.

The resulting JLD2 file is compatible with QuantumControl.load_optimization.

controls = get_controls(problem)

extracts the controls from problem.trajectories.

controls = get_controls(trajectories)

extracts the controls from a list of trajectories (i.e., from each trajectory's generator). Controls that occur multiple times in the different trajectories will occur only once in the result.

parameters = get_parameters(problem)

extracts the parameters from problem.trajectories.

parameters = get_parameters(trajectories)

collects and combines get parameter arrays from all the generators in trajectories. Note that this allows any custom generator type to define a custom get_parameters method to override the default of obtaining the parameters recursively from the controls inside the generator.

trajectory = substitute(trajectory::Trajectory, replacements)
trajectories = substitute(trajectories::Vector{<:Trajectory}, replacements)

recursively substitutes the initial_state, generator, and target_state.

Conditionally apply multi-threading to for loops.

This is a variation on Base.Threads.@threads that adds a run-time boolean flag to enable or disable threading. It is intended for internal use in packages building on QuantumControlBase.

Usage:

using QuantumControlBase: @threadsif

function optimize(trajectories; use_threads=true)
    @threadsif use_threads for k = 1:length(trajectories)
    # ...
    end
end