QuantumControlBase Package

Run optimize and store the result, or load the result if it exists.

result = @optimize_or_load(
    file,
    problem;
    method,
    suffix="jld2",
    tag=DrWatson.readenv("DRWATSON_TAG", true),
    gitpath=DrWatson.projectdir(),
    storepatch::Bool=DrWatson.readenv("DRWATSON_STOREPATCH", false),
    force=false,
    verbose=true,
    wsave_kwargs=Dict(),
    metadata=nothing,
    kwargs...
)

runs result = optimize(problem; method, kwargs...) and stores result in file. Note that the method keyword argument is mandatory. In addition to the result, the data in the output file may also contain some metadata, e.g. (automatically) "gitcommit" containing the git commit hash of the project that produced the file, and "script" with the file name and line number where @optimize_or_load was called, see load_optimization. If metadata is given as a dict on input, the data it contains will be included in the output file.

If file already exists (and force=false), load the result from that file instead of running the optimization.

The @optimize_or_load macro is intended to integrate well with the DrWatson framework for scientific projects and utilizes several configuration options and utility functions from DrWatson, see below. Note that even though DrWatson is recomended, you are not required to use if for your projects in order to use @optimize_or_load or any other part of QuantumControl.

I/O Keywords

The following keyword arguments determine how the result is stored:

• suffix. File extension of file, determining the output data format (see DrWatson Saving Tools). If file does not end with the given extension, it will be appended.
• tag: Whether to record the current "gitcommit" as metadata alongside the optimization result, via DrWatson.tagsave. If not given explicitly, determine automatically from suffix.
• gitpath, storepatch: Passed to DrWatson.tagsave if tag is true.
• force: If true, run and store the optimization regardless of whether file already exists.
• verbose: If true, print info about the process
• wsave_kwargs: Additional keyword arguments to pass to DrWatson.wsave, e.g., to enable compression

All other keyword arguments are passed directly to optimize.

Related Functions

• DrWatson.@produce_or_load: The lower-level backend implementing the functionality of @optimize_or_load.
• load_optimization: Function to load a file produced by @optimize_or_load

A full control problem with multiple objectives.

ControlProblem(;
   objectives,
   tlist,
   kwargs...
)

Note that the control problem can only be instantiated via keyword arguments.

The objectives are a list of Objective instances, each defining an initial state and a dynamical generator for the evolution of that state. Usually, the objective will also include a target state (see Objective) and possibly a weight.

The tlist is the time grid on which the time evolution of the initial states of each objective should be propagated.

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 simultaneously fulfill all objectives.

Optimization objective.

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

describes an optimization objective that is tracked by the time evolution of the given initial_state under the given generator, e.g., a time-dependent Hamiltonian or Liouvillian. Each objective represents a single propagated state on which an optimization functional may depend.

The most common control problems in quantum control (state-to-state, gate optimization) require that the initial_state evolves into a target_state, which should be given as a keyword argument.

An optimization functional usually depends on multiple forward-propagated states (i.e., multiple objectives). Sometimes, it is useful to weight the contributions of different objectives relative to each other, see, e.g., Goerz et al., New J. Phys. 16, 055012 (2014). To this end, a weight can be attached to each Objective as an optional keyword argument.

Any other keyword arguments are available to a custom functional as properties of the Objective .

Note that the Objective can only be instantiated via keyword arguments, with initial_state and generator being the only two mandatory keyword arguments.

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.

Construct a Liouvillian super-operator.

ℒ = liouvillian(Ĥ, c_ops=(); convention=:LvN)

calculates the sparse Liouvillian super-operator ℒ from the Hamiltonian Ĥ and a list c_ops of Lindblad operators.

With convention=:LvN, applying the resulting ℒ to a vectorized density matrix ρ⃗ calculates $\frac{d}{dt} \vec{\rho}(t) = ℒ \vec{\rho}(t)$ equivalent to the Liouville-von-Neumann equation for the density matrix $ρ̂$,

