Core Interface

Integrations

Everything you need to know how to integrate over time.

IntegrationType{F}

Encapsulates a time-grid, used to decide how to integrate over time.

Type Parameters

• F denotes the type for time values. Must be a real float.

Implementation

Any concrete sub-type G must implement the Prototypes interface.

In addition, the following methods must be implemented.

• nsteps(::G): the total number of finite jumps within this time grid.

  This is the maximum index of the timeat and stepat functions, but it is NOT the length of the lattice, since the minimum index is 0.

• timeat(::G, i::Int): time at index i (which may be zero)

• stepat(::G, i::Int): stepsize at index i (which may be zero)

  Let $r$ be the number of steps, $T$ be the duration, and $τ_i$ the step at index $i$. Then $∑_{i=0}^r τ_i = T$.

• Prototype(::Type{G}, r::Int; T, kwargs...): construct a prototypical grid of type G with r steps.

AbstractVector Interface

This type implements the AbstractVector interface, defined so that grid[i] == timeat(grid, i), where grid is the IntegrationType. Note that i here starts from 0. The collect function produces a concrete vector, in which i starts from 1. The lattice function does the same thing, but permits an allocation-free signature.

duration(grid::IntegrationType)

The total duration, ie. the last time point minus the first.

endtime(grid::IntegrationType)

Upper bound of a time integral.

integrate(grid::IntegrationType{F}, Φ)

Compute the integral $∫_0^T Φ(t)⋅dt$. Φ is a univariate scalar function of time, returning type F.

integrate(grid::IntegrationType, f̄::AbstractVector)

Treat the elements of f̄ as the function evaluations Φ(t) above. The length of f̄ must be length(grid) and its eltype must be real.

integrate(grid::IntegrationType, Φ, f̄s::AbstractVector...)

Compute the integral $∫_0^T Φ(t, f1, f2...)⋅dt$, where f1 is the value of f̄s[1] at the index corresponding to time t, etc. Φ is a multivariate scalar function, of time and one argument for each f̄, returning type F. The length of each element in f̄s must be length(grid).

lattice(grid::IntegrationType; result=nothing)

A vector of all time points.

This is equivalent to collect(grid) except for the result kwarg, which allows the caller to provide a pre-allocated array.

nsteps(grid::IntegrationType)::Int

The total number of finite jumps within this time grid.

This is the maximum index of the timeat and stepat functions, but it is NOT the length of the lattice, since the minimum index is 0.

starttime(grid::IntegrationType)

Lower bound of a time integral.

stepat(grid::IntegrationType{F}, i::Int)::F

The stepsize at index i (which may be zero).

Let $r$ be the number of steps, $T$ be the duration, and $τ_i$ the step at index $i$. Then $∑_{i=0}^r τ_i = T$.

stepsize(grid::IntegrationType)

The average step size, ie. duration divided by number of steps.

timeat(grid::IntegrationType{F}, i::Int)::F

The time at index i (which may be zero).

Signals

Time-dependent functions suitable for control signals with variational parameters.

The main motivation of this module is to provide a common interface for analytical gradients and optimization.

ConstrainedSignal(template::<:ParametricSignal, constraints::Vector{Symbol})

The parametric signal template, freezing all fields specified by constraints.

Frozen parameters are omitted from the Parameters interface. In other words, they do not appear in Parameters.names or Parameters.values, and they are not mutated by Parameters.bind!.

Example

Say you have a Trigonometric signal sub-typing ParametricSignal, where all fields amplitude A, phase ϕ, and frequency ν are registered as variational parameters. You want to run an optimization where the frequency ν is fixed to 4.608. Rather than implementing a whole new ParametricSignal identical except for a different implementation of the parameters function, use a ConstrainedSignal:

template = Trigonometric(0.0, 0.0, 4.608)   # Initialize A and ϕ to 0.0
signal = Constrained(template, :ν)

Note that the Constrained constructor above is just syntactic sugar for:

signal = ConstrainedSignal(template, [:ν])

ParametricSignal{P,R}

The super-type for simple user-defined SignalTypes.

Type Parameters

