Skip to content

Commit

Permalink
Merge branch 'master' into BBQLDPC
Browse files Browse the repository at this point in the history
  • Loading branch information
Fe-r-oz committed Sep 24, 2024
2 parents 7eeb82c + 05491d9 commit 8c4d158
Show file tree
Hide file tree
Showing 14 changed files with 396 additions and 172 deletions.
6 changes: 6 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,12 @@

# News

## dev

- 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).
Expand Down
33 changes: 26 additions & 7 deletions src/QuantumClifford.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ export
# Symbolic Clifford Ops
AbstractSymbolicOperator, AbstractSingleQubitOperator, AbstractTwoQubitOperator,
sHadamard, sPhase, sInvPhase, SingleQubitOperator, sId1, sX, sY, sZ,
sHadamardXY, sHadamardYZ, sSQRTX, sInvSQRTX, sSQRTY, sInvSQRTY, sCXYZ, sCZYX,
sCNOT, sCPHASE, sSWAP,
sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZCX, sZCY, sZCZ,
sZCrY, sInvZCrY,
Expand Down Expand Up @@ -149,11 +150,20 @@ function Tableau(paulis::AbstractVector{PauliOperator{Tₚ,Tᵥ}}) where {Tₚ<:
tab
end

Tableau(phases::AbstractVector{UInt8}, xs::AbstractMatrix{Bool}, zs::AbstractMatrix{Bool}) = Tableau(
phases, size(xs,2),
vcat(hcat((BitArray(xs[i,:]).chunks for i in 1:size(xs,1))...)::Matrix{UInt},
hcat((BitArray(zs[i,:]).chunks for i in 1:size(zs,1))...)::Matrix{UInt}) # type assertions to help Julia infer types
)
function Tableau(phases::AbstractVector{UInt8}, xs::AbstractMatrix{Bool}, zs::AbstractMatrix{Bool})
r_xs = size(xs, 1)
r_zs = size(zs, 1)
if length(phases) != r_xs || r_xs != r_zs
throw(DimensionMismatch(lazy"The length of phases ($(length(phases))), rows of xs ($r_xs), rows of zs ($r_zs) must all be equal."))
end
Tableau(
phases,size(xs, 2),
vcat(
hcat((BitArray(xs[i, :]).chunks for i in 1:r_xs)...)::Matrix{UInt},
hcat((BitArray(zs[i, :]).chunks for i in 1:r_zs)...)::Matrix{UInt} # type assertions to help Julia infer types
)
)
end

Tableau(phases::AbstractVector{UInt8}, xzs::AbstractMatrix{Bool}) = Tableau(phases, xzs[:,1:end÷2], xzs[:,end÷2+1:end])

Expand Down Expand Up @@ -758,10 +768,19 @@ function comm!(v, l::PauliOperator, r::Tableau)
v
end
comm!(v, l::Tableau, r::PauliOperator) = comm!(v, r, l)
@inline comm!(v, l::PauliOperator, r::Stabilizer, i::Int) = comm!(v, l, tab(r), i)
@inline comm!(v, l::Stabilizer, r::PauliOperator, i::Int) = comm!(v, tab(l), r, i)
@inline comm!(v, l::PauliOperator, r::Stabilizer) = comm!(v, l, tab(r))
@inline comm!(v, l::Stabilizer, r::PauliOperator) = comm!(v, tab(l), r)
function comm!(v, l::PauliOperator, r::Tableau, i)
v[i] = comm(l,r,i)
v
end
comm!(v, l::Tableau, r::PauliOperator, i) = comm!(v, r, l, i)
@inline comm!(v, l::PauliOperator, r::Stabilizer, i::Int) = comm!(v, l, tab(r), i)
@inline comm!(v, l::Stabilizer, r::PauliOperator, i::Int) = comm!(v, tab(l), r, i)
function comm!(v, s::Tableau, l::Int, r::Int)
v[l] = comm(s, l, r)
v
end
@inline comm!(v, s::Stabilizer, l::Int, r::Int) = comm!(v, tab(s), l, r)


Expand Down
2 changes: 1 addition & 1 deletion src/dense_cliffords.jl
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ function row_limit(str, limit=50)
end

digits_subchars = collect("₀₁₂₃₄₅₆₇₈₉")
digits_substr(n,nwidth) = join(([digits_subchars[d+1] for d in reverse(digits(n, pad=nwidth))]))
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))
Expand Down
12 changes: 9 additions & 3 deletions src/ecc/ECC.jl
Original file line number Diff line number Diff line change
@@ -1,14 +1,15 @@
module ECC