\[\frac{d}{dt} ρ̂(t) = -i [Ĥ, ρ̂(t)] + \sum_k\left( Â_k ρ̂ Â_k^\dagger - \frac{1}{2} A_k^\dagger Â_k ρ̂ - \frac{1}{2} ρ̂ Â_k^\dagger Â_k \right)\,,\]

where the Lindblad operators $Â_k$ are the elements of c_ops.

The Hamiltonian $Ĥ$ may be time-dependent, using a nested-tuple format by default, e.g., (Ĥ₀, (H₁, ϵ₁), (H₂, ϵ₂)), where ϵ₁ and ϵ₂ are functions of time. In this case, the resulting ℒ will also be in nested tuple format, ℒ = (ℒ₀, (ℒ₁, ϵ₁), (ℒ₂, ϵ₂)), where the initial element contains the superoperator ℒ₀ for the static component of the Liouvillian, i.e., the commutator with the drift Hamiltonian Ĥ₀, plus the dissipator (sum over $k$), as a sparse matrix. Time-dependent Lindblad operators are not supported. The remaining elements are tuples (ℒ₁, ϵ₁) and (ℒ₂, ϵ₂) corresponding to the commutators with the two control Hamiltonians, where ℒ₁ and ℒ₂ again are sparse matrices.

If $Ĥ$ is not time-dependent, the resulting ℒ will be a single-element tuple containing the Liouvillian as a sparse matrix, ℒ = (ℒ₀, ).

With convention=:TDSE, the Liouvillian will be constructed for the equation of motion $i \hbar \frac{d}{dt} \vec{\rho}(t) = ℒ \vec{\rho}(t)$ to match exactly the form of the time-dependent Schrödinger equation. While this notation is not standard in the literature of open quantum systems, it has the benefit that the resulting ℒ can be used in a numerical propagator for a (non-Hermitian) Schrödinger equation without any change. Thus, for numerical applications, convention=:TDSE is generally preferred. The returned ℒ between the two conventions differs only by a factor of $i$, since we generally assume $\hbar=1$.

The convention keyword argument is mandatory, to force a conscious choice.

See Goerz et. al. "Optimal control theory for a unitary operation under dissipative evolution", arXiv 1312.0111v2, Appendix B.2 for the explicit construction of the Liouvillian superoperator as a sparse matrix.

Load a previously stored optimization.

result = load_optimization(file; verbose=true, kwargs...)

recovers a result previously stored by @optimize_or_load.

result, metadata = load_optimization(file; return_metadata=true, kwargs...)

also obtains a metadata dict containing e.g., "gitcommit" or "script" depending on the options to @optimize_or_load.

Calling load_optimization with verbose=true (default) will show the metadata after loading the file.

Optimize a quantum control problem.

opt_result = optimize(problem; method=<method>, kwargs...)

optimizes towards a solution of given problem with the given optimization method. Any keyword argument temporarily overrides the corresponding keyword argument in problem.

opt_result = optimize(problem; method=:krotov, kwargs...)

optimizes problem using Krotov's method, see Krotov.optimize_krotov.

opt_result = optimize(problem; method=:GRAPE, kwargs...)

optimizes problem using GRadident Ascent Pulse Engineering (GRAPE), see GRAPE.optimize_grape.

Propagate with the dynamical generator of a control objective.

propagate_objective(obj, tlist; method=:auto, initial_state=obj.initial_state,
                    controls_map=IdDict(), kwargs...)

propagates initial_state under the dynamics described by obj.generator.

The optional dict control_map may be given to replace the controls in obj.generator (as obtained by getcontrols) with custom functions or vectors, e.g. with the controls resulting from optimization, see also substitute_controls.

If obj has a property/field prop_method or fw_prop_method, its value will be used as the default for method instead of :auto. An explicit keyword argument for method always overrides the default.

All other kwargs are forwarded to the underlying QuantumPropagators.propagate method for obj.initial_state.

