diff --git a/Project.toml b/Project.toml index 2a87a8b9..33ee19e2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,12 @@ name = "QuantumCollocation" uuid = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf" +version = "0.9.0" authors = ["Aaron Trowbridge and contributors"] -version = "0.8.1" [deps] DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08" @@ -20,13 +19,12 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8" [compat] -DirectTrajOpt = "0.4" +DirectTrajOpt = "0.5" Distributions = "0.25" ExponentialAction = "0.2" -JLD2 = "0.5" LinearAlgebra = "1.10, 1.11, 1.12" -NamedTrajectories = "0.5" -PiccoloQuantumObjects = "0.6" +NamedTrajectories = "0.6" +PiccoloQuantumObjects = "0.7" Random = "1.10, 1.11, 1.12" Reexport = "1.2" SparseArrays = "1.10, 1.11, 1.12" diff --git a/docs/literate/examples/multilevel_transmon.jl b/docs/literate/examples/multilevel_transmon.jl index ca837510..f59d074c 100644 --- a/docs/literate/examples/multilevel_transmon.jl +++ b/docs/literate/examples/multilevel_transmon.jl @@ -44,7 +44,7 @@ levels = 5 δ = 0.2 ## add a bound to the controls -a_bound = 0.2 +u_bound = 0.2 ## create the system sys = TransmonSystem(levels=levels, δ=δ) @@ -75,7 +75,7 @@ get_subspace_identity(op) |> sparse # We can then pass this embedded operator to the `UnitarySmoothPulseProblem` template to create the problem ## create the problem -prob = UnitarySmoothPulseProblem(sys, op, T, Δt; a_bound=a_bound) +prob = UnitarySmoothPulseProblem(sys, op, T, Δt; u_bound=u_bound) ## solve the problem solve!(prob; max_iter=50) @@ -98,8 +98,8 @@ plot_unitary_populations(prob.trajectory; fig_size=(900, 700)) ## create the a leakage suppression problem, initializing with the previous solution prob_leakage = UnitarySmoothPulseProblem(sys, op, T, Δt; - a_bound=a_bound, - a_guess=prob.trajectory.a[:, :], + u_bound=u_bound, + u_guess=prob.trajectory.u[:, :], piccolo_options=PiccoloOptions( leakage_constraint=true, leakage_constraint_value=1e-2, diff --git a/docs/literate/examples/two_qubit_gates.jl b/docs/literate/examples/two_qubit_gates.jl index 47d9b549..40aeaa1e 100644 --- a/docs/literate/examples/two_qubit_gates.jl +++ b/docs/literate/examples/two_qubit_gates.jl @@ -13,7 +13,7 @@ # Specifically, for a simple two-qubit system in a rotating frame, we have # ```math -# H = J_{12} \sigma_1^x \sigma_2^x + \sum_{i \in {1,2}} a_i^R(t) {\sigma^x_i \over 2} + a_i^I(t) {\sigma^y_i \over 2}. +# H = J_{12} \sigma_1^x \sigma_2^x + \sum_{i \in {1,2}} u_i^R(t) {\sigma^x_i \over 2} + u_i^I(t) {\sigma^y_i \over 2}. # ``` # where @@ -21,7 +21,7 @@ # ```math # \begin{align*} # J_{12} &= 0.001 \text{ GHz}, \\ -# |a_i^R(t)| &\leq 0.1 \text{ GHz} \\ +# |u_i^R(t)| &\leq 0.1 \text{ GHz} \\ # \end{align*} # ``` @@ -53,7 +53,7 @@ Id = GATES[:I] ## Define the parameters of the Hamiltonian J_12 = 0.001 # GHz -a_bound = 0.100 # GHz +u_bound = 0.100 # GHz ## Define the drift (coupling) Hamiltonian H_drift = J_12 * (σx ⊗ σx) @@ -62,9 +62,9 @@ H_drift = J_12 * (σx ⊗ σx) H_drives = [σx_1 / 2, σy_1 / 2, σx_2 / 2, σy_2 / 2] ## Define control (and higher derivative) bounds -a_bound = 0.1 -da_bound = 0.0005 -dda_bound = 0.0025 +u_bound = 0.1 +du_bound = 0.0005 +ddu_bound = 0.0025 ## Scale the Hamiltonians by 2π H_drift *= 2π @@ -95,11 +95,11 @@ prob = UnitarySmoothPulseProblem( U_goal, T, Δt; - a_bound=a_bound, - da_bound=da_bound, - dda_bound=dda_bound, - R_da=0.01, - R_dda=0.01, + u_bound=u_bound, + du_bound=du_bound, + ddu_bound=ddu_bound, + R_du=0.01, + R_ddu=0.01, Δt_max=Δt_max, piccolo_options=PiccoloOptions(bound_state=true), ) @@ -166,11 +166,11 @@ prob = UnitarySmoothPulseProblem( U_goal, T, Δt; - a_bound=a_bound, - da_bound=da_bound, - dda_bound=dda_bound, - R_da=0.01, - R_dda=0.01, + u_bound=u_bound, + du_bound=du_bound, + ddu_bound=ddu_bound, + R_du=0.01, + R_ddu=0.01, Δt_max=Δt_max, piccolo_options=PiccoloOptions(bound_state=true), ) diff --git a/docs/literate/man/ket_problem_templates.jl b/docs/literate/man/ket_problem_templates.jl index 9d2283ba..c0245988 100644 --- a/docs/literate/man/ket_problem_templates.jl +++ b/docs/literate/man/ket_problem_templates.jl @@ -39,7 +39,7 @@ solve!(state_prob, max_iter=100, verbose=true, print_level=1); println("After: ", rollout_fidelity(state_prob.trajectory, system)) # _extract the control pulses_ -state_prob.trajectory.a |> size +state_prob.trajectory.u |> size # ----- diff --git a/docs/literate/man/unitary_problem_templates.jl b/docs/literate/man/unitary_problem_templates.jl index 2121f883..6c2e140c 100644 --- a/docs/literate/man/unitary_problem_templates.jl +++ b/docs/literate/man/unitary_problem_templates.jl @@ -41,7 +41,7 @@ The `NamedTrajectory` object stores the control pulse, state variables, and the =# # _extract the control pulses_ -prob.trajectory.a |> size +prob.trajectory.u |> size # ----- diff --git a/docs/src/index.md b/docs/src/index.md index aaa32321..3aa93fcc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -36,26 +36,26 @@ Unitary 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, +For example, a *UnitarySmoothPulseProblem* is tasked with generating a *pulse* sequence $u_{1:T-1}$ in order 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 \\ + \qquad & U_{t+1} = \exp\{- i H(u_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 \\ + \qquad & U_{t+1} = \exp\{- i H(u_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 trajectory optimization*. +In each case, the dynamics between *knot points* $(U_t, u_t)$ and $(U_{t+1}, u_{t+1})$ are enforced as constraints on the states, which are free variables in the solver; this optimization framework is called *direct trajectory optimization*. ----- diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 2a263ab7..6656b57f 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -12,7 +12,6 @@ using DirectTrajOpt using PiccoloQuantumObjects using ExponentialAction -using JLD2 using LinearAlgebra using SparseArrays using TestItems @@ -53,7 +52,7 @@ function apply_piccolo_options!( end for (name, indices) ∈ zip(state_names, state_leakage_indices) - J += LeakageObjective(indices, name, traj, Qs=fill(piccolo_options.leakage_cost, traj.T)) + J += LeakageObjective(indices, name, traj, Qs=fill(piccolo_options.leakage_cost, traj.N)) push!(constraints, LeakageConstraint(val, indices, name, traj)) end end @@ -73,7 +72,7 @@ function apply_piccolo_options!( 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], + u -> [norm(u)^2 - piccolo_options.complex_control_norm_constraint_radius^2], piccolo_options.complex_control_norm_constraint_name, traj; equality=false, diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index b73e9eca..429f8686 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -89,14 +89,14 @@ end using NamedTrajectories using PiccoloQuantumObjects - T = 51 + N = 51 Δt = 0.2 - sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) ψ_init = Vector{ComplexF64}[[1.0, 0.0]] ψ_target = Vector{ComplexF64}[[0.0, 1.0]] prob = QuantumStateSmoothPulseProblem( - sys, ψ_init, ψ_target, T, Δt; + sys, ψ_init, ψ_target, N, Δt; piccolo_options=PiccoloOptions(verbose=false) ) initial = sum(get_timesteps(prob.trajectory)) diff --git a/src/problem_templates/quantum_state_sampling_problem.jl b/src/problem_templates/quantum_state_sampling_problem.jl index 49c4db81..4547b807 100644 --- a/src/problem_templates/quantum_state_sampling_problem.jl +++ b/src/problem_templates/quantum_state_sampling_problem.jl @@ -1,7 +1,41 @@ export QuantumStateSamplingProblem """ + QuantumStateSamplingProblem(systems, ψ_inits, ψ_goals, T, Δt; kwargs...) +Construct a quantum state sampling problem for multiple systems with shared controls. + +# Arguments +- `systems::AbstractVector{<:AbstractQuantumSystem}`: A vector of quantum systems. +- `ψ_inits::AbstractVector{<:AbstractVector{<:AbstractVector{<:ComplexF64}}}`: Initial states for each system. +- `ψ_goals::AbstractVector{<:AbstractVector{<:AbstractVector{<:ComplexF64}}}`: Target states for each system. +- `T::Int`: The number of time steps. +- `Δt::Union{Float64, AbstractVector{Float64}}`: The time step value or vector of time steps. + +# Keyword Arguments +- `ket_integrator=KetIntegrator`: The integrator to use for state dynamics. +- `system_weights=fill(1.0, length(systems))`: The weights for each system. +- `init_trajectory::Union{NamedTrajectory,Nothing}=nothing`: The initial trajectory. +- `state_name::Symbol=:ψ̃`: The name of the state variable. +- `control_name::Symbol=:u`: The name of the control variable. +- `timestep_name::Symbol=:Δt`: The name of the timestep variable. +- `u_bound::Float64=1.0`: The bound for the control amplitudes. +- `u_bounds=fill(u_bound, systems[1].n_drives)`: The bounds for the control amplitudes. +- `u_guess::Union{Matrix{Float64},Nothing}=nothing`: The initial guess for the control amplitudes. +- `du_bound::Float64=Inf`: The bound for the control first derivatives. +- `du_bounds=fill(du_bound, systems[1].n_drives)`: The bounds for the control first derivatives. +- `ddu_bound::Float64=1.0`: The bound for the control second derivatives. +- `ddu_bounds=fill(ddu_bound, systems[1].n_drives)`: The bounds for the control second derivatives. +- `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt)`: The minimum time step size. +- `Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt)`: The maximum time step size. +- `Q::Float64=100.0`: The fidelity weight. +- `R::Float64=1e-2`: The regularization weight. +- `R_u::Union{Float64,Vector{Float64}}=R`: The regularization weight for the control amplitudes. +- `R_du::Union{Float64,Vector{Float64}}=R`: The regularization weight for the control first derivatives. +- `R_ddu::Union{Float64,Vector{Float64}}=R`: The regularization weight for the control second derivatives. +- `state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing`: Indices of leakage states. +- `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: The constraints. +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. """ function QuantumStateSamplingProblem end @@ -15,22 +49,22 @@ function QuantumStateSamplingProblem( system_weights=fill(1.0, length(systems)), init_trajectory::Union{NamedTrajectory,Nothing}=nothing, state_name::Symbol=:ψ̃, - control_name::Symbol=:a, + control_name::Symbol=:u, timestep_name::Symbol=:Δt, - a_bound::Float64=1.0, - a_bounds=fill(a_bound, systems[1].n_drives), - a_guess::Union{Matrix{Float64},Nothing}=nothing, - da_bound::Float64=Inf, - da_bounds=fill(da_bound, systems[1].n_drives), - dda_bound::Float64=1.0, - dda_bounds=fill(dda_bound, systems[1].n_drives), - Δt_min::Float64=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_bound::Float64=1.0, + u_bounds=fill(u_bound, systems[1].n_drives), + u_guess::Union{Matrix{Float64},Nothing}=nothing, + du_bound::Float64=Inf, + du_bounds=fill(du_bound, systems[1].n_drives), + ddu_bound::Float64=1.0, + ddu_bounds=fill(ddu_bound, systems[1].n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt), + Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt), 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, + R::Float64=1e-2, + R_u::Union{Float64,Vector{Float64}}=R, + R_du::Union{Float64,Vector{Float64}}=R, + R_ddu::Union{Float64,Vector{Float64}}=R, state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), @@ -60,21 +94,21 @@ function QuantumStateSamplingProblem( T, Δt, sys.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_names=names, 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, bound_state=piccolo_options.bound_state, - a_guess=a_guess, + u_guess=u_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), timestep=timestep_name) + traj = merge(trajs, merge_names=(u=1, du=1, ddu=1, Δt=1), timestep=timestep_name) end control_names = [ @@ -83,9 +117,9 @@ function QuantumStateSamplingProblem( ] # Objective - J = QuadraticRegularizer(control_names[1], traj, R_a) - J += QuadraticRegularizer(control_names[2], traj, R_da) - J += QuadraticRegularizer(control_names[3], traj, R_dda) + J = QuadraticRegularizer(control_names[1], traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) for (weight, names) in zip(system_weights, state_names) for name in names @@ -154,15 +188,15 @@ end @testitem "Sample systems with single initial, target" begin using PiccoloQuantumObjects - T = 50 + N = 50 Δt = 0.2 - sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) - sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) + sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) + sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) ψ_init = Vector{ComplexF64}([1.0, 0.0]) ψ_target = Vector{ComplexF64}([0.0, 1.0]) prob = QuantumStateSamplingProblem( - [sys1, sys2], ψ_init, ψ_target, T, Δt; + [sys1, sys2], ψ_init, ψ_target, N, Δt; piccolo_options=PiccoloOptions(verbose=false) ) @@ -187,17 +221,17 @@ end @testitem "Sample systems with multiple initial, target" begin using PiccoloQuantumObjects - T = 50 + N = 50 Δt = 0.2 - sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) - sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]]) + sys1 = QuantumSystem(0.3 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) + sys2 = QuantumSystem(-0.3 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) # 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; + [sys1, sys2], ψ_inits, ψ_targets, N, Δt; piccolo_options=PiccoloOptions(verbose=false) ) @@ -219,4 +253,5 @@ end end end -# TODO: Test that a_guess can be used +# TODO: Test that u_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 519eda38..c2e629c5 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -1,14 +1,14 @@ export QuantumStateSmoothPulseProblem """ - QuantumStateSmoothPulseProblem(system, ψ_inits, ψ_goals, T, Δt; kwargs...) - QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, T, Δt; kwargs...) + QuantumStateSmoothPulseProblem(system, ψ_inits, ψ_goals, N, Δt; kwargs...) + QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, N, Δ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)` +Create a quantum state smooth pulse problem. The goal is to find a control pulse `u(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. +`ψ_goals` using `N` knot points of size `Δt`. This problem also controls the first and +second derivatives of the control pulse, `u̇(t)` and `ü(t)`, to ensure smoothness. # Arguments - `system::AbstractQuantumSystem`: The quantum system. @@ -22,30 +22,32 @@ or - `ψ_init::AbstractVector{<:ComplexF64}`: The initial state. - `ψ_goal::AbstractVector{<:ComplexF64}`: The target state. with -- `T::Int`: The number of timesteps. +- `N::Int`: The number of knot points. - `Δt::Float64`: The timestep size. # Keyword Arguments +- `ket_integrator=KetIntegrator`: the integrator to use for state dynamics +- `time_dependent_integrator=false`: whether to use time-dependent integrator - `state_name::Symbol=:ψ̃`: The name of the state variable. -- `control_name::Symbol=:a`: The name of the control variable. +- `control_name::Symbol=:u`: The name of the control variable. - `timestep_name::Symbol=:Δt`: The name of the timestep variable. - `init_trajectory::Union{NamedTrajectory, Nothing}=nothing`: The initial trajectory. -- `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 pulse. -- `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=fill(da_bound, length(system.G_drives))`: The bounds on the first derivative of the control pulse. -- `dda_bound::Float64=1.0`: The bound on the second derivative of the control pulse. -- `dda_bounds=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. -- `Δt_max::Float64=1.5 * Δt`: The maximum timestep size. -- `drive_derivative_σ::Float64=0.01`: The standard deviation of the drive derivative random initialization. +- `u_bound::Float64=1.0`: The bound on the control pulse. +- `u_bounds=fill(u_bound, sys.n_drives)`: The bounds on the control pulse. +- `u_guess::Union{AbstractMatrix{Float64}, Nothing}=nothing`: The initial guess for the control pulse. +- `du_bound::Float64=Inf`: The bound on the first derivative of the control pulse. +- `du_bounds=fill(du_bound, sys.n_drives)`: The bounds on the first derivative of the control pulse. +- `ddu_bound::Float64=1.0`: The bound on the second derivative of the control pulse. +- `ddu_bounds=fill(ddu_bound, sys.n_drives)`: The bounds on the second derivative of the control pulse. +- `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt)`: The minimum timestep size. +- `Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt)`: The maximum timestep size. - `Q::Float64=100.0`: The weight on the state objective. - `R=1e-2`: The weight on the control pulse and its derivatives. -- `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. +- `R_u::Union{Float64, Vector{Float64}}=R`: The weight on the control pulse. +- `R_du::Union{Float64, Vector{Float64}}=R`: The weight on the first derivative of the control pulse. +- `R_ddu::Union{Float64, Vector{Float64}}=R`: The weight on the second derivative of the control pulse. +- `state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing`: Indices of leakage states. - `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: The constraints. - `piccolo_options::PiccoloOptions=PiccoloOptions()`: The Piccolo options. """ @@ -55,27 +57,28 @@ function QuantumStateSmoothPulseProblem( sys::AbstractQuantumSystem, ψ_inits::Vector{<:AbstractVector{<:ComplexF64}}, ψ_goals::Vector{<:AbstractVector{<:ComplexF64}}, - T::Int, + N::Int, Δt::Union{Float64, <:AbstractVector{Float64}}; ket_integrator=KetIntegrator, + time_dependent_integrator=false, state_name::Symbol=:ψ̃, - control_name::Symbol=:a, + control_name::Symbol=:u, timestep_name::Symbol=:Δt, init_trajectory::Union{NamedTrajectory, Nothing}=nothing, - a_bound::Float64=1.0, - a_bounds=fill(a_bound, sys.n_drives), - a_guess::Union{AbstractMatrix{Float64}, Nothing}=nothing, - da_bound::Float64=Inf, - da_bounds=fill(da_bound, sys.n_drives), - dda_bound::Float64=1.0, - dda_bounds=fill(dda_bound, sys.n_drives), - Δt_min::Float64=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_bound::Float64=1.0, + u_bounds=fill(u_bound, sys.n_drives), + u_guess::Union{AbstractMatrix{Float64}, Nothing}=nothing, + du_bound::Float64=Inf, + du_bounds=fill(du_bound, sys.n_drives), + ddu_bound::Float64=1.0, + ddu_bounds=fill(ddu_bound, sys.n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt), + Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt), 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, + R_u::Union{Float64, Vector{Float64}}=R, + R_du::Union{Float64, Vector{Float64}}=R, + R_ddu::Union{Float64, Vector{Float64}}=R, state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), @@ -95,19 +98,20 @@ function QuantumStateSmoothPulseProblem( traj = initialize_trajectory( ψ_goals, ψ_inits, - T, + N, Δt, sys.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_name=state_name, control_name=control_name, timestep_name=timestep_name, zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, Δt_bounds=(Δt_min, Δt_max), bound_state=piccolo_options.bound_state, - a_guess=a_guess, + u_guess=u_guess, system=sys, rollout_integrator=piccolo_options.rollout_integrator, + store_times=time_dependent_integrator ) end @@ -123,9 +127,9 @@ function QuantumStateSmoothPulseProblem( ] # Objective - J = QuadraticRegularizer(control_names[1], traj, R_a) - J += QuadraticRegularizer(control_names[2], traj, R_da) - J += QuadraticRegularizer(control_names[3], traj, R_dda) + J = QuadraticRegularizer(control_names[1], traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) for name ∈ state_names J += KetInfidelityObjective(name, traj; Q=Q) @@ -146,10 +150,8 @@ function QuantumStateSmoothPulseProblem( state_integrators = [] for name ∈ state_names - push!( - state_integrators, - ket_integrator(sys, traj, name, control_name) - ) + integrator_ = ket_integrator(sys, traj, name, control_name) + push!(state_integrators, integrator_) end integrators = [ @@ -191,16 +193,16 @@ end @testitem "Test quantum state smooth pulse" begin using PiccoloQuantumObjects - T = 51 + N = 51 Δt = 0.2 - sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) ψ_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; + sys, ψ_init, ψ_target, N, Δt; piccolo_options=PiccoloOptions(verbose=false) ) initial = rollout_fidelity(prob.trajectory, sys) @@ -212,14 +214,14 @@ end @testitem "Test multiple quantum states smooth pulse" begin using PiccoloQuantumObjects - T = 50 + N = 50 Δt = 0.2 - sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) ψ_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; + sys, ψ_inits, ψ_targets, N, Δt; piccolo_options=PiccoloOptions(verbose=false) ) initial = rollout_fidelity(prob.trajectory, sys) diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 56649272..b4cc0619 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -2,12 +2,42 @@ export UnitaryFreePhaseProblem """ - UnitaryFreePhaseProblem(system::AbstractQuantumSystem, goal::Function, T, Δt; kwargs...) - + UnitaryFreePhaseProblem(system::AbstractQuantumSystem, goal::Function, T::Int, Δt::Union{Float64, AbstractVector{Float64}}; kwargs...) Construct a `DirectTrajOptProblem` for a free-time unitary gate problem with smooth control pulses enforced by constraining the second derivative of the pulse trajectory. The problem follows the same structure as `UnitarySmoothPulseProblem`, but allows for free global phases on the goal unitary, via cosines and sines parameterizing phase variables. The `goal` function should accept a vector of global phases `[cos(θ); sin(θ)]` and return an `AbstractPiccoloOperator`. + +# Arguments +- `system::AbstractQuantumSystem`: the system to be controlled +- `goal::Function`: function that takes phase vector and returns target unitary +- `T::Int`: the number of knot points +- `Δt::Union{Float64, AbstractVector{Float64}}`: the time step size or vector of time steps + +# Keyword Arguments +- `unitary_integrator=UnitaryIntegrator`: the integrator to use for unitary dynamics +- `state_name::Symbol = :Ũ⃗`: the name of the state +- `control_name::Symbol = :u`: the name of the control +- `timestep_name::Symbol = :Δt`: the name of the timestep +- `phase_name::Symbol = :θ`: the name of the phase variable +- `init_phases::Union{AbstractVector{Float64}, Nothing}=nothing`: initial phase values +- `init_trajectory::Union{NamedTrajectory, Nothing}=nothing`: an initial trajectory to use +- `u_guess::Union{Matrix{Float64}, Nothing}=nothing`: an initial guess for the control pulses +- `u_bound::Float64=1.0`: the bound on the control pulse +- `u_bounds=fill(u_bound, system.n_drives)`: the bounds on the control pulses +- `du_bound::Float64=Inf`: the bound on the control pulse derivative +- `du_bounds::Vector{Float64}=fill(du_bound, system.n_drives)`: the bounds on the control pulse derivatives +- `ddu_bound::Float64=1.0`: the bound on the control pulse second derivative +- `ddu_bounds::Vector{Float64}=fill(ddu_bound, system.n_drives)`: the bounds on the control pulse second derivatives +- `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt)`: the minimum time step size +- `Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt)`: the maximum time step size +- `Q::Float64=100.0`: the weight on the infidelity objective +- `R=1e-2`: the weight on the regularization terms +- `R_u::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulses +- `R_du::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse derivatives +- `R_ddu::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse second derivatives +- `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: the constraints to enforce +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: the options for the Piccolo solver """ function UnitaryFreePhaseProblem( system::AbstractQuantumSystem, @@ -15,26 +45,26 @@ function UnitaryFreePhaseProblem( T::Int, Δt::Union{Float64, AbstractVector{Float64}}; unitary_integrator=UnitaryIntegrator, - state_name::Symbol = :Ũ⃗, - control_name::Symbol = :a, + state_name::Symbol = :Ũ⃗, + control_name::Symbol = :u, timestep_name::Symbol = :Δt, phase_name::Symbol = :θ, init_phases::Union{AbstractVector{Float64}, Nothing}=nothing, 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), - 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=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_guess::Union{Matrix{Float64}, Nothing}=nothing, + u_bound::Float64=1.0, + u_bounds=fill(u_bound, system.n_drives), + du_bound::Float64=Inf, + du_bounds::Vector{Float64}=fill(du_bound, system.n_drives), + ddu_bound::Float64=1.0, + ddu_bounds::Vector{Float64}=fill(ddu_bound, system.n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt), + Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt), 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, + R_u::Union{Float64, Vector{Float64}}=R, + R_du::Union{Float64, Vector{Float64}}=R, + R_ddu::Union{Float64, Vector{Float64}}=R, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) @@ -68,7 +98,7 @@ function UnitaryFreePhaseProblem( T, Δt, system.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_name=state_name, control_name=control_name, timestep_name=timestep_name, @@ -76,7 +106,7 @@ function UnitaryFreePhaseProblem( 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, + u_guess=u_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, verbose=piccolo_options.verbose, @@ -92,9 +122,9 @@ function UnitaryFreePhaseProblem( if endswith(string(name), string(control_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) + J += QuadraticRegularizer(control_names[1], traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) # Optional Piccolo constraints and objectives J += apply_piccolo_options!( @@ -131,9 +161,9 @@ end @testitem "UnitaryFreePhaseProblem: basic construction" begin using PiccoloQuantumObjects - sys = QuantumSystem(0.3 * PAULIS.X, [PAULIS.Y]) + sys = QuantumSystem(0.3 * PAULIS.X, [PAULIS.Y], 10.0, [1.0]) U_goal = GATES.Z - T = 51 + N = 51 Δt = 0.2 function virtual_z(z::AbstractVector{<:Real}) @@ -146,7 +176,7 @@ end initial_phases = [pi/3] prob = UnitaryFreePhaseProblem( - sys, virtual_z, T, Δt, + sys, virtual_z, N, Δt, init_phases=initial_phases, piccolo_options=PiccoloOptions(verbose=false), phase_name=:ϕ, diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index dea11aed..eda5f60f 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -157,13 +157,13 @@ end H_drift = 0.1PAULIS[:Z] H_drives = [PAULIS[:X], PAULIS[:Y]] U_goal = GATES[:H] - T = 51 + N = 51 Δt = 0.2 - sys = QuantumSystem(H_drift, H_drives) + sys = QuantumSystem(H_drift, H_drives, 10.0, [1.0, 1.0]) prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, Δt_min=Δt * 0.01, + sys, U_goal, N, Δt_min=Δt * 0.01, piccolo_options=PiccoloOptions(verbose=false) ) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index f7e4d12f..5be46155 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -17,27 +17,28 @@ robust solution by including multiple systems reflecting the problem uncertainty - `Δt::Union{Float64, Vector{Float64}}`: The time step value or vector of time steps. # Keyword Arguments -- `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. +- `unitary_integrator=UnitaryIntegrator`: The integrator to use for unitary dynamics. +- `system_weights=fill(1.0, length(systems))`: The weights for each system. - `init_trajectory::Union{NamedTrajectory, Nothing} = nothing`: The initial trajectory. -- `state_name::Symbol = :Ũ⃗`: The name of the state variable. -- `control_name::Symbol = :a`: The name of the control variable. +- `piccolo_options::PiccoloOptions = PiccoloOptions()`: The Piccolo options. +- `state_name::Symbol = :Ũ⃗`: The name of the state variable. +- `control_name::Symbol = :u`: The name of the control variable. - `timestep_name::Symbol = :Δt`: The name of the timestep variable. - `constraints::Vector{<:AbstractConstraint} = AbstractConstraint[]`: The constraints. -- `a_bound::Float64 = 1.0`: The bound for the control amplitudes. -- `a_bounds = fill(a_bound, length(systems[1].G_drives))`: The bounds for the control amplitudes. -- `a_guess::Union{Matrix{Float64}, Nothing} = nothing`: The initial guess for the control amplitudes. -- `da_bound::Float64 = Inf`: The bound for the control first derivatives. -- `da_bounds = fill(da_bound, length(systems[1].G_drives))`: The bounds for the control first derivatives. -- `dda_bound::Float64 = 1.0`: The bound for the control second derivatives. -- `dda_bounds = fill(dda_bound, length(systems[1].G_drives))`: The bounds for the control second derivatives. -- `Δt_min::Float64 = 0.5 * Δt`: The minimum time step size. -- `Δt_max::Float64 = 1.5 * Δt`: The maximum time step size. +- `u_bound::Float64 = 1.0`: The bound for the control amplitudes. +- `u_bounds = fill(u_bound, systems[1].n_drives)`: The bounds for the control amplitudes. +- `u_guess::Union{Matrix{Float64}, Nothing} = nothing`: The initial guess for the control amplitudes. +- `du_bound::Float64 = Inf`: The bound for the control first derivatives. +- `du_bounds = fill(du_bound, systems[1].n_drives)`: The bounds for the control first derivatives. +- `ddu_bound::Float64 = 1.0`: The bound for the control second derivatives. +- `ddu_bounds = fill(ddu_bound, systems[1].n_drives)`: The bounds for the control second derivatives. +- `Δt_min::Float64 = Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt)`: The minimum time step size. +- `Δt_max::Float64 = Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt)`: The maximum time step size. - `Q::Float64 = 100.0`: The fidelity weight. - `R::Float64 = 1e-2`: The regularization weight. -- `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. +- `R_u::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control amplitudes. +- `R_du::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control first derivatives. +- `R_ddu::Union{Float64, Vector{Float64}} = R`: The regularization weight for the control second derivatives. - `piccolo_options::PiccoloOptions = PiccoloOptions()`: The Piccolo options. """ @@ -48,26 +49,26 @@ function UnitarySamplingProblem( Δt::Union{Float64,Vector{Float64}}; unitary_integrator=UnitaryIntegrator, system_weights=fill(1.0, length(systems)), - init_trajectory::Union{NamedTrajectory,Nothing}=nothing, - piccolo_options::PiccoloOptions=PiccoloOptions(), - state_name::Symbol=:Ũ⃗, - control_name::Symbol=:a, + init_trajectory::Union{NamedTrajectory, Nothing} = nothing, + piccolo_options::PiccoloOptions = PiccoloOptions(), + state_name::Symbol = :Ũ⃗, + control_name::Symbol=:u, timestep_name::Symbol=:Δt, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - a_bound::Float64=1.0, - a_bounds=fill(a_bound, systems[1].n_drives), - a_guess::Union{Matrix{Float64},Nothing}=nothing, - da_bound::Float64=Inf, - da_bounds=fill(da_bound, systems[1].n_drives), - dda_bound::Float64=1.0, - dda_bounds=fill(dda_bound, systems[1].n_drives), - Δt_min::Float64=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_bound::Float64=1.0, + u_bounds=fill(u_bound, systems[1].n_drives), + u_guess::Union{Matrix{Float64}, Nothing}=nothing, + du_bound::Float64=Inf, + du_bounds=fill(du_bound, systems[1].n_drives), + ddu_bound::Float64=1.0, + ddu_bounds=fill(ddu_bound, systems[1].n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt), + Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt), 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, + R::Float64=1e-2, + R_u::Union{Float64,Vector{Float64}}=R, + R_du::Union{Float64,Vector{Float64}}=R, + R_ddu::Union{Float64,Vector{Float64}}=R, kwargs... ) @assert length(systems) == length(operators) @@ -92,14 +93,14 @@ function UnitarySamplingProblem( T, Δt, sys.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_name=st, control_name=control_name, timestep_name=timestep_name, Δt_bounds=(Δt_min, Δt_max), geodesic=piccolo_options.geodesic, bound_state=piccolo_options.bound_state, - a_guess=a_guess, + u_guess=u_guess, system=sys, rollout_integrator=piccolo_options.rollout_integrator, verbose=false # loop @@ -107,7 +108,7 @@ function UnitarySamplingProblem( end traj = merge( - trajs, merge_names=(a=1, da=1, dda=1, Δt=1), timestep=timestep_name + trajs, merge_names=(u=1, du=1, ddu=1, Δt=1), timestep=timestep_name ) end @@ -117,9 +118,9 @@ function UnitarySamplingProblem( ] # Objective - J = QuadraticRegularizer(control_name, traj, R_a) - J += QuadraticRegularizer(control_names[2], traj, R_da) - J += QuadraticRegularizer(control_names[3], traj, R_dda) + J = QuadraticRegularizer(control_name, traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) for (weight, op, name) in zip(system_weights, operators, state_names) J += UnitaryInfidelityObjective(op, name, traj; Q=weight * Q) @@ -171,21 +172,21 @@ end @testitem "Sample robustness test" begin using PiccoloQuantumObjects - T = 50 + N = 50 Δt = 0.2 - timesteps = fill(Δt, T) + timesteps = fill(Δt, N) operator = GATES[:H] - systems(ζ) = QuantumSystem(ζ * GATES[:Z], [GATES[:X], GATES[:Y]]) + systems(ζ) = QuantumSystem(ζ * GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) samples = [0.0, 0.1] prob = UnitarySamplingProblem( - [systems(x) for x in samples], operator, T, Δt, + [systems(x) for x in samples], operator, N, Δt, piccolo_options=PiccoloOptions(verbose=false) ) solve!(prob, max_iter=100, print_level=1, verbose=false) base_prob = UnitarySmoothPulseProblem( - systems(samples[1]), operator, T, Δt, + systems(samples[1]), operator, N, piccolo_options=PiccoloOptions(verbose=false) ) solve!(base_prob, max_iter=100, verbose=false, print_level=1) diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 0ed24635..dc8de8c4 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -2,20 +2,20 @@ export UnitarySmoothPulseProblem @doc raw""" - UnitarySmoothPulseProblem(system::AbstractQuantumSystem, operator::AbstractPiccoloOperator, T::Int, Δt::Float64; kwargs...) - UnitarySmoothPulseProblem(H_drift, H_drives, operator, T, Δt; kwargs...) + UnitarySmoothPulseProblem(system::AbstractQuantumSystem, operator::AbstractPiccoloOperator, N::Int; kwargs...) + UnitarySmoothPulseProblem(H_drift, H_drives, operator, N; kwargs...) 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} -\underset{\vec{\tilde{U}}, a, \dot{a}, \ddot{a}, \Delta t}{\text{minimize}} & \quad -Q \cdot \ell\qty(\vec{\tilde{U}}_T, \vec{\tilde{U}}_{\text{goal}}) + \frac{1}{2} \sum_t \qty(R_a a_t^2 + R_{\dot{a}} \dot{a}_t^2 + R_{\ddot{a}} \ddot{a}_t^2) \\ -\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 \dot{a}_{t+1} - \dot{a}_t - \ddot{a}_t \Delta t_t = 0 \\ -& \quad |a_t| \leq a_{\text{bound}} \\ -& \quad |\ddot{a}_t| \leq \ddot{a}_{\text{bound}} \\ +\underset{\vec{\tilde{U}}, u, \dot{u}, \ddot{u}, \Delta t}{\text{minimize}} & \quad +Q \cdot \ell\qty(\vec{\tilde{U}}_N, \vec{\tilde{U}}_{\text{goal}}) + \frac{1}{2} \sum_t \qty(R_u u_t^2 + R_{\dot{u}} \dot{u}_t^2 + R_{\ddot{u}} \ddot{u}_t^2) \\ +\text{ subject to } & \quad \vb{P}^{(n)}\qty(\vec{\tilde{U}}_{t+1}, \vec{\tilde{U}}_t, u_t, \Delta t_t) = 0 \\ +& \quad u_{t+1} - u_t - \dot{u}_t \Delta t_t = 0 \\ +& \quad \dot{u}_{t+1} - \dot{u}_t - \ddot{u}_t \Delta t_t = 0 \\ +& \quad |u_t| \leq u_{\text{bound}} \\ +& \quad |\ddot{u}_t| \leq \ddot{u}_{\text{bound}} \\ & \quad \Delta t_{\text{min}} \leq \Delta t_t \leq \Delta t_{\text{max}} \\ \end{aligned} ``` @@ -23,11 +23,11 @@ Q \cdot \ell\qty(\vec{\tilde{U}}_T, \vec{\tilde{U}}_{\text{goal}}) + \frac{1}{2} 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)} } +\ell\qty(\vec{\tilde{U}}_N, \vec{\tilde{U}}_{\text{goal}}) = +\abs{1 - \frac{1}{N} \abs{ \tr \qty(U_{\text{goal}}, U_N)} } ``` -is the *infidelity* objective function, $Q$ is a weight, $R_a$, $R_{\dot{a}}$, and $R_{\ddot{a}}$ are weights on the regularization terms, and $\vb{P}^{(n)}$ is the $n$th-order Pade integrator. +is the *infidelity* objective function, $Q$ is a weight, $R_u$, $R_{\dot{u}}$, and $R_{\ddot{u}}$ are weights on the regularization terms, and $\vb{P}^{(n)}$ is the $n$th-order Pade integrator. # Arguments @@ -37,30 +37,28 @@ or - `H_drives::Vector{<:AbstractMatrix{<:Number}}`: the control hamiltonians with - `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 +- `N::Int`: the number of knot points # Keyword Arguments -- `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 +- `unitary_integrator=UnitaryIntegrator`: the integrator to use for unitary dynamics +- `state_name::Symbol = :Ũ⃗`: the name of the state +- `control_name::Symbol = :u`: 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=fill(a_bound, length(system.G_drives))`: the bounds on the control pulses, one for each drive -- `da_bound::Float64=Inf`: 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 -- `dda_bound::Float64=1.0`: the bound on the control pulse second derivative -- `dda_bounds=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 -- `Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt)`: the maximum time step size +- `u_guess::Union{Matrix{Float64}, Nothing}=nothing`: an initial guess for the control pulses +- `du_bound::Float64=Inf`: the bound on the control pulse derivative +- `du_bounds=fill(du_bound, system.n_drives)`: the bounds on the control pulse derivatives, one for each drive +- `ddu_bound::Float64=1.0`: the bound on the control pulse second derivative +- `ddu_bounds=fill(ddu_bound, system.n_drives)`: the bounds on the control pulse second derivatives, one for each drive +- `Δt_min::Float64=0.5 * system.T_max / N`: the minimum time step size +- `Δt_max::Float64=2.0 * system.T_max / N`: the maximum time step size - `Q::Float64=100.0`: the weight on the infidelity objective - `R=1e-2`: the weight on the regularization terms -- `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_u::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulses +- `R_du::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse derivatives +- `R_ddu::Union{Float64, Vector{Float64}}=R`: the weight on the regularization term for the control pulse second derivatives - `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: the constraints to enforce +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: the options for the Piccolo solver """ function UnitarySmoothPulseProblem end @@ -68,27 +66,24 @@ function UnitarySmoothPulseProblem end function UnitarySmoothPulseProblem( system::AbstractQuantumSystem, goal::AbstractPiccoloOperator, - T::Int, - Δt::Union{Float64, <:AbstractVector{Float64}}; + N::Int; unitary_integrator=UnitaryIntegrator, - state_name::Symbol = :Ũ⃗, - control_name::Symbol = :a, + state_name::Symbol = :Ũ⃗, + control_name::Symbol = :u, 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), - 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=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_guess::Union{Matrix{Float64}, Nothing}=nothing, + du_bound::Float64=Inf, + du_bounds=fill(du_bound, system.n_drives), + ddu_bound::Float64=1.0, + ddu_bounds=fill(ddu_bound, system.n_drives), + Δt_min::Float64=0.5 * system.T_max / N, + Δt_max::Float64=2.0 * system.T_max / N, 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, + R_u::Union{Float64, Vector{Float64}}=R, + R_du::Union{Float64, Vector{Float64}}=R, + R_ddu::Union{Float64, Vector{Float64}}=R, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) @@ -101,12 +96,16 @@ function UnitarySmoothPulseProblem( if !isnothing(init_trajectory) traj = init_trajectory else + u_bounds = ( + [system.drive_bounds[j][1] for j in 1:length(system.drive_bounds)], + [system.drive_bounds[j][2] for j in 1:length(system.drive_bounds)] + ) traj = initialize_trajectory( goal, - T, - Δt, + N, + system.T_max / N, system.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_name=state_name, control_name=control_name, timestep_name=timestep_name, @@ -114,7 +113,7 @@ function UnitarySmoothPulseProblem( 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, + u_guess=u_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, verbose=piccolo_options.verbose @@ -129,9 +128,9 @@ function UnitarySmoothPulseProblem( if endswith(string(name), string(control_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) + J += QuadraticRegularizer(control_names[1], traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) # Optional Piccolo constraints and objectives J += apply_piccolo_options!( @@ -148,6 +147,14 @@ function UnitarySmoothPulseProblem( DerivativeIntegrator(traj, control_names[2], control_names[3]), ] + # TODO: see if using a linear constraint here is the right choice + # derivative_constraints = [ + # DerivativeConstraint(control_name, control_names[2], traj), + # DerivativeConstraint(control_names[2], control_names[3], traj), + # ] + + # constraints = vcat(constraints, derivative_constraints) + return DirectTrajOptProblem( traj, J, @@ -171,14 +178,13 @@ end @testitem "Hadamard gate improvement" begin using PiccoloQuantumObjects - sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) + sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) U_goal = GATES[:H] - T = 51 - Δt = 0.2 + N = 51 prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt; - da_bound=1.0, + sys, U_goal, N; + du_bound=1.0, piccolo_options=PiccoloOptions(verbose=false) ) @@ -190,17 +196,16 @@ end @testitem "Bound states and control norm constraint" begin using PiccoloQuantumObjects - sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) + sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]], 10.0, [1.0, 1.0]) U_goal = GATES[:H] - T = 51 - Δt = 0.2 + N = 51 prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, + sys, U_goal, N; piccolo_options=PiccoloOptions( verbose=false, bound_state=true, - complex_control_norm_constraint_name=:a + complex_control_norm_constraint_name=:u ) ) @@ -213,14 +218,13 @@ end using PiccoloQuantumObjects a = annihilate(3) - sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) + sys = QuantumSystem([(a + a')/2, (a - a')/(2im)], 10.0, [1.0, 1.0]) U_goal = EmbeddedOperator(GATES[:H], sys) - T = 51 - Δt = 0.2 + N = 51 @testset "EmbeddedOperator: solve gate" begin prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, + sys, U_goal, N; piccolo_options=PiccoloOptions(verbose=false) ) @@ -231,8 +235,8 @@ end @testset "EmbeddedOperator: leakage constraint" begin prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt; - da_bound=1.0, + sys, U_goal, N; + du_bound=1.0, piccolo_options=PiccoloOptions( leakage_constraint=true, leakage_constraint_value=5e-2, diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl index 966f81b7..189c0044 100644 --- a/src/problem_templates/unitary_variational_problem.jl +++ b/src/problem_templates/unitary_variational_problem.jl @@ -23,29 +23,32 @@ Constructs a unitary variational problem for optimizing quantum control trajecto # Keyword Arguments -- `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: `:Ũ⃗ₐ`). -- `variational_scales::AbstractVector`: Scaling factors for the variational state variables (default: `1.0`). -- `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`: 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`: 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`: 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. +- `robust_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars]`: Times at which robustness to variations in the trajectory is enforced. +- `sensitive_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars]`: Times at which sensitivity to variations in the trajectory is enhanced. +- `variational_integrator=VariationalUnitaryIntegrator`: The integrator used for unitary evolution. +- `variational_scales::AbstractVector{<:Float64}=fill(1.0, length(system.G_vars))`: Scaling factors for the variational state variables. +- `state_name::Symbol = :Ũ⃗`: The name of the state variable in the trajectory. +- `variational_state_name::Symbol = :Ũ⃗ᵥ`: The name of the variational state variable. +- `control_name::Symbol = :u`: The name of the control variable. +- `timestep_name::Symbol = :Δt`: The name of the timestep variable. +- `init_trajectory::Union{NamedTrajectory, Nothing}=nothing`: An optional initial trajectory to start optimization. +- `u_bound::Float64=1.0`: The bound for the control variable `u`. +- `u_bounds=fill(u_bound, system.n_drives)`: Bounds for each control variable. +- `du_bound::Float64=Inf`: The bound for the derivative of the control variable. +- `du_bounds=fill(du_bound, system.n_drives)`: Bounds for each derivative of the control variable. +- `ddu_bound::Float64=1.0`: The bound for the second derivative of the control variable. +- `ddu_bounds=fill(ddu_bound, system.n_drives)`: Bounds for each second derivative of the control variable. +- `Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt)`: Minimum allowed timestep duration. +- `Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δt)`: Maximum allowed timestep duration. +- `Q::Float64=100.0`: Weight for the unitary infidelity objective. +- `Q_s::Float64=1e-2`: Weight for sensitivity objectives. +- `Q_r::Float64=100.0`: Weight for robustness objectives. +- `R=1e-2`: Regularization weight for control variables. +- `R_u::Union{Float64, Vector{Float64}}=R`: Regularization weights for control. +- `R_du::Union{Float64, Vector{Float64}}=R`: Regularization weights for control derivative. +- `R_ddu::Union{Float64, Vector{Float64}}=R`: Regularization weights for control second derivative. +- `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: Additional constraints for the optimization problem. +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: Options for configuring the Piccolo optimization framework. # Returns @@ -70,26 +73,26 @@ function UnitaryVariationalProblem( sensitive_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars], variational_integrator=VariationalUnitaryIntegrator, variational_scales::AbstractVector{<:Float64}=fill(1.0, length(system.G_vars)), - state_name::Symbol = :Ũ⃗, - variational_state_name::Symbol = :Ũ⃗ᵥ, - control_name::Symbol = :a, + state_name::Symbol = :Ũ⃗, + variational_state_name::Symbol = :Ũ⃗ᵥ, + control_name::Symbol = :u, 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=fill(da_bound, system.n_drives), - dda_bound::Float64=1.0, - dda_bounds=fill(dda_bound, system.n_drives), - Δt_min::Float64=0.5 * minimum(Δt), - Δt_max::Float64=2.0 * maximum(Δt), + u_bound::Float64=1.0, + u_bounds=fill(u_bound, system.n_drives), + du_bound::Float64=Inf, + du_bounds=fill(du_bound, system.n_drives), + ddu_bound::Float64=1.0, + ddu_bounds=fill(ddu_bound, system.n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * minimum(Δt), + Δt_max::Float64=Δt isa Float64 ? 2.0 * Δt : 2.0 * maximum(Δ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, + R_u::Union{Float64, Vector{Float64}}=R, + R_du::Union{Float64, Vector{Float64}}=R, + R_ddu::Union{Float64, Vector{Float64}}=R, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) @@ -114,7 +117,7 @@ function UnitaryVariationalProblem( T, Δt, system.n_drives, - (a_bounds, da_bounds, dda_bounds); + (u_bounds, du_bounds, ddu_bounds); state_name=state_name, control_name=control_name, timestep_name=timestep_name, @@ -156,11 +159,11 @@ function UnitaryVariationalProblem( if endswith(string(name), string(control_name)) ] - # objective + # 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) + J += QuadraticRegularizer(control_names[1], traj, R_u) + J += QuadraticRegularizer(control_names[2], traj, R_du) + J += QuadraticRegularizer(control_names[3], traj, R_ddu) # sensitivity for (name, scale, s, r) ∈ zip( @@ -211,30 +214,30 @@ end using LinearAlgebra using PiccoloQuantumObjects - system = QuantumSystem([PAULIS.X, PAULIS.Y]) - varsys = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X] ) - T = 50 + system = QuantumSystem([PAULIS.X, PAULIS.Y], 10.0, [1.0, 1.0]) + varsys = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X], 10.0, [1.0, 1.0]) + N = 50 Δt = 0.2 sense_scale = 8.0 sense_prob = UnitaryVariationalProblem( - varsys, GATES.X, T, Δt, + varsys, GATES.X, N, Δt, variational_scales=[sense_scale], - sensitive_times=[[T]], + sensitive_times=[[N]], 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, + varsys, GATES.X, N, Δt, variational_scales=[rob_scale], - robust_times=[[T]], + robust_times=[[N]], 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) + sense_n = norm(sense_scale * sense_prob.trajectory.Ũ⃗ᵥ1) + rob_n = norm(rob_scale * rob_prob.trajectory.Ũ⃗ᵥ1) @test sense_n > rob_n end \ No newline at end of file diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index dbb23fe2..a341b614 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -29,7 +29,7 @@ function FinalKetFidelityConstraint( ψ̃_name, traj, equality=false, - times=[traj.T] + times=[traj.N] ) end @@ -50,7 +50,7 @@ function FinalUnitaryFidelityConstraint( Ũ⃗_name, traj, equality=false, - times=[traj.T] + times=[traj.N] ) end @@ -73,7 +73,7 @@ function FinalUnitaryFidelityConstraint( θ_names, traj, equality=false, - times=[traj.T] + times=[traj.N] ) end @@ -92,7 +92,7 @@ function LeakageConstraint( indices::AbstractVector{Int}, name::Symbol, traj::NamedTrajectory; - times=1:traj.T, + times=1:traj.N, ) leakage_constraint(x) = abs2.(x[indices]) .- value diff --git a/src/quantum_integrators.jl b/src/quantum_integrators.jl index 8108a1c8..86eb6f4f 100644 --- a/src/quantum_integrators.jl +++ b/src/quantum_integrators.jl @@ -21,28 +21,29 @@ function KetIntegrator( sys::QuantumSystem, traj::NamedTrajectory, ψ̃::Symbol, - a::Symbol + u::Symbol ) - return BilinearIntegrator(sys.G, traj, ψ̃, a) + Ĝ = u_ -> sys.G(u_, 0.0) + return BilinearIntegrator(Ĝ, traj, ψ̃, u) end function UnitaryIntegrator( sys::QuantumSystem, traj::NamedTrajectory, - Ũ⃗::Symbol, - a::Symbol + Ũ⃗::Symbol, + u::Symbol ) - Ĝ = a_ -> I(sys.levels) ⊗ sys.G(a_) - return BilinearIntegrator(Ĝ, traj, Ũ⃗, a) + Ĝ = u_ -> I(sys.levels) ⊗ sys.G(u_, 0.0) + return BilinearIntegrator(Ĝ, traj, Ũ⃗, u) end function DensityMatrixIntegrator( sys::OpenQuantumSystem, traj::NamedTrajectory, ρ̃::Symbol, - a::Symbol + u::Symbol ) - return BilinearIntegrator(sys.𝒢, traj, ρ̃, a) + return BilinearIntegrator(sys.𝒢, traj, ρ̃, u) end # ----------------------------------------------------------------------------- # diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index f7da6c99..5e298aaa 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -162,7 +162,7 @@ function LeakageObjective( indices::AbstractVector{Int}, name::Symbol, traj::NamedTrajectory; - times=1:traj.T, + times=1:traj.N, Qs::AbstractVector{<:Float64}=fill(1.0, length(times)), ) leakage_objective(x) = sum(abs2, x[indices]) / length(indices) diff --git a/src/quantum_system_templates/cats.jl b/src/quantum_system_templates/cats.jl index 1d37a194..eb87ec0f 100644 --- a/src/quantum_system_templates/cats.jl +++ b/src/quantum_system_templates/cats.jl @@ -31,6 +31,8 @@ function CatSystem(; cat_levels::Int=13, buffer_levels::Int=3, prefactor::Real=1, + T_max::Real=1.0, + drive_bounds::Vector{<:Real}=[1.0, 1.0] ) params = Dict( :g2 => prefactor * g2, @@ -62,19 +64,21 @@ function CatSystem(; return OpenQuantumSystem( H_drift, H_drives, - L_dissipators; + T_max, + drive_bounds, + dissipation_operators=L_dissipators, params=params ) end -function get_cat_controls(sys::AbstractQuantumSystem, α::Real, T::Int) +function get_cat_controls(sys::AbstractQuantumSystem, α::Real, N::Int) @assert haskey(sys.params, :g2) "Requires photon transfer coupling between buffer and cat" @assert haskey(sys.params, :χ_aa) "Requires Kerr coupling for cat" buffer_drive = abs2(α) * sys.params[:g2] cat_kerr_correction = (2.0 * abs2(α) + 1.0) * sys.params[:χ_aa] return stack([ - fill(buffer_drive, T), - fill(cat_kerr_correction, T) + fill(buffer_drive, N), + fill(cat_kerr_correction, N) ], dims=1) end diff --git a/src/quantum_system_templates/rydberg.jl b/src/quantum_system_templates/rydberg.jl index 56a64e03..bdef256f 100644 --- a/src/quantum_system_templates/rydberg.jl +++ b/src/quantum_system_templates/rydberg.jl @@ -42,7 +42,9 @@ end local_detune::Bool=false, # If true, include one local detuning pattern. all2all::Bool=true, # If true, include all-to-all interactions. ignore_Y_drive::Bool=false, # If true, ignore the Y drive. (In the experiments, X&Y drives are implemented by Rabi amplitude and its phase.) - )::QuantumSystem + T_max::Float64=10.0, # Maximum evolution time + drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=[1.0, 1.0, 1.0], # Bounds for [Ωx, Ωy, Δ] or [Ωx, Δ] if ignore_Y_drive + ) -> QuantumSystem Returns a `QuantumSystem` object for the Rydberg atom chain in the spin basis |g⟩ = |0⟩ = [1, 0], |r⟩ = |1⟩ = [0, 1]. @@ -60,6 +62,8 @@ H = \sum_i 0.5*\Omega_i(t)\cos(\phi_i(t)) \sigma_i^x - 0.5*\Omega_i(t)\sin(\phi_ - `local_detune`: If true, include one local detuning pattern. - `all2all`: If true, include all-to-all interactions. - `ignore_Y_drive`: If true, ignore the Y drive. (In the experiments, X&Y drives are implemented by Rabi amplitude and its phase.) +- `T_max`: Maximum evolution time. +- `drive_bounds`: Bounds for drive amplitudes. """ function RydbergChainSystem(; N::Int=3, # number of atoms @@ -69,6 +73,8 @@ function RydbergChainSystem(; local_detune::Bool=false, all2all::Bool=true, ignore_Y_drive::Bool=false, + T_max::Float64=10.0, + drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=ignore_Y_drive ? [1.0, 1.0] : [1.0, 1.0, 1.0], ) PAULIS = ( I = ComplexF64[1 0; 0 1], @@ -117,7 +123,9 @@ function RydbergChainSystem(; return QuantumSystem( H_drift, - H_drives + H_drives, + T_max, + drive_bounds ) end @@ -126,5 +134,13 @@ end @testitem "Rydberg system test" begin using PiccoloQuantumObjects - @test RydbergChainSystem(N=3,cutoff_order=2,all2all=false) isa QuantumSystem + sys = RydbergChainSystem(N=3, cutoff_order=2, all2all=false) + @test sys isa QuantumSystem + @test sys.levels == 8 # 2^3 for 3 atoms + @test sys.n_drives == 3 # X, Y, and detuning + + # Test with ignore_Y_drive + sys2 = RydbergChainSystem(N=2, ignore_Y_drive=true) + @test sys2 isa QuantumSystem + @test sys2.n_drives == 2 # X and detuning only end diff --git a/src/quantum_system_templates/transmons.jl b/src/quantum_system_templates/transmons.jl index 7bc15071..a1b28087 100644 --- a/src/quantum_system_templates/transmons.jl +++ b/src/quantum_system_templates/transmons.jl @@ -31,6 +31,8 @@ where `a` is the annihilation operator. - `drives`: Whether to include drives in the Hamiltonian. """ function TransmonSystem(; + T_max::Float64=10.0, + drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}=fill(1.0, 2), ω::Float64=4.0, # GHz δ::Float64=0.2, # GHz levels::Int=3, @@ -39,8 +41,8 @@ function TransmonSystem(; mutiply_by_2π::Bool=true, lab_frame_type::Symbol=:duffing, drives::Bool=true, + kwargs... ) - @assert lab_frame_type ∈ (:duffing, :quartic, :cosine) "lab_frame_type must be one of (:duffing, :quartic, :cosine)" if lab_frame @@ -85,21 +87,11 @@ function TransmonSystem(; H_drives .*= 2π end - params = Dict{Symbol, Any}( - :ω => ω, - :δ => δ, - :levels => levels, - :lab_frame => lab_frame, - :frame_ω => frame_ω, - :mutiply_by_2π => mutiply_by_2π, - :lab_frame_type => lab_frame_type, - :drives => drives, - ) - return QuantumSystem( H_drift, - H_drives; - params=params + H_drives, + T_max, + drive_bounds ) end @@ -196,6 +188,8 @@ end ωs::AbstractVector{Float64}, δs::AbstractVector{Float64}, gs::AbstractMatrix{Float64}; + T_max::Float64=10.0, + drive_bounds::Union{Float64, Vector{<:Union{Tuple{Float64, Float64}, Float64}}}=1.0, levels_per_transmon::Int = 3, subsystem_levels::AbstractVector{Int} = fill(levels_per_transmon, length(ωs)), lab_frame=false, @@ -210,6 +204,8 @@ function MultiTransmonSystem( ωs::AbstractVector{Float64}, δs::AbstractVector{Float64}, gs::AbstractMatrix{Float64}; + T_max::Float64=10.0, + drive_bounds::Union{Float64, Vector{<:Union{Tuple{Float64, Float64}, Float64}}}=1.0, levels_per_transmon::Int = 3, subsystem_levels::AbstractVector{Int} = fill(levels_per_transmon, length(ωs)), lab_frame=false, @@ -221,11 +217,20 @@ function MultiTransmonSystem( @assert length(δs) == n_subsystems @assert size(gs) == (n_subsystems, n_subsystems) + # Convert drive_bounds to vector if scalar + if drive_bounds isa Float64 + drive_bounds_vec = fill(drive_bounds, 2) + else + drive_bounds_vec = drive_bounds + end + systems = QuantumSystem[] for (i, (ω, δ, levels)) ∈ enumerate(zip(ωs, δs, subsystem_levels)) if i ∈ subsystems sysᵢ = TransmonSystem( + T_max=T_max, + drive_bounds=drive_bounds_vec, levels=levels, ω=ω, δ=δ, @@ -251,7 +256,7 @@ function MultiTransmonSystem( levels = prod([sys.levels for sys in systems]) H_drift = sum(c -> c.op, couplings; init=zeros(ComplexF64, levels, levels)) - return CompositeQuantumSystem(H_drift, systems) + return CompositeQuantumSystem(H_drift, systems, T_max, drive_bounds_vec) end # *************************************************************************** # @@ -260,16 +265,12 @@ end using PiccoloQuantumObjects sys = TransmonSystem() @test sys isa QuantumSystem - @test haskey(sys.params, :ω) - @test haskey(sys.params, :δ) - @test sys.params[:levels] == 3 + @test sys.levels == 3 + @test sys.n_drives == 2 sys2 = TransmonSystem(ω=5.0, δ=0.3, levels=4, lab_frame=true, frame_ω=0.0, lab_frame_type=:duffing, drives=false) - @test sys2.params[:ω] == 5.0 - @test sys2.params[:δ] == 0.3 - @test sys2.params[:levels] == 4 - @test sys2.params[:lab_frame] == true - @test sys2.params[:drives] == false + @test sys2.levels == 4 + @test sys2.n_drives == 0 end @testitem "TransmonSystem: lab_frame_type variations" begin @@ -312,7 +313,7 @@ end comp = MultiTransmonSystem(ωs, δs, gs) @test comp isa CompositeQuantumSystem @test length(comp.subsystems) == 2 - @test !iszero(comp.H(zeros(comp.n_drives))) + @test !iszero(comp.H(zeros(comp.n_drives), 0.0)) comp2 = MultiTransmonSystem( ωs, δs, gs; @@ -323,7 +324,7 @@ end ) @test comp2 isa CompositeQuantumSystem @test length(comp2.subsystems) == 1 - @test !isapprox(norm(comp2.H(zeros(comp2.n_drives))), 0.0; atol=1e-12) + @test !isapprox(norm(comp2.H(zeros(comp2.n_drives), 0.0)), 0.0; atol=1e-12) end @testitem "MultiTransmonSystem: edge cases" begin @@ -340,6 +341,6 @@ end @test comp isa CompositeQuantumSystem @test length(comp.subsystems) == 2 # Only one drive - @test comp.subsystems[1].params[:drives] == false - @test comp.subsystems[2].params[:drives] == true + @test comp.subsystems[1].n_drives == 0 + @test comp.subsystems[2].n_drives == 2 end \ No newline at end of file diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 9b953c18..0b676edf 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -144,10 +144,37 @@ linear_interpolation(X::AbstractMatrix, Y::AbstractMatrix, n::Int) = # ============================================================================= # +""" + initialize_unitary_trajectory( + U_init::AbstractMatrix{<:Number}, + U_goal::AbstractPiccoloOperator, + N::Int; + geodesic::Bool=true, + system::Union{AbstractQuantumSystem, Nothing}=nothing + ) + +Generate an initial unitary trajectory from `U_init` to `U_goal`. + +# Arguments +- `U_init::AbstractMatrix{<:Number}`: Initial unitary operator +- `U_goal::AbstractPiccoloOperator`: Target unitary operator +- `N::Int`: Number of time steps + +# Keyword Arguments +- `geodesic::Bool=true`: Use geodesic interpolation (vs. linear interpolation) +- `system::Union{AbstractQuantumSystem, Nothing}=nothing`: System for drift Hamiltonian + +# Returns +- `Matrix{Float64}`: Trajectory of unitaries in iso-vec representation (column per timestep) + +# Notes +- Geodesic interpolation follows the shortest path on the unitary manifold +- If a system is provided, the drift Hamiltonian is used in geodesic computation +""" function initialize_unitary_trajectory( U_init::AbstractMatrix{<:Number}, U_goal::AbstractPiccoloOperator, - T::Int; + N::Int; geodesic::Bool=true, system::Union{AbstractQuantumSystem, Nothing}=nothing ) @@ -157,9 +184,9 @@ function initialize_unitary_trajectory( else H_drift = zeros(size(U_init)) end - Ũ⃗ = unitary_geodesic(U_init, U_goal, T, H_drift=H_drift) + Ũ⃗ = unitary_geodesic(U_init, U_goal, N, H_drift=H_drift) else - Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, T) + Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, N) end return Ũ⃗ end @@ -168,10 +195,38 @@ end # Initial controls # # ----------------------------------------------------------------------------- # +""" + initialize_control_trajectory( + n_drives::Int, + n_derivatives::Int, + N::Int, + bounds::VectorBound, + drive_derivative_σ::Float64 + ) + +Generate random initial control trajectories with derivatives. + +Creates smooth control trajectories by randomly sampling the base controls within bounds +and generating higher derivatives from a Gaussian distribution. + +# Arguments +- `n_drives::Int`: Number of independent control drives +- `n_derivatives::Int`: Number of derivatives to generate (e.g., 2 for u, du, ddu) +- `N::Int`: Number of time steps +- `bounds::VectorBound`: Bounds for the control amplitudes (vector or tuple of bounds) +- `drive_derivative_σ::Float64`: Standard deviation for random derivative initialization + +# Returns +- `Vector{Matrix{Float64}}`: Vector of control matrices [u, du, ddu, ...] + +# Notes +- Base controls are zero at initial and final timesteps for smooth boundary conditions +- Derivatives are sampled from N(0, drive_derivative_σ²) +""" function initialize_control_trajectory( n_drives::Int, n_derivatives::Int, - T::Int, + N::Int, bounds::VectorBound, drive_derivative_σ::Float64, ) @@ -187,18 +242,42 @@ function initialize_control_trajectory( a = hcat([ zeros(n_drives), - vcat([rand(a_dists[i], 1, T - 2) for i = 1:n_drives]...), + vcat([rand(a_dists[i], 1, N - 2) for i = 1:n_drives]...), zeros(n_drives) ]...) push!(controls, a) for _ in 1:n_derivatives - push!(controls, randn(n_drives, T) * drive_derivative_σ) + push!(controls, randn(n_drives, N) * drive_derivative_σ) end return controls end +""" + initialize_control_trajectory( + u::AbstractMatrix, + Δt::AbstractVecOrMat, + n_derivatives::Int + ) + +Generate control derivatives from a provided control trajectory. + +Takes a given control trajectory and computes its time derivatives using finite differences. +Ensures smooth transitions at boundaries to avoid constraint violations. + +# Arguments +- `u::AbstractMatrix`: Control trajectory (n_drives × N) +- `Δt::AbstractVecOrMat`: Time step size(s) +- `n_derivatives::Int`: Number of derivatives to compute + +# Returns +- `Vector{Matrix{Float64}}`: Vector of control matrices [u, du, ddu, ...] + +# Notes +- Uses finite difference approximation for derivatives +- Adjusts penultimate point to ensure smooth final derivative +""" function initialize_control_trajectory( a::AbstractMatrix, Δt::AbstractVecOrMat, @@ -227,32 +306,66 @@ initialize_control_trajectory(a::AbstractMatrix, Δt::Real, n_derivatives::Int) # ----------------------------------------------------------------------------- # """ - initialize_trajectory + initialize_trajectory( + state_data::Vector{<:AbstractMatrix{Float64}}, + state_inits::Vector{<:AbstractVector{Float64}}, + state_goals::Vector{<:AbstractVector{Float64}}, + state_names::AbstractVector{Symbol}, + N::Int, + Δt::Union{Float64, AbstractVecOrMat{<:Float64}}, + n_drives::Int, + control_bounds::Tuple{Vararg{VectorBound}}; + kwargs... + ) +Initialize a trajectory for a quantum control problem with custom state data. -Initialize a trajectory for a control problem. The trajectory is initialized with -data that should be consistently the same type (in this case, Float64). +# Arguments +- `state_data::Vector{<:AbstractMatrix{Float64}}`: Pre-computed state trajectories (one matrix per state) +- `state_inits::Vector{<:AbstractVector{Float64}}`: Initial state values +- `state_goals::Vector{<:AbstractVector{Float64}}`: Target state values +- `state_names::AbstractVector{Symbol}`: Names for each state component +- `N::Int`: Number of time steps +- `Δt::Union{Float64, AbstractVecOrMat{<:Float64}}`: Time step size(s) +- `n_drives::Int`: Number of control drives +- `control_bounds::Tuple{Vararg{VectorBound}}`: Bounds for controls and their derivatives +# Keyword Arguments +- `bound_state::Bool=false`: Whether to bound the state variables +- `control_name::Symbol=:u`: Name for the control variable +- `n_control_derivatives::Int=length(control_bounds) - 1`: Number of control derivatives +- `zero_initial_and_final_derivative::Bool=false`: Enforce zero derivatives at boundaries +- `timestep_name::Symbol=:Δt`: Name for the timestep variable +- `Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt)`: Bounds for the timestep +- `drive_derivative_σ::Float64=0.1`: Standard deviation for random control derivatives +- `u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing`: Initial guess for controls +- `global_component_data::NamedTuple=NamedTuple()`: Additional global trajectory components +- `verbose::Bool=false`: Print detailed initialization information +- `store_times::Bool=false`: Store cumulative time values in the trajectory + +# Returns +- `NamedTrajectory`: Initialized trajectory with states, controls, and timesteps """ function initialize_trajectory( state_data::Vector{<:AbstractMatrix{Float64}}, state_inits::Vector{<:AbstractVector{Float64}}, state_goals::Vector{<:AbstractVector{Float64}}, state_names::AbstractVector{Symbol}, - T::Int, + N::Int, Δt::Union{Float64, AbstractVecOrMat{<:Float64}}, n_drives::Int, control_bounds::Tuple{Vararg{VectorBound}}; bound_state=false, - control_name=:a, + control_name=:u, n_control_derivatives::Int=length(control_bounds) - 1, zero_initial_and_final_derivative=false, timestep_name=:Δt, Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt), drive_derivative_σ::Float64=0.1, - a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, global_component_data::NamedTuple{gname, <:Tuple{Vararg{AbstractVector{<:Real}}}} where gname=(;), verbose=false, + store_times=false, ) @assert length(state_data) == length(state_names) == length(state_inits) == length(state_goals) "state_data, state_names, state_inits, and state_goals must have the same length" @assert length(control_bounds) == n_control_derivatives + 1 "control_bounds must have $n_control_derivatives + 1 elements" @@ -274,12 +387,12 @@ function initialize_trajectory( # Timestep data if Δt isa Real - timestep_data = fill(Δt, 1, T) + timestep_data = fill(Δt, 1, N) elseif Δt isa AbstractVector timestep_data = reshape(Δt, 1, :) else timestep_data = Δt - @assert size(Δt) == (1, T) "Δt must be a Real, AbstractVector, or 1x$(T) AbstractMatrix" + @assert size(Δt) == (1, N) "Δt must be a Real, AbstractVector, or 1x$(N) AbstractMatrix" end timestep = timestep_name @@ -289,6 +402,11 @@ function initialize_trajectory( control_name => zeros(n_drives), ) + if store_times + initial = merge(initial, (; t=[0.0])) + t_data = cumsum(timestep_data, dims=2) + end + final = (; control_name => zeros(n_drives), ) @@ -313,24 +431,30 @@ function initialize_trajectory( end # Trajectory - if isnothing(a_guess) + if isnothing(u_guess) # Randomly sample controls control_data = initialize_control_trajectory( n_drives, n_control_derivatives, - T, + N, bounds[control_name], drive_derivative_σ ) else # Use provided controls and take derivatives - control_data = initialize_control_trajectory(a_guess, Δt, n_control_derivatives) + control_data = initialize_control_trajectory(u_guess, Δt, n_control_derivatives) end names = [state_names..., control_names..., timestep_name] values = [state_data..., control_data..., timestep_data] controls = (control_names[end], timestep_name) + if store_times + names = [names..., :t] + values = [values..., t_data] + controls = (controls..., :t) + end + return NamedTrajectory( (; (names .=> values)...), global_component_data; @@ -344,18 +468,52 @@ function initialize_trajectory( end """ - initialize_trajectory + initialize_trajectory( + U_goal::AbstractPiccoloOperator, + N::Int, + Δt::Union{Real, AbstractVecOrMat{<:Real}}, + n_drives::Int, + control_bounds::Tuple{Vararg{VectorBound}}; + kwargs... + ) -Trajectory initialization of unitaries. +Initialize a trajectory for unitary gate synthesis problems. + +Constructs a trajectory that evolves from an initial unitary (default: identity) to a target +unitary gate. The trajectory can use geodesic interpolation or rollout-based initialization. + +# Arguments +- `U_goal::AbstractPiccoloOperator`: Target unitary operator (can be `EmbeddedOperator`) +- `N::Int`: Number of time steps +- `Δt::Union{Real, AbstractVecOrMat{<:Real}}`: Time step size(s) +- `n_drives::Int`: Number of control drives +- `control_bounds::Tuple{Vararg{VectorBound}}`: Bounds for controls and their derivatives + +# Keyword Arguments +- `state_name::Symbol=:Ũ⃗`: Name for the unitary state variable (iso-vec representation) +- `U_init::AbstractMatrix{<:Number}=I`: Initial unitary operator +- `u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing`: Initial guess for controls +- `system::Union{AbstractQuantumSystem, Nothing}=nothing`: Quantum system for rollout +- `rollout_integrator::Function=expv`: Integrator for unitary dynamics +- `geodesic=true`: Use geodesic interpolation between unitaries +- Additional kwargs passed to the base `initialize_trajectory` method + +# Returns +- `NamedTrajectory`: Initialized trajectory with unitary states, controls, and timesteps + +# Notes +- If `u_guess` is provided, the trajectory is computed via rollout using the quantum system +- If `u_guess` is `nothing`, geodesic interpolation is used (requires `geodesic=true`) +- The unitary is stored in iso-vec representation for efficient optimization """ function initialize_trajectory( U_goal::AbstractPiccoloOperator, - T::Int, + N::Int, Δt::Union{Real, AbstractVecOrMat{<:Real}}, args...; - state_name::Symbol=:Ũ⃗, + state_name::Symbol=:Ũ⃗, U_init::AbstractMatrix{<:Number}=Matrix{ComplexF64}(I(size(U_goal, 1))), - a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, system::Union{AbstractQuantumSystem, Nothing}=nothing, rollout_integrator::Function=expv, geodesic=true, @@ -365,13 +523,13 @@ function initialize_trajectory( if Δt isa AbstractMatrix timesteps = vec(Δt) elseif Δt isa Float64 - timesteps = fill(Δt, T) + timesteps = fill(Δt, N) else timesteps = Δt end # Initial state and goal - Ũ⃗_init = operator_to_iso_vec(U_init) + Ũ⃗_init = operator_to_iso_vec(U_init) if U_goal isa EmbeddedOperator Ũ⃗_goal = operator_to_iso_vec(U_goal.operator) @@ -380,28 +538,28 @@ function initialize_trajectory( end # Construct state data - if isnothing(a_guess) - Ũ⃗_traj = initialize_unitary_trajectory( + if isnothing(u_guess) + Ũ⃗_traj = initialize_unitary_trajectory( U_init, U_goal, - T; + N; geodesic=geodesic, system=system ) else - @assert !isnothing(system) "System must be provided if a_guess is provided." - Ũ⃗_traj = unitary_rollout(Ũ⃗_init, a_guess, timesteps, system; integrator=rollout_integrator) + @assert !isnothing(system) "System must be provided if u_guess is provided." + Ũ⃗_traj = unitary_rollout(Ũ⃗_init, u_guess, timesteps, system; integrator=rollout_integrator) end return initialize_trajectory( - [Ũ⃗_traj], - [Ũ⃗_init], - [Ũ⃗_goal], + [Ũ⃗_traj], + [Ũ⃗_init], + [Ũ⃗_goal], [state_name], - T, + N, Δt, args...; - a_guess=a_guess, + u_guess=u_guess, kwargs... ) end @@ -409,21 +567,57 @@ end """ - initialize_trajectory + initialize_trajectory( + ψ_goals::AbstractVector{<:AbstractVector{ComplexF64}}, + ψ_inits::AbstractVector{<:AbstractVector{ComplexF64}}, + N::Int, + Δt::Union{Real, AbstractVector{<:Real}}, + n_drives::Int, + control_bounds::Tuple{Vararg{VectorBound}}; + kwargs... + ) + +Initialize a trajectory for quantum state transfer problems. -Trajectory initialization of quantum states. +Constructs a trajectory that evolves one or more quantum states from initial states to target +states. Supports multiple simultaneous state trajectories with shared controls. + +# Arguments +- `ψ_goals::AbstractVector{<:AbstractVector{ComplexF64}}`: Target quantum state(s) +- `ψ_inits::AbstractVector{<:AbstractVector{ComplexF64}}`: Initial quantum state(s) +- `N::Int`: Number of time steps +- `Δt::Union{Real, AbstractVector{<:Real}}`: Time step size(s) +- `n_drives::Int`: Number of control drives +- `control_bounds::Tuple{Vararg{VectorBound}}`: Bounds for controls and their derivatives + +# Keyword Arguments +- `state_name::Symbol=:ψ̃`: Base name for state variables (iso representation) +- `state_names::AbstractVector{<:Symbol}`: Explicit names for each state (auto-generated if not provided) +- `u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing`: Initial guess for controls +- `system::Union{AbstractQuantumSystem, Nothing}=nothing`: Quantum system for rollout +- `rollout_integrator::Function=expv`: Integrator for state dynamics +- Additional kwargs passed to the base `initialize_trajectory` method + +# Returns +- `NamedTrajectory`: Initialized trajectory with quantum states, controls, and timesteps + +# Notes +- States are stored in iso representation (real-valued vectors) for optimization +- If `u_guess` is provided, trajectories are computed via rollout +- If `u_guess` is `nothing`, states are linearly interpolated +- Multiple states share the same control trajectory """ function initialize_trajectory( ψ_goals::AbstractVector{<:AbstractVector{ComplexF64}}, ψ_inits::AbstractVector{<:AbstractVector{ComplexF64}}, - T::Int, + N::Int, Δt::Union{Real, AbstractVector{<:Real}}, args...; state_name=:ψ̃, state_names::AbstractVector{<:Symbol}=length(ψ_goals) == 1 ? [state_name] : [Symbol(string(state_name) * "$i") for i = 1:length(ψ_goals)], - a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, system::Union{AbstractQuantumSystem, Nothing}=nothing, rollout_integrator::Function=expv, kwargs... @@ -438,16 +632,16 @@ function initialize_trajectory( if Δt isa AbstractMatrix timesteps = vec(Δt) elseif Δt isa Float64 - timesteps = fill(Δt, T) + timesteps = fill(Δt, N) else timesteps = Δt end # Construct state data ψ̃_trajs = Matrix{Float64}[] - if isnothing(a_guess) + if isnothing(u_guess) for (ψ̃_init, ψ̃_goal) ∈ zip(ψ̃_inits, ψ̃_goals) - ψ̃_traj = linear_interpolation(ψ̃_init, ψ̃_goal, T) + ψ̃_traj = linear_interpolation(ψ̃_init, ψ̃_goal, N) push!(ψ̃_trajs, ψ̃_traj) end if system isa AbstractVector @@ -455,7 +649,7 @@ function initialize_trajectory( end else for ψ̃_init ∈ ψ̃_inits - ψ̃_traj = rollout(ψ̃_init, a_guess, timesteps, system; integrator=rollout_integrator) + ψ̃_traj = rollout(ψ̃_init, u_guess, timesteps, system; integrator=rollout_integrator) push!(ψ̃_trajs, ψ̃_traj) end end @@ -465,27 +659,65 @@ function initialize_trajectory( ψ̃_inits, ψ̃_goals, state_names, - T, + N, Δt, args...; - a_guess=a_guess, + u_guess=u_guess, kwargs... ) end """ - initialize_trajectory + initialize_trajectory( + ρ_init, + ρ_goal, + N::Int, + Δt::Union{Real, AbstractVecOrMat{<:Real}}, + args...; + state_name::Symbol=:ρ⃗̃, + u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + system::Union{OpenQuantumSystem, Nothing}=nothing, + rollout_integrator::Function=expv, + kwargs... + ) + +Initialize a trajectory for open quantum system density matrix evolution. + +Constructs a trajectory for evolving a density matrix from an initial state to a target +state, supporting both unitary and open system dynamics. Density matrices are stored in +iso-vectorized form for optimization. -Trajectory initialization of density matrices. +# Arguments +- `ρ_init`: Initial density matrix +- `ρ_goal`: Target density matrix +- `N::Int`: Number of time steps +- `Δt::Union{Real, AbstractVecOrMat{<:Real}}`: Time step size(s) +- `args...`: Additional arguments passed to base `initialize_trajectory` + +# Keyword Arguments +- `state_name::Symbol=:ρ⃗̃`: Name for the density matrix state variable +- `u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing`: Initial control guess +- `system::Union{OpenQuantumSystem, Nothing}=nothing`: Open quantum system for rollout +- `rollout_integrator::Function=expv`: Integrator for open system dynamics +- Additional kwargs passed to the base `initialize_trajectory` method + +# Returns +- `NamedTrajectory`: Initialized trajectory with density matrix state, controls, and timesteps + +# Notes +- Density matrices are stored in iso-vectorized representation +- If `u_guess` is provided, requires `system` and uses rollout for state trajectory +- If `u_guess` is `nothing`, uses linear interpolation between initial and target states +- Uses `open_rollout` for open system dynamics with Lindblad master equation """ function initialize_trajectory( ρ_init, ρ_goal, - T::Int, + N::Int, Δt::Union{Real, AbstractVecOrMat{<:Real}}, args...; state_name::Symbol=:ρ⃗̃, - a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + u_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, system::Union{OpenQuantumSystem, Nothing}=nothing, rollout_integrator::Function=expv, kwargs... @@ -494,7 +726,7 @@ function initialize_trajectory( if Δt isa AbstractMatrix timesteps = vec(Δt) elseif Δt isa Float64 - timesteps = fill(Δt, T) + timesteps = fill(Δt, N) else timesteps = Δt end @@ -504,14 +736,14 @@ function initialize_trajectory( ρ⃗̃_goal = density_to_iso_vec(ρ_goal) # Construct state data - if isnothing(a_guess) - ρ⃗̃_traj = linear_interpolation(ρ_init, ρ_goal, T) + if isnothing(u_guess) + ρ⃗̃_traj = linear_interpolation(ρ_init, ρ_goal, N) else - @assert !isnothing(system) "System must be provided if a_guess is provided." + @assert !isnothing(system) "System must be provided if u_guess is provided." ρ⃗̃_traj = open_rollout( ρ_init, - a_guess, + u_guess, timesteps, system; integrator=rollout_integrator @@ -523,10 +755,10 @@ function initialize_trajectory( [ρ⃗̃_init], [ρ⃗̃_goal], [state_name], - T, + N, Δt, args...; - a_guess=a_guess, + u_guess=u_guess, kwargs... ) end