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 233156167..06856e122 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,15 +5,30 @@ # News -## dev +## v0.9.13 - dev - 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` + +## v0.9.11 - 2024-09-27 + +- `hcat` of Tableaux objects +- `QuantumReedMuller` codes added to the ECC module +- **(breaking)** change the convention for how to provide a representation function in the constructor of `LPCode` -- strictly speaking a breaking change, but this is not an API that is publicly used in practice + +## v0.9.10 - 2024-09-26 + +- The lifted product class of quantum LDPC codes is implemented in the ECC submodule. +- **(fix)** `ECC.code_s` now gives the number of parity checks with redundancy. If you want the number of linearly independent parity checks, you can use `LinearAlgebra.rank`. - Implementing many more named single-qubit gates following naming convention similar to the stim package in python. - **(fix)** Bug fix to the `parity_checks(ReedMuller(r, m))` of classical Reed-Muller code (it was returning generator matrix). - `RecursiveReedMuller` code implementation as an alternative implementation of `ReedMuller`. + ## v0.9.9 - 2024-08-05 - `inv` is implemented for all Clifford operator types (symbolic, dense, sparse). diff --git a/Project.toml b/Project.toml index a5dc91bfb..4d300bc66 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.9" +version = "0.9.12" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -24,6 +24,7 @@ SumTypes = "8e1ec7a9-0e02-4297-b0fe-6433085c89f2" [weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -33,6 +34,7 @@ QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" [extensions] QuantumCliffordGPUExt = "CUDA" +QuantumCliffordHeckeExt = "Hecke" QuantumCliffordLDPCDecodersExt = "LDPCDecoders" QuantumCliffordMakieExt = "Makie" QuantumCliffordPlotsExt = "Plots" @@ -46,6 +48,7 @@ Combinatorics = "1.0" DataStructures = "0.18" 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" InteractiveUtils = "1.9" @@ -53,7 +56,7 @@ LDPCDecoders = "0.3.1" LinearAlgebra = "1.9" MacroTools = "0.5.9" Makie = "0.20, 0.21" -Nemo = "0.42, 0.43, 0.44, 0.45, 0.46" +Nemo = "0.42.1, 0.43, 0.44, 0.45, 0.46, 0.47" Plots = "1.38.0" PrecompileTools = "1.2" PyQDecoders = "0.2.1" diff --git a/docs/Project.toml b/docs/Project.toml index b414301a3..8ef02232a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,6 +5,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb" Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/make.jl b/docs/make.jl index 699924cda..b8d13fc30 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,6 +7,11 @@ using QuantumClifford using QuantumClifford.Experimental.NoisyCircuits using QuantumInterface +ENV["HECKE_PRINT_BANNER"] = "false" +import Hecke + +const QuantumCliffordHeckeExt = Base.get_extension(QuantumClifford, :QuantumCliffordHeckeExt) + #DocMeta.setdocmeta!(QuantumClifford, :DocTestSetup, :(using QuantumClifford); recursive=true) ENV["LINES"] = 80 # for forcing `displaysize(io)` to be big enough @@ -20,8 +25,9 @@ doctest = false, clean = true, sitename = "QuantumClifford.jl", format = Documenter.HTML(size_threshold_ignore = ["API.md"]), -modules = [QuantumClifford, QuantumClifford.Experimental.NoisyCircuits, QuantumClifford.ECC, QuantumInterface], +modules = [QuantumClifford, QuantumClifford.Experimental.NoisyCircuits, QuantumClifford.ECC, QuantumInterface, QuantumCliffordHeckeExt], warnonly = [:missing_docs], +linkcheck = true, authors = "Stefan Krastanov", pages = [ "QuantumClifford.jl" => "index.md", diff --git a/docs/src/ECC_API.md b/docs/src/ECC_API.md index d77b26097..c400a2bf4 100644 --- a/docs/src/ECC_API.md +++ b/docs/src/ECC_API.md @@ -4,3 +4,10 @@ Modules = [QuantumClifford.ECC] Private = false ``` + +## Implemented in an extension requiring `Hecke.jl` + +```@autodocs +Modules = [QuantumCliffordHeckeExt] +Private = true +``` \ No newline at end of file diff --git a/docs/src/references.bib b/docs/src/references.bib index 6a50a9194..29500a034 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -402,3 +402,88 @@ @inproceedings{brown2013short pages = {346--350}, doi = {10.1109/ISIT.2013.6620245} } + +@article{panteleev2021degenerate, + title = {Degenerate {{Quantum LDPC Codes With Good Finite Length Performance}}}, + author = {Panteleev, Pavel and Kalachev, Gleb}, + year = {2021}, + month = nov, + journal = {Quantum}, + volume = {5}, + eprint = {1904.02703}, + primaryclass = {quant-ph}, + pages = {585}, + issn = {2521-327X}, + doi = {10.22331/q-2021-11-22-585}, + archiveprefix = {arXiv} +} + + +@inproceedings{panteleev2022asymptotically, + title = {Asymptotically Good {{Quantum}} and Locally Testable Classical {{LDPC}} Codes}, + booktitle = {Proceedings of the 54th {{Annual ACM SIGACT Symposium}} on {{Theory}} of {{Computing}}}, + author = {Panteleev, Pavel and Kalachev, Gleb}, + year = {2022}, + month = jun, + pages = {375--388}, + publisher = {ACM}, + address = {Rome Italy}, + doi = {10.1145/3519935.3520017}, + isbn = {978-1-4503-9264-8} +} + +@article{roffe2023bias, + title = {Bias-Tailored Quantum {{LDPC}} Codes}, + author = {Roffe, Joschka and Cohen, Lawrence Z. and Quintavalle, Armanda O. and Chandra, Daryus and Campbell, Earl T.}, + year = {2023}, + month = may, + journal = {Quantum}, + volume = {7}, + pages = {1005}, + doi = {10.22331/q-2023-05-15-1005} +} + +@article{raveendran2022finite, + title = {Finite {{Rate QLDPC-GKP Coding Scheme}} That {{Surpasses}} the {{CSS Hamming Bound}}}, + author = {Raveendran, Nithin and Rengaswamy, Narayanan and Rozp{\k e}dek, Filip and Raina, Ankur and Jiang, Liang and Vasi{\'c}, Bane}, + year = {2022}, + month = jul, + journal = {Quantum}, + volume = {6}, + pages = {767}, + issn = {2521-327X}, + doi = {10.22331/q-2022-07-20-767}, +} + +@article{steane1999quantum, + title={Quantum reed-muller codes}, + author={Steane, Andrew M}, + journal={IEEE Transactions on Information Theory}, + volume={45}, + number={5}, + pages={1701--1703}, + year={1999}, + publisher={IEEE} +} + +@article{campbell2012magic, + title={Magic-state distillation in all prime dimensions using quantum reed-muller codes}, + author={Campbell, Earl T and Anwar, Hussain and Browne, Dan E}, + journal={Physical Review X}, + volume={2}, + number={4}, + pages={041021}, + year={2012}, + publisher={APS} +} + +@article{anderson2014fault, + title={Fault-tolerant conversion between the steane and reed-muller quantum codes}, + author={Anderson, Jonas T and Duclos-Cianci, Guillaume and Poulin, David}, + journal={Physical review letters}, + volume={113}, + number={8}, + pages={080501}, + year={2014}, + publisher={APS} +} diff --git a/docs/src/references.md b/docs/src/references.md index b37e94d4e..35e944a21 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -37,6 +37,9 @@ For quantum code construction routines: - [kitaev2003fault](@cite) - [fowler2012surface](@cite) - [knill1996concatenated](@cite) +- [steane1999quantum](@cite) +- [campbell2012magic](@cite) +- [anderson2014fault](@cite) For classical code construction routines: - [muller1954application](@cite) diff --git a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl new file mode 100644 index 000000000..29e9de8ce --- /dev/null +++ b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl @@ -0,0 +1,20 @@ +module QuantumCliffordHeckeExt + +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 +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 + +include("types.jl") +include("lifted.jl") +include("lifted_product.jl") + +end # module diff --git a/ext/QuantumCliffordHeckeExt/lifted.jl b/ext/QuantumCliffordHeckeExt/lifted.jl new file mode 100644 index 000000000..75964f983 --- /dev/null +++ b/ext/QuantumCliffordHeckeExt/lifted.jl @@ -0,0 +1,105 @@ +""" +$TYPEDEF + +Classical codes lifted over a group algebra, used for lifted product code construction ([panteleev2021degenerate](@cite), [panteleev2022asymptotically](@cite)) + +The parity-check matrix is constructed by applying `repr` to each element of `A`, +which is mathematically a linear map from a group algebra element to a binary matrix. +The size of the parity check matrix will enlarged with each element of `A` being inflated into a matrix. +The procedure is called a lift [panteleev2022asymptotically](@cite). + +## Constructors + +A lifted code can be constructed via the following approaches: + +1. A matrix of group algebra elements. + +2. A matrix of group elements, where a group element will be considered as a group algebra element by assigning a unit coefficient. + +3. A matrix of integers, where each integer represent the shift of a cyclic permutation. The order of the cyclic permutation should be specified. + +The default `GA` is the group algebra of `A[1, 1]`, the default representation `repr` is the permutation representation. + +## The representation function `repr` + +We use the default representation function `Hecke.representation_matrix` to convert a `GF(2)`-group algebra element to a binary matrix. +The default representation, provided by `Hecke`, is the permutation representation. + +We also accept a custom representation function (the `repr` field of the constructor). +Whatever the representation, the matrix elements need to be convertible to Integers (e.g. permit `lift(ZZ, ...)`). +Such a customization would be useful to reduce the number of bits required by the code construction. + +For example, if we use a D4 group for lifting, our default representation will be `8×8` permutation matrices, +where 8 is the group's order. +However, we can find a `4×4` matrix representation for the group, +e.g. by using the typical [`2×2` representation](https://en.wikipedia.org/wiki/Dihedral_group) +and converting it into binary representation by replacing "1" with the Pauli I, and "-1" with the Pauli X matrix. + +See also: [`LPCode`](@ref). + +$TYPEDFIELDS +""" +struct LiftedCode <: ClassicalCode + """the base matrix of the code, whose elements are in a group algebra.""" + A::GroupAlgebraElemMatrix + """the group algebra for which elements in `A` are from.""" + GA::GroupAlgebra + """ + a function that converts a group algebra element to a binary matrix; + default to be the permutation representation for GF(2)-algebra.""" + repr::Function + + function LiftedCode(A::GroupAlgebraElemMatrix; GA::GroupAlgebra=parent(A[1, 1]), repr::Function) + all(elem.parent == GA for elem in A) || error("The base ring of all elements in the code must be the same as the group algebra") + new(A, GA, repr) + end +end + +""" +`LiftedCode` constructor using the default `GF(2)` representation (coefficients converted to a permutation matrix by `representation_matrix` provided by Hecke). +""" # TODO doctest example +function LiftedCode(A::Matrix{GroupAlgebraElem{FqFieldElem, <: GroupAlgebra}}; GA::GroupAlgebra=parent(A[1,1])) + !(characteristic(base_ring(A[1, 1])) == 2) && error("The default permutation representation applies only to GF(2) group algebra; otherwise, a custom representation function should be provided") + LiftedCode(A; GA=GA, repr=representation_matrix) +end + +# TODO document and doctest example +function LiftedCode(group_elem_array::Matrix{<: GroupOrAdditiveGroupElem}; GA::GroupAlgebra=group_algebra(GF(2), parent(group_elem_array[1,1])), repr::Union{Function, Nothing}=nothing) + A = zeros(GA, size(group_elem_array)...) + for i in axes(group_elem_array, 1), j in axes(group_elem_array, 2) + A[i, j] = GA[A[i, j]] + end + if repr === nothing + return LiftedCode(A; GA=GA, repr=representation_matrix) + else + return LiftedCode(A; GA=GA, repr=repr) + end +end + +# TODO document and doctest example +function LiftedCode(shift_array::Matrix{Int}, l::Int; GA::GroupAlgebra=group_algebra(GF(2), abelian_group(l))) + A = zeros(GA, size(shift_array)...) + for i in 1:size(shift_array, 1) + for j in 1:size(shift_array, 2) + A[i, j] = GA[shift_array[i, j]%l+1] + end + end + return LiftedCode(A; GA=GA, repr=representation_matrix) +end + +lift_to_bool(x) = Bool(Int(lift(ZZ,x))) + +function concat_lift_repr(repr, mat) + x = repr.(mat) + y = hvcat(size(x,2), transpose(x)...) + z = Matrix(lift_to_bool.(y)) + return z +end + +function parity_checks(c::LiftedCode) + return lift(c.repr, c.A) +end + +code_n(c::LiftedCode) = size(c.A, 2) * size(zero(c.GA), 2) + +code_s(c::LiftedCode) = size(c.A, 1) * size(zero(c.GA), 1) diff --git a/ext/QuantumCliffordHeckeExt/lifted_product.jl b/ext/QuantumCliffordHeckeExt/lifted_product.jl new file mode 100644 index 000000000..ef09eddbb --- /dev/null +++ b/ext/QuantumCliffordHeckeExt/lifted_product.jl @@ -0,0 +1,198 @@ +""" +$TYPEDEF + +Lifted product codes ([panteleev2021degenerate](@cite), [panteleev2022asymptotically](@cite)) + +A lifted product code is defined by the hypergraph product of a base matrices `A` and the conjugate of another base matrix `B'`. +Here, the hypergraph product is taken over a group algebra, of which the base matrices are consisting. + +The binary parity check matrix is obtained by applying `repr` to each element of the matrix resulted from the hypergraph product, which is mathematically a linear map from each group algebra element to a binary matrix. + +## Constructors + +Multiple constructors are available: + +1. Two base matrices of group algebra elements. + +2. Two lifted codes, whose base matrices are for quantum code construction. + +3. Two base matrices of group elements, where each group element will be considered as a group algebra element by assigning a unit coefficient. + +4. Two base matrices of integers, where each integer represent the shift of a cyclic permutation. The order of the cyclic permutation should be specified. + +## Examples + +A [[882, 24, d ≤ 24]] code from Appendix B of [roffe2023bias](@cite). +We use the 1st constructor to generate the code and check its length and dimension. +During the construction, we do arithmetic operations to get the group algebra elements in base matrices `A` and `B`. +Here `x` is the generator of the group algebra, i.e., offset-1 cyclic permutation, and `GA(1)` is the unit element. + +```jldoctest +julia> import Hecke: group_algebra, GF, abelian_group, gens; import LinearAlgebra: diagind; + +julia> l = 63; GA = group_algebra(GF(2), abelian_group(l)); x = gens(GA)[]; + +julia> A = zeros(GA, 7, 7); + +julia> A[diagind(A)] .= x^27; + +julia> A[diagind(A, -1)] .= x^54; + +julia> A[diagind(A, 6)] .= x^54; + +julia> A[diagind(A, -2)] .= GA(1); + +julia> A[diagind(A, 5)] .= GA(1); + +julia> B = reshape([1 + x + x^6], (1, 1)); + +julia> c1 = LPCode(A, B); + +julia> code_n(c1), code_k(c1) +(882, 24) +``` + +A [[175, 19, d ≤ 0]] code from Eq. (18) in Appendix A of [raveendran2022finite](@cite), +following the 4th constructor. + +```jldoctest +julia> base_matrix = [0 0 0 0; 0 1 2 5; 0 6 3 1]; l = 7; + +julia> c2 = LPCode(base_matrix, l .- base_matrix', l); + +julia> code_n(c2), code_k(c2) +(175, 19) +``` + +## Code subfamilies and convenience constructors for them + +- When the base matrices of the `LPCode` are 1×1, the code is called a two-block group-algebra code [`two_block_group_algebra_codes`](@ref). +- When the base matrices of the `LPCode` are 1×1 and their elements are sums of cyclic permutations, the code is called a generalized bicycle code [`generalized_bicycle_codes`](@ref). +- When the two matrices are adjoint to each other, the code is called a bicycle code [`bicycle_codes`](@ref). + +## The representation function + +We use the default representation function `Hecke.representation_matrix` to convert a `GF(2)`-group algebra element to a binary matrix. +The default representation, provided by `Hecke`, is the permutation representation. + +We also accept a custom representation function as detailed in [`LiftedCode`](@ref). + +See also: [`LiftedCode`](@ref), [`two_block_group_algebra_codes`](@ref), [`generalized_bicycle_codes`](@ref), [`bicycle_codes`](@ref). + +$TYPEDFIELDS +""" +struct LPCode <: AbstractECC + """the first base matrix of the code, whose elements are in a group algebra.""" + A::GroupAlgebraElemMatrix + """the second base matrix of the code, whose elements are in the same group algebra as `A`.""" + B::GroupAlgebraElemMatrix + """the group algebra for which elements in `A` and `B` are from.""" + GA::GroupAlgebra + """ + a function that converts a group algebra element to a binary matrix; + default to be the permutation representation for GF(2)-algebra.""" + repr::Function + + function LPCode(A::GroupAlgebraElemMatrix, B::GroupAlgebraElemMatrix; GA::GroupAlgebra=parent(A[1,1]), repr::Function) + all(elem.parent == GA for elem in A) && all(elem.parent == GA for elem in B) || error("The base rings of all elements in both matrices must be the same as the group algebra") + new(A, B, GA, repr) + end + + function LPCode(c₁::LiftedCode, c₂::LiftedCode; GA::GroupAlgebra=c₁.GA, repr::Function=c₁.repr) + # we are using the group algebra and the representation function of the first lifted code + c₁.GA == GA && c₂.GA == GA || error("The base rings of both lifted codes must be the same as the group algebra") + new(c₁.A, c₂.A, GA, repr) + end +end + +# TODO document and doctest example +function LPCode(A::FqFieldGroupAlgebraElemMatrix, B::FqFieldGroupAlgebraElemMatrix; GA::GroupAlgebra=parent(A[1,1])) + LPCode(LiftedCode(A; GA=GA, repr=representation_matrix), LiftedCode(B; GA=GA, repr=representation_matrix); GA=GA, repr=representation_matrix) +end + +# TODO document and doctest example +function LPCode(group_elem_array1::Matrix{<: GroupOrAdditiveGroupElem}, group_elem_array2::Matrix{<: GroupOrAdditiveGroupElem}; GA::GroupAlgebra=group_algebra(GF(2), parent(group_elem_array1[1,1]))) + LPCode(LiftedCode(group_elem_array1; GA=GA), LiftedCode(group_elem_array2; GA=GA); GA=GA, repr=representation_matrix) +end + +# TODO document and doctest example +function LPCode(shift_array1::Matrix{Int}, shift_array2::Matrix{Int}, l::Int; GA::GroupAlgebra=group_algebra(GF(2), abelian_group(l))) + LPCode(LiftedCode(shift_array1, l; GA=GA), LiftedCode(shift_array2, l; GA=GA); GA=GA, repr=representation_matrix) +end + +iscss(::Type{LPCode}) = true + +function hgp(h₁::GroupAlgebraElemMatrix, h₂::GroupAlgebraElemMatrix) + r₁, n₁ = size(h₁) + r₂, n₂ = size(h₂) + # here we use `permutdims` instead of `transpose` to avoid recursive call + # convert LinearAlgebra.I to Matrix to fix incompatibility with Julia 1.11.1 + # TODO the performance may be affected by this workaround for large codes + hx = hcat(kron(h₁, Matrix(LinearAlgebra.I(n₂))), kron(Matrix(LinearAlgebra.I(r₁)), permutedims(group_algebra_conj.(h₂)))) + hz = hcat(kron(Matrix(LinearAlgebra.I(n₁)), h₂), kron(permutedims(group_algebra_conj.(h₁)), Matrix(LinearAlgebra.I(r₂)))) + hx, hz +end + +function parity_checks_xz(c::LPCode) + hx, hz = hgp(c.A, permutedims(group_algebra_conj.(c.B))) + hx, hz = concat_lift_repr(c.repr,hx), concat_lift_repr(c.repr,hz) + return hx, hz +end + +parity_checks_x(c::LPCode) = parity_checks_xz(c)[1] + +parity_checks_z(c::LPCode) = parity_checks_xz(c)[2] + +parity_checks(c::LPCode) = parity_checks(CSS(parity_checks_xz(c)...)) + +code_n(c::LPCode) = size(c.repr(zero(c.GA)), 2) * (size(c.A, 2) * size(c.B, 1) + size(c.A, 1) * size(c.B, 2)) + +code_s(c::LPCode) = size(c.repr(zero(c.GA)), 1) * (size(c.A, 1) * size(c.B, 1) + size(c.A, 2) * size(c.B, 2)) + +""" +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. + +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) +end + +""" +Generalized bicycle codes, which are a special case of 2GBA codes (and therefore of lifted product codes). +Here the group is chosen as the cyclic group of order `l`, +and the base matrices `a` and `b` are the sum of the group algebra elements corresponding to the shifts `a_shifts` and `b_shifts`. + +See also: [`two_block_group_algebra_codes`](@ref), [`bicycle_codes`](@ref). + +A [[254, 28, 14 ≤ d ≤ 20]] code from (A1) in Appendix B of [panteleev2021degenerate](@cite). + +```jldoctest +julia> c = generalized_bicycle_codes([0, 15, 20, 28, 66], [0, 58, 59, 100, 121], 127); + +julia> code_n(c), code_k(c) +(254, 28) +``` +""" +function generalized_bicycle_codes(a_shifts::Array{Int}, b_shifts::Array{Int}, l::Int) + GA = group_algebra(GF(2), abelian_group(l)) + a = sum(GA[n%l+1] for n in a_shifts) + b = sum(GA[n%l+1] for n in b_shifts) + two_block_group_algebra_codes(a, b) +end + +""" +Bicycle codes are a special case of generalized bicycle codes, +where `a` and `b` are conjugate to each other. +The order of the cyclic group is `l`, and the shifts `a_shifts` and `b_shifts` are reverse to each other. + +See also: [`two_block_group_algebra_codes`](@ref), [`generalized_bicycle_codes`](@ref). +""" # TODO doctest example +function bicycle_codes(a_shifts::Array{Int}, l::Int) + GA = group_algebra(GF(2), abelian_group(l)) + a = sum(GA[n÷l+1] for n in a_shifts) + two_block_group_algebra_codes(a, group_algebra_conj(a)) +end diff --git a/ext/QuantumCliffordHeckeExt/types.jl b/ext/QuantumCliffordHeckeExt/types.jl new file mode 100644 index 000000000..ecd0d7a68 --- /dev/null +++ b/ext/QuantumCliffordHeckeExt/types.jl @@ -0,0 +1,35 @@ +const GroupOrAdditiveGroupElem = Union{GroupElem,AdditiveGroupElem} + +const GroupAlgebraElemMatrix = Union{ + Matrix{<:GroupAlgebraElem}, + LinearAlgebra.Adjoint{<:GroupAlgebraElem,<:Matrix{<:GroupAlgebraElem}} +} + +const FqFieldGroupAlgebraElemMatrix = Union{ + Matrix{<:GroupAlgebraElem{FqFieldElem,<:GroupAlgebra}}, + LinearAlgebra.Adjoint{<:GroupAlgebraElem{FqFieldElem,<:GroupAlgebra},<:Matrix{<:GroupAlgebraElem{FqFieldElem,<:GroupAlgebra}}} +} + +""" +Compute the conjugate of a group algebra element. +The conjugate is defined by inversing elements in the associated group. +""" +function group_algebra_conj(a::GroupAlgebraElem{T}) where T + A = parent(a) + d = dim(A) + v = Vector{T}(undef, d) + for i in 1:d + v[i] = zero(base_ring(A)) + end + id_index = findfirst(x -> x == 1, one(A).coeffs) + # t = zero(base_ring(A)) + mt = multiplication_table(A, copy = false) + acoeff = coefficients(a, copy = false) + for i in 1:d + if acoeff[i] != 0 + k = findfirst(x -> x==id_index, mt[i, :]) # find the inverse of i-th element in the group + v[k] += acoeff[i] + end + end + return A(v) +end diff --git a/ext/QuantumCliffordLDPCDecodersExt/QuantumCliffordLDPCDecodersExt.jl b/ext/QuantumCliffordLDPCDecodersExt/QuantumCliffordLDPCDecodersExt.jl index 8c369cf82..e80c13af6 100644 --- a/ext/QuantumCliffordLDPCDecodersExt/QuantumCliffordLDPCDecodersExt.jl +++ b/ext/QuantumCliffordLDPCDecodersExt/QuantumCliffordLDPCDecodersExt.jl @@ -74,16 +74,16 @@ parity_checks(d::BeliefPropDecoder) = d.H parity_checks(d::BitFlipDecoder) = d.H function decode(d::BeliefPropDecoder, syndrome_sample) - row_x = syndrome_sample[1:d.cx] - row_z = syndrome_sample[d.cx+1:d.cx+d.cz] + row_x = @view syndrome_sample[1:d.cx] + row_z = @view syndrome_sample[d.cx+1:d.cx+d.cz] guess_z, success = LDPCDecoders.decode!(d.bpdecoderx, row_x) guess_x, success = LDPCDecoders.decode!(d.bpdecoderz, row_z) return vcat(guess_x, guess_z) end function decode(d::BitFlipDecoder, syndrome_sample) - row_x = syndrome_sample[1:d.cx] - row_z = syndrome_sample[d.cx+1:d.cx+d.cz] + row_x = @view syndrome_sample[1:d.cx] + row_z = @view syndrome_sample[d.cx+1:d.cx+d.cz] guess_z, success = LDPCDecoders.decode!(d.bfdecoderx, row_x) guess_x, success = LDPCDecoders.decode!(d.bfdecoderz, row_z) return vcat(guess_x, guess_z) diff --git a/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl b/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl index bfa557d05..e3496a733 100644 --- a/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl +++ b/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl @@ -69,8 +69,8 @@ end parity_checks(d::PyBP) = d.H function decode(d::PyBP, syndrome_sample) - row_x = syndrome_sample[1:d.nx] # TODO These copies and indirections might be costly! - row_z = syndrome_sample[d.nx+1:end] + row_x = @view syndrome_sample[1:d.nx] + row_z = @view syndrome_sample[d.nx+1:end] guess_z_errors = PythonCall.PyArray(d.pyx.decode(np.array(row_x))) guess_x_errors = PythonCall.PyArray(d.pyz.decode(np.array(row_z))) vcat(guess_x_errors, guess_z_errors) @@ -106,18 +106,19 @@ end parity_checks(d::PyMatchingDecoder) = d.H function decode(d::PyMatchingDecoder, syndrome_sample) - row_x = syndrome_sample[1:d.nx] # TODO This copy is costly! - row_z = syndrome_sample[d.nx+1:end] + row_x = @view syndrome_sample[1:d.nx] + row_z = @view syndrome_sample[d.nx+1:end] guess_z_errors = PythonCall.PyArray(d.pyx.decode(row_x)) guess_x_errors = PythonCall.PyArray(d.pyz.decode(row_z)) vcat(guess_x_errors, guess_z_errors) end function batchdecode(d::PyMatchingDecoder, syndrome_samples) - row_x = syndrome_samples[:,1:d.nx] # TODO This copy is costly! - row_z = syndrome_samples[:,d.nx+1:end] + row_x = @view syndrome_samples[:,1:d.nx] + row_z = @view syndrome_samples[:,d.nx+1:end] guess_z_errors = PythonCall.PyArray(d.pyx.decode_batch(row_x)) guess_x_errors = PythonCall.PyArray(d.pyz.decode_batch(row_z)) + n_cols_x = size(guess_x_errors, 2) hcat(guess_x_errors, guess_z_errors) end diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index 4137b00a3..6837483cb 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 @@ -174,7 +174,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 @@ -230,18 +230,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 @@ -379,7 +379,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)) @@ -498,11 +498,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 @@ -581,7 +581,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 @@ -589,7 +589,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 @@ -600,9 +600,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 @@ -611,17 +611,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)) @@ -936,6 +936,19 @@ function check_allrowscommute(stabilizer::Tableau) end check_allrowscommute(stabilizer::Stabilizer)=check_allrowscommute(tab(stabilizer)) +""" +Vertically concatenates tableaux. + +```jldoctest +julia> vcat(ghz(2), ghz(2)) ++ XX ++ ZZ ++ XX ++ ZZ +``` + +See also: [`hcat`](@ref) +""" function Base.vcat(tabs::Tableau...) Tableau( vcat((s.phases for s in tabs)...), @@ -943,7 +956,41 @@ function Base.vcat(tabs::Tableau...) hcat((s.xzs for s in tabs)...)) end -Base.vcat(stabs::Stabilizer...) = Stabilizer(vcat((tab(s) for s in stabs)...)) +Base.vcat(stabs::Stabilizer{T}...) where {T} = Stabilizer(vcat((tab(s) for s in stabs)...)) + +""" +Horizontally concatenates tableaux. + +```jldoctest +julia> hcat(ghz(2), ghz(2)) ++ XXXX ++ ZZZZ +``` + +See also: [`vcat`](@ref) +""" +function Base.hcat(tabs::Tableau...) # TODO this implementation is slow as it unpacks each bitvector into bits and repacks them -- reuse the various tableau inset functionality we have to speed this up + rows = size(tabs[1], 1) + cols = sum(map(nqubits, tabs)) + newtab = zero(Tableau, rows, cols) + cols_idx = 1 + for tab in tabs + rows_tab, cols_tab = size(tab) + if rows_tab != rows + throw(ArgumentError("All input Tableaux/Stabilizers must have the same number of rows.")) + end + for i in 1:rows + for j in 1:cols_tab + newtab[i, cols_idx+j-1]::Tuple{Bool,Bool} = tab[i, j]::Tuple{Bool,Bool} + end + newtab.phases[i] = (newtab.phases[i]+tab.phases[i])%4 + end + cols_idx += cols_tab + end + return newtab +end + +Base.hcat(stabs::Stabilizer{T}...) where {T} = Stabilizer(hcat((tab(s) for s in stabs)...)) ############################## # Unitary Clifford Operations @@ -990,6 +1037,16 @@ function _apply!(stab::AbstractStabilizer, p::PauliOperator, indices; phases::Va stab end +############################## +# Conversion and promotion +############################## + +Base.promote_rule(::Type{<:Destabilizer{T}} , ::Type{<:MixedDestabilizer{T}}) where {T<:Tableau} = MixedDestabilizer{T} +Base.promote_rule(::Type{<:MixedStabilizer{T}}, ::Type{<:MixedDestabilizer{T}}) where {T<:Tableau} = MixedDestabilizer{T} +Base.promote_rule(::Type{<:Stabilizer{T}} , ::Type{<:S} ) where {T<:Tableau, S<:Union{MixedStabilizer{T}, Destabilizer{T}, MixedDestabilizer{T}}} = S + +Base.convert(::Type{<:MixedDestabilizer{T}}, x::Union{Destabilizer{T}, MixedStabilizer{T}, Stabilizer{T}}) where {T <: Tableau} = MixedDestabilizer(x) + ############################## # Helpers for binary codes ############################## @@ -1102,8 +1159,14 @@ end """ Check basic consistency requirements of a stabilizer. Used in tests. """ -function stab_looks_good(s) - c = tab(canonicalize!(copy(s))) +function stab_looks_good(s; remove_redundant_rows=false) + # first remove redundant rows + c = if remove_redundant_rows + s1, _, rank = canonicalize!(copy(s), ranks=true) + tab(s1[1:rank]) + else + tab(canonicalize!(copy(s))) + end nrows, ncols = size(c) all((c.phases .== 0x0) .| (c.phases .== 0x2)) || return false H = stab_to_gf2(c) 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 955ec5818..cdda7742e 100644 --- a/src/ecc/ECC.jl +++ b/src/ecc/ECC.jl @@ -17,10 +17,11 @@ export parity_checks, parity_checks_x, parity_checks_z, iscss, code_n, code_s, code_k, rate, distance, isdegenerate, faults_matrix, naive_syndrome_circuit, shor_syndrome_circuit, naive_encoding_circuit, - RepCode, + RepCode, LiftedCode, CSS, Shor9, Steane7, Cleve8, Perfect5, Bitflip3, - Toric, Gottesman, Surface, Concat, CircuitCode, + Toric, Gottesman, Surface, Concat, CircuitCode, QuantumReedMuller, + LPCode, two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes, random_brickwork_circuit_code, random_all_to_all_circuit_code, evaluate_decoder, CommutationCheckECCSetup, NaiveSyndromeECCSetup, ShorSyndromeECCSetup, @@ -87,14 +88,22 @@ code_n(c::AbstractECC) = code_n(parity_checks(c)) code_n(s::Stabilizer) = nqubits(s) -"""The number of stabilizer checks in a code.""" +"""The number of stabilizer checks in a code. They might not be all linearly independent, thus `code_s >= code_n-code_k`. For the number of linearly independent checks you can use `LinearAlgebra.rank`.""" function code_s end - -code_s(c::AbstractECC) = code_s(parity_checks(c)) code_s(s::Stabilizer) = length(s) +code_s(c::AbstractECC) = code_s(parity_checks(c)) -"""The number of logical qubits in a code.""" -code_k(c) = code_n(c) - code_s(c) +""" +The number of logical qubits in a code. + +Note that when redundant rows exist in the parity check matrix, the number of logical qubits `code_k(c)` will be greater than `code_n(c) - code_s(c)`, where the difference equals the redundancy. +""" +function code_k(s::Stabilizer) + _, _, r = canonicalize!(Base.copy(s), ranks=true) + return code_n(s) - r +end + +code_k(c::AbstractECC) = code_k(parity_checks(c)) """The rate of a code.""" function rate(c) @@ -354,6 +363,7 @@ include("circuits.jl") include("decoder_pipeline.jl") include("codes/util.jl") + include("codes/classical_codes.jl") include("codes/css.jl") include("codes/bitflipcode.jl") @@ -367,6 +377,12 @@ include("codes/surface.jl") include("codes/concat.jl") include("codes/random_circuit.jl") include("codes/classical/reedmuller.jl") -include("codes/classical/bch.jl") include("codes/classical/recursivereedmuller.jl") +include("codes/classical/bch.jl") +include("codes/quantumreedmuller.jl") + +# qLDPC +include("codes/classical/lifted.jl") +include("codes/lifted_product.jl") + end #module 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/classical/lifted.jl b/src/ecc/codes/classical/lifted.jl new file mode 100644 index 000000000..a7c20fa04 --- /dev/null +++ b/src/ecc/codes/classical/lifted.jl @@ -0,0 +1,10 @@ +"""Classical codes lifted over a group algebra, used for lifted product code construction ([panteleev2021degenerate](@cite), [panteleev2022asymptotically](@cite)) + +Implemented as a package extension with Hecke. Check the [QuantumClifford documentation](http://qc.quantumsavory.org/stable/ECC_API/) for more details on that extension.""" +function LiftedCode(args...; kwargs...) + ext = Base.get_extension(QuantumClifford, :QuantumCliffordHeckeExt) + if isnothing(ext) + throw("The `LiftedCode` depends on the package `Hecke` but you have not installed or imported it yet. Immediately after you import `Hecke`, the `LiftedCode` will be available.") + end + return ext.LiftedCode(args...; kwargs...) +end diff --git a/src/ecc/codes/lifted_product.jl b/src/ecc/codes/lifted_product.jl new file mode 100644 index 000000000..338880702 --- /dev/null +++ b/src/ecc/codes/lifted_product.jl @@ -0,0 +1,19 @@ +"""Lifted product codes ([panteleev2021degenerate](@cite), [panteleev2022asymptotically](@cite)) + +Implemented as a package extension with Hecke. Check the [QuantumClifford documentation](http://qc.quantumsavory.org/stable/ECC_API/) for more details on that extension.""" +function LPCode(args...; kwargs...) + ext = Base.get_extension(QuantumClifford, :QuantumCliffordHeckeExt) + if isnothing(ext) + throw("The `LPCode` depends on the package `Hecke` but you have not installed or imported it yet. Immediately after you import `Hecke`, the `LPCode` will be available.") + end + return ext.LPCode(args...; kwargs...) +end + +"""Implemented in a package extension with Hecke.""" +function two_block_group_algebra_codes end + +"""Implemented in a package extension with Hecke.""" +function generalized_bicycle_codes end + +"""Implemented in a package extension with Hecke.""" +function bicycle_codes end diff --git a/src/ecc/codes/quantumreedmuller.jl b/src/ecc/codes/quantumreedmuller.jl new file mode 100644 index 000000000..6813d6646 --- /dev/null +++ b/src/ecc/codes/quantumreedmuller.jl @@ -0,0 +1,41 @@ +""" +The family of `[[2ᵐ - 1, 1, 3]]` CSS Quantum-Reed-Muller codes, as discovered by Steane in his 1999 paper [steane1999quantum](@cite). + +Quantum codes are constructed from shortened Reed-Muller codes `RM(1, m)`, by removing the first row and column of the generator matrix `Gₘ`. Similarly, we can define truncated dual codes `RM(m - 2, m)` using the generator matrix `Hₘ` [anderson2014fault](@cite). The quantum Reed-Muller codes `QRM(m)` derived from `RM(1, m)` are CSS codes. + +Given that the stabilizers of the quantum code are defined through the generator matrix of the classical code, the minimum distance of the quantum code corresponds to the minimum distance of the dual classical code, which is `d = 3`, thus it can correct any single qubit error. Since one stabilizer from the original and one from the dual code are removed in the truncation process, the code parameters are `[[2ᵐ - 1, 1, 3]]`. + +You might be interested in consulting [anderson2014fault](@cite) and [campbell2012magic](@cite) as well. + +The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/quantum_reed_muller). +""" +struct QuantumReedMuller <: AbstractECC + m::Int + function QuantumReedMuller(m) + if m < 3 + throw(DomainError("Invalid parameters: m must be bigger than 2 in order to have a valid code.")) + end + new(m) + end +end + +function iscss(::Type{QuantumReedMuller}) + return true +end + +function parity_checks(c::QuantumReedMuller) + RM₁₋ₘ = generator(RecursiveReedMuller(1, c.m)) + RM₍ₘ₋₂₎₋ₘ₎ = generator(RecursiveReedMuller(c.m-2, c.m)) + QRM = CSS(RM₁₋ₘ[2:end, 2:end], RM₍ₘ₋₂₎₋ₘ₎[2:end, 2:end]) + Stabilizer(QRM) +end + +code_n(c::QuantumReedMuller) = 2^c.m - 1 + +code_k(c::QuantumReedMuller) = 1 + +distance(c::QuantumReedMuller) = 3 + +parity_checks_x(c::QuantumReedMuller) = stab_to_gf2(parity_checks(QuantumReedMuller(c.m)))[1:c.m, 1:end÷2] + +parity_checks_z(c::QuantumReedMuller) = stab_to_gf2(parity_checks(QuantumReedMuller(c.m)))[end-(code_n(c::QuantumReedMuller)-2-c.m):end, end÷2+1:end] diff --git a/src/ecc/decoder_pipeline.jl b/src/ecc/decoder_pipeline.jl index 81fd1af30..aeb481851 100644 --- a/src/ecc/decoder_pipeline.jl +++ b/src/ecc/decoder_pipeline.jl @@ -101,7 +101,7 @@ function physical_ECC_circuit(H, setup::NaiveSyndromeECCSetup) end mem_error_circ = [PauliError(i, setup.mem_noise) for i in 1:nqubits(H)] - circ = [mem_error_circ..., noisy_syndrome_circ...] + circ = vcat(mem_error_circ, noisy_syndrome_circ) return circ, syndrome_bits, n_anc end @@ -119,7 +119,7 @@ function physical_ECC_circuit(H, setup::ShorSyndromeECCSetup) mem_error_circ = [PauliError(i, setup.mem_noise) for i in 1:nqubits(H)] - circ = [prep_anc..., mem_error_circ..., noisy_syndrome_circ...] + circ = vcat(prep_anc, mem_error_circ, noisy_syndrome_circ) circ, syndrome_bits, n_anc end @@ -141,12 +141,12 @@ function evaluate_decoder(d::AbstractSyndromeDecoder, setup::AbstractECCSetup, n # Evaluate the probability for X logical error (the Z-observable part of the faults matrix is used) X_error = evaluate_decoder( d, nsamples, - [encoding_circ..., physical_noisy_circ..., logZ_circ...], + vcat(encoding_circ, physical_noisy_circ, logZ_circ), syndrome_bits, logZ_bits, O[end÷2+1:end,:]) # Evaluate the probability for Z logical error (the X-observable part of the faults matrix is used) Z_error = evaluate_decoder( d, nsamples, - [preX..., encoding_circ..., physical_noisy_circ..., logX_circ...], + vcat(preX, encoding_circ, physical_noisy_circ, logX_circ), syndrome_bits, logX_bits, O[1:end÷2,:]) return (X_error, Z_error) end @@ -171,12 +171,21 @@ end function evaluate_guesses(measured_faults, guesses, faults_matrix) nsamples = size(guesses, 1) - guess_faults = (faults_matrix * guesses') .% 2 # TODO this can be faster and non-allocating by turning it into a loop - decoded = 0 - for i in 1:nsamples # TODO this can be faster and non-allocating by having the loop and the matrix multiplication on the line above work together and not store anything - (@view guess_faults[:,i]) == (@view measured_faults[i,:]) && (decoded += 1) + fails = 0 + for i in 1:nsamples + for j in 1:size(faults_matrix, 1) + sum_mod = 0 + @inbounds @simd for k in 1:size(faults_matrix, 2) + sum_mod += faults_matrix[j, k] * guesses[i, k] + end + sum_mod %= 2 + if sum_mod != measured_faults[i, j] + fails += 1 + break + end + end end - return (nsamples - decoded) / nsamples + return fails / nsamples end function evaluate_decoder(d::AbstractSyndromeDecoder, setup::CommutationCheckECCSetup, nsamples::Int) @@ -242,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..208f1ddfb 100644 --- a/src/grouptableaux.jl +++ b/src/grouptableaux.jl @@ -1,6 +1,3 @@ -using Graphs -using LinearAlgebra - """ Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)). @@ -15,8 +12,8 @@ julia> groupify(S"XZ ZX") ``` """ 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 + # 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)) @@ -27,7 +24,7 @@ function groupify(s::Stabilizer) end end end - return group + return group end @@ -44,8 +41,8 @@ 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 @@ -60,7 +57,7 @@ function minimal_generating_set(s::Stabilizer) end """ -Return the full Pauli group of a given length. Phases are ignored by default, +Return the full Pauli group of a given length. Phases are ignored by default, but can be included by setting `phases=true`. ```jldoctest @@ -131,8 +128,8 @@ julia> normalizer(T"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 + # 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) @@ -161,7 +158,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 +172,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,7 +180,7 @@ function centralizer(t::Tableau) end """ -Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to +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. ```jldoctest @@ -196,9 +193,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 +205,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 16b4eaf52..ad252e95f 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -161,8 +161,10 @@ julia> tensor(s, s) See also [`tensor_pow`](@ref).""" function tensor end -function tensor(ops::AbstractStabilizer...) # TODO optimize this by doing conversion to common type to enable preallocation - foldl(⊗, ops[2:end], init=ops[1]) +function tensor(ops::AbstractStabilizer...) # TODO optimize by pre-allocating one large tableau instead of the current quadratic fold + ct = promote_type(map(typeof, ops)...) + conv_ops = map(x -> convert(ct, x), ops) + return foldl(⊗, conv_ops) end """Repeated tensor product of an operators or a tableau. @@ -258,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..9803ed2b5 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) diff --git a/src/pauli_operator.jl b/src/pauli_operator.jl index e556628bc..a843ef10b 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,7 +87,7 @@ 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 +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)::Bool, ((p.xz[end÷2+_div(Tᵥₑ,i-1)+1] & Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))!=0x0)::Bool 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ᵥₑ}} @@ -176,7 +176,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..4b2e6e2e5 100644 --- a/src/project_trace_reset.jl +++ b/src/project_trace_reset.jl @@ -339,7 +339,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 +377,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 +497,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 +647,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] @@ -886,4 +886,4 @@ See also: [`traceout!`](@ref) """ function delete_columns(𝒮::Stabilizer, subset) 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 8f5857330..98f0e5d61 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) @@ -333,8 +333,8 @@ end @qubitop2 YCY (x1⊻z2⊻x2, z1⊻x2⊻z2 , x1⊻x2⊻z1 , x1⊻z1⊻z2, ~iszero( (x1 & ~z1 & ~x2 & z2) | (~x1 & z1 & x2 & ~z2))) @qubitop2 YCZ (x1⊻x2 , x2⊻z1 , x2 , z2⊻x1⊻z1, ~iszero( (x2 & (x1 ⊻ z1) & (z2 ⊻ x1)) )) -@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 & ~z2) | (x1 & ~z1 & ~x2 & z2) | (x1 & z1 & ~x2 & z2) | (x1 & z1 & x2 & z2))) +@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))) diff --git a/src/tableau_show.jl b/src/tableau_show.jl index 52e1adc07..303549573 100644 --- a/src/tableau_show.jl +++ b/src/tableau_show.jl @@ -76,15 +76,15 @@ function Base.show(io::IO, d::Destabilizer) print(io, "Destablizer $r×$q") elseif get(io, :limit, false) h,w = displaysize(io) - println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(d.tab,2)-4),0)) - _show(io, destabilizerview(d).tab, w, h÷2) - println(io, "\n𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(d.tab,2)-2),0)) - _show(io, stabilizerview(d).tab, w, h÷2) + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(tab(d),2)-4),0)) + _show(io, tab(destabilizerview(d)), w, h÷2) + println(io, "\n𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(tab(d),2)-2),0)) + _show(io, tab(stabilizerview(d)), w, h÷2) else - println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(d.tab,2)-4,0)) - _show(io, destabilizerview(d).tab, missing, missing) - println(io, "\n𝒮𝓉𝒶𝒷" * "━"^max(size(d.tab,2)-2,0)) - _show(io, stabilizerview(d).tab, missing, missing) + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(tab(d),2)-4,0)) + _show(io, tab(destabilizerview(d)), missing, missing) + println(io, "\n𝒮𝓉𝒶𝒷" * "━"^max(size(tab(d),2)-2,0)) + _show(io, tab(stabilizerview(d)), missing, missing) end end @@ -95,36 +95,36 @@ function Base.show(io::IO, d::MixedDestabilizer) print(io, "MixedDestablizer $r×$q") elseif get(io, :limit, false) h,w = displaysize(io) - println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(d.tab,2)-4),0)) - _show(io, destabilizerview(d).tab, w, h÷4) + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(tab(d),2)-4),0)) + _show(io, tab(destabilizerview(d)), w, h÷4) if r != q println(io) - println(io, "𝒳ₗ" * "━"^max(min(w-5,size(d.tab,2)),0)) - _show(io, logicalxview(d).tab, w, h÷4) + println(io, "𝒳ₗ" * "━"^max(min(w-5,size(tab(d),2)),0)) + _show(io, tab(logicalxview(d)), w, h÷4) end println(io) - println(io, "𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(d.tab,2)-2),0)) - _show(io, stabilizerview(d).tab, w, h÷4) + println(io, "𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(tab(d),2)-2),0)) + _show(io, tab(stabilizerview(d)), w, h÷4) if r != q println(io) - println(io, "𝒵ₗ" * "━"^max(min(w-5,size(d.tab,2)),0)) - _show(io, logicalzview(d).tab, w, h÷4) + println(io, "𝒵ₗ" * "━"^max(min(w-5,size(tab(d),2)),0)) + _show(io, tab(logicalzview(d)), w, h÷4) end else - println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(d.tab,2)-4,0)) - _show(io, destabilizerview(d).tab, missing, missing) + println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(tab(d),2)-4,0)) + _show(io, tab(destabilizerview(d)), missing, missing) if r != q println(io) - println(io, "𝒳ₗ" * "━"^max(size(d.tab,2),0)) - _show(io, logicalxview(d).tab, missing, missing) + println(io, "𝒳ₗ" * "━"^max(size(tab(d),2),0)) + _show(io, tab(logicalxview(d)), missing, missing) end println(io) - println(io, "𝒮𝓉𝒶𝒷" * "━"^max(size(d.tab,2)-2,0)) - _show(io, stabilizerview(d).tab, missing, missing) + println(io, "𝒮𝓉𝒶𝒷" * "━"^max(size(tab(d),2)-2,0)) + _show(io, tab(stabilizerview(d)), missing, missing) if r != q println(io) - println(io, "𝒵ₗ" * "━"^max(size(d.tab,2)),0) - _show(io, logicalzview(d).tab, missing, missing) + println(io, "𝒵ₗ" * "━"^max(size(tab(d),2)),0) + _show(io, tab(logicalzview(d)), missing, missing) end end end @@ -144,7 +144,7 @@ function _show(io::IO, c::CliffordOperator, limit=50, limit_vertical=20) continue end print(io, "X"*digits_substr(i,nwidth)*" ⟼ ") - _show(io, c.tab[i], _limit) + _show(io, tab(c)[i], _limit) println(io) end for i in range @@ -153,7 +153,7 @@ function _show(io::IO, c::CliffordOperator, limit=50, limit_vertical=20) continue end print(io, "Z"*digits_substr(i,nwidth)*" ⟼ ") - _show(io, c.tab[i+n], _limit) + _show(io, tab(c)[i+n], _limit) i!=n && println(io) end end diff --git a/test/Project.toml b/test/Project.toml index 67977b460..eb254c893 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -7,6 +7,7 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" HostCPUFeatures = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0" ILog2 = "2cd5bd5f-40a1-5050-9e10-fc8cdb6109f5" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -24,6 +25,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SIMD = "fdea26ae-647d-5447-a871-4b548cad5224" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" diff --git a/test/test_ecc.jl b/test/test_ecc.jl index 64531986d..075d27ac5 100644 --- a/test/test_ecc.jl +++ b/test/test_ecc.jl @@ -9,7 +9,7 @@ function test_naive_syndrome(c::AbstractECC, e::Bool) # create a random logical state unencoded_qubits = random_stabilizer(code_k(c)) - bufferqubits = one(Stabilizer,code_s(c)) + bufferqubits = one(Stabilizer, code_n(c) - code_k(c)) logicalqubits = bufferqubits⊗unencoded_qubits mctrajectory!(logicalqubits, naive_encoding_circuit(c)) if e @@ -50,7 +50,7 @@ ancqubits = code_s(code) regbits = ancqubits frames = PauliFrame(nframes, dataqubits+ancqubits, regbits) - circuit = [ecirc..., scirc...] + circuit = vcat(ecirc, scirc) pftrajectories(frames, circuit) @test sum(pfmeasurements(frames)) == 0 end diff --git a/test/test_ecc_base.jl b/test/test_ecc_base.jl index 9b4dac227..9f4a5ec9e 100644 --- a/test/test_ecc_base.jl +++ b/test/test_ecc_base.jl @@ -3,6 +3,10 @@ using QuantumClifford using QuantumClifford.ECC using InteractiveUtils +import Nemo: GF +import LinearAlgebra +import Hecke: group_algebra, abelian_group, gens + # generate instances of all implemented codes to make sure nothing skips being checked # We do not include smaller random circuit code because some of them has a bad distance and fails the TableDecoder test @@ -14,19 +18,59 @@ random_circuit_code_args = vcat( [map(f -> getfield(random_all_to_all_circuit_code(c...), f), fieldnames(CircuitCode)) for c in random_all_to_all_circuit_args] ) +# test codes LP04 and LP118 from Appendix A of [raveendran2022finite](@cite), +B04 = Dict( + 7 => [0 0 0 0; 0 1 2 5; 0 6 3 1], + 9 => [0 0 0 0; 0 1 6 7; 0 4 5 2], + 17 => [0 0 0 0; 0 1 2 11; 0 8 12 13], + 19 => [0 0 0 0; 0 2 6 9; 0 16 7 11] +) + +B118 = Dict( + 16 => [0 0 0 0 0; 0 2 4 7 11; 0 3 10 14 15], + 21 => [0 0 0 0 0; 0 4 5 7 17; 0 14 18 12 11], + 30 => [0 0 0 0 0; 0 2 14 24 25; 0 16 11 14 13], +) + +LP04 = [LPCode(base_matrix, l .- base_matrix', l) for (l, base_matrix) in B04] +LP118 = [LPCode(base_matrix, l .- base_matrix', l) for (l, base_matrix) in B118] + +# generalized bicyle codes from (A1) and (A2) Appendix B of [panteleev2021degenerate](@cite). +test_gb_codes = [ + generalized_bicycle_codes([0, 15, 20, 28, 66], [0, 58, 59, 100, 121], 127), # (A1) [[254, 28, 14≤d≤20]] + generalized_bicycle_codes([0, 1, 14, 16, 22], [0, 3, 13, 20, 42], 63), # (A2) [[126, 28, 8]] +] + +other_lifted_product_codes = [] + +# [[882, 24, d≤24]] code from (B1) in Appendix B of [panteleev2021degenerate](@cite) +l = 63 +GA = group_algebra(GF(2), abelian_group(l)) +A = zeros(GA, 7, 7) +x = gens(GA)[] +A[LinearAlgebra.diagind(A)] .= x^27 +A[LinearAlgebra.diagind(A, -1)] .= x^54 +A[LinearAlgebra.diagind(A, 6)] .= x^54 +A[LinearAlgebra.diagind(A, -2)] .= GA(1) +A[LinearAlgebra.diagind(A, 5)] .= GA(1) +B = reshape([1 + x + x^6], (1, 1)) +push!(other_lifted_product_codes, LPCode(A, B)) + const code_instance_args = Dict( - Toric => [(3,3), (4,4), (3,6), (4,3), (5,5)], - Surface => [(3,3), (4,4), (3,6), (4,3), (5,5)], - Gottesman => [3, 4, 5], - CSS => (c -> (parity_checks_x(c), parity_checks_z(c))).([Shor9(), Steane7(), Toric(4, 4)]), - Concat => [(Perfect5(), Perfect5()), (Perfect5(), Steane7()), (Steane7(), Cleve8()), (Toric(2, 2), Shor9())], - CircuitCode => random_circuit_code_args + :Toric => [(3,3), (4,4), (3,6), (4,3), (5,5)], + :Surface => [(3,3), (4,4), (3,6), (4,3), (5,5)], + :Gottesman => [3, 4, 5], + :CSS => (c -> (parity_checks_x(c), parity_checks_z(c))).([Shor9(), Steane7(), Toric(4, 4)]), + :Concat => [(Perfect5(), Perfect5()), (Perfect5(), Steane7()), (Steane7(), Cleve8()), (Toric(2, 2), Shor9())], + :CircuitCode => random_circuit_code_args, + :LPCode => (c -> (c.A, c.B)).(vcat(LP04, LP118, test_gb_codes, other_lifted_product_codes)), + :QuantumReedMuller => [3, 4, 5] ) function all_testablable_code_instances(;maxn=nothing) codeinstances = [] for t in subtypes(QuantumClifford.ECC.AbstractECC) - for c in get(code_instance_args, t, []) + for c in get(code_instance_args, t.name.name, []) codeinstance = t(c...) !isnothing(maxn) && nqubits(codeinstance) > maxn && continue push!(codeinstances, codeinstance) diff --git a/test/test_ecc_codeproperties.jl b/test/test_ecc_codeproperties.jl index 404592763..ec529d3e4 100644 --- a/test/test_ecc_codeproperties.jl +++ b/test/test_ecc_codeproperties.jl @@ -30,11 +30,10 @@ H = parity_checks(code) @test nqubits(code) == size(H, 2) == code_n(code) @test size(H, 1) == code_s(code) - @test code_s(code) + code_k(code) == code_n(code) - @test size(H, 1) < size(H, 2) + @test code_s(code) + code_k(code) >= code_n(code) # possibly exist redundant checks _, _, rank = canonicalize!(copy(H), ranks=true) - @test rank == size(H, 1) # TODO maybe weaken this if we want to permit codes with redundancies - @test QuantumClifford.stab_looks_good(copy(H)) + @test rank <= size(H, 1) + @test QuantumClifford.stab_looks_good(copy(H), remove_redundant_rows=true) end end end diff --git a/test/test_ecc_decoder_all_setups.jl b/test/test_ecc_decoder_all_setups.jl index f6ceee981..ecd747da3 100644 --- a/test/test_ecc_decoder_all_setups.jl +++ b/test/test_ecc_decoder_all_setups.jl @@ -26,7 +26,7 @@ #@show c #@show s #@show e - @assert max(e...) < noise/4 + @test max(e...) < noise/4 end end end @@ -35,8 +35,39 @@ ## @testset "belief prop decoders, good for sparse codes" begin + codes = vcat(LP04, LP118, test_gb_codes, other_lifted_product_codes) + + noise = 0.001 + + setups = [ + CommutationCheckECCSetup(noise), + NaiveSyndromeECCSetup(noise, 0), + ShorSyndromeECCSetup(noise, 0), + ] + + for c in codes + for s in setups + for d in [c -> PyBeliefPropOSDecoder(c, maxiter=2)] + nsamples = 10000 + if true + @test_broken false # TODO these are too slow to test in CI + continue + end + e = evaluate_decoder(d(c), s, nsamples) + # @show c + # @show s + # @show e + @test max(e...) <= noise + end + end + end + end + + + @testset "BitFlipDecoder decoder, good for sparse codes" begin codes = [ - # TODO + QuantumReedMuller(3), + QuantumReedMuller(4) ] noise = 0.001 @@ -49,12 +80,12 @@ for c in codes for s in setups - for d in [c->PyBeliefPropOSDecoder(c, maxiter=10)] + for d in [c->BitFlipDecoder(c, maxiter=10)] e = evaluate_decoder(d(c), s, 100000) - @show c - @show s - @show e - @assert max(e...) < noise/4 + #@show c + #@show s + #@show e + @test max(e...) < noise/4 end end end @@ -91,7 +122,7 @@ #@show c #@show s #@show e - @assert max(e...) < noise/5 + @test max(e...) < noise/5 end end end diff --git a/test/test_ecc_encoding.jl b/test/test_ecc_encoding.jl index cab0562ff..78312c57d 100644 --- a/test/test_ecc_encoding.jl +++ b/test/test_ecc_encoding.jl @@ -31,7 +31,7 @@ pre_encₖ = one(Stabilizer, code_k(code)) # n-k ancillary qubits in state zero prepended - pre_encₙ = one(Stabilizer, code_s(code)) ⊗ pre_encₖ + pre_encₙ = one(Stabilizer, code_n(code) - code_k(code)) ⊗ pre_encₖ # running the encoding circuit encodedₙ = mctrajectory!(pre_encₙ, circ)[1] |> canonicalize! diff --git a/test/test_ecc_quantumreedmuller.jl b/test/test_ecc_quantumreedmuller.jl new file mode 100644 index 000000000..efb76918f --- /dev/null +++ b/test/test_ecc_quantumreedmuller.jl @@ -0,0 +1,50 @@ +@testitem "Quantum Reed-Muller" begin + using Test + using Nemo: echelon_form, matrix, GF + using LinearAlgebra + using QuantumClifford + using QuantumClifford: canonicalize!, Stabilizer, stab_to_gf2 + using QuantumClifford.ECC + using QuantumClifford.ECC: AbstractECC, QuantumReedMuller, Steane7, CSS + + function designed_distance(mat) + dist = 3 + for row in eachrow(mat) + count = sum(row) + if count < dist + return false + end + end + return true + end + + @testset "Test QuantumReedMuller(r,m) properties" begin + for m in 3:10 + stab = parity_checks(QuantumReedMuller(m)) + H = stab_to_gf2(stab) + @test designed_distance(H) == true + # QuantumReedMuller(3) is the Steane7 code. + @test canonicalize!(parity_checks(Steane7())) == parity_checks(QuantumReedMuller(3)) + @test code_n(QuantumReedMuller(m)) == 2^m - 1 + @test code_k(QuantumReedMuller(m)) == 1 + @test distance(QuantumReedMuller(m)) == 3 + @test H == stab_to_gf2(parity_checks(CSS(parity_checks_x(QuantumReedMuller(m)), parity_checks_z(QuantumReedMuller(m))))) + # [[15,1,3]] qrm code from table 1 of https://arxiv.org/pdf/1705.0010 + qrm₁₅₁₃ = S"ZIZIZIZIZIZIZIZ + IZZIIZZIIZZIIZZ + IIIZZZZIIIIZZZZ + IIIIIIIZZZZZZZZ + IIZIIIZIIIZIIIZ + IIIIZIZIIIIIZIZ + IIIIIZZIIIIIIZZ + IIIIIIIIIZZIIZZ + IIIIIIIIIIIZZZZ + IIIIIIIIZIZIZIZ + XIXIXIXIXIXIXIX + IXXIIXXIIXXIIXX + IIIXXXXIIIIXXXX + IIIIIIIXXXXXXXX" + @test canonicalize!(parity_checks(qrm₁₅₁₃)) == canonicalize!(parity_checks(QuantumReedMuller(4))) + end + end +end diff --git a/test/test_ecc_syndromes.jl b/test/test_ecc_syndromes.jl index 45f4e0f89..0f6f67ea0 100644 --- a/test/test_ecc_syndromes.jl +++ b/test/test_ecc_syndromes.jl @@ -17,8 +17,8 @@ # no noise naive_frames = PauliFrame(nframes, naive_qubits, syndromebits) shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits)) - naive_circuit = [ecirc..., naive_scirc...] - shor_circuit = [ecirc..., shor_cat_scirc..., shor_scirc...] + naive_circuit = vcat(ecirc, naive_scirc) + shor_circuit = vcat(ecirc, shor_cat_scirc, shor_scirc) pftrajectories(naive_frames, naive_circuit) pftrajectories(shor_frames, shor_circuit) @test pfmeasurements(naive_frames) == pfmeasurements(shor_frames)[:,shor_bits] @@ -27,7 +27,7 @@ naive_frames = PauliFrame(nframes, naive_qubits, syndromebits) shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits)) pftrajectories(naive_frames, ecirc) - pftrajectories(shor_frames, [ecirc..., shor_cat_scirc...]) + pftrajectories(shor_frames, vcat(ecirc, shor_cat_scirc)) # manually injecting the same type of noise in the frames -- not really a user accessible API p = random_pauli(dataqubits, realphase=true) pₙ = embed(naive_qubits, 1:dataqubits, p) diff --git a/test/test_group_tableaux.jl b/test/test_group_tableaux.jl index afec47e90..d208afdff 100644 --- a/test/test_group_tableaux.jl +++ b/test/test_group_tableaux.jl @@ -1,11 +1,10 @@ -@testitem "Classical" begin +@testitem "group theory routines" begin using Test - using Random using QuantumClifford # Including sizes that would test off-by-one errors in the bit encoding. - test_sizes = [1, 2, 3, 4, 5, 7, 8, 9, 15, 16, 17] + test_sizes = [1, 2, 3, 4, 5, 7, 8, 9, 15, 16, 17] # Zero function(in groupify) slows down around 2^30(n=30),eventually breaks small_test_sizes = [1, 2, 3, 4, 5, 7] # Pauligroup slows around n = 8 @@ -45,7 +44,7 @@ end end #Test pauligroup - for n in [1, small_test_sizes...] + for n in [1, small_test_sizes...] @test length(QuantumClifford.pauligroup(n, phases=false)) == 4^n @test length(QuantumClifford.pauligroup(n, phases=true)) == 4^(n+1) end @@ -98,7 +97,7 @@ end c = contractor(s, subset) count = 0 - for stabilizer in s + for stabilizer in s contractable = true for i in subset if stabilizer[i] != (false, false) contractable = false end @@ -111,9 +110,9 @@ p = zero(PauliOperator, nqubits(s)) index = 0 for i in 1:nqubits(s) - if !(i in subset) - index+=1 - p[i] = contracted[index] + if !(i in subset) + index+=1 + p[i] = contracted[index] end end @test p in s || -1* p in s || 1im * p in s || -1im * p in s @@ -123,4 +122,4 @@ end end end -end \ No newline at end of file +end diff --git a/test/test_jet.jl b/test/test_jet.jl index 304f3b78c..7f0248ab7 100644 --- a/test/test_jet.jl +++ b/test/test_jet.jl @@ -8,6 +8,8 @@ using LinearAlgebra using Nemo using AbstractAlgebra + using Hecke + using StaticArrays rep = report_package("QuantumClifford"; ignored_modules=( @@ -19,9 +21,11 @@ AnyFrameModule(LinearAlgebra), AnyFrameModule(Nemo), AnyFrameModule(AbstractAlgebra), + AnyFrameModule(Hecke), + AnyFrameModule(StaticArrays), )) @show rep - @test_broken length(JET.get_reports(rep)) == 0 - @test length(JET.get_reports(rep)) <= 11 + @show length(JET.get_reports(rep)) + @test length(JET.get_reports(rep)) == 0 end diff --git a/test/test_noisycircuits.jl b/test/test_noisycircuits.jl index 8ae674161..4f5e92c3e 100644 --- a/test/test_noisycircuits.jl +++ b/test/test_noisycircuits.jl @@ -228,7 +228,7 @@ state = Register(MixedDestabilizer(S"ZZ"), zeros(Bool, 1)) meas = PauliMeasurement(P"ZI", 1) state, flag = applywstatus!(state, meas) - @test state.stab.rank == 2 + @test rank(state.stab) == 2 tab(state.stab).phases .= 0 @test stabilizerview(state.stab) == S"ZZ ZI" diff --git a/test/test_projections.jl b/test/test_projections.jl index 03d56e1e2..beae0ab70 100644 --- a/test/test_projections.jl +++ b/test/test_projections.jl @@ -76,12 +76,12 @@ @test_throws BadDataStructure pds, a, r = project!(copy(ds),p) pms, a, r = project!(copy(ms),p) @test mixed_stab_looks_good(pms) - @test pms.rank==3 - @test a==pms.rank && isnothing(r) + @test rank(pms)==3 + @test a==rank(pms) && isnothing(r) pmds, a, r = project!(copy(mds),p) @test mixed_destab_looks_good(pmds) - @test pmds.rank==3 - @test a==pmds.rank && isnothing(r) + @test rank(pmds)==3 + @test a==rank(pmds) && isnothing(r) p = P"ZZI" ps, a, r = project!(copy(s),p) @@ -90,11 +90,11 @@ @test_throws BadDataStructure pds, a, r = project!(copy(ds),p) pms, a, r = project!(copy(ms),p) @test mixed_stab_looks_good(pms) - @test pms.rank==2 + @test rank(pms)==2 @test a==0 && r==0x2 pmds, a, r = project!(copy(mds),p) @test mixed_destab_looks_good(pmds) - @test pmds.rank==2 + @test rank(pmds)==2 @test a==0 && r==0x2 @test canonicalize!(ps)==canonicalize!(stabilizerview(pms))==canonicalize!(stabilizerview(pmds)) @@ -107,11 +107,11 @@ @test a==2 && isnothing(r) pms, a, r = project!(copy(ms),p) @test mixed_stab_looks_good(pms) - @test pms.rank==2 + @test rank(pms)==2 @test a==2 && isnothing(r) pmds, a, r = project!(copy(mds),p) @test mixed_destab_looks_good(pmds) - @test pmds.rank==2 + @test rank(pmds)==2 @test a==2 && isnothing(r) @test canonicalize!(ps)==canonicalize!(stabilizerview(pms))==canonicalize!(stabilizerview(pds))==canonicalize!(stabilizerview(pmds)) end @@ -162,10 +162,10 @@ s = MixedStabilizer(s, 2) ms, a, r = project!(copy(s), P"IZI") @test (a, r) == (0, 0x0) # on commuting operator in the stabilizer - @test ms.rank == 2 + @test rank(ms) == 2 ms, a, r = project!(copy(s), P"IIZ") @test (a, r) == (3, nothing) # on commuting operator out of the stabilizer - @test ms.rank == 3 + @test rank(ms) == 3 s = S"ZII IZI" s = Destabilizer(s) @test_throws BadDataStructure project!(copy(s), P"IZI"; keep_result=true) # on comm @@ -182,16 +182,16 @@ s = MixedDestabilizer(s) mds, a, r = project!(copy(s), P"IZI"; keep_result=true) @test (a, r) == (0, 0x0) # on commuting operator in the stabilizer - @test mds.rank == 2 + @test rank(mds) == 2 mds, a, r = project!(copy(s), P"IIZ"; keep_result=true) @test (a, r) == (3, nothing) # on commuting operator out of the stabilizer - @test mds.rank == 3 + @test rank(mds) == 3 mds, a, r = project!(copy(s), P"IZI"; keep_result=false) @test (a, r) == (0, nothing) # on commuting operator in the stabilizer - @test mds.rank == 2 + @test rank(mds) == 2 mds, a, r = project!(copy(s), P"IIZ"; keep_result=false) @test (a, r) == (3, nothing) # on commuting operator out of the stabilizer - @test mds.rank == 3 + @test rank(mds) == 3 end @testset "Results from canonicalization vs from destabilizer" begin @test generate!(P"_Z", S"XZ") === nothing # for bug fixed in 4b536231c3ee4e6446262fcc61ba8da669415bc8 @@ -206,7 +206,7 @@ _, ams, rms = project!(ms,p) _, amd, rmd = project!(md,p) @test rs == rms == rmd - @test (md.rank!=r) || (canonicalize!(s) == canonicalize!(stabilizerview(ms))) + @test (rank(md)!=r) || (canonicalize!(s) == canonicalize!(stabilizerview(ms))) @test canonicalize!(stabilizerview(ms)) == canonicalize!(stabilizerview(md)) if as == 0 @test ams == amd == 0 diff --git a/test/test_stabs.jl b/test/test_stabs.jl index 63a5c562e..b479934c2 100644 --- a/test/test_stabs.jl +++ b/test/test_stabs.jl @@ -55,6 +55,12 @@ stabs = [s[1:i] for s in [random_stabilizer(n) for n in [32,16,16,64,63,65,129,128,127]] for i in rand(1:10)]; mdstabs = MixedDestabilizer.(stabs); @test canonicalize!(⊗(stabs...)) == canonicalize!(stabilizerview(⊗(mdstabs...))) + md = MixedDestabilizer(random_destabilizer(n)) + s = random_stabilizer(n) + mds = md⊗s + @test mixed_destab_looks_good(mds) + estab = stabilizerview(md)⊗s + @test canonicalize!(copy(stabilizerview(mds))) == canonicalize!(estab) end end @@ -90,4 +96,15 @@ @test stab_to_gf2(s2a) == stab_to_gf2(s2b) end end + + @testset "horizontal concatenation" begin + @test hcat(ghz(2), ghz(2)) == S"XXXX ZZZZ" + s1 = S"YZ -XX" + s2 = S"-ZY -YX" + @test hcat(copy(s1), copy(s2)) == S"-YZZY XXYX" + @test hcat(copy(s1), copy(s2), copy(s1), copy(s2)) == S"YZZYYZZY XXYXXXYX" + @test_throws ArgumentError hcat(copy(s1), random_stabilizer(3)) + @test hcat(copy(tab(s1)), copy(tab(s2))) == T"-YZZY XXYX" + @test hcat(copy(tab(s1)), copy(tab(s2)), copy(tab(s1)), copy(tab(s2))) == T"YZZYYZZY XXYXXXYX" + end end diff --git a/test/test_symcontrolled.jl b/test/test_symcontrolled.jl index 256f62be3..4fbd86a87 100644 --- a/test/test_symcontrolled.jl +++ b/test/test_symcontrolled.jl @@ -125,4 +125,21 @@ end end end + + @testset "Ket-based definition for Y and rY" begin + for control in (:Z,) + target = :Y + for r in (true, false) + s = Stabilizer(QuantumClifford._T_str(string(control))) + k1 = Ket(s) + s.tab.phases[1] = 0x2 + k2 = Ket(s) + i = Operator(tId1) + o = Operator(CliffordOperator(eval(Symbol(:s,target,))(1),1)) + gate = projector(k1)⊗i + (r ? -1 : -im) * projector(k2)⊗o # XXX there is a -1 here because global phase is arbitrary and thus our Operator(tY) is NOT im*Operator(X)*Operator(Z) + implemented_gate = Operator(CliffordOperator(eval(Symbol(:s,control,:C,(r ? :rY : :Y)))(1,2),2)) + @test gate≈implemented_gate + end + end + end end diff --git a/test/test_trace.jl b/test/test_trace.jl index 424399ce1..6ed1022e2 100644 --- a/test/test_trace.jl +++ b/test/test_trace.jl @@ -101,7 +101,7 @@ @test canonicalize!(copy(ssr2v)[:,perm])[1:z] == canonicalize!(copy(newstate)) @test canonicalize!(msr2v) == c[1:z] # Compare different datastractures - @test canonicalize!(copy(stabilizerview(mdr1)))==canonicalize!(copy(stabilizerview(msr1)))==canonicalize!(ssr1[1:mdr1.rank]) + @test canonicalize!(copy(stabilizerview(mdr1)))==canonicalize!(copy(stabilizerview(msr1)))==canonicalize!(ssr1[1:rank(mdr1)]) end end end