Propagate multiple objectives in parallel.

result = propagate_objectives(objectives, tlist; use_threads=true, kwargs...)

runs propagate_objective for every objective in objectives, collects and returns a vector of results. The propagation happens in parallel if use_threads=true (default). All keyword parameters are passed to propagate_objective, except that if initial_state is given, it must be a vector of initial states, one for each objective. 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.

Extended state-vector for the dynamic gradient.

Ψ̃ = GradVector(Ψ, num_controls)

for an initial state Ψ and num_controls control fields.

The GradVector conceptually corresponds to a direct-sum (block) column-vector $Ψ̃ = (|Ψ̃₁⟩, |Ψ̃₂⟩, … |Ψ̃ₙ⟩, |Ψ⟩)^T$, where $n$ is num_controls. With a matching $G̃$ as in the documentation of TimeDependentGradGenerator, we have

\[G̃ Ψ̃ = \begin{pmatrix} Ĥ |Ψ̃₁⟩ + Ĥ₁|Ψ⟩ \\ \vdots \\ Ĥ |Ψ̃ₙ⟩ + Ĥₙ|Ψ⟩ \\ Ĥ |Ψ⟩ \end{pmatrix}\]

and

\[e^{-i G̃ dt} \begin{pmatrix} 0 \\ \vdots \\ 0 \\ |Ψ⟩ \end{pmatrix} = \begin{pmatrix} \frac{∂}{∂ϵ₁} e^{-i Ĥ dt} |Ψ⟩ \\ \vdots \\ \frac{∂}{∂ϵₙ} e^{-i Ĥ dt} |Ψ⟩ \\ e^{-i Ĥ dt} |Ψ⟩ \end{pmatrix}.\]

Reset the given gradient vector for a new gradient evaluation.

resetgradvec!(Ψ̃::GradVector)

zeroes out Ψ̃.grad_states but leaves Ψ̃.state unaffected.

resetgradvec!(Ψ̃::GradVector, Ψ)

additionally sets Ψ̃.state to Ψ.

Static generator for the standard dynamic gradient.

G::GradGenerator = evalcontrols(G_of_t::TimeDependentGradGenerator, vals_dict)

is the result of plugging in specific values for all controls in a TimeDependentGradGenerator. See evalcontrols and evalcontrols!.

The resulting object can be multiplied directly with a GradVector, e.g., in the process of evaluating a piecewise-constant time propagation.

Extended generator for the standard dynamic gradient.

G̃ = TimeDependentGradGenerator(G)

contains the original time-dependent generator G (a Hamiltonian or Liouvillian) in G̃.G, a vector of control derivatives $∂G/∂ϵₗ(t)$ in G̃.control_derivs, and the controls in G̃.controls.

For a generator $G = Ĥ(t) = Ĥ₀ + ϵ₁(t) Ĥ₁ + … + ϵₙ(t) Ĥₙ$, this extended generator encodes the block-matrix

\[G̃ = \begin{pmatrix} Ĥ(t) & 0 & \dots & 0 & Ĥ₁ \\ 0 & Ĥ(t) & \dots & 0 & Ĥ₂ \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \dots & Ĥ(t) & Ĥₙ \\ 0 & 0 & \dots & 0 & Ĥ(t) \end{pmatrix}\]

Note that the $∂G/∂ϵₗ(t)$ ($Ĥₗ$ in the above example) are functions, to account for the possibility of non-linear control terms, see getcontrolderiv.

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.ConditionalThreads: @threadsif

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

Real-part functional.

J_T_re(ϕ, objectives; τ=nothing)

calculates

\[J_{T,\text{re}} = 1 - F_{\text{re}} \quad\in \begin{cases} [0, 2] & \text{in Hilbert space} \\ [0, 1] & \text{in Liouville space.} \end{cases}\]

All arguments are passed to f_tau while evaluating $F_{\text{re}}$ in F_re.

Square-modulus functional.

