minCoeff, maxCoeff = assign_thresholds(params, D1, maxpar [, maxpar_unc])

Build vector of frequency independent min/max parameter constraints for each control function. Here, minCoeff = -maxCoeff.

Arguments

• params:: objparams: Struct containing problem definition.
• D1:: Int64: Number of basis functions in each segment.
• maxpar::Array{Float64,1}: Maximum parameter value for each coupled control.

minCoeff, maxCoeff = assign_thresholds_ctrl_freq(params, D1, maxamp)

Build vector of parameter min/max constraints that can depend on the control function and carrier wave frequency, with minCoeff = -maxCoeff.

Arguments

• params:: objparams: Struct containing problem definition.
• D1:: Int64: Number of basis functions in each segment.
• maxamp:: Matrix{Float64}: maxamp[c,f] is the maximum parameter value for ctrl c and frequency f

minCoeff, maxCoeff = assign_thresholds_freq(maxamp, Ncoupled, Nfreq, D1)

Build vector of frequency dependent min/max parameter constraints, with minCoeff = -maxCoeff, when there are no uncoupled control functions.

Arguments

• maxamp::Array{Float64,1}: Maximum parameter value for each frequency
• Ncoupled::Int64: Number of coupled controls in the simulation
• Nfreq::Int64: Number of carrier wave frequencies used in the controls
• D1:: Int64: Number of basis functions in each control function

f = bcarrier2(t, params, func)

Evaluate a B-spline function with carrier waves. See also the bcparams constructor.

Arguments

• t::Float64: Evaluate spline at parameter t ∈ [0, param.T]
• param::params: Parameters for the spline
• func::Int64: Spline function index ∈ [0, param.Nseg-1]

f = bspline2(t, splineparam, splinefunc)

Evaluate a B-spline function. See also the splineparams constructor.

Arguments

• t::Float64: Evaluate spline at parameter t ∈ [0, param.T]
• param::splineparams: Parameters for the spline
• splinefunc::Int64: Spline function index ∈ [0, param.Nseg-1]

nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, maxpar [, Pmin = 40])

Estimate the number of time steps needed for the simulation, for the case without uncoupled controls.

Arguments

• T:: Float64: Final simulation
• H0::Array{Float64,2}: Time-independent part of the Hamiltonian matrix
• Hsym_ops:: Array{Float64,2}: Array of symmetric control Hamiltonians
• Hanti_ops:: Array{Float64,2}: Array of symmetric control Hamiltonians
• maxpar:: Array{Float64,1}: Maximum parameter value for each subsystem
• Pmin:: Int64: Number of time steps per shortest period (assuming a slowly varying Hamiltonian).

nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, Hunc_ops, 
                                          maxpar, max_flux[, Pmin = 40])

Estimate the number of time steps needed for the simulation, when there are uncoupled controls.

Arguments

• T:: Float64: Final simulation
• H0::Array{Float64,2}: Time-independent part of the Hamiltonian matrix
• Hsym_ops:: Array{Float64,2}: Array of symmetric control Hamiltonians
• Hanti_ops:: Array{Float64,2}: Array of symmetric control Hamiltonians
• Hunc_ops:: Array{Float64,2}: Array of uncoupled control Hamiltonians
• maxpar:: Array{Float64,1}: Maximum parameter value for each coupled control
• max_flux:: Array{Float64,1}: Maximum parameter value for each uncoupled control
• Pmin:: Int64: Number of time steps per shortest period (assuming a slowly varying Hamiltonian).

nsteps = calculate_timestep(T, H0, Hunc_ops, maxpar [, Pmin = 40])

Estimate the number of time steps needed for an accurate simulation, when there are no coupled controls

Arguments

• T:: Float64: Final simulation
• H0::Array{Float64,2}: Time-independent part of the Hamiltonian matrix
• Hunc_ops:: Array{Float64,2}: Array of uncoupled control Hamiltonians
• max_unc:: Array{Float64,1}: Maximum parameter value for each subsystem (uncoupled)
• Pmin:: Int64: Sample rate for accuracy (assuming a slowly varying Hamiltonian)