• P denotes the type of all variational parameters. Must be a real float.
• R denotes the type of $Ω(t)$ itself. May be any number type.

Implementation

Concrete sub-types S must be mutable structs, and all of its variational parameters (as indicated by the parameters function) must have type P.

The following methods must be implemented:

• Signals.parameters(Ω::S): returns a tuple of the fields in S treated as variational parameters.

• Signals.valueat(Ω::S, t::Real): the actual function $Ω(t)$. Must return a number of type R.

• Signals.partial(k::Int, Ω::S, t::Real): the partial derivative $∂Ω/∂x_k$ evaluated at time t, where $x_k$ is Ω's k-th variational parameter (ie. Parameters.names(Ω)[k]). Must return a number of type R.

• Base.string(Ω::S, names::AbstractVector{String}): a human-readable description of the signal, inserting each element of names in the place of the corresponding parameter. For example, complex constant signals may return a description like "A + i B", where A and B are the "names" given by the names argument.

SignalType{P,R}

Encapsulates a parametric and differntiable scalar function $Ω(t)$.

Type Parameters

• P denotes the type of all variational parameters. Must be a real float.
• R denotes the type of $Ω(t)$ itself. May be any number type.

Implementation

Any concrete sub-type S must implement all the Parameters interface.

• In particular, Parameters.values(Ω::S) must return a vector of type P.
• If you are trying to create your own signal type, you probably want to implement a ParametricSignal, which already has an implementation for the Parameters interface.

In addition, the following methods must be implemented:

• valueat(Ω::S, t::Real): the actual function $Ω(t)$. Must return a number of type R.

• partial(i::Int, Ω::S, t::Real): the partial derivative $∂Ω/∂x_i$ evaluated at time t, where $x_i$ is Ω's i-th variational parameter (ie. Parameters.names(Ω)[i]). Must return a number of type R.

• Base.string(Ω::S, names::AbstractVector{String}): a human-readable description of the signal, inserting each element of names in the place of the corresponding parameter. For example, a complex constant signal might return a description like "A + i B", where A and B are the "names" given by the names argument.

(signal::SignalType{P,R})(t)

Syntactic sugar: if Ω is a SignalType, then Ω(t) gives valueat(Ω,t).

The time t may be a scalar time, (abstract) vector of times, or an IntegrationType.

Base.string(Ω::SignalType, names::AbstractVector{String})::String

Substitutes the default name of each variational parameter for the ones in names.

Note that this is not the usual signature for Base.string!

Base.string(Ω::SignalType)

A human-readable string description of the signal.

Constrained(template::ParametricSignal, constraints::Symbol...)

Construct a ConstrainedSignal from a ParametricSignal and the fields to freeze.

parameters(::S<:ParametricSignal{P,R})::Tuple{Symbol...}

Returns a tuple of the fields in S treated as variational parameters.

parametertype(signal)

Returns the number type for parameters in signal.

partial(k::Int, signal::SignalType{P,R}, t::Real)::R

The partial derivative $∂Ω/∂x_k|_t$.

Here $x_k$ is the signal's k-th variational parameter (ie. Parameters.names(signal)[k]).

partial(i, signal, t̄; result=nothing)

Evaluate the partial at each point in the time lattice defined by an integration.

Parameters

• i: indexes which parameter to take the partial derivative with respect to
• signal: the SignalType.
• grid: an IntegrationType.

Keyword Arguments

• result: if provided, a pre-allocated array to store the returned vector

returntype(signal)

Returns the number type for function values of signal.

valueat(signal::SignalType{P,R}, t::Real)::R

The signal at time t, ie. $Ω(t)$.

valueat(signal, grid; result=nothing)

Evaluate the signal at each point in the time lattice defined by an integration.

Parameters

• signal: the SignalType.
• grid: an IntegrationType.

Keyword Arguments

• result: if provided, a pre-allocated array to store the returned vector

Devices

In silico representation of quantum devices, in which quantum states evolve in time.

In this package, the "static" components (e.g., qubit frequencies, couplings, etc.) and the "drive" components (e.g., control signal, variational parameters, etc.) are all integrated into a single DeviceType object. All you need to know how a quantum state ψ evolves up time T is in the device.

