From a838da5bdb6c5353b5006ec91e015a25974247a2 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 17 Jun 2024 13:02:27 -0500 Subject: [PATCH 01/35] Objectives scalar multiplication --- src/objectives.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/objectives.jl b/src/objectives.jl index 4dbe90c0..261e429a 100644 --- a/src/objectives.jl +++ b/src/objectives.jl @@ -89,6 +89,21 @@ function Objective(terms::Vector{Dict}) return +(Objective.(terms)...) end +function Base.:*(num::Real, obj::Objective) + L = (Z⃗, Z) -> num * obj.L(Z⃗, Z) + ∇L = (Z⃗, Z) -> num * obj.∇L(Z⃗, Z) + if isnothing(obj.∂²L) + ∂²L = nothing + ∂²L_structure = nothing + else + ∂²L = (Z⃗, Z) -> num * obj.∂²L(Z⃗, Z) + ∂²L_structure = obj.∂²L_structure + end + return Objective(L, ∇L, ∂²L, ∂²L_structure, obj.terms) +end + +Base.:*(obj::Objective, num::Real) = num * obj + function Objective(term::Dict) return eval(term[:type])(; delete!(term, :type)...) end From dce1c0bc10aa70a213260ba4ec8d7dc5cd40cbcb Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 17 Jun 2024 13:03:09 -0500 Subject: [PATCH 02/35] Plotting layout helper --- src/plotting.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/plotting.jl b/src/plotting.jl index 3f3594a9..bd2d7eed 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -11,6 +11,11 @@ using ..Problems pretty_print(X::AbstractMatrix) = Base.show(stdout,"text/plain", X) +function get_layout(index::Int, n::Int) + √n = isqrt(n) + 1 + return ((index - 1) ÷ √n + 1, (index - 1) % √n + 1) +end + # """ # plot_unitary_populations( # traj::NamedTrajectory; From 4771e74df7038cf1cb4ede5c14489e5c4fca0b95 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 17 Jun 2024 13:03:57 -0500 Subject: [PATCH 03/35] bug fix: verbose passed to constrain --- src/problems.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/problems.jl b/src/problems.jl index 9ed2ef78..3f405d52 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -98,7 +98,8 @@ function QuantumControlProblem( linear_constraints, n_dynamics_constraints, nonlinear_constraints, - n_variables + n_variables, + verbose=verbose ) variables = reshape(variables, traj.dim, traj.T) @@ -228,7 +229,8 @@ function initialize_optimizer!( linear_constraints::Vector{LinearConstraint}, n_dynamics_constraints::Int, nonlinear_constraints::Vector{NonlinearConstraint}, - n_variables::Int + n_variables::Int; + verbose::Bool=true ) nl_cons = fill( MOI.NLPBoundsPair(0.0, 0.0), @@ -258,7 +260,7 @@ function initialize_optimizer!( variables = MOI.add_variables(optimizer, n_variables) # add linear constraints - constrain!(optimizer, variables, linear_constraints, trajectory, verbose=true) + constrain!(optimizer, variables, linear_constraints, trajectory, verbose=verbose) return variables end From f8a3fa8a0e59e36ca2412ce0c02b2a4bf5946ec8 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 17 Jun 2024 13:04:31 -0500 Subject: [PATCH 04/35] quantum system is_reachable and is_reachable arg swap --- src/quantum_system_utils.jl | 33 +++++++++++++++++------- test/quantum_system_utils_test.jl | 43 ++++++++++++++++++++++++------- 2 files changed, 58 insertions(+), 18 deletions(-) diff --git a/src/quantum_system_utils.jl b/src/quantum_system_utils.jl index be9a2946..8c6e1eb5 100644 --- a/src/quantum_system_utils.jl +++ b/src/quantum_system_utils.jl @@ -4,6 +4,7 @@ export operator_algebra export is_reachable using ..EmbeddedOperators +using ..QuantumSystems using LinearAlgebra using SparseArrays @@ -154,8 +155,8 @@ function operator_algebra( end function fit_gen_to_basis( - basis::AbstractVector{<:AbstractMatrix{<:T}}, - gen::AbstractMatrix{<:T} + gen::AbstractMatrix{<:T}, + basis::AbstractVector{<:AbstractMatrix{<:T}} ) where T<:Number A = stack(vec.(basis)) b = vec(gen) @@ -163,8 +164,8 @@ function fit_gen_to_basis( end function is_in_span( - basis::AbstractVector{<:AbstractMatrix}, - gen::AbstractMatrix; + gen::AbstractMatrix, + basis::AbstractVector{<:AbstractMatrix}; subspace_indices::AbstractVector{<:Int}=1:size(gen, 1), atol=eps(Float32), return_effective_gen=false, @@ -172,7 +173,7 @@ function is_in_span( g_basis = [deepcopy(b[subspace_indices, subspace_indices]) for b ∈ basis] linearly_independent_subset!(g_basis) # Traceless basis needs traceless fit - x = fit_gen_to_basis(g_basis, gen) + x = fit_gen_to_basis(gen, g_basis) g_eff = sum(x .* g_basis) ε = norm(g_eff - gen, 2) if return_effective_gen @@ -194,8 +195,8 @@ is_reachable(hamiltonians, gate, subspace_levels, levels; kwargs...) - `levels::Vector{<:Int}`: levels of the system """ function is_reachable( - hamiltonians::AbstractVector{<:AbstractMatrix}, - gate::AbstractMatrix; + gate::AbstractMatrix, + hamiltonians::AbstractVector{<:AbstractMatrix}; subspace_indices::AbstractVector{<:Int}=1:size(gate, 1), compute_basis=true, remove_trace=true, @@ -215,12 +216,26 @@ function is_reachable( end return is_in_span( - basis, generator, + basis, subspace_indices=subspace_indices, atol=atol ) end - +function is_reachable( + gate::AbstractMatrix, + system::QuantumSystem; + use_drift::Bool=true, + kwargs... +) + if !iszero(system.H_drift) && use_drift + hamiltonians = [system.H_drift, system.H_drives...] + else + hamiltonians = system.H_drives + end + return is_reachable(gate, hamiltonians; kwargs...) end + + +end \ No newline at end of file diff --git a/test/quantum_system_utils_test.jl b/test/quantum_system_utils_test.jl index 83f1994c..650e5170 100644 --- a/test/quantum_system_utils_test.jl +++ b/test/quantum_system_utils_test.jl @@ -63,7 +63,17 @@ end # Check 1 qubit with complete basis gen = map(A -> kron_from_dict(A, H_ops), ["X", "Y"]) target = H_ops["Z"] - @test is_reachable(gen, target, compute_basis=true, verbose=false) + @test is_reachable(target, gen, compute_basis=true, verbose=false) + + # System + sys = QuantumSystem([GATES[:X], GATES[:Y], GATES[:Z]]) + target = GATES[:Z] + @test is_reachable(target, sys) + + # System with drift + sys = QuantumSystem(GATES[:Z], [GATES[:X]]) + target = GATES[:Z] + @test is_reachable(target, sys) # Check 2 qubit with complete basis XI = GATES[:X] ⊗ GATES[:I] @@ -85,16 +95,31 @@ end CX = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0] # Pass - @test is_reachable(complete_gen, R2) - @test is_reachable(complete_gen, CZ) - @test is_reachable(complete_gen, CX) - @test is_reachable(complete_gen, XI) + @test is_reachable(R2, complete_gen) + @test is_reachable(CZ, complete_gen) + @test is_reachable(CX, complete_gen) + @test is_reachable(XI, complete_gen) + + # Mostly fail + @test !is_reachable(R2, incomplete_gen) + @test !is_reachable(CZ, incomplete_gen) + @test !is_reachable(CX, incomplete_gen) + @test is_reachable(XI, incomplete_gen) + + # QuantumSystems + complete_gen_sys = QuantumSystem(complete_gen) + incomplete_gen_sys = QuantumSystem(incomplete_gen) + # Pass + @test is_reachable(R2, complete_gen_sys) + @test is_reachable(CZ, complete_gen_sys) + @test is_reachable(CX, complete_gen_sys) + @test is_reachable(XI, complete_gen_sys) # Mostly fail - @test !is_reachable(incomplete_gen, R2) - @test !is_reachable(incomplete_gen, CZ) - @test !is_reachable(incomplete_gen, CX) - @test is_reachable(incomplete_gen, XI) + @test !is_reachable(R2, incomplete_gen_sys) + @test !is_reachable(CZ, incomplete_gen_sys) + @test !is_reachable(CX, incomplete_gen_sys) + @test is_reachable(XI, incomplete_gen_sys) end @testitem "Lie Algebra subspace reachability" begin From a71f36029f242dfde9fd4c509b5f7a2f8b97431d Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 19:52:59 -0500 Subject: [PATCH 05/35] add trajectory initilization helpers --- .../unitary_smooth_pulse_problem.jl | 215 ++++++------------ src/trajectory_initialization.jl | 152 +++++++++++++ 2 files changed, 221 insertions(+), 146 deletions(-) diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 0db87bd1..f4892372 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -117,14 +117,6 @@ function UnitarySmoothPulseProblem( blas_multithreading=true, kwargs... ) - if operator isa EmbeddedOperator - U_goal = operator.operator - U_init = get_subspace_identity(operator) - else - U_goal = Matrix{ComplexF64}(operator) - U_init = Matrix{ComplexF64}(I(size(U_goal, 1))) - end - if !blas_multithreading BLAS.set_num_threads(1) end @@ -133,127 +125,43 @@ function UnitarySmoothPulseProblem( ipopt_options.hessian_approximation = "limited-memory" end + # Goal + if operator isa EmbeddedOperator + U_goal = operator.operator + U_init = get_subspace_identity(operator) + else + U_goal = Matrix{ComplexF64}(operator) + U_init = Matrix{ComplexF64}(I(size(U_goal, 1))) + end + Ũ⃗_init = operator_to_iso_vec(U_init) + Ũ⃗_goal = operator_to_iso_vec(U_goal) n_drives = length(system.G_drives) + # Trajectory if !isnothing(init_trajectory) traj = init_trajectory else - if free_time - if Δt isa Float64 - Δt = fill(Δt, 1, T) - end - end - - if isnothing(a_guess) - if geodesic - if operator isa EmbeddedOperator - Ũ⃗ = unitary_geodesic(operator, T) - else - Ũ⃗ = unitary_geodesic(U_goal, T) - end - else - Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, T) - end - - if a_bounds isa AbstractVector - a_dists = [Uniform(-a_bounds[i], a_bounds[i]) for i = 1:n_drives] - elseif a_bounds isa Tuple - a_dists = [Uniform(aᵢ_lb, aᵢ_ub) for (aᵢ_lb, aᵢ_ub) ∈ zip(a_bounds...)] - else - error("a_bounds must be a Vector or Tuple") - end - - a = hcat([ - zeros(n_drives), - vcat([rand(a_dists[i], 1, T - 2) for i = 1:n_drives]...), - zeros(n_drives) - ]...) - - da = randn(n_drives, T) * drive_derivative_σ - dda = randn(n_drives, T) * drive_derivative_σ - else - Ũ⃗ = unitary_rollout( - operator_to_iso_vec(U_init), - a_guess, - Δt, - system; - integrator=rollout_integrator - ) - a = a_guess - da = derivative(a, Δt) - dda = derivative(da, Δt) - - # to avoid constraint violation error at initial iteration - da[:, end] = da[:, end-1] + Δt[end-1] * dda[:, end-1] - end - - initial = ( - Ũ⃗ = operator_to_iso_vec(U_init), - a = zeros(n_drives), - ) - - final = ( - a = zeros(n_drives), + traj = initialize_trajectory( + Ũ⃗_init, + Ũ⃗_goal, + T, + Δt, + n_drives, + a_bounds, + dda_bounds; + free_time=free_time, + Δt_bounds=(Δt_min, Δt_max), + geodesic=geodesic, + bound_unitary=bound_unitary, + drive_derivative_σ=drive_derivative_σ, + a_guess=a_guess, + system=system, + rollout_integrator=rollout_integrator, ) - - goal = ( - Ũ⃗ = operator_to_iso_vec(U_goal), - ) - - bounds = ( - a = a_bounds, - dda = dda_bounds, - ) - - if bound_unitary - Ũ⃗_dim = size(Ũ⃗, 1) - Ũ⃗_bound = ( - Ũ⃗ = (-ones(Ũ⃗_dim), ones(Ũ⃗_dim)), - ) - bounds = merge(bounds, Ũ⃗_bound) - end - - if free_time - components = ( - Ũ⃗ = Ũ⃗, - a = a, - da = da, - dda = dda, - Δt = Δt, - ) - - bounds = merge(bounds, (Δt = (Δt_min, Δt_max),)) - - traj = NamedTrajectory( - components; - controls=(:dda, :Δt), - timestep=:Δt, - bounds=bounds, - initial=initial, - final=final, - goal=goal - ) - else - components = ( - Ũ⃗ = Ũ⃗, - a = a, - da = da, - dda = dda, - ) - - traj = NamedTrajectory( - components; - controls=(:dda,), - timestep=Δt, - bounds=bounds, - initial=initial, - final=final, - goal=goal - ) - end end + # Objective J = UnitaryInfidelityObjective(:Ũ⃗, traj, Q; subspace=operator isa EmbeddedOperator ? operator.subspace_indices : nothing, ) @@ -261,39 +169,17 @@ function UnitarySmoothPulseProblem( J += QuadraticRegularizer(:da, traj, R_da) J += QuadraticRegularizer(:dda, traj, R_dda) + # Constraints if leakage_suppression if operator isa EmbeddedOperator - leakage_indices = get_unitary_isomorphism_leakage_indices(operator) - J_leakage, slack_con = L1Regularizer( - :Ũ⃗, - traj; - R_value=R_leakage, - indices=leakage_indices + set_leakage_suppression!( + J, constraints, traj, operator, Ũ⃗_symb=Ũ⃗ᵢ, R=R_leakage ) - push!(constraints, slack_con) - J += J_leakage else @warn "leakage_suppression is not supported for non-embedded operators, ignoring." end end - if integrator == :pade - unitary_integrator = - UnitaryPadeIntegrator(system, :Ũ⃗, :a; order=pade_order, autodiff=autodiff) - elseif integrator == :exponential - unitary_integrator = - UnitaryExponentialIntegrator(system, :Ũ⃗, :a) - else - error("integrator must be one of (:pade, :exponential)") - end - - - integrators = [ - unitary_integrator, - DerivativeIntegrator(:a, :da, traj), - DerivativeIntegrator(:da, :dda, traj), - ] - if free_time if timesteps_all_equal push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) @@ -312,6 +198,23 @@ function UnitarySmoothPulseProblem( push!(constraints, norm_con) end + # Integrators + if integrator == :pade + unitary_integrator = + UnitaryPadeIntegrator(system, :Ũ⃗, :a; order=pade_order, autodiff=autodiff) + elseif integrator == :exponential + unitary_integrator = + UnitaryExponentialIntegrator(system, :Ũ⃗, :a) + else + error("integrator must be one of (:pade, :exponential)") + end + + integrators = [ + unitary_integrator, + DerivativeIntegrator(:a, :da, traj), + DerivativeIntegrator(:da, :dda, traj), + ] + return QuantumControlProblem( system, traj, @@ -330,7 +233,6 @@ function UnitarySmoothPulseProblem( end - """ UnitarySmoothPulseProblem( H_drift::AbstractMatrix{<:Number}, @@ -353,6 +255,27 @@ function UnitarySmoothPulseProblem( return UnitarySmoothPulseProblem(system, args...; kwargs...) end + +function set_leakage_suppression!( + J::Objective, + constraints::Vector{<:AbstractConstraint}, + traj::NamedTrajectory, + operator::EmbeddedOperator; + R::Real=1e-1, + Ũ⃗_symb::Symbol=:Ũ⃗, +) + # leakage + leakage_indices = get_unitary_isomorphism_leakage_indices(operator) + J_leakage, slack_con = L1Regularizer( + traj[Ũ⃗_symb], + traj; + R_value=R, + indices=leakage_indices + ) + push!(constraints, slack_con) + J += J_leakage +end + # *************************************************************************** # @testitem "Hadamard gate" begin diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 7bde8e9b..b143b07b 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -3,8 +3,11 @@ module TrajectoryInitialization export unitary_linear_interpolation export unitary_geodesic export linear_interpolation +export initialize_trajectory +using NamedTrajectories using LinearAlgebra +using Distributions using ..QuantumUtils using ..EmbeddedOperators @@ -128,4 +131,153 @@ end linear_interpolation(x::AbstractVector, y::AbstractVector, n::Int) = hcat(range(x, y, n)...) +# ============================================================================= + +VectorBound = Union{AbstractVector{R}, Tuple{AbstractVector{R}, AbstractVector{R}}} where R <: Real +ScalarBound = Union{R, Tuple{R, R}} where R <: Real + +function initialize_unitaries( + Ũ⃗_init::AbstractVector{<:Number}, + Ũ⃗_goal::AbstractVector{<:Number}, + T::Int; + geodesic=true +) + U_init = iso_vec_to_operator(Ũ⃗_init) + U_goal = iso_vec_to_operator(Ũ⃗_goal) + if geodesic + Ũ⃗ = unitary_geodesic(U_goal, T) + else + Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, T) + end + return Ũ⃗ +end + +function initialize_controls( + n_drives::Int, + T::Int, + a_bounds::VectorBound, + drive_derivative_σ::Float64 +) + if a_bounds isa AbstractVector + a_dists = [Uniform(-a_bounds[i], a_bounds[i]) for i = 1:n_drives] + elseif a_bounds isa Tuple + a_dists = [Uniform(aᵢ_lb, aᵢ_ub) for (aᵢ_lb, aᵢ_ub) ∈ zip(a_bounds...)] + else + error("a_bounds must be a Vector or Tuple") + end + + a = hcat([ + zeros(n_drives), + vcat([rand(a_dists[i], 1, T - 2) for i = 1:n_drives]...), + zeros(n_drives) + ]...) + + da = randn(n_drives, T) * drive_derivative_σ + dda = randn(n_drives, T) * drive_derivative_σ + return a, da, dda +end + +function initialize_controls( + a_guess::AbstractMatrix, +) + a = a_guess + da = derivative(a, Δt) + dda = derivative(da, Δt) + + # to avoid constraint violation error at initial iteration + da[:, end] = da[:, end-1] + Δt[end-1] * dda[:, end-1] + + return a, da, dda +end + +function initialize_trajectory( + Ũ⃗_init::AbstractVector{<:Number}, + Ũ⃗_goal::AbstractVector{<:Number}, + T::Int, + Δt::Real, + n_drives::Int, + a_bounds::VectorBound, + dda_bounds::VectorBound; + geodesic=true, + bound_unitary=false, + free_time=false, + Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt), + drive_derivative_σ::Float64=0.1, + a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + system=Union{AbstractQuantumSystem, Nothing}=nothing, + rollout_integrator::Function=exp, + Ũ⃗_keys::AbstractVector{<:Symbol}=[:Ũ⃗], +) + if free_time + if Δt isa Float64 + Δt = fill(Δt, 1, T) + end + end + + # Initial state and controls + if isnothing(a_guess) + Ũ⃗ = initialize_unitaries(Ũ⃗_init, Ũ⃗_goal, T; geodesic=geodesic) + a, da, dda = initialize_controls(n_drives, T, a_bounds, drive_derivative_σ) + else + @assert !isnothing(system) "system must be provided if a_guess is provided" + Ũ⃗ = unitary_rollout( + Ũ⃗_init, + a_guess, + Δt, + system; + integrator=rollout_integrator + ) + a, da, dda = initialize_controls(a_guess) + end + + # Constraints + Ũ⃗_inits = repeat([Ũ⃗_init], length(Ũ⃗_keys)) + initial = (; + (Ũ⃗_keys .=> Ũ⃗_inits)..., + a = zeros(n_drives), + ) + + final = ( + a = zeros(n_drives), + ) + + Ũ⃗_goals = repeat([Ũ⃗_goal], length(Ũ⃗_keys)) + goal = (; (Ũ⃗_keys .=> Ũ⃗_goals)...) + + # Bounds + bounds = (a = a_bounds, dda = dda_bounds,) + + if bound_unitary + Ũ⃗_dim = length(Ũ⃗_init) + Ũ⃗_bounds = repeat([(-ones(Ũ⃗_dim), ones(Ũ⃗_dim))], length(Ũ⃗_keys)) + merge!(bounds, (; (Ũ⃗_keys .=> Ũ⃗_bounds)...)) + end + + # Trajectory + Ũ⃗_values = repeat([Ũ⃗], length(Ũ⃗_keys)) + keys = [Ũ⃗_keys..., :a, :da, :dda] + values = [Ũ⃗_values..., a, da, dda] + + if free_time + push!(keys, :Δt) + push!(values, Δt) + controls = (:dda, :Δt) + timestep = :Δt + bounds = merge(bounds, (Δt = Δt_bounds,)) + else + controls = (:dda,) + timestep = Δt + end + + return NamedTrajectory( + (; (keys .=> values)...); + controls=controls, + timestep=timestep, + bounds=bounds, + initial=initial, + final=final, + goal=goal + ) +end + end From 6c6feb1eb56e0ff8cfd9256684ba3d23058fa2ee Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 20:20:06 -0500 Subject: [PATCH 06/35] pass verbose to constrain --- src/problems.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/problems.jl b/src/problems.jl index 9ed2ef78..8e4a867d 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -98,7 +98,8 @@ function QuantumControlProblem( linear_constraints, n_dynamics_constraints, nonlinear_constraints, - n_variables + n_variables, + verbose=verbose ) variables = reshape(variables, traj.dim, traj.T) @@ -228,7 +229,8 @@ function initialize_optimizer!( linear_constraints::Vector{LinearConstraint}, n_dynamics_constraints::Int, nonlinear_constraints::Vector{NonlinearConstraint}, - n_variables::Int + n_variables::Int; + verbose=true ) nl_cons = fill( MOI.NLPBoundsPair(0.0, 0.0), @@ -258,7 +260,7 @@ function initialize_optimizer!( variables = MOI.add_variables(optimizer, n_variables) # add linear constraints - constrain!(optimizer, variables, linear_constraints, trajectory, verbose=true) + constrain!(optimizer, variables, linear_constraints, trajectory, verbose=verbose) return variables end From e05619d6c19069293b2269f33c50914ba9b27963 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 21:32:39 -0500 Subject: [PATCH 07/35] fix embed op usage in geodesic --- src/embedded_operators.jl | 8 +++++ src/trajectory_initialization.jl | 54 ++++++++++++++------------------ 2 files changed, 32 insertions(+), 30 deletions(-) diff --git a/src/embedded_operators.jl b/src/embedded_operators.jl index 2f08d388..17c3235e 100644 --- a/src/embedded_operators.jl +++ b/src/embedded_operators.jl @@ -54,10 +54,18 @@ end EmbeddedOperator(op::Matrix{<:Number}, subspace_indices::AbstractVector{Int}, levels::Int) = EmbeddedOperator(op, subspace_indices, [levels]) +function embed(matrix::Matrix{ComplexF64}, op::EmbeddedOperator) + return embed(matrix, op.subspace_indices, prod(op.subsystem_levels)) +end + function unembed(op::EmbeddedOperator)::Matrix{ComplexF64} return op.operator[op.subspace_indices, op.subspace_indices] end +function unembed(matrix::AbstractMatrix, op::EmbeddedOperator) + return matrix[op.subspace_indices, op.subspace_indices] +end + Base.size(op::EmbeddedOperator) = size(op.operator) Base.size(op::EmbeddedOperator, dim::Union{Int, Nothing}) = size(op.operator, dim) diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index b143b07b..b08d9e16 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -65,38 +65,27 @@ Compute a geodesic connecting two unitary operators. function unitary_geodesic end function unitary_geodesic( - operator::EmbeddedOperator, + U_init::AbstractMatrix{<:Number}, + U_goal::EmbeddedOperator, samples::Int; kwargs... -) - U_goal = unembed(operator) - U_init = Matrix{ComplexF64}(I(size(U_goal, 1))) - Ũ⃗ = unitary_geodesic(U_init, U_goal, samples; 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(Ũ⃗ₜ), - operator.subspace_indices, - prod(operator.subsystem_levels)) - ) for Ũ⃗ₜ ∈ eachcol(Ũ⃗)]...) + 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... -) - N = size(U_goal, 1) - U₀ = Matrix{ComplexF64}(I(N)) - return unitary_geodesic(U₀, U_goal, samples; kwargs...) -end - -function unitary_geodesic( - U₀::AbstractMatrix{<:Number}, - U₁::AbstractMatrix{<:Number}, samples::Number; kwargs... ) - return unitary_geodesic(U₀, U₁, range(0, 1, samples); kwargs...) + return unitary_geodesic(U_init, U_goal, range(0, 1, samples); kwargs...) end function unitary_geodesic( @@ -137,15 +126,13 @@ VectorBound = Union{AbstractVector{R}, Tuple{AbstractVector{R}, AbstractVector{R ScalarBound = Union{R, Tuple{R, R}} where R <: Real function initialize_unitaries( - Ũ⃗_init::AbstractVector{<:Number}, - Ũ⃗_goal::AbstractVector{<:Number}, + U_init::AbstractMatrix{<:Number}, + U_goal::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, T::Int; geodesic=true ) - U_init = iso_vec_to_operator(Ũ⃗_init) - U_goal = iso_vec_to_operator(Ũ⃗_goal) if geodesic - Ũ⃗ = unitary_geodesic(U_goal, T) + Ũ⃗ = unitary_geodesic(U_init, U_goal, T) else Ũ⃗ = unitary_linear_interpolation(U_init, U_goal, T) end @@ -191,13 +178,13 @@ function initialize_controls( end function initialize_trajectory( - Ũ⃗_init::AbstractVector{<:Number}, - Ũ⃗_goal::AbstractVector{<:Number}, + U_goal::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, T::Int, Δt::Real, n_drives::Int, a_bounds::VectorBound, dda_bounds::VectorBound; + U_init::AbstractMatrix{<:Number}=Matrix{ComplexF64}(I(size(U_goal, 1))), geodesic=true, bound_unitary=false, free_time=false, @@ -214,9 +201,16 @@ function initialize_trajectory( end end + Ũ⃗_init = operator_to_iso_vec(U_init) + if U_goal isa EmbeddedOperator + Ũ⃗_goal = operator_to_iso_vec(U_goal.operator) + else + Ũ⃗_goal = operator_to_iso_vec(U_goal) + end + # Initial state and controls if isnothing(a_guess) - Ũ⃗ = initialize_unitaries(Ũ⃗_init, Ũ⃗_goal, T; geodesic=geodesic) + Ũ⃗ = initialize_unitaries(U_init, U_goal, T, geodesic=geodesic) a, da, dda = initialize_controls(n_drives, T, a_bounds, drive_derivative_σ) else @assert !isnothing(system) "system must be provided if a_guess is provided" From e1c04d89d360912b8cce3cc4187060bb9947e578 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 21:32:50 -0500 Subject: [PATCH 08/35] leakage test --- src/objectives.jl | 11 +++ .../unitary_smooth_pulse_problem.jl | 95 ++++++++++--------- 2 files changed, 62 insertions(+), 44 deletions(-) diff --git a/src/objectives.jl b/src/objectives.jl index 4dbe90c0..c4b07d5d 100644 --- a/src/objectives.jl +++ b/src/objectives.jl @@ -15,6 +15,7 @@ export QuadraticRegularizer export PairwiseQuadraticRegularizer export QuadraticSmoothnessRegularizer export L1Regularizer +export L1Regularizer! using TrajectoryIndexingUtils using ..QuantumUtils @@ -849,6 +850,16 @@ function L1Regularizer( return J, slack_con end +function L1Regularizer!( + constraints::Vector{<:AbstractConstraint}, + name::Symbol, + traj::NamedTrajectory; + kwargs... +) + J, slack_con = L1Regularizer(name, traj; kwargs...) + push!(constraints, slack_con) + return J +end function MinimumTimeObjective(; D::Float64=1.0, diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index f4892372..85707dbb 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -125,26 +125,13 @@ function UnitarySmoothPulseProblem( ipopt_options.hessian_approximation = "limited-memory" end - # Goal - if operator isa EmbeddedOperator - U_goal = operator.operator - U_init = get_subspace_identity(operator) - else - U_goal = Matrix{ComplexF64}(operator) - U_init = Matrix{ComplexF64}(I(size(U_goal, 1))) - end - - Ũ⃗_init = operator_to_iso_vec(U_init) - Ũ⃗_goal = operator_to_iso_vec(U_goal) - n_drives = length(system.G_drives) - # Trajectory if !isnothing(init_trajectory) traj = init_trajectory else + n_drives = length(system.G_drives) traj = initialize_trajectory( - Ũ⃗_init, - Ũ⃗_goal, + operator, T, Δt, n_drives, @@ -172,8 +159,12 @@ function UnitarySmoothPulseProblem( # Constraints if leakage_suppression if operator isa EmbeddedOperator - set_leakage_suppression!( - J, constraints, traj, operator, Ũ⃗_symb=Ũ⃗ᵢ, R=R_leakage + leakage_indices = get_unitary_isomorphism_leakage_indices(operator) + J += L1Regularizer!( + constraints, :Ũ⃗, traj, + R_value=R_leakage, + indices=leakage_indices, + eval_hessian=!hessian_approximation ) else @warn "leakage_suppression is not supported for non-embedded operators, ignoring." @@ -255,41 +246,57 @@ function UnitarySmoothPulseProblem( return UnitarySmoothPulseProblem(system, args...; kwargs...) end - -function set_leakage_suppression!( - J::Objective, - constraints::Vector{<:AbstractConstraint}, - traj::NamedTrajectory, - operator::EmbeddedOperator; - R::Real=1e-1, - Ũ⃗_symb::Symbol=:Ũ⃗, -) - # leakage - leakage_indices = get_unitary_isomorphism_leakage_indices(operator) - J_leakage, slack_con = L1Regularizer( - traj[Ũ⃗_symb], - traj; - R_value=R, - indices=leakage_indices - ) - push!(constraints, slack_con) - J += J_leakage -end - # *************************************************************************** # @testitem "Hadamard gate" begin - H_drift = GATES[:Z] - H_drives = [GATES[:X], GATES[:Y]] + sys = QuantumSystem(GATES[:Z], [GATES[:X], GATES[:Y]]) U_goal = GATES[:H] T = 51 Δt = 0.2 + + prob = UnitarySmoothPulseProblem( + sys, U_goal, T, Δt, + ipopt_options=Options(print_level=1), verbose=false + ) + + initial = unitary_fidelity(prob) + solve!(prob, max_iter=20) + final = unitary_fidelity(prob) + @test final > initial +end + +@testitem "EmbeddedOperator Hadamard gate" begin + a = annihilate(3) + sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) + U_goal = EmbeddedOperator(GATES[:H], sys) + T = 51 + Δt = 0.2 + + prob = UnitarySmoothPulseProblem( + sys, U_goal, T, Δt, + ipopt_options=Options(print_level=1), verbose=false + ) + + initial = unitary_fidelity(prob, subspace=U_goal.subspace_indices) + solve!(prob, max_iter=20) + final = unitary_fidelity(prob, subspace=U_goal.subspace_indices) + @test final > initial + + # Test leakage suppression + a = annihilate(4) + sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) + U_goal = EmbeddedOperator(GATES[:H], sys) + T = 50 + Δt = 0.2 prob = UnitarySmoothPulseProblem( - H_drift, H_drives, U_goal, T, Δt, - ipopt_options=Options(print_level=4) + sys, U_goal, T, Δt, + leakage_suppression=true, R_leakage=1e-1, + ipopt_options=Options(print_level=1), verbose=false ) - solve!(prob, max_iter=100) - @test unitary_fidelity(prob) > 0.99 + initial = unitary_fidelity(prob, subspace=U_goal.subspace_indices) + solve!(prob, max_iter=20) + final = unitary_fidelity(prob, subspace=U_goal.subspace_indices) + @test final > initial end From 40d9bede4aac9f0fae093ba7a6c7dba28911f171 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 22:04:39 -0500 Subject: [PATCH 09/35] scalar obj multiplication --- src/objectives.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/objectives.jl b/src/objectives.jl index c4b07d5d..f8be8db7 100644 --- a/src/objectives.jl +++ b/src/objectives.jl @@ -90,6 +90,21 @@ function Objective(terms::Vector{Dict}) return +(Objective.(terms)...) end +function Base.:*(num::Real, obj::Objective) + L = (Z⃗, Z) -> num * obj.L(Z⃗, Z) + ∇L = (Z⃗, Z) -> num * obj.∇L(Z⃗, Z) + if isnothing(obj.∂²L) + ∂²L = nothing + ∂²L_structure = nothing + else + ∂²L = (Z⃗, Z) -> num * obj.∂²L(Z⃗, Z) + ∂²L_structure = obj.∂²L_structure + end + return Objective(L, ∇L, ∂²L, ∂²L_structure, obj.terms) +end + +Base.:*(obj::Objective, num::Real) = num * obj + function Objective(term::Dict) return eval(term[:type])(; delete!(term, :type)...) end From 9f07cdd3791fde95224da51243d9195ce57fa0f7 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 22:29:58 -0500 Subject: [PATCH 10/35] sampling problem template --- src/problem_templates/_problem_templates.jl | 2 + .../unitary_sampling_problem.jl | 180 ++++++++++++++++++ 2 files changed, 182 insertions(+) create mode 100644 src/problem_templates/unitary_sampling_problem.jl diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 4d90e87f..65aa9801 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -4,6 +4,7 @@ export UnitarySmoothPulseProblem export UnitaryMinimumTimeProblem export UnitaryRobustnessProblem export UnitaryDirectSumProblem +export UnitarySamplingProblem export QuantumStateSmoothPulseProblem export QuantumStateMinimumTimeProblem @@ -32,6 +33,7 @@ include("unitary_smooth_pulse_problem.jl") include("unitary_minimum_time_problem.jl") include("unitary_robustness_problem.jl") include("unitary_direct_sum_problem.jl") +include("unitary_sampling_problem.jl") include("quantum_state_smooth_pulse_problem.jl") include("quantum_state_minimum_time_problem.jl") diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl new file mode 100644 index 00000000..be551d91 --- /dev/null +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -0,0 +1,180 @@ +@doc raw""" + UnitarySamplingProblem + + TODO: systems might need flexible bounds. + +""" +function UnitarySamplingProblem( + systems::AbstractVector{<:AbstractQuantumSystem}, + operator::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, + T::Int, + Δt::Union{Float64, Vector{Float64}}; + system_labels=string.(1:length(systems)), + system_weights=fill(1.0, length(systems)), + free_time=true, + init_trajectory::Union{NamedTrajectory, Nothing}=nothing, + a_bound::Float64=1.0, + a_bounds=fill(a_bound, length(systems[1].G_drives)), + a_guess::Union{Matrix{Float64}, Nothing}=nothing, + dda_bound::Float64=1.0, + dda_bounds=fill(dda_bound, length(systems[1].G_drives)), + Δt_min::Float64=0.5 * Δt, + Δt_max::Float64=1.5 * Δt, + drive_derivative_σ::Float64=0.01, + Q::Float64=100.0, + R=1e-2, + R_a::Union{Float64, Vector{Float64}}=R, + R_da::Union{Float64, Vector{Float64}}=R, + R_dda::Union{Float64, Vector{Float64}}=R, + leakage_suppression=false, + R_leakage=1e-1, + max_iter::Int=1000, + linear_solver::String="mumps", + ipopt_options::Options=Options(), + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], + timesteps_all_equal::Bool=true, + verbose::Bool=false, + integrator::Symbol=:pade, + rollout_integrator=exp, + bound_unitary=integrator == :exponential, + control_norm_constraint=false, + control_norm_constraint_components=nothing, + control_norm_R=nothing, + geodesic=true, + pade_order=4, + autodiff=pade_order != 4, + jacobian_structure=true, + hessian_approximation=false, + blas_multithreading=true, + kwargs... +) + if !blas_multithreading + BLAS.set_num_threads(1) + end + + if hessian_approximation + ipopt_options.hessian_approximation = "limited-memory" + end + + # Create keys for multiple systems + Ũ⃗_keys = [add_suffix(:Ũ⃗, ℓ) for ℓ ∈ system_labels] + + # Trajectory + if !isnothing(init_trajectory) + traj = init_trajectory + else + n_drives = length(systems[1].G_drives) + # TODO: Initial system? + traj = initialize_trajectory( + operator, + T, + Δt, + n_drives, + a_bounds, + dda_bounds; + free_time=free_time, + Δt_bounds=(Δt_min, Δt_max), + geodesic=geodesic, + bound_unitary=bound_unitary, + drive_derivative_σ=drive_derivative_σ, + a_guess=a_guess, + system=systems[1], + rollout_integrator=rollout_integrator, + Ũ⃗_keys=Ũ⃗_keys + ) + end + + # Objective + J = NullObjective() + for (wᵢ, Ũ⃗_key) in zip(system_weights, Ũ⃗_keys) + J += wᵢ * UnitaryInfidelityObjective( + Ũ⃗_key, traj, Q; + subspace=operator isa EmbeddedOperator ? operator.subspace_indices : nothing + ) + end + J += QuadraticRegularizer(:a, traj, R_a) + J += QuadraticRegularizer(:da, traj, R_da) + J += QuadraticRegularizer(:dda, traj, R_dda) + + # Constraints + if leakage_suppression + if operator isa EmbeddedOperator + leakage_indices = get_unitary_isomorphism_leakage_indices(operator) + for Ũ⃗_key in Ũ⃗_keys + J += L1Regularizer!( + constraints, Ũ⃗_key, traj, + R_value=R_leakage, + indices=leakage_indices, + eval_hessian=!hessian_approximation + ) + end + else + @warn "leakage_suppression is not supported for non-embedded operators, ignoring." + end + end + + if free_time + if timesteps_all_equal + push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) + end + end + + if control_norm_constraint + @assert !isnothing(control_norm_constraint_components) "control_norm_constraint_components must be provided" + @assert !isnothing(control_norm_R) "control_norm_R must be provided" + norm_con = ComplexModulusContraint( + :a, + control_norm_R, + traj; + name_comps=control_norm_constraint_components, + ) + push!(constraints, norm_con) + end + + # Integrators + unitary_integrators = AbstractIntegrator[] + for (sys, Ũ⃗_key) in zip(systems, Ũ⃗_keys) + if integrator == :pade + push!( + unitary_integrators, + UnitaryPadeIntegrator(sys, Ũ⃗_key, :a; order=pade_order, autodiff=autodiff) + ) + elseif integrator == :exponential + push!( + unitary_integrators, + UnitaryExponentialIntegrator(sys, Ũ⃗_key, :a) + ) + else + error("integrator must be one of (:pade, :exponential)") + end + end + + integrators = [ + unitary_integrators..., + DerivativeIntegrator(:a, :da, traj), + DerivativeIntegrator(:da, :dda, traj), + ] + + return QuantumControlProblem( + direct_sum(systems), + traj, + J, + integrators; + constraints=constraints, + max_iter=max_iter, + linear_solver=linear_solver, + verbose=verbose, + ipopt_options=ipopt_options, + jacobian_structure=jacobian_structure, + hessian_approximation=hessian_approximation, + eval_hessian=!hessian_approximation, + kwargs... + ) +end + + +# ζs = range(-.05, .05, length=5) +# systems = [system(ζ) for ζ ∈ ζs]; +# system_labels = string.(1:length(systems)) +# system_weights = fill(1.0, length(systems)) + From 2c6de399cd41f9a7ef436980a5e206036a204739 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 23:23:45 -0500 Subject: [PATCH 11/35] sample from distribution --- .../unitary_sampling_problem.jl | 24 +++++++++++++++---- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index be551d91..faea4008 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -172,9 +172,23 @@ function UnitarySamplingProblem( ) end - -# ζs = range(-.05, .05, length=5) -# systems = [system(ζ) for ζ ∈ ζs]; -# system_labels = string.(1:length(systems)) -# system_weights = fill(1.0, length(systems)) +function UnitarySamplingProblem( + system::Function, + distribution::Sampleable, + num_samples::Int, + operator::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, + T::Int, + Δt::Union{Float64, Vector{Float64}}; + kwargs... +) + samples = rand(distribution, num_samples) + systems = [system(x) for x in samples] + return UnitarySamplingProblem( + systems, + operator, + T, + Δt; + kwargs... + ) +end From 3e091643205773a0657730f184bd4e6e742da894 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 19 Jun 2024 23:50:44 -0500 Subject: [PATCH 12/35] test sampling robust problem --- .../unitary_sampling_problem.jl | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index faea4008..22c7ca6f 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -192,3 +192,50 @@ function UnitarySamplingProblem( ) end +# ============================================================================= + +@testitem "Sample robustness test" begin + using Distributions + + n_samples = 3 + T = 50 + Δt = 0.2 + timesteps = fill(Δt, T) + operator = GATES[:H] + systems(ζ) = QuantumSystem(ζ * GATES[:Z], [GATES[:X], GATES[:Y]]) + + prob = UnitarySamplingProblem( + systems, + Normal(0, 0.1), + n_samples, + operator, + T, + Δt; + verbose=false, + ipopt_options=Options(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + ) + + solve!(prob, max_iter=20) + + d_prob = UnitarySmoothPulseProblem( + systems(0), + operator, + T, + Δt; + verbose=false, + ipopt_options=Options(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + ) + solve!(prob, max_iter=20) + + # Check that the solution improves over the default + ζ_test = 0.02 + Ũ⃗_goal = operator_to_iso_vec(operator) + + Ũ⃗_end = unitary_rollout(prob.trajectory.a, timesteps, systems(ζ_test))[:, end] + fid = unitary_fidelity(Ũ⃗_end, Ũ⃗_goal) + + d_Ũ⃗_end = unitary_rollout(d_prob.trajectory.a, timesteps, systems(ζ_test))[:, end] + default_fid = unitary_fidelity(d_Ũ⃗_end, Ũ⃗_goal) + + @test fid > default_fid +end From 3c998803b8e84d82750ac3658deb17fb1eb3c599 Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 21 Jun 2024 13:12:35 -0500 Subject: [PATCH 13/35] =?UTF-8?q?bug=20fix:=20=CE=94t=20vector=20in=20kwar?= =?UTF-8?q?gs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/problem_templates/unitary_smooth_pulse_problem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 0db87bd1..b5ff4c8c 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -86,8 +86,8 @@ function UnitarySmoothPulseProblem( a_guess::Union{Matrix{Float64}, Nothing}=nothing, dda_bound::Float64=1.0, dda_bounds=fill(dda_bound, length(system.G_drives)), - Δt_min::Float64=0.5 * Δt, - Δt_max::Float64=1.5 * Δt, + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), + Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), drive_derivative_σ::Float64=0.01, Q::Float64=100.0, R=1e-2, From 6c0e7849226f42976a0f91b572660b68104a2413 Mon Sep 17 00:00:00 2001 From: andy Date: Sun, 23 Jun 2024 13:40:23 -0500 Subject: [PATCH 14/35] sampling docstring --- .../unitary_sampling_problem.jl | 48 ++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index 22c7ca6f..5520415a 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -1,7 +1,53 @@ @doc raw""" UnitarySamplingProblem - TODO: systems might need flexible bounds. +A `UnitarySamplingProblem` is a quantum control problem where the goal is to find a control pulse that generates a target unitary operator for a set of quantum systems. +The controls are shared among all systems, and the optimization seeks to find the control pulse that achieves the operator for each system. The idea is to enforce a +robust solution by including multiple systems reflecting the problem uncertainty. + +# Arguments +- `systems::AbstractVector{<:AbstractQuantumSystem}`: A vector of quantum systems. +- `operator::Union{EmbeddedOperator, AbstractMatrix{<:Number}}`: The target unitary operator. +- `T::Int`: The number of time steps. +- `Δt::Union{Float64, Vector{Float64}}`: The time step size. +- `system_labels::Vector{String}`: The labels for each system. +- `system_weights::Vector{Float64}`: The weights for each system. +- `free_time::Bool`: Whether to optimize the time steps. +- `init_trajectory::Union{NamedTrajectory, Nothing}`: The initial trajectory. +- `a_bound::Float64`: The bound for the control amplitudes. +- `a_bounds::Vector{Float64}`: The bounds for the control amplitudes. +- `a_guess::Union{Matrix{Float64}, Nothing}`: The initial guess for the control amplitudes. +- `dda_bound::Float64`: The bound for the control second derivatives. +- `dda_bounds::Vector{Float64}`: The bounds for the control second derivatives. +- `Δt_min::Float64`: The minimum time step size. +- `Δt_max::Float64`: The maximum time step size. +- `drive_derivative_σ::Float64`: The standard deviation for the drive derivative noise. +- `Q::Float64`: The fidelity weight. +- `R::Float64`: The regularization weight. +- `R_a::Union{Float64, Vector{Float64}}`: The regularization weight for the control amplitudes. +- `R_da::Union{Float64, Vector{Float64}}`: The regularization weight for the control first derivatives. +- `R_dda::Union{Float64, Vector{Float64}}`: The regularization weight for the control second derivatives. +- `leakage_suppression::Bool`: Whether to suppress leakage. +- `R_leakage::Float64`: The regularization weight for the leakage. +- `max_iter::Int`: The maximum number of iterations. +- `linear_solver::String`: The linear solver. +- `ipopt_options::Options`: The IPOPT options. +- `constraints::Vector{<:AbstractConstraint}`: The constraints. +- `timesteps_all_equal::Bool`: Whether to enforce equal time steps. +- `verbose::Bool`: Whether to print verbose output. +- `integrator::Symbol`: The integrator to use. +- `rollout_integrator`: The integrator for the rollout. +- `bound_unitary::Bool`: Whether to bound the unitary. +- `control_norm_constraint::Bool`: Whether to enforce a control norm constraint. +- `control_norm_constraint_components`: The components for the control norm constraint. +- `control_norm_R`: The regularization weight for the control norm constraint. +- `geodesic::Bool`: Whether to use the geodesic. +- `pade_order::Int`: The order of the Pade approximation. +- `autodiff::Bool`: Whether to use automatic differentiation. +- `jacobian_structure::Bool`: Whether to evaluate the Jacobian structure. +- `hessian_approximation::Bool`: Whether to approximate the Hessian. +- `blas_multithreading::Bool`: Whether to use BLAS multithreading. +- `kwargs...`: Additional keyword arguments. """ function UnitarySamplingProblem( From b01db3880c8acbdc4f8135d6c46b06191cc0ac65 Mon Sep 17 00:00:00 2001 From: andy Date: Sun, 23 Jun 2024 14:30:39 -0500 Subject: [PATCH 15/35] update tests --- .../unitary_sampling_problem.jl | 29 +++++++++++- src/trajectory_initialization.jl | 45 ++++++++++--------- 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index 5520415a..f1e779a7 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -124,7 +124,7 @@ function UnitarySamplingProblem( bound_unitary=bound_unitary, drive_derivative_σ=drive_derivative_σ, a_guess=a_guess, - system=systems[1], + system=systems, rollout_integrator=rollout_integrator, Ũ⃗_keys=Ũ⃗_keys ) @@ -284,4 +284,29 @@ end default_fid = unitary_fidelity(d_Ũ⃗_end, Ũ⃗_goal) @test fid > default_fid -end + + # Check initial guess initialization + a_guess = prob.trajectory.a + + g1_prob = UnitarySamplingProblem( + [systems(0), systems(0)], + operator, + T, + Δt; + verbose=false, + a_guess=a_guess, + ) + + @test g1_prob.trajectory.Ũ⃗1 ≈ g1_prob.trajectory.Ũ⃗2 + + g2_prob = UnitarySamplingProblem( + [systems(0), systems(0.1)], + operator, + T, + Δt; + verbose=false, + a_guess=a_guess, + ) + + @test ~(g2_prob.trajectory.Ũ⃗1 ≈ g2_prob.trajectory.Ũ⃗2) +end \ No newline at end of file diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index b08d9e16..02e9e3cb 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -10,6 +10,8 @@ using LinearAlgebra using Distributions using ..QuantumUtils +using ..QuantumSystems +using ..Rollouts using ..EmbeddedOperators """ @@ -164,10 +166,7 @@ function initialize_controls( return a, da, dda end -function initialize_controls( - a_guess::AbstractMatrix, -) - a = a_guess +function initialize_controls(a::AbstractMatrix, Δt::AbstractVecOrMat) da = derivative(a, Δt) dda = derivative(da, Δt) @@ -191,7 +190,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, - system=Union{AbstractQuantumSystem, Nothing}=nothing, + system::Union{AbstractQuantumSystem, AbstractVector{<:AbstractQuantumSystem}, Nothing}=nothing, rollout_integrator::Function=exp, Ũ⃗_keys::AbstractVector{<:Symbol}=[:Ũ⃗], ) @@ -208,22 +207,6 @@ function initialize_trajectory( Ũ⃗_goal = operator_to_iso_vec(U_goal) end - # Initial state and controls - if isnothing(a_guess) - Ũ⃗ = initialize_unitaries(U_init, U_goal, T, geodesic=geodesic) - a, da, dda = initialize_controls(n_drives, T, a_bounds, drive_derivative_σ) - else - @assert !isnothing(system) "system must be provided if a_guess is provided" - Ũ⃗ = unitary_rollout( - Ũ⃗_init, - a_guess, - Δt, - system; - integrator=rollout_integrator - ) - a, da, dda = initialize_controls(a_guess) - end - # Constraints Ũ⃗_inits = repeat([Ũ⃗_init], length(Ũ⃗_keys)) initial = (; @@ -247,8 +230,26 @@ function initialize_trajectory( merge!(bounds, (; (Ũ⃗_keys .=> Ũ⃗_bounds)...)) end + # Initial state and control values + if isnothing(a_guess) + Ũ⃗ = initialize_unitaries(U_init, U_goal, T, geodesic=geodesic) + Ũ⃗_values = repeat([Ũ⃗], length(Ũ⃗_keys)) + a, da, dda = initialize_controls(n_drives, T, a_bounds, drive_derivative_σ) + else + @assert !isnothing(system) "system must be provided if a_guess is provided" + if system isa AbstractVector + @assert length(system) == length(Ũ⃗_keys) "systems must have the same length as Ũ⃗_keys" + Ũ⃗_values = map(system) do sys + unitary_rollout(Ũ⃗_init, a_guess, Δt, sys; integrator=rollout_integrator) + end + else + unitary_rollout(Ũ⃗_init, a_guess, Δt, system; integrator=rollout_integrator) + Ũ⃗_values = repeat([Ũ⃗], length(Ũ⃗_keys)) + end + a, da, dda = initialize_controls(a_guess, Δt) + end + # Trajectory - Ũ⃗_values = repeat([Ũ⃗], length(Ũ⃗_keys)) keys = [Ũ⃗_keys..., :a, :da, :dda] values = [Ũ⃗_values..., a, da, dda] From 6dad0c0e812a11d36d6d385f664ed98749d2e7df Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 09:03:49 -0500 Subject: [PATCH 16/35] rename initialize_trajectory --- src/problem_solvers.jl | 4 ++-- src/problem_templates/unitary_sampling_problem.jl | 2 +- src/problem_templates/unitary_smooth_pulse_problem.jl | 2 +- src/problems.jl | 8 ++++---- src/trajectory_initialization.jl | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/problem_solvers.jl b/src/problem_solvers.jl index 860b58c0..651b8b1b 100644 --- a/src/problem_solvers.jl +++ b/src/problem_solvers.jl @@ -25,9 +25,9 @@ function solve!( set!(prob.optimizer, prob.options) if !isnothing(init_traj) - initialize_trajectory!(prob, init_traj) + set_trajectory!(prob, init_traj) else - initialize_trajectory!(prob) + set_trajectory!(prob) end MOI.optimize!(prob.optimizer) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index f1e779a7..4fb8f2fe 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -111,7 +111,7 @@ function UnitarySamplingProblem( else n_drives = length(systems[1].G_drives) # TODO: Initial system? - traj = initialize_trajectory( + traj = initialize_unitary_trajectory( operator, T, Δt, diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 85707dbb..cd425a82 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -130,7 +130,7 @@ function UnitarySmoothPulseProblem( traj = init_trajectory else n_drives = length(system.G_drives) - traj = initialize_trajectory( + traj = initialize_unitary_trajectory( operator, T, Δt, diff --git a/src/problems.jl b/src/problems.jl index 8e4a867d..82bc728f 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -4,7 +4,7 @@ export AbstractProblem export FixedTimeProblem export QuantumControlProblem -export initialize_trajectory! +export set_trajectory! export update_trajectory! export get_traj_data export get_datavec @@ -265,7 +265,7 @@ function initialize_optimizer!( return variables end -function initialize_trajectory!( +function set_trajectory!( prob::QuantumControlProblem, traj::NamedTrajectory ) @@ -277,8 +277,8 @@ function initialize_trajectory!( ) end -initialize_trajectory!(prob::QuantumControlProblem) = - initialize_trajectory!(prob, prob.trajectory) +set_trajectory!(prob::QuantumControlProblem) = + set_trajectory!(prob, prob.trajectory) function get_datavec(prob::QuantumControlProblem) Z⃗ = MOI.get( diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 02e9e3cb..f26659f9 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -3,7 +3,7 @@ module TrajectoryInitialization export unitary_linear_interpolation export unitary_geodesic export linear_interpolation -export initialize_trajectory +export initialize_unitary_trajectory using NamedTrajectories using LinearAlgebra @@ -176,7 +176,7 @@ function initialize_controls(a::AbstractMatrix, Δt::AbstractVecOrMat) return a, da, dda end -function initialize_trajectory( +function initialize_unitary_trajectory( U_goal::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, T::Int, Δt::Real, From da48d9470c3fade5c0bcfe97322abf9c8d092409 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 09:32:34 -0500 Subject: [PATCH 17/35] add state initialization --- .../quantum_state_smooth_pulse_problem.jl | 174 +++++------------- src/trajectory_initialization.jl | 78 ++++++++ 2 files changed, 123 insertions(+), 129 deletions(-) diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index b93ac72e..7d2be9b6 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -35,6 +35,7 @@ function QuantumStateSmoothPulseProblem( a_bound::Float64=Inf, a_bounds::Vector{Float64}=fill(a_bound, length(system.G_drives)), a_guess::Union{Matrix{Float64}, Nothing}=nothing, + rollout_integrator=exp, dda_bound::Float64=Inf, dda_bounds::Vector{Float64}=fill(dda_bound, length(system.G_drives)), Δt_min::Float64=0.5 * Δt, @@ -75,108 +76,28 @@ function QuantumStateSmoothPulseProblem( n_drives = length(system.G_drives) + # Trajectory if !isnothing(init_trajectory) traj = init_trajectory else - if free_time - Δt = fill(Δt, T) - end - - if isnothing(a_guess) - ψ̃s = NamedTuple([ - Symbol("ψ̃$i") => linear_interpolation(ψ̃_init, ψ̃_goal, T) - for (i, (ψ̃_init, ψ̃_goal)) in enumerate(zip(ψ̃_inits, ψ̃_goals)) - ]) - a_dists = [Uniform(-a_bounds[i], a_bounds[i]) for i = 1:n_drives] - a = hcat([ - zeros(n_drives), - vcat([rand(a_dists[i], 1, T - 2) for i = 1:n_drives]...), - zeros(n_drives) - ]...) - da = randn(n_drives, T) * drive_derivative_σ - dda = randn(n_drives, T) * drive_derivative_σ - else - ψ̃s = NamedTuple([ - Symbol("ψ̃$i") => rollout(ψ̃_init, a_guess, Δt, system) - for (i, ψ̃_init) in enumerate(ψ̃_inits) - ]) - a = a_guess - da = derivative(a, Δt) - dda = derivative(da, Δt) - end - - ψ̃_initial = NamedTuple([ - Symbol("ψ̃$i") => ψ̃_init - for (i, ψ̃_init) in enumerate(ψ̃_inits) - ]) - - control_initial = ( - a = zeros(n_drives), + traj = initialize_state_trajectory( + ψ̃_goals, + ψ̃_inits, + T, + Δt, + n_drives, + a_bounds, + dda_bounds; + free_time=free_time, + Δt_bounds=(Δt_min, Δt_max), + drive_derivative_σ=drive_derivative_σ, + a_guess=a_guess, + system=system, + rollout_integrator=rollout_integrator, ) - - initial = merge(ψ̃_initial, control_initial) - - final = ( - a = zeros(n_drives), - ) - - goal = NamedTuple([ - Symbol("ψ̃$i") => ψ̃_goal - for (i, ψ̃_goal) in enumerate(ψ̃_goals) - ]) - - if free_time - - control_components = ( - a = a, - da = da, - dda = dda, - Δt = Δt, - ) - - components = merge(ψ̃s, control_components) - - bounds = ( - a = a_bounds, - dda = dda_bounds, - Δt = (Δt_min, Δt_max), - ) - - traj = NamedTrajectory( - components; - controls=(:dda, :Δt), - timestep=:Δt, - bounds=bounds, - initial=initial, - final=final, - goal=goal - ) - else - control_components = ( - a = a, - da = da, - dda = dda, - ) - - components = merge(ψ̃s, control_components) - - bounds = ( - a = a_bounds, - dda = dda_bounds, - ) - - traj = NamedTrajectory( - components; - controls=(:dda,), - timestep=Δt, - bounds=bounds, - initial=initial, - final=final, - goal=goal - ) - end end + # Objective J = QuadraticRegularizer(:a, traj, R_a) J += QuadraticRegularizer(:da, traj, R_da) J += QuadraticRegularizer(:dda, traj, R_dda) @@ -185,43 +106,22 @@ function QuantumStateSmoothPulseProblem( J += QuantumStateObjective(Symbol("ψ̃$i"), traj, Q) end - L1_slack_constraints = [] - + # Constraints for name in L1_regularized_names if name in keys(L1_regularized_indices) - J_L1, slack_con = L1Regularizer(name, traj; R_value=R_L1, indices=L1_regularized_indices[name]) + J += L1Regularizer!( + constraints, name, traj, + R_value=R_L1, + indices=L1_regularized_indices[name], + eval_hessian=!hessian_approximation + ) else - J_L1, slack_con = L1Regularizer(name, traj; R_value=R_L1) + J += L1Regularizer!( + constraints, name, traj; + R_value=R_L1, + eval_hessian=!hessian_approximation + ) end - J += J_L1 - push!(L1_slack_constraints, slack_con) - end - - append!(constraints, L1_slack_constraints) - - if free_time - - ψ̃_integrators = [ - QuantumStatePadeIntegrator(system, Symbol("ψ̃$i"), :a) - for i = 1:length(ψ_inits) - ] - - integrators = [ - ψ̃_integrators..., - DerivativeIntegrator(:a, :da, traj), - DerivativeIntegrator(:da, :dda, traj) - ] - else - ψ̃_integrators = [ - QuantumStatePadeIntegrator(system, Symbol("ψ̃$i"), :a) - for i = 1:length(ψ_inits) - ] - - integrators = [ - ψ̃_integrators..., - DerivativeIntegrator(:a, :da, traj), - DerivativeIntegrator(:da, :dda, traj) - ] end if free_time @@ -230,6 +130,18 @@ function QuantumStateSmoothPulseProblem( end end + # Integrators + ψ̃_integrators = [ + QuantumStatePadeIntegrator(system, Symbol("ψ̃$i"), :a) + for i = 1:length(ψ_inits) + ] + + integrators = [ + ψ̃_integrators..., + DerivativeIntegrator(:a, :da, traj), + DerivativeIntegrator(:da, :dda, traj) + ] + return QuantumControlProblem( system, traj, @@ -253,3 +165,7 @@ function QuantumStateSmoothPulseProblem( system = QuantumSystem(H_drift, H_drives) return QuantumStateSmoothPulseProblem(system, args...; kwargs...) end + +# *************************************************************************** # + +# TODO: Tests diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index f26659f9..f38ebfd5 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -275,4 +275,82 @@ function initialize_unitary_trajectory( ) end +function initialize_state_trajectory( + ψ̃_goals::AbstractVector{<:AbstractVector{<:Real}}, + ψ̃_inits::AbstractVector{<:AbstractVector{<:Real}}, + T::Int, + Δt::Real, + n_drives::Int, + a_bounds::VectorBound, + dda_bounds::VectorBound; + free_time=false, + Δt_bounds::ScalarBound=(0.5 * Δt, 1.5 * Δt), + drive_derivative_σ::Float64=0.1, + a_guess::Union{AbstractMatrix{<:Float64}, Nothing}=nothing, + system::Union{AbstractQuantumSystem, AbstractVector{<:AbstractQuantumSystem}, Nothing}=nothing, + rollout_integrator::Function=exp, + ψ̃_keys::AbstractVector{<:Symbol}=[Symbol("ψ̃$i") for i = 1:length(ψ̃_goals)] +) + @assert length(ψ̃_inits) == length(ψ̃_goals) "ψ̃_inits and ψ̃_goals must have the same length" + @assert length(ψ̃_keys) == length(ψ̃_goals) "ψ̃_keys and ψ̃_goals must have the same length" + + if free_time + if Δt isa Float64 + Δt = fill(Δt, 1, T) + end + end + + # Constraints + state_initial = (; (ψ̃_keys .=> ψ̃_inits)...) + control_initial = (a = zeros(n_drives),) + initial = merge(state_initial, control_initial) + + final = (a = zeros(n_drives),) + + goal = (; (ψ̃_keys .=> ψ̃_goals)...) + + # Bounds + bounds = (a = a_bounds, dda = dda_bounds,) + + # Initial state and control values + if isnothing(a_guess) + ψ̃s = NamedTuple([ + k => linear_interpolation(ψ̃_init, ψ̃_goal, T) + for (k, ψ̃_init, ψ̃_goal) in zip(ψ̃_keys, ψ̃_inits, ψ̃_goals) + ]) + a, da, dda = initialize_controls(n_drives, T, a_bounds, drive_derivative_σ) + else + ψ̃s = NamedTuple([ + k => rollout(ψ̃_init, a_guess, Δt, system, integrator=rollout_integrator) + for (i, ψ̃_init) in zip(ψ̃_keys, ψ̃_inits) + ]) + a, da, dda = initialize_controls(a_guess, Δt) + end + + # Trajectory + keys = [ψ̃_keys..., :a, :da, :dda] + values = [ψ̃s..., a, da, dda] + + if free_time + push!(keys, :Δt) + push!(values, Δt) + controls = (:dda, :Δt) + timestep = :Δt + bounds = merge(bounds, (Δt = Δt_bounds,)) + else + controls = (:dda,) + timestep = Δt + end + + return NamedTrajectory( + (; (keys .=> values)...); + controls=controls, + timestep=timestep, + bounds=bounds, + initial=initial, + final=final, + goal=goal + ) +end + end From 39bc244ae172c9e7c5b5a23b113e830f779bae12 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 12:21:12 -0500 Subject: [PATCH 18/35] G fn in rollouts and add state integrator tests --- src/integrators/_integrator_utils.jl | 4 +- src/integrators/pade_integrators.jl | 53 ++++++++++--------- .../quantum_state_smooth_pulse_problem.jl | 41 +++++++++++--- src/rollouts.jl | 51 ++++++++++++------ src/trajectory_initialization.jl | 1 + 5 files changed, 99 insertions(+), 51 deletions(-) diff --git a/src/integrators/_integrator_utils.jl b/src/integrators/_integrator_utils.jl index 955679e7..f91d99c0 100644 --- a/src/integrators/_integrator_utils.jl +++ b/src/integrators/_integrator_utils.jl @@ -1,5 +1,5 @@ -# G(a) helper function -function G( +# G_bilinear(a) helper function +function G_bilinear( a::AbstractVector, G_drift::AbstractMatrix, G_drives::AbstractVector{<:AbstractMatrix} diff --git a/src/integrators/pade_integrators.jl b/src/integrators/pade_integrators.jl index 25720765..aafd2ba2 100644 --- a/src/integrators/pade_integrators.jl +++ b/src/integrators/pade_integrators.jl @@ -53,7 +53,7 @@ struct UnitaryPadeIntegrator <: QuantumPadeIntegrator dim::Int order::Int autodiff::Bool - G::Union{Function, Nothing} + G::Function """ UnitaryPadeIntegrator( @@ -111,6 +111,7 @@ struct UnitaryPadeIntegrator <: QuantumPadeIntegrator G_drift = sys.G_drift G_drives = sys.G_drives + G = isnothing(G) ? G_bilinear : G drive_anticomms, drift_anticomms = order == 4 ? build_anticomms(G_drift, G_drives, n_drives) : (nothing, nothing) @@ -193,13 +194,14 @@ struct QuantumStatePadeIntegrator <: QuantumPadeIntegrator @assert order ∈ [4, 6, 8, 10] "order must be in [4, 6, 8, 10]" @assert !isnothing(state_symb) "state_symb must be specified" @assert !isnothing(drive_symb) "drive_symb must be specified" - n_drives = length(sys.H_drives_real) - N = size(sys.H_drift_real, 1) + n_drives = length(sys.H_drives) + N = size(sys.H_drift, 1) dim = 2N I_2N = sparse(I(2N)) G_drift = sys.G_drift G_drives = sys.G_drives + G = isnothing(G) ? G_bilinear : G drive_anticomms, drift_anticomms = order == 4 ? build_anticomms(G_drift, G_drives, n_drives) : (nothing, nothing) @@ -234,7 +236,7 @@ function nth_order_pade( ) Ũₜ₊₁ = iso_vec_to_iso_operator(Ũ⃗ₜ₊₁) Ũₜ = iso_vec_to_iso_operator(Ũ⃗ₜ) - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) n = P.order ÷ 2 Gₜ_powers = compute_powers(Gₜ, n) B = P.I_2N + sum([ @@ -278,7 +280,7 @@ function nth_order_pade( aₜ::AbstractVector, Δt::Real ) - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) n = P.order ÷ 2 Gₜ_powers = compute_powers(Gₜ, n) B = P.I_2N + sum([ @@ -341,7 +343,7 @@ function ∂aₜ( for j = 1:n_drives Gʲ = P.G_drives[j] Gʲ_anticomm_Gₜ = - G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) for i = 0:P.N-1 ψ̃ⁱₜ₊₁ = @view Ũ⃗ₜ₊₁[i * isodim .+ (1:isodim)] ψ̃ⁱₜ = @view Ũ⃗ₜ[i * isodim .+ (1:isodim)] @@ -381,7 +383,7 @@ function ∂aₜ( for j = 1:n_drives Gʲ = P.G_drives[j] Gʲ_anticomm_Gₜ = - G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂aP[:, j] = -Δtₜ / 2 * Gʲ * (ψ̃ₜ₊₁ + ψ̃ₜ) + Δtₜ^2 / 12 * Gʲ_anticomm_Gₜ * (ψ̃ₜ₊₁ - ψ̃ₜ) @@ -402,8 +404,7 @@ function ∂Δtₜ( aₜ::AbstractVector, Δtₜ::Real ) - - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) Ũₜ₊₁ = iso_vec_to_iso_operator(Ũ⃗ₜ₊₁) Ũₜ = iso_vec_to_iso_operator(Ũ⃗ₜ) if P.order == 4 @@ -435,7 +436,7 @@ function ∂Δtₜ( Δtₜ::Real ) - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) if P.order==4 ∂ΔtₜP = -1/2 * Gₜ * (ψ̃ₜ₊₁ + ψ̃ₜ) + 1/6 * Δtₜ * Gₜ^2 * (ψ̃ₜ₊₁ - ψ̃ₜ) else @@ -481,7 +482,7 @@ end ∂Ũ⃗ₜP = spzeros(T, P.dim, P.dim) ∂Ũ⃗ₜ₊₁P = spzeros(T, P.dim, P.dim) - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) n = P.order ÷ 2 # can memoize this chunk of code, prly memoize G powers @@ -531,7 +532,7 @@ end ∂ΔtₜP = ∂Δtₜ(P, ψ̃ₜ₊₁, ψ̃ₜ, aₜ, Δtₜ) end - Gₜ = isnothing(P.G) ? G(aₜ, P.G_drift, P.G_drives) : P.G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) n = P.order ÷ 2 Gₜ_powers = compute_powers(Gₜ, n) B = P.I_2N + sum([(-1)^k * PADE_COEFFICIENTS[P.order][k] * Δtₜ^k * Gₜ_powers[k] for k = 1:n]) @@ -562,14 +563,14 @@ function μ∂aₜ∂Ũ⃗ₜ( n_drives = P.n_drives - if P.autodiff || !isnothing(P.G) + if P.autodiff elseif P.order == 4 μ∂aₜ∂Ũ⃗ₜP = Array{T}(undef, P.dim, n_drives) for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂aₜ∂Ũ⃗ₜ_block_i = -(Δtₜ / 2 * Gʲ + Δtₜ^2 / 12 * Ĝʲ) # sparse is necessary since blockdiag doesn't accept dense matrices ∂aₜ∂Ũ⃗ₜ = blockdiag(fill(sparse(∂aₜ∂Ũ⃗ₜ_block_i), P.N)...) @@ -593,7 +594,7 @@ function μ∂Ũ⃗ₜ₊₁∂aₜ( for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂Ũ⃗ₜ₊₁∂aₜ_block_i = -Δtₜ / 2 * Gʲ + Δtₜ^2 / 12 * Ĝʲ # sparse is necessary since blockdiag doesn't accept dense matrices ∂Ũ⃗ₜ₊₁∂aₜ = blockdiag(fill(sparse(∂Ũ⃗ₜ₊₁∂aₜ_block_i), P.N)...) @@ -615,7 +616,7 @@ function μ∂aₜ∂ψ̃ₜ( for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂aₜ∂ψ̃ₜP = -(Δtₜ / 2 * Gʲ + Δtₜ^2 / 12 * Ĝʲ) μ∂aₜ∂ψ̃ₜP[:, j] = ∂aₜ∂ψ̃ₜP' * μₜ end @@ -635,7 +636,7 @@ function μ∂ψ̃ₜ₊₁∂aₜ( for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂ψ̃ₜ₊₁∂aₜP = -Δtₜ / 2 * Gʲ + Δtₜ^2 / 12 * Ĝʲ μ∂ψ̃ₜ₊₁∂aₜP[j, :] = μₜ' * ∂ψ̃ₜ₊₁∂aₜP end @@ -677,7 +678,7 @@ function μ∂²aₜ( μₜ::AbstractVector{T}, ) where T <: Real - n_drives = length(drive_indices) + n_drives = P.n_drives μ∂²aₜP = Array{T}(undef, n_drives, n_drives) if P.order==4 @@ -707,7 +708,7 @@ function μ∂Δtₜ∂aₜ( if P.order == 4 for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) B = blockdiag(fill(sparse(-1/2 * Gʲ + 1/6 * Δtₜ * Ĝʲ), P.N)...) F = blockdiag(fill(sparse(1/2 * Gʲ + 1/6 * Δtₜ * Ĝʲ), P.N)...) ∂Δtₜ∂aₜ_j = B*Ũ⃗ₜ₊₁ - F*Ũ⃗ₜ @@ -732,7 +733,7 @@ function μ∂Δtₜ∂aₜ( if P.order == 4 for j = 1:n_drives Gʲ = P.G_drives[j] - Ĝʲ = G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) + Ĝʲ = P.G(aₜ, P.G_drift_anticomms[j], P.G_drive_anticomms[:, j]) ∂Δt∂aʲP = -1 / 2 * Gʲ * (ψ̃ₜ₊₁ + ψ̃ₜ) + 1 / 6 * Δtₜ * Ĝʲ * (ψ̃ₜ₊₁ - ψ̃ₜ) @@ -748,7 +749,7 @@ function μ∂Δtₜ∂Ũ⃗ₜ( Δtₜ::Real, μₜ::AbstractVector ) - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) minus_F = -(1/2 * Gₜ + 1/6 * Δtₜ * Gₜ^2) big_minus_F = blockdiag(fill(sparse(minus_F), P.N)...) return big_minus_F' * μₜ @@ -760,7 +761,7 @@ function μ∂Ũ⃗ₜ₊₁∂Δtₜ( Δtₜ::Real, μₜ::AbstractVector ) - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) B = -1/2 * Gₜ + 1/6 * Δtₜ * Gₜ^2 big_B = blockdiag(fill(sparse(B), P.N)...) return μₜ' * big_B @@ -773,7 +774,7 @@ function μ∂Δtₜ∂ψ̃ₜ( μₜ::AbstractVector ) # memoize the calc here - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) minus_F = -(1/2 * Gₜ + 1/6 * Δtₜ * Gₜ^2) return minus_F' * μₜ end @@ -784,7 +785,7 @@ function μ∂ψ̃ₜ₊₁∂Δtₜ( Δtₜ::Real, μₜ::AbstractVector ) - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) B = -1/2 * Gₜ + 1/6 * Δtₜ * Gₜ^2 return μₜ' * B end @@ -796,7 +797,7 @@ function μ∂²Δtₜ( aₜ::AbstractVector, μₜ::AbstractVector ) - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) ∂²Δtₜ_gen_block = 1/6 * Gₜ^2 ∂²Δtₜ_gen = blockdiag(fill(sparse(∂²Δtₜ_gen_block), P.N)...) ∂²Δtₜ = ∂²Δtₜ_gen * (Ũ⃗ₜ₊₁ - Ũ⃗ₜ) @@ -810,7 +811,7 @@ function μ∂²Δtₜ( aₜ::AbstractVector, μₜ::AbstractVector ) - Gₜ = G(aₜ, P.G_drift, P.G_drives) + Gₜ = P.G(aₜ, P.G_drift, P.G_drives) ∂²Δtₜ = 1/6 * Gₜ^2 * (ψ̃ₜ₊₁ - ψ̃ₜ) return μₜ' * ∂²Δtₜ end diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 7d2be9b6..739db6d3 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -32,11 +32,11 @@ function QuantumStateSmoothPulseProblem( Δt::Float64; free_time=true, init_trajectory::Union{NamedTrajectory, Nothing}=nothing, - a_bound::Float64=Inf, + a_bound::Float64=1.0, a_bounds::Vector{Float64}=fill(a_bound, length(system.G_drives)), a_guess::Union{Matrix{Float64}, Nothing}=nothing, rollout_integrator=exp, - dda_bound::Float64=Inf, + dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, length(system.G_drives)), Δt_min::Float64=0.5 * Δt, Δt_max::Float64=1.5 * Δt, @@ -68,11 +68,8 @@ function QuantumStateSmoothPulseProblem( ψ_goals = ψ_goal end - ψ_inits = Vector{ComplexF64}.(ψ_init) - ψ̃_inits = ket_to_iso.(ψ_init) - - ψ_goals = Vector{ComplexF64}.(ψ_goal) - ψ̃_goals = ket_to_iso.(ψ_goal) + ψ̃_inits = ket_to_iso.(Vector{ComplexF64}.(ψ_inits)) + ψ̃_goals = ket_to_iso.(Vector{ComplexF64}.(ψ_goals)) n_drives = length(system.G_drives) @@ -168,4 +165,32 @@ end # *************************************************************************** # -# TODO: Tests +@testitem "Test quantum state smooth pulse" begin + # System + T = 50 + Δt = 0.2 + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + ψ_init = [1.0, 0.0] + ψ_target = [0.0, 1.0] + + prob = QuantumStateSmoothPulseProblem( + sys, ψ_init, ψ_target, T, Δt; + ipopt_options=Options(print_level=1), verbose=false + ) + initial = fidelity(prob) + solve!(prob, max_iter=20) + final = fidelity(prob) + @test final > initial + + # Multiple initial and target states + ψ_inits = [[1.0, 0.0], [0.0, 1.0]] + ψ_targets = [[0.0, 1.0], [1.0, 0.0]] + prob = QuantumStateSmoothPulseProblem( + sys, ψ_inits, ψ_targets, T, Δt; + ipopt_options=Options(print_level=1), verbose=false + ) + initial = fidelity(prob) + solve!(prob, max_iter=20) + final = fidelity(prob) + @test all(final .> initial) +end diff --git a/src/rollouts.jl b/src/rollouts.jl index c6fd2e4a..16df25c9 100644 --- a/src/rollouts.jl +++ b/src/rollouts.jl @@ -21,7 +21,8 @@ function rollout( controls::AbstractMatrix{Float64}, Δt::Union{AbstractVector{Float64}, AbstractMatrix{Float64}, Float64}, system::AbstractQuantumSystem; - integrator=Integrators.fourth_order_pade + integrator=exp, + G=Integrators.G_bilinear ) if Δt isa AbstractMatrix @assert size(Δt, 1) == 1 @@ -41,11 +42,7 @@ function rollout( for t = 2:T aₜ₋₁ = controls[:, t - 1] - Gₜ = Integrators.G( - aₜ₋₁, - G_drift, - G_drives - ) + Gₜ = G(aₜ₋₁, G_drift, G_drives) Ψ̃[:, t] .= integrator(Gₜ * Δt[t - 1]) * Ψ̃[:, t - 1] end @@ -81,16 +78,46 @@ function QuantumUtils.fidelity( system::AbstractQuantumSystem; kwargs... ) - return fidelity(ket_to_iso(ψ₁), ket_to_iso(ψ_goal), args...; kwargs...) + return fidelity(ket_to_iso(ψ₁), ket_to_iso(ψ_goal), controls, Δt, system; kwargs...) end +function QuantumUtils.fidelity( + trajectory::NamedTrajectory, + system::AbstractQuantumSystem; + state_symb::Symbol=:ψ̃, + control_symb=:a, + kwargs... +) + fids = [] + for name in trajectory.names + if startswith(string(name), string(state_symb)) + init = trajectory.initial[name] + goal = trajectory.goal[name] + push!( + fids, + fidelity(init, goal, trajectory[control_symb], get_timesteps(trajectory), system; kwargs...) + ) + end + end + return fids +end + +function QuantumUtils.fidelity( + prob::QuantumControlProblem; + kwargs... +) + return fidelity(prob.trajectory, prob.system; kwargs...) +end + +# ============================================================================= # function unitary_rollout( Ũ⃗₁::AbstractVector{<:Real}, controls::AbstractMatrix{Float64}, Δt::Union{AbstractVector{Float64}, AbstractMatrix{Float64}, Float64}, system::AbstractQuantumSystem; - integrator=exp + integrator=exp, + G=Integrators.G_bilinear ) if Δt isa AbstractMatrix @assert size(Δt, 1) == 1 @@ -110,11 +137,7 @@ function unitary_rollout( for t = 2:T aₜ₋₁ = controls[:, t - 1] - Gₜ = Integrators.G( - aₜ₋₁, - G_drift, - G_drives - ) + Gₜ = G(aₜ₋₁, G_drift, G_drives) Ũ⃗ₜ₋₁ = Ũ⃗[:, t - 1] Ũₜ₋₁ = iso_vec_to_iso_operator(Ũ⃗ₜ₋₁) Ũₜ = integrator(Gₜ * Δt[t - 1]) * Ũₜ₋₁ @@ -125,8 +148,6 @@ function unitary_rollout( return Ũ⃗ end - - function unitary_rollout( controls::AbstractMatrix{Float64}, Δt::Union{AbstractVector{Float64}, AbstractMatrix{Float64}, Float64}, diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index f38ebfd5..1e863865 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -4,6 +4,7 @@ export unitary_linear_interpolation export unitary_geodesic export linear_interpolation export initialize_unitary_trajectory +export initialize_state_trajectory using NamedTrajectories using LinearAlgebra From 425fc1c3bc040134c6a2b4fb872fc94f49fa5c17 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 12:41:03 -0500 Subject: [PATCH 19/35] add geodesic method for goal, sample --- src/trajectory_initialization.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 1e863865..926b8702 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -85,12 +85,23 @@ end function unitary_geodesic( U_init::AbstractMatrix{<:Number}, U_goal::AbstractMatrix{<:Number}, - samples::Number; + samples::Int; kwargs... ) return unitary_geodesic(U_init, U_goal, range(0, 1, samples); kwargs...) end +function unitary_geodesic( + U_goal::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, + 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 + function unitary_geodesic( U₀::AbstractMatrix{<:Number}, U₁::AbstractMatrix{<:Number}, From a8a63f2f9b13c51d62d2fca08e4e77835b6b9290 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 13:18:36 -0500 Subject: [PATCH 20/35] exp integrators G_bilinear --- src/integrators/exponential_integrators.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/integrators/exponential_integrators.jl b/src/integrators/exponential_integrators.jl index 3390e554..9c357345 100644 --- a/src/integrators/exponential_integrators.jl +++ b/src/integrators/exponential_integrators.jl @@ -15,6 +15,7 @@ struct UnitaryExponentialIntegrator <: QuantumExponentialIntegrator n_drives::Int ketdim::Int dim::Int + G::Function function UnitaryExponentialIntegrator( sys::AbstractQuantumSystem, @@ -32,7 +33,8 @@ struct UnitaryExponentialIntegrator <: QuantumExponentialIntegrator drive_names, n_drives, ketdim, - dim + dim, + G_bilinear ) end end @@ -60,7 +62,7 @@ controls(integrator::UnitaryExponentialIntegrator) = integrator.drive_names aₜ = zₜ[traj.components[ℰ.drive_names]] end - Gₜ = G(aₜ, ℰ.G_drift, ℰ.G_drives) + Gₜ = ℰ.G(aₜ, ℰ.G_drift, ℰ.G_drives) return Ũ⃗ₜ₊₁ - expv(Δtₜ, I(ℰ.ketdim) ⊗ Gₜ, Ũ⃗ₜ) end @@ -98,7 +100,7 @@ end aₜ = zₜ[inds] - Gₜ = G(aₜ, ℰ.G_drift, ℰ.G_drives) + Gₜ = ℰ.G(aₜ, ℰ.G_drift, ℰ.G_drives) Id = I(ℰ.ketdim) @@ -108,7 +110,7 @@ end ∂Ũ⃗ₜℰ = -expĜₜ ∂aₜℰ = -ForwardDiff.jacobian( - a -> expv(Δtₜ, Id ⊗ G(a, ℰ.G_drift, ℰ.G_drives), Ũ⃗ₜ), + a -> expv(Δtₜ, Id ⊗ ℰ.G(a, ℰ.G_drift, ℰ.G_drives), Ũ⃗ₜ), aₜ ) @@ -128,6 +130,7 @@ struct QuantumStateExponentialIntegrator <: QuantumExponentialIntegrator n_drives::Int ketdim::Int dim::Int + G::Function function QuantumStateExponentialIntegrator( sys::AbstractQuantumSystem, @@ -145,7 +148,8 @@ struct QuantumStateExponentialIntegrator <: QuantumExponentialIntegrator drive_names, n_drives, ketdim, - dim + dim, + G_bilinear ) end end @@ -173,7 +177,7 @@ controls(integrator::QuantumStateExponentialIntegrator) = integrator.drive_names aₜ = zₜ[traj.components[ℰ.drive_names]] end - Gₜ = G(aₜ, ℰ.G_drift, ℰ.G_drives) + Gₜ = ℰ.G(aₜ, ℰ.G_drift, ℰ.G_drives) return ψ̃ₜ₊₁ - expv(Δtₜ, Gₜ, ψ̃ₜ) end @@ -203,7 +207,7 @@ end aₜ = zₜ[inds] - Gₜ = G(aₜ, ℰ.G_drift, ℰ.G_drives) + Gₜ = ℰ.G(aₜ, ℰ.G_drift, ℰ.G_drives) expGₜ = hermitian_exp(Δtₜ * Gₜ) @@ -211,7 +215,7 @@ end ∂ψ̃ₜℰ = -expGₜ ∂aₜℰ = -ForwardDiff.jacobian( - a -> expv(Δtₜ, G(a, ℰ.G_drift, ℰ.G_drives), ψ̃ₜ), + a -> expv(Δtₜ, ℰ.G(a, ℰ.G_drift, ℰ.G_drives), ψ̃ₜ), aₜ ) From 69110c230bca790008c9c147c21b145afe78f97e Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 13:25:56 -0500 Subject: [PATCH 21/35] break out exp integ tests into testitems --- test/integrators_test.jl | 4 ++-- test/runtests.jl | 2 -- test/{ => scripts}/integrator_test_1qubit.jl | 0 test/{ => scripts}/integrator_test_script.jl | 0 4 files changed, 2 insertions(+), 4 deletions(-) rename test/{ => scripts}/integrator_test_1qubit.jl (100%) rename test/{ => scripts}/integrator_test_script.jl (100%) diff --git a/test/integrators_test.jl b/test/integrators_test.jl index 693cd05b..df1a23d3 100644 --- a/test/integrators_test.jl +++ b/test/integrators_test.jl @@ -1,4 +1,4 @@ -@testset "testing UnitaryExponentialIntegrator" begin +@testitem "testing UnitaryExponentialIntegrator" begin using NamedTrajectories using ForwardDiff @@ -45,7 +45,7 @@ @test ∂Δtₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.Δt] end -@testset "testing QuantumStateExponentialIntegrator" begin +@testitem "testing QuantumStateExponentialIntegrator" begin using NamedTrajectories using ForwardDiff diff --git a/test/runtests.jl b/test/runtests.jl index 4206b8e6..804eaae3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,8 +17,6 @@ include("test_utils.jl") @testset "QuantumCollocation.jl" begin # include("objectives_test.jl") # include("dynamics_test.jl") - include("integrators_test.jl") - include("quantum_utils_test.jl") include("quantum_system_templates_test.jl") end diff --git a/test/integrator_test_1qubit.jl b/test/scripts/integrator_test_1qubit.jl similarity index 100% rename from test/integrator_test_1qubit.jl rename to test/scripts/integrator_test_1qubit.jl diff --git a/test/integrator_test_script.jl b/test/scripts/integrator_test_script.jl similarity index 100% rename from test/integrator_test_script.jl rename to test/scripts/integrator_test_script.jl From 42fe0caf11fba565e041e20847c6fa5c12073b9c Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 14:46:15 -0500 Subject: [PATCH 22/35] embedded op linear interp --- src/trajectory_initialization.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 926b8702..57e866e2 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -1,8 +1,8 @@ module TrajectoryInitialization -export unitary_linear_interpolation export unitary_geodesic export linear_interpolation +export unitary_linear_interpolation export initialize_unitary_trajectory export initialize_state_trajectory @@ -36,6 +36,14 @@ function unitary_linear_interpolation( return Ũ⃗ end +function unitary_linear_interpolation( + U_init::AbstractMatrix{<:Number}, + U_goal::EmbeddedOperator, + samples::Int +) + return unitary_linear_interpolation(U_init, U_goal.operator, samples) +end + """ unitary_geodesic( operator::EmbeddedOperator, From a9ed0c6b66155547684eae589c328b748924f7c5 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 15:01:29 -0500 Subject: [PATCH 23/35] remove println --- .../unitary_direct_sum_problem.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/problem_templates/unitary_direct_sum_problem.jl b/src/problem_templates/unitary_direct_sum_problem.jl index 7c632423..c367f7f8 100644 --- a/src/problem_templates/unitary_direct_sum_problem.jl +++ b/src/problem_templates/unitary_direct_sum_problem.jl @@ -122,8 +122,6 @@ function UnitaryDirectSumProblem( # Build the direct sum system - println("...trajectories...") - # merge suffix trajectories traj = direct_sum([add_suffix(p.trajectory, ℓ) for (p, ℓ) ∈ zip(probs, prob_labels)]) @@ -154,8 +152,6 @@ function UnitaryDirectSumProblem( build_trajectory_constraints = true constraints = AbstractConstraint[] - println("...constraints...") - # Add goal constraints for each problem for (p, ℓ) ∈ zip(probs, prob_labels) goal_symb, = keys(p.trajectory.goal) @@ -280,11 +276,11 @@ end T = 50 Δt = 0.2 ops = Options(print_level=1) - prob1 = UnitarySmoothPulseProblem(sys, U_goal1, T, Δt, free_time=false, ipopt_options=ops) - prob2 = UnitarySmoothPulseProblem(sys, U_goal2, T, Δt, free_time=false, ipopt_options=ops) + prob1 = UnitarySmoothPulseProblem(sys, U_goal1, T, Δt, free_time=false, verbose=false, ipopt_options=ops) + prob2 = UnitarySmoothPulseProblem(sys, U_goal2, T, Δt, free_time=false, verbose=false, ipopt_options=ops) # Test default - direct_sum_prob1 = UnitaryDirectSumProblem([prob1, prob2], 0.99, ipopt_options=ops) + direct_sum_prob1 = UnitaryDirectSumProblem([prob1, prob2], 0.99, verbose=false, ipopt_options=ops) state_names = vcat( add_suffix(prob1.trajectory.state_names, "1")..., add_suffix(prob2.trajectory.state_names, "2")... @@ -302,6 +298,7 @@ end 0.99, prob_labels=["a", "b"], graph=[("a", "b")], + verbose=false, ipopt_options=ops) state_names_ab = vcat( add_suffix(prob1.trajectory.state_names, "a")..., @@ -320,6 +317,7 @@ end 0.99, prob_labels=["a", "b"], graph=[("x", "b")], + verbose=false, ipopt_options=ops ) @@ -328,6 +326,7 @@ end [prob1, prob2], 0.99, graph=[(:a1, :a2)], + verbose=false, ipopt_options=ops ) @test issetequal(direct_sum_prob3.trajectory.state_names, state_names) @@ -348,6 +347,7 @@ end R_b=1e3, Q_symb=:dda, boundary_values=Dict("x"=>copy(prob1.trajectory[:dda])), + verbose=false, ipopt_options=ops ) # # TODO: Check for objectives? From 75c924d2af9d200e5a103bb0d32862163cc343f3 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 16:44:38 -0500 Subject: [PATCH 24/35] refactor tests into testitems and folders --- src/integrators/_integrators.jl | 1 + src/integrators/exponential_integrators.jl | 95 +++++++++++++++++++ .../_quantum_system_templates.jl | 1 + .../quantum_optics.jl | 20 +++- test/integrators_test.jl | 92 ------------------ test/runtests.jl | 7 +- 6 files changed, 117 insertions(+), 99 deletions(-) delete mode 100644 test/integrators_test.jl diff --git a/src/integrators/_integrators.jl b/src/integrators/_integrators.jl index f17db1dd..590086ff 100644 --- a/src/integrators/_integrators.jl +++ b/src/integrators/_integrators.jl @@ -35,6 +35,7 @@ using TrajectoryIndexingUtils using LinearAlgebra using SparseArrays using ForwardDiff +using TestItemRunner abstract type AbstractIntegrator end diff --git a/src/integrators/exponential_integrators.jl b/src/integrators/exponential_integrators.jl index 9c357345..68386db3 100644 --- a/src/integrators/exponential_integrators.jl +++ b/src/integrators/exponential_integrators.jl @@ -226,3 +226,98 @@ end return ∂ψ̃ₜℰ, ∂ψ̃ₜ₊₁ℰ, ∂aₜℰ end end + +# ******************************************************************************* # + +@testitem "testing UnitaryExponentialIntegrator" begin + using NamedTrajectories + using ForwardDiff + + T = 100 + H_drift = GATES[:Z] + H_drives = [GATES[:X], GATES[:Y]] + n_drives = length(H_drives) + + system = QuantumSystem(H_drift, H_drives) + + U_init = GATES[:I] + U_goal = GATES[:X] + + Ũ⃗_init = operator_to_iso_vec(U_init) + Ũ⃗_goal = operator_to_iso_vec(U_goal) + + dt = 0.1 + + Z = NamedTrajectory( + ( + Ũ⃗ = unitary_geodesic(U_goal, T), + a = randn(n_drives, T), + da = randn(n_drives, T), + Δt = fill(dt, 1, T), + ), + controls=(:da,), + timestep=:Δt, + goal=(Ũ⃗ = Ũ⃗_goal,) + ) + + ℰ = UnitaryExponentialIntegrator(system, :Ũ⃗, :a) + + + ∂Ũ⃗ₜℰ, ∂Ũ⃗ₜ₊₁ℰ, ∂aₜℰ, ∂Δtₜℰ = jacobian(ℰ, Z[1].data, Z[2].data, Z) + + ∂ℰ_forwarddiff = ForwardDiff.jacobian( + zz -> ℰ(zz[1:Z.dim], zz[Z.dim+1:end], Z), + [Z[1].data; Z[2].data] + ) + + @test ∂Ũ⃗ₜℰ ≈ ∂ℰ_forwarddiff[:, 1:ℰ.dim] + @test ∂Ũ⃗ₜ₊₁ℰ ≈ ∂ℰ_forwarddiff[:, Z.dim .+ (1:ℰ.dim)] + @test ∂aₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.a] + @test ∂Δtₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.Δt] +end + +@testitem "testing QuantumStateExponentialIntegrator" begin + using NamedTrajectories + using ForwardDiff + + T = 100 + H_drift = GATES[:Z] + H_drives = [GATES[:X], GATES[:Y]] + n_drives = length(H_drives) + + system = QuantumSystem(H_drift, H_drives) + + U_init = GATES[:I] + U_goal = GATES[:X] + + ψ̃_init = ket_to_iso(quantum_state("g", [2])) + ψ̃_goal = ket_to_iso(quantum_state("e", [2])) + + dt = 0.1 + + Z = NamedTrajectory( + ( + ψ̃ = linear_interpolation(ψ̃_init, ψ̃_goal, T), + a = randn(n_drives, T), + da = randn(n_drives, T), + Δt = fill(dt, 1, T), + ), + controls=(:da,), + timestep=:Δt, + goal=(ψ̃ = ψ̃_goal,) + ) + + ℰ = QuantumStateExponentialIntegrator(system, :ψ̃, :a) + + ∂ψ̃ₜℰ, ∂ψ̃ₜ₊₁ℰ, ∂aₜℰ, ∂Δtₜℰ = jacobian(ℰ, Z[1].data, Z[2].data, Z) + + ∂ℰ_forwarddiff = ForwardDiff.jacobian( + zz -> ℰ(zz[1:Z.dim], zz[Z.dim+1:end], Z), + [Z[1].data; Z[2].data] + ) + + @test ∂ψ̃ₜℰ ≈ ∂ℰ_forwarddiff[:, 1:ℰ.dim] + @test ∂ψ̃ₜ₊₁ℰ ≈ ∂ℰ_forwarddiff[:, Z.dim .+ (1:ℰ.dim)] + @test ∂aₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.a] + @test ∂Δtₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.Δt] +end diff --git a/src/quantum_system_templates/_quantum_system_templates.jl b/src/quantum_system_templates/_quantum_system_templates.jl index cec2a352..9080dc18 100644 --- a/src/quantum_system_templates/_quantum_system_templates.jl +++ b/src/quantum_system_templates/_quantum_system_templates.jl @@ -10,6 +10,7 @@ using ..QuantumUtils using ..QuantumSystems using LinearAlgebra +using TestItemRunner include("transmons.jl") include("rydberg.jl") diff --git a/src/quantum_system_templates/quantum_optics.jl b/src/quantum_system_templates/quantum_optics.jl index 56f3ee00..01199d27 100644 --- a/src/quantum_system_templates/quantum_optics.jl +++ b/src/quantum_system_templates/quantum_optics.jl @@ -24,4 +24,22 @@ function QuantumOpticsSystem( H_drives; constructor=QuantumOpticsSystem, ) -end \ No newline at end of file +end + +# ******************************************************************************* # + +@testitem "Quantum Optics System" begin + using QuantumOpticsBase + N = rand(1:5); + T = ComplexF64; + b = FockBasis(N); + a = QuantumOpticsBase.create(T, b); + H = a + a'; + sys = QuantumOpticsSystem(H, [H, H]); + @test typeof(sys) == QuantumSystem + @test sys.constructor == QuantumOpticsSystem + @test sys.H_drift == H.data + + # creation with non-Hermitian operators + @test_broken QuantumOpticsSystem(a, [a]) +end diff --git a/test/integrators_test.jl b/test/integrators_test.jl deleted file mode 100644 index df1a23d3..00000000 --- a/test/integrators_test.jl +++ /dev/null @@ -1,92 +0,0 @@ -@testitem "testing UnitaryExponentialIntegrator" begin - using NamedTrajectories - using ForwardDiff - - T = 100 - H_drift = GATES[:Z] - H_drives = [GATES[:X], GATES[:Y]] - n_drives = length(H_drives) - - system = QuantumSystem(H_drift, H_drives) - - U_init = GATES[:I] - U_goal = GATES[:X] - - Ũ⃗_init = operator_to_iso_vec(U_init) - Ũ⃗_goal = operator_to_iso_vec(U_goal) - - dt = 0.1 - - Z = NamedTrajectory( - ( - Ũ⃗ = unitary_geodesic(U_goal, T), - a = randn(n_drives, T), - da = randn(n_drives, T), - Δt = fill(dt, 1, T), - ), - controls=(:da,), - timestep=:Δt, - goal=(Ũ⃗ = Ũ⃗_goal,) - ) - - ℰ = UnitaryExponentialIntegrator(system, :Ũ⃗, :a) - - - ∂Ũ⃗ₜℰ, ∂Ũ⃗ₜ₊₁ℰ, ∂aₜℰ, ∂Δtₜℰ = jacobian(ℰ, Z[1].data, Z[2].data, Z) - - ∂ℰ_forwarddiff = ForwardDiff.jacobian( - zz -> ℰ(zz[1:Z.dim], zz[Z.dim+1:end], Z), - [Z[1].data; Z[2].data] - ) - - @test ∂Ũ⃗ₜℰ ≈ ∂ℰ_forwarddiff[:, 1:ℰ.dim] - @test ∂Ũ⃗ₜ₊₁ℰ ≈ ∂ℰ_forwarddiff[:, Z.dim .+ (1:ℰ.dim)] - @test ∂aₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.a] - @test ∂Δtₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.Δt] -end - -@testitem "testing QuantumStateExponentialIntegrator" begin - using NamedTrajectories - using ForwardDiff - - T = 100 - H_drift = GATES[:Z] - H_drives = [GATES[:X], GATES[:Y]] - n_drives = length(H_drives) - - system = QuantumSystem(H_drift, H_drives) - - U_init = GATES[:I] - U_goal = GATES[:X] - - ψ̃_init = ket_to_iso(quantum_state("g", [2])) - ψ̃_goal = ket_to_iso(quantum_state("e", [2])) - - dt = 0.1 - - Z = NamedTrajectory( - ( - ψ̃ = linear_interpolation(ψ̃_init, ψ̃_goal, T), - a = randn(n_drives, T), - da = randn(n_drives, T), - Δt = fill(dt, 1, T), - ), - controls=(:da,), - timestep=:Δt, - goal=(ψ̃ = ψ̃_goal,) - ) - - ℰ = QuantumStateExponentialIntegrator(system, :ψ̃, :a) - - ∂ψ̃ₜℰ, ∂ψ̃ₜ₊₁ℰ, ∂aₜℰ, ∂Δtₜℰ = jacobian(ℰ, Z[1].data, Z[2].data, Z) - - ∂ℰ_forwarddiff = ForwardDiff.jacobian( - zz -> ℰ(zz[1:Z.dim], zz[Z.dim+1:end], Z), - [Z[1].data; Z[2].data] - ) - - @test ∂ψ̃ₜℰ ≈ ∂ℰ_forwarddiff[:, 1:ℰ.dim] - @test ∂ψ̃ₜ₊₁ℰ ≈ ∂ℰ_forwarddiff[:, Z.dim .+ (1:ℰ.dim)] - @test ∂aₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.a] - @test ∂Δtₜℰ ≈ ∂ℰ_forwarddiff[:, Z.components.Δt] -end diff --git a/test/runtests.jl b/test/runtests.jl index 804eaae3..71aad2f1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,10 +14,5 @@ using NamedTrajectories include("test_utils.jl") -@testset "QuantumCollocation.jl" begin - # include("objectives_test.jl") - # include("dynamics_test.jl") - include("quantum_system_templates_test.jl") -end - +# Run testitem @run_package_tests From 5889cfe71475cefb2ea94338ca1fee685bad447b Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 16:49:12 -0500 Subject: [PATCH 25/35] update readme --- README.md | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/README.md b/README.md index 3392484d..dde910a1 100644 --- a/README.md +++ b/README.md @@ -48,3 +48,41 @@ See the example script [examples/scripts/single_qubit_gate.jl](examples/scripts/ ![Single Qubit X-Gate](images/T_100_Q_1000_iter_1000_00004_fidelity_0.9999999999994745.png) +## Quickstart developers guide + +__Install Julia__ [Juliaup](https://github.com/JuliaLang/juliaup) is an installer and version manager. This is one useful way to manage Julia versions and keep up with the latest changes. After installing, run `julia` to obtain the Julia _REPL_. + +__Julia environments__ +[(Documentation)](https://pkgdocs.julialang.org/v1/environments/#Using-someone-else's-project) Your project's environment is stored in _Project.toml_. You can interactively add packages to an environment by using the Julia command line _REPL_ and _package manager_. Start Julia in the project folder. Type `]` to enter the package manager. Type `activate .` to activate or create an environment specified by _Project.toml_ located in the current folder. Separately, you generate a manifest (solving the versions to create a valid environment) by running `instantiate`; instantiate will check that the environment is correct after you add all the packages you want. + +__Adding packages__ +[(Documentation)](https://pkgdocs.julialang.org/v1/managing-packages/#Adding-packages) The initial cell for a Piccolo notebook might look something like the following: +```Julia +using QuantumCollocation +using NamedTrajectories +using TrajectoryIndexingUtils + +using LinearAlgebra +using CairoMakie +``` + +First, let's install a standard library packages. From the package manager and in the current environment (type `julia`, `]`, and `activate .`), you can type `add LinearAlgebra` to install and precompile _LinearAlgebra_. Same with `CairoMakie`. These are like Numpy and Matplotlib. + +Second, let's install _Piccolo_. The first three packages (_QuantumCollocation_, _NamedTrajetories_, _TrajectoryIndexingUtils_) are the core of [Piccolo](https://docs.juliahub.com/General/Piccolo/stable/). We could do `add Piccolo` to get the three as a bundle from the Julia repository, which requires only `using Piccolo`. + +As a developer, we want to use the git repositories directly from [our Github page](https://github.com/aarontrowbridge). Clone, then add to the Project file with e.g. `dev ../relative/path/to/repo/QuantumCollocation` to obtain a development version of _QuantumCollocation_ pointing to the local Github code instead of the package repository. You can repeat this for the others, also. + +__Developing__ +[Revise.jl](https://timholy.github.io/Revise.jl/stable/) will let you edit source code, update packages, and reload the changes in your notebook---automatically! This is a great tool for development. `add Revise` from the REPL and then include it before any packages you intend to edit: +```Julia +using Revise +using QuantumCollocation +``` + +### Tips for Visual Studio Code +__Julia extension__ You can run Julia notebooks and much more with [the Julia extension](https://code.visualstudio.com/docs/languages/julia). Upon opening your project folder in VS code and attempting to run an `.ipynb`, you will see that VS Code finds the interpreters managed by juliaup and defaults to using the environment based on the _Project.toml_ in the project directory. + +__Fonts__ VS Code will not display all characters allowed by Julia. You can change the editor font family in the settings to `'JuliaMono'` to get full support. If you don't want to mix and mash, you can create a new VS Code settings profile for working in Julia at _File>Preferences>Profile_. + +__Tests__ Tests should automatically populate in VS Code when working with a Piccolo package. For example, just by adding the `QuantumCollocation.jl` folder to your workspace, you should see tests appear if you click on the _Testing_ sidebar icon. If you run one of these tests, a new Julia kernel is spawned for the test. You can find the kernel if you click on the _Julia_ sidebar icon (after installing the Julia extensions). Sometimes, for the tests to recognize new changes, you may need to manually kill this kernel to see your changes reflected. + From a522cd2cf6f79deb55477fb6ce3effadceba34d9 Mon Sep 17 00:00:00 2001 From: Andy Goldschmidt Date: Mon, 24 Jun 2024 18:49:20 -0500 Subject: [PATCH 26/35] Update README.md repos link --- README.md | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index dde910a1..e5ee2436 100644 --- a/README.md +++ b/README.md @@ -58,22 +58,24 @@ __Julia environments__ __Adding packages__ [(Documentation)](https://pkgdocs.julialang.org/v1/managing-packages/#Adding-packages) The initial cell for a Piccolo notebook might look something like the following: ```Julia +# Standard packages +using LinearAlgebra +using CairoMakie + +# Piccolo packages using QuantumCollocation using NamedTrajectories using TrajectoryIndexingUtils - -using LinearAlgebra -using CairoMakie ``` -First, let's install a standard library packages. From the package manager and in the current environment (type `julia`, `]`, and `activate .`), you can type `add LinearAlgebra` to install and precompile _LinearAlgebra_. Same with `CairoMakie`. These are like Numpy and Matplotlib. +First, let's install some standard packages (these are like Numpy and Matplotlib). Open the package manager in the current environment (type `julia`, `]`, and `activate .`), type `add LinearAlgebra` to install and precompile _LinearAlgebra_. Same with `CairoMakie`. -Second, let's install _Piccolo_. The first three packages (_QuantumCollocation_, _NamedTrajetories_, _TrajectoryIndexingUtils_) are the core of [Piccolo](https://docs.juliahub.com/General/Piccolo/stable/). We could do `add Piccolo` to get the three as a bundle from the Julia repository, which requires only `using Piccolo`. +Second, let's install _Piccolo_. There are three packages (_QuantumCollocation_, _NamedTrajetories_, _TrajectoryIndexingUtils_) inside [Piccolo](https://docs.juliahub.com/General/Piccolo/stable/). We could do `add Piccolo` to get the three as a bundle from the Julia repository. Instead of individually calling `using ...` for each, this approach only requires `using Piccolo` at the start of a file or notebook. -As a developer, we want to use the git repositories directly from [our Github page](https://github.com/aarontrowbridge). Clone, then add to the Project file with e.g. `dev ../relative/path/to/repo/QuantumCollocation` to obtain a development version of _QuantumCollocation_ pointing to the local Github code instead of the package repository. You can repeat this for the others, also. +As a developer, we want to use the git repositories directly from [the Kestrel Quantum Github page](https://github.com/kestrelquantum). Clone, then add the local packages to the Project file with e.g. `dev ../relative/path/to/repo/QuantumCollocation`. This command installs the development version of _QuantumCollocation_ pointing to the local Github code instead of the package repository. You can repeat this for the others, also. __Developing__ -[Revise.jl](https://timholy.github.io/Revise.jl/stable/) will let you edit source code, update packages, and reload the changes in your notebook---automatically! This is a great tool for development. `add Revise` from the REPL and then include it before any packages you intend to edit: +[Revise.jl](https://timholy.github.io/Revise.jl/stable/) will let you edit source code, update packages, and reload the changes in a notebook---automatically! This is a great tool for development. `add Revise` from the REPL and then include it before any packages you intend to edit: ```Julia using Revise using QuantumCollocation From d4f504247993a3729bbb221f3fe0df9ebbae05f0 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 20:57:34 -0500 Subject: [PATCH 27/35] change smooth pulse template to use piccolo options --- Manifest.toml | 6 +- docs/src/contribution_guide.md | 2 +- src/QuantumCollocation.jl | 4 +- src/dynamics.jl | 29 +++--- src/evaluators.jl | 4 +- src/integrators/pade_integrators.jl | 8 +- src/{ipopt_options.jl => options.jl} | 30 +++++- src/problem_solvers.jl | 12 +-- src/problem_templates/_problem_templates.jl | 2 +- .../quantum_state_minimum_time_problem.jl | 2 +- .../quantum_state_smooth_pulse_problem.jl | 6 +- .../unitary_direct_sum_problem.jl | 4 +- .../unitary_minimum_time_problem.jl | 8 +- .../unitary_robustness_problem.jl | 10 +- .../unitary_sampling_problem.jl | 6 +- .../unitary_smooth_pulse_problem.jl | 76 ++++++-------- src/problems.jl | 99 ++++++++++--------- test/problems_test.jl | 4 +- 18 files changed, 160 insertions(+), 152 deletions(-) rename src/{ipopt_options.jl => options.jl} (71%) diff --git a/Manifest.toml b/Manifest.toml index 4fd4a8d6..8ac99798 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.3" +julia_version = "1.10.0" manifest_format = "2.0" project_hash = "0af13b9126caf21a7349df842fd536f87be35b2f" @@ -338,7 +338,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.1+0" +version = "1.0.5+1" [[deps.CompositeTypes]] git-tree-sha1 = "bce26c3dab336582805503bed209faab1c279768" @@ -1542,7 +1542,7 @@ version = "0.3.24+0" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+4" +version = "0.3.23+2" [[deps.OpenEXR]] deps = ["Colors", "FileIO", "OpenEXR_jll"] diff --git a/docs/src/contribution_guide.md b/docs/src/contribution_guide.md index 8abd880f..2341baff 100644 --- a/docs/src/contribution_guide.md +++ b/docs/src/contribution_guide.md @@ -63,7 +63,7 @@ Tests are implemented using the [`TestItemRunner.jl`](https://github.com/julia-v prob = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=Options(print_level=1) + ipopt_options=IpoptOptions(print_level=1) ) solve!(prob, max_iter=100) diff --git a/src/QuantumCollocation.jl b/src/QuantumCollocation.jl index 9bc43028..8dfed574 100644 --- a/src/QuantumCollocation.jl +++ b/src/QuantumCollocation.jl @@ -38,8 +38,8 @@ include("dynamics.jl") include("evaluators.jl") @reexport using .Evaluators -include("ipopt_options.jl") -@reexport using .IpoptOptions +include("options.jl") +@reexport using .Options include("problems.jl") @reexport using .Problems diff --git a/src/dynamics.jl b/src/dynamics.jl index f9c07691..17db4bf5 100644 --- a/src/dynamics.jl +++ b/src/dynamics.jl @@ -199,7 +199,7 @@ end function QuantumDynamics( integrators::Vector{<:AbstractIntegrator}, traj::NamedTrajectory; - hessian_approximation=false, + eval_hessian=true, jacobian_structure=true, verbose=false ) @@ -235,10 +235,10 @@ function QuantumDynamics( ∂f = dynamics_jacobian(integrators, traj) - if hessian_approximation - μ∂²f = nothing - else + if eval_hessian μ∂²f = dynamics_hessian_of_lagrangian(integrators, traj) + else + μ∂²f = nothing end if verbose @@ -247,13 +247,7 @@ function QuantumDynamics( dynamics_dim = dim(integrators) - if hessian_approximation - ∂f_structure, ∂F_structure = dynamics_structure(∂f, traj, dynamics_dim; - verbose=verbose, - jacobian=jacobian_structure, - ) - μ∂²F_structure = nothing - else + if eval_hessian ∂f_structure, ∂F_structure, μ∂²f_structure, μ∂²F_structure = dynamics_structure(∂f, μ∂²f, traj, dynamics_dim; verbose=verbose, @@ -263,6 +257,12 @@ function QuantumDynamics( ) ) μ∂²f_nnz = length(μ∂²f_structure) + else + ∂f_structure, ∂F_structure = dynamics_structure(∂f, traj, dynamics_dim; + verbose=verbose, + jacobian=jacobian_structure, + ) + μ∂²F_structure = nothing end ∂f_nnz = length(∂f_structure) @@ -294,9 +294,7 @@ function QuantumDynamics( return ∂s end - if hessian_approximation - μ∂²F = nothing - else + if eval_hessian @views μ∂²F = (Z⃗::AbstractVector{<:Real}, μ⃗::AbstractVector{<:Real}) -> begin μ∂²s = Vector{eltype(Z⃗)}(undef, length(μ∂²F_structure)) Threads.@threads for t = 1:traj.T-1 @@ -310,7 +308,10 @@ function QuantumDynamics( end return μ∂²s end + else + μ∂²F = nothing end + return QuantumDynamics( integrators, F, diff --git a/src/evaluators.jl b/src/evaluators.jl index 50cba1fe..7f84882c 100644 --- a/src/evaluators.jl +++ b/src/evaluators.jl @@ -26,8 +26,8 @@ mutable struct PicoEvaluator <: MOI.AbstractNLPEvaluator trajectory::NamedTrajectory, objective::Objective, dynamics::QuantumDynamics, - nonlinear_constraints::Vector{<:NonlinearConstraint}, - eval_hessian::Bool + nonlinear_constraints::Vector{<:NonlinearConstraint}; + eval_hessian::Bool=true ) n_dynamics_constraints = dynamics.dim * (trajectory.T - 1) n_nonlinear_constraints = sum(con.dim for con ∈ nonlinear_constraints; init=0) diff --git a/src/integrators/pade_integrators.jl b/src/integrators/pade_integrators.jl index aafd2ba2..6eab2f03 100644 --- a/src/integrators/pade_integrators.jl +++ b/src/integrators/pade_integrators.jl @@ -96,8 +96,8 @@ struct UnitaryPadeIntegrator <: QuantumPadeIntegrator unitary_symb::Union{Symbol,Nothing}=nothing, drive_symb::Union{Symbol,Tuple{Vararg{Symbol}},Nothing}=nothing; order::Int=4, - autodiff::Bool=false, - G::Union{Function, Nothing}=nothing, + autodiff::Bool=order != 4, + G::Function=G_bilinear, ) @assert order ∈ [4, 6, 8, 10] "order must be in [4, 6, 8, 10]" @assert !isnothing(unitary_symb) "must specify unitary symbol" @@ -111,7 +111,6 @@ struct UnitaryPadeIntegrator <: QuantumPadeIntegrator G_drift = sys.G_drift G_drives = sys.G_drives - G = isnothing(G) ? G_bilinear : G drive_anticomms, drift_anticomms = order == 4 ? build_anticomms(G_drift, G_drives, n_drives) : (nothing, nothing) @@ -189,7 +188,7 @@ struct QuantumStatePadeIntegrator <: QuantumPadeIntegrator drive_symb::Union{Symbol,Tuple{Vararg{Symbol}},Nothing}=nothing; order::Int=4, autodiff::Bool=order != 4, - G::Union{Function, Nothing}=nothing, + G::Function=G_bilinear, ) @assert order ∈ [4, 6, 8, 10] "order must be in [4, 6, 8, 10]" @assert !isnothing(state_symb) "state_symb must be specified" @@ -201,7 +200,6 @@ struct QuantumStatePadeIntegrator <: QuantumPadeIntegrator G_drift = sys.G_drift G_drives = sys.G_drives - G = isnothing(G) ? G_bilinear : G drive_anticomms, drift_anticomms = order == 4 ? build_anticomms(G_drift, G_drives, n_drives) : (nothing, nothing) diff --git a/src/ipopt_options.jl b/src/options.jl similarity index 71% rename from src/ipopt_options.jl rename to src/options.jl index be41e736..db455e82 100644 --- a/src/ipopt_options.jl +++ b/src/options.jl @@ -1,18 +1,21 @@ -module IpoptOptions +module Options -export Options +export IpoptOptions +export PiccoloOptions export set! using Ipopt using Libdl using Base: @kwdef +abstract type AbstractOptions end + """ Solver options for Ipopt https://coin-or.github.io/Ipopt/OPTIONS.html#OPT_print_options_documentation """ -@kwdef mutable struct Options{T} +@kwdef mutable struct IpoptOptions{T} <: AbstractOptions tol::T = 1e-8 s_max::T = 100.0 max_iter::Int = 1_000 @@ -49,7 +52,26 @@ using Base: @kwdef watchdog_trial_iter_max = 3 end -function set!(optimizer::Ipopt.Optimizer, options::Options) + +""" + Solver settings for Quantum Collocation. +""" +@kwdef mutable struct PiccoloOptions <: AbstractOptions + verbose::Bool = true + free_time::Bool = true + timesteps_all_equal::Bool = true + integrator::Symbol = :pade + pade_order::Int = 4 + eval_hessian::Bool = false + jacobian_structure::Bool = true + rollout_integrator::Function = exp + geodesic = true + blas_multithreading::Bool = true + build_trajectory_constraints::Bool = true +end + + +function set!(optimizer::Ipopt.Optimizer, options::AbstractOptions) for name in fieldnames(typeof(options)) value = getfield(options, name) if name == :linear_solver diff --git a/src/problem_solvers.jl b/src/problem_solvers.jl index 651b8b1b..8b47ac06 100644 --- a/src/problem_solvers.jl +++ b/src/problem_solvers.jl @@ -5,7 +5,7 @@ export solve! using ..Constraints using ..Problems using ..SaveLoadUtils -using ..IpoptOptions +using ..Options using NamedTrajectories using MathOptInterface @@ -16,13 +16,13 @@ function solve!( init_traj=nothing, save_path=nothing, controls_save_path=nothing, - max_iter::Int=prob.options.max_iter, - linear_solver::String=prob.options.linear_solver, + max_iter::Int=prob.ipopt_options.max_iter, + linear_solver::String=prob.ipopt_options.linear_solver, ) - prob.options.max_iter = max_iter - prob.options.linear_solver = linear_solver + prob.ipopt_options.max_iter = max_iter + prob.ipopt_options.linear_solver = linear_solver - set!(prob.optimizer, prob.options) + set!(prob.optimizer, prob.ipopt_options) if !isnothing(init_traj) set_trajectory!(prob, init_traj) diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index 65aa9801..dff08dd4 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -19,7 +19,7 @@ using ..Objectives using ..Constraints using ..Integrators using ..Problems -using ..IpoptOptions +using ..Options using Distributions using NamedTrajectories diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index 18073f68..47510204 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -14,7 +14,7 @@ function QuantumStateMinimumTimeProblem( state_symbol::Symbol=:ψ̃, D=1.0, verbose::Bool=false, - ipopt_options::Options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), kwargs... ) @assert state_symbol ∈ trajectory.names diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 739db6d3..30464327 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -49,7 +49,7 @@ function QuantumStateSmoothPulseProblem( R_L1::Float64=20.0, max_iter::Int=1000, linear_solver::String="mumps", - ipopt_options::Options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], timesteps_all_equal::Bool=true, L1_regularized_names=Symbol[], @@ -175,7 +175,7 @@ end prob = QuantumStateSmoothPulseProblem( sys, ψ_init, ψ_target, T, Δt; - ipopt_options=Options(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), verbose=false ) initial = fidelity(prob) solve!(prob, max_iter=20) @@ -187,7 +187,7 @@ end ψ_targets = [[0.0, 1.0], [1.0, 0.0]] prob = QuantumStateSmoothPulseProblem( sys, ψ_inits, ψ_targets, T, Δt; - ipopt_options=Options(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), verbose=false ) initial = fidelity(prob) solve!(prob, max_iter=20) diff --git a/src/problem_templates/unitary_direct_sum_problem.jl b/src/problem_templates/unitary_direct_sum_problem.jl index c367f7f8..6971ad97 100644 --- a/src/problem_templates/unitary_direct_sum_problem.jl +++ b/src/problem_templates/unitary_direct_sum_problem.jl @@ -40,7 +40,7 @@ between each neighbor of the provided `probs`. - `hessian_approximation=true`: whether or not to use L-BFGS hessian approximation in Ipopt - `jacobian_structure=true`: whether or not to use the jacobian structure in Ipopt - `blas_multithreading=true`: whether or not to use multithreading in BLAS -- `ipopt_options::Options=Options()`: the options for the Ipopt solver +- `ipopt_options::IpoptOptions=IpoptOptions()`: the options for the Ipopt solver """ function UnitaryDirectSumProblem( @@ -66,7 +66,7 @@ function UnitaryDirectSumProblem( hessian_approximation=true, jacobian_structure=true, blas_multithreading=true, - ipopt_options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), kwargs... ) N = length(probs) diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index 9853b5f8..9b912782 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -42,7 +42,7 @@ J(\vec{\tilde{U}}, a, \dot{a}, \ddot{a}) + D \sum_t \Delta t_t \\ - `final_fidelity::Float64=0.99`: The final fidelity. - `D=1.0`: The weight for the minimum-time objective. - `verbose::Bool=false`: Whether to print additional information. -- `ipopt_options::Options=Options()`: The options for the Ipopt solver. +- `ipopt_options::IpoptOptions=IpoptOptions()`: The options for the Ipopt solver. - `kwargs...`: Additional keyword arguments to pass to `QuantumControlProblem`. """ function UnitaryMinimumTimeProblem end @@ -57,7 +57,7 @@ function UnitaryMinimumTimeProblem( final_fidelity::Float64=0.99, D=1.0, verbose::Bool=false, - ipopt_options::Options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), subspace=nothing, kwargs... ) @@ -145,7 +145,7 @@ end prob = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=Options(print_level=1) + ipopt_options=IpoptOptions(print_level=1) ) solve!(prob, max_iter=100) @@ -157,7 +157,7 @@ end mintime_prob = UnitaryMinimumTimeProblem( prob, final_fidelity=final_fidelity, - ipopt_options=Options(print_level=1) + ipopt_options=IpoptOptions(print_level=1) ) solve!(mintime_prob; max_iter=100) diff --git a/src/problem_templates/unitary_robustness_problem.jl b/src/problem_templates/unitary_robustness_problem.jl index 9dc576ae..0d75368a 100644 --- a/src/problem_templates/unitary_robustness_problem.jl +++ b/src/problem_templates/unitary_robustness_problem.jl @@ -5,7 +5,7 @@ subspace=nothing, eval_hessian=false, verbose=false, - ipopt_options=Options(), + ipopt_options=IpoptOptions(), kwargs... ) @@ -28,7 +28,7 @@ function UnitaryRobustnessProblem( subspace::AbstractVector{<:Integer}=1:size(Hₑ, 1), eval_hessian::Bool=false, verbose::Bool=false, - ipopt_options::Options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), kwargs... ) @assert unitary_symbol ∈ trajectory.names @@ -118,7 +118,7 @@ end sys, U_goal, T, Δt, geodesic=false, verbose=false, - ipopt_options=Options(print_level=1) + ipopt_options=IpoptOptions(print_level=1) ) solve!(probs["transmon"], max_iter=200) @@ -134,7 +134,7 @@ end final_fidelity=0.99, subspace=subspace, verbose=false, - ipopt_options=Options(recalc_y="yes", recalc_y_feas_tol=1.0, print_level=1) + ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=1.0, print_level=1) ) solve!(probs["robust"], max_iter=200) @@ -161,7 +161,7 @@ end H_error, trajectory, system, objective, integrators, constraints, final_fidelity=0.99, subspace=subspace, - ipopt_options=Options(recalc_y="yes", recalc_y_feas_tol=1e-1, print_level=4) + ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=1e-1, print_level=4) ) solve!(probs["unconstrained"]; max_iter=100) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index 4fb8f2fe..bdbf38f6 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -76,7 +76,7 @@ function UnitarySamplingProblem( R_leakage=1e-1, max_iter::Int=1000, linear_solver::String="mumps", - ipopt_options::Options=Options(), + ipopt_options::IpoptOptions=IpoptOptions(), constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], timesteps_all_equal::Bool=true, verbose::Bool=false, @@ -258,7 +258,7 @@ end T, Δt; verbose=false, - ipopt_options=Options(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + ipopt_options=IpoptOptions(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) ) solve!(prob, max_iter=20) @@ -269,7 +269,7 @@ end T, Δt; verbose=false, - ipopt_options=Options(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + ipopt_options=IpoptOptions(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) ) solve!(prob, max_iter=20) diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index a2803740..1e884e70 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -56,7 +56,7 @@ with - `R_leakage=1e-1`: the weight on the leakage suppression term - `max_iter::Int=1000`: the maximum number of iterations for the solver - `linear_solver::String="mumps"`: the linear solver to use -- `ipopt_options::Options=Options()`: the options for the Ipopt solver +- `ipopt_options::IpoptOptions=IpoptOptions()`: the options for the Ipopt solver - `constraints::Vector{<:AbstractConstraint}=AbstractConstraint[]`: additional constraints to add to the problem - `timesteps_all_equal::Bool=true`: whether or not to enforce that all time steps are equal - `verbose::Bool=false`: whether or not to print constructor output @@ -72,14 +72,20 @@ with - `subspace=nothing`: the subspace to use for the unitary integrator - `jacobian_structure=true`: whether or not to use the jacobian structure - `hessian_approximation=false`: whether or not to use L-BFGS hessian approximation in Ipopt -- `blas_multithreading=true`: whether or not to use multithreading in BLAS +- `blas_multithreading=true`: whether or not to use multithreading in BLAS\\ + + +TODO: control modulus norm, advanced feature, needs documentation + """ function UnitarySmoothPulseProblem( system::AbstractQuantumSystem, operator::Union{EmbeddedOperator, AbstractMatrix{<:Number}}, T::Int, Δt::Union{Float64, Vector{Float64}}; - free_time=true, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], init_trajectory::Union{NamedTrajectory, Nothing}=nothing, a_bound::Float64=1.0, a_bounds=fill(a_bound, length(system.G_drives)), @@ -96,35 +102,12 @@ function UnitarySmoothPulseProblem( R_dda::Union{Float64, Vector{Float64}}=R, leakage_suppression=false, R_leakage=1e-1, - max_iter::Int=1000, - linear_solver::String="mumps", - ipopt_options::Options=Options(), - constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - timesteps_all_equal::Bool=true, - verbose::Bool=false, - integrator::Symbol=:pade, - rollout_integrator=exp, - bound_unitary=integrator == :exponential, - # TODO: control modulus norm, advanced feature, needs documentation control_norm_constraint=false, control_norm_constraint_components=nothing, control_norm_R=nothing, - geodesic=true, - pade_order=4, - autodiff=pade_order != 4, - jacobian_structure=true, - hessian_approximation=false, - blas_multithreading=true, + bound_unitary=piccolo_options.integrator == :exponential, kwargs... ) - if !blas_multithreading - BLAS.set_num_threads(1) - end - - if hessian_approximation - ipopt_options.hessian_approximation = "limited-memory" - end - # Trajectory if !isnothing(init_trajectory) traj = init_trajectory @@ -137,14 +120,14 @@ function UnitarySmoothPulseProblem( n_drives, a_bounds, dda_bounds; - free_time=free_time, + free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), - geodesic=geodesic, + geodesic=piccolo_options.geodesic, bound_unitary=bound_unitary, drive_derivative_σ=drive_derivative_σ, a_guess=a_guess, system=system, - rollout_integrator=rollout_integrator, + rollout_integrator=piccolo_options.rollout_integrator, ) end @@ -164,15 +147,15 @@ function UnitarySmoothPulseProblem( constraints, :Ũ⃗, traj, R_value=R_leakage, indices=leakage_indices, - eval_hessian=!hessian_approximation + eval_hessian=piccolo_options.eval_hessian ) else @warn "leakage_suppression is not supported for non-embedded operators, ignoring." end end - if free_time - if timesteps_all_equal + if piccolo_options.free_time + if piccolo_options.timesteps_all_equal push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) end end @@ -190,10 +173,10 @@ function UnitarySmoothPulseProblem( end # Integrators - if integrator == :pade + if piccolo_options.integrator == :pade unitary_integrator = - UnitaryPadeIntegrator(system, :Ũ⃗, :a; order=pade_order, autodiff=autodiff) - elseif integrator == :exponential + UnitaryPadeIntegrator(system, :Ũ⃗, :a; order=piccolo_options.pade_order) + elseif piccolo_options.integrator == :exponential unitary_integrator = UnitaryExponentialIntegrator(system, :Ũ⃗, :a) else @@ -212,13 +195,8 @@ function UnitarySmoothPulseProblem( J, integrators; constraints=constraints, - max_iter=max_iter, - linear_solver=linear_solver, - verbose=verbose, ipopt_options=ipopt_options, - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - eval_hessian=!hessian_approximation, + piccolo_options=piccolo_options, kwargs... ) end @@ -256,7 +234,8 @@ end prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - ipopt_options=Options(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) initial = unitary_fidelity(prob) @@ -272,17 +251,21 @@ end T = 51 Δt = 0.2 + # Test embedded operator + # ---------------------- prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - ipopt_options=Options(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) initial = unitary_fidelity(prob, subspace=U_goal.subspace_indices) solve!(prob, max_iter=20) final = unitary_fidelity(prob, subspace=U_goal.subspace_indices) @test final > initial - + # Test leakage suppression + # ------------------------ a = annihilate(4) sys = QuantumSystem([(a + a')/2, (a - a')/(2im)]) U_goal = EmbeddedOperator(GATES[:H], sys) @@ -292,7 +275,8 @@ end prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, leakage_suppression=true, R_leakage=1e-1, - ipopt_options=Options(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) initial = unitary_fidelity(prob, subspace=U_goal.subspace_indices) diff --git a/src/problems.jl b/src/problems.jl index 82bc728f..522637ec 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -13,7 +13,7 @@ export get_objective using ..QuantumSystems using ..Integrators using ..Evaluators -using ..IpoptOptions +using ..Options using ..Constraints using ..Dynamics using ..Objectives @@ -42,7 +42,7 @@ mutable struct QuantumControlProblem <: AbstractProblem system::AbstractQuantumSystem trajectory::NamedTrajectory integrators::Union{Nothing,Vector{<:AbstractIntegrator}} - options::Options + ipopt_options::IpoptOptions params::Dict{Symbol, Any} end @@ -51,43 +51,52 @@ function QuantumControlProblem( traj::NamedTrajectory, obj::Objective, dynamics::QuantumDynamics; + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), + scale_factor_objective::Float64=1.0, additional_objective::Union{Nothing, Objective}=nothing, - eval_hessian::Bool=true, - options::Options=Options(), constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], params::Dict{Symbol, Any}=Dict{Symbol, Any}(), - max_iter::Int=options.max_iter, - linear_solver::String=options.linear_solver, - verbose=false, - build_trajectory_constraints::Bool=true, kwargs... ) - options.max_iter = max_iter - options.linear_solver = linear_solver + if !piccolo_options.blas_multithreading + BLAS.set_num_threads(1) + end + + if !piccolo_options.eval_hessian + ipopt_options.hessian_approximation = "limited-memory" + end nonlinear_constraints = NonlinearConstraint[con for con ∈ constraints if con isa NonlinearConstraint] + if scale_factor_objective != 1 + obj = scale_factor_objective * obj + end + if !isnothing(additional_objective) obj += additional_objective end - if verbose + if piccolo_options.verbose println(" building evaluator...") end - evaluator = PicoEvaluator(traj, obj, dynamics, nonlinear_constraints, eval_hessian) + + evaluator = PicoEvaluator( + traj, obj, dynamics, nonlinear_constraints, eval_hessian=piccolo_options.eval_hessian + ) n_dynamics_constraints = dynamics.dim * (traj.T - 1) n_variables = traj.dim * traj.T linear_constraints = LinearConstraint[con for con ∈ constraints if con isa LinearConstraint] - if build_trajectory_constraints + if piccolo_options.build_trajectory_constraints linear_constraints = LinearConstraint[trajectory_constraints(traj); linear_constraints] end optimizer = Ipopt.Optimizer() - if verbose + if piccolo_options.verbose println(" initializing optimizer...") end @@ -99,15 +108,15 @@ function QuantumControlProblem( n_dynamics_constraints, nonlinear_constraints, n_variables, - verbose=verbose + verbose=piccolo_options.verbose ) variables = reshape(variables, traj.dim, traj.T) params = merge(kwargs, params) - params[:eval_hessian] = eval_hessian - params[:options] = options + params[:ipopt_options] = ipopt_options + params[:piccolo_options] = piccolo_options params[:linear_constraints] = linear_constraints params[:nonlinear_constraints] = [ nl_constraint.params for nl_constraint ∈ nonlinear_constraints @@ -120,7 +129,7 @@ function QuantumControlProblem( system, traj, dynamics.integrators, - options, + ipopt_options, params ) end @@ -131,28 +140,26 @@ function QuantumControlProblem( obj::Objective, integrators::Vector{<:AbstractIntegrator}; params::Dict{Symbol,Any}=Dict{Symbol, Any}(), - ipopt_options::Options=Options(), - verbose=false, - jacobian_structure=true, - hessian_approximation=false, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), kwargs... ) - if verbose + if piccolo_options.verbose println(" building dynamics from integrators...") end dynamics = QuantumDynamics(integrators, traj; - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - verbose=verbose + jacobian_structure=piccolo_options.jacobian_structure, + eval_hessian=piccolo_options.eval_hessian, + verbose=piccolo_options.verbose ) return QuantumControlProblem( system, traj, obj, dynamics; - options=ipopt_options, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, params=params, - verbose=verbose, kwargs... ) end @@ -164,28 +171,26 @@ function QuantumControlProblem( obj::Objective, integrator::AbstractIntegrator; params::Dict{Symbol,Any}=Dict{Symbol, Any}(), - ipopt_options::Options=Options(), - verbose=false, - jacobian_structure=true, - hessian_approximation=false, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), kwargs... ) - if verbose + if piccolo_options.verbose println(" building dynamics from integrator...") end dynamics = QuantumDynamics(integrator, traj; - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - verbose=verbose + jacobian_structure=piccolo_options.jacobian_structure, + eval_hessian=piccolo_options.eval_hessian, + verbose=piccolo_options.verbose ) return QuantumControlProblem( system, traj, obj, dynamics; - options=ipopt_options, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, params=params, - verbose=verbose, kwargs... ) end @@ -196,28 +201,26 @@ function QuantumControlProblem( obj::Objective, f::Function; params::Dict{Symbol,Any}=Dict{Symbol, Any}(), - ipopt_options::Options=Options(), - verbose=false, - jacobian_structure=true, - hessian_approximation=false, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), kwargs... ) - if verbose + if piccolo_options.verbose println(" building dynamics from function...") end dynamics = QuantumDynamics(f, traj; - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - verbose=verbose + jacobian_structure=piccolo_options.jacobian_structure, + eval_hessian=piccolo_options.eval_hessian, + verbose=piccolo_options.verbose ) return QuantumControlProblem( system, traj, obj, dynamics; - options=ipopt_options, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, params=params, - verbose=verbose, kwargs... ) end diff --git a/test/problems_test.jl b/test/problems_test.jl index 08399db3..ec13115c 100644 --- a/test/problems_test.jl +++ b/test/problems_test.jl @@ -30,14 +30,14 @@ end prob_vanilla = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=Options(print_level=4) + ipopt_options=IpoptOptions(print_level=4) ) J_extra = QuadraticSmoothnessRegularizer(:dda, prob_vanilla.trajectory, 10.0) prob_additional = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=Options(print_level=4), + ipopt_options=IpoptOptions(print_level=4), additional_objective=J_extra, ) From ed376c138a60361bbd11c39a94e0e67d4389f517 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 21:01:30 -0500 Subject: [PATCH 28/35] state smooth pulse piccolo options --- .../quantum_state_smooth_pulse_problem.jl | 28 ++++++++----------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 30464327..cb5fe00f 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -30,12 +30,12 @@ function QuantumStateSmoothPulseProblem( ψ_goal::Union{AbstractVector{<:Number}, Vector{<:AbstractVector{<:Number}}}, T::Int, Δt::Float64; - free_time=true, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), init_trajectory::Union{NamedTrajectory, Nothing}=nothing, a_bound::Float64=1.0, a_bounds::Vector{Float64}=fill(a_bound, length(system.G_drives)), a_guess::Union{Matrix{Float64}, Nothing}=nothing, - rollout_integrator=exp, dda_bound::Float64=1.0, dda_bounds::Vector{Float64}=fill(dda_bound, length(system.G_drives)), Δt_min::Float64=0.5 * Δt, @@ -47,14 +47,9 @@ function QuantumStateSmoothPulseProblem( R_da::Union{Float64, Vector{Float64}}=R, R_dda::Union{Float64, Vector{Float64}}=R, R_L1::Float64=20.0, - max_iter::Int=1000, - linear_solver::String="mumps", - ipopt_options::IpoptOptions=IpoptOptions(), constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - timesteps_all_equal::Bool=true, L1_regularized_names=Symbol[], L1_regularized_indices::NamedTuple=NamedTuple(), - verbose=false, kwargs... ) @assert all(name ∈ L1_regularized_names for name in keys(L1_regularized_indices) if !isempty(L1_regularized_indices[name])) @@ -85,12 +80,12 @@ function QuantumStateSmoothPulseProblem( n_drives, a_bounds, dda_bounds; - free_time=free_time, + free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), drive_derivative_σ=drive_derivative_σ, a_guess=a_guess, system=system, - rollout_integrator=rollout_integrator, + rollout_integrator=piccolo_options.rollout_integrator, ) end @@ -110,19 +105,19 @@ function QuantumStateSmoothPulseProblem( constraints, name, traj, R_value=R_L1, indices=L1_regularized_indices[name], - eval_hessian=!hessian_approximation + eval_hessian=piccolo_options.eval_hessian ) else J += L1Regularizer!( constraints, name, traj; R_value=R_L1, - eval_hessian=!hessian_approximation + eval_hessian=piccolo_options.eval_hessian ) end end - if free_time - if timesteps_all_equal + if piccolo_options.free_time + if piccolo_options.timesteps_all_equal push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) end end @@ -145,10 +140,8 @@ function QuantumStateSmoothPulseProblem( J, integrators; constraints=constraints, - max_iter=max_iter, - linear_solver=linear_solver, - verbose=verbose, ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end @@ -173,6 +166,8 @@ end ψ_init = [1.0, 0.0] ψ_target = [0.0, 1.0] + # Single initial and target states + # -------------------------------- prob = QuantumStateSmoothPulseProblem( sys, ψ_init, ψ_target, T, Δt; ipopt_options=IpoptOptions(print_level=1), verbose=false @@ -183,6 +178,7 @@ end @test final > initial # Multiple initial and target states + # ---------------------------------- ψ_inits = [[1.0, 0.0], [0.0, 1.0]] ψ_targets = [[0.0, 1.0], [1.0, 0.0]] prob = QuantumStateSmoothPulseProblem( From aad39159bbffcee7a0f877638738d2687df1bea8 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 21:55:46 -0500 Subject: [PATCH 29/35] state minimum time refactor for options --- src/constraints.jl | 31 ++---- src/objectives.jl | 4 +- .../quantum_state_minimum_time_problem.jl | 97 ++++++++++++------- .../quantum_state_smooth_pulse_problem.jl | 5 +- src/problems.jl | 18 +++- 5 files changed, 91 insertions(+), 64 deletions(-) diff --git a/src/constraints.jl b/src/constraints.jl index bc8c8af9..c6542e3e 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -130,7 +130,7 @@ function FinalFidelityConstraint(; zdim::Union{Int,Nothing}=nothing, T::Union{Int,Nothing}=nothing, subspace::Union{AbstractVector{<:Integer}, Nothing}=nothing, - hessian::Bool=true + eval_hessian::Bool=true ) @assert !isnothing(fidelity_function) "must provide a fidelity function" @assert !isnothing(value) "must provide a fidelity value" @@ -194,8 +194,8 @@ function FinalFidelityConstraint(; end end - if hessian - ∂²ℱ(x) = ForwardDiff.hessian(fid, x) + if eval_hessian + ∂²ℱ(x) = ForwardDiff.eval_hessian(fid, x) ∂²ℱ_structure = hessian_of_lagrangian_structure(∂²ℱ, statedim, 1) @@ -242,7 +242,7 @@ function FinalUnitaryFidelityConstraint( val::Float64, traj::NamedTrajectory; subspace::Union{AbstractVector{<:Integer}, Nothing}=nothing, - hessian::Bool=true + eval_hessian::Bool=true ) @assert statesymb ∈ traj.names return FinalFidelityConstraint(; @@ -254,7 +254,7 @@ function FinalUnitaryFidelityConstraint( zdim=traj.dim, T=traj.T, subspace=subspace, - hessian=hessian + eval_hessian=eval_hessian ) end @@ -268,7 +268,8 @@ is the NamedTrajectory symbol representing the unitary. function FinalQuantumStateFidelityConstraint( statesymb::Symbol, val::Float64, - traj::NamedTrajectory, + traj::NamedTrajectory; + kwargs... ) @assert statesymb ∈ traj.names return FinalFidelityConstraint(; @@ -278,26 +279,12 @@ function FinalQuantumStateFidelityConstraint( goal=traj.goal[statesymb], statedim=traj.dims[statesymb], zdim=traj.dim, - T=traj.T + T=traj.T, + kwargs... ) end - -# function FinalStateFidelityConstraint( -# val::Float64, -# statesymb::Symbol, -# statedim::Int; -# fidelity_function::Function=fidelity -# ) -# return FinalFidelityConstraint(; -# fidelity_function=fidelity_function, -# value=val, -# statesymb=statesymb, -# statedim=statedim -# ) -# end - """ ComplexModulusContraint() diff --git a/src/objectives.jl b/src/objectives.jl index f8be8db7..613439b3 100644 --- a/src/objectives.jl +++ b/src/objectives.jl @@ -910,10 +910,10 @@ function MinimumTimeObjective(; return Objective(L, ∇L, ∂²L, ∂²L_structure, Dict[params]) end -function MinimumTimeObjective(traj::NamedTrajectory; D=1.0) +function MinimumTimeObjective(traj::NamedTrajectory; D=1.0, kwargs...) @assert traj.timestep isa Symbol "trajectory does not have a dynamical timestep" Δt_indices = [index(t, traj.components[traj.timestep][1], traj.dim) for t = 1:traj.T] - return MinimumTimeObjective(; D=D, Δt_indices=Δt_indices) + return MinimumTimeObjective(; D=D, Δt_indices=Δt_indices, kwargs...) end @doc raw""" diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index 47510204..0ecbae9c 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -6,63 +6,94 @@ TODO: Add documentation function QuantumStateMinimumTimeProblem end function QuantumStateMinimumTimeProblem( - trajectory::NamedTrajectory, - system::QuantumSystem, - objective::Objective, + traj::NamedTrajectory, + sys::QuantumSystem, + obj::Objective, integrators::Vector{<:AbstractIntegrator}, constraints::Vector{<:AbstractConstraint}; state_symbol::Symbol=:ψ̃, + final_fidelity::Union{Real, Nothing}=nothing, D=1.0, - verbose::Bool=false, ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), kwargs... ) - @assert state_symbol ∈ trajectory.names + state_names = [name for name in traj.names if startswith(name, state_symbol)] + @assert length(state_names) ≥ 1 "No matching states found in trajectory" - objective += MinimumTimeObjective(trajectory; D=D) + obj += MinimumTimeObjective(traj; D=D, eval_hessian=piccolo_options.eval_hessian) - final_fidelity = fidelity(trajectory[end][state_symbol], trajectory.goal[state_symbol]) + # Default to average state fidelity + if isnothing(final_fidelity) + vals = [fidelity(traj[n][:, end], traj.goal[n]) for n ∈ state_names] + final_fidelity = sum(vals) / length(vals) + end - fidelity_constraint = FinalQuantumStateFidelityConstraint( - state_symbol, - final_fidelity, - trajectory - ) + for state_name in state_names + fidelity_constraint = FinalQuantumStateFidelityConstraint( + state_name, + final_fidelity, + traj, + eval_hessian=piccolo_options.eval_hessian + ) - push!(constraints, fidelity_constraint) + push!(constraints, fidelity_constraint) + end return QuantumControlProblem( - system, - trajectory, - objective, + sys, + traj, + obj, integrators; constraints=constraints, - verbose=verbose, ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end function QuantumStateMinimumTimeProblem( - data_path::String; + prob::QuantumControlProblem; + obj::Objective=get_objective(prob), + constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), + build_trajectory_constraints=false, kwargs... ) - data = load(data_path) - system = data["system"] - trajectory = data["trajectory"] - objective = Objective(data["params"][:objective_terms]) - integrators = data["params"][:dynamics] - constraints = AbstractConstraint[ - data["params"][:linear_constraints]..., - NonlinearConstraint.(data["params"][:nonlinear_constraints])... - ] + new_piccolo_options = deepcopy(prob.piccolo_options) + new_piccolo_options.build_trajectory_constraints = build_trajectory_constraints + return QuantumStateMinimumTimeProblem( - trajectory, - system, - objective, - integrators, - constraints; - build_trajectory_constraints=false, + prob.trajectory, + prob.system, + obj, + prob.integrators, + constraints, + ipopt_options=prob.ipopt_options, + piccolo_options=new_piccolo_options, kwargs... ) end + +# *************************************************************************** # + +@testitem "Test quantum state minimum time" begin + using NamedTrajectories + + # System + T = 50 + Δt = 0.2 + sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]]) + ψ_init = [1.0, 0.0] + ψ_target = [0.0, 1.0] + + prob = QuantumStateSmoothPulseProblem( + sys, ψ_init, ψ_target, T, Δt; + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) + ) + initial = sum(get_timesteps(prob.trajectory)) + mintime_prob = QuantumStateMinimumTimeProblem(prob) + solve!(mintime_prob, max_iter=20) + final = sum(get_timesteps(mintime_prob.trajectory)) + @test final < initial +end diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index cb5fe00f..bd2bfcdb 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -17,10 +17,7 @@ Create a quantum control problem for smooth pulse optimization of a quantum state trajectory. -# Keyword Arguments - -# TODO: clean up this whole constructor - +TODO: Document args """ function QuantumStateSmoothPulseProblem end diff --git a/src/problems.jl b/src/problems.jl index 522637ec..6acb8f9f 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -9,6 +9,7 @@ export update_trajectory! export get_traj_data export get_datavec export get_objective +export get_constraints using ..QuantumSystems using ..Integrators @@ -43,6 +44,7 @@ mutable struct QuantumControlProblem <: AbstractProblem trajectory::NamedTrajectory integrators::Union{Nothing,Vector{<:AbstractIntegrator}} ipopt_options::IpoptOptions + piccolo_options::PiccoloOptions params::Dict{Symbol, Any} end @@ -113,10 +115,8 @@ function QuantumControlProblem( variables = reshape(variables, traj.dim, traj.T) + # Container for saving constraints and objectives params = merge(kwargs, params) - - params[:ipopt_options] = ipopt_options - params[:piccolo_options] = piccolo_options params[:linear_constraints] = linear_constraints params[:nonlinear_constraints] = [ nl_constraint.params for nl_constraint ∈ nonlinear_constraints @@ -130,6 +130,7 @@ function QuantumControlProblem( traj, dynamics.integrators, ipopt_options, + piccolo_options, params ) end @@ -306,5 +307,16 @@ function get_objective(prob::QuantumControlProblem) return Objective(prob.params[:objective_terms]) end +""" + get_constraints(prob::QuantumControlProblem) + +Return the constraints of the `prob::QuantumControlProblem`. +""" +function get_constraints(prob::QuantumControlProblem) + return AbstractConstraint[ + prob.params[:linear_constraints]..., + NonlinearConstraint.(prob.params[:nonlinear_constraints])... + ] +end end From b8f62bde508cc1d1ab9e4c60fb4906c07dadd4be Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 22:15:53 -0500 Subject: [PATCH 30/35] refactor tests to testitems; move traj init to new file from rollout --- .../quantum_state_smooth_pulse_problem.jl | 6 +- src/rollouts.jl | 9 +- test/objectives_test.jl | 92 +++++++++---------- test/problems_test.jl | 9 +- test/quantum_system_templates_test.jl | 27 ------ test/rollouts_test.jl | 56 +---------- test/runtests.jl | 15 +-- test/trajectory_initialization_test.jl | 59 ++++++++++++ 8 files changed, 119 insertions(+), 154 deletions(-) delete mode 100644 test/quantum_system_templates_test.jl create mode 100644 test/trajectory_initialization_test.jl diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index bd2bfcdb..9362417f 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -167,7 +167,8 @@ end # -------------------------------- prob = QuantumStateSmoothPulseProblem( sys, ψ_init, ψ_target, T, Δt; - ipopt_options=IpoptOptions(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) initial = fidelity(prob) solve!(prob, max_iter=20) @@ -180,7 +181,8 @@ end ψ_targets = [[0.0, 1.0], [1.0, 0.0]] prob = QuantumStateSmoothPulseProblem( sys, ψ_inits, ψ_targets, T, Δt; - ipopt_options=IpoptOptions(print_level=1), verbose=false + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) initial = fidelity(prob) solve!(prob, max_iter=20) diff --git a/src/rollouts.jl b/src/rollouts.jl index 16df25c9..74dba666 100644 --- a/src/rollouts.jl +++ b/src/rollouts.jl @@ -12,6 +12,7 @@ using ..QuantumSystems using ..EmbeddedOperators using ..Integrators using ..Problems +using ..DirectSums using LinearAlgebra using NamedTrajectories @@ -89,10 +90,10 @@ function QuantumUtils.fidelity( kwargs... ) fids = [] - for name in trajectory.names - if startswith(string(name), string(state_symb)) - init = trajectory.initial[name] - goal = trajectory.goal[name] + for symb in trajectory.names + if startswith(symb, state_symb) + init = trajectory.initial[symb] + goal = trajectory.goal[symb] push!( fids, fidelity(init, goal, trajectory[control_symb], get_timesteps(trajectory), system; kwargs...) diff --git a/test/objectives_test.jl b/test/objectives_test.jl index 9aba2dcc..0d163618 100644 --- a/test/objectives_test.jl +++ b/test/objectives_test.jl @@ -2,61 +2,45 @@ Testing objective struct functionality """ -#TODO: Depends on TestUtils - - -# @testset "Objectives" begin - -# initializing test trajectory -# T = 10 -# H_drift = GATES[:Z] -# H_drives = [GATES[:X], GATES[:Y]] - -# system = QuantumSystem(H_drift, H_drives) - -# P = FourthOrderPade(system) - -# function f(zₜ, zₜ₊₁) -# ψ̃ₜ₊₁ = zₜ₊₁[Z.components.ψ̃] -# ψ̃ₜ = zₜ[Z.components.ψ̃] -# uₜ = zₜ[Z.components.u] -# δψ̃ = P(ψ̃ₜ₊₁, ψ̃ₜ, uₜ, Z.timestep) -# return δψ̃ -# end - -# dynamics = QuantumDynamics(f, Z) -# evaluator = PicoEvaluator(Z, J, dynamics, true) - -# @testset "Quantum State Objective" begin - -# Z = NamedTrajectory( -# (ψ̃ = randn(4, T), u = randn(2, T)), -# controls=:u, -# dt=0.1, -# goal=(ψ̃ = [1, 0, 0, 0],) -# ) +@testitem "Quantum State Objective" begin + using LinearAlgebra + using NamedTrajectories + using ForwardDiff + include("test_utils.jl") + + T = 10 + Z = NamedTrajectory( + (ψ̃ = randn(4, T), u = randn(2, T)), + controls=:u, + timestep=0.1, + goal=(ψ̃ = [1., 0., 0., 0.],) + ) + loss = :InfidelityLoss + Q = 100.0 -# loss = :InfidelityLoss -# Q = 100.0 + J = QuantumObjective(:ψ̃, Z, loss, Q) -# J = QuantumObjective(:ψ̃, Z, loss, Q) + L = Z⃗ -> J.L(Z⃗, Z) + ∇L = Z⃗ -> J.∇L(Z⃗, Z) + ∂²L = Z⃗ -> J.∂²L(Z⃗, Z) + ∂²L_structure = J.∂²L_structure(Z) -# L = Z⃗ -> J.L(Z⃗, Z) -# ∇L = Z⃗ -> J.∇L(Z⃗, Z) -# ∂²L = Z⃗ -> J.∂²L(Z⃗, Z) -# ∂²L_structure = J.∂²L_structure(Z) + # test objective function gradient + @test ForwardDiff.gradient(L, Z.datavec) ≈ ∇L(Z.datavec) -# # test objective function gradient -# @test ForwardDiff.gradient(L, Z.datavec) ≈ ∇L(Z.datavec) + # test objective function hessian + shape = (Z.dim * Z.T, Z.dim * Z.T) + @test ForwardDiff.hessian(L, Z.datavec) ≈ dense(∂²L(Z.datavec), ∂²L_structure, shape) +end -# # test objective function hessian -# shape = (Z.dim * Z.T, Z.dim * Z.T) -# @test ForwardDiff.hessian(L, Z.datavec) ≈ dense(∂²L(Z.datavec), ∂²L_structure, shape) -# end +@testitem "Quadratic Regularizer Objective" begin + using LinearAlgebra + using NamedTrajectories + using ForwardDiff + include("test_utils.jl") -@testset "Quadratic Regularizer Objective" begin T = 10 Z = NamedTrajectory( @@ -87,7 +71,12 @@ )) end -@testset "Quadratic Smoothness Regularizer Objective" begin +@testitem "Quadratic Smoothness Regularizer Objective" begin + using LinearAlgebra + using NamedTrajectories + using ForwardDiff + include("test_utils.jl") + T = 10 Z = NamedTrajectory( @@ -118,7 +107,12 @@ end )) end -@testset "Unitary Objective" begin +@testitem "Unitary Objective" begin + using LinearAlgebra + using NamedTrajectories + using ForwardDiff + include("test_utils.jl") + T = 10 U_init = GATES[:I] diff --git a/test/problems_test.jl b/test/problems_test.jl index ec13115c..6586b2e2 100644 --- a/test/problems_test.jl +++ b/test/problems_test.jl @@ -1,13 +1,14 @@ """ - Testing problems + Testing problem features + TODO: - test problem creation - - test problem solving + - test problem iterations - test problem saving - test problem loading """ -@testitem "Problems" begin +@testitem "System creation" begin # initializing test system T = 5 H_drift = GATES[:Z] @@ -16,7 +17,7 @@ system = QuantumSystem(H_drift, H_drives) - # test problem creation + # test system creation end diff --git a/test/quantum_system_templates_test.jl b/test/quantum_system_templates_test.jl deleted file mode 100644 index 4837eaa1..00000000 --- a/test/quantum_system_templates_test.jl +++ /dev/null @@ -1,27 +0,0 @@ -""" -Test: Quantum System Templates [RydbergChainSystem, TransmonSystem, QuantumOpticsSystem] -""" - -@testitem "Rydberg Chain System" begin - -end - -@testitem "Transmon System" begin - -end - -@testitem "Quantum Optics System" begin - using QuantumOpticsBase - N = rand(1:5); - T = ComplexF64; - b = FockBasis(N); - a = QuantumOpticsBase.create(T, b); - H = a + a'; - sys = QuantumOpticsSystem(H, [H, H]); - @test typeof(sys) == QuantumSystem - @test sys.constructor == QuantumOpticsSystem - @test sys.H_drift == H.data - - # creation with non-Hermitian operators - @test_broken QuantumOpticsSystem(a, [a]) -end \ No newline at end of file diff --git a/test/rollouts_test.jl b/test/rollouts_test.jl index 0cdc9d6a..053685d6 100644 --- a/test/rollouts_test.jl +++ b/test/rollouts_test.jl @@ -1,59 +1,7 @@ """ Testing rollouts -""" - -@testitem "Geodesic" begin - using LinearAlgebra - - ## Group 1: identity to X (π rotation) - - # Test π rotation - U_α = GATES[:I] - U_ω = GATES[:X] - Us, H = unitary_geodesic( - U_α, U_ω, range(0, 1, 4), return_generator=true - ) - - @test size(Us, 2) == 4 - @test Us[:, 1] ≈ operator_to_iso_vec(U_α) - @test Us[:, end] ≈ operator_to_iso_vec(U_ω) - @test H' - H ≈ zeros(2, 2) - @test norm(H) ≈ π - - # Test modified timesteps (10x) - Us10, H10 = unitary_geodesic( - U_α, U_ω, range(-5, 5, 4), return_generator=true - ) - @test size(Us10, 2) == 4 - @test Us10[:, 1] ≈ operator_to_iso_vec(U_α) - @test Us10[:, end] ≈ operator_to_iso_vec(U_ω) - @test H10' - H10 ≈ zeros(2, 2) - @test norm(H10) ≈ π/10 - - # Test wrapped call - Us_wrap, H_wrap = unitary_geodesic(U_ω, 10, return_generator=true) - @test Us_wrap[:, 1] ≈ operator_to_iso_vec(GATES[:I]) - @test Us_wrap[:, end] ≈ operator_to_iso_vec(U_ω) - rollout = [exp(-im * H_wrap * t) for t ∈ range(0, 1, 10)] - Us_test = stack(operator_to_iso_vec.(rollout), dims=2) - @test isapprox(Us_wrap, Us_test) - - - ## Group 2: √X to X (π/2 rotation) - - # Test geodesic not at identity - U₀ = sqrt(GATES[:X]) - U₁ = GATES[:X] - Us, H = unitary_geodesic(U₀, U₁, 10, return_generator=true) - @test Us[:, 1] ≈ operator_to_iso_vec(U₀) - @test Us[:, end] ≈ operator_to_iso_vec(U_ω) + TODO: Add tests for rollouts +""" - rollout = [exp(-im * H * t) * U₀ for t ∈ range(0, 1, 10)] - Us_test = stack(operator_to_iso_vec.(rollout), dims=2) - @test isapprox(Us, Us_test) - Us_wrap = unitary_geodesic(U_ω, 4) - @test Us_wrap[:, 1] ≈ operator_to_iso_vec(GATES[:I]) - @test Us_wrap[:, end] ≈ operator_to_iso_vec(U_ω) -end diff --git a/test/runtests.jl b/test/runtests.jl index 71aad2f1..a9e6f4ae 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,18 +1,5 @@ using Test using TestItemRunner -using LinearAlgebra -using ForwardDiff -using SparseArrays -using Random; Random.seed!(1234) - -using QuantumCollocation -using NamedTrajectories - - - - -include("test_utils.jl") - -# Run testitem +# Run all testitem tests in package @run_package_tests diff --git a/test/trajectory_initialization_test.jl b/test/trajectory_initialization_test.jl new file mode 100644 index 00000000..1f7c604d --- /dev/null +++ b/test/trajectory_initialization_test.jl @@ -0,0 +1,59 @@ +""" + Testing trajectory initialization +""" + +@testitem "Geodesic" begin + using LinearAlgebra + + ## Group 1: identity to X (π rotation) + + # Test π rotation + U_α = GATES[:I] + U_ω = GATES[:X] + Us, H = unitary_geodesic( + U_α, U_ω, range(0, 1, 4), return_generator=true + ) + + @test size(Us, 2) == 4 + @test Us[:, 1] ≈ operator_to_iso_vec(U_α) + @test Us[:, end] ≈ operator_to_iso_vec(U_ω) + @test H' - H ≈ zeros(2, 2) + @test norm(H) ≈ π + + # Test modified timesteps (10x) + Us10, H10 = unitary_geodesic( + U_α, U_ω, range(-5, 5, 4), return_generator=true + ) + + @test size(Us10, 2) == 4 + @test Us10[:, 1] ≈ operator_to_iso_vec(U_α) + @test Us10[:, end] ≈ operator_to_iso_vec(U_ω) + @test H10' - H10 ≈ zeros(2, 2) + @test norm(H10) ≈ π/10 + + # Test wrapped call + Us_wrap, H_wrap = unitary_geodesic(U_ω, 10, return_generator=true) + @test Us_wrap[:, 1] ≈ operator_to_iso_vec(GATES[:I]) + @test Us_wrap[:, end] ≈ operator_to_iso_vec(U_ω) + rotation = [exp(-im * H_wrap * t) for t ∈ range(0, 1, 10)] + Us_test = stack(operator_to_iso_vec.(rotation), dims=2) + @test isapprox(Us_wrap, Us_test) + + + ## Group 2: √X to X (π/2 rotation) + + # Test geodesic not at identity + U₀ = sqrt(GATES[:X]) + U₁ = GATES[:X] + Us, H = unitary_geodesic(U₀, U₁, 10, return_generator=true) + @test Us[:, 1] ≈ operator_to_iso_vec(U₀) + @test Us[:, end] ≈ operator_to_iso_vec(U_ω) + + rotation = [exp(-im * H * t) * U₀ for t ∈ range(0, 1, 10)] + Us_test = stack(operator_to_iso_vec.(rotation), dims=2) + @test isapprox(Us, Us_test) + Us_wrap = unitary_geodesic(U_ω, 4) + @test Us_wrap[:, 1] ≈ operator_to_iso_vec(GATES[:I]) + @test Us_wrap[:, end] ≈ operator_to_iso_vec(U_ω) + +end \ No newline at end of file From e0c3800522593c84fb2401b0432a400a360887a5 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 22:37:27 -0500 Subject: [PATCH 31/35] unitary min time with piccolo params --- .../quantum_state_minimum_time_problem.jl | 6 +- .../unitary_minimum_time_problem.jl | 91 +++++++------------ src/problems.jl | 12 ++- 3 files changed, 44 insertions(+), 65 deletions(-) diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index 0ecbae9c..d542b972 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -63,12 +63,12 @@ function QuantumStateMinimumTimeProblem( new_piccolo_options.build_trajectory_constraints = build_trajectory_constraints return QuantumStateMinimumTimeProblem( - prob.trajectory, + copy(prob.trajectory), prob.system, obj, prob.integrators, - constraints, - ipopt_options=prob.ipopt_options, + constraints; + ipopt_options=deepcopy(prob.ipopt_options), piccolo_options=new_piccolo_options, kwargs... ) diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index 9b912782..16568a6c 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -41,8 +41,8 @@ J(\vec{\tilde{U}}, a, \dot{a}, \ddot{a}) + D \sum_t \Delta t_t \\ - `unitary_symbol::Symbol=:Ũ⃗`: The symbol for the unitary control. - `final_fidelity::Float64=0.99`: The final fidelity. - `D=1.0`: The weight for the minimum-time objective. -- `verbose::Bool=false`: Whether to print additional information. - `ipopt_options::IpoptOptions=IpoptOptions()`: The options for the Ipopt solver. +- `piccolo_options::PiccoloOptions=PiccoloOptions()`: The options for the Piccolo solver. - `kwargs...`: Additional keyword arguments to pass to `QuantumControlProblem`. """ function UnitaryMinimumTimeProblem end @@ -56,23 +56,24 @@ function UnitaryMinimumTimeProblem( unitary_symbol::Symbol=:Ũ⃗, final_fidelity::Float64=0.99, D=1.0, - verbose::Bool=false, ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), subspace=nothing, kwargs... ) @assert unitary_symbol ∈ trajectory.names - objective += MinimumTimeObjective(trajectory; D=D) + objective += MinimumTimeObjective(trajectory; D=D, eval_hessian=piccolo_options.eval_hessian) fidelity_constraint = FinalUnitaryFidelityConstraint( unitary_symbol, final_fidelity, trajectory; - subspace=subspace + subspace=subspace, + eval_hessian=piccolo_options.eval_hessian ) - constraints = AbstractConstraint[constraints..., fidelity_constraint] + constraints = push!(constraints, fidelity_constraint) return QuantumControlProblem( system, @@ -80,56 +81,30 @@ function UnitaryMinimumTimeProblem( objective, integrators; constraints=constraints, - verbose=verbose, ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end function UnitaryMinimumTimeProblem( prob::QuantumControlProblem; + objective::Objective=get_objective(prob), + constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), + build_trajectory_constraints=false, kwargs... ) - params = deepcopy(prob.params) - trajectory = copy(prob.trajectory) - system = prob.system - objective = Objective(params[:objective_terms]) - integrators = prob.integrators - constraints = [ - params[:linear_constraints]..., - NonlinearConstraint.(params[:nonlinear_constraints])... - ] + new_piccolo_options = deepcopy(prob.piccolo_options) + new_piccolo_options.build_trajectory_constraints = build_trajectory_constraints + return UnitaryMinimumTimeProblem( - trajectory, - system, - objective, - integrators, - constraints; - build_trajectory_constraints=false, - kwargs... - ) -end - -function UnitaryMinimumTimeProblem( - data_path::String; - kwargs... -) - data = load(data_path) - system = data["system"] - trajectory = data["trajectory"] - objective = Objective(data["params"][:objective_terms]) - integrators = data["integrators"] - constraints = AbstractConstraint[ - data["params"][:linear_constraints]..., - NonlinearConstraint.(data["params"][:nonlinear_constraints])... - ] - return UnitaryMinimumTimeProblem( - trajectory, - system, + copy(prob.trajectory), + prob.system, objective, - integrators, + prob.integrators, constraints; - build_trajectory_constraints=false, + ipopt_options=deepcopy(prob.ipopt_options), + piccolo_options=new_piccolo_options, kwargs... ) end @@ -137,6 +112,8 @@ end # *************************************************************************** # @testitem "Minimum time Hadamard gate" begin + using NamedTrajectories + H_drift = GATES[:Z] H_drives = [GATES[:X], GATES[:Y]] U_goal = GATES[:H] @@ -145,24 +122,20 @@ end prob = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=1) - ) - - solve!(prob, max_iter=100) - - @test unitary_fidelity(prob) > 0.99 - - final_fidelity = 0.99 - - mintime_prob = UnitaryMinimumTimeProblem( - prob, - final_fidelity=final_fidelity, - ipopt_options=IpoptOptions(print_level=1) + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) - solve!(mintime_prob; max_iter=100) + before = unitary_fidelity(prob) + solve!(prob, max_iter=20) + after = unitary_fidelity(prob) + @test after > before - @test unitary_fidelity(mintime_prob) > final_fidelity + # Soft fidelity constraint + final_fidelity = minimum([0.99, after]) + mintime_prob = UnitaryMinimumTimeProblem(prob, final_fidelity=final_fidelity) + solve!(mintime_prob; max_iter=40) - @test sum(mintime_prob.trajectory[:Δt]) < sum(prob.trajectory[:Δt]) + @test unitary_fidelity(mintime_prob) ≥ after + @test sum(get_timesteps(mintime_prob.trajectory)) < sum(get_timesteps(prob.trajectory)) end diff --git a/src/problems.jl b/src/problems.jl index 6acb8f9f..1e4f1139 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -302,20 +302,26 @@ end get_objective(prob::QuantumControlProblem) Return the objective function of the `prob::QuantumControlProblem`. + +TODO: Is deepcopy necessary? """ function get_objective(prob::QuantumControlProblem) - return Objective(prob.params[:objective_terms]) + params = deepcopy(prob.params) + return Objective(params[:objective_terms]) end """ get_constraints(prob::QuantumControlProblem) Return the constraints of the `prob::QuantumControlProblem`. + +TODO: Is deepcopy necessary? """ function get_constraints(prob::QuantumControlProblem) + params = deepcopy(prob.params) return AbstractConstraint[ - prob.params[:linear_constraints]..., - NonlinearConstraint.(prob.params[:nonlinear_constraints])... + params[:linear_constraints]..., + NonlinearConstraint.(params[:nonlinear_constraints])... ] end From 5d3e61833c265f928c71863303a4c6f741cf78cc Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 23:03:58 -0500 Subject: [PATCH 32/35] unitary robustness update to piccolo options --- src/constraints.jl | 2 +- .../quantum_state_minimum_time_problem.jl | 12 +- .../unitary_minimum_time_problem.jl | 26 ++-- .../unitary_robustness_problem.jl | 129 +++++++----------- 4 files changed, 75 insertions(+), 94 deletions(-) diff --git a/src/constraints.jl b/src/constraints.jl index c6542e3e..ff1981e3 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -195,7 +195,7 @@ function FinalFidelityConstraint(; end if eval_hessian - ∂²ℱ(x) = ForwardDiff.eval_hessian(fid, x) + ∂²ℱ(x) = ForwardDiff.hessian(fid, x) ∂²ℱ_structure = hessian_of_lagrangian_structure(∂²ℱ, statedim, 1) diff --git a/src/problem_templates/quantum_state_minimum_time_problem.jl b/src/problem_templates/quantum_state_minimum_time_problem.jl index d542b972..6e9a22a1 100644 --- a/src/problem_templates/quantum_state_minimum_time_problem.jl +++ b/src/problem_templates/quantum_state_minimum_time_problem.jl @@ -56,11 +56,12 @@ function QuantumStateMinimumTimeProblem( prob::QuantumControlProblem; obj::Objective=get_objective(prob), constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), + ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), + piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), build_trajectory_constraints=false, kwargs... ) - new_piccolo_options = deepcopy(prob.piccolo_options) - new_piccolo_options.build_trajectory_constraints = build_trajectory_constraints + piccolo_options.build_trajectory_constraints = build_trajectory_constraints return QuantumStateMinimumTimeProblem( copy(prob.trajectory), @@ -68,8 +69,8 @@ function QuantumStateMinimumTimeProblem( obj, prob.integrators, constraints; - ipopt_options=deepcopy(prob.ipopt_options), - piccolo_options=new_piccolo_options, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end @@ -96,4 +97,7 @@ end solve!(mintime_prob, max_iter=20) final = sum(get_timesteps(mintime_prob.trajectory)) @test final < initial + + # Test with final fidelity + QuantumStateMinimumTimeProblem(prob, final_fidelity=0.99) end diff --git a/src/problem_templates/unitary_minimum_time_problem.jl b/src/problem_templates/unitary_minimum_time_problem.jl index 16568a6c..acb62d83 100644 --- a/src/problem_templates/unitary_minimum_time_problem.jl +++ b/src/problem_templates/unitary_minimum_time_problem.jl @@ -13,11 +13,6 @@ kwargs... ) - UnitaryMinimumTimeProblem( - data_path::String; - kwargs... - ) - Create a minimum-time problem for unitary control. ```math @@ -54,7 +49,7 @@ function UnitaryMinimumTimeProblem( integrators::Vector{<:AbstractIntegrator}, constraints::Vector{<:AbstractConstraint}; unitary_symbol::Symbol=:Ũ⃗, - final_fidelity::Float64=0.99, + final_fidelity::Union{Real, Nothing}=nothing, D=1.0, ipopt_options::IpoptOptions=IpoptOptions(), piccolo_options::PiccoloOptions=PiccoloOptions(), @@ -63,6 +58,12 @@ function UnitaryMinimumTimeProblem( ) @assert unitary_symbol ∈ trajectory.names + if isnothing(final_fidelity) + final_fidelity = unitary_fidelity( + trajectory[unitary_symbol][:, end], trajectory.goal[unitary_symbol] + ) + end + objective += MinimumTimeObjective(trajectory; D=D, eval_hessian=piccolo_options.eval_hessian) fidelity_constraint = FinalUnitaryFidelityConstraint( @@ -91,11 +92,12 @@ function UnitaryMinimumTimeProblem( prob::QuantumControlProblem; objective::Objective=get_objective(prob), constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), + ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), + piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), build_trajectory_constraints=false, kwargs... ) - new_piccolo_options = deepcopy(prob.piccolo_options) - new_piccolo_options.build_trajectory_constraints = build_trajectory_constraints + piccolo_options.build_trajectory_constraints = build_trajectory_constraints return UnitaryMinimumTimeProblem( copy(prob.trajectory), @@ -103,8 +105,8 @@ function UnitaryMinimumTimeProblem( objective, prob.integrators, constraints; - ipopt_options=deepcopy(prob.ipopt_options), - piccolo_options=new_piccolo_options, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end @@ -138,4 +140,8 @@ end @test unitary_fidelity(mintime_prob) ≥ after @test sum(get_timesteps(mintime_prob.trajectory)) < sum(get_timesteps(prob.trajectory)) + + # Test without final fidelity + UnitaryMinimumTimeProblem(prob) + end diff --git a/src/problem_templates/unitary_robustness_problem.jl b/src/problem_templates/unitary_robustness_problem.jl index 0d75368a..ef19822d 100644 --- a/src/problem_templates/unitary_robustness_problem.jl +++ b/src/problem_templates/unitary_robustness_problem.jl @@ -1,11 +1,11 @@ @doc raw""" - UnitaryRobustnessProblem(Hₑ, trajectory, system, objective, integrators, constraints; + UnitaryRobustnessProblem( + H_error, trajectory, system, objective, integrators, constraints; unitary_symbol=:Ũ⃗, - final_fidelity=unitary_fidelity(trajectory[end][unitary_symbol], trajectory.goal[unitary_symbol]), + final_fidelity=nothing, subspace=nothing, - eval_hessian=false, - verbose=false, ipopt_options=IpoptOptions(), + piccolo_options=PiccoloOptions(), kwargs... ) @@ -17,30 +17,31 @@ function UnitaryRobustnessProblem end function UnitaryRobustnessProblem( - Hₑ::AbstractMatrix{<:Number}, + H_error::AbstractMatrix{<:Number}, trajectory::NamedTrajectory, system::QuantumSystem, objective::Objective, integrators::Vector{<:AbstractIntegrator}, constraints::Vector{<:AbstractConstraint}; unitary_symbol::Symbol=:Ũ⃗, - final_fidelity::Float64=unitary_fidelity(trajectory[end][unitary_symbol], trajectory.goal[unitary_symbol]), - subspace::AbstractVector{<:Integer}=1:size(Hₑ, 1), - eval_hessian::Bool=false, - verbose::Bool=false, + final_fidelity::Union{Real, Nothing}=nothing, + subspace::AbstractVector{<:Integer}=1:size(H_error, 1), ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), kwargs... ) @assert unitary_symbol ∈ trajectory.names - if !eval_hessian - ipopt_options.hessian_approximation = "limited-memory" + if isnothing(final_fidelity) + final_fidelity = unitary_fidelity( + trajectory[unitary_symbol][:, end], trajectory.goal[unitary_symbol] + ) end objective += InfidelityRobustnessObjective( - Hₑ, + H_error, trajectory, - eval_hessian=eval_hessian, + eval_hessian=piccolo_options.eval_hessian, subspace=subspace ) @@ -50,8 +51,7 @@ function UnitaryRobustnessProblem( trajectory; subspace=subspace ) - - constraints = AbstractConstraint[constraints..., fidelity_constraint] + push!(constraints, fidelity_constraint) return QuantumControlProblem( system, @@ -59,36 +59,33 @@ function UnitaryRobustnessProblem( objective, integrators; constraints=constraints, - verbose=verbose, ipopt_options=ipopt_options, - eval_hessian=eval_hessian, + piccolo_options=piccolo_options, kwargs... ) end function UnitaryRobustnessProblem( - Hₑ::AbstractMatrix{<:Number}, + H_error::AbstractMatrix{<:Number}, prob::QuantumControlProblem; + objective::Objective=get_objective(prob), + constraints::AbstractVector{<:AbstractConstraint}=get_constraints(prob), + ipopt_options::IpoptOptions=deepcopy(prob.ipopt_options), + piccolo_options::PiccoloOptions=deepcopy(prob.piccolo_options), + build_trajectory_constraints=false, kwargs... ) - params = deepcopy(prob.params) - trajectory = copy(prob.trajectory) - system = prob.system - objective = Objective(params[:objective_terms]) - integrators = prob.integrators - constraints = [ - params[:linear_constraints]..., - NonlinearConstraint.(params[:nonlinear_constraints])... - ] + piccolo_options.build_trajectory_constraints = build_trajectory_constraints return UnitaryRobustnessProblem( - Hₑ, - trajectory, - system, + H_error, + copy(prob.trajectory), + prob.system, objective, - integrators, + prob.integrators, constraints; - build_trajectory_constraints=false, + ipopt_options=ipopt_options, + piccolo_options=piccolo_options, kwargs... ) end @@ -108,66 +105,40 @@ end subspace = U_goal.subspace_indices T = 51 Δt = 0.2 - probs = Dict() - # -------------------------------------------- - # 1. test UnitarySmoothPulseProblem with subspace - # - rely on linear interpolation of unitary - # -------------------------------------------- - probs["transmon"] = UnitarySmoothPulseProblem( + # test initial problem + # --------------------- + prob = UnitarySmoothPulseProblem( sys, U_goal, T, Δt, - geodesic=false, + geodesic=true, verbose=false, - ipopt_options=IpoptOptions(print_level=1) + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false) ) - solve!(probs["transmon"], max_iter=200) + before = unitary_fidelity(prob, subspace=subspace) + solve!(prob, max_iter=20) + after = unitary_fidelity(prob, subspace=subspace) # Subspace gate success - @test unitary_fidelity(probs["transmon"], subspace=subspace) > 0.99 + @test after > before - # -------------------------------------------- - # 2. test UnitaryRobustnessProblem from previous problem - # -------------------------------------------- - probs["robust"] = UnitaryRobustnessProblem( - H_error, probs["transmon"], - final_fidelity=0.99, + # test robustness from previous problem + # -------------------------------------- + final_fidelity = 0.99 + rob_prob = UnitaryRobustnessProblem( + H_error, prob, + final_fidelity=final_fidelity, subspace=subspace, - verbose=false, - ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=1.0, print_level=1) + ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=10.0, print_level=1), ) - solve!(probs["robust"], max_iter=200) + solve!(rob_prob, max_iter=50) - eval_loss(problem, Loss) = Loss(vec(problem.trajectory.data), problem.trajectory) - loss = InfidelityRobustnessObjective(H_error, probs["transmon"].trajectory).L + loss(Z⃗) = InfidelityRobustnessObjective(H_error, prob.trajectory).L(Z⃗, prob.trajectory) # Robustness improvement over default - @test eval_loss(probs["robust"], loss) < eval_loss(probs["transmon"], loss) - - # Fidelity constraint approximately satisfied - @test isapprox(unitary_fidelity(probs["robust"]; subspace=subspace), 0.99, atol=0.025) - - # -------------------------------------------- - # 3. test UnitaryRobustnessProblem from default struct - # -------------------------------------------- - params = deepcopy(probs["transmon"].params) - trajectory = copy(probs["transmon"].trajectory) - system = probs["transmon"].system - objective = QuadraticRegularizer(:dda, trajectory, 1e-4) - integrators = probs["transmon"].integrators - constraints = AbstractConstraint[] - - probs["unconstrained"] = UnitaryRobustnessProblem( - H_error, trajectory, system, objective, integrators, constraints, - final_fidelity=0.99, - subspace=subspace, - ipopt_options=IpoptOptions(recalc_y="yes", recalc_y_feas_tol=1e-1, print_level=4) - ) - solve!(probs["unconstrained"]; max_iter=100) - - # Additonal robustness improvement after relaxed objective - @test eval_loss(probs["unconstrained"], loss) < eval_loss(probs["transmon"], loss) + @test loss(rob_prob.trajectory.datavec) < loss(prob.trajectory.datavec) # Fidelity constraint approximately satisfied - @test isapprox(unitary_fidelity(probs["unconstrained"]; subspace=subspace), 0.99, atol=0.025) + @test isapprox(unitary_fidelity(rob_prob; subspace=subspace), 0.99, atol=0.025) end From 21d246fcd461b4a7694944dcaabd9758da3eba7a Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 23:32:16 -0500 Subject: [PATCH 33/35] adapt direct sum to piccolo_options --- .../quantum_state_smooth_pulse_problem.jl | 2 +- .../unitary_direct_sum_problem.jl | 132 +++--------------- src/trajectory_initialization.jl | 16 +-- test/trajectory_initialization_test.jl | 14 ++ 4 files changed, 45 insertions(+), 119 deletions(-) diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index 9362417f..cc1fcf36 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -69,7 +69,7 @@ function QuantumStateSmoothPulseProblem( if !isnothing(init_trajectory) traj = init_trajectory else - traj = initialize_state_trajectory( + traj = initialize_quantum_state_trajectory( ψ̃_goals, ψ̃_inits, T, diff --git a/src/problem_templates/unitary_direct_sum_problem.jl b/src/problem_templates/unitary_direct_sum_problem.jl index 6971ad97..bcb89af9 100644 --- a/src/problem_templates/unitary_direct_sum_problem.jl +++ b/src/problem_templates/unitary_direct_sum_problem.jl @@ -60,39 +60,15 @@ function UnitaryDirectSumProblem( drive_reset_ratio::Float64=0.50, fidelity_cost::Bool=false, subspace::Union{AbstractVector{<:Integer}, Nothing}=nothing, - max_iter::Int=1000, - linear_solver::String="mumps", - verbose::Bool=false, - hessian_approximation=true, - jacobian_structure=true, - blas_multithreading=true, - ipopt_options::IpoptOptions=IpoptOptions(), + ipopt_options::IpoptOptions=deepcopy(probs[1].ipopt_options), + piccolo_options::PiccoloOptions=deepcopy(probs[1].piccolo_options), kwargs... ) N = length(probs) - if length(prob_labels) != N - throw(ArgumentError("Length of prob_labels must match length of probs")) - end - - if N < 2 - throw(ArgumentError("At least two problems are required")) - end - - if drive_reset_ratio < 0 || drive_reset_ratio > 1 - throw(ArgumentError("drive_reset_σ must be in [0, 1]")) - end - - if !isempty(intersect(keys(boundary_values), prob_labels)) - throw(ArgumentError("Boundary value keys cannot be in prob_labels")) - end - - if hessian_approximation - ipopt_options.hessian_approximation = "limited-memory" - end - - if !blas_multithreading - BLAS.set_num_threads(1) - end + @assert length(prob_labels) == N "prob_labels must match length of probs" + @assert N ≥ 2 "At least two problems are required" + @assert 0 ≤ drive_reset_ratio ≤ 1 "drive_reset_ratio must be in [0, 1]" + @assert isempty(intersect(keys(boundary_values), prob_labels)) "Boundary value keys cannot be in prob_labels" # Default chain graph and boundary boundary = Tuple{Symbol, Array}[] @@ -127,13 +103,11 @@ function UnitaryDirectSumProblem( # add noise to control data (avoid restoration) if drive_reset_ratio > 0 - σs = repeat([drive_derivative_σ], 2) for ℓ in prob_labels - a_symb = add_suffix(:a, ℓ) - da_symb = add_suffix(:da, ℓ) - dda_symb = add_suffix(:dda, ℓ) - a_bounds_lower, a_bounds_upper = traj.bounds[a_symb] - a, da, dda = randomly_fill_drives(traj.T, a_bounds_lower, a_bounds_upper, σs) + a_symb, da_symb, dda_symb = add_suffix(:a, ℓ), add_suffix(:da, ℓ), add_suffix(:dda, ℓ) + a, da, dda = TrajectoryInitialization.initialize_controls( + length(traj.components[a_symb]), traj.T, traj.bounds[a_symb], drive_derivative_σ + ) update!(traj, a_symb, (1 - drive_reset_ratio) * traj[a_symb] + drive_reset_ratio * a) update!(traj, da_symb, (1 - drive_reset_ratio) * traj[da_symb] + drive_reset_ratio * da) update!(traj, dda_symb, (1 - drive_reset_ratio) * traj[dda_symb] + drive_reset_ratio * dda) @@ -149,7 +123,7 @@ function UnitaryDirectSumProblem( system = direct_sum([add_suffix(p.system, ℓ) for (p, ℓ) ∈ zip(probs, prob_labels)]) # Rebuild trajectory constraints - build_trajectory_constraints = true + piccolo_options.build_trajectory_constraints = true constraints = AbstractConstraint[] # Add goal constraints for each problem @@ -160,7 +134,7 @@ function UnitaryDirectSumProblem( final_fidelity, traj, subspace=subspace, - hessian=!hessian_approximation + eval_hessian=piccolo_options.eval_hessian ) push!(constraints, fidelity_constraint) end @@ -168,8 +142,8 @@ function UnitaryDirectSumProblem( # Build the objective function J = PairwiseQuadraticRegularizer(traj, Q, graph) - for (symb, s₀) ∈ boundary - J += QuadraticRegularizer(symb, traj, R_b; baseline=s₀) + for (symb, val) ∈ boundary + J += QuadraticRegularizer(symb, traj, R_b; baseline=val) end for ℓ ∈ prob_labels @@ -185,7 +159,7 @@ function UnitaryDirectSumProblem( J += UnitaryInfidelityObjective( add_suffix(:Ũ⃗, ℓ), traj, Q_fid, subspace=subspace, - eval_hessian=!hessian_approximation + eval_hessian=piccolo_options.eval_hessian ) end end @@ -196,78 +170,14 @@ function UnitaryDirectSumProblem( J, integrators; constraints=constraints, - max_iter=max_iter, - linear_solver=linear_solver, - verbose=verbose, ipopt_options=ipopt_options, - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - eval_hessian=!hessian_approximation, - build_trajectory_constraints=build_trajectory_constraints + piccolo_options=piccolo_options, ) end -# *************************************************************************** # - -function randomly_fill_drives( - T::Int, - drive_bounds_lower::AbstractVector{<:Real}, - drive_bounds_upper::AbstractVector{<:Real}, - drive_derivatives_σ::AbstractVector{<:Real} -) - data = Matrix{Float64}[] - - # Drive - n_drives = length(drive_bounds_lower) - a_dists = [Uniform(drive_bounds_lower[i], drive_bounds_upper[i]) for i = 1:n_drives] - push!( - data, - hcat([ - zeros(n_drives), - vcat([rand(a_dists[i], 1, T - 2) for i = 1:n_drives]...), - zeros(n_drives) - ]...) - ) - - # Drive derivatives - for σ ∈ drive_derivatives_σ - push!( - data, - randn(n_drives, T) * σ - ) - end - return data -end - -function randomly_fill_drives( - T::Int, - drive_bounds::AbstractVector{<:Real}, - drive_derivatives_σ::AbstractVector{<:Real} -) - return randomly_fill_drives( - T, - -drive_bounds, - drive_bounds, - drive_derivatives_σ - ) -end # *************************************************************************** # -@testitem "Random drive initialization" begin - T = 10 - n_drives = 2 - drive_bounds = [1.0, 2.0] - drive_derivatives_σ = repeat([0.01], 2) - - a, da, dda = ProblemTemplates.randomly_fill_drives(T, drive_bounds, drive_derivatives_σ) - - @test size(a) == (n_drives, T) - @test size(da) == (n_drives, T) - @test size(dda) == (n_drives, T) - @test all([-drive_bounds[i] < minimum(a[i, :]) < drive_bounds[i] for i in 1:n_drives]) -end - @testitem "Construct direct sum problem" begin sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:X], GATES[:Y]]) U_goal1 = GATES[:X] @@ -275,12 +185,14 @@ end U_goal2 = U_ε'GATES[:X]*U_ε T = 50 Δt = 0.2 - ops = Options(print_level=1) - prob1 = UnitarySmoothPulseProblem(sys, U_goal1, T, Δt, free_time=false, verbose=false, ipopt_options=ops) - prob2 = UnitarySmoothPulseProblem(sys, U_goal2, T, Δt, free_time=false, verbose=false, ipopt_options=ops) + ops = IpoptOptions(print_level=1) + pops = PiccoloOptions(verbose=false, free_time=false) + + prob1 = UnitarySmoothPulseProblem(sys, U_goal1, T, Δt, ipopt_options=ops, piccolo_options=pops) + prob2 = UnitarySmoothPulseProblem(sys, U_goal2, T, Δt, ipopt_options=ops, piccolo_options=pops) # Test default - direct_sum_prob1 = UnitaryDirectSumProblem([prob1, prob2], 0.99, verbose=false, ipopt_options=ops) + direct_sum_prob1 = UnitaryDirectSumProblem([prob1, prob2], 0.99) state_names = vcat( add_suffix(prob1.trajectory.state_names, "1")..., add_suffix(prob2.trajectory.state_names, "2")... diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 57e866e2..0fe8b7b8 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -4,7 +4,7 @@ export unitary_geodesic export linear_interpolation export unitary_linear_interpolation export initialize_unitary_trajectory -export initialize_state_trajectory +export initialize_quantum_state_trajectory using NamedTrajectories using LinearAlgebra @@ -111,15 +111,15 @@ function unitary_geodesic( end function unitary_geodesic( - U₀::AbstractMatrix{<:Number}, - U₁::AbstractMatrix{<:Number}, + U_init::AbstractMatrix{<:Number}, + U_goal::AbstractMatrix{<:Number}, timesteps::AbstractVector{<:Number}; return_generator=false ) """ Compute the effective generator of the geodesic connecting U₀ and U₁. - U₁ = exp(-im * H * T) U₀ - log(U₁ * U₀') = -im * H * T + U_goal = exp(-im * H * T) U_init + log(U_goal * U_init') = -im * H * T Allow for the possibiltiy of unequal timesteps and ranges outside [0,1]. @@ -128,9 +128,9 @@ function unitary_geodesic( """ t₀ = timesteps[1] T = timesteps[end] - t₀ - H = im * log(U₁ * U₀') / T + H = im * log(U_goal * U_init') / T # -im prefactor is not included in H - U_geo = [exp(-im * H * (t - t₀)) * U₀ for t ∈ timesteps] + U_geo = [exp(-im * H * (t - t₀)) * U_init for t ∈ timesteps] Ũ⃗_geo = stack(operator_to_iso_vec.(U_geo), dims=2) if return_generator return Ũ⃗_geo, H @@ -295,7 +295,7 @@ function initialize_unitary_trajectory( ) end -function initialize_state_trajectory( +function initialize_quantum_state_trajectory( ψ̃_goals::AbstractVector{<:AbstractVector{<:Real}}, ψ̃_inits::AbstractVector{<:AbstractVector{<:Real}}, T::Int, diff --git a/test/trajectory_initialization_test.jl b/test/trajectory_initialization_test.jl index 1f7c604d..63aa48f8 100644 --- a/test/trajectory_initialization_test.jl +++ b/test/trajectory_initialization_test.jl @@ -2,6 +2,20 @@ Testing trajectory initialization """ +@testitem "Random drive initialization" begin + T = 10 + n_drives = 2 + drive_bounds = [1.0, 2.0] + drive_derivative_σ = 0.01 + + a, da, dda = TrajectoryInitialization.initialize_controls(n_drives, T, drive_bounds, drive_derivative_σ) + + @test size(a) == (n_drives, T) + @test size(da) == (n_drives, T) + @test size(dda) == (n_drives, T) + @test all([-drive_bounds[i] < minimum(a[i, :]) < drive_bounds[i] for i in 1:n_drives]) +end + @testitem "Geodesic" begin using LinearAlgebra From a6d349ad40812914e3b7a95b52b95a016b8f27fc Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 23:48:02 -0500 Subject: [PATCH 34/35] unitary sampling problem uses piccolo options --- .../unitary_sampling_problem.jl | 114 ++++++------------ 1 file changed, 40 insertions(+), 74 deletions(-) diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index bdbf38f6..4a9ede7e 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -57,8 +57,8 @@ function UnitarySamplingProblem( Δt::Union{Float64, Vector{Float64}}; system_labels=string.(1:length(systems)), system_weights=fill(1.0, length(systems)), - free_time=true, init_trajectory::Union{NamedTrajectory, Nothing}=nothing, + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], a_bound::Float64=1.0, a_bounds=fill(a_bound, length(systems[1].G_drives)), a_guess::Union{Matrix{Float64}, Nothing}=nothing, @@ -74,34 +74,14 @@ function UnitarySamplingProblem( R_dda::Union{Float64, Vector{Float64}}=R, leakage_suppression=false, R_leakage=1e-1, - max_iter::Int=1000, - linear_solver::String="mumps", - ipopt_options::IpoptOptions=IpoptOptions(), - constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], - timesteps_all_equal::Bool=true, - verbose::Bool=false, - integrator::Symbol=:pade, - rollout_integrator=exp, - bound_unitary=integrator == :exponential, control_norm_constraint=false, control_norm_constraint_components=nothing, control_norm_R=nothing, - geodesic=true, - pade_order=4, - autodiff=pade_order != 4, - jacobian_structure=true, - hessian_approximation=false, - blas_multithreading=true, + ipopt_options::IpoptOptions=IpoptOptions(), + piccolo_options::PiccoloOptions=PiccoloOptions(), + bound_unitary=piccolo_options.integrator == :exponential, kwargs... ) - if !blas_multithreading - BLAS.set_num_threads(1) - end - - if hessian_approximation - ipopt_options.hessian_approximation = "limited-memory" - end - # Create keys for multiple systems Ũ⃗_keys = [add_suffix(:Ũ⃗, ℓ) for ℓ ∈ system_labels] @@ -110,7 +90,6 @@ function UnitarySamplingProblem( traj = init_trajectory else n_drives = length(systems[1].G_drives) - # TODO: Initial system? traj = initialize_unitary_trajectory( operator, T, @@ -118,14 +97,14 @@ function UnitarySamplingProblem( n_drives, a_bounds, dda_bounds; - free_time=free_time, + free_time=piccolo_options.free_time, Δt_bounds=(Δt_min, Δt_max), - geodesic=geodesic, + geodesic=piccolo_options.geodesic, bound_unitary=bound_unitary, drive_derivative_σ=drive_derivative_σ, a_guess=a_guess, system=systems, - rollout_integrator=rollout_integrator, + rollout_integrator=piccolo_options.rollout_integrator, Ũ⃗_keys=Ũ⃗_keys ) end @@ -159,8 +138,8 @@ function UnitarySamplingProblem( end end - if free_time - if timesteps_all_equal + if piccolo_options.free_time + if piccolo_options.timesteps_all_equal push!(constraints, TimeStepsAllEqualConstraint(:Δt, traj)) end end @@ -180,12 +159,12 @@ function UnitarySamplingProblem( # Integrators unitary_integrators = AbstractIntegrator[] for (sys, Ũ⃗_key) in zip(systems, Ũ⃗_keys) - if integrator == :pade + if piccolo_options.integrator == :pade push!( unitary_integrators, - UnitaryPadeIntegrator(sys, Ũ⃗_key, :a; order=pade_order, autodiff=autodiff) + UnitaryPadeIntegrator(sys, Ũ⃗_key, :a; order=piccolo_options.pade_order) ) - elseif integrator == :exponential + elseif piccolo_options.integrator == :exponential push!( unitary_integrators, UnitaryExponentialIntegrator(sys, Ũ⃗_key, :a) @@ -207,13 +186,8 @@ function UnitarySamplingProblem( J, integrators; constraints=constraints, - max_iter=max_iter, - linear_solver=linear_solver, - verbose=verbose, ipopt_options=ipopt_options, - jacobian_structure=jacobian_structure, - hessian_approximation=hessian_approximation, - eval_hessian=!hessian_approximation, + piccolo_options=piccolo_options, kwargs... ) end @@ -242,70 +216,62 @@ end @testitem "Sample robustness test" begin using Distributions + using Random + Random.seed!(1234) - n_samples = 3 + n_samples = 5 T = 50 Δt = 0.2 timesteps = fill(Δt, T) operator = GATES[:H] systems(ζ) = QuantumSystem(ζ * GATES[:Z], [GATES[:X], GATES[:Y]]) + ip_ops = IpoptOptions(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + pi_ops = PiccoloOptions(verbose=false) + prob = UnitarySamplingProblem( - systems, - Normal(0, 0.1), - n_samples, - operator, - T, - Δt; - verbose=false, - ipopt_options=IpoptOptions(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + systems, Normal(0, 0.05), n_samples, operator, T, Δt, + ipopt_options=ip_ops, piccolo_options=pi_ops ) - solve!(prob, max_iter=20) d_prob = UnitarySmoothPulseProblem( - systems(0), - operator, - T, - Δt; - verbose=false, - ipopt_options=IpoptOptions(print_level=1, recalc_y = "yes", recalc_y_feas_tol = 1e1) + systems(0), operator, T, Δt, + ipopt_options=ip_ops, piccolo_options=pi_ops ) solve!(prob, max_iter=20) # Check that the solution improves over the default - ζ_test = 0.02 + # ------------------------------------------------- + ζ_tests = -0.05:0.01:0.05 Ũ⃗_goal = operator_to_iso_vec(operator) - - Ũ⃗_end = unitary_rollout(prob.trajectory.a, timesteps, systems(ζ_test))[:, end] - fid = unitary_fidelity(Ũ⃗_end, Ũ⃗_goal) - - d_Ũ⃗_end = unitary_rollout(d_prob.trajectory.a, timesteps, systems(ζ_test))[:, end] - default_fid = unitary_fidelity(d_Ũ⃗_end, Ũ⃗_goal) - - @test fid > default_fid + fids = [] + default_fids = [] + for ζ in ζ_tests + Ũ⃗_end = unitary_rollout(prob.trajectory.a, timesteps, systems(ζ))[:, end] + push!(fids, unitary_fidelity(Ũ⃗_end, Ũ⃗_goal)) + + d_Ũ⃗_end = unitary_rollout(d_prob.trajectory.a, timesteps, systems(ζ))[:, end] + push!(default_fids, unitary_fidelity(d_Ũ⃗_end, Ũ⃗_goal)) + end + @test sum(fids) > sum(default_fids) # Check initial guess initialization + # ---------------------------------- a_guess = prob.trajectory.a g1_prob = UnitarySamplingProblem( - [systems(0), systems(0)], - operator, - T, - Δt; - verbose=false, + [systems(0), systems(0)], operator, T, Δt, a_guess=a_guess, + piccolo_options=pi_ops ) @test g1_prob.trajectory.Ũ⃗1 ≈ g1_prob.trajectory.Ũ⃗2 g2_prob = UnitarySamplingProblem( - [systems(0), systems(0.1)], - operator, - T, - Δt; - verbose=false, + [systems(0), systems(0.05)], operator, T, Δt; a_guess=a_guess, + piccolo_options=pi_ops ) @test ~(g2_prob.trajectory.Ũ⃗1 ≈ g2_prob.trajectory.Ũ⃗2) From 667a3e71c67a58197496ee7a068532edbff1a6b4 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 24 Jun 2024 23:56:08 -0500 Subject: [PATCH 35/35] fix tests --- test/direct_sums_test.jl | 25 +++++++++++++++---------- test/problems_test.jl | 6 ++++-- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/test/direct_sums_test.jl b/test/direct_sums_test.jl index 0ee902ba..0e2081b3 100644 --- a/test/direct_sums_test.jl +++ b/test/direct_sums_test.jl @@ -106,13 +106,14 @@ end sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:X], GATES[:Y]]) T = 50 Δt = 0.2 - ops = Options(print_level=1) - prob1 = UnitarySmoothPulseProblem(sys, GATES[:X], T, Δt, free_time=false, ipopt_options=ops) - prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, free_time=false, ipopt_options=ops) + ip_ops = IpoptOptions(print_level=1) + pi_ops = PiccoloOptions(verbose=false, free_time=false) + prob1 = UnitarySmoothPulseProblem(sys, GATES[:X], T, Δt, piccolo_options=pi_ops, ipopt_options=ip_ops) + prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, piccolo_options=pi_ops, ipopt_options=ip_ops) # Direct sum problem with suffix extraction # Note: Turn off control reset - direct_sum_prob = UnitaryDirectSumProblem([prob1, prob2], 0.99, drive_reset_ratio=0.0, ipopt_options=ops) + direct_sum_prob = UnitaryDirectSumProblem([prob1, prob2], 0.99, drive_reset_ratio=0.0, ipopt_options=ip_ops) prob1_got = get_suffix(direct_sum_prob, "1") @test prob1_got.trajectory == add_suffix(prob1.trajectory, "1") @@ -129,9 +130,11 @@ end sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:Y]]) T = 50 Δt = 0.2 - ops = Options(print_level=1) - prob1 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, free_time=false, ipopt_options=ops) - prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, free_time=true, ipopt_options=ops) + ip_ops = IpoptOptions(print_level=1) + pi_false_ops = PiccoloOptions(verbose=false, free_time=false) + pi_true_ops = PiccoloOptions(verbose=false, free_time=true) + prob1 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, piccolo_options=pi_false_ops, ipopt_options=ip_ops) + prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, piccolo_options=pi_true_ops, ipopt_options=ip_ops) suffix = "_new" # UnitaryPadeIntegrator @@ -157,18 +160,20 @@ end sys = QuantumSystem(0.01 * GATES[:Z], [GATES[:Y]]) T = 50 Δt = 0.2 - ops = Options(print_level=1) + ops = IpoptOptions(print_level=1) + pi_false_ops = PiccoloOptions(verbose=false, free_time=false) + pi_true_ops = PiccoloOptions(verbose=false, free_time=true) suffix = "_new" timestep_symbol = :Δt - prob1 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, free_time=false, ipopt_options=ops) + prob1 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, piccolo_options=pi_false_ops, ipopt_options=ops) traj1 = direct_sum(prob1.trajectory, add_suffix(prob1.trajectory, suffix), free_time=true) # Direct sum (shared timestep name) @test get_suffix(traj1, suffix).timestep == timestep_symbol @test get_suffix(traj1, suffix, remove=true).timestep == timestep_symbol - prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, free_time=true, ipopt_options=ops) + prob2 = UnitarySmoothPulseProblem(sys, GATES[:Y], T, Δt, ipopt_options=ops, piccolo_options=pi_true_ops) traj2 = add_suffix(prob2.trajectory, suffix) # Trajectory (unique timestep name) diff --git a/test/problems_test.jl b/test/problems_test.jl index 6586b2e2..35c9d199 100644 --- a/test/problems_test.jl +++ b/test/problems_test.jl @@ -31,14 +31,16 @@ end prob_vanilla = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=4) + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false), ) J_extra = QuadraticSmoothnessRegularizer(:dda, prob_vanilla.trajectory, 10.0) prob_additional = UnitarySmoothPulseProblem( H_drift, H_drives, U_goal, T, Δt, - ipopt_options=IpoptOptions(print_level=4), + ipopt_options=IpoptOptions(print_level=1), + piccolo_options=PiccoloOptions(verbose=false), additional_objective=J_extra, )