QuantumControlBase Package

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 [20], 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.

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.

using Krotov
result = optimize(problem; method=Krotov, kwargs...)

optimizes the given control problem using Krotov's method, returning a KrotovResult.

Keyword arguments that control the optimization are taken from the keyword arguments used in the instantiation of problem; any of these can be overriden with explicit keyword arguments to optimize.

Required problem keyword arguments

• J_T: A function J_T(ϕ, trajectories) that evaluates the final time functional from a list ϕ of forward-propagated states and problem.trajectories.

Recommended problem keyword arguments

• lambda_a=1.0: The inverse Krotov step width λ_a for every pulse.
• update_shape=(t->1.0): A function S(t) for the "update shape" that scales the update for every pulse

If different controls require different lambda_a or update_shape, a dict pulse_options must be given instead of a global lambda_a and update_shape, see below.

Optional problem keyword arguments

The following keyword arguments are supported (with default values):

• pulse_options: A dictionary that maps every control (as obtained by get_controls from the problem.trajectories) to the following dict:

  • :lambda_a: The value for inverse Krotov step width λₐ
  • :update_shape: A function S(t) for the "update shape" that scales the Krotov pulse update.

  This overrides the global lambda_a and update_shape arguments.

• chi: A function chi!(χ, ϕ, trajectories) what receives a list ϕ of the forward propagated states and must set $|χₖ⟩ = -∂J_T/∂⟨ϕₖ|$. If not given, it will be automatically determined from J_T via make_chi with the default parameters.

• sigma=nothing: Function that calculate the second-order contribution. If not given, the first-order Krotov method is used.

• iter_start=0: the initial iteration number

• iter_stop=5000: the maximum iteration number

• prop_method: The propagation method to use for each trajectory, see below.

• update_hook: A function that receives the Krotov workspace, the iteration number, the list of updated pulses and the list of guess pulses as positional arguments. The function may mutate any of its arguments. This may be used e.g. to apply a spectral filter to the updated pulses or to perform similar manipulations.

• info_hook: A function (or tuple of functions) that receives the same arguments as update_hook, in order to write information about the current iteration to the screen or to a file. The default info_hook prints a table with convergence information to the screen. Runs after update_hook. The info_hook function may return a tuple, which is stored in the list of records inside the KrotovResult object.

• check_convergence: a function to check whether convergence has been reached. Receives a KrotovResult object result, and should set result.converged to true and result.message to an appropriate string in case of convergence. Multiple convergence checks can be performed by chaining functions with ∘. The convergence check is performed after any calls to update_hook and info_hook.

• verbose=false: If true, print information during initialization

Trajectory propagation

Krotov's method involves the forward and backward propagation for every Trajectory in the problem. The keyword arguments for each propagation (see propagate) are determined from any properties of each Trajectory that have a prop_ prefix, cf. init_prop_trajectory.

In situations where different parameters are required for the forward and backward propagation, instead of the prop_ prefix, the fw_prop_ and bw_prop_ prefix can be used, respectively. These override any setting with the prop_ prefix. This applies both to the properties of each Trajectory and the problem keyword arguments.

Note that the propagation method for each propagation must be specified. In most cases, it is sufficient (and recommended) to pass a global prop_method problem keyword argument.

using GRAPE
result = optimize(problem; method=GRAPE, kwargs...)

optimizes the given control problem via the GRAPE method, by minimizing the functional

\[J(\{ϵ_{nl}\}) = J_T(\{|ϕ_k(T)⟩\}) + λ_a J_a(\{ϵ_{nl}\})\]

where the final time functional $J_T$ depends explicitly on the forward-propagated states and the running cost $J_a$ depends explicitly on pulse values $ϵ_{nl}$ of the l'th control discretized on the n'th interval of the time grid.

Returns a GrapeResult.

Keyword arguments that control the optimization are taken from the keyword arguments used in the instantiation of problem; any of these can be overridden with explicit keyword arguments to optimize.

Required problem keyword arguments

• J_T: A function J_T(ϕ, trajectories; τ=τ) that evaluates the final time functional from a vector ϕ of forward-propagated states and problem.trajectories. For all trajectories that define a target_state, the element τₖ of the vector τ will contain the overlap of the state ϕₖ with the target_state of the k'th trajectory, or NaN otherwise.

Optional problem keyword arguments

• chi: A function chi!(χ, ϕ, trajectories) what receives a list ϕ of the forward propagated states and must set $|χₖ⟩ = -∂J_T/∂⟨ϕₖ|$. If not given, it will be automatically determined from J_T via make_chi with the default parameters.

