Enable more Aqua tests by solving piracy and ambiguities issues (#143)

Co-authored-by: Stefan Krastanov <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>

CHANGELOG.md

@@ -1,7 +1,8 @@
 # News
 
-## v0.4.2 - dev
+## v0.5.0 - 2024-08-11
 
+- `observable` now takes a default value as a kwarg, i.e., you need to make the substitution `observable(regs, obs, 0.0; time)` ↦ `observable(regs, obs; something=0.0, time)`
 - Bump QuantumSymbolics and QuantumOpticsBase compat bound and bump julia compat to 1.10.
 
 ## v0.4.1 - 2024-06-05

Project.toml

@@ -49,10 +49,10 @@ Makie = "0.20, 0.21"
 NetworkLayout = "0.4.4"
 PrecompileTools = "1"
 Printf = "1"
-QuantumClifford = "0.8.20, 0.9"
-QuantumInterface = "0.3.4"
-QuantumOptics = "1.0.5"
-QuantumOpticsBase = "0.4.22, 0.5"
+QuantumClifford = "0.9.9"
+QuantumInterface = "0.3.5"
+QuantumOptics = "1.1.0"
+QuantumOpticsBase = "0.5.3"
 QuantumSymbolics = "0.3"
 Random = "1"
 Reexport = "1.2.2"

docs/src/register_interface.md

@@ -138,13 +138,13 @@ Measure a quantum observable. The dispatch down the call three is very similar t
 observable
 ```
 
-#### `observable(refs::Tuple{Vararg{RegRef, N}}, obs, something=nothing; time)`
+#### `observable(refs::Tuple{Vararg{RegRef, N}}, obs; something=nothing, time=nothing)`
 
 Calculate the value of an observable on the state in the sequence of [`RegRef`](@ref)s at a specified `time`. If these registers are not instantiated, return `something`.
 
 `refs` can also be `Tuple{Vararg{RegRef, N}}` or a single [`RegRef`](@ref).
 
-#### `observable(regs::Vector{Register}, indices, obs, something=nothing; time)`
+#### `observable(regs::Vector{Register}, indices, obs; something=nothing, time=nothing)`
 
 `indices` refers to the slots inside of the given `regs`.
 
@@ -167,8 +167,8 @@ If `operation<:Symbolic`, then an `express(obs, repr, ::UseAsObservable)` call i
 ```@raw html
 <div class="mermaid">
 flowchart TB
-  A["<code>observable(refs::Vector{RegRef}, obs, something=nothing; time)</code>"]
-  B["<code>observable(regs::Vector{Register}, indices, obs, something=nothing; time)</code>"]
+  A["<code>observable(refs::Vector{RegRef}, obs; something=nothing, time)</code>"]
+  B["<code>observable(regs::Vector{Register}, indices, obs; something=nothing, time)</code>"]
   subgraph TOP [lower from registers to states]
     direction LR
     B1["<code>uptotime!</code>"]

docs/src/symbolics.md

@@ -203,4 +203,4 @@ express(MixedState(X1)/2+SProjector(Z1)/2, CliffordRepr())
 
 !!! warning "Stabilizer state expressions"
 
-    The state written as $\frac{|Z₁⟩⊗|Z₁⟩+|Z₂⟩⊗|Z₂⟩}{√2}$ is a well known stabilizer state, namely a Bell state. However, automatically expressing it as a stabilizer is a prohibitively expensive computational operation in general. We do not perform that computation automatically. If you want to ensure that states you define can be automatically converted to tableaux for Clifford simulations, avoid using sumation of kets. On the other hand, in all of our Clifford Monte-Carlo simulations, `⊗` is fully supported, as well as [`SProjector`](@ref), [`MixedState`](@ref), [`StabilizerState`](@ref), and sumation of density matrices.
+    The state written as $\frac{|Z₁⟩⊗|Z₁⟩+|Z₂⟩⊗|Z₂⟩}{√2}$ is a well known stabilizer state, namely a Bell state. However, automatically expressing it as a stabilizer is a prohibitively expensive computational operation in general. We do not perform that computation automatically. If you want to ensure that states you define can be automatically converted to tableaux for Clifford simulations, avoid using summation of kets. On the other hand, in all of our Clifford Monte-Carlo simulations, `⊗` is fully supported, as well as [`SProjector`](@ref), [`MixedState`](@ref), [`StabilizerState`](@ref), and summation of density matrices.

examples/colorcentermodularcluster/2_real_time_visualization.jl

@@ -68,7 +68,7 @@ record(F, "colorcentermodularcluster-02.simdashboard.mp4", step_ts, framerate=FR
         l = length(neighs)
         obs = observables[l]
         regs = [net[i, 2] for i in [v, neighs...]]
-        real(observable(regs, obs, 0.0; time=now(sim)))
+        real(observable(regs, obs; something=0.0, time=now(sim)))
     end
     linkcolors[] .= fid
     push!(fids[],mean(fid))

examples/colorcentermodularcluster/3_makie_interactive.jl

@@ -209,7 +209,7 @@ function continue_singlerun!(sim, net,
             l = length(neighs)
             obs = observables[l]
             regs = [net[i, 2] for i in [v, neighs...]]
-            real(observable(regs, obs, 0.0; time=now(sim)))
+            real(observable(regs, obs; something=0.0, time=now(sim)))
         end)
         linkcolors[] .= fid
         push!(fids[],mean(fid))

examples/congestionchain/setup.jl

@@ -191,8 +191,8 @@ end
         reg2 = network[node2]
         @yield request(reg1[q1]) & request(reg2[q2])
         uptotime!((reg1[q1], reg2[q2]), now(sim))
-        fXX = real(observable((reg1[q1],reg2[q2]), XX, 0.0; time=now(sim)))
-        fZZ = real(observable((reg1[q1],reg2[q2]), ZZ, 0.0; time=now(sim)))
+        fXX = real(observable((reg1[q1],reg2[q2]), XX; something=0.0, time=now(sim)))
+        fZZ = real(observable((reg1[q1],reg2[q2]), ZZ; something=0.0, time=now(sim)))
         push!(fidelityXXlog[], fXX)
         push!(fidelityZZlog[], fZZ)
         push!(timelog[], now(sim)-last_success)

examples/firstgenrepeater/4_visualization.jl

@@ -56,8 +56,8 @@ step_ts = range(0, 100, step=0.1)
 record(fig, "firstgenrepeater-07.observable.mp4", step_ts; framerate=10, visible=true) do t
     run(sim, t)
 
-    fXX = real(observable(registers[[1,last]], [2,2], XX, 0.0; time=t))
-    fZZ = real(observable(registers[[1,last]], [2,2], ZZ, 0.0; time=t))
+    fXX = real(observable(registers[[1,last]], [2,2], XX; something=0.0, time=t))
+    fZZ = real(observable(registers[[1,last]], [2,2], ZZ; something=0.0, time=t))
     push!(fidXX[],fXX)
     push!(fidZZ[],fZZ)
     push!(ts[],t)

examples/firstgenrepeater/5_clifford_full_example.jl

@@ -56,8 +56,8 @@ step_ts = range(0, 100, step=0.1)
 record(fig, "firstgenrepeater-08.clifford.mp4", step_ts, framerate=10, visible=true) do t
     run(sim, t)
 
-    fXX = real(observable(registers[[1,last]], [2,2], XX, 0.0; time=t))
-    fZZ = real(observable(registers[[1,last]], [2,2], ZZ, 0.0; time=t))
+    fXX = real(observable(registers[[1,last]], [2,2], XX; something=0.0, time=t))
+    fZZ = real(observable(registers[[1,last]], [2,2], ZZ; something=0.0, time=t))
     push!(fidXX[],fXX)
     push!(fidZZ[],fZZ)
     push!(ts[],t)

examples/firstgenrepeater/6.1_compare_formalisms_noplot.jl

@@ -47,8 +47,8 @@ function monte_carlo_trajectory(;
     for t in sampled_times
         run(sim, t)
 
-        fXX = real(observable(registers[[1,len]], [2,2], XX, 0.0; time=t))
-        fZZ = real(observable(registers[[1,len]], [2,2], ZZ, 0.0; time=t))
+        fXX = real(observable(registers[[1,len]], [2,2], XX; something=0.0, time=t))
+        fZZ = real(observable(registers[[1,len]], [2,2], ZZ; something=0.0, time=t))
         push!(fidXX,fXX)
         push!(fidZZ,fZZ)
     end

examples/firstgenrepeater/6_compare_formalisms.jl

@@ -49,8 +49,8 @@ function monte_carlo_trajectory(;
     for t in sampled_times
         run(sim, t)
 
-        fXX = real(observable(registers[[1,len]], [2,2], XX, 0.0; time=t))
-        fZZ = real(observable(registers[[1,len]], [2,2], ZZ, 0.0; time=t))
+        fXX = real(observable(registers[[1,len]], [2,2], XX; something=0.0, time=t))
+        fZZ = real(observable(registers[[1,len]], [2,2], ZZ; something=0.0, time=t))
         push!(fidXX,fXX)
         push!(fidZZ,fZZ)
     end

src/QuantumSavory.jl

@@ -34,7 +34,7 @@ export project_traceout! #TODO should move to QuantumInterface
 
 using QuantumSymbolics:
     AbstractRepresentation, AbstractUse,
-    CliffordRepr, QuantumOpticsRepr, QuantumMCRepr,
+    CliffordRepr, consistent_representation, QuantumOpticsRepr, QuantumMCRepr,
     metadata, istree, operation, arguments, Symbolic, # from Symbolics
     HGate, XGate, YGate, ZGate, CPHASEGate, CNOTGate,
     XBasisState, YBasisState, ZBasisState,
@@ -102,6 +102,8 @@ include("noninstant.jl")
 include("backends/quantumoptics/quantumoptics.jl")
 include("backends/clifford/clifford.jl")
 
+include("ambiguity_fix.jl")
+
 include("concurrentsim.jl")
 
 include("plots.jl")

src/ambiguity_fix.jl

@@ -0,0 +1,6 @@
+
+Register(traits, ::Union{Tuple{}, AbstractVector{Union{}}}) = Register(traits)
+subsystemcompose() = error()
+project_traceout!(::Union{Ket,Operator}, ::Int, ::Union{Tuple{}, AbstractVector{Union{}}}) = error()
+QuantumClifford.apply!(::QuantumClifford.MixedDestabilizer, ::QuantumSavory.QCGateSequence, ::Type{<:QuantumClifford.AbstractSymbolicOperator}) = error()
+observable(::Union{Tuple{}, AbstractVector{Union{}}}, ::Base.AbstractVecOrTuple{Int}) = error()

src/backends/clifford/clifford.jl

@@ -3,12 +3,8 @@ import QuantumClifford: MixedDestabilizer
 
 subsystemcompose(states::QuantumClifford.MixedDestabilizer...) = QuantumClifford.tensor(states...)
 
-nsubsystems(state::QuantumClifford.MixedDestabilizer) = QuantumClifford.nqubits(state)
-
 default_repr(::QuantumClifford.MixedDestabilizer) = CliffordRepr()
 
-apply!(state::QuantumClifford.MixedDestabilizer, indices, operation::Type{<:QuantumClifford.AbstractSymbolicOperator}) = QuantumClifford.apply!(state, operation, indices)
-
 ispadded(::QuantumClifford.MixedDestabilizer) = false
 
 const _qc_l = copy(express(Z1, CliffordRepr()))

src/backends/clifford/express.jl

@@ -1,4 +1,5 @@
-function project_traceout!(state::QuantumClifford.MixedDestabilizer,stateindex,basis::Symbolic{AbstractOperator})
+
+function project_traceout!(state::QuantumClifford.MixedDestabilizer,stateindex::Int,basis::Symbolic{AbstractOperator})
     # do this if ispadded() = true
     #state, res = express_qc_proj(basis)(state, stateindex)
     # do this if ispadded() = false
@@ -12,14 +13,14 @@ express_qc_proj(::XGate) = QuantumClifford.projectX! # should be project*rand! i
 express_qc_proj(::YGate) = QuantumClifford.projectY! # should be project*rand! if ispadded()=true
 express_qc_proj(::ZGate) = QuantumClifford.projectZ! # should be project*rand! if ispadded()=true
 
-function observable(state::QuantumClifford.MixedDestabilizer, indices, operation)
+function observable(state::QuantumClifford.MixedDestabilizer, indices::Base.AbstractVecOrTuple{Int}, operation)
     operation = express(operation, CliffordRepr(), UseAsObservable())
     op = embed(QuantumClifford.nqubits(state), indices, operation)
     QuantumClifford.expect(op, state)
 end
 
 # This is a bit of a hack to work specifically with SProjector. If you start needing more of these for other types, consider doing a bit of a redesign. This all should pass through `express(...,::UseAsObservable)`.
-function observable(state::QuantumClifford.MixedDestabilizer, indices, operation::SProjector)
+function observable(state::QuantumClifford.MixedDestabilizer, indices::Base.AbstractVecOrTuple{Int}, operation::SProjector)
     pstate = express(operation.ket, CliffordRepr())
     QuantumClifford.nqubits(state)==length(indices)==QuantumClifford.nqubits(pstate) || error("An attempt was made to measure a projection observable while using Clifford representation for the qubits. However, the qubits that are being observed are entangled with other qubits. Currently this is not supported. Consider tracing out the extra qubits or using Pauli observables that do not suffer from this embedding limitation.")
     dot(pstate, state)

src/backends/clifford/should_upstream.jl

@@ -1,4 +1,3 @@
-QuantumClifford.apply!(state::QuantumClifford.MixedDestabilizer, op::Type{<:QuantumClifford.AbstractSymbolicOperator}, indices) = QuantumClifford.apply!(state, op(indices...)) # TODO piracy to be moved to QuantumClifford
 
 struct QCGateSequence <: QuantumClifford.AbstractSymbolicOperator
     gates # TODO constructor that flattens nested QCGateSequence

src/backends/quantumoptics/express.jl

@@ -4,10 +4,9 @@ _overlap(l::Symbolic{AbstractKet}, r::Operator) = _overlap(express(l, QOR), r)
 _project_and_drop(state::Ket, project_on::Symbolic{AbstractKet}, basis_index) = _project_and_drop(state, express(project_on, QOR), basis_index)
 _project_and_drop(state::Operator, project_on::Symbolic{AbstractKet}, basis_index) = _project_and_drop(state, express(project_on, QOR), basis_index)
 
-function project_traceout!(state::Union{Ket,Operator},stateindex,basis::Symbolic{AbstractOperator})
+function project_traceout!(state::Union{Ket,Operator},stateindex::Int,basis::Symbolic{AbstractOperator})
     project_traceout!(state,stateindex,eigvecs(basis))
 end
-
-function project_traceout!(state::Union{Ket,Operator},stateindex,basis::Base.AbstractVecOrTuple{<:Symbolic{AbstractKet}})
+function project_traceout!(state::Union{Ket,Operator},stateindex::Int,basis::Base.AbstractVecOrTuple{<:Symbolic{AbstractKet}})
     project_traceout!(state,stateindex,express.(basis,(QOR,)))
 end

src/backends/quantumoptics/quantumoptics.jl

@@ -17,21 +17,17 @@ subsystemcompose(op::Operator, k::Ket) = tensor(op,dm(k)) # TODO this should be
 default_repr(::StateVector) = QOR
 default_repr(::Operator) = QOR
 
-traceout!(s::StateVector, i) = ptrace(s,i)
-traceout!(s::Operator, i) = ptrace(s,i)
-
 ispadded(::StateVector) = false
 ispadded(::Operator) = false
 
-function observable(state::Union{<:Ket,<:Operator}, indices, operation)
+function observable(state::Union{<:Ket,<:Operator}, indices::Base.AbstractVecOrTuple{Int}, operation)
     operation = express(operation, QOR)
     e = basis(state)==basis(operation)
     op = e ? operation : embed(basis(state), indices, operation)
     expect(op, state)
 end
 
-
-function project_traceout!(state::Union{Ket,Operator},stateindex,psis::Base.AbstractVecOrTuple{<:Ket})
+function project_traceout!(state::Union{Ket,Operator},stateindex::Int,psis::Base.AbstractVecOrTuple{Ket})
     if nsubsystems(state) == 1 # TODO is there a way to do this in a single function, instead of _overlap vs _project_and_drop
         _overlaps = [_overlap(psi,state) for psi in psis]
         branch_probs = cumsum(_overlaps)

src/baseops/apply.jl

@@ -7,7 +7,7 @@ The appropriate representation of the gate is used,
 depending on the formalism under which a quantum state is stored in the given registers.
 The Hilbert spaces of the registers are automatically joined if necessary.
 """
-function apply!(regs::Vector{Register}, indices::Vector{Int}, operation; time=nothing)
+function apply!(regs::Vector{Register}, indices::Base.AbstractVecOrTuple{Int}, operation; time=nothing)
     max_time = maximum((r.accesstimes[i] for (r,i) in zip(regs,indices)))
     if !isnothing(time)
         time<max_time && error("The simulation was commanded to apply $(operation) at time t=$(time) although the current simulation time is higher at t=$(max_time). Consider using locks around the offending operations.")

src/baseops/initialize.jl

@@ -15,7 +15,7 @@ set the state of the given slots in the given registers to `state`.
 e.g., kets or density matrices from `QuantumOptics.jl`
 or tableaux from `QuantumClifford.jl`.
 """
-function initialize!(regs::Vector{Register},indices::Vector{Int},state; time=nothing)
+function initialize!(regs::Vector{Register},indices::Base.AbstractVecOrTuple{Int},state; time=nothing)
     length(regs)==length(indices)==nsubsystems(state) || throw(DimensionMismatch(lazy"Attempting to initialize a set of registers with a state that does not have the correct number of subsystems."))
     stateref = StateRef(state, collect(regs), collect(indices))
     for (si,(reg,ri)) in enumerate(zip(regs,indices))
@@ -34,7 +34,7 @@ initialize!(reg::Register,i::Int,state; time=nothing) = initialize!([reg],[i],st
 initialize!(r::RegRef, state; time=nothing) = initialize!(r.reg, r.idx, state; time)
 initialize!(r::Vector{Register},i::Vector{Int},state::Symbolic; time=nothing) = initialize!(r,i,express(state,consistent_representation(r,i,state)); time)
 
-function QuantumSymbolics.consistent_representation(regs,idx,state)
+function QuantumSymbolics.consistent_representation(regs::Base.AbstractVecOrTuple{Register},idx,state)
     reprs = Set([r.reprs[i] for (r,i) in zip(regs,idx)])
     consistent_representation(reprs,state)
 end

src/baseops/observable.jl

@@ -5,7 +5,7 @@ Calculate the expectation value of a quantum observable on the given register an
 of the `obs` observable (using the appropriate formalism, depending on the state
 representation in the given registers).
 """
-function observable(regs::Base.AbstractVecOrTuple{Register}, indices::Base.AbstractVecOrTuple{Int}, obs, something=nothing; time=nothing) # TODO weird split between positional and keyword arguments
+function observable(regs::Base.AbstractVecOrTuple{Register}, indices::Base.AbstractVecOrTuple{Int}, obs; something=nothing, time=nothing)
     staterefs = [r.staterefs[i] for (r,i) in zip(regs, indices)]
     # TODO it should still work even if they are not represented in the same state
     (any(isnothing, staterefs) || !all(s->s===staterefs[1], staterefs)) && return something
@@ -14,5 +14,5 @@ function observable(regs::Base.AbstractVecOrTuple{Register}, indices::Base.Abstr
     state_indices = [r.stateindices[i] for (r,i) in zip(regs, indices)]
     observable(state, state_indices, obs)
 end
-observable(refs::Base.AbstractVecOrTuple{RegRef}, obs, something=nothing; time=nothing) = observable(map(r->r.reg, refs), map(r->r.idx, refs), obs, something; time)
-observable(ref::RegRef, obs, something=nothing; time=nothing) = observable([ref.reg], [ref.idx], obs, something; time)
+observable(refs::Base.AbstractVecOrTuple{RegRef}, obs; something=nothing, time=nothing) = observable(map(r->r.reg, refs), map(r->r.idx, refs), obs; something, time)
+observable(ref::RegRef, obs; something=nothing, time=nothing) = observable([ref.reg], [ref.idx], obs; something, time)

src/baseops/subsystemcompose.jl

@@ -2,7 +2,6 @@ nsubsystems(s::StateRef) = length(s.registers) # nsubsystems(s.state[]) TODO thi
 nsubsystems_padded(s::StateRef) = nsubsystems(s.state[])
 nsubsystems(r::Register) = length(r.staterefs)
 nsubsystems(r::RegRef) = 1
-nsubsystems(::Nothing) = 1 # TODO consider removing this and reworking the functions that depend on it. E.g., a reason to have it when performing a project_traceout measurement on a state that contains only one subsystem
 
 function swap!(reg1::Register, reg2::Register, i1::Int, i2::Int; time=nothing)
     if reg1===reg2 && i1==i2

src/noninstant.jl

@@ -6,7 +6,7 @@ struct NonInstantGate <: AbstractNoninstantOperation
     duration # TODO assert larger than zero
 end
 
-function apply!(regs::Vector{Register}, indices::Vector{Int}, operation::NonInstantGate; time=nothing)
+function apply!(regs::Vector{Register}, indices::Base.AbstractVecOrTuple{Int}, operation::NonInstantGate; time=nothing)
     _, new_time = apply!(regs, indices, operation.gate; time)
     uptotime!(regs, indices, new_time+operation.duration)
     regs, new_time+operation.duration
@@ -18,7 +18,7 @@ struct ConstantHamiltonianEvolution <: AbstractNoninstantOperation
     duration # TODO assert larger than zero
 end
 
-function apply!(regs::Vector{Register}, indices::Vector{Int}, operation::ConstantHamiltonianEvolution; time=nothing) # TODO very significant code repetition with the general purpose apply!
+function apply!(regs::Vector{Register}, indices::Base.AbstractVecOrTuple{Int}, operation::ConstantHamiltonianEvolution; time=nothing) # TODO very significant code repetition with the general purpose apply!
     max_time = maximum((r.accesstimes[i] for (r,i) in zip(regs,indices)))
     if !isnothing(time)
         time<max_time && error("The simulation was commanded to apply $(operation) at time t=$(time) although the current simulation time is higher at t=$(max_time). Consider using locks around the offending operations.")

src/quantumchannel.jl

@@ -46,10 +46,11 @@ struct QuantumChannel{T}
     trait::T
     queue::ConcurrentSim.DelayQueue{Register}
     background::Any
+    function QuantumChannel(queue::ConcurrentSim.DelayQueue{Register}, background=nothing, trait::T=Qubit()) where T
+        new{T}(trait, queue, background)
+    end
 end
 
-QuantumChannel(queue::ConcurrentSim.DelayQueue{Register}, background=nothing, trait=Qubit()) = QuantumChannel(trait, queue, background)
-
 QuantumChannel(env::ConcurrentSim.Simulation, delay, background=nothing, trait=Qubit()) = QuantumChannel(ConcurrentSim.DelayQueue{Register}(env, delay), background, trait)
 Register(qc::QuantumChannel) = Register([qc.trait], [qc.background])