J_T_sm(ϕ, objectives; τ=nothing)

calculates

\[J_{T,\text{sm}} = 1 - F_{\text{sm}} \quad\in [0, 1].\]

All arguments are passed to f_tau while evaluating $F_{\text{sm}}$ in F_sm.

State-to-state phase-insensitive functional.

J_T_ss(ϕ, objectives; τ=nothing)

calculates

\[J_{T,\text{ss}} = 1 - F_{\text{ss}} \in [0, 1].\]

All arguments are passed to F_ss.

Running cost for the pulse fluence.

J_a = J_a_fluence(pulsevals, tlist)

calculates

\[J_a = \sum_l \int_0^T |ϵ_l(t)|^2 dt = \left(\sum_{ln} |ϵ_{ln}|^2 \right) dt\]

where $ϵ_{ln}$ are the values in the (vectorized) pulsevals, n is the index of the intervals of the time grid, and $dt$ is the time step, taken from the first time interval of tlist and assumed to be uniform.

Convert a functional from acting on a gate to acting on propagated states.

J_T = gate_functional(J_T_U; kwargs...)

constructs a functional J_T that meets the requirements for for Krotov/GRAPE and make_chi. That is, the output J_T takes positional positional arguments ϕ and objectives. The input functional J_T_U is assumed to have the signature J_T_U(U; kwargs...) where U is a matrix with elements $U_{ij} = ⟨Ψ_i|ϕ_j⟩$, where $|Ψ_i⟩$ is the initial_state of the i'th objectives (assumed to be the i'th canonical basis state) and $|ϕ_j⟩$ is the result of forward-propagating $|Ψ_j⟩$. That is, U is the projection of the time evolution operator into the subspace defined by the basis in the initial_states of the objectives.

See also

• make_gate_chi — create a corresponding chi function that acts more efficiently than the general make_chi.

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

chi! = make_chi(
    J_T,
    objectives;
    force_zygote=false,
    via=(any(isnothing(obj.target_state) for obj in objectives) ? :phi : :tau),
    use_finite_differences=false
)

creates a function chi!(χ, ϕ, objectives; τ) 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 objectives as positional parameters, and a vector τ as a keyword argument, see e.g. J_T_sm. If all objectives 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 τ.

For functionals where $-∂J_T/∂⟨ϕ_k|$ is known analytically, that analytic derivative will be returned, e.g.,

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

Otherwise, or if force_zygote=true or use_finite_differences=true, the derivative to calculate $|χ_k⟩$ will be evaluated automatically, via automatic differentiation with Zygote, or via finite differences (which primarily serves for testing the Zygote gradient).

When evaluating $|χ_k⟩$ automatically, if via=:phi is given , $|χ_k(T)⟩$ is calculated directly as defined a 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 objectives 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)\,.\]

make_analytic_chi(::typeof(J_T_sm), objectives) = chi_sm!

Return a function to evaluate $|χ_k⟩ = -∂J_T(Û)/∂⟨ϕ_k|$ via the chain rule.

chi! = make_gate_chi(J_T_U, objectives; use_finite_differences=false, kwargs...)

returns a function equivalent to

chi! = make_chi(gate_functional(J_T_U; kwargs...), objectives)

\[\begin{split} |χ_k⟩ &= -\frac{∂}{∂⟨ϕ_k|} J_T \\ &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} \frac{∂ U_{ik}}{∂⟨ϕ_k|} \\ &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} |Ψ_i⟩ \end{split}\]

where $|Ψ_i⟩$ is the basis state stored as the initial_state of the i'th objective, see gate_functional.

The gradient $∇_U J_T$ is obtained via automatic differentiation, or via finite differences if use_finite_differences=true.

Compared to the more general make_chi, make_gate_chi will generally have a slightly smaller numerical overhead, as it pushes the use of automatic differentiation down by one level.

With use_finite_differences=true, this routine serves to test and debug gradients for gate functionals obtained by automatic differentiation.

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

