diff --git a/src/PiccoloQuantumObjects.jl b/src/PiccoloQuantumObjects.jl index 1aa4e5f..2283020 100644 --- a/src/PiccoloQuantumObjects.jl +++ b/src/PiccoloQuantumObjects.jl @@ -14,4 +14,13 @@ include("quantum_systems.jl") include("embedded_operators.jl") @reexport using .EmbeddedOperators +include("quantum_object_utils.jl") +@reexport using .QuantumObjectUtils + +include("quantum_system_utils.jl") +@reexport using .QuantumSystemUtils + +include("quantum_system_templates/_quantum_system_templates.jl") +@reexport using .QuantumSystemTemplates + end diff --git a/src/quantum_object_utils.jl b/src/quantum_object_utils.jl new file mode 100644 index 0000000..492de6a --- /dev/null +++ b/src/quantum_object_utils.jl @@ -0,0 +1,192 @@ +module QuantumObjectUtils + +export operator_from_string +export ket_from_string +export ket_from_bitstring + +export haar_random +export haar_identity + +export create +export annihilate + +using ..Gates + +using LinearAlgebra +using TestItemRunner + +# TODO: +# [ ] Remove need for oscillator operators (used by tests) +# [ ] Allow multi-character symbols for operators_from_string +# [ ] Remove need for otimes symbol or avoid import conflicts with other packages + + +@doc raw""" +operator_from_string(operator::String; lookup::Dict{Symbol, AbstractMatrix}=PAULIS) + + Reduce the string (each character is one key) via operators from a dictionary. + +""" +function operator_from_string( + operator::String; + lookup::Dict{Symbol, <:AbstractMatrix}=PAULIS +)::Matrix{ComplexF64} + # TODO: allow multi-character keys, ['(', ')'] + + # split string into keys and replace with operators + characters = [Symbol(c) for c ∈ operator] + operators = replace(characters, lookup...) + + return foldr(kron, operators) +end + +function cavity_state(state::Int, levels::Int)::Vector{ComplexF64} + @assert state ≤ levels - 1 "Level $state is not allowed for $levels levels" + ket = zeros(levels) + ket[state + 1] = 1 + return ket +end + +@doc raw""" + ket_from_string( + ket::String, + levels::Vector{Int}; + level_dict=Dict(:g => 0, :e => 1, :f => 2, :h => 2), + return_states=false + ) + +Construct a quantum state from a string ket representation. + +# Example + +# TODO: add example +""" +function ket_from_string( + ket::String, + levels::Vector{Int}; + level_dict=Dict(:g => 0, :e => 1, :f => 2, :h => 2), + return_states=false +)::Vector{ComplexF64} + kets = [] + + for x ∈ split(ket, ['(', ')']) + if x == "" + continue + elseif all(Symbol(xᵢ) ∈ keys(level_dict) for xᵢ ∈ x) + append!(kets, x) + elseif occursin("+", x) + superposition = split(x, '+') + @assert all(all(Symbol(xᵢ) ∈ keys(level_dict) for xᵢ ∈ x) for x ∈ superposition) "Invalid ket: $x" + @assert length(superposition) == 2 "Only two states can be superposed for now" + push!(kets, x) + else + error("Invalid ket: $x") + end + end + + states = [] + + for (ψᵢ, l) ∈ zip(kets, levels) + if ψᵢ isa AbstractString && occursin("+", ψᵢ) + superposition = split(ψᵢ, '+') + superposition_states = [level_dict[Symbol(x)] for x ∈ superposition] + @assert all(state ≤ l - 1 for state ∈ superposition_states) "Level $ψᵢ is not allowed for $l levels" + superposition_state = sum([ + cavity_state(state, l) for state ∈ superposition_states + ]) + normalize!(superposition_state) + push!(states, superposition_state) + else + state = level_dict[Symbol(ψᵢ)] + @assert state ≤ l - 1 "Level $ψᵢ is not allowed for $l levels" + push!(states, cavity_state(state, l)) + end + end + + if return_states + return states + else + return kron([1.0], states...) + end +end + +@doc raw""" + ket_from_bitstring(ket::String) + +Get the state vector for a qubit system given a ket string `ket` of 0s and 1s. +""" +function ket_from_bitstring(ket::String)::Vector{ComplexF64} + cs = [c for c ∈ ket] + @assert all(c ∈ "01" for c ∈ cs) + states = [c == '0' ? [1, 0] : [0, 1] for c ∈ cs] + return foldr(kron, states) +end + +### +### Random operators +### + +@doc raw""" + haar_random(n::Int) + +Generate a random unitary matrix using the Haar measure for an `n`-dimensional system. +""" +function haar_random(n::Int) + # Ginibre matrix + Z = (randn(n, n) + im * randn(n, n)) / √2 + F = qr(Z) + # QR correction (R main diagonal is real, strictly positive) + Λ = diagm(diag(F.R) ./ abs.(diag(F.R))) + return F.Q * Λ +end + +@doc raw""" + haar_identity(n::Int, radius::Number) + +Generate a random unitary matrix close to the identity matrix using the Haar measure for an `n`-dimensional system with a given `radius`. +""" +function haar_identity(n::Int, radius::Number) + # Ginibre matrix + Z = (I + radius * (randn(n, n) + im * randn(n, n)) / √2) / (1 + radius) + F = qr(Z) + # QR correction (R main diagonal is real, strictly positive) + Λ = diagm(diag(F.R) ./ abs.(diag(F.R))) + return F.Q * Λ +end + +### +### Oscillator operators +### + +@doc raw""" + annihilate(levels::Int) + +Get the annihilation operator for a system with `levels` levels. +""" +function annihilate(levels::Int)::Matrix{ComplexF64} + return diagm(1 => map(sqrt, 1:levels - 1)) +end + +@doc raw""" + create(levels::Int) + +Get the creation operator for a system with `levels` levels. +""" +function create(levels::Int) + return collect(annihilate(levels)') +end + +# ============================================================================= # + +@testitem "Test ket_from_bitstring function" begin + using LinearAlgebra + @test ket_from_bitstring("0") == [1, 0] + @test ket_from_bitstring("1") == [0, 1] + @test ket_from_bitstring("00") == [1, 0, 0, 0] + @test ket_from_bitstring("01") == [0, 1, 0, 0] + @test ket_from_bitstring("10") == [0, 0, 1, 0] + @test ket_from_bitstring("11") == [0, 0, 0, 1] +end + + +end diff --git a/src/quantum_system_templates/_quantum_system_templates.jl b/src/quantum_system_templates/_quantum_system_templates.jl new file mode 100644 index 0000000..cca2e3e --- /dev/null +++ b/src/quantum_system_templates/_quantum_system_templates.jl @@ -0,0 +1,18 @@ +module QuantumSystemTemplates + +export TransmonSystem +export TransmonDipoleCoupling +export MultiTransmonSystem +export RydbergChainSystem +export QuantumOpticsSystem + +using ..QuantumSystems +using ..QuantumObjectUtils + +using LinearAlgebra +using TestItemRunner + +include("transmons.jl") +include("rydberg.jl") + +end diff --git a/src/quantum_system_templates/rydberg.jl b/src/quantum_system_templates/rydberg.jl new file mode 100644 index 0000000..98b6175 --- /dev/null +++ b/src/quantum_system_templates/rydberg.jl @@ -0,0 +1,126 @@ + +function generate_pattern(N::Int, i::Int) + # Create an array filled with 'I' + qubits = fill('I', N) + # Insert 'n' at position i and i+1, ensuring it doesn't exceed the array bounds + if i <= N && i+1 <= N + qubits[i] = 'n' + qubits[i+1] = 'n' + end + return join(qubits) +end +function generate_pattern_with_gap(N::Int, i::Int, gap::Int) + # Create an array filled with 'I' + qubits = fill('I', N) + # Insert 'n' at position i and i+gap+1, ensuring it doesn't exceed the array bounds + if i <= N && i+gap+1 <= N + qubits[i] = 'n' + qubits[i+gap+1] = 'n' + end + return join(qubits) +end + +""" +Embed a character into a string at a specific position. +""" +function lift(x::Char,i::Int, N::Int) + qubits = fill('I', N) + qubits[i] = x + return join(qubits) +end +@doc raw""" + RydbergChainSystem(; + N::Int=3, # number of atoms + C::Float64=862690*2π, + distance::Float64=10.0, # μm + cutoff_order::Int=2, # 1 is nearest neighbor, 2 is next-nearest neighbor, etc. + local_detune::Bool=false, # If true, include one local detuning pattern. + all2all::Bool=true, # If true, include all-to-all interactions. + ignore_Y_drive::Bool=false, # If true, ignore the Y drive. (In the experiments, X&Y drives are implemented by Rabi amplitude and its phase.) + ) -> QuantumSystem + +Returns a `QuantumSystem` object for the Rydberg atom chain in the spin basis + |g⟩ = |0⟩ = [1, 0], |r⟩ = |1⟩ = [0, 1]. + +```math +H = \sum_i 0.5*\Omega_i(t)\cos(\phi_i(t)) \sigma_i^x - 0.5*\Omega_i(t)\sin(\phi_i(t)) \sigma_i^y - \sum_i \Delta_i(t)n_i + \sum_{i [1 0; 0 1], :X => [0 1; 1 0], :Y => [0 -im; im 0], :Z => [1 0; 0 -1], :n => [0 0; 0 1]) + + if all2all + H_drift = zeros(ComplexF64, 2^N, 2^N) + for gap in 0:N-2 + for i in 1:N-gap-1 + H_drift += C * operator_from_string( + generate_pattern_with_gap(N, i, gap), + lookup = PAULIS + ) / ((gap + 1) * distance)^6 + end + end + else + if cutoff_order == 1 + H_drift = sum([C*operator_from_string(generate_pattern(N,i),lookup=PAULIS)/(distance^6) for i in 1:N-1]) + elseif cutoff_order == 2 + H_drift = sum([C*operator_from_string(generate_pattern(N,i),lookup=PAULIS)/(distance^6) for i in 1:N-1]) + H_drift += sum([C*operator_from_string(generate_pattern_with_gap(N,i,1),lookup=PAULIS)/((2*distance)^6) for i in 1:N-2]) + else + error("Higher cutoff order not supported") + end + end + + H_drives = Matrix{ComplexF64}[] + + # Add global X drive + Hx = sum([0.5*operator_from_string(lift('X',i,N), lookup=PAULIS) for i in 1:N]) + push!(H_drives, Hx) + + if !ignore_Y_drive + # Add global Y drive + Hy = sum([0.5*operator_from_string(lift('Y',i,N), lookup=PAULIS) for i in 1:N]) + push!(H_drives, Hy) + end + + # Add global detuning + H_detune = -sum([operator_from_string(lift('n',i,N), lookup=PAULIS) for i in 1:N]) + push!(H_drives, H_detune) + + params = Dict{Symbol, Any}( + :N => N, + :C => C, + :distance => distance, + :cutoff_order => cutoff_order, + :local_detune => local_detune, + :all2all => all2all, + ) + + return QuantumSystem( + H_drift, + H_drives; + constructor=RydbergChainSystem, + params=params, + ) +end + +@testitem "Rydberg system test" begin + @test RydbergChainSystem(N=3,cutoff_order=2,all2all=false) isa QuantumSystem +end diff --git a/src/quantum_system_templates/transmons.jl b/src/quantum_system_templates/transmons.jl new file mode 100644 index 0000000..7c75405 --- /dev/null +++ b/src/quantum_system_templates/transmons.jl @@ -0,0 +1,239 @@ +@doc raw""" + TransmonSystem(; + ω::Float64=4.4153, # GHz + δ::Float64=0.17215, # GHz + levels::Int=3, + lab_frame::Bool=false, + frame_ω::Float64=ω, + ) -> QuantumSystem + +Returns a `QuantumSystem` object for a transmon qubit, with the Hamiltonian + +```math +H = \omega a^\dagger a - \frac{\delta}{2} a^\dagger a^\dagger a a +``` + +where `a` is the annihilation operator. + +# Keyword Arguments +- `ω`: The frequency of the transmon, in GHz. +- `δ`: The anharmonicity of the transmon, in GHz. +- `levels`: The number of levels in the transmon. +- `lab_frame`: Whether to use the lab frame Hamiltonian, or an ω-rotating frame. +- `frame_ω`: The frequency of the rotating frame, in GHz. +- `mutiply_by_2π`: Whether to multiply the Hamiltonian by 2π, set to true by default because the frequency is in GHz. +- `lab_frame_type`: The type of lab frame Hamiltonian to use, one of (:duffing, :quartic, :cosine). +- `drives`: Whether to include drives in the Hamiltonian. +""" +function TransmonSystem(; + ω::Float64=4.0, # GHz + δ::Float64=0.2, # GHz + levels::Int=3, + lab_frame::Bool=false, + frame_ω::Float64=lab_frame ? 0.0 : ω, + mutiply_by_2π::Bool=true, + lab_frame_type::Symbol=:duffing, + drives::Bool=true, +) + + @assert lab_frame_type ∈ (:duffing, :quartic, :cosine) "lab_frame_type must be one of (:duffing, :quartic, :cosine)" + + if lab_frame + if frame_ω ≉ 0.0 + frame_ω = 0.0 + end + end + + if frame_ω ≉ 0.0 + lab_frame = false + end + + a = annihilate(levels) + + if lab_frame + if lab_frame_type == :duffing + H_drift = ω * a' * a - δ / 2 * a' * a' * a * a + elseif lab_frame_type == :quartic + ω₀ = ω + δ + H_drift = ω₀ * a' * a - δ / 12 * (a + a')^4 + elseif lab_frame_type == :cosine + ω₀ = ω + δ + E_C = δ + E_J = ω₀^2 / 8E_C + n̂ = im / 2 * (E_J / 2E_C)^(1/4) * (a - a') + φ̂ = (2E_C / E_J)^(1/4) * (a + a') + H_drift = 4 * E_C * n̂^2 - E_J * cos(φ̂) + # H_drift = 4 * E_C * n̂^2 - E_J * (I - φ̂^2 / 2 + φ̂^4 / 24) + end + else + H_drift = (ω - frame_ω) * a' * a - δ / 2 * a' * a' * a * a + end + + if drives + H_drives = [a + a', 1.0im * (a - a')] + else + H_drives = Matrix{ComplexF64}[] + end + + if mutiply_by_2π + H_drift *= 2π + H_drives .*= 2π + end + + params = Dict{Symbol, Any}( + :ω => ω, + :δ => δ, + :levels => levels, + :lab_frame => lab_frame, + :frame_ω => frame_ω, + :mutiply_by_2π => mutiply_by_2π, + :lab_frame_type => lab_frame_type, + :drives => drives, + ) + + return QuantumSystem( + H_drift, + H_drives; + constructor=TransmonSystem, + params=params, + ) +end + +@doc raw""" + TransmonDipoleCoupling( + g_ij::Float64, + pair::Tuple{Int, Int}, + subsystem_levels::Vector{Int}; + lab_frame::Bool=false, + ) -> QuantumSystemCoupling + + TransmonDipoleCoupling( + g_ij::Float64, + pair::Tuple{Int, Int}, + sub_systems::Vector{QuantumSystem}; + kwargs... + ) -> QuantumSystemCoupling + +Returns a `QuantumSystemCoupling` object for a transmon qubit. In the lab frame, the Hamiltonian coupling term is + +```math +H = g_{ij} (a_i + a_i^\dagger) (a_j + a_j^\dagger) +``` + +In the rotating frame, the Hamiltonian coupling term is + +```math +H = g_{ij} (a_i a_j^\dagger + a_i^\dagger a_j) +``` + +where `a_i` is the annihilation operator for the `i`th transmon. + +""" +function TransmonDipoleCoupling end + +function TransmonDipoleCoupling( + g_ij::Float64, + pair::Tuple{Int, Int}, + subsystem_levels::Vector{Int}; + lab_frame::Bool=false, + mulitply_by_2π::Bool=true, +) + i, j = pair + a_i = lift(annihilate(subsystem_levels[i]), i, subsystem_levels) + a_j = lift(annihilate(subsystem_levels[j]), j, subsystem_levels) + + if lab_frame + op = g_ij * (a_i + a_i') * (a_j + a_j') + else + op = g_ij * (a_i * a_j' + a_i' * a_j) + end + + if mulitply_by_2π + op *= 2π + end + + params = Dict{Symbol, Any}( + :lab_frame => lab_frame, + :mulitply_by_2π => mulitply_by_2π, + ) + + return QuantumSystemCoupling( + op, + g_ij, + pair, + subsystem_levels, + TransmonDipoleCoupling, + params + ) +end + +function TransmonDipoleCoupling( + g_ij::Float64, + pair::Tuple{Int, Int}, + sub_systems::Vector{QuantumSystem}; + kwargs... +) + subsystem_levels = [sys.levels for sys ∈ sub_systems] + return TransmonDipoleCoupling(g_ij, pair, subsystem_levels; kwargs...) +end + +""" + MultiTransmonSystem( + ωs::AbstractVector{Float64}, + δs::AbstractVector{Float64}, + gs::AbstractMatrix{Float64}; + levels_per_transmon::Int = 3, + subsystem_levels::AbstractVector{Int} = fill(levels_per_transmon, length(ωs)), + lab_frame=false, + subsystems::AbstractVector{Int} = 1:length(ωs), + subsystem_drive_indices::AbstractVector{Int} = 1:length(ωs), + kwargs... + ) -> CompositeQuantumSystem + +Returns a `CompositeQuantumSystem` object for a multi-transmon system. +""" +function MultiTransmonSystem( + ωs::AbstractVector{Float64}, + δs::AbstractVector{Float64}, + gs::AbstractMatrix{Float64}; + levels_per_transmon::Int = 3, + subsystem_levels::AbstractVector{Int} = fill(levels_per_transmon, length(ωs)), + lab_frame=false, + subsystems::AbstractVector{Int} = 1:length(ωs), + subsystem_drive_indices::AbstractVector{Int} = 1:length(ωs), + kwargs... +) + n_subsystems = length(ωs) + @assert length(δs) == n_subsystems + @assert size(gs) == (n_subsystems, n_subsystems) + + systems = QuantumSystem[] + + for (i, (ω, δ, levels)) ∈ enumerate(zip(ωs, δs, subsystem_levels)) + if i ∈ subsystems + sysᵢ = TransmonSystem( + levels=levels, + ω=ω, + δ=δ, + lab_frame=lab_frame, + drives=i ∈ subsystem_drive_indices + ) + push!(systems, sysᵢ) + end + end + + couplings = QuantumSystemCoupling[] + + for i = 1:n_subsystems-1 + for j = i+1:n_subsystems + if i ∈ subsystems && j ∈ subsystems + push!( + couplings, + TransmonDipoleCoupling(gs[i, j], (i, j), systems; lab_frame=lab_frame) + ) + end + end + end + + return CompositeQuantumSystem(systems, couplings) +end diff --git a/src/quantum_system_utils.jl b/src/quantum_system_utils.jl new file mode 100644 index 0000000..0b5a71e --- /dev/null +++ b/src/quantum_system_utils.jl @@ -0,0 +1,375 @@ +module QuantumSystemUtils + +export operator_algebra +export is_reachable + +using ..EmbeddedOperators +using ..QuantumObjectUtils +using ..QuantumSystems +using ..Gates + +using LinearAlgebra +using SparseArrays +using TestItemRunner + + + +commutator(A::AbstractMatrix, B::AbstractMatrix) = A * B - B * A + +# TODO: Any difference between LinearAlgebra::ishermitian? +is_hermitian(H::AbstractMatrix; atol=eps(Float32)) = + all(isapprox.(H - H', 0.0, atol=atol)) + +function is_linearly_dependent(basis::Vector{<:AbstractMatrix}; kwargs...) + return is_linearly_dependent(stack(vec.(basis)); kwargs...) +end + +function is_linearly_dependent(M::AbstractMatrix; eps=eps(Float32), verbose=true) + if size(M, 2) > size(M, 1) + if verbose + println("Linearly dependent because columns > rows, $(size(M, 2)) > $(size(M, 1)).") + end + return true + end + # QR decomposition has a zero R on diagonal if linearly dependent + val = minimum(abs.(diag(qr(M).R))) + return isapprox(val, 0.0, atol=eps) +end + +function linearly_independent_indices( + basis::Vector{<:AbstractMatrix}; + order=1:length(basis), + kwargs... +) + @assert issetequal(order, 1:length(basis)) "Order must enumerate entire basis." + bᵢ = Int[] + for i ∈ order + if !is_linearly_dependent([basis[bᵢ]..., basis[i]]; kwargs...) + push!(bᵢ, i) + end + end + return bᵢ +end + +function linearly_independent_subset(basis::Vector{<:AbstractMatrix}; kwargs...) + bᵢ = linearly_independent_indices(basis; kwargs...) + return deepcopy(basis[bᵢ]) +end + +function linearly_independent_subset!(basis::Vector{<:AbstractMatrix}; kwargs...) + bᵢ = linearly_independent_indices(basis; kwargs...) + deleteat!(basis, setdiff(1:length(basis), bᵢ)) + return nothing +end + +traceless(M::AbstractMatrix) = M - tr(M) * I / size(M, 1) + +""" +operator_algebra(generators; return_layers=false, normalize=false, verbose=false, remove_trace=true) + + Compute the Lie algebra basis for the given generators. + If return_layers is true, the Lie tree layers are also returned. + + Can normalize the basis and enforce traceless generators. + + # Arguments + - `generators::Vector{<:AbstractMatrix}`: generators of the Lie algebra + - `return_layers::Bool=false`: return the Lie tree layers + - `normalize::Bool=false`: normalize the basis + - `verbose::Bool=false`: print debug information + - `remove_trace::Bool=true`: remove trace from generators +""" +function operator_algebra( + generators::Vector{<:AbstractMatrix{T}}; + return_layers=false, + normalize=false, + verbose=false, + remove_trace=true +) where T<:Number + # Initialize basis (traceless, normalized) + basis = normalize ? [g / norm(g) for g ∈ generators] : deepcopy(generators) + if remove_trace + @. basis = traceless(basis) + end + + # Initialize layers + current_layer = deepcopy(basis) + if return_layers + all_layers = Vector{Matrix{T}}[deepcopy(basis)] + end + + ℓ = 1 + if verbose + println("ℓ = $ℓ") + end + if is_linearly_dependent(basis) + println("Linearly dependent generators.") + else + # Note: Use left normalized commutators + # Note: Jacobi identity is not checked + need_basis = true + algebra_dim = size(first(generators), 1)^2 - 1 + while need_basis + layer = Matrix{T}[] + # Repeat commutators until no new operators are found + for op ∈ current_layer + for gen ∈ generators + if !need_basis + continue + end + + test = commutator(gen, op) + if all(test .≈ 0) + continue + # Current basis is assumed to be linearly independent + elseif is_linearly_dependent([basis..., test], verbose=verbose) + continue + else + test .= is_hermitian(test) ? test : im * test + test .= normalize ? test / norm(test) : test + push!(basis, test) + push!(layer, test) + need_basis = length(basis) < algebra_dim ? true : false + end + end + end + + if isempty(layer) + if verbose + println("Subspace termination.") + end + break + else + current_layer = layer + ℓ += 1 + if verbose + println("ℓ = $ℓ") + end + end + + if return_layers + append!(all_layers, [current_layer]) + end + end + end + + if return_layers + return basis, all_layers + else + return basis + end +end + +function fit_gen_to_basis( + gen::AbstractMatrix{<:T}, + basis::AbstractVector{<:AbstractMatrix{<:T}} +) where T<:Number + A = stack(vec.(basis)) + b = vec(gen) + return A \ b +end + +function is_in_span( + gen::AbstractMatrix, + basis::AbstractVector{<:AbstractMatrix}; + subspace::AbstractVector{<:Int}=1:size(gen, 1), + atol=eps(Float32), + return_effective_gen=false, +) + g_basis = [deepcopy(b[subspace, subspace]) for b ∈ basis] + linearly_independent_subset!(g_basis) + # Traceless basis needs traceless fit + x = fit_gen_to_basis(gen, g_basis) + g_eff = sum(x .* g_basis) + ε = norm(g_eff - gen, 2) + if return_effective_gen + return ε < atol, g_eff + else + return ε < atol + end +end + +""" + is_reachable(gate, hamiltonians; kwargs...) + +Check if the gate is reachable from the given generators. If subspace_indices are provided, +then the gate should be given in the subspace. + +# Arguments +- `gate::AbstractMatrix`: target gate +- `hamiltonians::AbstractVector{<:AbstractMatrix}`: generators of the Lie algebra + +# Keyword Arguments +- `subspace::AbstractVector{<:Int}=1:size(gate, 1)`: subspace indices +- `compute_basis::Bool=true`: compute the basis +- `remove_trace::Bool=true`: remove trace from generators +- `verbose::Bool=false`: print debug information +- `atol::Float32=eps(Float32)`: absolute tolerance +""" +function is_reachable( + gate::AbstractMatrix, + hamiltonians::AbstractVector{<:AbstractMatrix}; + subspace::AbstractVector{<:Int}=1:size(gate, 1), + compute_basis=true, + remove_trace=true, + verbose=false, + atol=eps(Float32) +) + @assert size(gate, 1) == length(subspace) "Gate must be given in the subspace." + generator = im * log(gate) + + if remove_trace + generator = traceless(generator) + end + + if compute_basis + basis = operator_algebra(hamiltonians, remove_trace=remove_trace, verbose=verbose) + else + basis = hamiltonians + end + + return is_in_span( + generator, + basis, + subspace=subspace, + 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 + +function is_reachable( + gate::EmbeddedOperator, + system::QuantumSystem; + kwargs... +) + return is_reachable( + unembed(gate), + system; + subspace=gate.subspace_indices, + kwargs... + ) +end + +# ============================================================================= # + +@testitem "Lie algebra basis" begin + # Check 1 qubit with complete basis + gen = operator_from_string.(["X", "Y"]) + basis = operator_algebra(gen, return_layers=false, verbose=false) + @test length(basis) == size(first(gen), 1)^2-1 + + # Check 1 qubit with complete basis and layers + basis, layers = operator_algebra(gen, return_layers=true, verbose=false) + @test length(basis) == size(first(gen), 1)^2-1 + + # Check 1 qubit with subspace + gen = operator_from_string.(["X"]) + basis = operator_algebra(gen, verbose=false) + @test length(basis) == 1 + + # Check 2 qubit with complete basis + gen = operator_from_string.(["XX", "YY", "XI", "YI", "IY", "IX"]) + basis = operator_algebra(gen, verbose=false) + @test length(basis) == size(first(gen), 1)^2-1 + + # Check 2 qubit with linearly dependent basis + gen = operator_from_string.(["XX", "YY", "XI", "XI", "IY", "IX"]) + basis = operator_algebra(gen, verbose=false) + @test length(basis) == length(gen) + + # Check 2 qubit with pair of 1-qubit subspaces + gen = operator_from_string.(["XI", "YI", "IY", "IX"]) + basis = operator_algebra(gen, verbose=false) + @test length(basis) == 2 * (2^2 - 1) +end + + +@testitem "Lie Algebra reachability" begin + using LinearAlgebra + + H_ops = Dict( + "X" => GATES[:X], + "Y" => GATES[:Y], + "Z" => GATES[:Z] + ) + + # Check 1 qubit with complete basis + gen = operator_from_string.(["X", "Y"]) + target = H_ops["Z"] + @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] + IX = GATES[:I] ⊗ GATES[:X] + YI = GATES[:Y] ⊗ GATES[:I] + IY = GATES[:I] ⊗ GATES[:Y] + XX = GATES[:X] ⊗ GATES[:X] + YY = GATES[:Y] ⊗ GATES[:Y] + ZI = GATES[:Z] ⊗ GATES[:I] + IZ = GATES[:I] ⊗ GATES[:Z] + ZZ = GATES[:Z] ⊗ GATES[:Z] + + complete_gen = [XX+YY, XI, YI, IX, IY] + incomplete_gen = [XI, ZZ] + r = [0, 1, 2, 3, 4] + r /= norm(r) + R2 = exp(-im * sum([θ * H for (H, θ) in zip(complete_gen, r)])) + CZ = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 -1] + CX = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0] + + # Pass + @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(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 + # TODO: implement tests +end + +end