using LinearAlgebra
using LinearAlgebra: I
using QuantumClifford
using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
using DocStringExtensions
using Combinatorics: combinations
using SparseArrays: sparse
using Statistics: std
using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible, nullspace, hom, free_module, kernel, domain, map, gens
using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible, echelon_form, nullspace, hom, free_module, kernel, domain, map, gens

Check warning on line 12 in src/ecc/ECC.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"hom" should be "home".

abstract type AbstractECC end

Expand Down Expand Up @@ -70,6 +71,9 @@ In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the gener
"""
function generator_polynomial end

"""The generator matrix of a code."""
function generator end

parity_checks(s::Stabilizer) = s
Stabilizer(c::AbstractECC) = parity_checks(c)
MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
Expand Down Expand Up @@ -335,8 +339,9 @@ isdegenerate(c::AbstractECC, args...) = isdegenerate(parity_checks(c), args...)
isdegenerate(c::AbstractStabilizer, args...) = isdegenerate(stabilizerview(c), args...)

function isdegenerate(H::Stabilizer, errors) # Described in https://quantumcomputing.stackexchange.com/questions/27279
syndromes = comm.((H,), errors) # TODO This can be optimized by having something that always returns bitvectors
return length(Set(syndromes)) != length(errors)
syndromes = map(e -> comm(H,e), errors) # TODO This can be optimized by having something that always returns bitvectors
syndrome_set = Set(syndromes)
return length(syndrome_set) != length(errors)
end

function isdegenerate(H::Stabilizer, d::Int=1)
Expand Down Expand Up @@ -364,4 +369,5 @@ include("codes/random_circuit.jl")
include("codes/bbqldpc.jl")
include("codes/classical/reedmuller.jl")
include("codes/classical/bch.jl")
include("codes/classical/recursivereedmuller.jl")
end #module
64 changes: 64 additions & 0 deletions src/ecc/codes/classical/recursivereedmuller.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
"""
A construction of the Reed-Muller class of codes using the recursive definition.
The Plotkin `(u, u + v)` construction defines a recursive relation between generator matrices of Reed-Muller `(RM)` codes [abbe2020reed](@cite). To derive the generator matrix `G(m, r)` for `RM(r, m)`, the generator matrices of lower-order codes are utilized:
- `G(r - 1, m - 1)`: Generator matrix of `RM(r - 1, m - 1)`
- `G(r, m - 1)`: Generator matrix of `RM(r, m - 1)`
The generator matrix `G(m, r)` of `RM(m, r)` is formulated as follows in matrix notation:
```math
G(m, r) = \begin{bmatrix}
G(r, m - 1) & G(r, m - 1) \\
0 & G(r - 1, m - 1)
\end{bmatrix}
```
Here, the matrix 0 denotes an all-zero matrix with dimensions matching `G(r - 1, m - 1)`. This recursive approach facilitates the construction of higher-order Reed-Muller codes based on the generator matrices of lower-order codes.
In addition, the dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
See also: `ReedMuller`
"""
struct RecursiveReedMuller <: ClassicalCode
r::Int
m::Int

function RecursiveReedMuller(r, m)
if r < 0 || r > m
throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
end
new(r, m)
end
end

function _recursiveReedMuller(r::Int, m::Int)
if r == 1 && m == 1
return Matrix{Int}([1 1; 0 1])
elseif r == m
return Matrix{Int}(I, 2^m, 2^m)
elseif r == 0
return Matrix{Int}(ones(1, 2^m))
else
Gᵣₘ₋₁ = _recursiveReedMuller(r, m - 1)
Gᵣ₋₁ₘ₋₁ = _recursiveReedMuller(r - 1, m - 1)
return vcat(hcat(Gᵣₘ₋₁, Gᵣₘ₋₁), hcat(zeros(Int, size(Gᵣ₋₁ₘ₋₁)...), Gᵣ₋₁ₘ₋₁))
end
end

function generator(c::RecursiveReedMuller)
return _recursiveReedMuller(c.r, c.m)
end

function parity_checks(c::RecursiveReedMuller)
H = generator(RecursiveReedMuller(c.m - c.r - 1, c.m))
return H
end

code_n(c::RecursiveReedMuller) = 2 ^ c.m

code_k(c::RecursiveReedMuller) = sum(binomial.(c.m, 0:c.r))

distance(c::RecursiveReedMuller) = 2 ^ (c.m - c.r)

rate(c::RecursiveReedMuller) = code_k(c::RecursiveReedMuller) / code_n(c::RecursiveReedMuller)
52 changes: 37 additions & 15 deletions src/ecc/codes/classical/reedmuller.jl
Original file line number Diff line number Diff line change
@@ -1,38 +1,60 @@
"""The family of Reed-Muller codes, as discovered by Muller in his 1954 paper [muller1954application](@cite) and Reed who proposed the first efficient decoding algorithm [reed1954class](@cite).
Let `m` be a positive integer and `r` a nonnegative integer with `r ≤ m`. These linear codes, denoted as `RM(r, m)`, have order `r` (where `0 ≤ r ≤ m`) and codeword length `n` of `2ᵐ`.
Two special cases of generator(RM(r, m)) exist:
1. `generator(RM(0, m))`: This is the `0ᵗʰ`-order `RM` code, similar to the binary repetition code with length `2ᵐ`. It's characterized by a single basis vector containing all ones.
2. `generator(RM(m, m))`: This is the `mᵗʰ`-order `RM` code. It encompasses the entire field `F(2ᵐ)`, representing all possible binary strings of length `2ᵐ`.
You might be interested in consulting [raaphorst2003reed](@cite), [abbe2020reed](@cite), and [djordjevic2021quantum](@cite) as well.
The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller)
"""
The dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
abstract type ClassicalCode end
The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller).
See also: `RecursiveReedMuller`
"""
struct ReedMuller <: ClassicalCode
r::Int
m::Int

