Basic Implementations

TrapezoidalIntegration(t0::F, tf::F, r::Int)

First-order grid using a trapezoidal rule time integration.

Arguments

• t0, tf: lower and upper bounds of the integral
• r: the number of time steps (remember number of lattice points is r+1)

Explanation

This grid gives a first-order quadrature, meaning every step takes you forward in time.

Additionally, this grid uses uniform spacing. You'd think that would mean stepat(grid, i) would give $τ_i = τ$ for every i, where $τ ≡ T / r$ and $T ≡ t_f - t_0$. But careful! The sum of all τ_i must match the length of the integral, ie. T. But there are r+1 points, and (r+1)⋅τ > T. How do we reconcile this? A "Left Hand Sum" would omit t=T from the time points; a "Right Hand Sum" would omit t=0. The trapezoidal rule omits half a point from each.

That sounds awfully strange, but it's mathematically sound! We only integrate through half of each boundary time point t=0 and t=T. Thus, those points, and only those points, have a spacing of τ/2.

TemporalLattice(T::Real, r::Int)

Convenience constructor to make a trapezoidal time grid from t=0.0 to T.

julia> grid = TemporalLattice(20.0, 400);

julia> validate(grid);

julia> Integrations.duration(grid)
20.0
julia> Integrations.nsteps(grid)
400

ComplexConstantSignal(A::F, B::F) where {F<:AbstractFloat}

The constant complex signal $Ω(t) = A + iB$.

ConstantSignal(A::AbstractFloat)

The constant real signal $Ω(t) = A$.

Constant(A)
Constant(A, B)

Convenience constructors for a constant signal.

The single argument form constructs a ConstantSignal when given a real number, or a ComplexConstantSignal when given a complex number. The two-argument form constructs a ComplexConstantSignal, taking A as the real part and B as the imaginary part.

julia> grid = TemporalLattice(20.0, 400);

julia> realonly = Constant(2.0);

julia> validate(realonly; grid=grid, t=10.0, rms=1e-6);

julia> valueat(realonly, 10.0)
2.0
julia> Parameters.names(realonly)
1-element Vector{String}:
 "A"

julia> complex = Constant(2.0, 1.0);

julia> validate(complex; grid=grid, t=10.0, rms=1e-6);

julia> valueat(complex, 10.0)
2.0 + 1.0im
julia> Parameters.names(complex)
2-element Vector{String}:
 "A"
 "B"

julia> alsocomplex = Constant(2.0 + 1.0im);

julia> validate(alsocomplex; grid=grid, t=10.0, rms=1e-6);

julia> typeof(alsocomplex) == typeof(complex)
true
julia> valueat(alsocomplex, 10.0) == valueat(complex, 10.0)
true
julia> Parameters.names(alsocomplex) == Parameters.names(complex)
true

julia> imagfrozen = Constrained(complex, :B);

julia> validate(imagfrozen; grid=grid, t=10.0, rms=1e-6);

julia> valueat(imagfrozen, 10.0) == valueat(complex, 10.0)
true
julia> Parameters.names(imagfrozen) == Parameters.names(realonly)
true

CompositeSignal(components::AbstractVector{<:SignalType{P,R}})

A signal which is the sum of each sub-signal in components.

Each component should be the same type of signal, for type stability. If you need to compose different types of signals, you should probably implement your own custom SignalType.

Composed(components::SignalType...)

Convenience constructor to combine multiple signals into a CompositeSignal.

julia> grid = TemporalLattice(20.0, 400);

julia> realpart = Constrained(Constant(2.0+0.0im), :B);

julia> imagpart = Constrained(Constant(0.0+1.0im), :A);

julia> signal = Composed(realpart, imagpart);

julia> validate(signal; grid=grid, t=10.0, rms=1e-6);

julia> Parameters.names(signal)
2-element Vector{String}:
 "A.1"
 "B.2"

WindowedSignal(windows, starttimes)

A signal which applies a different function for each window.

Arguments

• windows: a vector of signals
• starttimes: a vector of times transitioning each window

Both windows and starttimes have the same length; starttimes[i] indicates when windows[i] begins.

This signal is undefined for times t < starttimes[1]. Normally, starttimes[1] == 0.

Windows with dynamically changing numbers of parameters are unsupported.

Windowed(signal, starttimes)
Windowed(signal, T, W)

Convenience constructors to segment a single signal into a WindowedSignal.

Feed in starttimes directly, or make W uniformly spaced windows up to maximum time T.