grad_J_a! = make_grad_J_a(
    J_a,
    tlist;
    force_zygote=false,
    use_finite_differences=false
)

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.

If force_zygote=true, automatic differentiation with Zygote will be used to calculate the derivative.

If use_finite_differences=true, the derivative will be calculated via finite differences. This may be used to verify Zygote gradients.

By default, for functionals J_a that have a known analytic derivative, that analytic derivative will be used. For unknown functions, the derivative will be calculated via Zygote.

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

State-to-state phase-insensitive fidelity.

F_ss(ϕ, objectives; τ=nothing)

calculates

\[F_{\text{ss}} = \frac{1}{N} \sum_{k=1}^{N} w_k |τ_k|^2 \quad\in [0, 1]\]

with $N$, $w_k$ and $τ_k$ as in f_tau.

Backward boundary states $|χ⟩$ for functional J_T_re.

chi_re!(χ, ϕ, objectives; τ=nothing)

sets the elements of χ according to

\[|χ_k⟩ = -\frac{∂ J_{T,\text{re}}}{∂ ⟨ϕ_k(T)|} = \frac{1}{2N} w_k |ϕ^{\tgt}_k⟩\]

with $|ϕ^{\tgt}_k⟩$ and $w_k$ as defined in f_tau.

Note: this function can be obtained with make_chi(J_T_re, objectives).

Analytic derivative for J_a_fluence.

grad_J_a_fluence!(∇J_a, pulsevals, tlist)

sets the (vectorized) elements of ∇J_a to $2 ϵ_{ln} dt$, where $ϵ_{ln}$ are the (vectorized) elements of pulsevals and $dt$ is the time step, taken from the first time interval of tlist and assumed to be uniform.

Real-part fidelity.

F_re(ϕ, objectives; τ=nothing)

calculates

\[F_{\text{re}} = \Re[f_{τ}] = \Re\left[ \frac{1}{N} \sum_{k=1}^{N} w_k τ_k \right] \quad\in \begin{cases} [-1, 1] & \text{in Hilbert space} \\ [0, 1] & \text{in Liouville space.} \end{cases}\]

with $w_k$ the weight for the k'th objective and $τ_k$ the overlap of the k'th propagated state with the k'th target state, and $N$ the number of objectives.

All arguments are passed to f_tau to evaluate $f_τ$.

Average complex overlap of the target states with forward-propagated states.

f_tau(ϕ, objectives; τ=nothing)

calculates

\[f_τ = \frac{1}{N} \sum_{k=1}^{N} w_k τ_k\]

with

\[τ_k = ⟨ϕ_k^\tgt|ϕ_k(T)⟩\]

in Hilbert space, or

\[τ_k = \tr[ρ̂_k^{\tgt\,\dagger} ρ̂_k(T)]\]

in Liouville space, where $|ϕ_k⟩$ or $ρ̂_k$ are the elements of ϕ, and $|ϕ_k^\tgt⟩$ or $ρ̂_k^\tgt$ are the target states from the target_state field of the objectives. If τ is given as a keyword argument, it must contain the values τ_k according to the above definition. Otherwise, the $τ_k$ values will be calculated internally.

$N$ is the number of objectives, and $w_k$ is the weight attribute for each objective. The weights are not automatically normalized, they are assumed to have values such that the resulting $f_τ$ lies in the unit circle of the complex plane. Usually, this means that the weights should sum to $N$.

Backward boundary states $|χ⟩$ for functional J_T_ss.

chi_ss!(χ, ϕ, objectives; τ=nothing)

sets the elements of χ according to

\[|χ_k⟩ = -\frac{∂ J_{T,\text{ss}}}{∂ ⟨ϕ_k(T)|} = \frac{1}{N} w_k τ_k |ϕ^{\tgt}_k⟩\,,\]

with $|ϕ^{\tgt}_k⟩$, $τ_k$ and $w_k$ as defined in f_tau.

Note: this function can be obtained with make_chi(J_T_ss, objectives).

