diff --git a/Project.toml b/Project.toml index 939f1d8f..8aff5bdd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumCollocation" uuid = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf" authors = ["Aaron Trowbridge and contributors"] -version = "0.7.0" +version = "0.7.2" [deps] DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" @@ -21,14 +21,14 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8" [compat] -DirectTrajOpt = "0.1.0" +DirectTrajOpt = "0.2" Distributions = "0.25" ExponentialAction = "0.2" Interpolations = "0.15" JLD2 = "0.5" LinearAlgebra = "1.10, 1.11" NamedTrajectories = "0.3" -PiccoloQuantumObjects = "0.3" +PiccoloQuantumObjects = "0.4" Random = "1.10, 1.11" Reexport = "1.2" SparseArrays = "1.10, 1.11" diff --git a/docs/Project.toml b/docs/Project.toml deleted file mode 100644 index 8df91abc..00000000 --- a/docs/Project.toml +++ /dev/null @@ -1,8 +0,0 @@ -[deps] -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" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/docs/dev/ipopt_callbacks.jl b/docs/dev/ipopt_callbacks.jl deleted file mode 100644 index fa072deb..00000000 --- a/docs/dev/ipopt_callbacks.jl +++ /dev/null @@ -1,88 +0,0 @@ -# ```@meta -# CollapsedDocStrings = true -# ``` -# # 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). - -# ## Callbacks - -using QuantumCollocation -using NamedTrajectories - -import ..QuantumStateSmoothPulseProblem -import ..Callbacks - -# 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, - obj_value::Float64, - inf_pr::Float64, - inf_du::Float64, - mu::Float64, - d_norm::Float64, - regularization_size::Float64, - alpha_du::Float64, - alpha_pr::Float64, - ls_trials::Cint, -) - return true -end - -# This gives the user access to some of the optimization state internals at each iteration. -# A callback function with any subset of these arguments can be passed into the `solve!` function via the `callback` keyword argument see below. - -# The callback function can be used to stop the optimization early by returning `false`. The following callback when passed to `solve!` will stop the optimization after the first iteration: -my_callback = (kwargs...) -> false - -# Single initial and target states -# -------------------------------- -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) -) - - -# The callback function can be used to monitor the optimization progress, save intermediate results, or modify the optimization process. -# For example, the following callback function saves the optimization trajectory at each iteration - this can be useful for debugging or plotting the optimization progress. -# `trajectory_history_callback` from the `Callbacks` module -callback, trajectory_history = QuantumCollocation.Callbacks.trajectory_history_callback(prob) -solve!(prob, max_iter=20, callback=callback) - -# Save trajectory images into files which can be used to create a gif like the following: -for (iter, traj) in enumerate(trajectory_history) - str_index = lpad(iter, length(string(length(trajectory_history))), "0") - plot("./iteration-$str_index-trajectory.png", traj, [:ψ̃, :a], xlims=(-Δt, (T+5)*Δt), ylims=(ψ̃1 = (-2, 2), a = (-1.1, 1.1))) -end - -# ![pulse optimization animation](../../assets/animation.gif) - -# Using a callback to get the best trajectory from all the optimization iterations -sys2 = QuantumSystem(0.15 * GATES[:Z], [GATES[:X], GATES[:Y]]) -ψ_init2 = Vector{ComplexF64}([0.0, 1.0]) -ψ_target2 = Vector{ComplexF64}([1.0, 0.0]) - -# Using other callbacks from the callback library -# -------------------------------- -# Callback used here is `best_rollout_fidelity_callback` which appends the best trajectories based on monotonically increasing fidelity of the rollout -prob2 = QuantumStateSmoothPulseProblem( - sys2, ψ_init2, ψ_target2, T, Δt; - ipopt_options=IpoptOptions(print_level=1), - piccolo_options=PiccoloOptions(verbose=false) -) - -best_trajectory_callback, best_trajectory_list = best_rollout_fidelity_callback(prob2) -solve!(prob2, max_iter=20, callback=best_trajectory_callback) -# fidelity of the last iterate -@show Losses.fidelity(prob2) - -# fidelity of the best iterate -@show QuantumCollocation.fidelity(best_trajectory_list[end], prob2.system) diff --git a/docs/literate/man/problem_templates.jl b/docs/literate/man/problem_templates.jl deleted file mode 100644 index 75107fd0..00000000 --- a/docs/literate/man/problem_templates.jl +++ /dev/null @@ -1,197 +0,0 @@ -# ```@meta -# CollapsedDocStrings = true -# ``` -using NamedTrajectories -using PiccoloQuantumObjects -using QuantumCollocation - -#= -# 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 -`DirectTrajOptProblem` object, which stores all the parts of the optimal control problem. - -This page provides a brief overview of each problem template, broken down by the state of -the problem being solved. - -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)]_ - - diff --git a/docs/make.jl b/docs/make.jl deleted file mode 100644 index 7323c6e8..00000000 --- a/docs/make.jl +++ /dev/null @@ -1,72 +0,0 @@ -using QuantumCollocation -using Documenter -using Literate - -push!(LOAD_PATH, joinpath(@__DIR__, "..", "src")) - -@info "Building Documenter site for QuantumCollocation.jl" -open(joinpath(@__DIR__, "src", "index.md"), write = true) do io - for line in eachline(joinpath(@__DIR__, "..", "README.md")) - if occursin("", line) - comment_content = match(r"", line).captures[1] - write(io, comment_content * "\n") - else - write(io, line * "\n") - end - end -end - -# TODO: Callbacks are currently broken -pages = [ - "Home" => "index.md", - "Manual" => [ - "Problem Templates" => "generated/man/problem_templates.md", - ], - "Library" => "lib.md", -] - -format = Documenter.HTML(; - prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://docs.harmoniqs.co/QuantumCollocation.jl", - edit_link="main", - assets=String[], - mathengine = MathJax3(Dict( - :loader => Dict("load" => ["[tex]/physics"]), - :tex => Dict( - "inlineMath" => [["\$","\$"], ["\\(","\\)"]], - "tags" => "ams", - "packages" => [ - "base", - "ams", - "autoload", - "physics" - ], - ), - )), -) - -src = joinpath(@__DIR__, "src") -lit = joinpath(@__DIR__, "literate") - -lit_output = joinpath(src, "generated") - -for (root, _, files) ∈ walkdir(lit), file ∈ files - splitext(file)[2] == ".jl" || continue - ipath = joinpath(root, file) - opath = splitdir(replace(ipath, lit=>lit_output))[1] - Literate.markdown(ipath, opath) -end - -makedocs(; - modules=[QuantumCollocation], - authors="Aaron Trowbridge and contributors", - sitename="QuantumCollocation.jl", - format=format, - pages=pages, - warnonly=true, -) - -deploydocs(; - repo="github.com/harmoniqs/QuantumCollocation.jl.git", - devbranch="main", -) diff --git a/docs/src/assets/animation.gif b/docs/src/assets/animation.gif deleted file mode 100644 index 69e1adf7..00000000 Binary files a/docs/src/assets/animation.gif and /dev/null differ diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg deleted file mode 100644 index 9ee565d5..00000000 --- a/docs/src/assets/logo.svg +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/docs/src/assets/piccolo_logo.svg b/docs/src/assets/piccolo_logo.svg deleted file mode 100644 index d91bee46..00000000 --- a/docs/src/assets/piccolo_logo.svg +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md deleted file mode 100644 index 83d0ef82..00000000 --- a/docs/src/index.md +++ /dev/null @@ -1,116 +0,0 @@ -```@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 deleted file mode 100644 index dc0310d9..00000000 --- a/docs/src/lib.md +++ /dev/null @@ -1,26 +0,0 @@ -# Library - -## Problem Templates -```@autodocs -Modules = [QuantumCollocation.ProblemTemplates] -``` - -## Options -```@autodocs -Modules = [QuantumCollocation.Options] -``` - -## Trajectory Initialization -```@autodocs -Modules = [QuantumCollocation.TrajectoryInitialization] -``` - -## Trajectory Interpolations -```@autodocs -Modules = [QuantumCollocation.TrajectoryInterpolations] -``` - -## Quantum System Templates -```@autodocs -Modules = [QuantumCollocation.QuantumSystemTemplates] -``` \ No newline at end of file diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index fd6e512d..ce6e09ab 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -1,5 +1,7 @@ module ProblemTemplates +export apply_piccolo_options! + using ..TrajectoryInitialization using ..QuantumObjectives using ..QuantumConstraints @@ -17,6 +19,7 @@ using JLD2 using TestItems include("unitary_smooth_pulse_problem.jl") +include("unitary_variational_problem.jl") include("unitary_minimum_time_problem.jl") include("unitary_sampling_problem.jl") diff --git a/src/problem_templates/quantum_state_sampling_problem.jl b/src/problem_templates/quantum_state_sampling_problem.jl index 0859e3d1..710cded6 100644 --- a/src/problem_templates/quantum_state_sampling_problem.jl +++ b/src/problem_templates/quantum_state_sampling_problem.jl @@ -88,7 +88,7 @@ function QuantumStateSamplingProblem( for (weight, names) in zip(system_weights, state_names) for name in names - J += KetInfidelityLoss(name, traj, Q=weight * Q) + J += KetInfidelityObjective(name, traj, Q=weight * Q) end end diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index f5d9798a..6cd18912 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -128,7 +128,7 @@ function QuantumStateSmoothPulseProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) for name ∈ state_names - J += KetInfidelityLoss(name, traj; Q=Q) + J += KetInfidelityObjective(name, traj; Q=Q) end # Optional Piccolo constraints and objectives diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index 7e334079..f89451da 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -120,7 +120,7 @@ function UnitarySamplingProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) for (weight, op, name) in zip(system_weights, operators, state_names) - J += UnitaryInfidelityLoss(op, name, traj; Q=weight * Q) + J += UnitaryInfidelityObjective(op, name, traj; Q=weight * Q) end # Optional Piccolo constraints and objectives diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 107d0de4..23bd639c 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -122,7 +122,7 @@ function UnitarySmoothPulseProblem( end # Objective - J = UnitaryInfidelityLoss(goal, state_name, traj; Q=Q) + J = UnitaryInfidelityObjective(goal, state_name, traj; Q=Q) control_names = [ name for name ∈ traj.names @@ -229,4 +229,4 @@ end @test final > initial end -# TODO: Test changing names of control, state, and timestep +# TODO: Test changing names of control, state, and timestep \ No newline at end of file diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl new file mode 100644 index 00000000..44d3d3ec --- /dev/null +++ b/src/problem_templates/unitary_variational_problem.jl @@ -0,0 +1,207 @@ +export UnitaryVariationalProblem + + +""" + UnitaryVariationalProblem( + system::VariationalQuantumSystem, + goal::AbstractPiccoloOperator, + T::Int, + Δt::Union{Float64, <:AbstractVector{Float64}}; + robust_times::AbstractVector{<:Union{AbstractVector{Int}, Nothing}}=[nothing for s ∈ system.G_vars], + sensitive_times::AbstractVector{<:Union{AbstractVector{Int}, Nothing}}=[nothing for s ∈ system.G_vars], + kwargs... + ) + +Constructs a unitary variational problem for optimizing quantum control trajectories. + +# Arguments +- `system::VariationalQuantumSystem`: The quantum system to be controlled, containing variational parameters. +- `goal::AbstractPiccoloOperator`: The target operator or state to achieve at the end of the trajectory. +- `T::Int`: The total number of timesteps in the trajectory. +- `Δt::Union{Float64, <:AbstractVector{Float64}}`: The timestep duration or a vector of timestep durations. +- `robust_times::AbstractVector`: Times at which robustness to variations in the trajectory is enforced. +- `sensitive_times::AbstractVector`: Times at which sensitivity to variations in the trajectory is enhanced. +- `unitary_integrator`: The integrator used for unitary evolution (default: `VariationalUnitaryIntegrator`). +- `state_name::Symbol`: The name of the state variable in the trajectory (default: `:Ũ⃗`). +- `variational_state_name::Symbol`: The name of the variational state variable (default: `:Ũ⃗ₐ`). +- `control_name::Symbol`: The name of the control variable (default: `:a`). +- `timestep_name::Symbol`: The name of the timestep variable (default: `:Δt`). +- `init_trajectory::Union{NamedTrajectory, Nothing}`: An optional initial trajectory to start optimization. +- `a_bound::Float64`: The bound for the control variable `a` (default: `1.0`). +- `a_bounds::Vector`: Bounds for each control variable (default: filled with `a_bound`). +- `da_bound::Float64`: The bound for the derivative of the control variable (default: `Inf`). +- `da_bounds::Vector`: Bounds for each derivative of the control variable. +- `dda_bound::Float64`: The bound for the second derivative of the control variable (default: `1.0`). +- `dda_bounds::Vector`: Bounds for each second derivative of the control variable. +- `Δt_min::Float64`: Minimum allowed timestep duration. +- `Δt_max::Float64`: Maximum allowed timestep duration. +- `Q::Float64`: Weight for the unitary infidelity objective (default: `100.0`). +- `Q_v::Float64`: Weight for sensitivity objectives (default: `1.0`). +- `R`: Regularization weight for control variables (default: `1e-2`). +- `R_a`, `R_da`, `R_dda`: Regularization weights for control, its derivative, and second derivative. +- `constraints::Vector`: Additional constraints for the optimization problem. +- `piccolo_options::PiccoloOptions`: Options for configuring the Piccolo optimization framework. + +# Returns +A `DirectTrajOptProblem` object representing the optimization problem, including the +trajectory, objective, integrators, and constraints. + +# Notes +This function constructs a trajectory optimization problem for quantum control using +variational principles. It supports robust and sensitive trajectory design, regularization, +and optional constraints. The problem is solved using the Piccolo optimization framework. + +""" +function UnitaryVariationalProblem end + +function UnitaryVariationalProblem( + system::VariationalQuantumSystem, + goal::AbstractPiccoloOperator, + T::Int, + Δt::Union{Float64, <:AbstractVector{Float64}}; + robust_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars], + sensitive_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars], + variational_integrator=VariationalUnitaryIntegrator, + variational_scale::Float64=1.0, + state_name::Symbol = :Ũ⃗, + variational_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), + da_bound::Float64=Inf, + da_bounds::Vector{Float64}=fill(da_bound, system.n_drives), + 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), + Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), + Q::Float64=100.0, + Q_s::Float64=1e-2, + Q_r::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, + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], + piccolo_options::PiccoloOptions=PiccoloOptions(), +) + if piccolo_options.verbose + println(" constructing UnitaryVariationalProblem...") + println("\tusing integrator: $(typeof(variational_integrator))") + println("\ttotal variational parameters: $(length(system.G_vars))") + end + + # Trajectory + if !isnothing(init_trajectory) + traj = deepcopy(init_trajectory) + else + traj = initialize_trajectory( + goal, + T, + Δt, + system.n_drives, + (a_bounds, da_bounds, dda_bounds); + state_name=state_name, + control_name=control_name, + timestep_name=timestep_name, + Δt_bounds=(Δt_min, Δt_max), + zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, + geodesic=piccolo_options.geodesic, + bound_state=piccolo_options.bound_state, + rollout_integrator=piccolo_options.rollout_integrator, + verbose=piccolo_options.verbose + ) + end + + _, Ũ⃗_vars = variational_unitary_rollout( + traj, + system; + unitary_name=state_name, + drive_name=control_name + ) + + variational_state_names = [ + Symbol(string(variational_state_name) * "$(i)") for i in eachindex(system.G_vars) + ] + + for (name, Ũ⃗_v) in zip(variational_state_names, Ũ⃗_vars) + add_component!(traj, name, Ũ⃗_v / variational_scale; type=:state) + traj.initial = merge(traj.initial, (name => Ũ⃗_v[:, 1] / variational_scale, )) + end + + control_names = [ + name for name ∈ traj.names + if endswith(string(name), string(control_name)) + ] + + # objective + J = UnitaryInfidelityObjective(goal, state_name, traj; Q=Q) + J += QuadraticRegularizer(control_names[1], traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) + + # sensitivity + for (name, s, r) ∈ zip(variational_state_names, sensitive_times, robust_times) + @assert isdisjoint(s, r) + J += UnitarySensitivityObjective( + name, + traj, + [s; r]; + Qs=[fill(-Q_s, length(s)); fill(Q_r, length(r))], + scale=variational_scale) + end + + # Optional Piccolo constraints and objectives + apply_piccolo_options!( + J, constraints, piccolo_options, traj, state_name, timestep_name; + state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing + ) + + integrators = [ + variational_integrator(system, traj, state_name, variational_state_names, control_name, scale=variational_scale), + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), + ] + + return DirectTrajOptProblem( + traj, + J, + integrators; + constraints=constraints + ) +end + +# *************************************************************************** # + +@testitem "Sensitive and robust" begin + using LinearAlgebra + using PiccoloQuantumObjects + + system = QuantumSystem([PAULIS.X, PAULIS.Y]) + varsys = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X] ) + T = 50 + Δt = 0.2 + + sense_scale = 8.0 + sense_prob = UnitaryVariationalProblem( + varsys, GATES.X, T, Δt, + variational_scale=sense_scale, + sensitive_times=[[T]], + piccolo_options=PiccoloOptions(verbose=false) + ) + solve!(sense_prob, max_iter=20, print_level=1, verbose=false) + + rob_scale = 1 / 8.0 + rob_prob = UnitaryVariationalProblem( + varsys, GATES.X, T, Δt, + variational_scale=rob_scale, + robust_times=[[T]], + piccolo_options=PiccoloOptions(verbose=false) + ) + solve!(rob_prob, max_iter=20, print_level=1, verbose=false) + + sense_n = norm(sense_scale * sense_prob.trajectory.Ũ⃗ᵥ1) + rob_n = norm(rob_scale * rob_prob.trajectory.Ũ⃗ᵥ1) + @assert sense_n > rob_n +end \ No newline at end of file diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index 7e3a988b..2e77dafd 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -21,7 +21,7 @@ function FinalKetFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = ψ̃ -> [ - ket_infidelity_loss(ψ̃, ψ_goal) - abs(1 - final_fidelity) + abs(QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal) - final_fidelity) ] return NonlinearKnotPointConstraint( @@ -44,7 +44,7 @@ function FinalUnitaryFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = Ũ⃗ -> [ - unitary_infidelity_loss(Ũ⃗, U_goal) - abs(1 - final_fidelity) + abs(QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) - final_fidelity) ] return NonlinearKnotPointConstraint( diff --git a/src/quantum_integrators.jl b/src/quantum_integrators.jl index 4c8e144a..c0172973 100644 --- a/src/quantum_integrators.jl +++ b/src/quantum_integrators.jl @@ -3,14 +3,20 @@ module QuantumIntegrators export KetIntegrator export UnitaryIntegrator export DensityMatrixIntegrator +export VariationalUnitaryIntegrator using LinearAlgebra using NamedTrajectories using DirectTrajOpt using PiccoloQuantumObjects +using SparseArrays const ⊗ = kron +# ----------------------------------------------------------------------------- # +# Default Integrators +# ----------------------------------------------------------------------------- # + function KetIntegrator( sys::QuantumSystem, traj::NamedTrajectory, @@ -39,4 +45,40 @@ function DensityMatrixIntegrator( return BilinearIntegrator(sys.𝒢, traj, ρ̃, a) end +# ----------------------------------------------------------------------------- # +# Variational Integrators +# ----------------------------------------------------------------------------- # + +function VariationalKetIntegrator( + sys::VariationalQuantumSystem, + traj::NamedTrajectory, + ψ̃::Symbol, + ψ̃_variations::AbstractVector{Symbol}, + a::Symbol; + scale::Float64=1.0, +) + var_ψ̃ = vcat(ψ̃, ψ̃_variations...) + G = a -> Isomorphisms.var_G(sys.G(a), [G(a) / scale for G in sys.G_vars]) + return BilinearIntegrator(G, traj, var_ψ̃, a) +end + +function VariationalUnitaryIntegrator( + sys::VariationalQuantumSystem, + traj::NamedTrajectory, + Ũ⃗::Symbol, + Ũ⃗_variations::AbstractVector{Symbol}, + a::Symbol; + scale::Float64=1.0, +) + var_Ũ⃗ = vcat(Ũ⃗, Ũ⃗_variations...) + + function Ĝ(a) + G0 = sys.G(a) + Gs = typeof(G0)[I(sys.levels) ⊗ G(a) / scale for G in sys.G_vars] + return Isomorphisms.var_G(I(sys.levels) ⊗ G0, Gs) + end + return BilinearIntegrator(Ĝ, traj, var_Ũ⃗, a) +end + + end \ No newline at end of file diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index 7c085bdd..21ed1cc8 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -1,11 +1,9 @@ module QuantumObjectives -export ket_infidelity_loss -export unitary_infidelity_loss -export density_matrix_pure_state_infidelity_loss -export KetInfidelityLoss -export UnitaryInfidelityLoss -export DensityMatrixPureStateInfidelityLoss +export KetInfidelityObjective +export UnitaryInfidelityObjective +export DensityMatrixPureStateInfidelityObjective +export UnitarySensitivityObjective using LinearAlgebra using NamedTrajectories @@ -16,23 +14,22 @@ using DirectTrajOpt # Kets # --------------------------------------------------------- -function ket_infidelity_loss( +function ket_fidelity_loss( ψ̃::AbstractVector, ψ_goal::AbstractVector{<:Complex{Float64}} ) ψ = iso_to_ket(ψ̃) - ℱ = abs2(ψ_goal' * ψ) - return abs(1 - ℱ) + return abs2(ψ_goal' * ψ) end -function KetInfidelityLoss( +function KetInfidelityObjective( ψ̃_name::Symbol, traj::NamedTrajectory; Q=100.0 ) ψ_goal = iso_to_ket(traj.goal[ψ̃_name]) - ℓ = ψ̃ -> ket_infidelity_loss(ψ̃, ψ_goal) - return TerminalLoss(ℓ, ψ̃_name, traj; Q=Q) + ℓ = ψ̃ -> abs(1 - ket_fidelity_loss(ψ̃, ψ_goal)) + return TerminalObjective(ℓ, ψ̃_name, traj; Q=Q) end @@ -40,36 +37,34 @@ end # Unitaries # --------------------------------------------------------- -function unitary_infidelity_loss( - Ũ⃗::AbstractVector, +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real}, U_goal::AbstractMatrix{<:Complex{Float64}} ) U = iso_vec_to_operator(Ũ⃗) n = size(U, 1) - ℱ = abs2(tr(U_goal' * U)) / n^2 - return abs(1 - ℱ) + return abs2(tr(U_goal' * U)) / n^2 end -function unitary_infidelity_loss( - Ũ⃗::AbstractVector, +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real}, 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 - ℱ) + return 1 / (n * (n + 1)) * (abs(tr(M'M)) + abs2(tr(M))) end -function UnitaryInfidelityLoss( +function UnitaryInfidelityObjective( U_goal::AbstractPiccoloOperator, Ũ⃗_name::Symbol, traj::NamedTrajectory; Q=100.0 ) - ℓ = Ũ⃗ -> unitary_infidelity_loss(Ũ⃗, U_goal) - return TerminalLoss(ℓ, Ũ⃗_name, traj; Q=Q) + ℓ = Ũ⃗ -> abs(1 - unitary_fidelity_loss(Ũ⃗, U_goal)) + return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q) end # --------------------------------------------------------- @@ -85,14 +80,45 @@ function density_matrix_pure_state_infidelity_loss( return abs(1 - ℱ) end -function DensityMatrixPureStateInfidelityLoss( +function DensityMatrixPureStateInfidelityObjective( ρ̃_name::Symbol, ψ_goal::AbstractVector{<:Complex{Float64}}, traj::NamedTrajectory; Q=100.0 ) ℓ = ρ̃ -> density_matrix_pure_state_infidelity_loss(ρ̃, ψ_goal) - return TerminalLoss(ℓ, ρ̃_name, traj; Q=Q) + return TerminalObjective(ℓ, ρ̃_name, traj; Q=Q) +end + +# --------------------------------------------------------- +# Sensitivity +# --------------------------------------------------------- + +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real} +) + U = iso_vec_to_operator(Ũ⃗) + n = size(U, 1) + return abs2(tr(U' * U)) / n^2 end +function UnitarySensitivityObjective( + name::Symbol, + traj::NamedTrajectory, + times::AbstractVector{Int}; + Qs::AbstractVector{<:Float64}=fill(1.0, length(times)), + scale::Float64=1.0, +) + ℓ = Ũ⃗ -> scale^4 * unitary_fidelity_loss(Ũ⃗) + + return KnotPointObjective( + ℓ, + name, + traj; + Qs=Qs, + times=times + ) +end + + end \ No newline at end of file diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 3c329801..7ab0cdd3 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -424,6 +424,7 @@ function initialize_trajectory( # Construct phase data phase_data = isnothing(phase_operators) ? nothing : π * randn(length(phase_operators)) + return initialize_trajectory( [Ũ⃗_traj], [Ũ⃗_init], @@ -438,6 +439,8 @@ function initialize_trajectory( ) end + + """ initialize_trajectory