julia> grid = TemporalLattice(20.0, 400);

julia> signal = Windowed(Constant(2.0), 20.0, 5);

julia> validate(signal; grid=grid, t=10.0, rms=1e-6);

julia> Parameters.names(signal)
5-element Vector{String}:
 "A.1"
 "A.2"
 "A.3"
 "A.4"
 "A.5"

get_window_from_parameter(signal::WindowedSignal, i::Int)

Identify the window index given a parameter index (by counting parameters in windows).

get_window_from_time(signal::WindowedSignal, t::Real)

Identify the window index given the time (by inspecting starttimes).

RealWindowedResonantTransmonDevice{m}(ω, δ, g, quples, T, Ω)

A minimalist transmon device with real-valued constant windows driven on resonance.

Windows are always equally spaced. Drives are approximated with RWA.

Type Parameters

• F: (inferred from arguments) the float type of this device
• m: the integer number of levels per transmon

Parameters

• ω: an (abstract) vector of qubit resonance frequencies
• δ: an (abstract) vector of qubit anharmonicities
• g: an (abstract) vector of qubit coupling strengths
• quples: an (abstract) vector of quples identifying each coupling
• T: the total pulse duration applied on this device
• Ω: a matrix of complex amplitudes. Ω[w,q] is the amplitude applied on qubit q in the time-window indexed by w.

RealWindowedResonantTransmonDevice{m}(ω, δ, g, quples, T, W::Int)

Convenience constructor where Ω is initialized to zero, with W time windows.

Prototype(::Type{RealWindowedResonantTransmonDevice{F,m}}, n::Int; kwargs...)

A prototypical RealWindowedResonantTransmonDevice with the following decisions:

• All anharmonicities are constant.
• Couplings are linear.
• Each coupling strength equals the difference in resonance frequencies of the coupled qubits.
• By default, all resonance frequencies are equally spaced (so, coupling strengths are constant) but this can be controlled through kwargs.

Default parameters are vaguely reminiscent of IBM devices circa 2021, although the default behavior of linearly-spaced resonance frequencies is not realistic and should be avoided outside of testing/benchmarking.

Keyword Arguments

• ω0=4.82: resonance frequency of first qubit.
• Δω=0.02: the spacing in resonance frequencies between adjacent qubits. When passed as a float (including the default), resonance frequencies are linearly spaced. Instead, you can pass this as an explicit vector with n-1 elements.
• δ0=0.30: the constant anharmonicity.
• T=10.0: total pulse duration.
• W=1: number of window segments.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(RWRTDevice{Float64,3}; n=2);

julia> validate(device; grid=grid, t=10.0);

julia> nlevels(device)
3
julia> nqubits(device)
2

__windowindex(device, t)

Infer the column index of device.Ω i which time t falls.

WindowedResonantTransmonDevice{m}(ω, δ, g, quples, T, Ω)

A minimalist transmon device with complex-constant windows driven on resonance.

Windows are always equally spaced. Drives are approximated with RWA.

Type Parameters

• F: (inferred from arguments) the float type of this device
• m: the integer number of levels per transmon

Parameters

• ω: an (abstract) vector of qubit resonance frequencies
• δ: an (abstract) vector of qubit anharmonicities
• g: an (abstract) vector of qubit coupling strengths
• quples: an (abstract) vector of quples identifying each coupling
• T: the total pulse duration applied on this device
• Ω: a matrix of complex amplitudes. Ω[w,q] is the amplitude applied on qubit q in the time-window indexed by w.

WindowedResonantTransmonDevice{m}(ω, δ, g, quples, T, W::Int)

Convenience constructor where Ω is initialized to zero, with W time windows.

Prototype(::Type{WindowedResonantTransmonDevice{F,m}}, n::Int; kwargs...)

A prototypical WindowedResonantTransmonDevice with the following decisions:

• All anharmonicities are constant.
• Couplings are linear.
• Each coupling strength equals the difference in resonance frequencies of the coupled qubits.
• By default, all resonance frequencies are equally spaced (so, coupling strengths are constant) but this can be controlled through kwargs.

Default parameters are vaguely reminiscent of IBM devices circa 2021, although the default behavior of linearly-spaced resonance frequencies is not realistic and should be avoided outside of testing/benchmarking.

Keyword Arguments