Square-modulus fidelity.

F_sm(ϕ, objectives; τ=nothing)

calculates

\[F_{\text{sm}} = |f_τ|^2 = \left\vert\frac{1}{N} \sum_{k=1}^{N} w_k τ_k\right\vert^2 = \frac{1}{N^2} \sum_{k=1}^{N} \sum_{j=1}^{N} w_k w_j τ̄_k τ_j \quad\in [0, 1]\,,\]

with $w_k$ the weight for the k'th objective and $τ_k$ the overlap of the k'th propagated state with the k'th target state, $τ̄_k$ the complex conjugate of $τ_k$, and $N$ the number of objectives.

All arguments are passed to f_tau to evaluate $f_τ$.

Backward boundary states $|χ⟩$ for functional J_T_sm.

chi_sm!(χ, ϕ, objectives; τ=nothing)

sets the elements of χ according to

\[|χ_k⟩ = -\frac{\partial J_{T,\text{sm}}}{\partial ⟨ϕ_k(T)|} = \frac{1}{N^2} w_k \sum_{j}^{N} w_j τ_j |ϕ_k^{\tgt}⟩\]

with $|ϕ^{\tgt}_k⟩$, $τ_j$ and $w_k$ as defined in f_tau.

Note: this function can be obtained with make_chi(J_T_sm, objectives).

Blackman window shape.

blackman(t, t₀, T; a=0.16)

calculates

\[B(t; t_0, T) = \frac{1}{2}\left( 1 - a - \cos\left(2π \frac{t - t_0}{T - t_0}\right) + a \cos\left(4π \frac{t - t_0}{T - t_0}\right) \right)\,,\]

for a scalar t, with $a$ = 0.16.

See http://en.wikipedia.org/wiki/Window_function#Blackman_windows

A Blackman shape looks nearly identical to a Gaussian with a 6-sigma interval between t₀ and T. Unlike the Gaussian, however, it will go exactly to zero at the edges. Thus, Blackman pulses are often preferable to Gaussians.

Box shape (Theta-function).

box(t, t₀, T)

evaluates the Heaviside (Theta-) function $\Theta(t) = 1$ for $t_0 \le t \le T$; and $\Theta(t) = 0$ otherwise.

Flat shape (amplitude 1.0) with a switch-on/switch-off from zero.

flattop(t; T, t_rise, t₀=0.0, t_fall=t_rise, func=:blackman)

evaluates a shape function that starts at 0 at $t=t₀$, and ramps to to 1 during the t_rise interval. The function then remains at value 1, before ramping down to 0 again during the interval t_fall before T. For $t < t₀$ and $t > T$, the shape is zero.

The default switch-on/-off shape is half of a Blackman window (see blackman).

For func=:sinsq, the switch-on/-off shape is a sine-squared curve.

Set up a dummy control problem.

problem = dummy_control_problem(;
    N=10, n_objectives=1, n_controls=1, n_steps=50, dt=1.0, sparsity=0.5,
    complex_operators=true, hermitian=true, kwargs...)

Sets up a control problem with random (sparse) Hermitian matrices.

Arguments

• N: The dimension of the Hilbert space
• n_objectives: The number of objectives in the optimization. All objectives will have the same Hamiltonian, but random initial and target states.
• n_controls: The number of controls, that is, the number of control terms in the control Hamiltonian. Each control is an array of random values, normalized on the intervals of the time grid.
• n_steps: The number of time steps (intervals of the time grid)
• dt: The time step
• sparsity: The sparsity of the Hamiltonians, as a number between 0.0 and 1.0. For sparsity=1.0, the Hamiltonians will be dense matrices.
• complex_operators: Whether or not the drift/control operators will be complex-valued or real-valued.
• hermitian: Whether or not all drift/control operators will be Hermitian matrices.
• kwargs: All other keyword arguments are passed on to ControlProblem

Generate an HTML report for existing coverage data.