function ReedMuller(r, m)
if r < 0 || m < 1 || m >= 11
throw(ArgumentError("Invalid parameters: r must be non-negative and m must be positive and < 11 in order to obtain a valid code and to remain tractable"))
if r < 0 || r > m
throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
end
new(r, m)
end
end

function variables_xi(m, i)
return repeat([fill(1, 2^(m - i - 1)); fill(0, 2^(m - i - 1))], outer = 2^i)
function _variablesₓᵢ_rm(m, i)
return repeat([fill(1, 2 ^ (m - i - 1)); fill(0, 2 ^ (m - i - 1))], outer = 2 ^ i)
end

function vmult(vecs...)
function _vmult_rm(vecs...)
return [reduce(*, a, init=1) for a in zip(vecs...)]
end

function parity_checks(c::ReedMuller)
r=c.r
m=c.m
xi = [variables_xi(m, i) for i in 0:m - 1]
row_matrices = [reduce(vmult, [xi[i + 1] for i in S], init = ones(Int, 2^m)) for s in 0:r for S in combinations(0:m - 1, s)]
function generator(c::ReedMuller)
r = c.r
m = c.m
xᵢ = [_variablesₓᵢ_rm(m, i) for i in 0:m - 1]
row_matrices = [reduce(_vmult_rm, [xᵢ[i + 1] for i in S], init = ones(Int, 2 ^ m)) for s in 0:r for S in combinations(0:m - 1, s)]
rows = length(row_matrices)
cols = length(row_matrices[1])
H = reshape(vcat(row_matrices...), cols, rows)'
end
G = reshape(vcat(row_matrices...), cols, rows)'
G = Matrix{Bool}(G)
return G
end

function parity_checks(c::ReedMuller)
H = generator(ReedMuller(c.m - c.r - 1, c.m))
return H
end

code_n(c::ReedMuller) = 2 ^ c.m

code_k(c::ReedMuller) = sum(binomial.(c.m, 0:c.r))

