JET related fixes (#63)

* Fix boxed variable in closure (a.k.a. THE Julia closure problem)

* Add Project.toml to benchmarks

* Relax statistical test for noisy circuits

* Destabilizer(view) bug fix and zero! JET fix

* Forbid non-square tableaux in CliffordOperator constructor

* Improve _project! on (Mixed)Destabilizer as guided by JET

* version bump

CHANGELOG.md

@@ -1,5 +1,11 @@
 # News
 
+## v0.5.5 - 2022-07-05
+
+- **(breaking fix)** `CliffordOperator` constructor called on a square tableau occasionally resulted in incorrect results due to row-reordering during cannonicalization.
+- Continuing static analysis fixes thanks to JET.
+- Optimization of `canonicalize_clip!`, now resulting in much fewer allocations.
+
 ## v0.5.4 - 2022-07-03
 
 - Start using `JET.jl` for static analysis during CI.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]>"]
-version = "0.5.4"
+version = "0.5.5"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

benchmark/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+QuantumClifford = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"

src/QuantumClifford.jl

@@ -204,8 +204,8 @@ Base.zero(::Type{<:PauliOperator}, q) = PauliOperator(zeros(UInt8), q, zeros(UIn
 Base.zero(p::PauliOperator) = zero(PauliOperator, nqubits(p))
 
 """Zero-out the phases and single-qubit operators in a [`PauliOperator`](@ref)"""
-@inline function zero!(p::PauliOperator)
-    fill!(p.xz, zero(eltype(p.xz)))
+@inline function zero!(p::PauliOperator{Tz,Tv}) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}}
+    fill!(p.xz, zero(Tve))
     p.phase[] = 0x0
     p
 end
