From 1a61fe296984d619fc358ec4e5b2a269211e6364 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Thu, 1 May 2025 11:02:54 -0500 Subject: [PATCH 01/12] drift compensation in init --- src/problem_templates/_problem_templates.jl | 6 +- src/trajectory_initialization.jl | 151 +++++++++----------- 2 files changed, 72 insertions(+), 85 deletions(-) diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 89d91eab..f6722691 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -57,10 +57,10 @@ function apply_piccolo_options!( end if free_time - if piccolo_options.verbose - println("\tapplying timesteps_all_equal constraint: $(traj.timestep)") - end if piccolo_options.timesteps_all_equal + if piccolo_options.verbose + println("\tapplying timesteps_all_equal constraint: $(traj.timestep)") + end push!( constraints, TimeStepsAllEqualConstraint(traj) diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 7ab0cdd3..1d5939fb 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -18,6 +18,8 @@ using TestItems # Initial states # # ----------------------------------------------------------------------------- # +linear_interpolation(x::AbstractVector, y::AbstractVector, n::Int) = hcat(range(x, y, n)...) + """ unitary_linear_interpolation( U_init::AbstractMatrix, @@ -47,81 +49,10 @@ function unitary_linear_interpolation( return unitary_linear_interpolation(U_init, U_goal.operator, samples) end -""" - unitary_geodesic( - operator::EmbeddedOperator, - samples::Int; - kwargs... - ) - - unitary_geodesic( - U_goal::AbstractMatrix{<:Number}, - samples::Int; - kwargs... - ) - - unitary_geodesic( - U₀::AbstractMatrix{<:Number}, - U₁::AbstractMatrix{<:Number}, - samples::Number; - kwargs... - ) - - unitary_geodesic( - U₀::AbstractMatrix{<:Number}, - U₁::AbstractMatrix{<:Number}, - timesteps::AbstractVector{<:Number}; - return_generator=false - ) - -Compute a geodesic connecting two unitary operators. -""" -function unitary_geodesic end - -function unitary_geodesic( - U_init::AbstractMatrix{<:Number}, - U_goal::EmbeddedOperator, - samples::Int; - kwargs... -) - U1 = unembed(U_init, U_goal) - U2 = unembed(U_goal) - Ũ⃗ = unitary_geodesic(U1, U2, samples; kwargs...) - return hcat([ - operator_to_iso_vec(embed(iso_vec_to_operator(Ũ⃗ₜ), U_goal)) - for Ũ⃗ₜ ∈ eachcol(Ũ⃗) - ]...) -end - -function unitary_geodesic( - U_init::AbstractMatrix{<:Number}, - U_goal::AbstractMatrix{<:Number}, - samples::Int; - kwargs... -) - return unitary_geodesic(U_init, U_goal, range(0, 1, samples); kwargs...) -end - -function unitary_geodesic( - U_goal::AbstractPiccoloOperator, - samples::Int; - kwargs... -) - if U_goal isa EmbeddedOperator - U_goal = U_goal.operator - end - return unitary_geodesic( - Matrix{ComplexF64}(I(size(U_goal, 1))), - U_goal, - samples; - kwargs... - ) -end - """ unitary_geodesic(U_init, U_goal, times; kwargs...) -Compute the geodesic connecting U_init and U_goal at the specified times. Allows for the possibility of unequal times and ranges outside [0,1]. +Compute the geodesic connecting U_init and U_goal at the specified times. # Arguments - `U_init::AbstractMatrix{<:Number}`: The initial unitary operator. @@ -129,21 +60,29 @@ Compute the geodesic connecting U_init and U_goal at the specified times. Allows - `times::AbstractVector{<:Number}`: The times at which to evaluate the geodesic. # Keyword Arguments -- `return_unitary_isos::Bool=true`: If true returns a matrix where each column is a unitary isovec, i.e. vec(vcat(real(U), imag(U))). If false, returns a vector of unitary matrices. -- `return_generator::Bool=false`: If true, returns the effective Hamiltonian generating the geodesic. +- `return_unitary_isos::Bool=true`: If true returns a matrix where each column is a unitary + isovec, i.e. vec(vcat(real(U), imag(U))). If false, returns a vector of unitary matrices. +- `return_generator::Bool=false`: If true, returns the effective Hamiltonian generating + the geodesic. """ +function unitary_geodesic end + function unitary_geodesic( U_init::AbstractMatrix{<:Number}, U_goal::AbstractMatrix{<:Number}, times::AbstractVector{<:Number}; return_unitary_isos=true, - return_generator=false + return_generator=false, + H_drift::AbstractMatrix{<:Number}=zeros(size(U_init)), ) t₀ = times[1] T = times[end] - t₀ - H = im * log(U_goal * U_init') / T + + U_drift(t) = exp(-im * H_drift * t) + H = im * log(U_drift(T)' * (U_goal * U_init')) / T # -im prefactor is not included in H - U_geo = [exp(-im * H * (t - t₀)) * U_init for t ∈ times] + U_geo = [U_drift(t) * exp(-im * H * (t - t₀)) * U_init for t ∈ times] + if !return_unitary_isos if return_generator return U_geo, H @@ -160,8 +99,44 @@ function unitary_geodesic( end end -linear_interpolation(x::AbstractVector, y::AbstractVector, n::Int) = - hcat(range(x, y, n)...) +function unitary_geodesic( + U_goal::AbstractPiccoloOperator, + samples::Int; + kwargs... +) + return unitary_geodesic( + I(size(U_goal)), + U_goal, + samples; + kwargs... + ) +end + +function unitary_geodesic( + U_init::AbstractMatrix{<:Number}, + U_goal::EmbeddedOperator, + samples::Int; + H_drift::AbstractMatrix{<:Number}=zeros(size(U_init)), + kwargs... +) + H_drift = unembed(H_drift, U_goal) + U1 = unembed(U_init, U_goal) + U2 = unembed(U_goal) + Ũ⃗ = unitary_geodesic(U1, U2, samples; H_drift=H_drift, kwargs...) + return hcat([ + operator_to_iso_vec(embed(iso_vec_to_operator(Ũ⃗ₜ), U_goal)) + for Ũ⃗ₜ ∈ eachcol(Ũ⃗) + ]...) +end + +function unitary_geodesic( + U_init::AbstractMatrix{<:Number}, + U_goal::AbstractMatrix{<:Number}, + samples::Int; + kwargs... +) + return unitary_geodesic(U_init, U_goal, range(0, 1, samples); kwargs...) +end # ============================================================================= # @@ -172,10 +147,16 @@ function initialize_unitary_trajectory( U_init::AbstractMatrix{<:Number}, U_goal::AbstractPiccoloOperator, T::Int; - geodesic=true + geodesic::Bool=true, + system::Union{AbstractQuantumSystem, Nothing}=nothing ) if geodesic - Ũ⃗ = unitary_geodesic(U_init, U_goal, T) + if system isa AbstractQuantumSystem + H_drift = Matrix(get_drift(system)) + else + H_drift = zeros(size(U_init)) + end + Ũ⃗ = unitary_geodesic(U_init, U_goal, T, H_drift=H_drift) else Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, T) end @@ -415,7 +396,13 @@ function initialize_trajectory( # Construct state data if isnothing(a_guess) - Ũ⃗_traj = initialize_unitary_trajectory(U_init, U_goal, T; geodesic=geodesic) + Ũ⃗_traj = initialize_unitary_trajectory( + U_init, + U_goal, + T; + 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) From 8daf3e1447289038f0ab35722a4d52a3bf00fc8b Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Fri, 2 May 2025 22:38:53 -0500 Subject: [PATCH 02/12] free phase template --- src/problem_templates/_problem_templates.jl | 1 + .../unitary_free_phase_problem.jl | 100 ++++++++++++++++++ src/quantum_objectives.jl | 18 +++- 3 files changed, 118 insertions(+), 1 deletion(-) create mode 100644 src/problem_templates/unitary_free_phase_problem.jl diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index f6722691..5407fc5f 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -20,6 +20,7 @@ include("unitary_smooth_pulse_problem.jl") include("unitary_variational_problem.jl") include("unitary_minimum_time_problem.jl") include("unitary_sampling_problem.jl") +include("unitary_free_phase_problem.jl") include("quantum_state_smooth_pulse_problem.jl") include("quantum_state_minimum_time_problem.jl") diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl new file mode 100644 index 00000000..41f80a9f --- /dev/null +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -0,0 +1,100 @@ +export UnitaryFreePhaseProblem + + +function UnitaryFreePhaseProblem( + system::AbstractQuantumSystem, + goal::Function, + T::Int, + Δt::Union{Float64, <:AbstractVector{Float64}}, + initial_phases::AbstractVector{Float64}; + unitary_integrator=UnitaryIntegrator, + state_name::Symbol = :Ũ⃗, + control_name::Symbol = :a, + timestep_name::Symbol = :Δt, + phase_name::Symbol = :θ, + 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=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), + Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), + Q::Float64=100.0, + R=1e-2, + R_a::Union{Float64, Vector{Float64}}=R, + R_da::Union{Float64, Vector{Float64}}=R, + R_dda::Union{Float64, Vector{Float64}}=R, + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], + piccolo_options::PiccoloOptions=PiccoloOptions(), +) + if piccolo_options.verbose + println(" constructing UnitarySmoothPulseProblem...") + println("\tusing integrator: $(typeof(unitary_integrator))") + println("\tinitial free phases: $(phase_name) = $(initial_phases)") + end + + # Trajectory + if !isnothing(init_trajectory) + traj = init_trajectory + else + traj = initialize_trajectory( + goal(initial_phases), + T, + Δt, + system.n_drives, + (a_bounds, da_bounds, dda_bounds); + state_name=state_name, + control_name=control_name, + timestep_name=timestep_name, + phase_name=phase_name, + phase_data=initial_phases, + Δ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, + a_guess=a_guess, + system=system, + rollout_integrator=piccolo_options.rollout_integrator, + verbose=piccolo_options.verbose + ) + end + + # Objective + J = UnitaryFreePhaseInfidelityObjective( + goal, + state_name, + phase_name, + traj + ) + + control_names = [ + name for name ∈ traj.names + 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) + + # Optional Piccolo constraints and objectives + ProblemTemplates.apply_piccolo_options!( + J, constraints, piccolo_options, traj, state_name, timestep_name; + state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing + ) + + integrators = [ + unitary_integrator(system, traj, state_name, control_name), + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), + ] + + return DirectTrajOptProblem( + traj, + J, + integrators; + constraints=constraints + ) +end \ No newline at end of file diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index 21ed1cc8..8f3718b5 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -4,6 +4,7 @@ export KetInfidelityObjective export UnitaryInfidelityObjective export DensityMatrixPureStateInfidelityObjective export UnitarySensitivityObjective +export UnitaryFreePhaseInfidelityObjective using LinearAlgebra using NamedTrajectories @@ -39,7 +40,7 @@ end function unitary_fidelity_loss( Ũ⃗::AbstractVector{<:Real}, - U_goal::AbstractMatrix{<:Complex{Float64}} + U_goal::AbstractMatrix{<:Complex{<:Real}} ) U = iso_vec_to_operator(Ũ⃗) n = size(U, 1) @@ -67,6 +68,21 @@ function UnitaryInfidelityObjective( return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q) end +function UnitaryFreePhaseInfidelityObjective( + U_goal::Function, + Ũ⃗_name::Symbol, + θ_name::Symbol, + traj::NamedTrajectory; + Q=100.0 +) + d = traj.global_dims[θ_name] + function ℓ(x) + Ũ⃗, θ = x[1:end-d], x[end-d+1:end] + return abs(1 - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ))) + end + return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q, global_names=[θ_name]) +end + # --------------------------------------------------------- # Density Matrices # --------------------------------------------------------- From 02802ae7b9494fbf034b3155add4ef5c3501f9dc Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Mon, 5 May 2025 12:29:53 -0500 Subject: [PATCH 03/12] update init for free phase problem, use global constraint and sin, cos --- .../unitary_free_phase_problem.jl | 36 +++++++++++++++---- src/quantum_objectives.jl | 20 ++++++++--- src/trajectory_initialization.jl | 12 +------ 3 files changed, 45 insertions(+), 23 deletions(-) diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 41f80a9f..33fb2d01 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -1,6 +1,8 @@ export UnitaryFreePhaseProblem +""" +""" function UnitaryFreePhaseProblem( system::AbstractQuantumSystem, goal::Function, @@ -36,12 +38,25 @@ function UnitaryFreePhaseProblem( println("\tinitial free phases: $(phase_name) = $(initial_phases)") end + # Construct global phases + x = Symbol("cos$(phase_name)") + y = Symbol("sin$(phase_name)") + phase_names = [x, y] + phase_data = (; x => cos.(initial_phases), y => sin.(initial_phases)) + trig_phases = [phase_data[x]; phase_data[y]] + if piccolo_options.verbose + println("\tusing global names: ", phase_names) + end + + @assert goal(trig_phases) isa AbstractPiccoloOperator "expected goal([cos(θ); sin(θ)])" + eval_goal = goal(trig_phases) + # Trajectory if !isnothing(init_trajectory) traj = init_trajectory else traj = initialize_trajectory( - goal(initial_phases), + eval_goal, T, Δt, system.n_drives, @@ -49,8 +64,6 @@ function UnitaryFreePhaseProblem( state_name=state_name, control_name=control_name, timestep_name=timestep_name, - phase_name=phase_name, - phase_data=initial_phases, Δt_bounds=(Δt_min, Δt_max), zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, geodesic=piccolo_options.geodesic, @@ -58,7 +71,8 @@ function UnitaryFreePhaseProblem( a_guess=a_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, - verbose=piccolo_options.verbose + verbose=piccolo_options.verbose, + global_data=phase_data, ) end @@ -66,8 +80,9 @@ function UnitaryFreePhaseProblem( J = UnitaryFreePhaseInfidelityObjective( goal, state_name, - phase_name, - traj + phase_names, + traj; + Q=Q ) control_names = [ @@ -82,8 +97,15 @@ function UnitaryFreePhaseProblem( # Optional Piccolo constraints and objectives ProblemTemplates.apply_piccolo_options!( J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing + state_leakage_indices=eval_goal isa EmbeddedOperator ? get_leakage_indices(eval_goal) : nothing ) + + # Phase constraint + function phase_norm(z) + x, y = z[1:length(z) ÷ 2], z[length(z) ÷ 2 + 1:end] + return x .^ 2 + y .^2 .- 1 + end + push!(constraints, GlobalNonlinearConstraint(phase_norm, phase_names, traj)) integrators = [ unitary_integrator(system, traj, state_name, control_name), diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index 8f3718b5..54bd340a 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -71,16 +71,26 @@ end function UnitaryFreePhaseInfidelityObjective( U_goal::Function, Ũ⃗_name::Symbol, - θ_name::Symbol, + θ_names::AbstractVector{Symbol}, traj::NamedTrajectory; Q=100.0 ) - d = traj.global_dims[θ_name] - function ℓ(x) - Ũ⃗, θ = x[1:end-d], x[end-d+1:end] + d = sum(traj.global_dims[n] for n in θ_names) + function ℓ(z) + Ũ⃗, θ = z[1:end-d], z[end-d+1:end] return abs(1 - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ))) end - return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q, global_names=[θ_name]) + return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q, global_names=θ_names) +end + +function UnitaryFreePhaseInfidelityObjective( + U_goal::Function, + Ũ⃗_name::Symbol, + θ_name::Symbol, + traj::NamedTrajectory; + kwargs... +) + return UnitaryFreePhaseInfidelityObjective(U_goal, Ũ⃗_name, [θ_name], traj; kwargs...) end # --------------------------------------------------------- diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 1d5939fb..eafd7a30 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -251,8 +251,7 @@ function initialize_trajectory( Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt), drive_derivative_σ::Float64=0.1, a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, - phase_name::Symbol=:ϕ, - phase_data::Union{AbstractVector{<:Real}, Nothing}=nothing, + global_data::NamedTuple{gname, <:Tuple{Vararg{AbstractVector{<:Real}}}} where gname=(;), verbose=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" @@ -342,9 +341,6 @@ function initialize_trajectory( controls = (control_names[end],) end - # Construct global data for free phases - global_data = isnothing(phase_data) ? (;) : (; phase_name => phase_data) - return NamedTrajectory( (; (names .=> values)...); controls=controls, @@ -373,7 +369,6 @@ function initialize_trajectory( system::Union{AbstractQuantumSystem, Nothing}=nothing, rollout_integrator::Function=expv, geodesic=true, - phase_operators::Union{AbstractVector{<:AbstractMatrix}, Nothing}=nothing, kwargs... ) # Construct timesteps @@ -407,10 +402,6 @@ function initialize_trajectory( @assert !isnothing(system) "System must be provided if a_guess is provided." Ũ⃗_traj = unitary_rollout(Ũ⃗_init, a_guess, timesteps, system; integrator=rollout_integrator) end - - # Construct phase data - phase_data = isnothing(phase_operators) ? nothing : π * randn(length(phase_operators)) - return initialize_trajectory( [Ũ⃗_traj], @@ -420,7 +411,6 @@ function initialize_trajectory( T, Δt, args...; - phase_data=phase_data, a_guess=a_guess, kwargs... ) From 156a89330c4faa11e1c6d35c30bda7b14a4d87a9 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Fri, 9 May 2025 19:28:22 -0500 Subject: [PATCH 04/12] update fidelity constraints for globals --- .../unitary_free_phase_problem.jl | 4 +- .../unitary_minimum_time_problem.jl | 62 +++++++++++++++++++ src/quantum_constraints.jl | 29 ++++++++- src/quantum_objectives.jl | 2 +- 4 files changed, 92 insertions(+), 5 deletions(-) diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 33fb2d01..2112f1b1 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -33,7 +33,7 @@ function UnitaryFreePhaseProblem( piccolo_options::PiccoloOptions=PiccoloOptions(), ) if piccolo_options.verbose - println(" constructing UnitarySmoothPulseProblem...") + println(" constructing UnitaryFreePhaseProblem...") println("\tusing integrator: $(typeof(unitary_integrator))") println("\tinitial free phases: $(phase_name) = $(initial_phases)") end @@ -105,7 +105,7 @@ function UnitaryFreePhaseProblem( x, y = z[1:length(z) ÷ 2], z[length(z) ÷ 2 + 1:end] return x .^ 2 + y .^2 .- 1 end - push!(constraints, GlobalNonlinearConstraint(phase_norm, phase_names, traj)) + push!(constraints, NonlinearGlobalConstraint(phase_norm, phase_names, traj)) integrators = [ unitary_integrator(system, traj, state_name, control_name), diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index 49a3f82d..9fb1586e 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -91,6 +91,68 @@ function UnitaryMinimumTimeProblem( ) end +# --------------------------------------------------------------------------- # +# Free phases +# --------------------------------------------------------------------------- # + +function UnitaryMinimumTimeProblem( + trajectory::NamedTrajectory, + goal::Function, + objective::Objective, + dynamics::TrajectoryDynamics, + constraints::AbstractVector{<:AbstractConstraint}; + unitary_name::Symbol=:Ũ⃗, + phase_name::Symbol=:θ, + final_fidelity::Float64=1.0, + D::Float64=100.0, + piccolo_options::PiccoloOptions=PiccoloOptions(), +) + # Collect phase names + phase_names = [ + n for n in keys(trajectory.global_data) if endswith(string(n), string(phase_name)) + ] + + if piccolo_options.verbose + println(" constructing UnitaryMinimumTimeProblem...") + println("\tfinal fidelity: $(final_fidelity)") + println("\tphase names: $(phase_names)") + end + + objective += MinimumTimeObjective( + trajectory; D=D, timesteps_all_equal=piccolo_options.timesteps_all_equal + ) + + fidelity_constraint = FinalUnitaryFidelityConstraint( + goal, unitary_name, phase_names, final_fidelity, trajectory + ) + + constraints = push!(constraints, fidelity_constraint) + + return DirectTrajOptProblem( + trajectory, + objective, + dynamics, + constraints + ) +end + +function UnitaryMinimumTimeProblem( + prob::DirectTrajOptProblem, + goal::Function; + objective::Objective=prob.objective, + constraints::AbstractVector{<:AbstractConstraint}=prob.constraints, + kwargs... +) + return UnitaryMinimumTimeProblem( + prob.trajectory, + goal, + objective, + prob.dynamics, + constraints; + kwargs... + ) +end + # *************************************************************************** # @testitem "Minimum time Hadamard gate" begin diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index 2e77dafd..48c4152d 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -21,7 +21,7 @@ function FinalKetFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = ψ̃ -> [ - abs(QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal) - final_fidelity) + final_fidelity - QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal) ] return NonlinearKnotPointConstraint( @@ -44,7 +44,7 @@ function FinalUnitaryFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = Ũ⃗ -> [ - abs(QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) - final_fidelity) + final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) ] return NonlinearKnotPointConstraint( @@ -56,4 +56,29 @@ function FinalUnitaryFidelityConstraint( ) end +function FinalUnitaryFidelityConstraint( + U_goal::Function, + Ũ⃗_name::Symbol, + θ_names::AbstractVector{Symbol}, + final_fidelity::Float64, + traj::NamedTrajectory +) + d = sum(traj.global_dims[n] for n in θ_names) + function terminal_constraint(z) + Ũ⃗, θ = z[1:end-d], z[end-d+1:end] + return [ + final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ)) + ] + end + + return NonlinearGlobalKnotPointConstraint( + terminal_constraint, + Ũ⃗_name, + θ_names, + traj, + equality=false, + times=[traj.T] + ) +end + end \ No newline at end of file diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index 54bd340a..e345c769 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -80,7 +80,7 @@ function UnitaryFreePhaseInfidelityObjective( Ũ⃗, θ = z[1:end-d], z[end-d+1:end] return abs(1 - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ))) end - return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q, global_names=θ_names) + return TerminalObjective(ℓ, Ũ⃗_name, θ_names, traj; Q=Q) end function UnitaryFreePhaseInfidelityObjective( From b01c9f34c810a20219eec854106218a15aaf2b2f Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Fri, 9 May 2025 20:08:27 -0500 Subject: [PATCH 05/12] deep copy constraints in min time --- src/problem_templates/unitary_minimum_time_problem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index 9fb1586e..dea98b72 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -78,7 +78,7 @@ function UnitaryMinimumTimeProblem( prob::DirectTrajOptProblem, goal::AbstractPiccoloOperator; objective::Objective=prob.objective, - constraints::AbstractVector{<:AbstractConstraint}=prob.constraints, + constraints::AbstractVector{<:AbstractConstraint}=deepcopy(prob.constraints), kwargs... ) return UnitaryMinimumTimeProblem( @@ -140,7 +140,7 @@ function UnitaryMinimumTimeProblem( prob::DirectTrajOptProblem, goal::Function; objective::Objective=prob.objective, - constraints::AbstractVector{<:AbstractConstraint}=prob.constraints, + constraints::AbstractVector{<:AbstractConstraint}=deepcopy(prob.constraints), kwargs... ) return UnitaryMinimumTimeProblem( From 87bbdc4ad1bfd83a011ee13539811fcb30c7cea9 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Sun, 18 May 2025 13:44:25 -0500 Subject: [PATCH 06/12] log fid constraint --- src/quantum_constraints.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index 48c4152d..f84b4135 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -21,7 +21,7 @@ function FinalKetFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = ψ̃ -> [ - final_fidelity - QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal) + log(final_fidelity / QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal)) ] return NonlinearKnotPointConstraint( @@ -44,7 +44,8 @@ function FinalUnitaryFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = Ũ⃗ -> [ - final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) + # final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) + log(final_fidelity / QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal)) ] return NonlinearKnotPointConstraint( @@ -67,7 +68,8 @@ function FinalUnitaryFidelityConstraint( function terminal_constraint(z) Ũ⃗, θ = z[1:end-d], z[end-d+1:end] return [ - final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ)) + # final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ)) + log(final_fidelity / QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ))) ] end From 5d843e9de98196b3ee3e8fbddef0077aa95f05f0 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Sun, 18 May 2025 13:44:44 -0500 Subject: [PATCH 07/12] goal from phases --- .../unitary_free_phase_problem.jl | 23 +++++++++++-------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 2112f1b1..3a55bd01 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -1,5 +1,8 @@ export UnitaryFreePhaseProblem +mean(x::AbstractVector) = sum(x) / length(x) + + """ """ @@ -7,13 +10,13 @@ function UnitaryFreePhaseProblem( system::AbstractQuantumSystem, goal::Function, T::Int, - Δt::Union{Float64, <:AbstractVector{Float64}}, - initial_phases::AbstractVector{Float64}; + Δt::Union{Float64, <:AbstractVector{Float64}}; unitary_integrator=UnitaryIntegrator, state_name::Symbol = :Ũ⃗, control_name::Symbol = :a, 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, @@ -35,26 +38,28 @@ function UnitaryFreePhaseProblem( if piccolo_options.verbose println(" constructing UnitaryFreePhaseProblem...") println("\tusing integrator: $(typeof(unitary_integrator))") - println("\tinitial free phases: $(phase_name) = $(initial_phases)") + println("\tinitial free phases: $(phase_name) = $(init_phases)") end # Construct global phases x = Symbol("cos$(phase_name)") y = Symbol("sin$(phase_name)") phase_names = [x, y] - phase_data = (; x => cos.(initial_phases), y => sin.(initial_phases)) - trig_phases = [phase_data[x]; phase_data[y]] if piccolo_options.verbose println("\tusing global names: ", phase_names) end - @assert goal(trig_phases) isa AbstractPiccoloOperator "expected goal([cos(θ); sin(θ)])" - eval_goal = goal(trig_phases) - # Trajectory if !isnothing(init_trajectory) + trig_phases = [init_trajectory.global_data[x]; init_trajectory.global_data[y]] + @assert goal(trig_phases) isa AbstractPiccoloOperator "expected goal([cos(θ); sin(θ)])" + eval_goal = goal(trig_phases) traj = init_trajectory else + phase_data = (; x => cos.(init_phases), y => sin.(init_phases)) + trig_phases = [phase_data[x]; phase_data[y]] + @assert goal(trig_phases) isa AbstractPiccoloOperator "expected goal([cos(θ); sin(θ)])" + eval_goal = goal(trig_phases) traj = initialize_trajectory( eval_goal, T, @@ -67,7 +72,7 @@ function UnitaryFreePhaseProblem( Δ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, + bound_state=false, # TODO: Hardcoded a_guess=a_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, From 6ae3726e1e7e257965d56f332d92b8a84eb7ce58 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Mon, 19 May 2025 11:32:28 -0500 Subject: [PATCH 08/12] only free time version of min time obj --- src/problem_templates/unitary_minimum_time_problem.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index dea98b72..ae40eba3 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -53,9 +53,8 @@ function UnitaryMinimumTimeProblem( println("\tfinal fidelity: $(final_fidelity)") end - objective += MinimumTimeObjective( - trajectory; D=D, timesteps_all_equal=piccolo_options.timesteps_all_equal - ) + objective += MinimumTimeObjective(trajectory, D=D) + # timesteps_all_equal=piccolo_options.timesteps_all_equal fidelity_constraint = FinalUnitaryFidelityConstraint( goal, unitary_name, final_fidelity, trajectory @@ -118,9 +117,8 @@ function UnitaryMinimumTimeProblem( println("\tphase names: $(phase_names)") end - objective += MinimumTimeObjective( - trajectory; D=D, timesteps_all_equal=piccolo_options.timesteps_all_equal - ) + objective += MinimumTimeObjective(trajectory, D=D) + # timesteps_all_equal=piccolo_options.timesteps_all_equal fidelity_constraint = FinalUnitaryFidelityConstraint( goal, unitary_name, phase_names, final_fidelity, trajectory From e983815ed50aea1599a3504914b0b78f531bf5e1 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Fri, 13 Jun 2025 13:39:11 -0500 Subject: [PATCH 09/12] update for new NT --- Project.toml | 2 +- src/problem_templates/unitary_variational_problem.jl | 3 +++ src/trajectory_initialization.jl | 4 +--- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 7458db6a..c036b5b7 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ Interpolations = "0.15" JLD2 = "0.5" LinearAlgebra = "1.10, 1.11" NamedTrajectories = "0.3" -PiccoloQuantumObjects = "0.4" +PiccoloQuantumObjects = "0.5" Random = "1.10, 1.11" Reexport = "1.2" SparseArrays = "1.10, 1.11" diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl index 8bc013b7..1d702549 100644 --- a/src/problem_templates/unitary_variational_problem.jl +++ b/src/problem_templates/unitary_variational_problem.jl @@ -87,6 +87,8 @@ function UnitaryVariationalProblem( constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) + @assert true "Broken by NT refactor" + if piccolo_options.verbose println(" constructing UnitaryVariationalProblem...") println("\tusing integrator: $(typeof(variational_integrator))") @@ -126,6 +128,7 @@ function UnitaryVariationalProblem( Symbol(string(variational_state_name) * "$(i)") for i in eachindex(system.G_vars) ] + # TODO: Must refactor for (name, scale, Ũ⃗_v) in zip(variational_state_names, variational_scales, Ũ⃗_vars) add_component!(traj, name, Ũ⃗_v / scale; type=:state) traj.initial = merge(traj.initial, (name => Ũ⃗_v[:, 1] / scale,)) diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index eafd7a30..d708a16c 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -6,6 +6,7 @@ export unitary_linear_interpolation export initialize_trajectory using NamedTrajectories +import NamedTrajectories: ScalarBound, VectorBound using PiccoloQuantumObjects using Distributions @@ -140,9 +141,6 @@ end # ============================================================================= # -const VectorBound = Union{AbstractVector{R}, Tuple{AbstractVector{R}, AbstractVector{R}}} where R <: Real -const ScalarBound = Union{R, Tuple{R, R}} where R <: Real - function initialize_unitary_trajectory( U_init::AbstractMatrix{<:Number}, U_goal::AbstractPiccoloOperator, From 0af3ab310c01b2426615f83c059bde42edc0a772 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Tue, 17 Jun 2025 10:08:06 -0500 Subject: [PATCH 10/12] parametric NT, free phase and variational PT tests --- Project.toml | 5 +- src/problem_templates/_problem_templates.jl | 10 ++-- .../unitary_free_phase_problem.jl | 46 ++++++++++++++++-- .../unitary_minimum_time_problem.jl | 9 ++-- .../unitary_smooth_pulse_problem.jl | 2 +- .../unitary_variational_problem.jl | 43 ++++++++++++----- src/quantum_constraints.jl | 3 +- src/trajectory_initialization.jl | 47 +++++++------------ 8 files changed, 103 insertions(+), 62 deletions(-) diff --git a/Project.toml b/Project.toml index c036b5b7..71b233b6 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.2" +version = "0.8.0" [deps] DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" @@ -21,14 +21,11 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8" [compat] -DirectTrajOpt = "0.2.3" Distributions = "0.25" ExponentialAction = "0.2" Interpolations = "0.15" JLD2 = "0.5" LinearAlgebra = "1.10, 1.11" -NamedTrajectories = "0.3" -PiccoloQuantumObjects = "0.5" Random = "1.10, 1.11" Reexport = "1.2" SparseArrays = "1.10, 1.11" diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 5407fc5f..a04852f9 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -10,12 +10,16 @@ using TrajectoryIndexingUtils using NamedTrajectories using DirectTrajOpt using PiccoloQuantumObjects -using LinearAlgebra -using SparseArrays + using ExponentialAction using JLD2 +using LinearAlgebra +using SparseArrays using TestItems +# using Statistics +mean(x) = sum(x) / length(x) + include("unitary_smooth_pulse_problem.jl") include("unitary_variational_problem.jl") include("unitary_minimum_time_problem.jl") @@ -26,7 +30,7 @@ include("quantum_state_smooth_pulse_problem.jl") include("quantum_state_minimum_time_problem.jl") include("quantum_state_sampling_problem.jl") - +# TODO: refactor function apply_piccolo_options!( J::Objective, constraints::AbstractVector{<:AbstractConstraint}, diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 3a55bd01..41d75212 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -1,10 +1,13 @@ export UnitaryFreePhaseProblem -mean(x::AbstractVector) = sum(x) / length(x) - """ + UnitaryFreePhaseProblem(system::AbstractQuantumSystem, goal::Function, T, Δt; 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`. """ function UnitaryFreePhaseProblem( system::AbstractQuantumSystem, @@ -72,12 +75,12 @@ function UnitaryFreePhaseProblem( Δt_bounds=(Δt_min, Δt_max), zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, geodesic=piccolo_options.geodesic, - bound_state=false, # TODO: Hardcoded + bound_state=piccolo_options.bound_state, a_guess=a_guess, system=system, rollout_integrator=piccolo_options.rollout_integrator, verbose=piccolo_options.verbose, - global_data=phase_data, + global_component_data=phase_data, ) end @@ -124,4 +127,39 @@ function UnitaryFreePhaseProblem( integrators; constraints=constraints ) +end + +# *************************************************************************** # + +@testitem "UnitaryFreePhaseProblem: basic construction" begin + using PiccoloQuantumObjects + + sys = QuantumSystem(0.3 * PAULIS.X, [PAULIS.Y]) + U_goal = GATES.Z + T = 51 + Δt = 0.2 + + function virtual_z(z::AbstractVector{<:Real}) + x, y = z[1:length(z)÷2], z[1+length(z)÷2:end] + # U_goal ≈ R * U + R = reduce(kron, [xᵢ * PAULIS.I + im * yᵢ * PAULIS.Z for (xᵢ, yᵢ) in zip(x, y)]) + return R'U_goal + end + + initial_phases = [pi/3] + + prob = UnitaryFreePhaseProblem( + sys, virtual_z, T, Δt, + init_phases=initial_phases, + piccolo_options=PiccoloOptions(verbose=false), + phase_name=:ϕ, + ) + + @test prob isa DirectTrajOptProblem + @test length(prob.trajectory.global_data) == 2length(initial_phases) + @test prob.trajectory.global_names == (:cosϕ, :sinϕ) + + before = copy(prob.trajectory.global_data) + solve!(prob, max_iter=10, verbose=false, print_level=1) + @test prob.trajectory.global_data != before end \ No newline at end of file diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index dcbbaa7a..b8ca0f8b 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -54,7 +54,6 @@ function UnitaryMinimumTimeProblem( end objective += MinimumTimeObjective(trajectory, D=D) - # timesteps_all_equal=piccolo_options.timesteps_all_equal fidelity_constraint = FinalUnitaryFidelityConstraint( goal, unitary_name, final_fidelity, trajectory @@ -81,7 +80,7 @@ function UnitaryMinimumTimeProblem( kwargs... ) return UnitaryMinimumTimeProblem( - prob.trajectory, + deepcopy(prob.trajectory), goal, objective, prob.dynamics, @@ -107,9 +106,7 @@ function UnitaryMinimumTimeProblem( piccolo_options::PiccoloOptions=PiccoloOptions(), ) # Collect phase names - phase_names = [ - n for n in keys(trajectory.global_data) if endswith(string(n), string(phase_name)) - ] + phase_names = [n for n in trajectory.global_names if endswith(string(n), string(phase_name))] if piccolo_options.verbose println(" constructing UnitaryMinimumTimeProblem...") @@ -142,7 +139,7 @@ function UnitaryMinimumTimeProblem( kwargs... ) return UnitaryMinimumTimeProblem( - prob.trajectory, + deepcopy(prob.trajectory), goal, objective, prob.dynamics, diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 2a560f44..f7e85f98 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -2,7 +2,7 @@ export UnitarySmoothPulseProblem @doc raw""" - UnitarySmoothPulseProblem(system::AbstractQuantumSystem, operator, T, Δt; kwargs...) + UnitarySmoothPulseProblem(system::AbstractQuantumSystem, operator::AbstractPiccoloOperator, T::Int, Δt::Float64; kwargs...) UnitarySmoothPulseProblem(H_drift, H_drives, operator, T, Δt; 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., diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl index 1d702549..43dd666c 100644 --- a/src/problem_templates/unitary_variational_problem.jl +++ b/src/problem_templates/unitary_variational_problem.jl @@ -15,10 +15,14 @@ export UnitaryVariationalProblem 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. + +# 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`). @@ -44,10 +48,12 @@ Constructs a unitary variational problem for optimizing quantum control trajecto - `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. @@ -87,12 +93,16 @@ function UnitaryVariationalProblem( constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) - @assert true "Broken by NT refactor" - if piccolo_options.verbose println(" constructing UnitaryVariationalProblem...") println("\tusing integrator: $(typeof(variational_integrator))") println("\ttotal variational parameters: $(length(system.G_vars))") + if !isempty(robust_times) + println("\trobust knot points: $(robust_times)") + end + if !isempty(sensitive_times) + println("\tsensitive knot points: $(sensitive_times)") + end end # Trajectory @@ -124,15 +134,22 @@ function UnitaryVariationalProblem( drive_name=control_name ) - variational_state_names = [ + # Add variational components to the trajectory + var_state_names = Tuple( Symbol(string(variational_state_name) * "$(i)") for i in eachindex(system.G_vars) - ] - - # TODO: Must refactor - for (name, scale, Ũ⃗_v) in zip(variational_state_names, variational_scales, Ũ⃗_vars) - add_component!(traj, name, Ũ⃗_v / scale; type=:state) - traj.initial = merge(traj.initial, (name => Ũ⃗_v[:, 1] / scale,)) - end + ) + var_comps_inits = NamedTuple{var_state_names}( + Ũ⃗_v[:, 1] / scale for (scale, Ũ⃗_v) in zip(variational_scales, Ũ⃗_vars) + ) + var_comps_data = NamedTuple{var_state_names}( + Ũ⃗_v / scale for (scale, Ũ⃗_v) in zip(variational_scales, Ũ⃗_vars) + ) + traj = add_components( + traj, + var_comps_data; + type=:state, + initial=merge(traj.initial, var_comps_inits) + ) control_names = [ name for name ∈ traj.names @@ -147,7 +164,7 @@ function UnitaryVariationalProblem( # sensitivity for (name, scale, s, r) ∈ zip( - variational_state_names, + var_state_names, variational_scales, sensitive_times, robust_times @@ -170,7 +187,7 @@ function UnitaryVariationalProblem( integrators = [ variational_integrator( - system, traj, state_name, variational_state_names, control_name, + system, traj, state_name, [var_state_names...], control_name, scales=variational_scales ), DerivativeIntegrator(traj, control_name, control_names[2]), @@ -216,5 +233,5 @@ end sense_n = norm(sense_scale * sense_prob.trajectory.Ũ⃗ᵥ1) rob_n = norm(rob_scale * rob_prob.trajectory.Ũ⃗ᵥ1) - @assert sense_n > rob_n + @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 d9d30427..48c4152d 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -67,8 +67,7 @@ function FinalUnitaryFidelityConstraint( function terminal_constraint(z) Ũ⃗, θ = z[1:end-d], z[end-d+1:end] return [ - # final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ)) - log(final_fidelity / QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ))) + final_fidelity - QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal(θ)) ] end diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index d708a16c..19e0424b 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -6,7 +6,7 @@ export unitary_linear_interpolation export initialize_trajectory using NamedTrajectories -import NamedTrajectories: ScalarBound, VectorBound +import NamedTrajectories.StructNamedTrajectory: ScalarBound, VectorBound using PiccoloQuantumObjects using Distributions @@ -241,7 +241,6 @@ function initialize_trajectory( n_drives::Int, control_bounds::Tuple{Vararg{VectorBound}}; bound_state=false, - free_time=true, control_name=:a, n_control_derivatives::Int=length(control_bounds) - 1, zero_initial_and_final_derivative=false, @@ -249,7 +248,7 @@ function initialize_trajectory( Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt), drive_derivative_σ::Float64=0.1, a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, - global_data::NamedTuple{gname, <:Tuple{Vararg{AbstractVector{<:Real}}}} where gname=(;), + global_component_data::NamedTuple{gname, <:Tuple{Vararg{AbstractVector{<:Real}}}} where gname=(;), verbose=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" @@ -271,19 +270,15 @@ function initialize_trajectory( control_bounds = NamedTuple{control_names}(control_bounds) # Timestep data - if free_time - if Δt isa Real - Δt = fill(Δt, 1, T) - elseif Δt isa AbstractVector - Δt = reshape(Δt, 1, :) - else - @assert size(Δt) == (1, T) "Δt must be a Real, AbstractVector, or 1x$(T) AbstractMatrix" - end - timestep = timestep_name + if Δt isa Real + timestep_data = fill(Δt, 1, T) + elseif Δt isa AbstractVector + timestep_data = reshape(Δt, 1, :) else - @assert Δt isa Real "Δt must be a Real if free_time is false" - timestep = Δt + timestep_data = Δt + @assert size(Δt) == (1, T) "Δt must be a Real, AbstractVector, or 1x$(T) AbstractMatrix" end + timestep = timestep_name # Constraints initial = (; @@ -305,6 +300,8 @@ function initialize_trajectory( # Bounds bounds = control_bounds + bounds = merge(bounds, (; timestep_name => Δt_bounds,)) + # Put unit box bounds on the state if bound_state is true if bound_state state_dim = length(state_inits[1]) @@ -315,7 +312,7 @@ function initialize_trajectory( # Trajectory if isnothing(a_guess) # Randomly sample controls - a_values = initialize_control_trajectory( + control_data = initialize_control_trajectory( n_drives, n_control_derivatives, T, @@ -324,30 +321,22 @@ function initialize_trajectory( ) else # Use provided controls and take derivatives - a_values = initialize_control_trajectory(a_guess, Δt, n_control_derivatives) + control_data = initialize_control_trajectory(a_guess, Δt, n_control_derivatives) end - names = [state_names..., control_names...] - values = [state_data..., a_values...] - - if free_time - push!(names, timestep_name) - push!(values, Δt) - controls = (control_names[end], timestep_name) - bounds = merge(bounds, (; timestep_name => Δt_bounds,)) - else - controls = (control_names[end],) - end + names = [state_names..., control_names..., timestep_name] + values = [state_data..., control_data..., timestep_data] + controls = (control_names[end], timestep_name) return NamedTrajectory( - (; (names .=> values)...); + (; (names .=> values)...), + global_component_data; controls=controls, timestep=timestep, bounds=bounds, initial=initial, final=final, goal=goal, - global_data=global_data ) end From e793c6a4f3579aaa072e827b44019de64f541f03 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Mon, 23 Jun 2025 12:11:05 -0500 Subject: [PATCH 11/12] add leakage constraint to piccolo options --- src/piccolo_options.jl | 8 +- src/problem_templates/_problem_templates.jl | 99 +++++++------------ .../quantum_state_sampling_problem.jl | 14 +-- .../quantum_state_smooth_pulse_problem.jl | 12 ++- .../unitary_free_phase_problem.jl | 15 +-- .../unitary_sampling_problem.jl | 13 ++- .../unitary_smooth_pulse_problem.jl | 63 +++++++----- .../unitary_variational_problem.jl | 13 ++- src/quantum_constraints.jl | 4 +- src/quantum_objectives.jl | 2 +- 10 files changed, 120 insertions(+), 123 deletions(-) diff --git a/src/piccolo_options.jl b/src/piccolo_options.jl index db9c3a80..d51aa647 100644 --- a/src/piccolo_options.jl +++ b/src/piccolo_options.jl @@ -19,7 +19,7 @@ Options for the Piccolo quantum optimal control library. - `complex_control_norm_constraint_radius::Float64 = 1.0`: Radius of the complex control norm constraint. - `bound_state::Bool = false`: Bound the state. - `leakage_suppression::Bool = false`: Suppress leakage. -- `R_leakage::Float64 = 1.0`: Leakage suppression parameter. +- `leakage_cost::Float64 = 1.0`: Leakage suppression parameter. """ @kwdef mutable struct PiccoloOptions verbose::Bool = true @@ -30,9 +30,9 @@ Options for the Piccolo quantum optimal control library. complex_control_norm_constraint_name::Union{Nothing, Symbol} = nothing complex_control_norm_constraint_radius::Float64 = 1.0 bound_state::Bool = false - leakage_suppression::Bool = false - state_leakage_indices::AbstractVector{Int} = Int[] - R_leakage::Float64 = 1.0 + leakage_constraint::Bool = false + leakage_constraint_value::Float64 = 1e-2 + leakage_cost::Float64 = 1e-2 end end \ No newline at end of file diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index a04852f9..2a263ab7 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -17,9 +17,6 @@ using LinearAlgebra using SparseArrays using TestItems -# using Statistics -mean(x) = sum(x) / length(x) - include("unitary_smooth_pulse_problem.jl") include("unitary_variational_problem.jl") include("unitary_minimum_time_problem.jl") @@ -30,47 +27,45 @@ include("quantum_state_smooth_pulse_problem.jl") include("quantum_state_minimum_time_problem.jl") include("quantum_state_sampling_problem.jl") -# TODO: refactor + function apply_piccolo_options!( - J::Objective, - constraints::AbstractVector{<:AbstractConstraint}, piccolo_options::PiccoloOptions, - traj::NamedTrajectory, - state_names::AbstractVector{Symbol}, - timestep_name::Symbol; - state_leakage_indices::Union{Nothing, AbstractVector{<:AbstractVector{Int}}}=nothing, - free_time::Bool=true, + constraints::AbstractVector{<:AbstractConstraint}, + traj::NamedTrajectory; + state_names::Union{Nothing, Symbol, AbstractVector{Symbol}}=nothing, + state_leakage_indices::Union{Nothing, AbstractVector{Int}, AbstractVector{<:AbstractVector{Int}}}=nothing, ) - if piccolo_options.leakage_suppression - throw(error("L1 is not implemented.")) - # if piccolo_options.verbose - # println("\tapplying leakage suppression: $(state_names)") - # end + J = NullObjective(traj) + + if piccolo_options.leakage_constraint + val = piccolo_options.leakage_constraint_value + if piccolo_options.verbose + println("\tapplying leakage suppression: $(state_names) < $(val)") + end + + if isnothing(state_leakage_indices) + throw(ValueError("Leakage indices are required for leakage suppression.")) + end + + if state_names isa Symbol + state_names = [state_names] + state_leakage_indices = [state_leakage_indices] + end - # if isnothing(state_leakage_indices) - # error("You must provide leakage indices for leakage suppression.") - # end - # for (state_name, leakage_indices) ∈ zip(state_names, state_leakage_indices) - # J += L1Regularizer!( - # constraints, - # state_name, - # traj; - # R_value=piccolo_options.R_leakage, - # indices=leakage_indices, - # ) - # end + for (name, indices) ∈ zip(state_names, state_leakage_indices) + J += LeakageObjective(indices, name, traj, Qs=fill(piccolo_options.leakage_cost, traj.T)) + push!(constraints, LeakageConstraint(val, indices, name, traj)) + end end - if free_time - if piccolo_options.timesteps_all_equal - if piccolo_options.verbose - println("\tapplying timesteps_all_equal constraint: $(traj.timestep)") - end - push!( - constraints, - TimeStepsAllEqualConstraint(traj) - ) + if piccolo_options.timesteps_all_equal + if piccolo_options.verbose + println("\tapplying timesteps_all_equal constraint: $(traj.timestep)") end + push!( + constraints, + TimeStepsAllEqualConstraint(traj) + ) end if !isnothing(piccolo_options.complex_control_norm_constraint_name) @@ -86,35 +81,7 @@ function apply_piccolo_options!( push!(constraints, norm_con) end - return + return J end -function apply_piccolo_options!( - J::Objective, - constraints::AbstractVector{<:AbstractConstraint}, - piccolo_options::PiccoloOptions, - traj::NamedTrajectory, - state_name::Symbol, - timestep_name::Symbol; - state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, - kwargs... -) - state_names = [ - name for name ∈ traj.names - if startswith(string(name), string(state_name)) - ] - - return apply_piccolo_options!( - J, - constraints, - piccolo_options, - traj, - state_names, - timestep_name; - state_leakage_indices=isnothing(state_leakage_indices) ? nothing : fill(state_leakage_indices, length(state_names)), - kwargs... - ) -end - - end diff --git a/src/problem_templates/quantum_state_sampling_problem.jl b/src/problem_templates/quantum_state_sampling_problem.jl index a80fcb3d..2d1b4a7c 100644 --- a/src/problem_templates/quantum_state_sampling_problem.jl +++ b/src/problem_templates/quantum_state_sampling_problem.jl @@ -10,7 +10,7 @@ function QuantumStateSamplingProblem( ψ_inits::AbstractVector{<:AbstractVector{<:AbstractVector{<:ComplexF64}}}, ψ_goals::AbstractVector{<:AbstractVector{<:AbstractVector{<:ComplexF64}}}, T::Int, - Δt::Union{Float64,Vector{Float64}}; + Δt::Union{Float64, AbstractVector{Float64}}; ket_integrator=KetIntegrator, system_weights=fill(1.0, length(systems)), init_trajectory::Union{NamedTrajectory,Nothing}=nothing, @@ -24,13 +24,14 @@ function QuantumStateSamplingProblem( da_bounds::Vector{Float64}=fill(da_bound, systems[1].n_drives), dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, systems[1].n_drives), - Δt_min::Float64=0.5 * Δt, - Δt_max::Float64=1.5 * Δt, + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=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, + state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) @@ -93,9 +94,10 @@ function QuantumStateSamplingProblem( end # Optional Piccolo constraints and objectives - apply_piccolo_options!( - J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=piccolo_options.state_leakage_indices + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=vcat(state_names...), + state_leakage_indices=isnothing(state_leakage_indices) ? nothing : fill(state_leakage_indices, length(vcat(state_names...))), ) # Integrators diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 821b680b..eacba975 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -69,13 +69,14 @@ function QuantumStateSmoothPulseProblem( da_bounds::Vector{Float64}=fill(da_bound, sys.n_drives), dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, sys.n_drives), - Δt_min::Float64=0.001 * Δt, - Δt_max::Float64=2.0 * Δt, + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=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, + state_leakage_indices::Union{Nothing, AbstractVector{Int}}=nothing, constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], piccolo_options::PiccoloOptions=PiccoloOptions(), ) @@ -131,9 +132,10 @@ function QuantumStateSmoothPulseProblem( end # Optional Piccolo constraints and objectives - apply_piccolo_options!( - J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=piccolo_options.state_leakage_indices + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=state_names, + state_leakage_indices= isnothing(state_leakage_indices) ? nothing : fill(state_leakage_indices, length(state_names)) ) state_names = [ diff --git a/src/problem_templates/unitary_free_phase_problem.jl b/src/problem_templates/unitary_free_phase_problem.jl index 41d75212..d7fa22d1 100644 --- a/src/problem_templates/unitary_free_phase_problem.jl +++ b/src/problem_templates/unitary_free_phase_problem.jl @@ -13,7 +13,7 @@ function UnitaryFreePhaseProblem( system::AbstractQuantumSystem, goal::Function, T::Int, - Δt::Union{Float64, <:AbstractVector{Float64}}; + Δt::Union{Float64, AbstractVector{Float64}}; unitary_integrator=UnitaryIntegrator, state_name::Symbol = :Ũ⃗, control_name::Symbol = :a, @@ -28,8 +28,8 @@ function UnitaryFreePhaseProblem( 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), + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=2.0 * maximum(Δt), Q::Float64=100.0, R=1e-2, R_a::Union{Float64, Vector{Float64}}=R, @@ -103,9 +103,12 @@ function UnitaryFreePhaseProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) # Optional Piccolo constraints and objectives - ProblemTemplates.apply_piccolo_options!( - J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=eval_goal isa EmbeddedOperator ? get_leakage_indices(eval_goal) : nothing + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=state_name, + state_leakage_indices=eval_goal isa EmbeddedOperator ? + get_iso_vec_leakage_indices(eval_goal) : + nothing ) # Phase constraint diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index f8a7abbc..6387cf80 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -61,8 +61,8 @@ function UnitarySamplingProblem( da_bounds::Vector{Float64}=fill(da_bound, systems[1].n_drives), dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, systems[1].n_drives), - Δt_min::Float64=0.5 * Δt, - Δt_max::Float64=1.5 * Δt, + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=2.0 * maximum(Δt), Q::Float64=100.0, R=1e-2, R_a::Union{Float64,Vector{Float64}}=R, @@ -126,9 +126,12 @@ function UnitarySamplingProblem( end # Optional Piccolo constraints and objectives - apply_piccolo_options!( - J, constraints, piccolo_options, traj, state_names, timestep_name; - state_leakage_indices=all(op -> op isa EmbeddedOperator, operators) ? get_leakage_indices.(operators) : nothing + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=state_names, + state_leakage_indices=all(op -> op isa EmbeddedOperator, operators) ? + get_iso_vec_leakage_indices.(operators) : + nothing ) # Integrators diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index f7e85f98..4bb25e8b 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -82,8 +82,8 @@ function UnitarySmoothPulseProblem( 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), + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=2.0 * maximum(Δt), Q::Float64=100.0, R=1e-2, R_a::Union{Float64, Vector{Float64}}=R, @@ -134,9 +134,12 @@ function UnitarySmoothPulseProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) # Optional Piccolo constraints and objectives - apply_piccolo_options!( - J, constraints, piccolo_options, traj, state_name, timestep_name; - state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=state_name, + state_leakage_indices=goal isa EmbeddedOperator ? + get_iso_vec_leakage_indices(goal) : + nothing ) integrators = [ @@ -165,7 +168,7 @@ end # *************************************************************************** # -@testitem "Hadamard gate" begin +@testitem "Hadamard gate improvement" begin using PiccoloQuantumObjects sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) @@ -181,11 +184,10 @@ end initial = unitary_rollout_fidelity(prob.trajectory, sys) solve!(prob, max_iter=100, verbose=false, print_level=1) - final = unitary_rollout_fidelity(prob.trajectory, sys) - @test final > initial + @test unitary_rollout_fidelity(prob.trajectory, sys) > initial end -@testitem "Hadamard gate with bound states and control norm constraint" begin +@testitem "Bound states and control norm constraint" begin using PiccoloQuantumObjects sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) @@ -202,14 +204,12 @@ end ) ) - initial = unitary_rollout_fidelity(prob.trajectory, sys) - solve!(prob, max_iter=100, verbose=false, print_level=1) - final = unitary_rollout_fidelity(prob.trajectory, sys) - @test_broken false - # @test final > initial + initial = copy(prob.trajectory.data) + solve!(prob, max_iter=5, verbose=false, print_level=1) + @test prob.trajectory.data != initial end -@testitem "EmbeddedOperator Hadamard gate" begin +@testitem "EmbeddedOperator tests" begin using PiccoloQuantumObjects a = annihilate(3) @@ -218,15 +218,32 @@ end T = 51 Δt = 0.2 - prob = UnitarySmoothPulseProblem( - sys, U_goal, T, Δt, - piccolo_options=PiccoloOptions(verbose=false) - ) + @testset "EmbeddedOperator: solve gate" begin + prob = UnitarySmoothPulseProblem( + sys, U_goal, T, Δt, + piccolo_options=PiccoloOptions(verbose=false) + ) - initial = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) - solve!(prob, max_iter=100, verbose=false, print_level=1) - final = unitary_rollout_fidelity(prob.trajectory, sys, subspace=U_goal.subspace) - @test final > initial + initial = copy(prob.trajectory.data) + solve!(prob, max_iter=5, verbose=false, print_level=1) + @test prob.trajectory.data != initial + end + + @testset "EmbeddedOperator: leakage constraint" begin + prob = UnitarySmoothPulseProblem( + sys, U_goal, T, Δt; + da_bound=1.0, + piccolo_options=PiccoloOptions( + leakage_constraint=true, + leakage_constraint_value=5e-2, + leakage_cost=1e-1, + verbose=false + ) + ) + initial = copy(prob.trajectory.data) + solve!(prob, max_iter=5, verbose=false, print_level=1) + @test prob.trajectory.data != initial + end end # TODO: Test changing names of control, state, and timestep diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl index 43dd666c..e17dc460 100644 --- a/src/problem_templates/unitary_variational_problem.jl +++ b/src/problem_templates/unitary_variational_problem.jl @@ -81,8 +81,8 @@ function UnitaryVariationalProblem( 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), + Δt_min::Float64=0.5 * minimum(Δt), + Δt_max::Float64=2.0 * maximum(Δt), Q::Float64=100.0, Q_s::Float64=1e-2, Q_r::Float64=100.0, @@ -180,9 +180,12 @@ function UnitaryVariationalProblem( 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 + J += apply_piccolo_options!( + piccolo_options, constraints, traj; + state_names=state_name, + state_leakage_indices=goal isa EmbeddedOperator ? + get_iso_vec_leakage_indices(goal) : + nothing ) integrators = [ diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index 7ffc4d99..dbb23fe2 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -94,8 +94,8 @@ function LeakageConstraint( traj::NamedTrajectory; times=1:traj.T, ) - leakage_constraint(x) = [sum(abs2.(x[indices])) - value] - + leakage_constraint(x) = abs2.(x[indices]) .- value + return NonlinearKnotPointConstraint( leakage_constraint, name, diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index d2cd761d..f7da6c99 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -165,7 +165,7 @@ function LeakageObjective( times=1:traj.T, Qs::AbstractVector{<:Float64}=fill(1.0, length(times)), ) - leakage_objective(x) = sum(abs2.(x[indices])) + leakage_objective(x) = sum(abs2, x[indices]) / length(indices) return KnotPointObjective( leakage_objective, From 456fc30ba725c907f168df4c7cc18555397f00fc Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Mon, 23 Jun 2025 12:32:25 -0500 Subject: [PATCH 12/12] piccolo options docstring --- src/piccolo_options.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/piccolo_options.jl b/src/piccolo_options.jl index d51aa647..5a654426 100644 --- a/src/piccolo_options.jl +++ b/src/piccolo_options.jl @@ -4,6 +4,9 @@ export PiccoloOptions using ExponentialAction + +# TODO: Add duration and symmetry options + """ PiccoloOptions @@ -17,9 +20,10 @@ Options for the Piccolo quantum optimal control library. - `zero_initial_and_final_derivative::Bool=false`: Zero the initial and final control pulse derivatives. - `complex_control_norm_constraint_name::Union{Nothing, Symbol} = nothing`: Name of the complex control norm constraint. - `complex_control_norm_constraint_radius::Float64 = 1.0`: Radius of the complex control norm constraint. -- `bound_state::Bool = false`: Bound the state. -- `leakage_suppression::Bool = false`: Suppress leakage. -- `leakage_cost::Float64 = 1.0`: Leakage suppression parameter. +- `bound_state::Bool = false`: Bound the state variables <= 1.0. +- `leakage_constraint::Bool = false`: Suppress leakage with constraint and cost. +- `leakage_constraint_value::Float64 = 1e-2`: Value for the leakage constraint. +- `leakage_cost::Float64 = 1e-2`: Leakage suppression parameter. """ @kwdef mutable struct PiccoloOptions verbose::Bool = true