DeviceType{F}

Encapuslates a device Hamiltonian, under which quantum computational states evolve.

Type Parameters

• F: the float type associated with a device

Implementation

Any concrete sub-type D must implement the Prototypes and Parameters interfaces.

• In particular, if any static operators in your device depend on variational parameters, you should consult the "Note on Caching" below.

In addition, all methods in the following sections must be implemented.

• Counting methods
• Algebra methods
• Operator methods
• Gradient methods
• Benchmarking methods

If your device's drive channels are all local (acting on one qubit at a time), you should implement a LocallyDrivenDevice, which has a few extra requirements.

Counting methods:

• nqubits(::D): the number of qubits in the device - call this n.
• nlevels(::D): the number of physical levels in each "qubit" - call this m.
• ndrives(::D): the number of distinct drive channels.
• ngrades(::D): the number of distinct gradient operators.
• noperators(::D): the number of operators used to define an algebra for each "qubit".

Each of these methods returns an integer.

Algebra methods:

• localalgebra(::D): a 4d array ā where ā[:,:,σ,q] is the σ'th algebraic operator on qubit q, in the bare basis.

This method should define result=nothing as a keyword argument; when passed, use it as the array to store your result in.

Operator methods:

• qubithamiltonian(::D, ā, q::Int): the static components of the device Hamiltonian local to qubit q.

• staticcoupling(::D, ā): the static components of the device Hamiltonian nonlocal to any one qubit.

• driveoperator(::D, ā, i::Int, t::Real): the distinct drive operator for channel i at time t

• gradeoperator(::D, ā, j::Int, t::Real): the distinct gradient operator indexed by j at time t

Each of these methods should define result=nothing as a keyword argument; when passed, use it as the array to store your result in. When result is not passed, you should return a new array of type Complex{eltype(D)}. Aside from this, do your best to minimize allocations.

Each of these methods takes a 4darray ā, indexed as described in localalgebra. For example, for an algebra defined in terms of bosonic ladder operators, $\hat a_q$ might be accessed with ā[:,:,1,q] and $\hat a_q^\dagger$ with ā[:,:,1,q]'. For an algebra defined in terms of Pauli operators, $X_q$ might be accessed with ā[:,:,1,q], $Y_q$ with ā[:,:,2,q], and $Z_q$ with ā[:,:,3,q].

Usually, each ā[:,:,σ,q] is defined on the full Hilbert space (ie. m^n × m^n), but sometimes the code exploits a simple tensor structure by passing in local m × m operators instead, so do not assume a specific size a priori. Do NOT modify these operators, as they are usually drawn from a cache.

Gradient methods:

• gradient(::D, grid::Integrations.IntegrationType, ϕ̄): the gradient vector for each variational parameter in the device.

Each partial is generally an integral over at least one gradient signal. The argument grid identifies the temporal lattice on which ϕ̄ is defined. The argument ϕ̄ is a 2d array; ϕ̄[:,:,j] contains the jth gradient signal $ϕ_j(t)$ evaluated at each point in grid.

This method should define result=nothing as a keyword argument; when passed, use it as the array to store your result in.

Notes on Caching

This module uses the Memoization package to cache some arrays as they are calculated.

This does not apply to any method which depends on an absolute time t, though it does apply to methods depending only on a relative time τ. For example, the propagator for a static Hamiltonian is cached, but not one for a drive Hamiltonian.

Usually, variational parameters only affect time-dependent methods, but if any of your device's static operators do depend on a variational parameter, you should be careful to empty the cache when Parameters.bind! is called.

You can completely clear everything in the cache with:

Memoization.empty_all_caches!()

Alternatively, selectively clear caches for affected functions via:

Memoization.empty_cache!(fn)

I don't know if it's possible to selectively clear cached values for specific methods. If it can be done, it would require obtaining the actual Dict being used as a cache for a particular function, figuring out exactly how that cache is indexed, and manually removing elements matching your targeted method signature.

LocallyDrivenDevice

Super-type for device objects whose drive channels act locally on individual qubits.