change_target!(params, new_Utarget)

Update the unitary target in the objparams object.

Arguments

• param::objparams: Object holding the problem definition
• new_Utarget::Array{ComplexF64,2}: New unitary target as a two-dimensional complex-valued array (matrix) of dimension Ntot x N

estimate_Neumann!(tol, params, maxpar[, maxunc])

Estimate the number of terms needed by the Neumann series approach for solving the linear system during the implicit steps of the Störmer-Verlet scheme. See also neumann!

Arguments

• tol:: Float64: Error tolerance in inverting implicit SV term
• params:: objparams: Struct containing problem definition
• maxpar:: Array{Float64,1}: Maximum parameter value for each coupled control
• maxunc:: Array{Float64,1}: (optional) Maximum parameter value for each uncoupled controls

pj [, qj] = evalctrl(params, pcof0, td, func)

Evaluate the control function with index func at an array of time levels td.

NOTE: the control function index func is 1-based.

NOTE: The return value(s) depend on func. For func∈[1,Ncoupled], pj, qj are returned. Otherwise, only pj is returned, corresponding to control number func.

Arguments

• params:: objparams: Struct with problem definition
• pcof0:: Array{Float64,1}): Vector of parameter values
• td:: Array{Float64,1}): Time values control is to be evaluated
• jFunc:: Int64: Index of the control signal desired

For comparing other solvers with the built-in linear solver using Gaussian elimination.

gamma, used_stages = getgamma(order, stages)

Obtain step size coefficients (gamma) for the composition method with the given order and number of stages.

Arguments:

• order: The desired order of the compositional method. Should be 2, 4, 6, 8,

or 10.

• stages: The desired number of stages of the compositional method. Default

value will be used if input is invalid. '0' is reserved for using a default number of stages.

gradbcarrier2!(t, params, func, g) -> g

Evaluate the gradient of a control function with respect to all coefficient.

NOTE: the index of the control functions is 0-based. For a set of coupled controls, mod(func,2)=0 corresponds to ∇ pj(t) and mod(func,2) = 1 corresponds to ∇ qj(t), where j = div(func,2).

Arguments

• t::Float64: Evaluate spline at parameter t ∈ [0, param.T]
• params::bcparams: Parameters for the spline
• func::Int64: Control function index ∈ [0, param.Nseg-1]
• g::Array{Float64,1}: Preallocated array to store calculated gradient

g = gradbspline2(t, param, splinefunc)

Evaluate the gradient of a spline function with respect to all coefficient. NOTE: the index of the spline functions are 0-based. For a set of coupled controls, mod(splinefunc,2)=0 corresponds to ∇ pj(t) and mod(splinefunc,2) = 1 corresponds to ∇ qj(t), where j = div(splinefunc,2).

Arguments

• t::Float64: Evaluate spline at parameter t ∈ [0, param.T]
• param::splineparams: Spline parameter object
• splinefunc::Int64: Spline function index ∈ [0, param.Nseg-1]

forbiddenlev = identify_forbidden_levels(params[, custom = 0])

Build a Bool array indicating which energy levels are forbidden levels in the state vector. The forbidden levels in a state vector are defined as thos corresponding to the highest energy level in at least one of its subsystems.

Arguments

• params:: objparams: Struct with problem definition
• custom:: Int64: For nonzero value special stirap pulses case

guardlev = identify_guard_levels(params[, custom = 0])

Build a Bool array indicating if a given energy level is a guard level in the simulation.

Arguments

• params:: objparams: Struct with problem definition
• custom:: Int64: A nonzero value gives a special stirap pulses case

u_init = initial_cond(Ne, Ng)

Setup a basis of canonical unit vectors that span the essential Hilbert space, setting all guard levels to zero

Arguments

• Ne:: Array{Int64}: Array holding the number of essential levels in each system
• Ng:: Array{Int64}: Array holding the number of guard levels in each system

