diff --git a/Project.toml b/Project.toml index 23d88e1d..939f1d8f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "QuantumCollocation" uuid = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf" authors = ["Aaron Trowbridge and contributors"] -version = "0.6.0" +version = "0.7.0" [deps] +DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" @@ -12,7 +13,6 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08" PiccoloQuantumObjects = "5a402ddf-f93c-42eb-975e-5582dcda653d" -QuantumCollocationCore = "2b384925-53cb-4042-a8d2-6faa627467e1" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -21,14 +21,14 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8" [compat] +DirectTrajOpt = "0.1.0" Distributions = "0.25" ExponentialAction = "0.2" Interpolations = "0.15" JLD2 = "0.5" LinearAlgebra = "1.10, 1.11" -NamedTrajectories = "0.2" +NamedTrajectories = "0.3" PiccoloQuantumObjects = "0.3" -QuantumCollocationCore = "0.3" Random = "1.10, 1.11" Reexport = "1.2" SparseArrays = "1.10, 1.11" diff --git a/docs/Project.toml b/docs/Project.toml index 7e07b254..8df91abc 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,9 +3,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08" +PiccoloQuantumObjects = "5a402ddf-f93c-42eb-975e-5582dcda653d" QuantumCollocation = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf" -QuantumCollocationCore = "2b384925-53cb-4042-a8d2-6faa627467e1" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" - -[compat] -NamedTrajectories = "0.2" \ No newline at end of file diff --git a/docs/literate/man/ipopt_callbacks.jl b/docs/dev/ipopt_callbacks.jl similarity index 95% rename from docs/literate/man/ipopt_callbacks.jl rename to docs/dev/ipopt_callbacks.jl index c69a1876..fa072deb 100644 --- a/docs/literate/man/ipopt_callbacks.jl +++ b/docs/dev/ipopt_callbacks.jl @@ -1,9 +1,9 @@ # ```@meta # CollapsedDocStrings = true # ``` -# # IpOpt Callbacks +# # IPOPT Callbacks -# This page describes the callback functions that can be used with the IpOpt solver (in the future, may describe more general callback behavior). +# This page describes the callback functions that can be used with the IPOPT solver (in the future, may describe more general callback behavior). # ## Callbacks @@ -13,7 +13,7 @@ using NamedTrajectories import ..QuantumStateSmoothPulseProblem import ..Callbacks -# By default, IpOpt callbacks are called at each optimization step with the following signature: +# By default, IPOPT callbacks are called at each optimization step with the following signature: function full_argument_list_callback( alg_mod::Cint, iter_count::Cint, diff --git a/docs/literate/man/problem_templates.jl b/docs/literate/man/problem_templates.jl index a5d12f9d..75107fd0 100644 --- a/docs/literate/man/problem_templates.jl +++ b/docs/literate/man/problem_templates.jl @@ -1,19 +1,197 @@ # ```@meta # CollapsedDocStrings = true # ``` +using NamedTrajectories +using PiccoloQuantumObjects +using QuantumCollocation -# # Problem Templates +#= +# Problem Templates -# We provide a number of problem templates for making it simple and easy to set up and solve certain types of quantum optimal control problems. These templates all construct a `QuantumControlProblem` object. The problem templates are: +We provide a number of problem templates for making it simple and easy to set up and solve +certain types of quantum optimal control problems. These templates all construct a +`DirectTrajOptProblem` object, which stores all the parts of the optimal control problem. -# ## Unitary Smooth Pulse Problem +This page provides a brief overview of each problem template, broken down by the state of +the problem being solved. -# ```@docs; canonical = false -# UnitarySmoothPulseProblem -# ``` +Ket Problem Templates: +- [Quantum State Smooth Pulse Problem](#Quantum-State-Smooth-Pulse-Problem) +- [Quantum State Minimum Time Problem](#Quantum-State-Minimum-Time-Problem) +- [Quantum State Sampling Problem](#Quantum-State-Sampling-Problem) + +Unitary Problem Templates: +- [Unitary Smooth Pulse Problem](#Unitary-Smooth-Pulse-Problem) +- [Unitary Minimum Time Problem](#Unitary-Minimum-Time-Problem) +- [Unitary Sampling Problem](#Unitary-Sampling-Problem) +=# + +# ----- +# ## Ket Problem Templates + +#= +#### Quantum State Smooth Pulse Problem +```@docs; canonical = false +QuantumStateSmoothPulseProblem +``` + +Each problem starts with a `QuantumSystem` object, which is used to define the system's +Hamiltonian and control operators. The goal is to find a control pulse that drives the +intial state, `ψ_init`, to a target state, `ψ_goal`. +=# + +# _define the quantum system_ +system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) +ψ_init = Vector{ComplexF64}([1.0, 0.0]) +ψ_goal = Vector{ComplexF64}([0.0, 1.0]) +T = 51 +Δt = 0.2 + +# _create the smooth pulse problem_ +state_prob = QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, T, Δt); + +# _check the fidelity before solving_ +println("Before: ", rollout_fidelity(state_prob.trajectory, system)) + +# _solve the problem_ +solve!(state_prob, max_iter=100, verbose=true, print_level=1); + +# _check the fidelity after solving_ +println("After: ", rollout_fidelity(state_prob.trajectory, system)) + +# _extract the control pulses_ +state_prob.trajectory.a |> size + +#= +#### Quantum State Minimum Time Problem +```@docs; canonical = false +QuantumStateMinimumTimeProblem +``` +=# + +# _create the minimum time problem_ +min_state_prob = QuantumStateMinimumTimeProblem(state_prob, ψ_goal); + +# _check the previous duration_ +println("Duration before: ", get_duration(state_prob.trajectory)) + +# _solve the minimum time problem_ +solve!(min_state_prob, max_iter=100, verbose=true, print_level=1); + +# _check the new duration_ +println("Duration after: ", get_duration(min_state_prob.trajectory)) + +# _the fidelity is preserved by a constraint_ +println("Fidelity after: ", rollout_fidelity(min_state_prob.trajectory, system)) + +#= +#### Quantum State Sampling Problem +```@docs; canonical = false +QuantumStateSamplingProblem +``` +=# + +# _create a sampling problem_ +driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) +sampling_state_prob = QuantumStateSamplingProblem([system, driftless_system], ψ_init, ψ_goal, T, Δt); + +# _new keys are added to the trajectory for the new states_ +println(sampling_state_prob.trajectory.state_names) + +# _solve the sampling problem for a few iterations_ +solve!(sampling_state_prob, max_iter=25, verbose=true, print_level=1); + +# _check the fidelity of the sampling problem (use the updated key to get the initial and goal)_ +println("After (original system): ", rollout_fidelity(sampling_state_prob.trajectory, system, state_name=:ψ̃1_system_1)) +println("After (new system): ", rollout_fidelity(sampling_state_prob.trajectory, driftless_system, state_name=:ψ̃1_system_1)) + +# _compare this to using the original problem on the new system_ +println("After (new system, original `prob`): ", rollout_fidelity(state_prob.trajectory, driftless_system)) + + +# ----- +# ## Unitary Problem Templates + +#= +#### Unitary Smooth Pulse Problem + +```@docs; canonical = false +UnitarySmoothPulseProblem +``` + +The `UnitarySmoothPulseProblem` is similar to the `QuantumStateSmoothPulseProblem`, but +instead of driving the system to a target state, the goal is to drive the system to a +target unitary operator, `U_goal`. + +=# + +system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) +U_goal = GATES.H +T = 51 +Δt = 0.2 + +prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt); + +# _check the fidelity before solving_ +println("Before: ", unitary_rollout_fidelity(prob.trajectory, system)) + +# _finding an optimal control is as simple as calling `solve!`_ +solve!(prob, max_iter=100, verbose=true, print_level=1); + +# _check the fidelity after solving_ +println("After: ", unitary_rollout_fidelity(prob.trajectory, system)) + +#= +The `NamedTrajectory` object stores the control pulse, state variables, and the time grid. +=# + +# _extract the control pulses_ +prob.trajectory.a |> size + +#= +#### Unitary Minimum Time Problem + +```@docs; canonical = false +UnitaryMinimumTimeProblem +``` + +The goal of this problem is to find the shortest time it takes to drive the system to a +target unitary operator, `U_goal`. The problem is solved by minimizing the sum of all of +the time steps. It is constructed from `prob` in the previous example. +=# + +min_prob = UnitaryMinimumTimeProblem(prob, U_goal); + +# _check the previous duration_ +println("Duration before: ", get_duration(prob.trajectory)) + +# _solve the minimum time problem_ +solve!(min_prob, max_iter=100, verbose=true, print_level=1); + +# _check the new duration_ +println("Duration after: ", get_duration(min_prob.trajectory)) + +# _the fidelity is preserved by a constraint_ +println("Fidelity after: ", unitary_rollout_fidelity(min_prob.trajectory, system)) + +#= +#### Unitary Sampling Problem + +```@docs; canonical = false +UnitarySamplingProblem +``` + +A sampling problem is used to solve over multiple quantum systems with the same control. +This can be useful for exploring robustness, for example. +=# + +# _create a sampling problem_ +driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) +sampling_prob = UnitarySamplingProblem([system, driftless_system], U_goal, T, Δt); + +# _new keys are addded to the trajectory for the new states_ +println(sampling_prob.trajectory.state_names) + +# _the `solve!` proceeds as in the [Quantum State Sampling Problem](#Quantum-State-Sampling-Problem)]_ -# ## Unitary Minimum Time Problem -# ```@docs; canonical = false -# UnitaryMinimumTimeProblem -# ``` diff --git a/docs/literate/man/rollouts.jl b/docs/literate/man/rollouts.jl deleted file mode 100644 index 87aacad6..00000000 --- a/docs/literate/man/rollouts.jl +++ /dev/null @@ -1,16 +0,0 @@ -# ```@meta -# CollapsedDocStrings = true -# ``` - -# # Rollouts - -using QuantumCollocation -using SparseArrays # for visualization - -#= - -Rollouts are a way to visualize the evolution of a quantum system. The various rollout -functions provided in this module allow for the validation of the solution to a quantum -optimal control problem. - -=# \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 54cde02f..1793984a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,12 +16,11 @@ open(joinpath(@__DIR__, "src", "index.md"), write = true) do io end end +# TODO: Callbacks are currently broken pages = [ "Home" => "index.md", "Manual" => [ "Problem Templates" => "generated/man/problem_templates.md", - "Rollouts" => "generated/man/rollouts.md", - "Callbacks" => "generated/man/ipopt_callbacks.md", ], "Library" => "lib.md", ] diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..83d0ef82 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,116 @@ +```@raw html +
+ + Piccolo.jl + +
+ +
+ + + + + + + +
+ Documentation +
+ + Stable + + + Dev + + + arXiv + + 		+ Build Status +
+ + Build Status + + + Coverage + + 		+ License +
+ + MIT License + + 		+ Support +
+ + Unitary Fund + +
+
+ +
+ Quickly set up and solve problem templates for quantum optimal control +
+
+``` + +# QuantumCollocation.jl + +**QuantumCollocation.jl** sets up and solves *quantum control problems* as nonlinear programs (NLPs). In this context, a generic quantum control problem looks like +```math +\begin{aligned} + \arg \min_{\mathbf{Z}}\quad & J(\mathbf{Z}) \\ + \nonumber \text{s.t.}\qquad & \mathbf{f}(\mathbf{Z}) = 0 \\ + \nonumber & \mathbf{g}(\mathbf{Z}) \le 0 +\end{aligned} +``` +where $\mathbf{Z}$ is a trajectory containing states and controls, from [NamedTrajectories.jl](https://github.com/harmoniqs/NamedTrajectories.jl). + +### Problem Templates + +*Problem Templates* are reusable design patterns for setting up and solving common quantum control problems. + +For example, a *UnitarySmoothPulseProblem* is tasked with generating a *pulse* sequence $a_{1:T-1}$ in orderd to minimize infidelity, subject to constraints from the Schroedinger equation, +```math + \begin{aligned} + \arg \min_{\mathbf{Z}}\quad & |1 - \mathcal{F}(U_T, U_\text{goal})| \\ + \nonumber \text{s.t.} + \qquad & U_{t+1} = \exp\{- i H(a_t) \Delta t_t \} U_t, \quad \forall\, t \\ + \end{aligned} +``` +while a *UnitaryMinimumTimeProblem* minimizes time and constrains fidelity, +```math + \begin{aligned} + \arg \min_{\mathbf{Z}}\quad & \sum_{t=1}^T \Delta t_t \\ + \qquad & U_{t+1} = \exp\{- i H(a_t) \Delta t_t \} U_t, \quad \forall\, t \\ + \nonumber & \mathcal{F}(U_T, U_\text{goal}) \ge 0.9999 + \end{aligned} +``` + +In each case, the dynamics between *knot points* $(U_t, a_t)$ and $(U_{t+1}, a_{t+1})$ are enforced as constraints on the states, which are free variables in the solver; this optimization framework is called *direct collocation*. For details of our implementation please see our award-winning IEEE QCE 2023 paper, [Direct Collocation for Quantum Optimal Control](https://arxiv.org/abs/2305.03261). If you use QuantumCollocation.jl in your work, please cite :raised_hands:! + +Problem templates give the user the ability to add other constraints and objective functions to this problem and solve it efficiently using [Ipopt.jl](https://github.com/jump-dev/Ipopt.jl) and [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) under the hood. + +## Installation + +This package is registered! To install, enter the Julia REPL, type `]` to enter pkg mode, and then run: +```julia +pkg> add QuantumCollocation +``` + +## Example + +### Single Qubit Hadamard Gate +```Julia +using QuantumCollocation + +T = 50 +Δt = 0.2 +system = QuantumSystem([PAULIS[:X], PAULIS[:Y]]) +U_goal = GATES.H + +# Hadamard Gate +prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt) +solve!(prob, max_iter=100) +``` diff --git a/docs/src/lib.md b/docs/src/lib.md index 835b8de1..11d9e133 100644 --- a/docs/src/lib.md +++ b/docs/src/lib.md @@ -10,16 +10,6 @@ Modules = [QuantumCollocation.ProblemTemplates] Modules = [QuantumCollocation.Options] ``` -## Rollouts -```@autodocs -Modules = [QuantumCollocation.Rollouts] -``` - -## Direct Sums -```@autodocs -Modules = [QuantumCollocation.DirectSums] -``` - ## Trajectory Initialization ```@autodocs Modules = [QuantumCollocation.TrajectoryInitialization] diff --git a/src/QuantumCollocation.jl b/src/QuantumCollocation.jl index 82168459..454e9d39 100644 --- a/src/QuantumCollocation.jl +++ b/src/QuantumCollocation.jl @@ -2,8 +2,10 @@ module QuantumCollocation using Reexport -@reexport using QuantumCollocationCore -@reexport using PiccoloQuantumObjects +@reexport using DirectTrajOpt + +include("piccolo_options.jl") +@reexport using .Options include("trajectory_initialization.jl") @reexport using .TrajectoryInitialization @@ -11,6 +13,15 @@ include("trajectory_initialization.jl") include("trajectory_interpolations.jl") @reexport using .TrajectoryInterpolations +include("quantum_objectives.jl") +@reexport using .QuantumObjectives + +include("quantum_constraints.jl") +@reexport using .QuantumConstraints + +include("quantum_integrators.jl") +@reexport using .QuantumIntegrators + include("problem_templates/_problem_templates.jl") @reexport using .ProblemTemplates diff --git a/src/piccolo_options.jl b/src/piccolo_options.jl new file mode 100644 index 00000000..d784f70d --- /dev/null +++ b/src/piccolo_options.jl @@ -0,0 +1,38 @@ +module Options + +export PiccoloOptions + +using ExponentialAction + +""" + PiccoloOptions + +Options for the Piccolo quantum optimal control library. + +# Fields +- `verbose::Bool = true`: Print verbose output +- `timesteps_all_equal::Bool = true`: Use equal timesteps +- `rollout_integrator::Function = expv`: Integrator to use for rollout +- `geodesic = true`: Use the geodesic to initialize the optimization. +- `zero_initial_and_final_derivative::Bool=false`: Zero the initial and final control pulse derivatives. +- `complex_control_norm_constraint_name::Union{Nothing, Symbol} = nothing`: Name of the complex control norm constraint. +- `complex_control_norm_constraint_radius::Float64 = 1.0`: Radius of the complex control norm constraint. +- `bound_state::Bool = false`: Bound the state. +- `leakage_suppression::Bool = false`: Suppress leakage. +- `R_leakage::Float64 = 1.0`: Leakage suppression parameter. +""" +@kwdef mutable struct PiccoloOptions + verbose::Bool = true + timesteps_all_equal::Bool = true + rollout_integrator::Function = expv + geodesic::Bool = true + zero_initial_and_final_derivative::Bool = true + complex_control_norm_constraint_name::Union{Nothing, Symbol} = nothing + complex_control_norm_constraint_radius::Float64 = 1.0 + bound_state::Bool = false + leakage_suppression::Bool = false + state_leakage_indices::AbstractVector{Int} = Int[] + R_leakage::Float64 = 1.0 +end + +end \ No newline at end of file diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 9beda1e5..fd6e512d 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -1,14 +1,14 @@ module ProblemTemplates -using ..DirectSums -using ..Rollouts using ..TrajectoryInitialization -using ..Losses +using ..QuantumObjectives +using ..QuantumConstraints +using ..QuantumIntegrators +using ..Options -using Distributions using TrajectoryIndexingUtils using NamedTrajectories -using QuantumCollocationCore +using DirectTrajOpt using PiccoloQuantumObjects using LinearAlgebra using SparseArrays @@ -18,16 +18,12 @@ using TestItems include("unitary_smooth_pulse_problem.jl") include("unitary_minimum_time_problem.jl") -include("unitary_robustness_problem.jl") -include("unitary_direct_sum_problem.jl") include("unitary_sampling_problem.jl") -include("unitary_bang_bang_problem.jl") include("quantum_state_smooth_pulse_problem.jl") include("quantum_state_minimum_time_problem.jl") include("quantum_state_sampling_problem.jl") -include("density_operator_smooth_pulse_problem.jl") function apply_piccolo_options!( J::Objective, @@ -37,37 +33,49 @@ function apply_piccolo_options!( state_names::AbstractVector{Symbol}, timestep_name::Symbol; state_leakage_indices::Union{Nothing, AbstractVector{<:AbstractVector{Int}}}=nothing, + free_time::Bool=true, ) if piccolo_options.leakage_suppression - if isnothing(state_leakage_indices) - error("You must provide leakage indices for leakage suppression.") - end - for (state_name, leakage_indices) ∈ zip(state_names, state_leakage_indices) - J += L1Regularizer!( - constraints, - state_name, - traj; - R_value=piccolo_options.R_leakage, - indices=leakage_indices, - eval_hessian=piccolo_options.eval_hessian - ) - end + throw(error("L1 is not implemented.")) + # if piccolo_options.verbose + # println("\tapplying leakage suppression: $(state_names)") + # end + + # if isnothing(state_leakage_indices) + # error("You must provide leakage indices for leakage suppression.") + # end + # for (state_name, leakage_indices) ∈ zip(state_names, state_leakage_indices) + # J += L1Regularizer!( + # constraints, + # state_name, + # traj; + # R_value=piccolo_options.R_leakage, + # indices=leakage_indices, + # ) + # end end - if piccolo_options.free_time + if free_time + if piccolo_options.verbose + println("\tapplying timesteps_all_equal constraint: $(traj.timestep)") + end if piccolo_options.timesteps_all_equal push!( constraints, - TimeStepsAllEqualConstraint(timestep_name, traj) + TimeStepsAllEqualConstraint(traj) ) end end if !isnothing(piccolo_options.complex_control_norm_constraint_name) - norm_con = ComplexModulusContraint( + if piccolo_options.verbose + println("\tapplying complex control norm constraint: $(piccolo_options.complex_control_norm_constraint_name)") + end + norm_con = NonlinearKnotPointConstraint( + a -> [norm(a)^2 - piccolo_options.complex_control_norm_constraint_radius^2], piccolo_options.complex_control_norm_constraint_name, - piccolo_options.complex_control_norm_constraint_radius, traj; + equality=false, ) push!(constraints, norm_con) end @@ -83,6 +91,7 @@ function apply_piccolo_options!( state_name::Symbol, timestep_name::Symbol; state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, + kwargs... ) state_names = [ name for name ∈ traj.names @@ -97,6 +106,7 @@ function apply_piccolo_options!( state_names, timestep_name; state_leakage_indices=isnothing(state_leakage_indices) ? nothing : fill(state_leakage_indices, length(state_names)), + kwargs... ) end diff --git a/src/problem_templates/density_operator_smooth_pulse_problem.jl b/src/problem_templates/density_operator_smooth_pulse_problem.jl deleted file mode 100644 index 90be0f05..00000000 --- a/src/problem_templates/density_operator_smooth_pulse_problem.jl +++ /dev/null @@ -1,123 +0,0 @@ -export DensityOperatorSmoothPulseProblem - -function DensityOperatorSmoothPulseProblem( - system::OpenQuantumSystem, - ρ_init::AbstractMatrix, - ψ_goal::AbstractVector, - T::Int, - Δt::Union{Float64, Vector{Float64}}; - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), - constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - init_trajectory::Union{NamedTrajectory, Nothing}=nothing, - a_bound::Float64=1.0, - a_bounds=fill(a_bound, system.n_drives), - a_guess::Union{Matrix{Float64}, Nothing}=nothing, - da_bound::Float64=Inf, - da_bounds=fill(da_bound, system.n_drives), - dda_bound::Float64=1.0, - dda_bounds=fill(dda_bound, system.n_drives), - Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), - Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), - drive_derivative_σ::Float64=0.01, - Q::Float64=100.0, - R=1e-2, - R_a::Union{Float64, Vector{Float64}}=R, - R_da::Union{Float64, Vector{Float64}}=R, - R_dda::Union{Float64, Vector{Float64}}=R, - leakage_suppression=false, - R_leakage=1e-1, - control_norm_constraint=false, - control_norm_constraint_components=nothing, - control_norm_R=nothing, - kwargs... -) - # Trajectory - if !isnothing(init_trajectory) - traj = init_trajectory - else - traj = initialize_trajectory( - ρ_init, - ψ_goal*ψ_goal', - T, - Δt, - system.n_drives, - (a_bounds, da_bounds, dda_bounds); - free_time=piccolo_options.free_time, - Δt_bounds=(Δt_min, Δt_max), - drive_derivative_σ=drive_derivative_σ, - a_guess=a_guess, - system=system, - ) - end - - # Objective - J = DensityOperatorPureStateInfidelityObjective(:ρ⃗̃, ψ_goal; Q=Q) - J += QuadraticRegularizer(:a, traj, R_a) - J += QuadraticRegularizer(:da, traj, R_da) - J += QuadraticRegularizer(:dda, traj, R_dda) - - # Constraints - if leakage_suppression - if operator isa EmbeddedOperator - leakage_indices = get_iso_vec_leakage_indices(operator) - J += L1Regularizer!( - constraints, :Ũ⃗, traj, - R_value=R_leakage, - indices=leakage_indices, - eval_hessian=piccolo_options.eval_hessian - ) - else - @warn "leakage_suppression is not supported for non-embedded operators, ignoring." - end - end - - if piccolo_options.free_time - if piccolo_options.timesteps_all_equal - push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) - end - end - - if control_norm_constraint - @assert !isnothing(control_norm_constraint_components) "control_norm_constraint_components must be provided" - @assert !isnothing(control_norm_R) "control_norm_R must be provided" - norm_con = ComplexModulusContraint( - :a, - control_norm_R, - traj; - name_comps=control_norm_constraint_components, - ) - push!(constraints, norm_con) - end - - # Integrators - # if piccolo_options.integrator == :pade - # unitary_integrator = - # UnitaryPadeIntegrator(system, :Ũ⃗, :a, traj; order=piccolo_options.pade_order) - # elseif piccolo_options.integrator == :exponential - # unitary_integrator = - # UnitaryExponentialIntegrator(system, :Ũ⃗, :a, traj) - # else - # error("integrator must be one of (:pade, :exponential)") - # end - - density_operator_integrator = DensityOperatorExponentialIntegrator( - :ρ⃗̃, :a, system, traj - ) - - integrators = [ - density_operator_integrator, - DerivativeIntegrator(:a, :da, traj), - DerivativeIntegrator(:da, :dda, traj), - ] - - return QuantumControlProblem( - traj, - J, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - kwargs... - ) -end diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index e3f74ff1..d0a6c8b2 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -5,124 +5,110 @@ export QuantumStateMinimumTimeProblem QuantumStateMinimumTimeProblem(traj, sys, obj, integrators, constraints; kwargs...) QuantumStateMinimumTimeProblem(prob; kwargs...) -Construct a `QuantumControlProblem` for the minimum time problem of reaching a target state. - -# Arguments -- `traj::NamedTrajectory`: The initial trajectory. -- `sys::QuantumSystem`: The quantum system. -- `obj::Objective`: The objective function. -- `integrators::Vector{<:AbstractIntegrator}`: The integrators. -- `constraints::Vector{<:AbstractConstraint}`: The constraints. -or -- `prob::QuantumControlProblem`: The quantum control problem. +Construct a `DirectTrajOptProblem` for the minimum time problem of reaching a target state. # Keyword Arguments - `state_name::Symbol=:ψ̃`: The symbol for the state variables. - `final_fidelity::Union{Real, Nothing}=nothing`: The final fidelity. - `D=1.0`: The cost weight on the time. -- `ipopt_options::IpoptOptions=IpoptOptions()`: The Ipopt options. - `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. -- `kwargs...`: Additional keyword arguments, passed to `QuantumControlProblem`. - """ function QuantumStateMinimumTimeProblem end function QuantumStateMinimumTimeProblem( - traj::NamedTrajectory, - obj::Objective, - integrators::Vector{<:AbstractIntegrator}, + trajectory::NamedTrajectory, + ψ_goals::AbstractVector{<:AbstractVector{<:ComplexF64}}, + objective::Objective, + dynamics::TrajectoryDynamics, constraints::Vector{<:AbstractConstraint}; state_name::Symbol=:ψ̃, - control_name::Symbol=:a, - final_fidelity::Union{Real, Nothing}=nothing, - D=1.0, - ipopt_options::IpoptOptions=IpoptOptions(), + final_fidelity::Float64=1.0, + D=100.0, piccolo_options::PiccoloOptions=PiccoloOptions(), - kwargs... ) - state_names = [name for name in traj.names if startswith(name, state_name)] + state_names = [name for name in trajectory.names if startswith(string(name), string(state_name))] @assert length(state_names) ≥ 1 "No matching states found in trajectory" + @assert length(state_names) == length(ψ_goals) "Number of states and goals must match" - obj += MinimumTimeObjective(traj; D=D, eval_hessian=piccolo_options.eval_hessian) - - # Default to average state fidelity - if isnothing(final_fidelity) - vals = [iso_fidelity(traj[n][:, end], traj.goal[n]) for n ∈ state_names] - final_fidelity = sum(vals) / length(vals) + if piccolo_options.verbose + println(" constructing QuantumStateMinimumTimeProblem...") + println("\tfinal fidelity: $(final_fidelity)") end - for state_name in state_names - fidelity_constraint = FinalQuantumStateFidelityConstraint( + objective += MinimumTimeObjective( + trajectory, D=D, timesteps_all_equal=piccolo_options.timesteps_all_equal + ) + + for (state_name, ψ_goal) in zip(state_names, ψ_goals) + fidelity_constraint = FinalKetFidelityConstraint( + ψ_goal, state_name, final_fidelity, - traj, - eval_hessian=piccolo_options.eval_hessian + trajectory, ) push!(constraints, fidelity_constraint) end - return QuantumControlProblem( - traj, - obj, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... + return DirectTrajOptProblem( + trajectory, + objective, + dynamics, + constraints ) end +# TODO: goals from trajectory? + function QuantumStateMinimumTimeProblem( - prob::QuantumControlProblem; - obj::Objective=get_objective(prob), - constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), - ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), - piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), - build_trajectory_constraints=false, + prob::DirectTrajOptProblem, + ψ_goals::AbstractVector{<:AbstractVector{<:ComplexF64}}; + objective::Objective=prob.objective, + constraints::AbstractVector{<:AbstractConstraint}=prob.constraints, kwargs... ) - piccolo_options.build_trajectory_constraints = build_trajectory_constraints - return QuantumStateMinimumTimeProblem( - copy(prob.trajectory), - obj, - prob.integrators, + prob.trajectory, + ψ_goals, + objective, + prob.dynamics, constraints; - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, kwargs... ) end +function QuantumStateMinimumTimeProblem( + prob::DirectTrajOptProblem, + ψ_goal::AbstractVector{<:ComplexF64}; + kwargs... +) + return QuantumStateMinimumTimeProblem(prob, [ψ_goal]; kwargs...) +end + # *************************************************************************** # @testitem "Test quantum state minimum time" begin - using NamedTrajectories - using PiccoloQuantumObjects - - # System - T = 50 - Δt = 0.2 - sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) - ψ_init = Vector{ComplexF64}[[1.0, 0.0]] - ψ_target = Vector{ComplexF64}[[0.0, 1.0]] - - prob = QuantumStateSmoothPulseProblem( - sys, ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions( - verbose=false, - integrator=:pade - ) - ) - initial = sum(get_timesteps(prob.trajectory)) - mintime_prob = QuantumStateMinimumTimeProblem(prob) - solve!(mintime_prob, max_iter=20) - final = sum(get_timesteps(mintime_prob.trajectory)) - @test final < initial - - # Test with final fidelity - QuantumStateMinimumTimeProblem(prob, final_fidelity=0.99) + using NamedTrajectories + using PiccoloQuantumObjects + + T = 51 + Δt = 0.2 + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + ψ_init = Vector{ComplexF64}[[1.0, 0.0]] + ψ_target = Vector{ComplexF64}[[0.0, 1.0]] + + prob = QuantumStateSmoothPulseProblem( + sys, ψ_init, ψ_target, T, Δt; + piccolo_options=PiccoloOptions(verbose=false) + ) + initial = sum(get_timesteps(prob.trajectory)) + + min_prob = QuantumStateMinimumTimeProblem( + prob, ψ_target, + piccolo_options=PiccoloOptions(verbose=false) + ) + solve!(min_prob, max_iter=50, verbose=false, print_level=1) + final = sum(get_timesteps(min_prob.trajectory)) + + @test final < initial end diff --git a/src/problem_templates/quantum_state_sampling_problem.jl b/src/problem_templates/quantum_state_sampling_problem.jl index 62e7f3f7..0859e3d1 100644 --- a/src/problem_templates/quantum_state_sampling_problem.jl +++ b/src/problem_templates/quantum_state_sampling_problem.jl @@ -1,5 +1,8 @@ export QuantumStateSamplingProblem +""" + +""" function QuantumStateSamplingProblem end function QuantumStateSamplingProblem( @@ -8,10 +11,9 @@ function QuantumStateSamplingProblem( ψ_goals::AbstractVector{<:AbstractVector{<:AbstractVector{<:ComplexF64}}}, T::Int, Δt::Union{Float64,Vector{Float64}}; + ket_integrator=KetIntegrator, system_weights=fill(1.0, length(systems)), init_trajectory::Union{NamedTrajectory,Nothing}=nothing, - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), state_name::Symbol=:ψ̃, control_name::Symbol=:a, timestep_name::Symbol=:Δt, @@ -29,12 +31,17 @@ function QuantumStateSamplingProblem( R_a::Union{Float64,Vector{Float64}}=R, R_da::Union{Float64,Vector{Float64}}=R, R_dda::Union{Float64,Vector{Float64}}=R, - leakage_indices::Union{Nothing, <:AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - kwargs... + piccolo_options::PiccoloOptions=PiccoloOptions(), ) @assert length(ψ_inits) == length(ψ_goals) + if piccolo_options.verbose + println(" constructing QuantumStateSamplingProblem...") + println("\tusing integrator: $(typeof(ket_integrator))") + println("\tusing $(length(ψ_inits)) initial state(s)") + end + # Outer dimension is the system, inner dimension is the initial state state_names = [ [Symbol(string(state_name) * "$(i)_system_$(j)") for i in eachindex(ψs)] @@ -56,21 +63,17 @@ function QuantumStateSamplingProblem( state_names=names, control_name=control_name, timestep_name=timestep_name, - free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), + zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, bound_state=piccolo_options.bound_state, a_guess=a_guess, system=sys, rollout_integrator=piccolo_options.rollout_integrator, - verbose=piccolo_options.verbose + verbose=false # loop ) end - traj = merge( - trajs, - merge_names=(; a=1, da=1, dda=1, Δt=1), - free_time=piccolo_options.free_time - ) + traj = merge(trajs, merge_names=(; a=1, da=1, dda=1, Δt=1),free_time=true) end control_names = [ @@ -79,57 +82,41 @@ function QuantumStateSamplingProblem( ] # Objective - J = QuadraticRegularizer(control_names[1], traj, R_a; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[2], traj, R_da; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[3], traj, R_dda; timestep_name=timestep_name) + J = QuadraticRegularizer(control_names[1], traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) for (weight, names) in zip(system_weights, state_names) for name in names - J += weight * QuantumStateObjective(name, traj, Q) + J += KetInfidelityLoss(name, traj, Q=weight * Q) end end # Optional Piccolo constraints and objectives apply_piccolo_options!( J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=leakage_indices + state_leakage_indices=piccolo_options.state_leakage_indices ) # Integrators state_integrators = [] - for (system, names) ∈ zip(systems, state_names) + for (sys, names) ∈ zip(systems, state_names) for name ∈ names - if piccolo_options.integrator == :pade - state_integrator = QuantumStatePadeIntegrator( - name, control_name, system, traj; - order=piccolo_options.pade_order - ) - elseif piccolo_options.integrator == :exponential - state_integrator = QuantumStateExponentialIntegrator( - name, control_name, system, traj - ) - else - error("integrator must be one of (:pade, :exponential)") - end - push!(state_integrators, state_integrator) + push!(state_integrators, ket_integrator(sys, traj, name, control_name)) end end integrators = [ state_integrators..., - DerivativeIntegrator(control_name, control_names[2], traj), - DerivativeIntegrator(control_names[2], control_names[3], traj), + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), ] - return QuantumControlProblem( + return DirectTrajOptProblem( traj, J, integrators; constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... ) end @@ -156,88 +143,78 @@ function QuantumStateSamplingProblem( args...; kwargs... ) + # Broadcast the initial and target states to all systems return QuantumStateSamplingProblem(systems, [ψ_init], [ψ_goal], args...; kwargs...) end # *************************************************************************** # @testitem "Sample systems with single initial, target" begin - # System + using PiccoloQuantumObjects + T = 50 Δt = 0.2 sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) - - # Single initial and target states - # -------------------------------- ψ_init = Vector{ComplexF64}([1.0, 0.0]) ψ_target = Vector{ComplexF64}([0.0, 1.0]) - + prob = QuantumStateSamplingProblem( [sys1, sys2], ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) - + state_name = :ψ̃ state_names = [n for n ∈ prob.trajectory.names if startswith(string(n), string(state_name))] - sys1_state_names = [n for n ∈ state_names if endswith(string(n), "1")] - + sys_state_names = [n for n ∈ state_names if endswith(string(n), "1")] + # Separately compute all unique initial and goal state fidelities - init = [rollout_fidelity(prob.trajectory, sys1; state_name=n) for n in sys1_state_names] - solve!(prob, max_iter=20) - final = [rollout_fidelity(prob.trajectory, sys1, state_name=n) for n in sys1_state_names] - @test all(final .> init) - - # Check that a_guess can be used - prob = QuantumStateSamplingProblem( - [sys1, sys2], ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - a_guess=prob.trajectory.a - ) - solve!(prob, max_iter=20) - - # Compare a (very smooth) solution without robustness - # ------------------------------------- - prob_default = QuantumStateSmoothPulseProblem( - sys2, ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - R_dda=1e2, - robustness=false - ) - solve!(prob_default, max_iter=20) - final_default = rollout_fidelity(prob_default.trajectory, sys1) - # Pick any initial state - final_robust = rollout_fidelity(prob.trajectory, sys1, state_name=state_names[1]) - @test final_robust > final_default + inits = [] + for sys in [sys1, sys2] + push!(inits, [rollout_fidelity(prob.trajectory, sys; state_name=n) for n in state_names]) + end + + solve!(prob, max_iter=20, print_level=1, verbose=false) + + for (init, sys) in zip(inits, [sys1, sys2]) + final = [rollout_fidelity(prob.trajectory, sys, state_name=n) for n in state_names] + @test all(final .> init) + end end @testitem "Sample systems with multiple initial, target" begin - # System + using PiccoloQuantumObjects + T = 50 Δt = 0.2 sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) - sys2 = QuantumSystem(0.0 * GATES[:Z], [GATES[:X], GATES[:Y]]) - + sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) + # Multiple initial and target states - # ---------------------------------- ψ_inits = Vector{ComplexF64}.([[1.0, 0.0], [0.0, 1.0]]) ψ_targets = Vector{ComplexF64}.([[0.0, 1.0], [1.0, 0.0]]) - + prob = QuantumStateSamplingProblem( [sys1, sys2], ψ_inits, ψ_targets, T, Δt; - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) - + state_name = :ψ̃ state_names = [n for n ∈ prob.trajectory.names if startswith(string(n), string(state_name))] - sys1_state_names = [n for n ∈ state_names if endswith(string(n), "1")] - - init = [rollout_fidelity(prob.trajectory, sys1, state_name=n) for n in sys1_state_names] - solve!(prob, max_iter=20) - final = [rollout_fidelity(prob.trajectory, sys1, state_name=n) for n in sys1_state_names] - @test all(final .> init) + sys_state_names = [n for n ∈ state_names if endswith(string(n), "1")] + + # Separately compute all unique initial and goal state fidelities + inits = [] + for sys in [sys1, sys2] + push!(inits, [rollout_fidelity(prob.trajectory, sys; state_name=n) for n in state_names]) + end + + solve!(prob, max_iter=20, print_level=1, verbose=false) + + for (init, sys) in zip(inits, [sys1, sys2]) + final = [rollout_fidelity(prob.trajectory, sys, state_name=n) for n in state_names] + @test all(final .> init) + end end + +# TODO: Test that a_guess can be used diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 836014f3..f5d9798a 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -6,9 +6,10 @@ export QuantumStateSmoothPulseProblem QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, T, Δt; kwargs...) QuantumStateSmoothPulseProblem(H_drift, H_drives, args...; kwargs...) -Create a quantum state smooth pulse problem. The goal is to find a control pulse -`a(t)` that drives all of the initial states `ψ_inits` to the corresponding -target states `ψ_goals` using `T` timesteps of size `Δt`. This problem also controls the first and second derivatives of the control pulse, `da(t)` and `dda(t)`, to ensure smoothness. +Create a quantum state smooth pulse problem. The goal is to find a control pulse `a(t)` +that drives all of the initial states `ψ_inits` to the corresponding target states +`ψ_goals` using `T` timesteps of size `Δt`. This problem also controls the first and +second derivatives of the control pulse, `da(t)` and `dda(t)`, to ensure smoothness. # Arguments - `system::AbstractQuantumSystem`: The quantum system. @@ -27,8 +28,6 @@ with # Keyword Arguments -- `ipopt_options::IpoptOptions=IpoptOptions()`: The IPOPT options. -- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. - `state_name::Symbol=:ψ̃`: The name of the state variable. - `control_name::Symbol=:a`: The name of the control variable. - `timestep_name::Symbol=:Δt`: The name of the timestep variable. @@ -38,7 +37,6 @@ with - `a_guess::Union{Matrix{Float64}, Nothing}=nothing`: The initial guess for the control pulse. - `da_bound::Float64=Inf`: The bound on the first derivative of the control pulse. - `da_bounds::Vector{Float64}=fill(da_bound, length(system.G_drives))`: The bounds on the first derivative of the control pulse. -- `zero_initial_and_final_derivative::Bool=false`: Whether to enforce zero initial and final derivative. - `dda_bound::Float64=1.0`: The bound on the second derivative of the control pulse. - `dda_bounds::Vector{Float64}=fill(dda_bound, length(system.G_drives))`: The bounds on the second derivative of the control pulse. - `Δt_min::Float64=0.5 * Δt`: The minimum timestep size. @@ -49,8 +47,8 @@ with - `R_a::Union{Float64, Vector{Float64}}=R`: The weight on the control pulse. - `R_da::Union{Float64, Vector{Float64}}=R`: The weight on the first derivative of the control pulse. - `R_dda::Union{Float64, Vector{Float64}}=R`: The weight on the second derivative of the control pulse. -- `leakage_operator::Union{Nothing, EmbeddedOperator}=nothing`: The leakage operator, if leakage suppression is desired. - `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: The constraints. +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. """ function QuantumStateSmoothPulseProblem end @@ -60,8 +58,7 @@ function QuantumStateSmoothPulseProblem( ψ_goals::Vector{<:AbstractVector{<:ComplexF64}}, T::Int, Δt::Union{Float64, <:AbstractVector{Float64}}; - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), + ket_integrator=KetIntegrator, state_name::Symbol=:ψ̃, control_name::Symbol=:a, timestep_name::Symbol=:Δt, @@ -71,23 +68,26 @@ function QuantumStateSmoothPulseProblem( a_guess::Union{Matrix{Float64}, Nothing}=nothing, da_bound::Float64=Inf, da_bounds::Vector{Float64}=fill(da_bound, sys.n_drives), - zero_initial_and_final_derivative::Bool=false, dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, sys.n_drives), Δt_min::Float64=0.001 * Δt, Δt_max::Float64=2.0 * Δt, - drive_derivative_σ::Float64=0.01, Q::Float64=100.0, R=1e-2, R_a::Union{Float64, Vector{Float64}}=R, R_da::Union{Float64, Vector{Float64}}=R, R_dda::Union{Float64, Vector{Float64}}=R, - leakage_indices::Union{Nothing, <:AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - kwargs... + piccolo_options::PiccoloOptions=PiccoloOptions(), ) @assert length(ψ_inits) == length(ψ_goals) + if piccolo_options.verbose + println(" constructing QuantumStateSmoothPulseProblem...") + println("\tusing integrator: $(typeof(ket_integrator))") + println("\tusing $(length(ψ_inits)) initial state(s)") + end + # Trajectory if !isnothing(init_trajectory) traj = init_trajectory @@ -102,11 +102,9 @@ function QuantumStateSmoothPulseProblem( state_name=state_name, control_name=control_name, timestep_name=timestep_name, - zero_initial_and_final_derivative=zero_initial_and_final_derivative, - free_time=piccolo_options.free_time, + zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, Δt_bounds=(Δt_min, Δt_max), bound_state=piccolo_options.bound_state, - drive_derivative_σ=drive_derivative_σ, a_guess=a_guess, system=sys, rollout_integrator=piccolo_options.rollout_integrator, @@ -125,85 +123,45 @@ function QuantumStateSmoothPulseProblem( ] # Objective - J = QuadraticRegularizer(control_names[1], traj, R_a; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[2], traj, R_da; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[3], traj, R_dda; timestep_name=timestep_name) + J = QuadraticRegularizer(control_names[1], traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) for name ∈ state_names - J += QuantumStateObjective(name, traj, Q) + J += KetInfidelityLoss(name, traj; Q=Q) end # Optional Piccolo constraints and objectives apply_piccolo_options!( J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=leakage_indices + state_leakage_indices=piccolo_options.state_leakage_indices ) - # Integrators + state_names = [ + name for name ∈ traj.names + if startswith(string(name), string(state_name)) + ] + state_integrators = [] - if length(ψ_inits) == 1 - if piccolo_options.integrator == :pade - state_integrators = [QuantumStatePadeIntegrator( - state_name, - control_name, - sys, - traj; - order=piccolo_options.pade_order - )] - elseif piccolo_options.integrator == :exponential - state_integrators = [QuantumStateExponentialIntegrator( - state_name, - control_name, - sys, - traj - )] - else - error("integrator must be one of (:pade, :exponential)") - end - else - state_names = [ - name for name ∈ traj.names - if startswith(string(name), string(state_name)) - ] - state_integrators = [] - for i = 1:length(ψ_inits) - if piccolo_options.integrator == :pade - state_integrator = QuantumStatePadeIntegrator( - state_names[i], - control_name, - sys, - traj; - order=piccolo_options.pade_order - ) - elseif piccolo_options.integrator == :exponential - state_integrator = QuantumStateExponentialIntegrator( - state_names[i], - control_name, - sys, - traj - ) - else - error("integrator must be one of (:pade, :exponential)") - end - push!(state_integrators, state_integrator) - end + + for name ∈ state_names + push!( + state_integrators, + ket_integrator(sys, traj, name, control_name) + ) end integrators = [ state_integrators..., - DerivativeIntegrator(control_name, control_names[2], traj), - DerivativeIntegrator(control_names[2], control_names[3], traj) + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]) ] - return QuantumControlProblem( + return DirectTrajOptProblem( traj, J, integrators; constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... ) end @@ -230,78 +188,29 @@ end # *************************************************************************** # @testitem "Test quantum state smooth pulse" begin - # System - T = 50 + using PiccoloQuantumObjects + + T = 51 Δt = 0.2 sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) ψ_init = Vector{ComplexF64}([1.0, 0.0]) ψ_target = Vector{ComplexF64}([0.0, 1.0]) - + # Single initial and target states # -------------------------------- prob = QuantumStateSmoothPulseProblem( sys, ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) initial = rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=50) + solve!(prob, max_iter=50, print_level=1, verbose=false) final = rollout_fidelity(prob.trajectory, sys) @test final > initial - - # Multiple initial and target states - # ---------------------------------- - ψ_inits = Vector{ComplexF64}.([[1.0, 0.0], [0.0, 1.0]]) - ψ_targets = Vector{ComplexF64}.([[0.0, 1.0], [1.0, 0.0]]) - prob = QuantumStateSmoothPulseProblem( - sys, ψ_inits, ψ_targets, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false) - ) - initial = rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=50) - final = rollout_fidelity(prob.trajectory, sys) - @test all(final .> initial) end -@testitem "Test quantum state smooth pulse w/ exponential integrator" begin - # System - T = 50 - Δt = 0.2 - sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) - ψ_init = Vector{ComplexF64}([1.0, 0.0]) - ψ_target = Vector{ComplexF64}([0.0, 1.0]) - integrator=:exponential +@testitem "Test multiple quantum states smooth pulse" begin + using PiccoloQuantumObjects - # Single initial and target states - # -------------------------------- - prob = QuantumStateSmoothPulseProblem( - sys, ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, integrator=integrator) - ) - initial = rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) - final = rollout_fidelity(prob.trajectory, sys) - @test final > initial - - # Multiple initial and target states - # ---------------------------------- - ψ_inits = Vector{ComplexF64}.([[1.0, 0.0], [0.0, 1.0]]) - ψ_targets = Vector{ComplexF64}.([[0.0, 1.0], [1.0, 0.0]]) - prob = QuantumStateSmoothPulseProblem( - sys, ψ_inits, ψ_targets, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, integrator=integrator) - ) - initial = rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) - final = rollout_fidelity(prob.trajectory, sys) - @test all(final .> initial) -end - -@testitem "Test quantum state with multiple initial states and final states" begin - # System T = 50 Δt = 0.2 sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) @@ -310,14 +219,11 @@ end prob = QuantumStateSmoothPulseProblem( sys, ψ_inits, ψ_targets, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - state_name=:psi, - control_name=:u, - timestep_name=:dt + piccolo_options=PiccoloOptions(verbose=false) ) initial = rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) + solve!(prob, max_iter=50, print_level=1, verbose=false) final = rollout_fidelity(prob.trajectory, sys) + final, initial @test all(final .> initial) end diff --git a/src/problem_templates/unitary_bang_bang_problem.jl b/src/problem_templates/unitary_bang_bang_problem.jl deleted file mode 100644 index 9669180f..00000000 --- a/src/problem_templates/unitary_bang_bang_problem.jl +++ /dev/null @@ -1,265 +0,0 @@ -export UnitaryBangBangProblem - - -@doc raw""" - UnitaryBangBangProblem(system::QuantumSystem, operator, T, Δt; kwargs...) - UnitaryBangBangProblem(H_drift, H_drives, operator, T, Δt; kwargs...) - -Construct a `QuantumControlProblem` for a free-time unitary gate problem with bang-bang control pulses. - -```math -\begin{aligned} -\underset{\vec{\tilde{U}}, a, \dot{a}, \Delta t}{\text{minimize}} & \quad -Q \cdot \ell\qty(\vec{\tilde{U}}_T, \vec{\tilde{U}}_{\text{goal}}) + R_{\text{bang-bang}} \cdot \sum_t |\dot{a}_t| \\ -\text{ subject to } & \quad \vb{P}^{(n)}\qty(\vec{\tilde{U}}_{t+1}, \vec{\tilde{U}}_t, a_t, \Delta t_t) = 0 \\ -& \quad a_{t+1} - a_t - \dot{a}_t \Delta t_t = 0 \\ -& \quad |a_t| \leq a_{\text{bound}} \\ -& \quad |\dot{a}_t| \leq da_{\text{bound}} \\ -& \quad \Delta t_{\text{min}} \leq \Delta t_t \leq \Delta t_{\text{max}} \\ -\end{aligned} -``` - -where, for $U \in SU(N)$, - -```math -\ell\qty(\vec{\tilde{U}}_T, \vec{\tilde{U}}_{\text{goal}}) = -\abs{1 - \frac{1}{N} \abs{ \tr \qty(U_{\text{goal}}, U_T)} } -``` - -is the *infidelity* objective function, $Q$ is a weight, $R_a$, and $R_{\dot{a}}$ are weights on the regularization terms, and $\vb{P}^{(n)}$ is the $n$th-order Pade integrator. - -TODO: Document bang-bang modification. - -# Arguments - -- `H_drift::AbstractMatrix{<:Number}`: the drift hamiltonian -- `H_drives::Vector{<:AbstractMatrix{<:Number}}`: the control hamiltonians -or -- `system::QuantumSystem`: the system to be controlled -with -- `operator::AbstractPiccoloOperator`: the target unitary, either in the form of an `EmbeddedOperator` or a `Matrix{ComplexF64} -- `T::Int`: the number of timesteps -- `Δt::Float64`: the (initial) time step size - -# Keyword Arguments -- `ipopt_options::IpoptOptions=IpoptOptions()`: the options for the Ipopt solver -- `piccolo_options::PiccoloOptions=PiccoloOptions()`: the options for the Piccolo solver -- `state_name::Symbol = :Ũ⃗`: the name of the state variable -- `control_name::Symbol = :a`: the name of the control variable -- `timestep_name::Symbol = :Δt`: the name of the timestep variable -- `init_trajectory::Union{NamedTrajectory, Nothing}=nothing`: an initial trajectory to use -- `a_bound::Float64=1.0`: the bound on the control pulse -- `a_bounds=fill(a_bound, length(system.G_drives))`: the bounds on the control pulses, one for each drive -- `a_guess::Union{Matrix{Float64}, Nothing}=nothing`: an initial guess for the control pulses -- `da_bound::Float64=1.0`: the bound on the control pulse derivative -- `da_bounds=fill(da_bound, length(system.G_drives))`: the bounds on the control pulse derivatives, one for each drive -- `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt)`: the minimum time step size -- `Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt)`: the maximum time step size -- `drive_derivative_σ::Float64=0.01`: the standard deviation of the initial guess for the control pulse derivatives -- `Q::Float64=100.0`: the weight on the infidelity objective -- `R=1e-2`: the weight on the regularization terms -- `quadratic_control_regularization=false`: whether or not to use quadratic regularization for the control pulses -- `R_a::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulses -- `R_da::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse derivatives -- `R_bang_bang::Union{Float64, Vector{Float64}}=1e-1`: the weight on the bang-bang regularization term -- `phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing`: the phase operators for free phase corrections -- `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: the constraints to enforce -""" -function UnitaryBangBangProblem end - -function UnitaryBangBangProblem( - system::AbstractQuantumSystem, - operator::AbstractPiccoloOperator, - T::Int, - Δt::Union{Float64, Vector{Float64}}; - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), - state_name::Symbol = :Ũ⃗, - control_name::Symbol = :a, - timestep_name::Symbol = :Δt, - init_trajectory::Union{NamedTrajectory, Nothing}=nothing, - a_bound::Float64=1.0, - a_bounds=fill(a_bound, system.n_drives), - a_guess::Union{Matrix{Float64}, Nothing}=nothing, - da_bound::Float64=1.0, - da_bounds=fill(da_bound, system.n_drives), - Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), - Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), - Q::Float64=100.0, - R=1e-2, - R_a::Union{Float64, Vector{Float64}}=R, - R_da::Union{Float64, Vector{Float64}}=0.0, - R_bang_bang::Union{Float64, Vector{Float64}}=1e-1, - phase_name::Symbol=:ϕ, - phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing, - constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - kwargs... -) - # Trajectory - if !isnothing(init_trajectory) - traj = init_trajectory - else - traj = initialize_trajectory( - operator, - T, - Δt, - system.n_drives, - (a_bounds, da_bounds); - state_name=state_name, - control_name=control_name, - timestep_name=timestep_name, - free_time=piccolo_options.free_time, - Δt_bounds=(Δt_min, Δt_max), - geodesic=piccolo_options.geodesic, - bound_state=piccolo_options.bound_state, - a_guess=a_guess, - system=system, - rollout_integrator=piccolo_options.rollout_integrator, - phase_name=phase_name, - phase_operators=phase_operators - ) - end - - subspace = operator isa EmbeddedOperator ? operator.subspace_indices : nothing - - # Objective - if isnothing(phase_operators) - J = UnitaryInfidelityObjective(state_name, traj, Q; - subspace=subspace, - ) - else - J = UnitaryFreePhaseInfidelityObjective( - state_name, phase_name, phase_operators, traj, Q; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian, - ) - end - - control_names = [ - name for name ∈ traj.names - if endswith(string(name), string(control_name)) - ] - - J += QuadraticRegularizer(control_names[1], traj, R_a; timestep_name=timestep_name) - # Default to R_da = 0.0 because L1 is turned on by this problem. - J += QuadraticRegularizer(control_names[2], traj, R_da; timestep_name=timestep_name) - - # Constraints - if R_bang_bang isa Float64 - R_bang_bang = fill(R_bang_bang, system.n_drives) - end - J += L1Regularizer!( - constraints, control_names[2], traj, - R=R_bang_bang, eval_hessian=piccolo_options.eval_hessian - ) - - # Optional Piccolo constraints and objectives - apply_piccolo_options!(J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=operator isa EmbeddedOperator ? get_leakage_indices(operator) : nothing - ) - - # Integrators - if piccolo_options.integrator == :pade - unitary_integrator = - UnitaryPadeIntegrator(state_name, control_names[1], system, traj; order=piccolo_options.pade_order) - elseif piccolo_options.integrator == :exponential - unitary_integrator = - UnitaryExponentialIntegrator(state_name, control_names[1], system, traj) - else - error("integrator must be one of (:pade, :exponential)") - end - - integrators = [ - unitary_integrator, - DerivativeIntegrator(control_names[1], control_names[2], traj), - ] - - return QuantumControlProblem( - traj, - J, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... - ) -end - -function UnitaryBangBangProblem( - H_drift::AbstractMatrix{<:Number}, - H_drives::Vector{<:AbstractMatrix{<:Number}}, - args...; - kwargs... -) - system = QuantumSystem(H_drift, H_drives) - return UnitaryBangBangProblem(system, args...; kwargs...) -end - -# *************************************************************************** # - -@testitem "Bang-bang Hadamard gate" begin - sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:X], GATES[:Y]]) - U_goal = GATES[:H] - T = 51 - Δt = 0.2 - - ipopt_options = IpoptOptions(print_level=1, max_iter=25) - piccolo_options = PiccoloOptions(verbose=false, pade_order=12) - - prob = UnitaryBangBangProblem( - sys, U_goal, T, Δt, - R_bang_bang=10., - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=:u - ) - - smooth_prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt; - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=:u - ) - initial = unitary_rollout_fidelity(prob.trajectory, sys; drive_name=:u) - solve!(prob) - final = unitary_rollout_fidelity(prob.trajectory, sys; drive_name=:u) - @test final > initial - solve!(smooth_prob) - threshold = 1e-3 - a_sparse = sum(prob.trajectory.du .> 5e-2) - a_dense = sum(smooth_prob.trajectory.du .> 5e-2) - @test a_sparse < a_dense -end - -@testitem "Free phase Y gate using X" begin - using Random - Random.seed!(1234) - phase_name = :ϕ - phase_operators = [PAULIS[:Z]] - - sys = QuantumSystem([PAULIS[:X]]) - - prob = UnitaryBangBangProblem( - sys, GATES[:Y], 51, 0.2; - phase_operators=phase_operators, - phase_name=phase_name, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, free_time=false) - ) - - before = prob.trajectory.global_data[phase_name] - solve!(prob, max_iter=20) - after = prob.trajectory.global_data[phase_name] - - @test before ≠ after - - @test unitary_rollout_fidelity( - prob.trajectory, - sys; - phases=prob.trajectory.global_data[phase_name], - phase_operators=phase_operators - ) > 0.9 - - @test unitary_rollout_fidelity(prob.trajectory, sys) < 0.9 -end diff --git a/src/problem_templates/unitary_direct_sum_problem.jl b/src/problem_templates/unitary_direct_sum_problem.jl deleted file mode 100644 index 1bf5b7bb..00000000 --- a/src/problem_templates/unitary_direct_sum_problem.jl +++ /dev/null @@ -1,273 +0,0 @@ -export UnitaryDirectSumProblem - - -@doc """ - UnitaryDirectSumProblem(probs, final_fidelity; kwargs...) - -Construct a `QuantumControlProblem` as a direct sum of unitary gate problems. The -purpose is to find solutions that are as close as possible with respect to one of -their components. In particular, this is useful for finding interpolatable control solutions. - -A graph of edges (specified by problem labels) will enforce a `PairwiseQuadraticRegularizer` between -the component trajectories of the problem in `probs` corresponding to the names of the edge in `edges` -with corresponding edge weight `Q`. - -Boundary values can be included to enforce a `QuadraticRegularizer` on edges where one of the nodes is -not optimized. The boundary values are specified as a dictionary with keys corresponding to the edge -labels and values corresponding to the boundary values. - -The default behavior is to use a 1D chain for the graph, i.e., enforce a `PairwiseQuadraticRegularizer` -between each neighbor of the provided `probs`. - -# Arguments -- `probs::AbstractVector{<:QuantumControlProblem}`: the problems to combine -- `final_fidelity::Real`: the fidelity to enforce between the component final unitaries and the component goal unitaries - -# Keyword Arguments -- `prob_labels::AbstractVector{<:String}}`: the labels for the problems -- graph::Union{Nothing, AbstractVector{<:Tuple{String, String}}, AbstractVector{<:Tuple{Symbol, Symbol}}}`: the graph of edges to enforce -- `boundary_values::Union{Nothing, AbstractDict{<:String, <:AbstractArray}, AbstractDict{<:Symbol, <:AbstractVector}}=nothing`: the boundary values for the problems -- `Q::Union{Float64, Vector{Float64}}=100.0`: the weights on the pairwise regularizers -- `Q_symb::Symbol=:Ũ⃗`: the symbol to use for the regularizer -- `R::Float64=1e-2`: the shared weight on all control terms (:a, :da, :dda is assumed) -- `R_a::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulses -- `R_da::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse derivatives -- `R_dda::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse second derivatives -- `R_b::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the boundary values -- `drive_derivative_σ::Float64=0.01`: the standard deviation of the initial guess for the control pulse derivatives -- `drive_reset_ratio::Float64=0.1`: amount of random noise to add to the control data (can help avoid hitting restoration if provided problems are converged) -- `fidelity_cost::Bool=false`: whether or not to include a fidelity cost in the objective -- `subspace::Union{AbstractVector{<:Integer}, Nothing}=nothing`: the subspace to use for the fidelity of each problem -- `ipopt_options::IpoptOptions=IpoptOptions()`: the options for the Ipopt solver -- `piccolo_options::PiccoloOptions=PiccoloOptions()`: the options for the Piccolo solver - -""" -function UnitaryDirectSumProblem( - probs::AbstractVector{<:QuantumControlProblem}, - final_fidelity::Real; - prob_labels::AbstractVector{<:String}=[string(i) for i ∈ 1:length(probs)], - graph::Union{Nothing, AbstractVector{<:Tuple{String, String}}, AbstractVector{<:Tuple{Symbol, Symbol}}}=nothing, - boundary_values::Union{AbstractDict{<:String, <:AbstractArray}, AbstractDict{<:Symbol, <:AbstractArray}}=Dict{String, Array}(), - control_name::Symbol=:a, - Q::Union{Float64, Vector{Float64}}=100.0, - Q_symb::Symbol=:dda, - R::Float64=1e-2, - R_a::Union{Float64, Vector{Float64}}=R, - R_da::Union{Float64, Vector{Float64}}=R, - R_dda::Union{Float64, Vector{Float64}}=R, - R_b::Union{Float64, Vector{Float64}}=R, - drive_derivative_σ::Float64=0.01, - drive_reset_ratio::Float64=0.50, - fidelity_cost::Bool=false, - subspace::Union{AbstractVector{<:Integer}, Nothing}=nothing, - ipopt_options::IpoptOptions=deepcopy(probs[1].ipopt_options), - piccolo_options::PiccoloOptions=deepcopy(probs[1].piccolo_options), - kwargs... -) - N = length(probs) - @assert length(prob_labels) == N "prob_labels must match length of probs" - @assert N ≥ 2 "At least two problems are required" - @assert 0 ≤ drive_reset_ratio ≤ 1 "drive_reset_ratio must be in [0, 1]" - @assert isempty(intersect(keys(boundary_values), prob_labels)) "Boundary value keys cannot be in prob_labels" - @assert all([:dda ∈ p.trajectory.names for p in probs]) "Only smooth pulse problems are supported." - n_derivatives = 2 - - # Default chain graph and boundary - boundary = Tuple{Symbol, Array}[] - if isnothing(graph) - graph = [ - (add_suffix(Q_symb, i), add_suffix(Q_symb, j)) - for (i, j) ∈ zip(prob_labels[1:N-1], prob_labels[2:N]) - ] - else - # Check that String edge labels are in prob_labels or boundary, and make Symbols - if eltype(eltype(graph)) == String - graph_symbols = Tuple{Symbol, Symbol}[] - for edge in graph - if edge[1] in prob_labels && edge[2] in prob_labels - push!(graph_symbols, (add_suffix(Q_symb, edge[1]), add_suffix(Q_symb, edge[2]))) - elseif edge[1] in keys(boundary_values) && edge[2] in prob_labels - push!(boundary, (add_suffix(Q_symb, edge[2]), boundary_values[edge[1]])) - elseif edge[1] in prob_labels && edge[2] in keys(boundary_values) - push!(boundary, (add_suffix(Q_symb, edge[1]), boundary_values[edge[2]])) - else - throw(ArgumentError("Edge labels must be in prob_labels or boundary_values")) - end - end - graph = graph_symbols - end - end - - # Build the direct sum system - - # merge suffix trajectories - traj = merge([add_suffix(p.trajectory, ℓ) for (p, ℓ) ∈ zip(probs, prob_labels)]) - - # add noise to control data (avoid restoration) - if drive_reset_ratio > 0 - for ℓ in prob_labels - a_symb, da_symb, dda_symb = add_suffix(:a, ℓ), add_suffix(:da, ℓ), add_suffix(:dda, ℓ) - n_drives = length(traj.components[a_symb]) - a, da, dda = TrajectoryInitialization.initialize_control_trajectory( - n_drives, - n_derivatives, - traj.T, - traj.bounds[a_symb], - drive_derivative_σ - ) - update!(traj, a_symb, (1 - drive_reset_ratio) * traj[a_symb] + drive_reset_ratio * a) - update!(traj, da_symb, (1 - drive_reset_ratio) * traj[da_symb] + drive_reset_ratio * da) - update!(traj, dda_symb, (1 - drive_reset_ratio) * traj[dda_symb] + drive_reset_ratio * dda) - end - end - - # Rebuild integrators - integrators = vcat([ - add_suffix(p.integrators, ℓ, p.trajectory, traj) - for (p, ℓ) ∈ zip(probs, prob_labels) - ]...) - - # Rebuild trajectory constraints - constraints = AbstractConstraint[] - - # Add goal constraints for each problem - for (p, ℓ) ∈ zip(probs, prob_labels) - goal_symb, = keys(p.trajectory.goal) - fidelity_constraint = FinalUnitaryFidelityConstraint( - add_suffix(goal_symb, ℓ), - final_fidelity, - traj, - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - push!(constraints, fidelity_constraint) - end - - # Build the objective function - J = PairwiseQuadraticRegularizer(traj, Q, graph) - - for (symb, val) ∈ boundary - J += QuadraticRegularizer(symb, traj, R_b; baseline=val) - end - - for ℓ ∈ prob_labels - J += QuadraticRegularizer(add_suffix(:a, ℓ), traj, R_a) - J += QuadraticRegularizer(add_suffix(:da, ℓ), traj, R_da) - J += QuadraticRegularizer(add_suffix(:dda, ℓ), traj, R_dda) - end - - # Add fidelity cost - if fidelity_cost - for ℓ ∈ prob_labels - Q_fid = isa(Q, Number) ? Q : Q[1] - J += UnitaryInfidelityObjective( - add_suffix(:Ũ⃗, ℓ), traj, Q_fid, - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - end - end - - return QuantumControlProblem( - traj, - J, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... - ) -end - - -# *************************************************************************** # - -@testitem "Construct direct sum problem" begin - using QuantumCollocationCore - using PiccoloQuantumObjects - sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:X], GATES[:Y]]) - U_goal1 = GATES[:X] - U_ε = haar_identity(2, 0.33) - U_goal2 = U_ε'GATES[:X]*U_ε - T = 50 - Δt = 0.2 - ops = IpoptOptions(print_level=1) - pops = PiccoloOptions(verbose=false, free_time=false) - - prob1 = UnitarySmoothPulseProblem(sys, U_goal1, T, Δt, ipopt_options=ops, piccolo_options=pops) - prob2 = UnitarySmoothPulseProblem(sys, U_goal2, T, Δt, ipopt_options=ops, piccolo_options=pops) - - # Test default - direct_sum_prob1 = UnitaryDirectSumProblem([prob1, prob2], 0.99) - state_names = vcat( - add_suffix(prob1.trajectory.state_names, "1")..., - add_suffix(prob2.trajectory.state_names, "2")... - ) - control_names = vcat( - add_suffix(prob1.trajectory.control_names, "1")..., - add_suffix(prob2.trajectory.control_names, "2")... - ) - @test issetequal(direct_sum_prob1.trajectory.state_names, state_names) - @test issetequal(direct_sum_prob1.trajectory.control_names, control_names) - - # Test label graph - direct_sum_prob2 = UnitaryDirectSumProblem( - [prob1, prob2], - 0.99, - prob_labels=["a", "b"], - graph=[("a", "b")], - verbose=false, - ipopt_options=ops) - state_names_ab = vcat( - add_suffix(prob1.trajectory.state_names, "a")..., - add_suffix(prob2.trajectory.state_names, "b")... - ) - control_names_ab = vcat( - add_suffix(prob1.trajectory.control_names, "a")..., - add_suffix(prob2.trajectory.control_names, "b")... - ) - @test issetequal(direct_sum_prob2.trajectory.state_names, state_names_ab) - @test issetequal(direct_sum_prob2.trajectory.control_names, control_names_ab) - - # Test bad graph - @test_throws ArgumentError UnitaryDirectSumProblem( - [prob1, prob2], - 0.99, - prob_labels=["a", "b"], - graph=[("x", "b")], - verbose=false, - ipopt_options=ops - ) - - # Test symbol graph - direct_sum_prob3 = UnitaryDirectSumProblem( - [prob1, prob2], - 0.99, - graph=[(:a1, :a2)], - verbose=false, - ipopt_options=ops - ) - @test issetequal(direct_sum_prob3.trajectory.state_names, state_names) - @test issetequal(direct_sum_prob3.trajectory.control_names, control_names) - - # Test triple - direct_sum_prob4 = UnitaryDirectSumProblem([prob1, prob2, prob1], 0.99, ipopt_options=ops) - state_names_triple = vcat(state_names..., add_suffix(prob1.trajectory.state_names, "3")...) - control_names_triple = vcat(control_names..., add_suffix(prob1.trajectory.control_names, "3")...) - @test issetequal(direct_sum_prob4.trajectory.state_names, state_names_triple) - @test issetequal(direct_sum_prob4.trajectory.control_names, control_names_triple) - - # Test boundary values - direct_sum_prob5 = UnitaryDirectSumProblem( - [prob1, prob2], - 0.99, - graph=[("x", "1"), ("1", "2")], - R_b=1e3, - Q_symb=:dda, - boundary_values=Dict("x"=>copy(prob1.trajectory[:dda])), - verbose=false, - ipopt_options=ops - ) - # # TODO: Check for objectives? -end diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index bcebdf22..49a3f82d 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -3,16 +3,17 @@ export UnitaryMinimumTimeProblem @doc raw""" UnitaryMinimumTimeProblem( + goal::AbstractPiccoloOperator, trajectory::NamedTrajectory, - system::AbstractQuantumSystem, objective::Objective, - integrators::Vector{<:AbstractIntegrator}, - constraints::Vector{<:AbstractConstraint}; + dynamics::TrajectoryDynamics, + constraints::AbstractVector{<:AbstractConstraint}; kwargs... ) UnitaryMinimumTimeProblem( - prob::QuantumControlProblem; + goal::AbstractPiccoloOperator, + prob::DirectTrajOptProblem; kwargs... ) @@ -28,109 +29,64 @@ J(\vec{\tilde{U}}, a, \dot{a}, \ddot{a}) + D \sum_t \Delta t_t \\ \end{aligned} ``` -# Arguments -- `trajectory::NamedTrajectory`: The initial trajectory. -- `system::AbstractQuantumSystem`: The quantum system. -- `objective::Objective`: The objective function (additional to the minimum-time objective). -- `integrators::Vector{<:AbstractIntegrator}`: The integrators. -- `constraints::Vector{<:AbstractConstraint}`: The constraints. - # Keyword Arguments -- `unitary_name::Symbol=:Ũ⃗`: The symbol for the unitary control. -- `final_fidelity::Float64=0.99`: The final fidelity. -- `D=1.0`: The weight for the minimum-time objective. -- `ipopt_options::IpoptOptions=IpoptOptions()`: The options for the Ipopt solver. -- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The options for the Piccolo solver. -- `kwargs...`: Additional keyword arguments to pass to `QuantumControlProblem`. +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. +- `unitary_name::Symbol=:Ũ⃗`: The name of the unitary for the goal. +- `final_fidelity::Float64=1.0`: The final fidelity constraint. +- `D::Float64=1.0`: The scaling factor for the minimum-time objective. """ function UnitaryMinimumTimeProblem end function UnitaryMinimumTimeProblem( trajectory::NamedTrajectory, - system::AbstractQuantumSystem, + goal::AbstractPiccoloOperator, objective::Objective, - integrators::Vector{<:AbstractIntegrator}, - constraints::Vector{<:AbstractConstraint}; + dynamics::TrajectoryDynamics, + constraints::AbstractVector{<:AbstractConstraint}; unitary_name::Symbol=:Ũ⃗, - control_name::Symbol=:a, - final_fidelity::Union{Real, Nothing}=nothing, - D=1.0, - ipopt_options::IpoptOptions=IpoptOptions(), + final_fidelity::Float64=1.0, + D::Float64=100.0, piccolo_options::PiccoloOptions=PiccoloOptions(), - phase_name::Symbol=:ϕ, - phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing, - subspace=nothing, - kwargs... ) - @assert unitary_name ∈ trajectory.names + if piccolo_options.verbose + println(" constructing UnitaryMinimumTimeProblem...") + println("\tfinal fidelity: $(final_fidelity)") + end objective += MinimumTimeObjective( - trajectory; D=D, eval_hessian=piccolo_options.eval_hessian + trajectory; D=D, timesteps_all_equal=piccolo_options.timesteps_all_equal ) - U_T = trajectory[unitary_name][:, end] - U_G = trajectory.goal[unitary_name] - subspace = isnothing(subspace) ? axes(iso_vec_to_operator(U_T), 1) : subspace - - if isnothing(phase_operators) - if isnothing(final_fidelity) - final_fidelity = iso_vec_unitary_fidelity(U_T, U_G, subspace=subspace) - end - - fidelity_constraint = FinalUnitaryFidelityConstraint( - unitary_name, final_fidelity, trajectory; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - else - if isnothing(final_fidelity) - phases = trajectory.global_data[phase_name] - final_fidelity = iso_vec_unitary_free_phase_fidelity( - U_T, U_G, phases, phase_operators; subspace=subspace - ) - end - - fidelity_constraint = FinalUnitaryFreePhaseFidelityConstraint( - unitary_name, phase_name, phase_operators, final_fidelity, trajectory; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - end + fidelity_constraint = FinalUnitaryFidelityConstraint( + goal, unitary_name, final_fidelity, trajectory + ) constraints = push!(constraints, fidelity_constraint) - return QuantumControlProblem( + return DirectTrajOptProblem( trajectory, objective, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... + dynamics, + constraints ) end +# TODO: how to handle apply_piccolo_options? +# TODO: goal from trajectory? + function UnitaryMinimumTimeProblem( - prob::QuantumControlProblem, - sys::AbstractQuantumSystem; - objective::Objective=get_objective(prob), - constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), - ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), - piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), - build_trajectory_constraints=false, + prob::DirectTrajOptProblem, + goal::AbstractPiccoloOperator; + objective::Objective=prob.objective, + constraints::AbstractVector{<:AbstractConstraint}=prob.constraints, kwargs... ) - piccolo_options.build_trajectory_constraints = build_trajectory_constraints - return UnitaryMinimumTimeProblem( - copy(prob.trajectory), - sys, + prob.trajectory, + goal, objective, - prob.integrators, + prob.dynamics, constraints; - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, kwargs... ) end @@ -139,6 +95,7 @@ end @testitem "Minimum time Hadamard gate" begin using NamedTrajectories + using PiccoloQuantumObjects H_drift = PAULIS[:Z] H_drives = [PAULIS[:X], PAULIS[:Y]] @@ -150,62 +107,31 @@ end prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) before = unitary_rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=50) + solve!(prob; max_iter=100, verbose=false, print_level=1) after = unitary_rollout_fidelity(prob.trajectory, sys) @test after > before - # Soft fidelity constraint - final_fidelity = minimum([0.99, after]) - mintime_prob = UnitaryMinimumTimeProblem(prob, sys, final_fidelity=final_fidelity) - solve!(mintime_prob; max_iter=100) - - # Test fidelity is approximatley staying above the constraint - @test unitary_rollout_fidelity(mintime_prob.trajectory, sys) ≥ (final_fidelity - 0.1 * final_fidelity) - duration_after = sum(get_timesteps(mintime_prob.trajectory)) - duration_before = sum(get_timesteps(prob.trajectory)) - @test duration_after < duration_before - - # Set up without a final fidelity - @test UnitaryMinimumTimeProblem(prob, sys) isa QuantumControlProblem - -end - -@testitem "Minimum time free phase" begin - using NamedTrajectories - - phase_operators = [PAULIS[:Z]] - sys = QuantumSystem([PAULIS[:X]]) - prob = UnitarySmoothPulseProblem( - sys, GATES[:H], 51, 0.2, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - phase_operators=phase_operators - ) - - # Soft fidelity constraint - final_fidelity = minimum([0.99, unitary_rollout_fidelity(prob.trajectory, sys)]) - mintime_prob = UnitaryMinimumTimeProblem( - prob, - sys; - final_fidelity=final_fidelity, - phase_operators=phase_operators + # soft fidelity constraint + min_prob = UnitaryMinimumTimeProblem( + prob, U_goal, + piccolo_options=PiccoloOptions(verbose=false) ) - solve!(mintime_prob; max_iter=100) + solve!(min_prob; max_iter=100, verbose=false, print_level=1) - duration_after = sum(get_timesteps(mintime_prob.trajectory)) + # test fidelity has stayed above the constraint + constraint_tol = 0.95 + final_fidelity = minimum([0.99, after]) + @test unitary_rollout_fidelity(min_prob.trajectory, sys) ≥ constraint_tol * final_fidelity + duration_after = sum(get_timesteps(min_prob.trajectory)) duration_before = sum(get_timesteps(prob.trajectory)) @test duration_after < duration_before +end - # Quick check for using default fidelity - @test UnitaryMinimumTimeProblem( - prob, - sys; - final_fidelity=final_fidelity, - phase_operators=phase_operators - ) isa QuantumControlProblem +@testitem "Test relaxed final_fidelity constraint" begin + final_fidelity = 0.95 + @test_broken false end diff --git a/src/problem_templates/unitary_robustness_problem.jl b/src/problem_templates/unitary_robustness_problem.jl deleted file mode 100644 index 4baea9c1..00000000 --- a/src/problem_templates/unitary_robustness_problem.jl +++ /dev/null @@ -1,264 +0,0 @@ -export UnitaryRobustnessProblem - - -@doc raw""" - UnitaryRobustnessProblem( - H_error, - trajectory, - system, - objective, - integrators, - constraints; - kwargs... - ) - - UnitaryRobustnessProblem(Hₑ, prob::QuantumControlProblem; kwargs...) - -Create a quantum control problem for robustness optimization of a unitary trajectory. - -# Keyword Arguments -- `unitary_name::Symbol=:Ũ⃗`: The symbol for the unitary trajectory in `trajectory`. -- `final_fidelity::Union{Real, Nothing}=nothing`: The target fidelity for the final unitary. -- `ipopt_options::IpoptOptions=IpoptOptions()`: Options for the Ipopt solver. -- `piccolo_options::PiccoloOptions=PiccoloOptions()`: Options for the Piccolo solver. -- `kwargs...`: Additional keyword arguments passed to `QuantumControlProblem`. -""" -function UnitaryRobustnessProblem end - -function UnitaryRobustnessProblem( - H_error::AbstractPiccoloOperator, - trajectory::NamedTrajectory, - system::QuantumSystem, - objective::Objective, - integrators::Vector{<:AbstractIntegrator}, - constraints::Vector{<:AbstractConstraint}; - unitary_name::Symbol=:Ũ⃗, - control_name::Symbol=:a, - final_fidelity::Union{Real, Nothing}=nothing, - phase_name::Symbol=:ϕ, - phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing, - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), - subspace=nothing, - kwargs... -) - @assert unitary_name ∈ trajectory.names - - objective += UnitaryRobustnessObjective( - H_error=H_error, - eval_hessian=piccolo_options.eval_hessian - ) - - U_T = trajectory[unitary_name][:, end] - U_G = trajectory.goal[unitary_name] - subspace = isnothing(subspace) ? axes(iso_vec_to_operator(U_T), 1) : subspace - - if isnothing(phase_operators) - if isnothing(final_fidelity) - final_fidelity = iso_vec_unitary_fidelity(U_T, U_G, subspace=subspace) - end - - fidelity_constraint = FinalUnitaryFidelityConstraint( - unitary_name, - final_fidelity, - trajectory; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - else - if isnothing(final_fidelity) - phases = trajectory.global_data[phase_name] - final_fidelity = iso_vec_unitary_free_phase_fidelity( - U_T, U_G, phases, phase_operators; subspace=subspace - ) - end - - fidelity_constraint = FinalUnitaryFreePhaseFidelityConstraint( - unitary_name, - phase_name, - phase_operators, - final_fidelity, - trajectory; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian - ) - end - - push!(constraints, fidelity_constraint) - - return QuantumControlProblem( - trajectory, - objective, - integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... - ) -end - -function UnitaryRobustnessProblem( - H_error::AbstractPiccoloOperator, - prob::QuantumControlProblem, - sys::AbstractQuantumSystem; - objective::Objective=get_objective(prob), - constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), - ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), - piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), - build_trajectory_constraints=false, - kwargs... -) - piccolo_options.build_trajectory_constraints = build_trajectory_constraints - - return UnitaryRobustnessProblem( - H_error, - copy(prob.trajectory), - sys, - objective, - prob.integrators, - constraints; - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - kwargs... - ) -end - -# *************************************************************************** # - -@testitem "Robust and Subspace Templates" begin - H_drift = zeros(3, 3) - H_drives = [create(3) + annihilate(3), im * (create(3) - annihilate(3))] - sys = QuantumSystem(H_drift, H_drives) - - U_goal = EmbeddedOperator(:X, sys) - H_embed = EmbeddedOperator(:Z, sys) - T = 51 - Δt = 0.2 - - # test initial problem - # --------------------- - prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, - R=1e-12, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false) - ) - before = unitary_rollout_fidelity(prob.trajectory, sys; subspace=U_goal.subspace) - solve!(prob, max_iter=15) - after = unitary_rollout_fidelity(prob.trajectory, sys; subspace=U_goal.subspace) - - # Subspace gate success - @test after > before - - - # set up without a final fidelity - # ------------------------------- - @test UnitaryRobustnessProblem(H_embed, prob, sys) isa QuantumControlProblem - - - # test robustness from previous problem - # -------------------------------------- - final_fidelity = 0.99 - rob_prob = UnitaryRobustnessProblem( - H_embed, prob, sys; - final_fidelity=final_fidelity, - ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=100.0, print_level=1), - ) - solve!(rob_prob, max_iter=50) - - loss(Z⃗) = UnitaryRobustnessObjective(H_error=H_embed).L(Z⃗, prob.trajectory) - - # Robustness improvement over default (or small initial) - # TODO: Can this test be improved? (might fail if unlucky) - after = loss(rob_prob.trajectory.datavec) - before = loss(prob.trajectory.datavec) - @test (after < before) || (before < 0.25) - - # TODO: Fidelity constraint approximately satisfied - @test_skip isapprox(unitary_rollout_fidelity(rob_prob; subspace=U_goal.subspace), 0.99, atol=0.05) -end - -@testitem "Set up a free phase problem" begin - using LinearAlgebra - δ1 = δ2 = -0.1 - T = 75 - Δt = 1.0 - n_levels = 3 - a = annihilate(n_levels) - id = I(n_levels) - a1 = kron(a, id) - a2 = kron(id, a) - H_drift = δ1 / 2 * a1' * a1' * a1 * a1 + δ2 / 2 * a2' * a2' * a2 * a2 - H_drives = [a1'a1, a2'a2, a1'a2 + a1*a2', im * (a1'a2 - a1 * a2')] - system = QuantumSystem(H_drift, H_drives) - U_goal = EmbeddedOperator( - GATES[:CZ], - get_subspace_indices([1:2, 1:2], [n_levels, n_levels]), - [n_levels, n_levels] - ) - - phase_operators = [PAULIS[:Z], PAULIS[:Z]] - prob = UnitarySmoothPulseProblem( - system, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, eval_hessian=false), - phase_operators=phase_operators, - ) - - ZZ = EmbeddedOperator( - reduce(kron, phase_operators), - get_subspace_indices([1:2, 1:2], [n_levels, n_levels]), - [n_levels, n_levels] - ) - - @test UnitaryRobustnessProblem( - ZZ, - prob, - system, - phase_operators=phase_operators, - subspace=U_goal.subspace, - ) isa QuantumControlProblem -end - -@testitem "Set up a free phase problem" begin - using LinearAlgebra - δ1 = δ2 = -0.1 - T = 75 - Δt = 1.0 - n_levels = 3 - a = annihilate(n_levels) - id = I(n_levels) - a1 = kron(a, id) - a2 = kron(id, a) - H_drift = δ1 / 2 * a1' * a1' * a1 * a1 + δ2 / 2 * a2' * a2' * a2 * a2 - H_drives = [a1'a1, a2'a2, a1'a2 + a1*a2', im * (a1'a2 - a1 * a2')] - system = QuantumSystem(H_drift, H_drives) - U_goal = EmbeddedOperator( - GATES[:CZ], - get_subspace_indices([1:2, 1:2], [n_levels, n_levels]), - [n_levels, n_levels] - ) - - phase_operators = [PAULIS[:Z], PAULIS[:Z]] - prob = UnitarySmoothPulseProblem( - system, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, eval_hessian=false), - phase_operators=phase_operators, - ) - - ZZ = EmbeddedOperator( - reduce(kron, phase_operators), - get_subspace_indices([1:2, 1:2], [n_levels, n_levels]), - [n_levels, n_levels] - ) - - @test UnitaryRobustnessProblem( - ZZ, - prob, - system; - phase_operators=phase_operators, - subspace=U_goal.subspace, - ) isa QuantumControlProblem -end diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index d6e41ae9..7e334079 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -4,13 +4,15 @@ export UnitarySamplingProblem @doc raw""" UnitarySamplingProblem(systemns, operator, T, Δt; kwargs...) -A `UnitarySamplingProblem` is a quantum control problem where the goal is to find a control pulse that generates a target unitary operator for a set of quantum systems. -The controls are shared among all systems, and the optimization seeks to find the control pulse that achieves the operator for each system. The idea is to enforce a +A `UnitarySamplingProblem` is a quantum control problem where the goal is to find a control +pulse that generates a target unitary operator for a set of quantum systems. +The controls are shared among all systems, and the optimization seeks to find the control +pulse that achieves the operator for each system. The idea is to enforce a robust solution by including multiple systems reflecting the problem uncertainty. # Arguments - `systems::AbstractVector{<:AbstractQuantumSystem}`: A vector of quantum systems. -- `operator::AbstractPiccoloOperator`: The target unitary operator. +- `operators::AbstractVector{<:AbstractPiccoloOperator}`: A vector of target operators. - `T::Int`: The number of time steps. - `Δt::Union{Float64, Vector{Float64}}`: The time step value or vector of time steps. @@ -18,8 +20,6 @@ robust solution by including multiple systems reflecting the problem uncertainty - `system_labels::Vector{String} = string.(1:length(systems))`: The labels for each system. - `system_weights::Vector{Float64} = fill(1.0, length(systems))`: The weights for each system. - `init_trajectory::Union{NamedTrajectory, Nothing} = nothing`: The initial trajectory. -- `ipopt_options::IpoptOptions = IpoptOptions()`: The IPOPT options. -- `piccolo_options::PiccoloOptions = PiccoloOptions()`: The Piccolo options. - `state_name::Symbol = :Ũ⃗`: The name of the state variable. - `control_name::Symbol = :a`: The name of the control variable. - `timestep_name::Symbol = :Δt`: The name of the timestep variable. @@ -38,7 +38,7 @@ robust solution by including multiple systems reflecting the problem uncertainty - `R_a::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control amplitudes. - `R_da::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control first derivatives. - `R_dda::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control second derivatives. -- `kwargs...`: Additional keyword arguments. +- `piccolo_options::PiccoloOptions = PiccoloOptions()`: The Piccolo options. """ function UnitarySamplingProblem( @@ -46,9 +46,9 @@ function UnitarySamplingProblem( operators::AbstractVector{<:AbstractPiccoloOperator}, T::Int, Δt::Union{Float64,Vector{Float64}}; + unitary_integrator=UnitaryIntegrator, system_weights=fill(1.0, length(systems)), init_trajectory::Union{NamedTrajectory,Nothing}=nothing, - ipopt_options::IpoptOptions=IpoptOptions(), piccolo_options::PiccoloOptions=PiccoloOptions(), state_name::Symbol=:Ũ⃗, control_name::Symbol=:a, @@ -71,6 +71,12 @@ function UnitarySamplingProblem( kwargs... ) @assert length(systems) == length(operators) + + if piccolo_options.verbose + println(" constructing UnitarySamplingProblem...") + println("\tusing integrator: $(typeof(unitary_integrator))") + println("\tusing $(length(systems)) systems") + end # Trajectory state_names = [ @@ -80,31 +86,27 @@ function UnitarySamplingProblem( if !isnothing(init_trajectory) traj = init_trajectory else - trajs = map(zip(systems, operators, state_names)) do (sys, op, s) + trajs = map(zip(systems, operators, state_names)) do (sys, op, st) initialize_trajectory( op, T, Δt, sys.n_drives, (a_bounds, da_bounds, dda_bounds); - state_name=s, + state_name=st, control_name=control_name, timestep_name=timestep_name, - free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), geodesic=piccolo_options.geodesic, bound_state=piccolo_options.bound_state, a_guess=a_guess, system=sys, rollout_integrator=piccolo_options.rollout_integrator, + verbose=false # loop ) end - traj = merge( - trajs, - merge_names=(; a=1, da=1, dda=1, Δt=1), - free_time=piccolo_options.free_time - ) + traj = merge(trajs, merge_names=(; a=1, da=1, dda=1, Δt=1), free_time=true) end control_names = [ @@ -113,15 +115,12 @@ function UnitarySamplingProblem( ] # Objective - J = QuadraticRegularizer(control_name, traj, R_a; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[2], traj, R_da; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[3], traj, R_dda; timestep_name=timestep_name) + J = QuadraticRegularizer(control_name, traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) for (weight, op, name) in zip(system_weights, operators, state_names) - J += weight * UnitaryInfidelityObjective( - name, traj, Q; - subspace=op isa EmbeddedOperator ? op.subspace_indices : nothing - ) + J += UnitaryInfidelityLoss(op, name, traj; Q=weight * Q) end # Optional Piccolo constraints and objectives @@ -131,38 +130,17 @@ function UnitarySamplingProblem( ) # Integrators - unitary_integrators = AbstractIntegrator[] - for (sys, Ũ⃗_name) in zip(systems, state_names) - if piccolo_options.integrator == :pade - push!( - unitary_integrators, - UnitaryPadeIntegrator(Ũ⃗_name, control_name, sys, traj; order=piccolo_options.pade_order) - ) - elseif piccolo_options.integrator == :exponential - push!( - unitary_integrators, - UnitaryExponentialIntegrator(Ũ⃗_name, control_name, sys, traj) - ) - else - error("integrator must be one of (:pade, :exponential)") - end - end - integrators = [ - unitary_integrators..., - DerivativeIntegrator(control_name, control_names[2], traj), - DerivativeIntegrator(control_names[2], control_names[3], traj), + [unitary_integrator(sys, traj, n, control_name) for (sys, n) in zip(systems, state_names)]..., + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), ] - return QuantumControlProblem( + return DirectTrajOptProblem( traj, J, integrators; constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... ) end @@ -183,83 +161,35 @@ function UnitarySamplingProblem( ) end -function UnitarySamplingProblem( - system::Function, - distribution::Sampleable, - num_samples::Int, - args...; - kwargs... -) - samples = rand(distribution, num_samples) - systems = [system(x) for x in samples] - return UnitarySamplingProblem( - systems, - args...; - kwargs... - ) -end - # *************************************************************************** # @testitem "Sample robustness test" begin - using Distributions - using Random - Random.seed!(1234) + using PiccoloQuantumObjects - n_samples = 5 T = 50 Δt = 0.2 timesteps = fill(Δt, T) operator = GATES[:H] systems(ζ) = QuantumSystem(ζ * GATES[:Z], [GATES[:X], GATES[:Y]]) - - ip_ops = IpoptOptions(print_level=1, recalc_y="yes", recalc_y_feas_tol=1e1) - pi_ops = PiccoloOptions(verbose=false) - + + samples = [0.0, 0.1] prob = UnitarySamplingProblem( - systems, Normal(0, 0.05), n_samples, operator, T, Δt, - ipopt_options=ip_ops, piccolo_options=pi_ops + [systems(x) for x in samples], operator, T, Δt, + piccolo_options=PiccoloOptions(verbose=false) ) - solve!(prob, max_iter=20) - - d_prob = UnitarySmoothPulseProblem( - systems(0), operator, T, Δt, - ipopt_options=ip_ops, piccolo_options=pi_ops + solve!(prob, max_iter=100, print_level=1, verbose=false) + + base_prob = UnitarySmoothPulseProblem( + systems(samples[1]), operator, T, Δt, + piccolo_options=PiccoloOptions(verbose=false) ) - solve!(prob, max_iter=20) - - # Check that the solution improves over the default - # ------------------------------------------------- - ζ_tests = -0.05:0.01:0.05 - Ũ⃗_goal = operator_to_iso_vec(operator) - fids = [] - default_fids = [] - for ζ in ζ_tests - Ũ⃗_end = unitary_rollout(prob.trajectory.a, timesteps, systems(ζ))[:, end] - push!(fids, iso_vec_unitary_fidelity(Ũ⃗_end, Ũ⃗_goal)) - - d_Ũ⃗_end = unitary_rollout(d_prob.trajectory.a, timesteps, systems(ζ))[:, end] - push!(default_fids, iso_vec_unitary_fidelity(d_Ũ⃗_end, Ũ⃗_goal)) + solve!(base_prob, max_iter=100, verbose=false, print_level=1) + + fid = [] + base_fid = [] + for x in range(samples[1], samples[1], length=5) + push!(fid, unitary_rollout_fidelity(prob.trajectory, systems(0.1), unitary_name=:Ũ⃗_system_1)) + push!(base_fid, unitary_rollout_fidelity(base_prob.trajectory, systems(0.1))) end - @test sum(fids) > sum(default_fids) - - # Check initial guess initialization - # ---------------------------------- - a_guess = prob.trajectory.a - - g1_prob = UnitarySamplingProblem( - [systems(0), systems(0)], operator, T, Δt, - a_guess=a_guess, - piccolo_options=pi_ops - ) - - @test g1_prob.trajectory.Ũ⃗_system_1 ≈ g1_prob.trajectory.Ũ⃗_system_2 - - g2_prob = UnitarySamplingProblem( - [systems(0), systems(0.05)], operator, T, Δt; - a_guess=a_guess, - piccolo_options=pi_ops - ) - - @test ~(g2_prob.trajectory.Ũ⃗_system_1 ≈ g2_prob.trajectory.Ũ⃗_system_2) + @test sum(fid) > sum(base_fid) end diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 4cee0579..107d0de4 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -5,7 +5,7 @@ export UnitarySmoothPulseProblem UnitarySmoothPulseProblem(system::AbstractQuantumSystem, operator, T, Δt; kwargs...) UnitarySmoothPulseProblem(H_drift, H_drives, operator, T, Δt; kwargs...) -Construct a `QuantumControlProblem` for a free-time unitary gate problem with smooth control pulses enforced by constraining the second derivative of the pulse trajectory, i.e., +Construct a `DirectTrajOptProblem` for a free-time unitary gate problem with smooth control pulses enforced by constraining the second derivative of the pulse trajectory, i.e., ```math \begin{aligned} @@ -36,23 +36,21 @@ or - `H_drift::AbstractMatrix{<:Number}`: the drift hamiltonian - `H_drives::Vector{<:AbstractMatrix{<:Number}}`: the control hamiltonians with -- `operator::AbstractPiccoloOperator`: the target unitary, either in the form of an `EmbeddedOperator` or a `Matrix{ComplexF64} +- `goal::AbstractPiccoloOperator`: the target unitary, either in the form of an `EmbeddedOperator` or a `Matrix{ComplexF64} - `T::Int`: the number of timesteps - `Δt::Float64`: the (initial) time step size # Keyword Arguments -- `ipopt_options::IpoptOptions=IpoptOptions()`: the options for the Ipopt solver - `piccolo_options::PiccoloOptions=PiccoloOptions()`: the options for the Piccolo solver - `state_name::Symbol = :Ũ⃗`: the name of the state - `control_name::Symbol = :a`: the name of the control - `timestep_name::Symbol = :Δt`: the name of the timestep - `init_trajectory::Union{NamedTrajectory, Nothing}=nothing`: an initial trajectory to use +- `a_guess::Union{Matrix{Float64}, Nothing}=nothing`: an initial guess for the control pulses - `a_bound::Float64=1.0`: the bound on the control pulse - `a_bounds::Vector{Float64}=fill(a_bound, length(system.G_drives))`: the bounds on the control pulses, one for each drive -- `a_guess::Union{Matrix{Float64}, Nothing}=nothing`: an initial guess for the control pulses - `da_bound::Float64=Inf`: the bound on the control pulse derivative - `da_bounds::Vector{Float64}=fill(da_bound, length(system.G_drives))`: the bounds on the control pulse derivatives, one for each drive -- `zero_initial_and_final_derivative::Bool=false`: whether to enforce zero initial and final control pulse derivatives - `dda_bound::Float64=1.0`: the bound on the control pulse second derivative - `dda_bounds::Vector{Float64}=fill(dda_bound, length(system.G_drives))`: the bounds on the control pulse second derivatives, one for each drive - `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt)`: the minimum time step size @@ -62,28 +60,26 @@ with - `R_a::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulses - `R_da::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse derivatives - `R_dda::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse second derivatives -- `phase_name::Symbol=:ϕ`: the name of the phase -- `phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing`: the phase operators for free phase corrections - `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: the constraints to enforce """ +function UnitarySmoothPulseProblem end + function UnitarySmoothPulseProblem( system::AbstractQuantumSystem, - operator::AbstractPiccoloOperator, + goal::AbstractPiccoloOperator, T::Int, Δt::Union{Float64, <:AbstractVector{Float64}}; - ipopt_options::IpoptOptions=IpoptOptions(), - piccolo_options::PiccoloOptions=PiccoloOptions(), + unitary_integrator=UnitaryIntegrator, state_name::Symbol = :Ũ⃗, control_name::Symbol = :a, timestep_name::Symbol = :Δt, init_trajectory::Union{NamedTrajectory, Nothing}=nothing, + a_guess::Union{Matrix{Float64}, Nothing}=nothing, a_bound::Float64=1.0, a_bounds=fill(a_bound, system.n_drives), - a_guess::Union{Matrix{Float64}, Nothing}=nothing, da_bound::Float64=Inf, da_bounds::Vector{Float64}=fill(da_bound, system.n_drives), - zero_initial_and_final_derivative::Bool=false, dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, system.n_drives), Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), @@ -93,18 +89,20 @@ function UnitarySmoothPulseProblem( R_a::Union{Float64, Vector{Float64}}=R, R_da::Union{Float64, Vector{Float64}}=R, R_dda::Union{Float64, Vector{Float64}}=R, - phase_name::Symbol=:ϕ, - phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - kwargs... + piccolo_options::PiccoloOptions=PiccoloOptions(), ) + if piccolo_options.verbose + println(" constructing UnitarySmoothPulseProblem...") + println("\tusing integrator: $(typeof(unitary_integrator))") + end # Trajectory if !isnothing(init_trajectory) traj = init_trajectory else traj = initialize_trajectory( - operator, + goal, T, Δt, system.n_drives, @@ -112,81 +110,46 @@ function UnitarySmoothPulseProblem( state_name=state_name, control_name=control_name, timestep_name=timestep_name, - free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), - zero_initial_and_final_derivative=zero_initial_and_final_derivative, + zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, geodesic=piccolo_options.geodesic, bound_state=piccolo_options.bound_state, a_guess=a_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, - phase_name=phase_name, - phase_operators=phase_operators + verbose=piccolo_options.verbose ) end - # Subspace - subspace = operator isa EmbeddedOperator ? operator.subspace : nothing - # Objective - if isnothing(phase_operators) - J = UnitaryInfidelityObjective( - state_name, traj, Q; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian, - ) - else - J = UnitaryFreePhaseInfidelityObjective( - state_name, phase_name, phase_operators, traj, Q; - subspace=subspace, - eval_hessian=piccolo_options.eval_hessian, - ) - end + J = UnitaryInfidelityLoss(goal, state_name, traj; Q=Q) control_names = [ name for name ∈ traj.names if endswith(string(name), string(control_name)) ] - J += QuadraticRegularizer(control_names[1], traj, R_a; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[2], traj, R_da; timestep_name=timestep_name) - J += QuadraticRegularizer(control_names[3], traj, R_dda; timestep_name=timestep_name) + J += QuadraticRegularizer(control_names[1], traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) # Optional Piccolo constraints and objectives apply_piccolo_options!( J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=operator isa EmbeddedOperator ? get_leakage_indices(operator) : nothing + state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing ) - # Integrators - - if piccolo_options.integrator == :pade - unitary_integrator = - UnitaryPadeIntegrator(state_name, control_name, system, traj; - order=piccolo_options.pade_order - ) - elseif piccolo_options.integrator == :exponential - unitary_integrator = - UnitaryExponentialIntegrator(state_name, control_name, system, traj) - else - error("integrator must be one of (:pade, :exponential)") - end - integrators = [ - unitary_integrator, - DerivativeIntegrator(control_name, control_names[2], traj), - DerivativeIntegrator(control_names[2], control_names[3], traj), + unitary_integrator(system, traj, state_name, control_name), + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), ] - return QuantumControlProblem( + return DirectTrajOptProblem( traj, J, integrators; - constraints=constraints, - ipopt_options=ipopt_options, - piccolo_options=piccolo_options, - control_name=control_name, - kwargs... + constraints=constraints ) end @@ -203,6 +166,8 @@ end # *************************************************************************** # @testitem "Hadamard gate" begin + using PiccoloQuantumObjects + sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) U_goal = GATES[:H] T = 51 @@ -211,163 +176,57 @@ end prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt; da_bound=1.0, - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) initial = unitary_rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) + solve!(prob, max_iter=100, verbose=false, print_level=1) final = unitary_rollout_fidelity(prob.trajectory, sys) @test final > initial end -@testitem "Hadamard gate with exponential integrator" begin - sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) - U_goal = GATES[:H] - T = 51 - Δt = 0.2 - - prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, integrator=:exponential) - ) - - initial = unitary_rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) - final = unitary_rollout_fidelity(prob.trajectory, sys) - @test final > initial -end +@testitem "Hadamard gate with bound states and control norm constraint" begin + using PiccoloQuantumObjects -@testitem "Hadamard gate with exponential integrator, bounded states, and control norm constraint" begin sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) U_goal = GATES[:H] T = 51 Δt = 0.2 - piccolo_options = PiccoloOptions( - verbose=false, - integrator=:exponential, - # jacobian_structure=false, - bound_state=true, - complex_control_norm_constraint_name=:a - ) - prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=piccolo_options + piccolo_options=PiccoloOptions( + verbose=false, + bound_state=true, + complex_control_norm_constraint_name=:a + ) ) initial = unitary_rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=20) + solve!(prob, max_iter=100, verbose=false, print_level=1) final = unitary_rollout_fidelity(prob.trajectory, sys) - @test final > initial + @test_broken false + # @test final > initial end - - @testitem "EmbeddedOperator Hadamard gate" begin + using PiccoloQuantumObjects + a = annihilate(3) sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) U_goal = EmbeddedOperator(GATES[:H], sys) T = 51 Δt = 0.2 - # Test embedded operator - # ---------------------- - prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false) - ) - - initial = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) - solve!(prob, max_iter=20) - final = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) - @test final > initial - - # Test leakage suppression - # ------------------------ - a = annihilate(4) - sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) - U_goal = EmbeddedOperator(GATES[:H], sys) - T = 50 - Δt = 0.2 - prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - leakage_suppression=true, R_leakage=1e-1, - ipopt_options=IpoptOptions(print_level=1), piccolo_options=PiccoloOptions(verbose=false) ) initial = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) - solve!(prob, max_iter=20) + solve!(prob, max_iter=100, verbose=false, print_level=1) final = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) @test final > initial end -@testitem "Additional Objective" begin - H_drift = GATES[:Z] - H_drives = [GATES[:X], GATES[:Y]] - U_goal = GATES[:H] - T = 50 - Δt = 0.2 - - prob_vanilla = UnitarySmoothPulseProblem( - H_drift, H_drives, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - ) - - J_extra = QuadraticSmoothnessRegularizer(:dda, prob_vanilla.trajectory, 10.0) - - prob_additional = UnitarySmoothPulseProblem( - H_drift, H_drives, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false), - additional_objective=J_extra, - ) - - J_prob_vanilla = Problems.get_objective(prob_vanilla) - - J_additional = Problems.get_objective(prob_additional) - - Z = prob_vanilla.trajectory - Z⃗ = vec(prob_vanilla.trajectory) - - @test J_prob_vanilla.L(Z⃗, Z) + J_extra.L(Z⃗, Z) ≈ J_additional.L(Z⃗, Z) -end -@testitem "Free phase Y gate using X" begin - using Random - # Random.seed!(1234) - phase_name = :ϕ - phase_operators = [PAULIS[:Z]] - sys = QuantumSystem([PAULIS[:X]]) - prob = UnitarySmoothPulseProblem( - sys, - GATES[:Y], - 51, - 0.2; - phase_operators=phase_operators, - phase_name=phase_name, - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false, free_time=false) - ) - - before = prob.trajectory.global_data[phase_name] - solve!(prob, max_iter=50) - after = prob.trajectory.global_data[phase_name] - - @test before ≠ after - - @test unitary_rollout_fidelity( - prob.trajectory, - sys; - phases=prob.trajectory.global_data[phase_name], - phase_operators=phase_operators - ) > 0.9 - - @test unitary_rollout_fidelity(prob.trajectory, sys) < 0.9 -end +# TODO: Test changing names of control, state, and timestep diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl new file mode 100644 index 00000000..7e3a988b --- /dev/null +++ b/src/quantum_constraints.jl @@ -0,0 +1,59 @@ +module QuantumConstraints + +using ..QuantumObjectives + +using NamedTrajectories +using PiccoloQuantumObjects +using DirectTrajOpt + +export FinalKetFidelityConstraint +export FinalUnitaryFidelityConstraint + + +# --------------------------------------------------------- +# Kets +# --------------------------------------------------------- + +function FinalKetFidelityConstraint( + ψ_goal::AbstractVector{<:Complex{Float64}}, + ψ̃_name::Symbol, + final_fidelity::Float64, + traj::NamedTrajectory +) + terminal_constraint = ψ̃ -> [ + ket_infidelity_loss(ψ̃, ψ_goal) - abs(1 - final_fidelity) + ] + + return NonlinearKnotPointConstraint( + terminal_constraint, + ψ̃_name, + traj, + equality=false, + times=[traj.T] + ) +end + +# --------------------------------------------------------- +# Unitaries +# --------------------------------------------------------- + +function FinalUnitaryFidelityConstraint( + U_goal::AbstractPiccoloOperator, + Ũ⃗_name::Symbol, + final_fidelity::Float64, + traj::NamedTrajectory +) + terminal_constraint = Ũ⃗ -> [ + unitary_infidelity_loss(Ũ⃗, U_goal) - abs(1 - final_fidelity) + ] + + return NonlinearKnotPointConstraint( + terminal_constraint, + Ũ⃗_name, + traj, + equality=false, + times=[traj.T] + ) +end + +end \ No newline at end of file diff --git a/src/quantum_integrators.jl b/src/quantum_integrators.jl new file mode 100644 index 00000000..4c8e144a --- /dev/null +++ b/src/quantum_integrators.jl @@ -0,0 +1,42 @@ +module QuantumIntegrators + +export KetIntegrator +export UnitaryIntegrator +export DensityMatrixIntegrator + +using LinearAlgebra +using NamedTrajectories +using DirectTrajOpt +using PiccoloQuantumObjects + +const ⊗ = kron + +function KetIntegrator( + sys::QuantumSystem, + traj::NamedTrajectory, + ψ̃::Symbol, + a::Symbol +) + return BilinearIntegrator(sys.G, traj, ψ̃, a) +end + +function UnitaryIntegrator( + sys::QuantumSystem, + traj::NamedTrajectory, + Ũ⃗::Symbol, + a::Symbol +) + Ĝ = a_ -> I(sys.levels) ⊗ sys.G(a_) + return BilinearIntegrator(Ĝ, traj, Ũ⃗, a) +end + +function DensityMatrixIntegrator( + sys::OpenQuantumSystem, + traj::NamedTrajectory, + ρ̃::Symbol, + a::Symbol +) + return BilinearIntegrator(sys.𝒢, traj, ρ̃, a) +end + +end \ No newline at end of file diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl new file mode 100644 index 00000000..7c085bdd --- /dev/null +++ b/src/quantum_objectives.jl @@ -0,0 +1,98 @@ +module QuantumObjectives + +export ket_infidelity_loss +export unitary_infidelity_loss +export density_matrix_pure_state_infidelity_loss +export KetInfidelityLoss +export UnitaryInfidelityLoss +export DensityMatrixPureStateInfidelityLoss + +using LinearAlgebra +using NamedTrajectories +using PiccoloQuantumObjects +using DirectTrajOpt + +# --------------------------------------------------------- +# Kets +# --------------------------------------------------------- + +function ket_infidelity_loss( + ψ̃::AbstractVector, + ψ_goal::AbstractVector{<:Complex{Float64}} +) + ψ = iso_to_ket(ψ̃) + ℱ = abs2(ψ_goal' * ψ) + return abs(1 - ℱ) +end + +function KetInfidelityLoss( + ψ̃_name::Symbol, + traj::NamedTrajectory; + Q=100.0 +) + ψ_goal = iso_to_ket(traj.goal[ψ̃_name]) + ℓ = ψ̃ -> ket_infidelity_loss(ψ̃, ψ_goal) + return TerminalLoss(ℓ, ψ̃_name, traj; Q=Q) +end + + +# --------------------------------------------------------- +# Unitaries +# --------------------------------------------------------- + +function unitary_infidelity_loss( + Ũ⃗::AbstractVector, + U_goal::AbstractMatrix{<:Complex{Float64}} +) + U = iso_vec_to_operator(Ũ⃗) + n = size(U, 1) + ℱ = abs2(tr(U_goal' * U)) / n^2 + return abs(1 - ℱ) +end + +function unitary_infidelity_loss( + Ũ⃗::AbstractVector, + op::EmbeddedOperator +) + U_goal = unembed(op) + U = iso_vec_to_operator(Ũ⃗)[op.subspace, op.subspace] + n = length(op.subspace) + M = U_goal'U + ℱ = 1 / (n * (n + 1)) * (abs(tr(M'M)) + abs2(tr(M))) + return abs(1 - ℱ) +end + +function UnitaryInfidelityLoss( + U_goal::AbstractPiccoloOperator, + Ũ⃗_name::Symbol, + traj::NamedTrajectory; + Q=100.0 +) + ℓ = Ũ⃗ -> unitary_infidelity_loss(Ũ⃗, U_goal) + return TerminalLoss(ℓ, Ũ⃗_name, traj; Q=Q) +end + +# --------------------------------------------------------- +# Density Matrices +# --------------------------------------------------------- + +function density_matrix_pure_state_infidelity_loss( + ρ̃::AbstractVector, + ψ::AbstractVector{<:Complex{Float64}} +) + ρ = iso_vec_to_density(ρ̃) + ℱ = real(ψ' * ρ * ψ) + return abs(1 - ℱ) +end + +function DensityMatrixPureStateInfidelityLoss( + ρ̃_name::Symbol, + ψ_goal::AbstractVector{<:Complex{Float64}}, + traj::NamedTrajectory; + Q=100.0 +) + ℓ = ρ̃ -> density_matrix_pure_state_infidelity_loss(ρ̃, ψ_goal) + return TerminalLoss(ℓ, ρ̃_name, traj; Q=Q) +end + +end \ No newline at end of file diff --git a/src/quantum_system_templates/cats.jl b/src/quantum_system_templates/cats.jl index ca7c0f2e..8956d332 100644 --- a/src/quantum_system_templates/cats.jl +++ b/src/quantum_system_templates/cats.jl @@ -2,7 +2,7 @@ export CatSystem export coherent_ket export get_cat_controls -function coherent_ket(α::Union{Real, Complex}, levels::Int) +function coherent_ket(α::Union{Real, Complex}, levels::Int)::Vector{ComplexF64} return [exp(-0.5 * abs2(α)) * α^n / sqrt(factorial(n)) for n in 0:levels-1] end diff --git a/src/quantum_system_templates/rydberg.jl b/src/quantum_system_templates/rydberg.jl index 761423c7..3abd06bf 100644 --- a/src/quantum_system_templates/rydberg.jl +++ b/src/quantum_system_templates/rydberg.jl @@ -120,5 +120,7 @@ function RydbergChainSystem(; end @testitem "Rydberg system test" begin + using PiccoloQuantumObjects + @test RydbergChainSystem(N=3,cutoff_order=2,all2all=false) isa QuantumSystem end diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index cb38de92..3c329801 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -6,18 +6,13 @@ export unitary_linear_interpolation export initialize_trajectory using NamedTrajectories +using PiccoloQuantumObjects using Distributions using ExponentialAction using LinearAlgebra using TestItems -using ..Isomorphisms -using ..QuantumSystems -using ..Rollouts -using ..EmbeddedOperators -using ..DirectSums - # ----------------------------------------------------------------------------- # # Initial states # @@ -267,7 +262,7 @@ function initialize_trajectory( n_drives::Int, control_bounds::Tuple{Vararg{VectorBound}}; bound_state=false, - free_time=false, + free_time=true, control_name=:a, n_control_derivatives::Int=length(control_bounds) - 1, zero_initial_and_final_derivative=false, @@ -290,7 +285,7 @@ function initialize_trajectory( Symbol("d"^i * string(control_name)) for i = 1:n_control_derivatives ] if verbose - println("control_derivative_names: $control_derivative_names") + println("\tcontrol derivative names: $control_derivative_names") end control_names = (control_name, control_derivative_names...) @@ -587,6 +582,7 @@ end @testitem "Geodesic" begin using LinearAlgebra + using PiccoloQuantumObjects ## Group 1: identity to X (π rotation) @@ -643,6 +639,8 @@ end @testitem "unitary trajectory initialization" begin using NamedTrajectories + using PiccoloQuantumObjects + U_goal = GATES[:X] T = 10 Δt = 0.1