@@ -350,9 +350,9 @@ Base.getindex(stab::Stabilizer, i::Int) = PauliOperator(stab.phases[i], nqubits(
      return (x,z)
 end
 Base.getindex(stab::Stabilizer, r) = Stabilizer(stab.phases[r], nqubits(stab), stab.xzs[:,r])
-Base.getindex(stab::Stabilizer, r, c) = Stabilizer([s[c] for s in stab[r]])
-Base.getindex(stab::Stabilizer, r, c::Int) = stab[r,[c]]
-Base.getindex(stab::Stabilizer, r::Int, c) = stab[r][c]
+Base.getindex(stab::Stabilizer, r::Union{Colon,AbstractVector}, c::Union{Colon,AbstractVector}) = Stabilizer([s[c] for s in stab[r]])
+Base.getindex(stab::Stabilizer, r::Union{Colon,AbstractVector}, c::Int) = stab[r,[c]]
+Base.getindex(stab::Stabilizer, r::Int, c::Union{Colon,AbstractVector}) = stab[r][c]
 Base.view(stab::Stabilizer, r) = Stabilizer(view(stab.phases, r), nqubits(stab), view(stab.xzs, :, r))
 
 Base.iterate(stab::Stabilizer, state=1) = state>length(stab) ? nothing : (stab[state], state+1)
@@ -491,7 +491,7 @@ struct Destabilizer{Tzv<:AbstractVector{UInt8},Tm<:AbstractMatrix{<:Unsigned}} <
         elseif r<=n
             mixed_destab = MixedDestabilizer(s)
             tab = vcat(destabilizerview(mixed_destab),stabilizerview(mixed_destab))
-            new{typeof(s.phases),typeof(s.xzs)}(tab)
+            new{typeof(tab.phases),typeof(tab.xzs)}(tab)
         else
             throw(DomainError("To construct a `Destabilizer`, either use a 2n×n tableau of destabilizer and stabilizer rows that is directly used or an r×n (r<n) tableau of stabilizers from which destabilizers are automatically computed. Or better, just use `MixedDestabilizer`."))
         end

src/dense_cliffords.jl

@@ -51,11 +51,11 @@ struct CliffordOperator{Tzv<:AbstractVector{UInt8},Tm<:AbstractMatrix{<:Unsigned
     function CliffordOperator(stab::Stabilizer{Tzv,Tm}) where {Tzv,Tm}
         if size(stab,1)==2*size(stab,2)
             new{Tzv,Tm}(stab)
-        elseif size(stab,1)==size(stab,2)
-            destab = tab(Destabilizer(stab))
-            new{typeof(destab.phases),typeof(destab.xzs)}(destab) # TODO be smarter about type signatures here... there should be a better way
+        #elseif size(stab,1)==size(stab,2) # TODO be able to work with squara tableaux (by reversing all row operations)
+        #    destab = tab(Destabilizer(stab))
+        #    new{typeof(destab.phases),typeof(destab.xzs)}(destab) # TODO be smarter about type signatures here... there should be a better way
         else
-            throw(DimensionMismatch("Input tableau should be square (in which case the destabilizers are calculated) or of size 2n×n (in which case it is used directly)."))
+            throw(DimensionMismatch("Input tableau should be of size 2n×n (top half is the X mappings and the bottom half are the Z mappings)."))
         end
     end
 end

src/entanglement.jl

@@ -34,22 +34,21 @@ function canonicalize_clip!(state::AbstractStabilizer; phases::Bool=true)
     _canonicalize_clip!(state; phases=Val(phases))
 end
 function  _canonicalize_clip!(state::AbstractStabilizer; phases::Val{B}=Val(true)) where B
-    # TODO: the function has spurious allocation, to be further optimized
     tab = stabilizerview(state)
     rows, columns = size(stabilizerview(state))
     # step 1: pregauge
     i = 1 # index to place used stab
     for j in 1:columns
         # find first row that is not I in col j
-        k1 = findfirst(k->|(tab[k,j]...), i:rows)
+        k1 = findfirst(let j=j; k->|(tab[k,j]...) end, i:rows)
         # find second row that is not I and not same as k1
-        if k1!==nothing
+        if !isnothing(k1)
             k1 += i-1
-            k2 = findfirst(k->
-                (|(tab[k,j]...) & # not identity
-                ((tab[k,j][1]⊻tab[k1,j][1]) | (tab[k,j][2]⊻tab[k1,j][2]))), # not same as k1
-                k1+1:columns)
-            if k2!==nothing
+            k2 = findfirst(let j=j, k1=k1; k->
+                           (|(tab[k,j]...) & # not identity
+                           (tab[k,j]!=tab[k1,j])) end, # not same as k1
+                           k1+1:columns)
+            if !isnothing(k2)
                 k2 += k1
                 # move k1 and k2 up to i and i+1
                 rowswap!(state, k1, i; phases=phases)
@@ -83,27 +82,27 @@ function  _canonicalize_clip!(state::AbstractStabilizer; phases::Val{B}=Val(true
     unfrozen_rows = Array(rows:-1:1)
     for j in columns:-1:1 # in reversed order to keep left ends
         # find first row that is not I in col j
-        k1 = findfirst(k->|(tab[k,j]...), unfrozen_rows)
+        k1 = findfirst(let j=j; k->|(tab[k,j]...) end, unfrozen_rows)
         # find second row that is not I and not same as k1
         if k1!==nothing
             k1_row = unfrozen_rows[k1]
-            k2 = findfirst(k->
-                (|(tab[k,j]...) & # not identity
-                ((tab[k,j][1]⊻tab[k1_row,j][1]) | (tab[k,j][2]⊻tab[k1_row,j][2]))), # not same as k1
-                unfrozen_rows[k1+1:end])
+            k2 = findfirst(let j=j, k1_row=k1_row; k->
+                           (|(tab[k,j]...) & # not identity
+                           (tab[k,j]!=tab[k1_row,j])) end, # not same as k1
+                           @view unfrozen_rows[k1+1:end])
 
             if k2!==nothing
                 k2 += k1
                 k2_row = unfrozen_rows[k2]
                 # use them to eliminate others
                 # for rows between k1 and k2, use k1
-                for m in unfrozen_rows[k1+1:k2-1]
+                for m in @view unfrozen_rows[k1+1:k2-1]
                     if !(tab[m,j][1]⊻tab[k1_row,j][1]) && !(tab[m,j][2]⊻tab[k1_row,j][2])
                         mul_left!(state, m, k1_row; phases=phases)
                     end
                 end
                 # for other rows, use both
-                for m in unfrozen_rows[k2+1:end]
+                for m in @view unfrozen_rows[k2+1:end]
                     if !(tab[m,j][1]⊻tab[k1_row,j][1]) && !(tab[m,j][2]⊻tab[k1_row,j][2])
                         mul_left!(state, m, k1_row; phases=phases)
                     elseif !(tab[m,j][1]⊻tab[k2_row,j][1]) && !(tab[m,j][2]⊻tab[k2_row,j][2])
@@ -116,7 +115,7 @@ function  _canonicalize_clip!(state::AbstractStabilizer; phases::Val{B}=Val(true
                 deleteat!(unfrozen_rows, (k1, k2))
             else # can only find k1
                 # use it to eliminate others
-                for m in unfrozen_rows[k1+1:end]
+                for m in @view unfrozen_rows[k1+1:end]
                     if !(tab[m,j][1]⊻tab[k1_row,j][1]) && !(tab[m,j][2]⊻tab[k1_row,j][2])
                         mul_left!(state, m, k1_row; phases=phases)
                     end
@@ -149,8 +148,11 @@ function bigram(state::AbstractStabilizer; clip::Bool=true)
     rows, columns = size(tab)
     bg = zeros(Int, rows, 2)
     for i in 1:rows
-        bg[i, 1] = findfirst(j->|(tab[i,j]...), 1:columns)
-        bg[i, 2] = findlast(j->|(tab[i,j]...), 1:columns)
+        l = findfirst(j->|(tab[i,j]...), 1:columns)
+        r = findlast(j->|(tab[i,j]...), 1:columns)
+        (isnothing(l) || isnothing(r)) && throw(DomainError("the tableau is inconsistent (check if it is clip-canonicalized and Hermitian)"))
+        bg[i, 1] = l
+        bg[i, 2] = r
     end
     bg
 end

src/linalg.jl

@@ -42,7 +42,7 @@ Otherwise the inner product is `2^(-logdot/2)`.
 The actual inner product can be computed with `LinearAlgebra.dot`.
 
 Based on [garcia2012efficient](@cite)."""
-function logdot(s1::AbstractStabilizer, s2::AbstractStabilizer) # TODO verify rank
+function logdot(s1::AbstractStabilizer, s2::AbstractStabilizer) # TODO verify rank # TODO this is currently very inefficient as we discard the destabilizers and then recreate them
     logdot(stabilizerview(s1),stabilizerview(s2))
 end
 
@@ -53,7 +53,7 @@ function logdot(s1::Stabilizer, s2::Stabilizer)
     if nqubits(s1)!=nqubits(s2)
         throw(DimensionMismatch("Inner product can be calculated only between states with the same number of qubits."))
     end
-    c1_inv = inv(CliffordOperator(copy(s1)))
+    c1_inv = inv(CliffordOperator(tab(MixedDestabilizer(copy(s1)))))
     s2_prime = canonicalize!(c1_inv*s2)
     canonicalize!(s2_prime)
     k = 0

src/project_trace_reset.jl

@@ -338,7 +338,7 @@ end
     end
 end
 
-function _project!(d::Union{Destabilizer,MixedDestabilizer},pauli::PauliOperator;keep_result::Val{Bkr}=Val(true),phases::Val{Bp}=Val(true)) where {Bkr, Bp}
+function _project!(d::Destabilizer,pauli::PauliOperator;keep_result::Val{Bkr}=Val(true),phases::Val{Bp}=Val(true)) where {Bkr, Bp} # repetition between Destabilizer and MixedDestabilizer, but the redundancy makes the two codes slightly simpler and easier to infer
     anticommutes = 0
     tab = d.tab
     stabilizer = stabilizerview(d)
@@ -357,6 +357,39 @@ function _project!(d::Union{Destabilizer,MixedDestabilizer},pauli::PauliOperator
                                    :project!,
                                    :Destabilizer))
         end
+        if Bkr
+            new_pauli = zero(pauli)
+            new_pauli.phase[] = pauli.phase[]
+            for i in 1:r
+                comm(pauli,destabilizer,i)!=0 && mul_left!(new_pauli, stabilizer, i, phases=phases)
+            end
+            result = new_pauli.phase[]
+        else
+            result = nothing
+        end
+    else
+        anticomm_update_rows(tab,pauli,r,n,anticommutes,phases)
+        destabilizer[anticommutes] = stabilizer[anticommutes]
+        stabilizer[anticommutes] = pauli
+        result = nothing
+    end
+    d, anticommutes, result
+end
+
+function _project!(d::MixedDestabilizer,pauli::PauliOperator;keep_result::Val{Bkr}=Val(true),phases::Val{Bp}=Val(true)) where {Bkr, Bp} # repetition between Destabilizer and MixedDestabilizer, but the redundancy makes the two codes slightly simpler and easier to infer
+    anticommutes = 0
+    tab = d.tab
+    stabilizer = stabilizerview(d)
+    destabilizer = destabilizerview(d)
+    r = trusted_rank(d)
+    n = length(d) # not `nqubits(d)` in case we have an incomplete tableau    # Check whether we anticommute with any of the stabilizer rows
+    for i in 1:r # The explicit loop is faster than anticommutes = findfirst(row->comm(pauli,stabilizer,row)!=0x0, 1:r); both do not allocate.
+        if comm(pauli,stabilizer,i)!=0x0
+            anticommutes = i
+            break
+        end
+    end
+    if anticommutes == 0
         anticomlog = 0
         # Check whether we anticommute with any of the logical X rows
         for i in r+1:n # The explicit loop is faster than findfirst.
@@ -383,7 +416,7 @@ function _project!(d::Union{Destabilizer,MixedDestabilizer},pauli::PauliOperator
                 rowswap!(tab, r+1+n, anticomlog)
             end
             anticomm_update_rows(tab,pauli,r+1,n,r+1,phases)
-            d.rank += 1
+            d.rank+=1
             anticommutes = d.rank
             tab[r+1] = tab[n+r+1]
             tab[n+r+1] = pauli

test/test_cliff.jl

@@ -14,6 +14,9 @@ function test_cliff()
         end
     end
     @testset "Clifford Operators" begin
+        @testset "Constructors" begin
+            @test_throws DimensionMismatch CliffordOperator(S"X")
+        end
         @testset "Permutations of qubits" begin
             for c in [tCNOT, tId1⊗tHadamard, tCNOT⊗tCNOT, tensor_pow(tCNOT,6), tensor_pow(tCNOT,7), tensor_pow(tCNOT,6)⊗tPhase, tensor_pow(tCNOT,7)⊗tPhase]
                 for rep in 1:5

test/test_hash.jl

@@ -2,7 +2,7 @@ function test_hash()
     @testset "Hashing" begin
         @test hash(P"X") == hash(P"X")
         @test hash(S"X") == hash(S"X")
-        @test hash(C"X") == hash(C"X")
+        @test hash(C"X Z") == hash(C"X Z")
     end
 end
 

test/test_jet.jl

@@ -33,7 +33,7 @@ function test_jet()
         @test isempty(JET.get_reports(@report_call QuantumClifford._canonicalize_gott!(s)))
         @test isempty(JET.get_reports(@report_call QuantumClifford._canonicalize_rref!(s,[1,3])))
 
-        @test_broken isempty(JET.get_reports(report_package("QuantumClifford")))
+        @test_broken isempty(JET.get_reports(report_package("QuantumClifford"; ignored_modules=(AnyFrameModule(Graphs.LinAlg),))))
     end
 end
 

test/test_noisycircuits.jl

@@ -17,7 +17,7 @@ function test_noisycircuits()
                 with_purification = mctrajectories(init, [n,g1,g2,m,v], trajectories=500)
                 @test with_purification[failure_stat] > 5
                 @test with_purification[false_success_stat] > 10
-                @test with_purification[true_success_stat] > 430
+                @test with_purification[true_success_stat] > 420
                 without_purification = mctrajectories(init, [n,v], trajectories=500)
                 @test get(without_purification,failure_stat,0) == 0
                 @test without_purification[false_success_stat] > 10
@@ -57,12 +57,12 @@ function test_noisycircuits()
                 @test compare(mc,pe,false_success_stat)
                 @test compare(mc,pe,true_success_stat)
             end
-            
+
             @testset "Symbolic" begin
                 for statetype in [MixedDestabilizer]
                     R, (e,) = AbstractAlgebra.PolynomialRing(AbstractAlgebra.RealField, ["e"])
                     unity = R(1);
-                    
+
                     good_bell_state = statetype(S"XX
                                                   ZZ")
                     initial_state = good_bell_state⊗good_bell_state
@@ -296,4 +296,4 @@ function test_noisycircuits()
     end
 end
 
-test_noisycircuits()
+test_noisycircuits()