diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e00b68835..1f7eba665 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -4,6 +4,14 @@ on: branches: [master, main] tags: ["*"] pull_request: + +concurrency: + # group by workflow and ref; the last slightly strange component ensures that for pull + # requests, we limit to 1 concurrent job, but for the master branch we don't + group: ${{ github.workflow }}-${{ github.ref }}-${{ github.ref != 'refs/heads/master' || github.run_number }} + # Cancel intermediate builds, but only if it is a pull request build. + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + env: PYTHON: ~ jobs: diff --git a/.gitignore b/.gitignore index 8bf2daafc..b159e1a46 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,6 @@ LocalPreferences.toml */.*swp scratch/ *.cov -.vscode \ No newline at end of file +.vscode +test/.CondaPkg/ +docs/.CondaPkg/ diff --git a/CHANGELOG.md b/CHANGELOG.md index e260d5c02..7bf3a570d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,22 @@ # News +## v0.9.14 - 2024-11-03 + +- **(fix)** `affectedqubits()` on `sMX`, `sMY`, and `sMR*` +- **(fix)** restrictive type-assert in `MixedDestabilizer` failing on views of tableaux +- Implementing additional named two-qubit gates: `sSQRTXX, sInvSQRTXX, sSQRTYY, sInvSQRTYY` + +## v0.9.13 - 2024-10-30 + +- New error-correction group theory tools: + - `canonicalize_noncomm` function to find a generating set with minimal anticommutivity + - `SubsystemCodeTableau` data structure to represent the output of `canonicalize_noncomm` + - `commutify` function to find a commutative version of a non-commutative set of Paulis with minimal changes + - `matroid_parent` to, for set of Paulis that doesn't represent a state, find a version + that does. +- Implementing additional named two-qubit gates: `sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP, sSQRTZZ, sInvSQRTZZ` + ## v0.9.12 - 2024-10-18 - Minor compat fixes for julia 1.11 in the handling of `hgp` diff --git a/Project.toml b/Project.toml index 4d300bc66..e49f22503 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumClifford" uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1" authors = ["Stefan Krastanov and QuantumSavory community members"] -version = "0.9.12" +version = "0.9.14" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -50,7 +50,7 @@ DocStringExtensions = "0.9" Graphs = "1.9" Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34" HostCPUFeatures = "0.1.6" -ILog2 = "0.2.3" +ILog2 = "0.2.3, 1, 2" InteractiveUtils = "1.9" LDPCDecoders = "0.3.1" LinearAlgebra = "1.9" diff --git a/docs/src/noisycircuits_mc.md b/docs/src/noisycircuits_mc.md index 37e900df4..d24585573 100644 --- a/docs/src/noisycircuits_mc.md +++ b/docs/src/noisycircuits_mc.md @@ -15,7 +15,7 @@ Import with `using QuantumClifford.Experimental.NoisyCircuits`. This module enables the simulation of noisy Clifford circuits through a Monte Carlo method where the same circuit is evaluated multiple times with random errors interspersed through it as prescribed by a given error model. -Below is an example of a purification circuit. We first prepare the circuit we desire to use, including a noise model. `Quantikz.jl` was is used to visualize the circuit. +Below is an example of a purification circuit. We first prepare the circuit we desire to use, including a noise model. `Quantikz.jl` is used to visualize the circuit. ```@example 1 using QuantumClifford # hide @@ -55,8 +55,8 @@ If you want to create a custom gate type (e.g. calling it `Operation`), you need The `Symbol` is the status of the operation. Predefined statuses are kept in the `registered_statuses` list, but you can add more. Be sure to expand this list if you want the trajectory simulators using your custom statuses to output all trajectories. -There is also [`applynoise!`](@ref) which is convenient wait to create a noise model that can then be plugged into the [`NoisyGate`](@ref) struct, +There is also [`applynoise!`](@ref) which is a convenient way to create a noise model that can then be plugged into the [`NoisyGate`](@ref) struct, letting you reuse the predefined perfect gates and measurements. However, you can also just make up your own noise operator simply by implementing [`applywstatus!`](@ref) for it. -You can also consult the [list of implemented operators](@ref noisycircuits_ops). \ No newline at end of file +You can also consult the [list of implemented operators](@ref noisycircuits_ops). diff --git a/docs/src/references.bib b/docs/src/references.bib index fdcc25706..002960a68 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -403,6 +403,31 @@ @inproceedings{brown2013short doi = {10.1109/ISIT.2013.6620245} } +@article{RevModPhys.87.307, + title = {Quantum error correction for quantum memories}, + author = {Terhal, Barbara M.}, + journal = {Rev. Mod. Phys.}, + volume = {87}, + issue = {2}, + pages = {307--346}, + numpages = {40}, + year = {2015}, + month = {Apr}, + publisher = {American Physical Society}, + doi = {10.1103/RevModPhys.87.307}, + url = {https://link.aps.org/doi/10.1103/RevModPhys.87.307} +} + +@misc{goodenough2024bipartiteentanglementnoisystabilizer, + title={Bipartite entanglement of noisy stabilizer states through the lens of stabilizer codes}, + author={Kenneth Goodenough and Aqil Sajjad and Eneet Kaur and Saikat Guha and Don Towsley}, + year={2024}, + eprint={2406.02427}, + archivePrefix={arXiv}, + primaryClass={quant-ph}, + url={https://arxiv.org/abs/2406.02427}, +} + @article{panteleev2021degenerate, title = {Degenerate {{Quantum LDPC Codes With Good Finite Length Performance}}}, author = {Panteleev, Pavel and Kalachev, Gleb}, @@ -488,10 +513,20 @@ @article{anderson2014fault publisher={APS} } +@article{lin2024quantum, + title={Quantum two-block group algebra codes}, + author={Lin, Hsiang-Ku and Pryadko, Leonid P}, + journal={Physical Review A}, + volume={109}, + number={2}, + pages={022407}, + year={2024}, + publisher={APS} +} + @article{wang2024coprime, title={Coprime Bivariate Bicycle Codes and their Properties}, author={Wang, Ming and Mueller, Frank}, journal={arXiv preprint arXiv:2408.10001}, year={2024} } - diff --git a/ext/QuantumCliffordGPUExt/apply_noise.jl b/ext/QuantumCliffordGPUExt/apply_noise.jl index 4aa99ed4b..201eb7c60 100644 --- a/ext/QuantumCliffordGPUExt/apply_noise.jl +++ b/ext/QuantumCliffordGPUExt/apply_noise.jl @@ -1,16 +1,11 @@ -using QuantumClifford: _div, _mod +using QuantumClifford: get_bitmask_idxs #according to https://github.com/JuliaGPU/CUDA.jl/blob/ac1bc29a118e7be56d9edb084a4dea4224c1d707/test/core/device/random.jl#L33 #CUDA.jl supports calling rand() inside kernel function applynoise!(frame::PauliFrameGPU{T},noise::UnbiasedUncorrelatedNoise,i::Int) where {T <: Unsigned} p = noise.p - lowbit = T(1) - ibig = _div(T,i-1)+1 - ismall = _mod(T,i-1) - ismallm = lowbit<<(ismall) - - stab = frame.frame - xzs = tab(stab).xzs + xzs = tab(frame.frame).xzs + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) rows = size(stab, 1) @run_cuda applynoise_kernel(xzs, p, ibig, ismallm, rows) rows diff --git a/ext/QuantumCliffordGPUExt/pauli_frames.jl b/ext/QuantumCliffordGPUExt/pauli_frames.jl index 34a0565fd..4d6e03654 100644 --- a/ext/QuantumCliffordGPUExt/pauli_frames.jl +++ b/ext/QuantumCliffordGPUExt/pauli_frames.jl @@ -1,3 +1,5 @@ +using QuantumClifford: get_bitmask_idxs + ############################## # sMZ ############################## @@ -21,10 +23,7 @@ function apply!(frame::PauliFrameGPU{T}, op::QuantumClifford.sMZ) where {T <: Un op.bit == 0 && return frame i = op.qubit xzs = frame.frame.tab.xzs - lowbit = T(1) - ibig = QuantumClifford._div(T,i-1)+1 - ismall = QuantumClifford._mod(T,i-1) - ismallm = lowbit<<(ismall) + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) (@run_cuda apply_sMZ_kernel!(xzs, frame.measurements, op, ibig, ismallm, length(frame)) length(frame)) return frame end @@ -55,10 +54,7 @@ end function apply!(frame::PauliFrameGPU{T}, op::QuantumClifford.sMRZ) where {T <: Unsigned} # TODO sMRX, sMRY i = op.qubit xzs = frame.frame.tab.xzs - lowbit = T(1) - ibig = QuantumClifford._div(T,i-1)+1 - ismall = QuantumClifford._mod(T,i-1) - ismallm = lowbit<<(ismall) + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) (@run_cuda apply_sMRZ_kernel!(xzs, frame.measurements, op, ibig, ismallm, length(frame)) length(frame)) return frame end diff --git a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl index 29e9de8ce..b8508ff66 100644 --- a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl +++ b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl @@ -5,14 +5,15 @@ using DocStringExtensions import QuantumClifford, LinearAlgebra import Hecke: Group, GroupElem, AdditiveGroup, AdditiveGroupElem, GroupAlgebra, GroupAlgebraElem, FqFieldElem, representation_matrix, dim, base_ring, - multiplication_table, coefficients, abelian_group, group_algebra + multiplication_table, coefficients, abelian_group, group_algebra, rand import Nemo import Nemo: characteristic, matrix_repr, GF, ZZ, lift import QuantumClifford.ECC: AbstractECC, CSS, ClassicalCode, hgp, code_k, code_n, code_s, iscss, parity_checks, parity_checks_x, parity_checks_z, parity_checks_xz, - two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes + two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes, check_repr_commutation_relation +include("util.jl") include("types.jl") include("lifted.jl") include("lifted_product.jl") diff --git a/ext/QuantumCliffordHeckeExt/lifted_product.jl b/ext/QuantumCliffordHeckeExt/lifted_product.jl index 69a10ccef..b4c174a08 100644 --- a/ext/QuantumCliffordHeckeExt/lifted_product.jl +++ b/ext/QuantumCliffordHeckeExt/lifted_product.jl @@ -182,12 +182,32 @@ code_s(c::LPCode) = size(c.repr(zero(c.GA)), 1) * (size(c.A, 1) * size(c.B, 1) + Two-block group algebra (2GBA) codes, which are a special case of lifted product codes from two group algebra elements `a` and `b`, used as `1x1` base matrices. +Here is an example of a [[56, 28, 2]] 2BGA code from Table 2 of [lin2024quantum](@cite) +with direct product of `C₄ x C₂`. + +```jldoctest +julia> import Hecke: group_algebra, GF, abelian_group, gens; + +julia> GA = group_algebra(GF(2), abelian_group([14,2])); + +julia> x = gens(GA)[1]; + +julia> s = gens(GA)[2]; + +julia> A = 1 + x^7 + +julia> B = 1 + x^7 + s + x^8 + s*x^7 + x + +julia> c = two_block_group_algebra_codes(A,B); + +julia> code_n(c), code_k(c) +(56, 28) +``` + See also: [`LPCode`](@ref), [`generalized_bicycle_codes`](@ref), [`bicycle_codes`](@ref) -""" # TODO doctest example +""" function two_block_group_algebra_codes(a::GroupAlgebraElem, b::GroupAlgebraElem) - A = reshape([a], (1, 1)) - B = reshape([b], (1, 1)) - LPCode(A, B) + LPCode([a;;], [b;;]) end """ diff --git a/ext/QuantumCliffordHeckeExt/util.jl b/ext/QuantumCliffordHeckeExt/util.jl new file mode 100644 index 000000000..d8fc65e85 --- /dev/null +++ b/ext/QuantumCliffordHeckeExt/util.jl @@ -0,0 +1,15 @@ +""" +Checks the commutation relation between the left and right representation matrices +for two randomly-sampled elements `a` and `b` in the group algebra `ℱ[G]` with a general group `G`. +It verifies the commutation relation that states, `L(a)·R(b) = R(b)·L(a)`. This +property shows that matrices from the left and right representation sets commute +with each other, which is an important property related to the CSS orthogonality +condition. +""" +function check_repr_commutation_relation(GA::GroupAlgebra) + a, b = rand(GA), rand(GA) + # Check commutation relation: L(a)R(b) = R(b)L(a) + L_a = representation_matrix(a) + R_b = representation_matrix(b, :right) + return L_a * R_b == R_b * L_a +end diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index 60e020440..6b90af7a8 100644 --- a/src/QuantumClifford.jl +++ b/src/QuantumClifford.jl @@ -26,7 +26,7 @@ export nqubits, stabilizerview, destabilizerview, logicalxview, logicalzview, phases, fastcolumn, fastrow, - bitview, quantumstate, tab, + bitview, quantumstate, tab, rank, BadDataStructure, affectedqubits, #TODO move to QuantumInterface? # GF2 @@ -51,7 +51,8 @@ export sHadamardXY, sHadamardYZ, sSQRTX, sInvSQRTX, sSQRTY, sInvSQRTY, sCXYZ, sCZYX, sCNOT, sCPHASE, sSWAP, sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZCX, sZCY, sZCZ, - sZCrY, sInvZCrY, + sZCrY, sInvZCrY, sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP, + sSQRTZZ, sInvSQRTZZ, sSQRTXX, sInvSQRTXX, sSQRTYY, sInvSQRTYY, # Misc Ops SparseGate, sMX, sMY, sMZ, PauliMeasurement, Reset, sMRX, sMRY, sMRZ, @@ -77,6 +78,7 @@ export graphstate, graphstate!, graph_gatesequence, graph_gate, # Group theory tools groupify, minimal_generating_set, pauligroup, normalizer, centralizer, contractor, delete_columns, + canonicalize_noncomm, commutify, matroid_parent, SubsystemCodeTableau, # Clipped Gauge canonicalize_clip!, bigram, entanglement_entropy, # mctrajectories @@ -173,7 +175,7 @@ Tableau(xzs::AbstractMatrix{Bool}) = Tableau(zeros(UInt8, size(xzs,1)), xzs[:,1: Tableau(t::Tableau) = t -function _T_str(a) # TODO this can be optimized by not creating intermediary PauliOperator objects +function _T_str(a::Union{String,SubString{String}}) # TODO this can be optimized by not creating intermediary PauliOperator objects paulis = [_P_str(strip(s.match)) for s in eachmatch(r"[+-]?\h*[i]?\h*[XYZI_]+", a)] Tableau(paulis) end @@ -229,18 +231,18 @@ end Base.firstindex(tab::Tableau) = 1 -Base.lastindex(tab::Tableau) = length(tab.phases) +Base.lastindex(tab::Tableau) = length(tab.phases)::Int Base.lastindex(tab::Tableau, i) = size(tab)[i] -Base.eachindex(tab::Tableau) = Base.OneTo(lastindex(tab.phases)) +Base.eachindex(tab::Tableau) = Base.OneTo(lastindex(tab.phases)::Int) Base.axes(tab::Tableau) = (Base.OneTo(lastindex(tab)), Base.OneTo(nqubits(tab))) Base.axes(tab::Tableau,i) = axes(tab)[i] -Base.size(tab::Tableau) = (length(tab.phases),nqubits(tab)) +Base.size(tab::Tableau) = (length(tab.phases)::Int, nqubits(tab)) Base.size(tab::Tableau,i) = size(tab)[i] -Base.length(tab::Tableau) = length(tab.phases) +Base.length(tab::Tableau) = length(tab.phases)::Int Base.:(==)(l::Tableau, r::Tableau) = r.nqubits==l.nqubits && r.phases==l.phases && r.xzs==l.xzs @@ -378,7 +380,7 @@ macro S_str(a) quote Stabilizer(_T_str($a)) end end Base.getindex(stab::Stabilizer, i::Int) = tab(stab)[i] -Base.getindex(stab::Stabilizer, i) = Stabilizer(tab(stab)[i]) +Base.getindex(stab::Stabilizer, i) = Stabilizer(tab(stab)[i]::Tableau) @inline Base.getindex(stab::Stabilizer, r::Int, c) = tab(stab)[r,c] Base.getindex(stab::Stabilizer, r, c) = Stabilizer(tab(stab)[r,c]) Base.view(stab::Stabilizer, r) = Stabilizer(view(tab(stab),r)) @@ -497,11 +499,11 @@ function MixedStabilizer(s::Stabilizer{T}) where {T} MixedStabilizer(spadded,zr) end -MixedStabilizer(s::Stabilizer,rank::Int) = MixedStabilizer(tab(s),rank) +MixedStabilizer(s::Stabilizer,rank::Int) = MixedStabilizer(tab(s), rank) -Base.length(d::MixedStabilizer) = length(d.tab) +Base.length(d::MixedStabilizer) = length(tab(d)) -Base.copy(ms::MixedStabilizer) = MixedStabilizer(copy(ms.tab),ms.rank) +Base.copy(ms::MixedStabilizer) = MixedStabilizer(copy(tab(ms)), rank(ms)) ############################## # Mixed Destabilizer states @@ -569,8 +571,8 @@ function MixedDestabilizer(stab::Stabilizer{T}; undoperm=true, reportperm=false) t[n+r+s+1:end] = sZ # The other logical set in the tableau end if undoperm - t = t[:,invperm(permx[permz])]::T - return MixedDestabilizer(t, r+s)::MixedDestabilizer{T} + t = t[:,invperm(permx[permz])] + return MixedDestabilizer(t, r+s) end if reportperm return (MixedDestabilizer(t, r+s)::MixedDestabilizer{T}, r, permx, permz) @@ -580,7 +582,7 @@ function MixedDestabilizer(stab::Stabilizer{T}; undoperm=true, reportperm=false) end function MixedDestabilizer(d::Destabilizer, r::Int) - l,n = size(d.tab) + l,n = size(tab(d)) if l==2n MixedDestabilizer(tab(d), r) else @@ -588,7 +590,7 @@ function MixedDestabilizer(d::Destabilizer, r::Int) end end function MixedDestabilizer(d::Destabilizer) - l,n = size(d.tab) + l,n = size(tab(d)) if l==2n MixedDestabilizer(d, nqubits(d)) else @@ -599,9 +601,9 @@ end MixedDestabilizer(d::MixedStabilizer) = MixedDestabilizer(stabilizerview(d)) MixedDestabilizer(d::MixedDestabilizer) = d -Base.length(d::MixedDestabilizer) = length(d.tab)÷2 +Base.length(d::MixedDestabilizer) = length(tab(d))÷2 -Base.copy(d::MixedDestabilizer) = MixedDestabilizer(copy(d.tab),d.rank) +Base.copy(d::MixedDestabilizer) = MixedDestabilizer(copy(tab(d)),rank(d)) ############################## # Subtableau views @@ -610,17 +612,17 @@ Base.copy(d::MixedDestabilizer) = MixedDestabilizer(copy(d.tab),d.rank) """A view of the subtableau corresponding to the stabilizer. See also [`tab`](@ref), [`destabilizerview`](@ref), [`logicalxview`](@ref), [`logicalzview`](@ref)""" @inline stabilizerview(s::Stabilizer) = s @inline stabilizerview(s::Destabilizer) = Stabilizer(@view tab(s)[end÷2+1:end]) -@inline stabilizerview(s::MixedStabilizer) = Stabilizer(@view tab(s)[1:s.rank]) -@inline stabilizerview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[end÷2+1:end÷2+s.rank]) +@inline stabilizerview(s::MixedStabilizer) = Stabilizer(@view tab(s)[1:rank(s)]) +@inline stabilizerview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[end÷2+1:end÷2+rank(s)]) """A view of the subtableau corresponding to the destabilizer. See also [`tab`](@ref), [`stabilizerview`](@ref), [`logicalxview`](@ref), [`logicalzview`](@ref)""" @inline destabilizerview(s::Destabilizer) = Stabilizer(@view tab(s)[1:end÷2]) -@inline destabilizerview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[1:s.rank]) +@inline destabilizerview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[1:rank(s)]) """A view of the subtableau corresponding to the logical X operators. See also [`tab`](@ref), [`stabilizerview`](@ref), [`destabilizerview`](@ref), [`logicalzview`](@ref)""" -@inline logicalxview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[s.rank+1:end÷2]) +@inline logicalxview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[rank(s)+1:end÷2]) """A view of the subtableau corresponding to the logical Z operators. See also [`tab`](@ref), [`stabilizerview`](@ref), [`destabilizerview`](@ref), [`logicalxview`](@ref)""" -@inline logicalzview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[end÷2+s.rank+1:end]) +@inline logicalzview(s::MixedDestabilizer) = Stabilizer(@view tab(s)[end÷2+rank(s)+1:end]) """The number of qubits of a given state.""" @inline nqubits(s::AbstractStabilizer) = nqubits(tab(s)) @@ -910,6 +912,30 @@ function unsafe_bitfindnext_(chunks::AbstractVector{T}, start::Int) where T<:Uns return nothing end +""" +$(TYPEDSIGNATURES) + +Computes bitmask indices for an unsigned integer at index `i` +within the binary structure of a `Tableau` or `PauliOperator`. + +For `Tableau`, the method operates on the `.xzs` field, while +for `PauliOperator`, it uses the `.xz` field. It calculates +the following values based on the index `i`: + +- `lowbit`, the lowest bit. +- `ibig`, the index of the word containing the bit. +- `ismall`, the position of the bit within the word. +- `ismallm`, a bitmask isolating the specified bit. +""" +@inline function get_bitmask_idxs(xzs::AbstractArray{<:Unsigned}, i::Int) + T = eltype(xzs) + lowbit = T(1) + ibig = _div(T, i-1) + 1 + ismall = _mod(T, i-1) + ismallm = lowbit << ismall + return lowbit, ibig, ismall, ismallm +end + """Permute the qubits (i.e., columns) of the tableau in place.""" function Base.permute!(s::Tableau, perm::AbstractVector) for r in 1:size(s,1) diff --git a/src/affectedqubits.jl b/src/affectedqubits.jl index da6314a65..3609710e3 100644 --- a/src/affectedqubits.jl +++ b/src/affectedqubits.jl @@ -15,6 +15,5 @@ affectedqubits(c::CliffordOperator) = 1:nqubits(c) affectedqubits(c::ClassicalXOR) = () affectedbits(o) = () -affectedbits(m::sMRZ) = (m.bit,) -affectedbits(m::sMZ) = (m.bit,) +affectedbits(m::Union{sMRZ,sMZ,sMRX,sMX,sMRY,sMY}) = m.bit==0 ? () : (m.bit,) affectedbits(c::ClassicalXOR) = (c.bits..., c.store) diff --git a/src/dense_cliffords.jl b/src/dense_cliffords.jl index baec10f6e..e4113e0e2 100644 --- a/src/dense_cliffords.jl +++ b/src/dense_cliffords.jl @@ -62,7 +62,7 @@ CliffordOperator(op::CliffordOperator) = op CliffordOperator(paulis::AbstractVector{<:PauliOperator}) = CliffordOperator(Tableau(paulis)) CliffordOperator(destab::Destabilizer) = CliffordOperator(tab(destab)) -Base.:(==)(l::CliffordOperator, r::CliffordOperator) = l.tab == r.tab +Base.:(==)(l::CliffordOperator, r::CliffordOperator) = tab(l) == tab(r) Base.hash(c::T, h::UInt) where {T<:CliffordOperator} = hash(T, hash(tab(c), h)) Base.getindex(c::CliffordOperator, args...) = getindex(tab(c), args...) @@ -85,16 +85,16 @@ digits_subchars = collect("₀₁₂₃₄₅₆₇₈₉") digits_substr(n::Int,nwidth::Int) = join(([digits_subchars[d+1] for d in reverse(digits(n, pad=nwidth))])) function Base.copy(c::CliffordOperator) - CliffordOperator(copy(c.tab)) + CliffordOperator(copy(tab(c))) end -@inline nqubits(c::CliffordOperator) = nqubits(c.tab) +@inline nqubits(c::CliffordOperator) = nqubits(tab(c)) -Base.zero(c::CliffordOperator) = CliffordOperator(zero(c.tab)) +Base.zero(c::CliffordOperator) = CliffordOperator(zero(tab(c))) Base.zero(::Type{<:CliffordOperator}, n) = CliffordOperator(zero(Tableau, 2n, n)) function Base.:(*)(l::AbstractCliffordOperator, r::CliffordOperator) - tab = copy(r.tab) + tab = copy(QuantumClifford.tab(r)) apply!(Stabilizer(tab),l) # TODO maybe not the most elegant way to perform apply!(::Tableau, gate) CliffordOperator(tab) end @@ -106,7 +106,7 @@ end # TODO create Base.permute! and getindex(..., permutation_array) function permute(c::CliffordOperator,p) # TODO this is a slow stupid implementation - CliffordOperator(Tableau([c.tab[i][p] for i in 1:2*nqubits(c)][vcat(p,p.+nqubits(c))])) + CliffordOperator(Tableau([tab(c)[i][p] for i in 1:2*nqubits(c)][vcat(p,p.+nqubits(c))])) end """Nonvectorized version of `apply!` used for unit tests.""" diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl index cdda7742e..ff4444694 100644 --- a/src/ecc/ECC.jl +++ b/src/ecc/ECC.jl @@ -1,9 +1,16 @@ module ECC -using LinearAlgebra -using LinearAlgebra: I -using QuantumClifford -using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer +using LinearAlgebra: LinearAlgebra, I, rank, tr +using QuantumClifford: QuantumClifford, AbstractOperation, AbstractStabilizer, + AbstractTwoQubitOperator, Stabilizer, PauliOperator, + random_brickwork_clifford_circuit, random_all_to_all_clifford_circuit, + canonicalize!, canonicalize_gott!, + logicalxview, logicalzview, stabilizerview, destabilizerview, tab, phases, + sCNOT, sSWAP, sHadamard, sPhase, sInvPhase, + sZCX, sZCY, sZCZ, sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZ, sX, sY, sMRZ, sMRX, + single_x, single_y, single_z, random_pauli!, PauliError, + apply!, comm, comm!, stab_to_gf2, embed, @S_str, affectedqubits, affectedbits, + pftrajectories, pfmeasurements, mctrajectories import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits using DocStringExtensions using Combinatorics: combinations diff --git a/src/ecc/circuits.jl b/src/ecc/circuits.jl index 45aabf1a7..72b3f338c 100644 --- a/src/ecc/circuits.jl +++ b/src/ecc/circuits.jl @@ -74,6 +74,7 @@ function perm_to_transpositions(perm) for i in n:-1:1 if perm[i]!=i j = findfirst(==(i), perm) + @assert !isnothing(j) push!(transpositions, (i, j)) perm[j] = perm[i] end diff --git a/src/ecc/codes/util.jl b/src/ecc/codes/util.jl index 1063785f9..b96a75120 100644 --- a/src/ecc/codes/util.jl +++ b/src/ecc/codes/util.jl @@ -6,3 +6,6 @@ function hgp(h₁,h₂) hz = hcat(kron(LinearAlgebra.I(n₁), h₂), kron(h₁', LinearAlgebra.I(r₂))) hx, hz end + +"""Implemented in a package extension with Hecke.""" +function check_repr_commutation_relation end diff --git a/src/ecc/decoder_pipeline.jl b/src/ecc/decoder_pipeline.jl index 3ac2680cf..aeb481851 100644 --- a/src/ecc/decoder_pipeline.jl +++ b/src/ecc/decoder_pipeline.jl @@ -251,7 +251,7 @@ function create_lookup_table(code::Stabilizer) for bit_to_be_flipped in 1:qubits for error_type in [single_x, single_y, single_z] # Generate e⃗ - error = error_type(qubits, bit_to_be_flipped) + error = error_type(qubits, bit_to_be_flipped)::PauliOperator{Array{UInt8, 0}, Vector{UInt}} # Calculate s⃗ # (check which stabilizer rows do not commute with the Pauli error) syndrome = comm(error, code) diff --git a/src/fastmemlayout.jl b/src/fastmemlayout.jl index f7019cd1a..67e2bbbd7 100644 --- a/src/fastmemlayout.jl +++ b/src/fastmemlayout.jl @@ -20,14 +20,14 @@ fastrow(t::Tableau{Tₚᵥ,Tₘ}) where {Tₚᵥ, Tₘ<:Adjoint} = Tableau(t.pha fastcolumn(t::Tableau{Tₚᵥ,Tₘ}) where {Tₚᵥ, Tₘ} = Tableau(t.phases, t.nqubits, collect(t.xzs')') fastcolumn(t::Tableau{Tₚᵥ,Tₘ}) where {Tₚᵥ, Tₘ<:Adjoint} = t -fastrow(s::Stabilizer) = Stabilizer(fastrow(s.tab)) -fastcolumn(s::Stabilizer) = Stabilizer(fastcolumn(s.tab)) +fastrow(s::Stabilizer) = Stabilizer(fastrow(tab(s))) +fastcolumn(s::Stabilizer) = Stabilizer(fastcolumn(tab(s))) -fastrow(s::Destabilizer) = Destabilizer(fastrow(s.tab)) -fastcolumn(s::Destabilizer) = Destabilizer(fastcolumn(s.tab)) +fastrow(s::Destabilizer) = Destabilizer(fastrow(tab(s))) +fastcolumn(s::Destabilizer) = Destabilizer(fastcolumn(tab(s))) -fastrow(s::MixedStabilizer) = MixedStabilizer(fastrow(s.tab), s.rank) -fastcolumn(s::MixedStabilizer) = MixedStabilizer(fastcolumn(s.tab), s.rank) +fastrow(s::MixedStabilizer) = MixedStabilizer(fastrow(tab(s)), rank(s)) +fastcolumn(s::MixedStabilizer) = MixedStabilizer(fastcolumn(tab(s)), rank(s)) -fastrow(s::MixedDestabilizer) = MixedDestabilizer(fastrow(s.tab), s.rank) -fastcolumn(s::MixedDestabilizer) = MixedDestabilizer(fastcolumn(s.tab), s.rank) +fastrow(s::MixedDestabilizer) = MixedDestabilizer(fastrow(tab(s)), rank(s)) +fastcolumn(s::MixedDestabilizer) = MixedDestabilizer(fastcolumn(tab(s)), rank(s)) diff --git a/src/grouptableaux.jl b/src/grouptableaux.jl index b0f76e1c2..6f76802cf 100644 --- a/src/grouptableaux.jl +++ b/src/grouptableaux.jl @@ -1,5 +1,109 @@ -using Graphs -using LinearAlgebra +""" +A tableau representation of the non-commutative canonical form of a set of Paulis, +which is used in [`commutify`](@ref). + +They are organized in the same form as [`MixedDestabilizer`](@ref) +with a stabilizer, destabilizer, logical X, and logical Z components. +""" +mutable struct SubsystemCodeTableau{T <: Tableau} <: AbstractStabilizer + tab::T + index::Int + r::Int + m::Int + k::Int +end + +function SubsystemCodeTableau(t::Tableau) + index = 1 + for i in range(1, stop=length(t), step=2) + if i + 1 > length(t) + break + end + if comm(t[i], t[i+1]) == 0x01 + index = i+2 # index to split t into non-commuting pairs and commuting operators + end + end + s = Stabilizer(t[index:length(t)]) + ind = 1 + if length(s)>nqubits(s)#if stabilizer is overdetermined, Destabilizer constructor throws error + m = length(s) + tab = zero(Tableau, length(t)+length(s), nqubits(t)) + for i in s + tab[ind] = zero(PauliOperator, nqubits(s)) + ind+=1 + end + else + d = Destabilizer(s) + m = length(d) + tab = zero(Tableau, length(t)+length(d), nqubits(t)) + for p in destabilizerview(d) + tab[ind] = p + ind+=1 + end + end + for i in range(1, stop=index-1, step=2) + tab[ind] = t[i] + ind+=1 + end + for p in s + tab[ind] = p + ind+=1 + end + for i in range(2, stop=index, step=2) + tab[ind] = t[i] + ind+=1 + end + return SubsystemCodeTableau(tab, index, length(s), m, Int((index-1)/2)) +end + + +Base.copy(t::SubsystemCodeTableau) = SubsystemCodeTableau(copy(t.tab), t.index, t.r, t.m, t.k) + +function Base.show(io::IO, t::SubsystemCodeTableau) + r = t.r+t.m+2*t.k + q = nqubits(t) + if get(io, :compact, false) | haskey(io, :typeinfo) + print(io, "MixedDestablizer $r×$q") + elseif get(io, :limit, false) + h,w = displaysize(io) + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(t.tab,2)-4),0)) + _show(io, destabilizerview(t).tab, w, h÷4) + if r != q + println(io) + println(io, "𝒳" * "━"^max(min(w-5,size(t.tab,2)),0)) + _show(io, logicalxview(t).tab, w, h÷4) + end + println(io) + println(io, "𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(t.tab,2)-2),0)) + _show(io, stabilizerview(t).tab, w, h÷4) + if r != q + println(io) + println(io, "𝒵" * "━"^max(min(w-5,size(t.tab,2)),0)) + _show(io, logicalzview(t).tab, w, h÷4) + end + else + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-4,0)) + _show(io, destabilizerview(t).tab, missing, missing) + if r != q + println(io) + println(io, "𝒳ₗ" * "━"^max(size(t.tab,2),0)) + _show(io, logicalxview(t).tab, missing, missing) + end + println(io) + println(io, "𝒮𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-2,0)) + _show(io, stabilizerview(t).tab, missing, missing) + if r != q + println(io) + println(io, "𝒵ₗ" * "━"^max(size(t.tab,2))) + _show(io, logicalzview(t).tab, missing, missing) + end + end +end + +@inline stabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[s.m+s.k+1:s.m+s.k+s.r]) +@inline destabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[1:s.r]) +@inline logicalxview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[s.m+1:s.m+s.k]) +@inline logicalzview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[length(s.tab)-s.k+1:length(s.tab)]) """ Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)). @@ -14,20 +118,24 @@ julia> groupify(S"XZ ZX") + YY ``` """ -function groupify(s::Stabilizer) - # Create a `Tableau` of 2ⁿ n-qubit identity Pauli operators(where n is the size of - # `Stabilizer` s), then multiply each one by a different subset of the elements in s to +function groupify(s::Stabilizer; phases=false) + # Create a `Tableau` of 2ⁿ n-qubit identity Pauli operators(where n is the size of + # `Stabilizer` s), then multiply each one by a different subset of the elements in s to # create all 2ⁿ unique elements in the group generated by s, then return the `Tableau`. - n = length(s)::Int - group = zero(Tableau, 2^n, nqubits(s)) + if phases == true + throw(ArgumentError("in groupify phases=true functionality not yet implemented")) + end + gen_set = minimal_generating_set(s) + n = length(gen_set)::Int + group = zero(Tableau, 2^n, nqubits(gen_set)) for i in 0:2^n-1 for (digit_order, j) in enumerate(digits(i, base=2, pad=n)) if j == 1 - group[i+1] *= s[digit_order] + mul_left!(group, i+1, tab(gen_set), digit_order, phases=Val(false)) end end end - return group + return group end @@ -44,14 +152,14 @@ julia> minimal_generating_set(S"__ XZ ZX YY") ``` """ function minimal_generating_set(s::Stabilizer) - # Canonicalize `Stabilizer` s, then return a `Stabilizer` with all non-identity Pauli operators - # in the result. If s consists of only identity operators, return the negative + # Canonicalize `Stabilizer` s, then return a `Stabilizer` with all non-identity Pauli operators + # in the result. If s consists of only identity operators, return the negative # identity operator if one is contained in s, and the positive identity operator otherwise. s, _, r = canonicalize!(copy(s), ranks=true) if r == 0 gs = zero(Stabilizer, 1, nqubits(s)) if 0x02 in phases(s) - gs[1] = -1 * gs[1] + phases(gs)[1] = 0x2 end return gs else @@ -60,7 +168,154 @@ function minimal_generating_set(s::Stabilizer) end """ -Return the full Pauli group of a given length. Phases are ignored by default, +For a not-necessarily commutative set of Paulis, return a generating set of the form +⟨A₁, A₂, ... Aₖ, Aₖ₊₁, ... Aₘ, B₁, B₂, ... Bₖ⟩ where pairs Aₖ, Bₖ anticommute and all other pairings commute. Based on [RevModPhys.87.307](@cite) + +Returns the generating set as a data structure of type [`SubsystemCodeTableau`](@ref). The +`logicalxview` function returns the ⟨A₁, A₂,... Aₖ⟩, and the `logicalzview` +function returns ⟨B₁, B₂, ... Bₖ⟩. `stabilizerview` returns ⟨Aₖ₊₁, ... Aₘ⟩ +as a Stabilizer, and `destabilizerview` returns the Destabilizer of that Stabilizer. + +Phases are zeroed-out in this canonicalization. + +```jldoctest +julia> canonicalize_noncomm(T"XX XZ XY") +𝒟ℯ𝓈𝓉𝒶𝒷 ++ Z_ +𝒳━━ ++ XX +𝒮𝓉𝒶𝒷 ++ X_ +𝒵━━ ++ XZ +``` +""" +function canonicalize_noncomm(t::Tableau) + for i in eachindex(t) phases(t)[i] = 0x00 end + loc = zero(Tableau, 2*length(t), nqubits(t)) + x_index = 0 + z_index = length(loc) + for i in eachindex(t) + for j in eachindex(t) + if comm(t[i], t[j]) == 0x01 + for k in eachindex(t) + if k !=i && k != j + if comm(t[k], t[i]) == 0x01 + mul_left!(t, k, j, phases=Val(false)) + end + if comm(t[k], t[j]) == 0x01 + mul_left!(t, k, i, phases=Val(false)) + end + end + end + if !(t[i] in loc) + x_index+=1 + loc[x_index]= t[i] + end + if !(t[j] in loc) + loc[z_index]= t[j] + z_index-=1 + end + end + end + end + ind = 2*x_index+1 #format tableau for SubsystemCodeTableau + k = x_index + m = length(t)-2*k + for i in range(k, stop =1, step=-1) + loc[i+m]= loc[i] + end + i = m+1 + j = m+k + while i < j #ensure anticommutative pairs are lined up correctly + loc[i], loc[j], = loc[j], loc[i] + i+=1 + j-=1 + end + s_index = m+k+1 + for i in eachindex(t) + if !(t[i] in loc) + loc[s_index]= t[i] + s_index+=1 + end + end + r = s_index-m-k-1 + if r <= nqubits(t) + d = destabilizerview(Destabilizer(Stabilizer(loc[m+k+1:s_index-1]))) + for i in eachindex(d) + loc[i] =d[i] + end + else for i in 1:m zero!(loc, i) end end + return SubsystemCodeTableau(loc, ind, r, m, k) +end + +""" +For a not-necessarily commutative set of Paulis S, +computed S', the [non-commutative canonical form](@ref canonicalize_noncomm) of of S. +For each pair Aₖ, Bₖ of anticommutative Paulis in S', add a qubit to each Pauli in the set: +X to Aₖ, Z to Bₖ, and I to each other operator to produce S'', a fully commutative set. Return +S'' as well as a list of the indices of the added qubits. + +The returned object is a Stabilizer that is also useful for the [`matroid_parent`](@ref) function. + +```jldoctest +julia> commutify(T"XX XZ XY")[1] ++ XXX ++ X__ ++ XZZ + +julia> commutify(T"XX XZ XY")[2] +3:3 +``` +""" +function commutify(t) + loc = canonicalize_noncomm(t) + commutative = zero(Stabilizer, 2*loc.k+loc.r, nqubits(loc)+loc.k) + puttableau!(commutative, logicalxview(loc), 0,0) + puttableau!(commutative, stabilizerview(loc), loc.k,0) + puttableau!(commutative, logicalzview(loc), loc.k+loc.r,0) + x= T"X" + z=T"Z" + for i in 0:(loc.k-1) + puttableau!(commutative, x, i, nqubits(loc)+i) + puttableau!(commutative, z, i+loc.r+loc.k, nqubits(loc)+i) + end + to_delete = nqubits(loc)+1:nqubits(loc)+loc.k + return commutative, to_delete +end + +""" +For a given set S of Paulis that does not necessarily represent a state, +return a set of Paulis S' that represents a state. +S' is a superset of [commutified](@ref commutify) S. +Additionally returns two arrays representing deletions needed to produce S. +Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite) + +By deleting the qubits in the first output array from S', taking the [`normalizer`](@ref) of S', then +deleting the qubits in the second returned array from the [`normalizer`](@ref) of S', S is reproduced. + +```jldoctest +julia> matroid_parent(T"XX")[1] ++ X_X ++ XX_ ++ ZZZ + +julia> matroid_parent(T"XX")[2] +3:3 + +julia> matroid_parent(T"XX")[3] +3:2 +``` +""" +function matroid_parent(t::Tableau) + com, d1= commutify(t) + norm = normalizer(com.tab) + state, d2 = commutify(norm) + return state, d2, d1 +end + +""" +Return the full Pauli group of a given length. Phases are ignored by default, but can be included by setting `phases=true`. ```jldoctest @@ -130,15 +385,21 @@ julia> normalizer(T"X") + X ``` """ -function normalizer(t::Tableau) - # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each - # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector +function normalizer(t::Tableau; phases=false) + # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each + # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector # converted to `Tableau`. n = nqubits(t) - pgroup = pauligroup(n, phases=false) - ptype = typeof(t[1]) - normalizer = ptype[] - for p in pgroup + ptype =typeof(P"I") + norm = ptype[] + + p = zero(PauliOperator, n) + paulis = ((false, false), (true, false), (false, true), (true, true)) + for i in Iterators.product(Iterators.repeated(paulis, n)...) + zero!(p) + for (j, k) in enumerate(i) + p[j] = k + end commutes = true for q in t if comm(p, q) == 0x01 @@ -146,11 +407,15 @@ function normalizer(t::Tableau) end end if commutes - push!(normalizer, p) + push!(norm, copy(p)) + end + if phases + for phase in [-1, 1im, -1im] + push!(norm, phase *p) + end end end - - return Tableau(normalizer) + return Tableau(norm) end """ @@ -161,7 +426,7 @@ julia> centralizer(T"XX ZZ _Z") + ZZ ``` """ -function centralizer(t::Tableau) +function centralizer(t::Tableau) center = typeof(t[1])[] for P in t commutes = 0 @@ -175,7 +440,7 @@ function centralizer(t::Tableau) push!(center, P) end end - if length(center) == 0 + if length(center) == 0 return Tableau(zeros(Bool, 1,1)) end c = Tableau(center) @@ -183,8 +448,8 @@ function centralizer(t::Tableau) end """ -Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to -the given subset, without the entries corresponding to subset. +Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to +the given subset, without the entries corresponding to subset. Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite) ```jldoctest julia> contractor(S"_X X_", [1]) @@ -196,9 +461,9 @@ function contractor(s::Stabilizer, subset) for p in s contractable = true for i in subset - if p[i] != (false, false) - contractable = false - break + if p[i] != (false, false) + contractable = false + break end end if contractable push!(result, p[setdiff(1:length(p), subset)]) end @@ -208,4 +473,4 @@ function contractor(s::Stabilizer, subset) else return Tableau(zeros(Bool, 1,1)) end -end \ No newline at end of file +end diff --git a/src/linalg.jl b/src/linalg.jl index b95fba888..ad252e95f 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -260,7 +260,7 @@ function tensor(ops::CliffordOperator...) # TODO implement \otimes for Destabili last_zrow = ntot last_xrow = 0 for op in ops - t = op.tab + t = QuantumClifford.tab(op) _, last_zrow, _ = puttableau!(tab, (@view t[end÷2+1:end]), last_zrow, last_xrow) _, last_xrow, _ = puttableau!(tab, (@view t[1:end÷2]), last_xrow, last_xrow) end diff --git a/src/mul_leftright.jl b/src/mul_leftright.jl index 122fb1e61..08f905664 100644 --- a/src/mul_leftright.jl +++ b/src/mul_leftright.jl @@ -56,12 +56,6 @@ function mul_ordered_lv!(r::AbstractVector{T}, l::AbstractVector{T}; phases::Val end =# -function mul_ordered!(r::SubArray{T,1,P,I1,L1}, l::SubArray{T,1,P,I2,L2}; phases::Val{B}=Val(true)) where {T<:Unsigned, B, I1, I2, L1, L2, P<:Adjoint} - # This method exists because SIMD.jl does not play well with Adjoint - # Delete it and try `QuantumClifford.mul_left!(fastcolumn(random_stabilizer(194)), 2, 1)` # works fine for 192 - _mul_ordered_nonvec!(r,l; phases=B) -end - function mul_ordered!(r::SubArray{T,1,P,I2,L2}, l::AbstractVector{T}; phases::Val{B}=Val(true)) where {T<:Unsigned, B, I2, L2, P<:Adjoint} # This method exists because SIMD.jl does not play well with Adjoint _mul_ordered_nonvec!(r,l; phases=B) @@ -163,6 +157,12 @@ end # On Tableaux ############################## +@inline function mul_left!(s::Tableau, m, t::Tableau, i; phases::Val{B}=Val(true)) where B + extra_phase = mul_left!((@view s.xzs[:,m]), (@view t.xzs[:,i]); phases=phases) + B && (s.phases[m] = (extra_phase+s.phases[m]+s.phases[i])&0x3) + s +end + @inline function mul_left!(s::Tableau, m, i; phases::Val{B}=Val(true)) where B extra_phase = mul_left!((@view s.xzs[:,m]), (@view s.xzs[:,i]); phases=phases) B && (s.phases[m] = (extra_phase+s.phases[m]+s.phases[i])&0x3) diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl index 45e9763c8..af84bd63b 100644 --- a/src/pauli_frames.jl +++ b/src/pauli_frames.jl @@ -82,12 +82,8 @@ end function apply!(frame::PauliFrame, op::sMZ) # TODO sMY, and faster sMX op.bit == 0 && return frame i = op.qubit - xzs = frame.frame.tab.xzs - T = eltype(xzs) - lowbit = T(1) - ibig = _div(T,i-1)+1 - ismall = _mod(T,i-1) - ismallm = lowbit<<(ismall) + xzs = tab(frame.frame).xzs + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) @inbounds @simd for f in eachindex(frame) should_flip = !iszero(xzs[ibig,f] & ismallm) @@ -99,12 +95,8 @@ end function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX i = op.qubit - xzs = frame.frame.tab.xzs - T = eltype(xzs) - lowbit = T(1) - ibig = _div(T,i-1)+1 - ismall = _mod(T,i-1) - ismallm = lowbit<<(ismall) + xzs = tab(frame.frame).xzs + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) if op.bit != 0 @inbounds @simd for f in eachindex(frame) @@ -122,45 +114,39 @@ end function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int) p = noise.p - T = eltype(frame.frame.tab.xzs) + xzs = tab(frame.frame).xzs - lowbit = T(1) - ibig = _div(T,i-1)+1 - ismall = _mod(T,i-1) - ismallm = lowbit<<(ismall) + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) p = p/3 @inbounds @simd for f in eachindex(frame) r = rand() if r < p # X error - frame.frame.tab.xzs[ibig,f] ⊻= ismallm + xzs[ibig,f] ⊻= ismallm elseif r < 2p # Z error - frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm + xzs[end÷2+ibig,f] ⊻= ismallm elseif r < 3p # Y error - frame.frame.tab.xzs[ibig,f] ⊻= ismallm - frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm + xzs[ibig,f] ⊻= ismallm + xzs[end÷2+ibig,f] ⊻= ismallm end end return frame end function applynoise!(frame::PauliFrame,noise::PauliNoise,i::Int) - T = eltype(frame.frame.tab.xzs) + xzs = tab(frame.frame).xzs - lowbit = T(1) - ibig = _div(T,i-1)+1 - ismall = _mod(T,i-1) - ismallm = lowbit<<(ismall) + lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i) @inbounds @simd for f in eachindex(frame) r = rand() if r < noise.px # X error - frame.frame.tab.xzs[ibig,f] ⊻= ismallm + xzs[ibig,f] ⊻= ismallm elseif r < noise.px+noise.pz # Z error - frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm + xzs[end÷2+ibig,f] ⊻= ismallm elseif r < noise.px+noise.pz+noise.py # Y error - frame.frame.tab.xzs[ibig,f] ⊻= ismallm - frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm + xzs[ibig,f] ⊻= ismallm + xzs[end÷2+ibig,f] ⊻= ismallm end end return frame diff --git a/src/pauli_operator.jl b/src/pauli_operator.jl index e556628bc..2246d4d5f 100644 --- a/src/pauli_operator.jl +++ b/src/pauli_operator.jl @@ -77,7 +77,7 @@ function zbit(p::PauliOperator) [(word>>s)&one==one for word in zview(p) for s in 0:size-1][begin:p.nqubits] end -function _P_str(a) +function _P_str(a::Union{String,SubString{String}}) letters = filter(x->occursin(x,"_IZXY"),a) phase = phasedict[strip(filter(x->!occursin(x,"_IZXY"),a))] PauliOperator(phase, [l=='X'||l=='Y' for l in letters], [l=='Z'||l=='Y' for l in letters]) @@ -87,19 +87,23 @@ macro P_str(a) quote _P_str($a) end end -Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, i::Int) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}} = (p.xz[_div(Tᵥₑ, i-1)+1] & Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))!=0x0, (p.xz[end÷2+_div(Tᵥₑ,i-1)+1] & Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))!=0x0 +function Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, i::Int) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}} + _, ibig, _, ismallm = get_bitmask_idxs(p.xz,i) + ((p.xz[ibig] & ismallm) != 0x0)::Bool, ((p.xz[end÷2+ibig] & ismallm) != 0x0)::Bool +end Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, r) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}} = PauliOperator(p.phase[], xbit(p)[r], zbit(p)[r]) function Base.setindex!(p::PauliOperator{Tₚ,Tᵥ}, (x,z)::Tuple{Bool,Bool}, i) where {Tₚ, Tᵥₑ, Tᵥ<:AbstractVector{Tᵥₑ}} + _, ibig, _, ismallm = get_bitmask_idxs(p.xz,i) if x - p.xz[_div(Tᵥₑ,i-1)+1] |= Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1) + p.xz[ibig] |= ismallm else - p.xz[_div(Tᵥₑ,i-1)+1] &= ~(Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1)) + p.xz[ibig] &= ~(ismallm) end if z - p.xz[end÷2+_div(Tᵥₑ,i-1)+1] |= Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1) + p.xz[end÷2+ibig] |= ismallm else - p.xz[end÷2+_div(Tᵥₑ,i-1)+1] &= ~(Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1)) + p.xz[end÷2+ibig] &= ~(ismallm) end p end @@ -176,7 +180,8 @@ end function embed(n::Int, indices, p::PauliOperator) if nqubits(p) == length(indices) pout = zero(typeof(p), n) - @inbounds @simd for i in eachindex(indices) + @inbounds @simd for i_ in eachindex(indices) + i = i_::Int pout[indices[i]] = p[i] end pout.phase[] = p.phase[] diff --git a/src/project_trace_reset.jl b/src/project_trace_reset.jl index 19e49521b..4f7fa6211 100644 --- a/src/project_trace_reset.jl +++ b/src/project_trace_reset.jl @@ -42,17 +42,15 @@ function _generate!(pauli::PauliOperator{Tₚ,Tᵥ}, stabilizer::Stabilizer{Tabl xzs = tab(stabilizer).xzs xs = @view xzs[1:end÷2,:] zs = @view xzs[end÷2+1:end,:] - lowbit = Tₘₑ(0x1) zerobit = Tₘₑ(0x0) px,pz = xview(pauli), zview(pauli) used_indices = Int[] used = 0 # remove Xs while (i=unsafe_bitfindnext_(px,1); i !== nothing) # TODO awkward notation due to https://github.com/JuliaLang/julia/issues/45499 - jbig = _div(Tₘₑ,i-1)+1 - jsmall = lowbit<<_mod(Tₘₑ,i-1) - candidate = findfirst(e->e&jsmall!=zerobit, # TODO some form of reinterpret might be faster than equality check - xs[jbig,used+1:end]) + _, ibig, _, ismallm = get_bitmask_idxs(xzs,i) + candidate = findfirst(e->e&ismallm!=zerobit, # TODO some form of reinterpret might be faster than equality check + xs[ibig,used+1:end]) if isnothing(candidate) return nothing else @@ -63,10 +61,9 @@ function _generate!(pauli::PauliOperator{Tₚ,Tᵥ}, stabilizer::Stabilizer{Tabl end # remove Zs while (i=unsafe_bitfindnext_(pz,1); i !== nothing) # TODO awkward notation due to https://github.com/JuliaLang/julia/issues/45499 - jbig = _div(Tₘₑ,i-1)+1 - jsmall = lowbit<<_mod(Tₘₑ,i-1) - candidate = findfirst(e->e&jsmall!=zerobit, # TODO some form of reinterpret might be faster than equality check - zs[jbig,used+1:end]) + _, ibig, _, ismallm = get_bitmask_idxs(xzs,i) + candidate = findfirst(e->e&ismallm!=zerobit, # TODO some form of reinterpret might be faster than equality check + zs[ibig,used+1:end]) if isnothing(candidate) return nothing else @@ -339,7 +336,7 @@ end 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 + tab = QuantumClifford.tab(d) stabilizer = stabilizerview(d) destabilizer = destabilizerview(d) r = trusted_rank(d) @@ -377,7 +374,7 @@ 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 + tab = QuantumClifford.tab(d) stabilizer = stabilizerview(d) destabilizer = destabilizerview(d) r = trusted_rank(d) @@ -497,7 +494,7 @@ end """Internal method used to implement [`projectX!`](@ref), [`projectZ!`](@ref), and [`projectY!`](@ref).""" function project_cond!(d::MixedDestabilizer,qubit::Int,cond::Val{IS},reset::Val{RESET};keep_result::Bool=true,phases::Val{PHASES}=Val(true)) where {IS,RESET,PHASES} anticommutes = 0 - tab = d.tab + tab = QuantumClifford.tab(d) stabilizer = stabilizerview(d) destabilizer = destabilizerview(d) r = d.rank @@ -647,7 +644,7 @@ function traceout!(s::Union{MixedStabilizer, MixedDestabilizer}, qubits; phases= if rank return (s, i) else return s end end -function _expand_pauli(pauli,qubits,n) # TODO rename and make public +function _expand_pauli(pauli::PauliOperator,qubits,n) # TODO rename and make public expanded = zero(PauliOperator,n) for (ii, i) in enumerate(qubits) expanded[i] = pauli[ii] @@ -885,5 +882,6 @@ julia> delete_columns(S"XYZ YZX ZXY", [1,3]) See also: [`traceout!`](@ref) """ function delete_columns(𝒮::Stabilizer, subset) + if length(𝒮) == 0 return 𝒮 end return 𝒮[:, setdiff(1:nqubits(𝒮), subset)] -end \ No newline at end of file +end diff --git a/src/randoms.jl b/src/randoms.jl index 1d4fff388..c87932fd4 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -184,7 +184,7 @@ function nemo_inv(a, n)::Matrix{UInt8} end """Sample (h, S) from the distribution P_n(h, S) from Bravyi and Maslov Algorithm 1.""" -function quantum_mallows(rng, n) # each one is benchmakred in benchmarks/quantum_mallows.jl +function quantum_mallows(rng::AbstractRNG, n::Int) # each one is benchmarked in benchmarks/quantum_mallows.jl arr = collect(1:n) hadamard = falses(n) perm = zeros(Int64, n) @@ -202,7 +202,7 @@ end """ This function samples a number from 1 to `n` where `n >= 1` probability of outputting `i` is proportional to `2^i`""" -function sample_geometric_2(rng, n::Integer) +function sample_geometric_2(rng::AbstractRNG, n::Integer) n < 1 && throw(DomainError(n)) if n<30 k = rand(rng, 2:UInt(2)^n) @@ -217,7 +217,7 @@ function sample_geometric_2(rng, n::Integer) end """Assign (symmetric) random ints to off diagonals of matrix.""" -function fill_tril(rng, matrix, n; symmetric::Bool=false) +function fill_tril(rng::AbstractRNG, matrix, n; symmetric::Bool=false) # Add (symmetric) random ints to off diagonals @inbounds for row in 1:n, col in 1:row-1 b = rand(rng, Bool) diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl index e73622217..fcff4cb38 100644 --- a/src/symbolic_cliffords.jl +++ b/src/symbolic_cliffords.jl @@ -203,7 +203,7 @@ SingleQubitOperator(p::sInvSQRTY) = SingleQubitOperator(p.q, false, tr SingleQubitOperator(o::SingleQubitOperator) = o function SingleQubitOperator(op::CliffordOperator, qubit) nqubits(op)==1 || throw(DimensionMismatch("You are trying to convert a multiqubit `CliffordOperator` into a symbolic `SingleQubitOperator`.")) - SingleQubitOperator(qubit,op.tab[1,1]...,op.tab[2,1]...,(~).(iszero.(op.tab.phases))...) + SingleQubitOperator(qubit,tab(op)[1,1]...,tab(op)[2,1]...,(~).(iszero.(tab(op).phases))...) end SingleQubitOperator(op::CliffordOperator) = SingleQubitOperator(op, 1) @@ -308,6 +308,16 @@ macro qubitop2(name, kernel) end # x1 z1 x2 z2 @qubitop2 SWAP (x2 , z2 , x1 , z1 , false) + +@qubitop2 SWAPCX (x2 , z2⊻z1 , x2⊻x1 , z1 , ~iszero((x1 & z1 & x2 & z2) | (~x1 & z1 & x2 & ~z2))) +@qubitop2 InvSWAPCX (x2⊻x1 , z2 , x1 , z2⊻z1 , ~iszero((x1 & z1 & x2 & z2) | ( x1 &~z1 &~x2 & z2))) + +@qubitop2 ISWAP (x2 , x1⊻z2⊻x2 , x1 , x1⊻x2⊻z1 , ~iszero((x1 & z1 & ~x2) | (~x1 & x2 & z2))) +@qubitop2 InvISWAP (x2 , x1⊻z2⊻x2 , x1 , x1⊻x2⊻z1 , ~iszero((x1 &~z1 & ~x2) | (~x1 & x2 &~z2))) + +@qubitop2 CZSWAP (x2 , z2⊻x1 , x1 , x2⊻z1 , ~iszero((x1 & ~z1 & x2 & z2) | (x1 & z1 & x2 & ~z2))) +@qubitop2 CXSWAP (x2⊻x1 , z2 , x1 , z2⊻z1 , ~iszero((x1 & ~z1 &~x2 & z2) | (x1 & z1 & x2 & z2))) + @qubitop2 CNOT (x1 , z1⊻z2 , x2⊻x1 , z2 , ~iszero( (x1 & z1 & x2 & z2) | (x1 & z2 &~(z1|x2)) )) @qubitop2 CPHASE (x1 , z1⊻x2 , x2 , z2⊻x1 , ~iszero( (x1 & z1 & x2 &~z2) | (x1 &~z1 & x2 & z2) )) @@ -326,6 +336,15 @@ end @qubitop2 ZCrY (x1, x1⊻z1⊻x2⊻z2, x1⊻x2, x1⊻z2, ~iszero((x1 &~z1 & x2) | (x1 & ~z1 & ~z2) | (x1 & x2 & ~z2))) @qubitop2 InvZCrY (x1, x1⊻z1⊻x2⊻z2, x1⊻x2, x1⊻z2, ~iszero((x1 & z1 &~x2) | (x1 & z1 & z2) | (x1 &~x2 & z2))) +@qubitop2 SQRTZZ (x1 , x1⊻x2⊻z1 , x2 , x1⊻z2⊻x2 , ~iszero((x1 & z1 & ~x2) | (~x1 & x2 & z2))) +@qubitop2 InvSQRTZZ (x1 , x1⊻x2⊻z1 , x2 , x1⊻z2⊻x2 , ~iszero((x1 &~z1 & ~x2) | (~x1 & x2 &~z2))) + +@qubitop2 SQRTXX (z1⊻z2⊻x1, z1 , z1⊻x2⊻z2, z2 , ~iszero((~x1 & z1 &~z2) | (~z1 &~x2 & z2))) +@qubitop2 InvSQRTXX (z1⊻z2⊻x1, z1 , z1⊻x2⊻z2, z2 , ~iszero(( x1 & z1 &~z2) | (~z1 & x2 & z2))) + +@qubitop2 SQRTYY (z1⊻x2⊻z2, x1⊻z2⊻x2, x1⊻z1⊻z2, x1⊻x2⊻z1, ~iszero((~x1 &~z1 & x2 &~z2) | ( x1 &~z1 &~x2 &~z2) | ( x1 &~z1 & x2 & z2) | ( x1 & z1 & x2 &~z2))) +@qubitop2 InvSQRTYY (z1⊻x2⊻z2, x1⊻z2⊻x2, x1⊻z1⊻z2, x1⊻x2⊻z1, ~iszero(( x1 & z1 &~x2 & z2) | (~x1 & z1 & x2 & z2) | (~x1 & z1 &~x2 &~z2) | (~x1 &~z1 &~x2 & z2))) + #= To get the boolean formulas for the phase, it is easiest to first write down the truth table for the phase: for i in 0:15 @@ -370,20 +389,32 @@ function Base.show(io::IO, op::AbstractTwoQubitOperator) end end -LinearAlgebra.inv(op::sSWAP) = sSWAP(op.q1, op.q2) -LinearAlgebra.inv(op::sCNOT) = sCNOT(op.q1, op.q2) -LinearAlgebra.inv(op::sCPHASE) = sCPHASE(op.q1, op.q2) -LinearAlgebra.inv(op::sZCX) = sZCX(op.q1, op.q2) -LinearAlgebra.inv(op::sZCY) = sZCY(op.q1, op.q2) -LinearAlgebra.inv(op::sZCZ) = sZCZ(op.q1, op.q2) -LinearAlgebra.inv(op::sXCX) = sXCX(op.q1, op.q2) -LinearAlgebra.inv(op::sXCY) = sXCY(op.q1, op.q2) -LinearAlgebra.inv(op::sXCZ) = sXCZ(op.q1, op.q2) -LinearAlgebra.inv(op::sYCX) = sYCX(op.q1, op.q2) -LinearAlgebra.inv(op::sYCY) = sYCY(op.q1, op.q2) -LinearAlgebra.inv(op::sYCZ) = sYCZ(op.q1, op.q2) -LinearAlgebra.inv(op::sZCrY) = sInvZCrY(op.q1, op.q2) -LinearAlgebra.inv(op::sInvZCrY) = sZCrY(op.q1, op.q2) +LinearAlgebra.inv(op::sSWAP) = sSWAP(op.q1, op.q2) +LinearAlgebra.inv(op::sCNOT) = sCNOT(op.q1, op.q2) +LinearAlgebra.inv(op::sCPHASE) = sCPHASE(op.q1, op.q2) +LinearAlgebra.inv(op::sZCX) = sZCX(op.q1, op.q2) +LinearAlgebra.inv(op::sZCY) = sZCY(op.q1, op.q2) +LinearAlgebra.inv(op::sZCZ) = sZCZ(op.q1, op.q2) +LinearAlgebra.inv(op::sXCX) = sXCX(op.q1, op.q2) +LinearAlgebra.inv(op::sXCY) = sXCY(op.q1, op.q2) +LinearAlgebra.inv(op::sXCZ) = sXCZ(op.q1, op.q2) +LinearAlgebra.inv(op::sYCX) = sYCX(op.q1, op.q2) +LinearAlgebra.inv(op::sYCY) = sYCY(op.q1, op.q2) +LinearAlgebra.inv(op::sYCZ) = sYCZ(op.q1, op.q2) +LinearAlgebra.inv(op::sZCrY) = sInvZCrY(op.q1, op.q2) +LinearAlgebra.inv(op::sInvZCrY) = sZCrY(op.q1, op.q2) +LinearAlgebra.inv(op::sSWAPCX) = sInvSWAPCX(op.q1, op.q2) +LinearAlgebra.inv(op::sInvSWAPCX) = sSWAPCX(op.q1, op.q2) +LinearAlgebra.inv(op::sCZSWAP) = sCZSWAP(op.q1, op.q2) +LinearAlgebra.inv(op::sCXSWAP) = sSWAPCX(op.q1, op.q2) +LinearAlgebra.inv(op::sISWAP) = sInvISWAP(op.q1, op.q2) +LinearAlgebra.inv(op::sInvISWAP) = sISWAP(op.q1, op.q2) +LinearAlgebra.inv(op::sSQRTZZ) = sInvSQRTZZ(op.q1, op.q2) +LinearAlgebra.inv(op::sInvSQRTZZ) = sSQRTZZ(op.q1, op.q2) +LinearAlgebra.inv(op::sSQRTXX) = sInvSQRTXX(op.q1, op.q2) +LinearAlgebra.inv(op::sInvSQRTXX) = sSQRTXX(op.q1, op.q2) +LinearAlgebra.inv(op::sSQRTYY) = sInvSQRTYY(op.q1, op.q2) +LinearAlgebra.inv(op::sInvSQRTYY) = sSQRTYY(op.q1, op.q2) ############################## # Functions that perform direct application of common operators without needing an operator instance @@ -394,12 +425,9 @@ LinearAlgebra.inv(op::sInvZCrY) = sZCrY(op.q1, op.q2) """Apply a Pauli Z to the `i`-th qubit of state `s`. You should use `apply!(stab,sZ(i))` instead of this.""" function apply_single_z!(stab::AbstractStabilizer, i) s = tab(stab) - Tₘₑ = eltype(s.xzs) - bigi = _div(Tₘₑ,i-1)+1 - smalli = _mod(Tₘₑ,i-1) - mask = Tₘₑ(0x1)<