Inherit from this type if your driveoperator and gradeoperator methods depend only on a single qubit, i.e. ā[:,:,:,q]. This enables more efficient propagation methods which exploit a tensor product structure.

Implementation

Any concrete sub-type D must implement everything required in the DeviceType interface, so consult the documentation for DeviceType carefully.

In addition, the following methods must be implemented:

• drivequbit(::D, i::Int): index of the qubit on which channel i is applied.
• gradequbit(::D, j::Int): index of the qubit associated with the jth gradient operator.

It's usually trivial to infer the channel index i associated with each gradient operator, in which case gradequbit(device, j) = drivequbit(device, i), but this is left as an implementation detail.

basisrotation(tgt::Bases.BasisType, src::Bases.BasisType, device::DeviceType)

Calculate the basis rotation U which transforms $|ψ_{src}⟩ → |ψ_{tgt}⟩ = U|ψ_{src}⟩$.

braket(op, device[, basis], ψ1, ψ2)

The braket of an operator describing a device with respect to states ψ1 and ψ2.

If $A$ is the operator specified by op, this method calculates $⟨ψ1|A|ψ2⟩$.

Arguments

• op::Operators.OperatorType: which operator to estimate (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the state ψ is represented in. Defaults to Bases.BARE when omitted.
• ψ1, ψ2: Statevectors defined over the full Hilbert space of the device.

dress(device::DeviceType)

Diagonalize the static Hamiltonian and apply some post-processing to define the so-called DRESSED basis.

Compute the vector of eigenvalues Λ and the rotation matrix U for a given basis.

U is an operator acting on the global Hilbert space of the device.

The result is packed into a LinearAlgebra.Eigen object. It may be unpacked directly into a vector eigenvalues Λ and a matrix of eigenvectors U by

Λ, U = dress(device)

Alternatively:

ΛU = dress(device)
Λ = ΛU.values
U = ΛU.vectors

driveoperator(device::DeviceType, ā, i::Int, t::Real; result=nothing)

The distinct drive operator for channel i at time t.

This method is a function of algebraic operators given by ā[:,:,σ,q], constructed the globalalgebra method. If device is a LocallyDrivenDevice, ā may also have been constructed from the localalgebra method.

The array is stored in result if provided. If result is not provided, the array is of type Complex{eltype(device)}.

drivequbit(device, i::Int)::Int

Index of the qubit on which channel i is applied.

evolve!(op, device[, basis], t, ψ)

Propagate a state ψ by a time t under the Hermitian op describing a device.

This function is identical to propagate!, except that the cache is not used for intermediate propagator matrices, and that it is undefined for time-dependent operators. Look to the Evolutions module for algorithms compatible with time-dependence!

Arguments

• op::Operators.OperatorType: which operator to evolve under (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the state ψ is represented in. Defaults to Bases.BARE when omitted.
• t::Real: the amount to move forward in time by.
• ψ: Either a vector or a matrix, defined over the full Hilbert space of the device.

evolver(op, device[, basis], t; kwargs...)

A unitary propagator describing evolution under a Hermitian operator for a time t.

This function is identical to propagator, except that the argument t is considered an absolute time so it is never cached, and that it is undefined for time-dependent operators. It exists solely to perform rotating-frame rotations at every time-step without worrying about over-caching.

Arguments

• op::Operators.OperatorType: which operator to evolve under (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the operator will be represented in. Defaults to Bases.BARE when omitted.
• t::Real: the amount to move forward in time by.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result returns an array of type Complex{eltype(device)}.

expectation(op, device[, basis], ψ)

The expectation value of an operator describing a device with respect to the state ψ.

If $A$ is the operator specified by op, this method calculates $⟨ψ|A|ψ⟩$.

Arguments

• op::Operators.OperatorType: which operator to estimate (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the state ψ is represented in. Defaults to Bases.BARE when omitted.
• ψ: A statevector defined over the full Hilbert space of the device.

globalalgebra(device::DeviceType[, basis::Bases.BasisType])

A globalized version of the matrix list defined by localalgebra.

That is, each operator ā[:,:,q,σ] is an N⨯N matrix acting on the whole Hilbert space, rather than an m⨯m matrix acting on the space of a single qubit.

(N is used throughout the docs for nstates(device) and m for nlevels(device).)

The array is stored in result or, if not provided, returned from a cache.

gradeoperator(device::DeviceType, ā, j::Int, t::Real; result=nothing)

The distinct gradient operator indexed by j at time t.

I have defined the "gradient operator" $Â_j$ as the Hermitian operator for which the jth gradient signal is $ϕ_j = ⟨λ|(iÂ_j)|ψ⟩ + h.t.$.

This method is a function of algebraic operators given by ā[:,:,σ,q], constructed the globalalgebra method. If device is a LocallyDrivenDevice, ā may also have been constructed from the localalgebra method.

The array is stored in result if provided. If result is not provided, the array is of type Complex{eltype(device)}.

gradequbit(device, j::Int)::Int

Index of the qubit associated with the jth gradient operator.

gradient(::DeviceType, grid::Integrations.IntegrationType, ϕ; result=nothing)

The gradient vector of partials for each variational parameter in the device.

Each partial is generally an integral over at least one gradient signal. The argument grid identifies the temporal lattice on which ϕ is defined. The argument ϕ is a 2d array; ϕ[:,j] contains the jth gradient signal $ϕ_j(t)$ evaluated at each point in grid.

The array is stored in result if provided.

localalgebra(::DeviceType; result=nothing)

Construct a 4d array ā representing all operators defining an algebra for all qubits.

The matrix ā[:,:,σ,q] is a local operator acting on the space of a single qubit (meaning it is an m⨯m matrix, if m is the result of nlevels(device)).

The array is stored in result or, if not provided, returned from a cache. If result is not provided, the array is of type Complex{eltype(device)}.

Implementation

To use the cache, simply include as the first line:

isnothing(result) && return _localalgebra(device)

If there is nothing yet in the cache, the _localalgebra function will simply call your method again, but with an empty array passed as result.

localdriveoperators(device, t; kwargs...)

A matrix list v̄, where v̄[:,:,q] represents a sum of all drives acting on qubit q in the bare basis.

Arguments

• device::LocallyDrivenDevice: which device is being described.
• t::Real: the time each drive operator is evaluated at.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result will return an array with type Complex{eltype(device)}.

localdrivepropagators(device, τ, t; kwargs...)

A matrix list ū, where ū[:,:,q] is the propagator for a local drive term.

Arguments

• device::LocallyDrivenDevice: which device is being described.
• τ::Real: the amount to move forward in time by. Note that the propagation is only approximate for time-dependent operators. The smaller τ is, the more accurate the approximation.
• t::Real: the time each drive operator is evaluated at.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result will return an array with type Complex{eltype(device)}.

localqubitevolvers(device, τ; kwargs...)

A matrix list ū, where each ū[:,:,q] is a propagator for a local qubit hamiltonian.

This function is identical to localqubitevolvers, except that the argument t is considered an absolute time so it is never cached.

Arguments

• device::DeviceType: which device is being described.
• τ::Real: the amount to move forward in time by. Note that the propagation is only approximate for time-dependent operators. The smaller τ is, the more accurate the approximation.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result will return an array of type Complex{eltype(device)}.

localqubitoperators(device[, basis]; kwargs...)

A matrix list h̄, where each h̄[:,:,q] represents a local qubit hamiltonian in the bare basis.

Arguments

• device::DeviceType: which device is being described.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result will return a cached result, with type Complex{eltype(device)}.

localqubitpropagators(device, τ; kwargs...)

A matrix list ū, where each ū[:,:,q] is a propagator for a local qubit hamiltonian in the bare basis.

Arguments

• device::DeviceType: which device is being described.
• τ::Real: the amount to move forward in time by.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result will return a cached result, with type Complex{eltype(device)}.

ndrives(device::DeviceType)::Int

The number of distinct drive channels.

ngrades(device::DeviceType)::Int

The number of distinct gradient operators.

nlevels(device::DeviceType)::Int

The number of physical levels in each "qubit".

noperators(device::DeviceType)::Int

The number of operators used to define an algebra for each "qubit".

nqubits(device::DeviceType)::Int

The number of qubits in the device.

nstates(device::DeviceType)

The total number of states in the physical Hilbert space of the device.

(This is as opposed to nlevels(device), the number of states in the physical Hilbert space of a single independent qubit.)

operator(op, device[, basis]; kwargs...)

A Hermitian operator describing a device, represented in the given basis.

For example, to construct the static Hamiltonian of a device in the dressed basis, call operator(Operators.STATIC, device, Bases.DRESSED).

Arguments

• op::Operators.OperatorType: which operator to construct (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the operator will be represented in. Defaults to Bases.BARE when omitted.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result returns an array of type Complex{eltype(device)}. For static operators only, omitting result will return a cached result.

propagate!(op, device[, basis], τ, ψ)

Propagate a state ψ by a small time τ under the Hermitian op describing a device.

Arguments

• op::Operators.OperatorType: which operator to evolve under (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the state ψ is represented in. Defaults to Bases.BARE when omitted.
• τ::Real: the amount to move forward in time by. Note that the propagation is only approximate for time-dependent operators. The smaller τ is, the more accurate the approximation.
• ψ: Either a vector or a matrix, defined over the full Hilbert space of the device.

propagator(op, device[, basis], τ; kwargs...)

A unitary propagator describing evolution under a Hermitian operator for a small time τ.

Arguments

• op::Operators.OperatorType: which operator to evolve under (eg. static, drive, etc.).
• device::DeviceType: which device is being described.
• basis::Bases.BasisType: which basis the operator will be represented in. Defaults to Bases.BARE when omitted.
• τ::Real: the amount to move forward in time by. Note that the propagation is only approximate for time-dependent operators. The smaller τ is, the more accurate the approximation.

Keyword Arguments

• result: a pre-allocated array of compatible type and shape, used to store the result.

  Omitting result returns an array of type Complex{eltype(device)}. For static operators only, omitting result will return a cached result.

qubithamiltonian(device::DeviceType, ā, q::Int; result=nothing)

The static components of the device Hamiltonian local to qubit q.

This method is a function of algebraic operators given by ā[:,:,σ,q], constructed by either the localalgebra or globalalgebra methods.

The array is stored in result if provided. If result is not provided, the array is of type Complex{eltype(device)}.

staticcoupling(device::DeviceType, ā, q::Int; result=nothing)

The static components of the device Hamiltonian nonlocal to any one qubit.

This method is a function of algebraic operators given by ā[:,:,σ,q], constructed the globalalgebra method.

The array is stored in result if provided. If result is not provided, the array is of type Complex{eltype(device)}.

Evolutions

Algorithms to run time evolution, and related constructs like gradient signals.

EvolutionType

Defines a particular algorithm for performing time evolution.

Implementation

Any concrete sub-type A must implement the following methods:

• workbasis(::A): which Bases.BasisType the evolution algorithm uses

• evolve!(::A, device, grid, ψ; callback=nothing): evolve a state ψ in-place on a time grid

If possible, it should also implement:

• `gradientsignals(::A, device, basis, grid, ψ0, r, Ō; kwargs...): compute the gradient signals of a device corresponding to multiple observables

evolve(evolution, device, grid, ψ0; result=nothing, kwargs...)
evolve(evolution, device, basis, grid, ψ0; result=nothing, kwargs...)

Evolve a quantum computational state under a device Hamiltonian.

If basis is provided, ψ0 is taken to be represented in that basis. Otherwise, the workbasis of evolution is assumed.

This method simply copies ψ0 (to result if provided, or else to a new array), then calls the mutating function evolve! on the copy. Please see evolve! for detailed documentation.

evolve!(evolution, device, grid, ψ; callback=nothing)
evolve!(evolution, device, basis, grid, ψ; callback=nothing)

Evolve a quantum computational state in-place under a device Hamiltonian.

This method both mutates and returns ψ.

Arguments

• evolution::EvolutionType: which evolution algorithm to use.
• device::Devices.DeviceType: specifies which Hamiltonian to evolve under.
• basis::Bases.BasisType: which basis ψ is represented in. Assumed to be the workbasis of evolution when omitted.
• grid::Integrations.IntegrationType: defines the time integration bounds (eg. from 0 to $T$)
• ψ::AbstractVector: the initial statevector, defined on the full Hilbert space of the device.

Keyword Arguments

• callback: a function which is called at each iteration of the time evolution. The function is passed three arguments: - i: indexes the iteration - t: the current time point - ψ: the current statevector, in the work basis The function is called after having evolved ψ into |ψ(t)⟩.

Implementation

Only the signature omitting basis need be implemented, so you can assume ψ is represented in whatever basis you return from workbasis.

The signature including basis will automatically rotate ψ into the workbasis, call the method you implement, and then rotate back.

gradientsignals(evolution, device, grid, ψ, Ō; kwargs...)
gradientsignals(evolution, device, basis, grid, ψ, Ō; kwargs...)

The gradient signals associated with a given device Hamiltonian, and an observable O.

Gradient signals are used to calculate analytical derivatives of a control pulse.

Arguments

• evolution::EvolutionType: which evolution algorithm to use.
• device::Devices.DeviceType: specifies which Hamiltonian to evolve under. Also identifies each of the gradient operators used to calculate gradient signals.
• basis::Bases.BasisType: which basis ψ is represented in. Assumed to be the workbasis of evolution when omitted.
• grid::Integrations.IntegrationType: defines the time integration bounds (eg. from 0 to $T$)
• ψ::AbstractVector: the initial statevector, defined on the full Hilbert space of the device.
• Ō::Union{LAT.MatrixList,AbstractMatrix}: a list of Hermitian observables, represented as matrices, or a single such matrix. Gradients are calculated with respect to the expectation ⟨O⟩ at time $T$.

This method signature assumes ψ, Ō are represented in the workbasis of evolution.

Keyword Arguments

• result: an (optional) pre-allocated array to store gradient signals
• callback: a function called at each iteration of the gradient signal calculation. The function is passed three arguments: - i: indexes the iteration - t: the current time point - ψ: the current statevector, in the BARE basis The function is called after having evolved ψ into |ψ(t)⟩, but before calculating ϕ̄[i,:]. Evolution here runs backwards.

Returns

A 3d array ϕ̄, where each ϕ̄[:,j,k] is the gradient signal $ϕ_j(t)$ defined with respect to the observable $Ô_k$, or a 2d array when Ō is just a single matrix rather than a matrix list.

Explanation

A gradient signal $ϕ_j(t)$ is defined with respect to a gradient operator $Â_j$, an observable $Ô$, a time-dependent state |ψ(t)⟩, and total pulse duration T.

Let us define the expectation value $E(T) ≡ ⟨ψ(T)|Ô|ψ(T)⟩$.

Define the co-state $|λ(t)⟩$ as the (un-normalized) statevector which satisfies $E(T)=⟨λ(t)|ψ(t)⟩$ for any time t∊[0,T]. The gradient signal is defined as $ϕ_j(t) ≡ ⟨λ(t)|(iÂ_j)|ψ(t)⟩ + h.t.$.

Implementation

Only the signature omitting basis need be implemented, so can assume ψ, Ō are represented in whatever basis you return from workbasis.

The signature including basis will automatically rotate ψ, Ō into the workbasis before calling the method you implement.

Similiarly, you should only implement the method where Ō is a MatrixList. When a single matrix, the (already-implemented) method will reshaped it into a MatrixList, call the method you implement, and then reshape the resulting array.

workbasis(evolution::EvolutionType)::Bases.BasisType

Which basis the evolution algorithm works in.

Also defines the default basis to interpret ψ as, in evolution methods.

CostFunctions

Detailed instructions for how to compute an energy from a set of variational parameters.

Note that the CostFunctionType defined here refers to a type, not a function. But you can call CostFunctions.cost_function(::CostFunctionType) to get a function. This is the thing you would feed into an optimizer.

CostFunctionType{F}

Encapsulates a cost function, the thing you plug into an optimization algorithm.

Implementation

Any concrete sub-type CF must implement the following methods:

• cost_function(::CF): returns a Function which takes a parameter vector and returns the output of the cost function ie. a callable expression (x::Vector) -> f(x)
• grad!function(::CF): returns a mutating Function which takes a gradient vector (to be mutated) and a parameter vector, and writes the gradient vector to the first argument. As a matter of habit, the resulting gradient vector should also be returned.
• Base.length(::CF): the number of parameters this cost function takes

If your cost function involves calculating the expectation value of a time-evolved state, you should implement an EnergyFunction (even if it isn't strictly an energy). This type has a couple extra requirements to allow energy trajectories over the evolution.

(fn::CostFunctionType)(x̄::AbstractVector)

Evaluate the value of fn at the point x̄.

This is syntactic sugar for constructing a dedicated cost_function and calling it. It will not normally take advantage of cached work variables, so avoid using it in high-performance code.

EnergyFunction{F}

Super-type for cost functions calculating the expectation value of a time-evolved state.

Implementation

Any concrete sub-type CF must implement everything required in the CostFunctionType interface, so consult the documentation for CostFunctionType carefully.

In additon, the following method must be implemented:

• trajectory_callback(::CF, E::AbstractVector; callback=nothing) returns a Function compatible with Evolutions.evolve callback ie. a callable expression (i::Int, t::Real, ψ::Vector) -> Nothing which sets E[i] to the energy of a partially evolved wavefunction ψ. If callback is provided, the function calls that callback afterwards.

Finally, the following methods must now accept a keyword argument:

• cost_function(::CF; callback=nothing): When callback is provided, the time evolution must call it at each timestep.

• grad!function(::CF; ϕ=nothing): When ϕ is provided, write the gradient signals to it.

Base.eltype(fn::CostFunctionType)

Gives the number type for parameters of this cost function.

Base.length(fn::CostFunctionType)::Int

Gives the number of parameters for this cost function.

cost_function(fn::CostFunctionType)::Function

Converts the struct fn into a literal function of a parameter vector.

The function accepts a parameter vector x̄ which should have the type and length given by eltype(fn) and length(fn). The function returns the value of fn at the point x̄.

cost_function(fn::EnergyFunction; callback=nothing)

Same as for CostFunctionType except that whenever the function is called, the time evolution calls callback (if provided) in each time step.

grad!function(fn::CostFunctionType)::Function

Constructs a (mutating) function to calculate the gradient of fn at a particular point.

The function accepts a gradient vector ∇f̄ (to be mutated) and a parameter vector x̄. Both should have the type and length given by eltype(fn) and length(fn). After the function is called, ∇f̄ contains the gradient of fn at the point x̄. As per the Julia guidelines on mutating functions, ∇f̄ itself should then be returned.

grad!function(fn::EnergyFunction; ϕ=nothing)

Same as for CostFunctionType except that whenever the function is called, ϕ (if provided) is updated to contain the gradient signals. The array ϕ should be a 3d array with shape (:,nG,nK), where nK is the number of observables in the energy function, nG is the number of gradient operators in the underlying device, and the remaining dimension is the size of the time grid.

grad_function(fn::CostFunctionType)

Constructs a Function to calculate the gradient of fn at a particular point.

The function accepts a parameter vector x̄ which should have the type and length given by eltype(fn) and length(fn). The function returns the gradient (a vector) of fn at the point x̄.

nobservables(fn::EnergyFunction)::Int

Identify the number of Hermitian observables needed for this energy function.

For example, to measure the normalized energy, separate observables are needed for both energy and normalization, and the results are combined in a non-linear way to produce the final outcome.

trajectory_callback(fn::EnergyFunction, E::AbstractVector; callback=nothing)::Function

Make a callback to write the energy at each step of a time evolution.

The callback function should be compatible with Evolutions.evolve (ie. a callable expression (i::Int, t::Real, ψ::Vector) -> Nothing), which sets E[1+i] to the energy of a partially evolved wavefunction ψ. Note the 1+, due to i indexing a time integration, which starts from 0.

If callback is provided, the function calls that callback afterwards (i.e. callback chaining).