bug in observable when nsubsystems(state)==nsubsystems(op) but the indices are disordered (#212)

Author: Krastanov

.typos.toml

@@ -1,5 +1,6 @@
 [default.extend-words]
 ket = "ket"
+BA = "BA"
 
 [files]
 extend-exclude = ["src/StatesZoo/zalm/*.nb"]

Project.toml

@@ -46,7 +46,7 @@ Makie = "0.21, 0.22"
 NetworkLayout = "0.4.4"
 PrecompileTools = "1"
 Printf = "1"
-QuantumClifford = "0.9.9"
+QuantumClifford = "0.9.19"
 QuantumInterface = "0.3.5"
 QuantumOptics = "1.1.0"
 QuantumOpticsBase = "0.5.3"

src/QuantumSavory.jl

@@ -26,8 +26,20 @@ using SumTypes: @sum_type, isvariant, @cases
 import Combinatorics
 using Combinatorics: powerset
 
-import QuantumInterface: basis, tensor, ⊗, apply!, traceout!, nsubsystems,
-    AbstractOperator, AbstractKet, AbstractSuperOperator, Basis, SpinBasis
+import QuantumClifford
+import QuantumClifford: MixedDestabilizer
+
+import QuantumOpticsBase
+import QuantumOpticsBase: StateVector, Ket, Operator,
+    basisstate, spinup, spindown, sigmap, sigmax, sigmay, sigmaz, destroy, spre, spost
+
+import QuantumOptics
+import QuantumOptics: timeevolution
+
+import QuantumInterface: basis, tensor, ⊗, apply!, traceout!, nsubsystems, permutesystems,
+    projector, identityoperator, embed, dm, expect, ptrace,
+    AbstractOperator, AbstractKet, AbstractSuperOperator,
+    Basis, GenericBasis, CompositeBasis, SpinBasis
 
 export apply!, traceout!, removebackref!, nsubsystems
 export project_traceout! #TODO should move to QuantumInterface

src/backends/clifford/clifford.jl

@@ -1,6 +1,3 @@
-import QuantumClifford
-import QuantumClifford: MixedDestabilizer
-
 subsystemcompose(states::QuantumClifford.MixedDestabilizer...) = QuantumClifford.tensor(states...)
 
 default_repr(::QuantumClifford.MixedDestabilizer) = CliffordRepr()

src/backends/clifford/express.jl

@@ -22,6 +22,7 @@ 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::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)
+    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. Message us on the issue tracker if you want this functionality implemented.")
+    cstate = permutesystems(state, indices)
+    dot(pstate, cstate)
 end

src/backends/quantumoptics/quantumoptics.jl

@@ -1,12 +1,3 @@
-import QuantumOpticsBase
-import QuantumOpticsBase: GenericBasis, CompositeBasis,
-    StateVector, AbstractSuperOperator, Ket, Operator,
-    basisstate, spinup, spindown, sigmap, sigmax, sigmay, sigmaz, destroy,
-    projector, identityoperator, embed, dm, expect, ptrace, spre, spost
-import QuantumOptics
-import QuantumOptics: timeevolution
-import QuantumInterface: nsubsystems
-
 const QOR = QuantumOpticsRepr()
 
 subsystemcompose(states::Ket...) = tensor(states...)
@@ -22,8 +13,9 @@ ispadded(::Operator) = false
 
 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)
+    # TODO if indices is ascending 1:n we can skip this embed -- such an improvement should be upstreamed to QuantumOpticsBase, so that embed is faster
+    # TODO if nsubsystems(state) == 1 the embed should still work and be a no-op -- this should be upstreamed to QuantumInterface
+    op = nsubsystems(state) == 1 ? operation : embed(basis(state), indices, operation)
     expect(op, state)
 end
 

test/test_observable.jl

@@ -1,17 +1,50 @@
-@testitem "Observable" tags=[:observable] begin
+@testitem "Observable" begin
 
-const bell = StabilizerState("XX ZZ")
-# or equivalently `const bell = (Z₁⊗Z₁+Z₂⊗Z₂)/√2`,
-# however converting to stabilizer state for Clifford simulations
-# is not implemented (and can not be done efficiently).
+@testset "entangled observable" begin
+    bell = StabilizerState("XX ZZ")
+    # or equivalently `const bell = (Z₁⊗Z₁+Z₂⊗Z₂)/√2`,
+    # however converting to stabilizer state for Clifford simulations
+    # is not implemented (and can not be done efficiently).
 
-for rep in [QuantumOpticsRepr(), CliffordRepr()]
-    a = Register(2,rep)
-    initialize!(a[1:2], bell)
-    @test observable(a[1:2], SProjector(bell)) ≈ 1.0
-    @test observable(a[1:2], σˣ⊗σˣ) ≈ 1.0
-    apply!(a[1], σʸ)
-    @test observable(a[1:2], SProjector(bell)) ≈ 0.0
-    @test observable(a[1:2], σˣ⊗σˣ) ≈ -1.0
+    for rep in [QuantumOpticsRepr(), CliffordRepr()]
+        a = Register(2,rep)
+        initialize!(a[1:2], bell)
+        @test observable(a[1:2], SProjector(bell)) ≈ 1.0
+        @test observable(a[1:2], σˣ⊗σˣ) ≈ 1.0
+        apply!(a[1], σʸ)
+        @test observable(a[1:2], SProjector(bell)) ≈ 0.0
+        @test observable(a[1:2], σˣ⊗σˣ) ≈ -1.0
+    end
 end
+
+@testset "separable observable with order flipping" begin
+    A = StabilizerState("X")
+    B = StabilizerState("Z")
+    AB = StabilizerState("XI IZ")
+    BA = StabilizerState("ZI IX")
+
+    for rep in [QuantumOpticsRepr(), CliffordRepr()]
+        r1 = Register(2,rep)
+        r2 = Register(2,rep)
+        r12 = Register(2,rep)
+        r21 = Register(2,rep)
+
+        initialize!(r1[1], A)
+        initialize!(r2[2], B)
+        initialize!(r12[1:2], AB)
+        initialize!((r21[2], r21[1]), BA)
+
+        @test observable(r12[1:2], SProjector(AB)) ≈ 1.0
+        if rep == CliffordRepr()
+            @test_throws "entangled with other qubits" observable(r12[1], SProjector(A)) ≈ 1.0
+        else
+            @test observable(r12[1], SProjector(A)) ≈ 1.0
+        end
+        @test observable(r21[1:2], SProjector(AB)) ≈ 1.0
+        @test_broken observable((r1[1], r2[2]), SProjector(AB)) ≈ 1.0
+        @test observable(r1[1], SProjector(A)) ≈ 1.0
+        @test observable(r2[2], SProjector(B)) ≈ 1.0
+    end
+end
+
 end