• ω0=4.82: resonance frequency of first qubit.
• Δω=0.02: the spacing in resonance frequencies between adjacent qubits. When passed as a float (including the default), resonance frequencies are linearly spaced. Instead, you can pass this as an explicit vector with n-1 elements.
• δ0=0.30: the constant anharmonicity.
• T=10.0: total pulse duration.
• W=1: number of window segments.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(CWRTDevice{Float64,3}; n=2);

julia> validate(device; grid=grid, t=10.0);

julia> nlevels(device)
3
julia> nqubits(device)
2

__windowindex(device, t)

Infer the column index of device.Ω i which time t falls.

TransmonDevice{m}(ω, δ, g, quples, q, Ω, Δ)

A transmon device where control signals and drive frequencies are signals.

Drives are approximated with RWA.

Type Parameters

• F: (inferred from arguments) the float type of this device
• m: the integer number of levels per transmon

Parameters

• ω: an (abstract) vector of qubit resonance frequencies
• δ: an (abstract) vector of qubit anharmonicities
• g: an (abstract) vector of qubit coupling strengths
• quples: an (abstract) vector of quples identifying each coupling
• q: the qubits corresponding to each drive.
• Ω: an (abstract) vector of control signals for each drive. May be real or complex.
• Δ: an (abstract) vector of detunings for each drive frequency. Must be real.

TransmonDevice{m}(ω, δ, g, quples; kwargs...)

Convenience constructor allowing a more semantic approach to inputting signals.

Keyword Arguments

• q: may be a an (abstract) vector of Int. Defaults to one drive for each qubit (i.e. 1:length(ω))
• Ω: may be a SignalType or an (abstract) vector of SignalTypes. If a single SignalType is provided, it is duplicated for each drive.
• Δ: may be a Bool or a SignalType or an (abstract) vector of SignalTypes. If a single SignalType is provided, it is duplicated for each drive. If Δ=true, drive frequencies are constant signals initialized on resonance. If Δ=false (the default), drive frequencies are frozen on resonance.

Prototype(::Type{TransmonDevice{F,m}}, n::Int; kwargs...)

A prototypical TransmonDevice with the following decisions:

• All anharmonicities are constant.
• Couplings are linear.
• Each coupling strength equals the difference in resonance frequencies of the coupled qubits.
• By default, all resonance frequencies are equally spaced (so, coupling strengths are constant) but this can be controlled through kwargs.
• Drives match those of the defaults when using the kwarg constructor.

Default parameters are vaguely reminiscent of IBM devices circa 2021, although the default behavior of linearly-spaced resonance frequencies is not realistic and should be avoided outside of testing/benchmarking.

Keyword Arguments

• ω0=4.82: resonance frequency of first qubit.
• Δω=0.02: the spacing in resonance frequencies between adjacent qubits. When passed as a float (including the default), resonance frequencies are linearly spaced. Instead, you can pass this as an explicit vector with n-1 elements.
• δ0=0.30: the constant anharmonicity.
• T=10.0: the pulse duration, but this has no effect since the default signals are constant.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,3}; n=2);

julia> validate(device; grid=grid, t=10.0);

julia> nlevels(device)
3
julia> nqubits(device)
2

ROTATING_FRAME

A Trotterization method calculating drive terms in the rotating frame of the device.

The work basis for this algorithm is Bases.DRESSED. This ensures the rotating-frame evolution $U_t ≡ exp(-itH_0)$ is quite cheap. Even so, this algorithm exponentiates the matrix $U_t' V(t) U_t$ at each time step, so it is not terribly efficient.

A gradientsignals method is not currently supported for this evolution algorithm.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,2}; n=2);

julia> evolution = ROTATING_FRAME;

julia> validate(evolution; grid=grid, device=device, skipgradient=true);

julia> workbasis(evolution)
CtrlVQE.Bases.Dressed()

QUBIT_FRAME

A Trotterization method alternately propagating static and drive terms.

The work basis for this algorithm is Bases.BARE. The static term propagator is expensive but only computed once. If the drive terms are local (as for a LocallyDrivenDevice), the drive propagator is relatively cheap.

Beware that this algorithm implicitly employs a trapezoidal rule, irrespective of the grid passed to evolution functions.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,2}; n=2);

julia> evolution = QUBIT_FRAME;

julia> validate(evolution; grid=grid, device=device, skipgradient=true);

julia> workbasis(evolution)
CtrlVQE.Bases.Bare()

DenseLeakage(reference, device, basis, frame, grid, evolution)

Calculate leakage of a reference state after evolution.

Leakage is the probability of finding any qubit outside the |0⟩,|1⟩ subspace.