juq2qis(params, pcof, samplerate, q_ind, fileName="ctrl.dat", node_loc="c")

Evaluate control functions and export them into a format that is readable by Qiskit.

Arguments

• params:: objparams: Struct with problem definition
• pcof:: Array{Float64,1}): Vector of parameter values
• samplerate:: Float64: Samplerate for quantum device (number of samples per ns for the IQ mixer)
• q_ind:: Int64: Index of the control function to output (1 <= q_ind <= Nctrl*Nfreq)
• fileName:: String: Name of output file containing controls to be handled by Qiskit
• node_loc:: String: Node location, "c" for cell centered, "n" for node-centered, default is cell-centered

marg_prob = marginalize3(params, unitaryhist)

Evaluate marginalized probabilities for the case of 3 subsystems.

Arguments

• param:: objparams: Struct with problem definition
• unitaryhist:: Array{Complex{Float64},3}): State vector history for each timestep

wmat = orig_wmatsetup(Ne, Ng)

NOTE: This an alternative version of wmatsetup() that uses slightly different coefficients.

Build the default positive semi-definite weighting matrix W to calculate the leakage into higher energy forbidden states

Arguments

• Ne::Array{Int64,1}: Number of essential energy levels for each subsystem
• Ng::Array{Int64,1}: Number of guard energy levels for each subsystem

pconv = plot_conv_hist(params [, convname:: String=""])

Plot the optimization convergence history, including history of the different terms in the objective function and the norm of the gradient.

Arguments

• param:: objparams: Struct with problem definition
• convname:: String: Name of plot file to be generated

plt =  plot_energy(unitaryhistory, params)

Plot the evolution of the expected energy for each initial condition.

Arguments

• unitaryhistory:: Array{ComplexF64,3}: Array holding the time evolution of the state for each initial condition
• params:: objparams: Struct with problem definition

plt =  plot_final_unitary(final_unitary, params)

Plot the essential levels of the solution operator at a fixed time and return a plot handle

Arguments

• final_unitary:: Array{ComplexF64,2}: Ntot by Ness array holding the final state for each initial condition
• params:: objparams: Struct with problem definition
• fid:: Float64: Gate fidelity (for plot title)

pl = plot_results(params, pcof; [casename = "test", savefiles = false, samplerate = 32])

Create array of plot objects that can be visualized by, e.g., display(pl[1]).

Arguments

• params::objparams: Object holding problem definition
• pcof::Array{Float64,1}: Parameter vector
• casename::String: Default: "test". String used in plot titles and in file names
• savefiles::Bool: Default: false.Set to true to save plots on files with automatically generated filenames
• samplerate:: Int64: Default: 32 samples per unit time (ns). Sample rate for generating plots.

plt = plotspecified(us, params, guardlev, specifiedlev)

Plot the evolution of the state vector for specified levels.

Arguments

• us:: Array{Complex{Float64},3}): State vector history for each timestep
• param:: objparams: Struct with problem definition
• guardlev:: Array{Bool,1}): Boolean array indicating if a certain level is a guard level
• specifiedlev:: Array{Bool,1}: Boolean array indicating which levels to be plotted

plt = plotunitary(us, params, guardlev)

Plot the evolution of the state vector.

Arguments

• us:: Array{Complex{Float64},3}): State vector history for each timestep
• param:: objparams: Struct with problem definition
• guardlev:: Array{Bool,1}): Boolean array indicating if a certain level is a guard level

pcof = read_pcof(refFileName)

Read the parameter vector pcof from a JLD2 formatted file

Arguments

• refFileName: String holding the name of the file.

pcof = run_optimizer(prob, pcof0 [, baseName:: String=""])

Call IPOPT to optimizize the control functions.

Arguments

• prob:: IpoptProblem: Struct containing Ipopt parameters callback functions
• pcof0:: Array{Float64, 1}: Initial guess for the parameter values
• baseName:: String: Name of file for saving the optimized parameters; extension ".jld2" is appended

