diff --git a/CHANGELOG.md b/CHANGELOG.md index b3d6dc762..91ae0b046 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/Project.toml b/Project.toml index fec5ef238..36ab133f7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumClifford" uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1" authors = ["Stefan Krastanov "] -version = "0.5.4" +version = "0.5.5" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 000000000..9005636ab --- /dev/null +++ b/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" diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index caf3c01bf..c41b1a739 100644 --- a/src/QuantumClifford.jl +++ b/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|(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 diff --git a/src/linalg.jl b/src/linalg.jl index 588cd5309..79c73af47 100644 --- a/src/linalg.jl +++ b/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 diff --git a/src/project_trace_reset.jl b/src/project_trace_reset.jl index b18e3bca3..7c9dcd1ef 100644 --- a/src/project_trace_reset.jl +++ b/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 diff --git a/test/test_cliff.jl b/test/test_cliff.jl index dcbbc169f..d9b0627ad 100644 --- a/test/test_cliff.jl +++ b/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 diff --git a/test/test_hash.jl b/test/test_hash.jl index ac805c459..e70e38192 100644 --- a/test/test_hash.jl +++ b/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 diff --git a/test/test_jet.jl b/test/test_jet.jl index f933b0224..c30977610 100644 --- a/test/test_jet.jl +++ b/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 diff --git a/test/test_noisycircuits.jl b/test/test_noisycircuits.jl index eb0bee73d..827580f16 100644 --- a/test/test_noisycircuits.jl +++ b/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() \ No newline at end of file +test_noisycircuits()