Parameters

• reference: the initial statevector before evolution.
• device: the device under which the state evolves.
• basis: the basis that reference is input as, and the basis which for which leakage is defined.
• frame: the rotating frame for which leakage is defined.
• grid: the time grid on which the state evolves.
• evolution: the algorithm to calculate the time evolution.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,2}; n=2);

julia> ψ0 = LAT.basisvector(Complex{eltype(device)}, nstates(device), 1);

julia> costfn = DenseLeakage(ψ0, device, Bases.DRESSED, Operators.STATIC, grid, QUBIT_FRAME);

julia> validate(costfn; rms=1e-6, grid=grid, device=device);

DenseObservable(observable, reference, device, basis, frame, grid, evolution)

Calculate the expectation value of an observable after evolution of a reference state.

Parameters

• observable: the matrix (defined on the whole Hilbert space) to measure.
• reference: the initial statevector before evolution.
• device: the device under which the state evolves.
• basis: the basis that observable and reference are input as.
• frame: the rotating frame in which observable is to be measured.
• grid: the time grid on which the state evolves.
• evolution: the algorithm to calculate the time evolution.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,2}; n=2);

julia> ψ0 = LAT.basisvector(Complex{eltype(device)}, nstates(device), 1);

julia> O = LAT.basisvectors(Complex{eltype(device)}, nstates(device));

julia> costfn = DenseObservable(O, ψ0, device, Bases.DRESSED, Operators.STATIC, grid, QUBIT_FRAME);

julia> validate(costfn; rms=1e-6, grid=grid, device=device);

CompositeCostFunction(components::AbstractVector{CostFunctionType})

The sum of several cost-functions, with matching length and float type.

Use this eg. to combine an energy function with one or more penalty functions.

This struct records the last values computed (value or gradient) for each cost function in its values and gradients fields.

Beware that the components field is a vector with abstract eltype. This is a relatively mild form of type instability; certain compiled code could hypothetically be sub-optimal. But it does not appear to cause any noticeable disadvantage; see Examples/CompositeCostFunctions.

CompositeCostFunction(components::CostFunctionType{F}...)

Alternate constructor, letting each function be passed as its own argument.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(CWRTDevice{Float64,2}; n=2);

julia> ψ0 = LAT.basisvector(Complex{eltype(device)}, nstates(device), 1);

julia> O = LAT.basisvectors(Complex{eltype(device)}, nstates(device));

julia> energyfn = DenseObservable(O, ψ0, device, Bases.DRESSED, Operators.STATIC, grid, QUBIT_FRAME);

julia> penaltyfn = WindowedResonantPenalty(device; A=0.8);

julia> costfn = CompositeCostFunction(energyfn, penaltyfn);

julia> validate(costfn; rms=1e-6);

WindowedResonantPenalty(device, A::Vector, σ::Vector, λ::Vector)

Compute penalties for exceeding the maximum amplitude of a device.

Let $u=rac{|Ω|-A}{σ}$. When $u>0$ the penalty is computed as $Λ = λ xp(u - 1/u)$.

Parameters

• device: a compatible device (see below).
• A: the maximum modulus of any amplitude, on each drive.
• σ: the steepness parameter of the penalty function, on each drive.
• λ: the strength parameter of the penalty function, on each drive

Device Compatibility

Intended for use with the RWRTDevice and CWRTDevice provided in the basics, this type should work if:

• device has a field Ω
• Ω is a matrix of (potentially complex) amplitudes, with n columns, where n = Devices.nqubits(device).
• Parameters.values(x) is a vectorization of Ω (and nothing else), reinterpreted as a vector of floats when Ω is complex.

WindowedResonantPenalty(device; A=1.0, σ=A, λ=1.0)

Alternate constructor, with kwargs and sensible defaults.

Each parameter may be provided as a vector or as a scalar, in which case it is automatically expanded to a vector where the given value is assigned to each qubit.

We have selected σ defaulting to A because heuristically it seems to work well. Don't feel too attached to that choice.

julia> device = Prototype(CWRTDevice{Float64,2}; n=2);

julia> costfn = WindowedResonantPenalty(device; A=0.8);

julia> x = collect(range(0.0, 1.0, length(costfn)))
4-element Vector{Float64}:
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0

julia> validate(costfn; x=x, rms=1e-6);

julia> costfn(x)
0.22571605879846993
julia> grad_function(costfn)(x)
4-element Vector{Float64}:
 0.0
 0.0
 0.7767775331625435
 1.1651662997438155