generate_coverage_html(path="./"; covdir="coverage", genhtml="genhtml")

creates a folder covdir and use the external genhtml program to write an HTML coverage report into that folder.

Construct a random complex matrix of size N×N with spectral radius ρ.

random_complex_matrix(N, ρ)

Construct a random sparse complex matrix.

random_complex_sparse_matrix(N, ρ, sparsity)

returns a matrix of size N×N with spectral radius ρ and the given sparsity (number between zero and one that is the approximate fraction of non-zero elements).

Construct a random Hermitian matrix of size N×N with spectral radius ρ.

random_hermitian_matrix(N, ρ)

Construct a random Hermitian real matrix of size N×N with spectral radius ρ.

random_hermitian_real_matrix(N, ρ)

Construct a random sparse Hermitian matrix.

random_hermitian_sparse_matrix(N, ρ, sparsity)

returns a matrix of size N×N with spectral radius ρ and the given sparsity (number between zero and one that is the approximate fraction of non-zero elements).

Construct a random sparse Hermitian real matrix.

random_hermitian_sparse_real_matrix(N, ρ, sparsity)

returns a matrix of size N×N with spectral radius ρ and the given sparsity (number between zero and one that is the approximate fraction of non-zero elements).

Construct a random real-valued matrix of size N×N with spectral radius ρ.

random_real_matrix(N, ρ)

Construct a random sparse real-valued matrix.

random_real_sparse_matrix(N, ρ, sparsity)

returns a matrix of size N×N with spectral radius ρ and the given sparsity (number between zero and one that is the approximate fraction of non-zero elements).

Return a random, normalized Hilbert space state vector of dimension N.

random_state_vector(N)

Print out a coverage summary from existing coverage data.

show_coverage(path="./src"; sort_by=nothing)

prints a a table showing the tracked files in path, the total number of tracked lines in that file ("Total"), the number of lines with coverage ("Hit"), the number of lines without coverage ("Missed") and the "Coverage" as a percentage.

The coverage data is collected from .cov files in path.

Optionally, the table can be sorted by passing the name of a column to sort_by, e..g. sort_py=:Missed.

Run a package test-suite in a subprocess.

test(
    file="test/runtests.jl";
    root=pwd(),
    project="test",
    code_coverage="user",
    show_coverage=(code_coverage == "user"),
    color=<inherit>,
    compiled_modules=<inherit>,
    startup_file=<inherit>,
    depwarn=<inherit>,
    inline=<inherit>,
    check_bounds="yes",
    track_allocation=<inherit>,
    threads=<inherit>,
    genhtml=false,
    covdir="coverage"
)

runs the test suite of the package located at root by running include(file) inside a new julia process.

This is similar to what Pkg.test() does, but differs in the "sandboxing" approach. While Pkg.test() creates a new temporary sandboxed environment, test() uses an existing environment in project (the test subfolder by default). This allows testing against the dev-versions of other packages. It requires that the test folder contains both a Project.toml and a Manifest.toml file.

The test() function also differs from directly including test/runtests.jl in the REPL in that it can generate coverage data and reports (this is only possible when running tests in a subprocess).

If show_coverage is passed as true (default), a coverage summary is shown. Further, if genhtml is true, a full HTML coverage report will be generated in covdir (relative to root). This requires the genhtml executable (part of the lcov package). Instead of true, it is also possible to pass the path to the genhtml exectuable.

All other keyword arguments correspond to the respective command line flag for the julia executable that is run as the subprocess.

This function is intended to be exposed in a project's development-REPL.

Result returned by optimize_with_dummy_method.

Run a dummy optimization.

result = optimize(problem, method=:dummymethod)

runs through and "optimization" of the given problem where in each iteration, the amplitude of the guess pulses is diminished by 10%. The (summed) vector norm of the the control serves as the value of the optimization functional.

Perfect-entanglers distance measure.

D = D_PE(U; unitarity_weight=0.0, absolute_square=false)