save_pcof(refFileName, pcof)

Save the parameter vector pcof on a JLD2 formatted file with handle pcof

Arguments

• refFileName: String holding the name of the file.
• pcof: Vector of floats holding the parameters.

set_adjoint_Sv_type!(params, new_sv_type)

For continuation only: update the sv_type in the objparams object.

Arguments

• param::objparams: Object holding the problem definition
• new_sv_type:: Int64: New value for sv_type. Must be 1, 2, or 3.

prob = setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff; maxIter=50, 
                        lbfgsMax=10, startFromScratch=true, ipTol=1e-5, acceptTol=1e-5, acceptIter=15,
                        nodes=[0.0], weights=[1.0])

Setup structure containing callback functions and convergence criteria for optimization via IPOPT. Note the last two inputs, nodes', andweights', are to be used when performing a simple risk-neutral optimization where the fundamental frequency is random.

Arguments

• params:: objparams: Struct with problem definition
• wa::Working_Arrays: Struct containing preallocated working arrays
• nCoeff:: Int64: Number of parameters in optimization
• minCoeff:: Array{Float64, 1}: Minimum allowable value for each parameter
• maxCoeff:: Array{Float64, 1}: Maximum allowable value for each parameter
• maxIter:: Int64: Maximum number of iterations to be taken by optimizer (keyword arg)
• lbfgsMax:: Int64: Maximum number of past iterates for Hessian approximation by L-BFGS (keyword arg)
• startFromScratch:: Bool: Specify whether the optimization is starting from file or not (keyword arg)
• ipTol:: Float64: Desired convergence tolerance (relative) (keyword arg)
• acceptTol:: Float64: Acceptable convergence tolerance (relative) (keyword arg)
• acceptIter:: Int64: Number of acceptable iterates before triggering termination (keyword arg)
• nodes:: Array{Float64, 1}: Risk-neutral opt: User specified quadrature nodes on the interval [-ϵ,ϵ] for some ϵ (keyword arg)
• weights:: Array{Float64, 1}: Risk-neutral opt: User specified quadrature weights on the interval [-ϵ,ϵ] for some ϵ (keyword arg)

omega1[, omega2, omega3] = setup_rotmatrices(Ne, Ng, fund_freq)

Build diagonal rotation matrices based on the |0⟩to |1⟩ transition frequency in each sub-system.

Arguments

• Ne::Array{Int64,1}: Number of essential energy levels for each subsystem
• Ng::Array{Int64,1}: Number of guard energy levels for each subsystem
• fund_freq::Array{Float64}: Transitions frequency [GHz] for each subsystem

objf = traceobjgrad(pcof0, params, wa[, verbose = false, evaladjoint = true])

Perform a forward and/or adjoint Schrödinger solve to evaluate the objective function and/or gradient.

Arguments

• pcof0::Array{Float64,1}: Array of parameter values defining the controls
• param::objparams: Struct with problem definition
• wa::Working_Arrays: Struct containing preallocated working arrays
• verbose::Bool = false: Run simulation with additional terminal output and store state history.
• evaladjoint::Bool = true: Solve the adjoint equation and calculate the gradient of the objective function.

wmat = wmatsetup(Ne, Ng)

Build the default positive semi-definite weighting matrix W to calculate the leakage into higher energy forbidden states

Arguments

• Ne::Array{Int64,1}: Number of essential energy levels for each subsystem
• Ng::Array{Int64,1}: Number of guard energy levels for each subsystem

zero_start_end!(params, D1, minCoeff, maxCoeff)

Force the control functions to start and end at zero by setting zero bounds for the first two and last two parameters in each B-spline segment.

Arguments

• params:: objparams: Struct containing problem definition.
• D1:: Int64: Number of basis functions in each segment.
• minCoeff:: Vector{Float64}: Lower parameter bounds to be modified
• maxCoeff:: Vector{Float64}: Upper parameter bounds to be modified