SignalStrengthPenalty(grid, signal, A, σ, λ)

Compute penalties for a signal exceeding a maximum absolute value over a time interval.

Let $u(t)=\frac{|Ω(t)|-A}{σ}$, where Ω(t) is defined by signal. When $u>0$ the instantaneous penalty is computed as $Λ(t) = λ⋅\exp[u(t) - 1/u(t)]$. The total penalty is $\frac{1}{T}⋅\int_0^T Λ(t)⋅dt$, where the integration is defined by grid.

Parameters

• grid: defines the integration (including time bounds) over which to penalize.
• signal: the signal to penalize.
• A: the maximum modulus of the signal.
• σ: the steepness parameter of the penalty function.
• λ: the strength parameter of the penalty function.

Note that, if you do not wish the penalty to be normalized per unit time, λ should be selected to be proportional to the duration of grid.

SignalStrengthPenalty(device; A=1.0, σ=A, λ=1.0)

Alternate constructor, with kwargs and sensible defaults.

We have selected σ defaulting to A because heuristically it seems to work well. Don't feel too attached to that choice.

julia> grid = TemporalLattice(20.0, 400);

julia> signal = Windowed(Constant(2.0), 20.0, 5);

julia> costfn = SignalStrengthPenalty(grid, signal; A=0.8);

julia> x = collect(range(0.0, 1.0, length(costfn)))
5-element Vector{Float64}:
 0.0
 0.25
 0.5
 0.75
 1.0

julia> validate(costfn; x=x, rms=1e-6);

julia> costfn(x)
0.00473294635352182
julia> grad_function(costfn)(x)
5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.10057511001233899

AmplitudePenalty(device::DeviceType, penalties::AbstractVector{<:CostFunctionType})

The sum of several penalty functions acting on a device's amplitude signals.

In order for this type to be usable, each penalty in penalties must act on parameters corresponding to the signals contained in device.Ω, and parameters in device must be disjoint, listing parameters for each device.Ω then each device.ν. This interface is suitable for the basic TransmonDevice.

This struct records the last values computed (value or gradient) for each cost function in its values and gradients fields.

Beware that the components field is a vector with abstract eltype. This is a relatively mild form of type instability; certain compiled code could hypothetically be sub-optimal. But it does not appear to cause any noticeable disadvantage; see Examples/CompositeCostFunctions.

julia> grid = TemporalLattice(20.0, 400);

julia> device = Prototype(TransmonDevice{Float64,2}; n=2);

julia> penalties = [SignalStrengthPenalty(grid, signal; A=0.8) for signal in device.Ω];

julia> costfn = AmplitudePenalty(device, penalties);

julia> x = collect(range(0.0, 1.0, length(costfn)))
4-element Vector{Float64}:
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0

julia> validate(costfn; x=x, rms=1e-6);

julia> costfn(x)
0.22571605879846895
julia> grad_function(costfn)(x)
4-element Vector{Float64}:
 0.0
 0.0
 0.7767775331625492
 1.165166299743822

DetuningPenalty(device::DeviceType, penalties::AbstractVector{<:CostFunctionType})

The sum of several penalty functions acting on a device's detuning signals.

In order for this type to be usable, each penalty in penalties must act on parameters corresponding to the signals contained in device.Δ, and parameters in device must be disjoint, listing parameters for each device.Ω then each device.Δ. This interface is suitable for the basic TransmonDevice.

This struct records the last values computed (value or gradient) for each cost function in its values and gradients fields.

Beware that the components field is a vector with abstract eltype. This is a relatively mild form of type instability; certain compiled code could hypothetically be sub-optimal. But it does not appear to cause any noticeable disadvantage; see Examples/CompositeCostFunctions.

julia> grid = TemporalLattice(20.0, 400);

julia> device = TransmonDevice{2}([4.82, 4.84], [0.00, 0.00], [0.02], [Quple(1,2)]; Ω=Constrained(Constant(zero(ComplexF64)), :B), Δ=Constant(zero(Float64)));

julia> penalties = [SignalStrengthPenalty(grid, signal; A=0.8) for signal in device.Δ];

julia> costfn = DetuningPenalty(device, penalties);

julia> x = collect(range(0.0, 1.0, length(costfn)))
4-element Vector{Float64}:
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0

julia> validate(costfn; x=x, rms=1e-6);

julia> costfn(x)
0.02351774585600897
julia> grad_function(costfn)(x)
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.4997520994401937