• J_a: A function J_a(pulsevals, tlist) that evaluates running costs over the pulse values, where pulsevals are the vectorized values $ϵ_{nl}$, where n are in indices of the time intervals and l are the indices over the controls, i.e., [ϵ₁₁, ϵ₂₁, …, ϵ₁₂, ϵ₂₂, …] (the pulse values for each control are contiguous). If not given, the optimization will not include a running cost.

• gradient_method=:gradgen: One of :gradgen (default) or :taylor. With gradient_method=:gradgen, the gradient is calculated using QuantumGradientGenerators. With gradient_method=:taylor, it is evaluated via a Taylor series, see Eq. (20) in Kuprov and Rogers, J. Chem. Phys. 131, 234108 (2009) [17].

• taylor_grad_max_order=100: If given with gradient_method=:taylor, the maximum number of terms in the Taylor series. If taylor_grad_check_convergence=true (default), if the Taylor series does not convergence within the given number of terms, throw an an error. With taylor_grad_check_convergence=true, this is the exact order of the Taylor series.

• taylor_grad_tolerance=1e-16: If given with gradient_method=:taylor and taylor_grad_check_convergence=true, stop the Taylor series when the norm of the term falls below the given tolerance. Ignored if taylor_grad_check_convergence=false.

• taylor_grad_check_convergence=true: If given as true (default), check the convergence after each term in the Taylor series an stop as soon as the norm of the term drops below the given number. If false, stop after exactly taylor_grad_max_order terms.

• lambda_a=1: A weight for the running cost J_a.

• grad_J_a: A function to calculate the gradient of J_a. If not given, it will be automatically determined.

• upper_bound: An upper bound for the value of any optimized control. Time-dependent upper bounds can be specified via pulse_options.

• lower_bound: A lower bound for the value of any optimized control. Time-dependent lower bounds can be specified via pulse_options.

• pulse_options: A dictionary that maps every control (as obtained by get_controls from the problem.trajectories) to a dict with the following possible keys:

  • :upper_bounds: A vector of upper bound values, one for each intervals of the time grid. Values of Inf indicate an unconstrained upper bound for that time interval, respectively the global upper_bound, if given.
  • :lower_bounds: A vector of lower bound values. Values of -Inf indicate an unconstrained lower bound for that time interval,

• update_hook: Not implemented

• info_hook: A function (or tuple of functions) that receives the same arguments as update_hook, in order to write information about the current iteration to the screen or to a file. The default info_hook prints a table with convergence information to the screen. Runs after update_hook. The info_hook function may return a tuple, which is stored in the list of records inside the GrapeResult object.

• check_convergence: A function to check whether convergence has been reached. Receives a GrapeResult object result, and should set result.converged to true and result.message to an appropriate string in case of convergence. Multiple convergence checks can be performed by chaining functions with ∘. The convergence check is performed after any calls to update_hook and info_hook.

• x_tol: Parameter for Optim.jl

• f_tol: Parameter for Optim.jl

• g_tol: Parameter for Optim.jl

• show_trace: Parameter for Optim.jl

• extended_trace: Parameter for Optim.jl

• show_every: Parameter for Optim.jl

• allow_f_increases: Parameter for Optim.jl

• optimizer: An optional Optim.jl optimizer (Optim.AbstractOptimizer instance). If not given, an L-BFGS-B optimizer will be used.

• prop_method: The propagation method to use for each trajectory, see below.

• verbose=false: If true, print information during initialization

Trajectory propagation

GRAPE may involve three types of propagation:

• A forward propagation for every Trajectory in the problem
• A backward propagation for every trajectory
• A backward propagation of a gradient generator for every trajectory.

The keyword arguments for each propagation (see propagate) are determined from any properties of each Trajectory that have a prop_ prefix, cf. init_prop_trajectory.

In situations where different parameters are required for the forward and backward propagation, instead of the prop_ prefix, the fw_prop_ and bw_prop_ prefix can be used, respectively. These override any setting with the prop_ prefix. Similarly, properties for the backward propagation of the gradient generators can be set with properties that have a grad_prop_ prefix. These prefixes apply both to the properties of each Trajectory and the problem keyword arguments.

Note that the propagation method for each propagation must be specified. In most cases, it is sufficient (and recommended) to pass a global prop_method problem keyword argument.

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.

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.

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!

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.

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.

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

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.

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.

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.

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!

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.