distance(c::ReedMuller) = 2 ^ (c.m - c.r)

rate(c::ReedMuller) = code_k(c::ReedMuller) / code_n(c::ReedMuller)
2 changes: 1 addition & 1 deletion src/ecc/decoder_pipeline.jl
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ function evaluate_decoder(d::AbstractSyndromeDecoder, setup::AbstractECCSetup, n

physical_noisy_circ, syndrome_bits, n_anc = physical_ECC_circuit(H, setup)
encoding_circ = naive_encoding_circuit(H)
preX = [sHadamard(i) for i in n-k+1:n]
preX = sHadamard[sHadamard(i) for i in n-k+1:n]

mdH = MixedDestabilizer(H)
logX_circ, _, logX_bits = naive_syndrome_circuit(logicalxview(mdH), n_anc+1, last(syndrome_bits)+1)
Expand Down
10 changes: 10 additions & 0 deletions src/misc_ops.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,11 +29,21 @@ end
struct SparseGate{T<:Tableau} <: AbstractCliffordOperator # TODO simplify type parameters (remove nesting)
cliff::CliffordOperator{T}
indices::Vector{Int}
function SparseGate(cliff::CliffordOperator{T}, indices::Vector{Int}) where T<:Tableau
if length(indices) != nqubits(cliff)
throw(ArgumentError("The number of target qubits (indices) must match the qubit count in the CliffordOperator."))
end
new{T}(cliff, indices)
end
end

SparseGate(c,t::Tuple) = SparseGate(c,collect(t))

function apply!(state::AbstractStabilizer, g::SparseGate; kwargs...)
m = maximum(g.indices)
if m > nqubits(state)
throw(ArgumentError(lazy"SparseGate was attempted on invalid qubit index $(m) when the state contains only $(nqubits(state)) qubits."))
end
apply!(state, g.cliff, g.indices; kwargs...)
end

Expand Down
65 changes: 44 additions & 21 deletions src/symbolic_cliffords.jl
Original file line number Diff line number Diff line change
Expand Up @@ -85,12 +85,20 @@ macro qubitop1(name, kernel)
end
end

@qubitop1 Hadamard (z , x , x!=0 && z!=0)
@qubitop1 Phase (x , xz , x!=0 && z!=0)
@qubitop1 InvPhase (x , xz , x!=0 && z==0)
@qubitop1 X (x , z , z!=0)
@qubitop1 Y (x , z , (xz)!=0)
@qubitop1 Z (x , z , x!=0)
@qubitop1 Hadamard (z ,x , x!=0 && z!=0)
@qubitop1 HadamardXY (x ,xz , x==0 && z!=0)
@qubitop1 HadamardYZ (xz ,z , x!=0 && z==0)
@qubitop1 Phase (x ,xz , x!=0 && z!=0)
@qubitop1 InvPhase (x ,xz , x!=0 && z==0)
@qubitop1 X (x ,z , z!=0)
@qubitop1 Y (x ,z , (xz)!=0)
@qubitop1 Z (x ,z , x!=0)
@qubitop1 SQRTX (xz ,z , x==0 && z!=0)
@qubitop1 InvSQRTX (xz ,z , x!=0 && z!=0)
@qubitop1 SQRTY (z ,x , z==0)
@qubitop1 InvSQRTY (z ,x , z!=0 && x==0)
@qubitop1 CXYZ (xz ,x , z==0 && x==0)
@qubitop1 CZYX (z ,xz , z==0 && x==0)

"""A "symbolic" single-qubit Identity operation.
Expand Down Expand Up @@ -177,13 +185,21 @@ function _apply!(stab::AbstractStabilizer, op::SingleQubitOperator; phases::Val{
stab
end

SingleQubitOperator(h::sHadamard) = SingleQubitOperator(h.q, false, true , true , false, false, false)
SingleQubitOperator(p::sPhase) = SingleQubitOperator(p.q, true , true , false, true , false, false)
SingleQubitOperator(p::sInvPhase) = SingleQubitOperator(p.q, true , true , false, true , true , false)
SingleQubitOperator(p::sId1) = SingleQubitOperator(p.q, true , false, false, true , false, false)
SingleQubitOperator(p::sX) = SingleQubitOperator(p.q, true , false, false, true , false, true)
SingleQubitOperator(p::sY) = SingleQubitOperator(p.q, true , false, false, true , true , true)
SingleQubitOperator(p::sZ) = SingleQubitOperator(p.q, true , false, false, true , true , false)
SingleQubitOperator(h::sHadamard) = SingleQubitOperator(h.q, false, true , true , false, false, false)
SingleQubitOperator(p::sPhase) = SingleQubitOperator(p.q, true , true , false, true , false, false)
SingleQubitOperator(p::sInvPhase) = SingleQubitOperator(p.q, true , true , false, true , true , false)
SingleQubitOperator(p::sId1) = SingleQubitOperator(p.q, true , false, false, true , false, false)
SingleQubitOperator(p::sX) = SingleQubitOperator(p.q, true , false, false, true , false, true)
SingleQubitOperator(p::sY) = SingleQubitOperator(p.q, true , false, false, true , true , true)
SingleQubitOperator(p::sZ) = SingleQubitOperator(p.q, true , false, false, true , true , false)
SingleQubitOperator(p::sCXYZ) = SingleQubitOperator(p.q, true , true , true , false, false, false)
SingleQubitOperator(p::sCZYX) = SingleQubitOperator(p.q, false, true , true , true , false, false)
SingleQubitOperator(p::sHadamardXY) = SingleQubitOperator(p.q, true , true , false, true , false, true)
SingleQubitOperator(p::sHadamardYZ) = SingleQubitOperator(p.q, true , false, true , true , true , false)
SingleQubitOperator(p::sSQRTX) = SingleQubitOperator(p.q, true , false, true , true , false, true)
SingleQubitOperator(p::sInvSQRTX) = SingleQubitOperator(p.q, true , false, true , true , false, false)
SingleQubitOperator(p::sSQRTY) = SingleQubitOperator(p.q, false, true , true , false, true , false)
SingleQubitOperator(p::sInvSQRTY) = SingleQubitOperator(p.q, false, true , true , false, false, true)
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`."))
Expand Down Expand Up @@ -232,14 +248,21 @@ function LinearAlgebra.inv(op::SingleQubitOperator)
return SingleQubitOperator(c, op.q)
end

LinearAlgebra.inv(h::sHadamard) = sHadamard(h.q)
LinearAlgebra.inv(p::sPhase) = sInvPhase(p.q)
LinearAlgebra.inv(p::sInvPhase) = sPhase(p.q)
LinearAlgebra.inv(p::sId1) = sId1(p.q)
LinearAlgebra.inv(p::sX) = sX(p.q)
LinearAlgebra.inv(p::sY) = sY(p.q)
LinearAlgebra.inv(p::sZ) = sZ(p.q)

LinearAlgebra.inv(h::sHadamard) = sHadamard(h.q)
LinearAlgebra.inv(p::sPhase) = sInvPhase(p.q)
LinearAlgebra.inv(p::sInvPhase) = sPhase(p.q)
LinearAlgebra.inv(p::sId1) = sId1(p.q)
LinearAlgebra.inv(p::sX) = sX(p.q)
LinearAlgebra.inv(p::sY) = sY(p.q)
LinearAlgebra.inv(p::sZ) = sZ(p.q)
LinearAlgebra.inv(p::sHadamardXY) = sHadamardXY(p.q)
LinearAlgebra.inv(p::sHadamardYZ) = sHadamardYZ(p.q)
LinearAlgebra.inv(p::sSQRTX) = sInvSQRTX(p.q)
LinearAlgebra.inv(p::sInvSQRTX) = sSQRTX(p.q)
LinearAlgebra.inv(p::sSQRTY) = sInvSQRTY(p.q)
LinearAlgebra.inv(p::sInvSQRTY) = sSQRTY(p.q)
LinearAlgebra.inv(p::sCZYX) = sCXYZ(p.q)
LinearAlgebra.inv(p::sCXYZ) = sCZYX(p.q)
##############################
# Two-qubit gates
##############################
Expand Down
Loading

0 comments on commit 8c4d158

Please sign in to comment.