For a given two-qubit gate $Û$, this is defined via the local_invariants $g_1$, $g_2$, $g_3$ as

\[D = g_3 \sqrt{g_1^2 + g_2^2} - g_1\]

This describes the geometric distance of the quantum gate from the polyhedron of perfect entanglers in the Weyl chamber.

This equation is only meaningful under the assumption that $Û$ is unitary. If the two-qubit level are a logical subspace embedded in a larger physical Hilbert space, loss of population from the logical subspace may lead to a non-unitary $Û$. In this case, the unitarity measure can be added to the functional by giving a unitary_weight ∈ [0, 1) that specifies the relative proportion of the $D$ term and the unitarity term.

By specifying absolute_square=true, the functional is modified as $D → |D|²$, optimizing specifically for the boundary of the perfect entanglers polyhedron. This accounts for the fact that $D$ can take negative values inside the polyhedron, as well as the W1 region of the Weyl chamber (the one adjacent to SWAP). This may be especially useful in a system with population loss (unitarity_weight > 0), as it avoids situations where the optimization goes deeper into the perfect entanglers while increasing population loss.

Construct the canonical gate for the given Weyl chamber coordinates.

Û = canonical_gate(c₁, c₂, c₃)

constructs the two qubit gate $Û$ as

\[Û = \exp\left[i\frac{π}{2} (c_1 σ̂_x σ̂_x + c_2 σ̂_y σ̂_y + c_3 σ̂_z σ̂_z)\right]\]

where $σ̂_{x,y,z}$ are the Pauli matrices.

Calculate the maximum gate concurrence.

C = gate_concurrence(U)
C = gate_concurrence(c₁, c₂, c₃)

calculates that maximum concurrence $C ∈ [0, 1]$ that the two two-qubit gate U, respectively the gate described by the Weyl chamber coordinates c₁, c₂, c₃ (see weyl_chamber_coordinates) can generate.

See Kraus, Cirac, Phys. Rev. A 63, 062309 (2001)

Check whether a given gate is in (a specific region of) the Weyl chamber.

in_weyl_chamber(c₁, c₂, c₃)

checks whether c₁, c₂, c₃ are valid Weyl chamber coordinates.

in_weyl_chamber(U; region="PE")
in_weyl_chamber(c₁, c₂, c₃; region="PE")

checks whether the two-qubit gate U, respectively the gate described by the Weyl chamber coordinates c₁, c₂, c₃ (see weyl_chamber_coordinates) is a perfect entangler. The region can be any other of the regions returned by weyl_chamber_region.

Calculate the local invariants g₁, g₂, g₃ for a two-qubit gate.

g₁, g₂, g₃ = local_invariants(U)

Unitarity of a matrix.

pop_loss = 1 - unitarity(U)

measures the loss of population from the subspace described by U. E.g., for a two-qubit gate, U is a 4×4 matrix. The unitarity is defined as $\Re[\tr(Û^†Û) / N]$ where $N$ is the dimension of $Û$.

Calculate the Weyl chamber coordinates c₁, c₂, c₃ for a two-qubit gate.

c₁, c₂, c₃ = weyl_chamber_coordinates(U)

calculates the Weyl chamber coordinates using the algorithm described in Childs et al., PRA 68, 052311 (2003).

Identify which region of the Weyl chamber a given gate is located in.

region = weyl_chamber_region(U)
region = weyl_chamber_region(c₁, c₂, c₃)

identifies which region of the Weyl chamber the given two-qubit gate U, respectively the gate identified by the Weyl chamber coordinates c₁, c₂, c₃ (see weyl_chamber_coordinates) is in, as a string. Possible outputs are:

• "PE": gate is in the polyhedron of perfect entanglers.
• "W0": gate is between the identity and the perfect entanglers.
• "W0*": gate is between CPHASE(2π) and the perfect entanglers.
• "W1": gate is between SWAP and the perfect entanglers.

For invalid Weyl chamber coordinates, an empty string is returned.