implement fault_matrix (#138)

* implement fault_matrix

* fix printing

CHANGELOG.md

@@ -5,6 +5,13 @@
 
 # News
 
+## v0.8.9 - 2023-07-04
+
+- In the unexported experimental ECC module:
+    - we now implement `fault_matrix` which gives the mapping between single-qubit physical errors and logical observables.
+    - `MixedDestabilizer` and `Stabilizer` now have constructors when acting on an ECC object.
+- `stab_to_gf2` now works with Pauli operators as well.
+
 ## v0.8.8 - 2023-06-23
 
 - Bump `QuantumInterface` compat.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org>"]
-version = "0.8.8"
+version = "0.8.9"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

docs/make.jl

@@ -9,6 +9,9 @@ using QuantumInterface
 
 #DocMeta.setdocmeta!(QuantumClifford, :DocTestSetup, :(using QuantumClifford); recursive=true)
 
+ENV["LINES"] = 80    # for forcing `displaysize(io)` to be big enough
+ENV["COLUMNS"] = 80
+
 bib = CitationBibliography(joinpath(@__DIR__,"src/references.bib"))
 
 makedocs(

src/QuantumClifford.jl

@@ -1063,7 +1063,7 @@ end
 ##############################
 
 """The F(2,2) matrix of a given tableau, represented as the concatenation of two binary matrices, one for X and one for Z."""
-function stab_to_gf2(s::Tableau)::Matrix{Bool}
+function stab_to_gf2(s::Tableau)
     r, n = size(s)
     H = zeros(Bool,r,2n)
     for iᵣ in 1:r
@@ -1073,6 +1073,14 @@ function stab_to_gf2(s::Tableau)::Matrix{Bool}
     end
     H
 end
+function stab_to_gf2(p::PauliOperator)
+    n = nqubits(p)
+    H = zeros(Bool,2n)
+    @inbounds @simd for i in 1:n
+        H[i], H[i+n] = p[i]
+    end
+    H
+end
 stab_to_gf2(s::Stabilizer) = stab_to_gf2(tab(s))
 
 """Gaussian elimination over the binary field."""

src/ecc/ECC.jl

@@ -2,6 +2,7 @@ module ECC
 
 using QuantumClifford
 using QuantumClifford: AbstractOperation
+import QuantumClifford: Stabilizer, MixedDestabilizer
 
 abstract type AbstractECC end
 
@@ -11,6 +12,9 @@ function encoding_circuit end
 """Parity check tableau of a code."""
 function parity_checks end
 
+Stabilizer(c::AbstractECC) = parity_checks(c)
+MixedDestabilizer(c::AbstractECC) = MixedDestabilizer(Stabilizer(c))
+
 """The number of physical qubits in a code."""
 function code_n end
 
@@ -78,14 +82,195 @@ end
 
 """Logical X operations of a code."""
 function logx_ops(c::AbstractECC)
-    MixedDest = MixedDestabilizer(parity_checks(c))
-    logicalxview(MixedDest)
+    md = MixedDestabilizer(parity_checks(c))
+    logicalxview(md)
 end
 
 """Logical Z operations of a code."""
 function logz_ops(c::AbstractECC)
-    MixedDest = MixedDestabilizer(parity_checks(c))
-    logicalzview(MixedDest)
+    md = MixedDestabilizer(parity_checks(c))
+    logicalzview(md)
+end
+
+"""Error-to-logical-observable map (a.k.a. fault matrix) of a code.
+
+For a code with n physical qubits and k logical qubits this function returns
+a 2k × 2n binary matrix O such that
+`O[i,j]` is true if the logical observable of index `i`
+is flipped by the single physical qubit error of index `j`.
+Indexing is such that:
+
+- `O[1:k,:]` is the error-to-logical-X map
+- `O[k+1:2k,:]` is the error-to-logical-Z map
+- `O[:,1:n]` is the X-physical-error-to-logical map
+- `O[n+1:2n,:]` is the Z-physical-error-to-logical map
+
+E.g. for `k=1`, `n=10`, then
+if `O[2,5]` is true, then the logical Z observable is flipped by a X₅ error;
+and if `O[1,12]` is true, then the logical X observable is flipped by a Z₂ error.
+
+Of note is that there is a lot of freedom in choosing the logical operations!
+A logical operator multiplied by a stabilizer operator is still a logical operator.
+Similarly there is a different fault matrix for each choice of logical operators.
+But once the logical operators are picked, the fault matrix is fixed.
+
+Below we show an example that uses the Shor code. While it is not the smallest code,
+it is a convenient choice to showcase the importance of the fault matrix when dealing
+with degenerate codes where a correction operation and an error do not need to be the same.
+
+First, consider a single-qubit error, potential correction operations, and their effect on the Shor code:
+
+```jldoctest faults_matrix
+julia> using QuantumClifford.ECC: faults_matrix, Shor9
+
+julia> state = MixedDestabilizer(Shor9())
+𝒟ℯ𝓈𝓉𝒶𝒷━━━━━
++ Z________
++ ___Z_____
++ _X_______
++ __X______
++ ____X____
++ _____X___
++ ______X__
++ _______X_
+𝒳ₗ━━━━━━━━━
++ ______XXX
+𝒮𝓉𝒶𝒷━━━━━━━
++ XXX___XXX
++ ___XXXXXX
++ ZZ_______
++ Z_Z______
++ ___ZZ____
++ ___Z_Z___
++ ______Z_Z
++ _______ZZ
+𝒵ₗ━━━━━━━━━
++ Z__Z____Z
+
+julia> err_Z₁ = single_z(9,1) # the error we will simulate
++ Z________
+
+julia> cor_Z₂ = single_z(9,2) # the correction operation we will perform
++ _Z_______
+
+julia> err_Z₁ * state # observe that one of the syndrome bits is now flipped
+𝒟ℯ𝓈𝓉𝒶𝒷━━━━━
++ Z________
++ ___Z_____
++ _X_______
++ __X______
++ ____X____
++ _____X___
++ ______X__
++ _______X_
+𝒳ₗ━━━━━━━━━
++ ______XXX
+𝒮𝓉𝒶𝒷━━━━━━━
+- XXX___XXX
++ ___XXXXXX
++ ZZ_______
++ Z_Z______
++ ___ZZ____
++ ___Z_Z___
++ ______Z_Z
++ _______ZZ
+𝒵ₗ━━━━━━━━━
++ Z__Z____Z
+
+julia> cor_Z₂ * err_Z₁ * state # we are back to a good code state
+𝒟ℯ𝓈𝓉𝒶𝒷━━━━━
++ Z________
++ ___Z_____
+- _X_______
++ __X______
++ ____X____
++ _____X___
++ ______X__
++ _______X_
+𝒳ₗ━━━━━━━━━
++ ______XXX
+𝒮𝓉𝒶𝒷━━━━━━━
++ XXX___XXX
++ ___XXXXXX
++ ZZ_______
++ Z_Z______
++ ___ZZ____
++ ___Z_Z___
++ ______Z_Z
++ _______ZZ
+𝒵ₗ━━━━━━━━━
++ Z__Z____Z
+
+julia> bad_Z₆Z₉ = single_z(9,6) * single_z(9,9) # a different "correction" operation
++ _____Z__Z
+
+julia> bad_Z₆Z₉ * err_Z₁ * state # the syndrome is trivial, but now we have a logical error
+𝒟ℯ𝓈𝓉𝒶𝒷━━━━━
++ Z________
++ ___Z_____
++ _X_______
++ __X______
++ ____X____
+- _____X___
++ ______X__
++ _______X_
+𝒳ₗ━━━━━━━━━
+- ______XXX
+𝒮𝓉𝒶𝒷━━━━━━━
++ XXX___XXX
++ ___XXXXXX
++ ZZ_______
++ Z_Z______
++ ___ZZ____
++ ___Z_Z___
++ ______Z_Z
++ _______ZZ
+𝒵ₗ━━━━━━━━━
++ Z__Z____Z
+```
+
+The success of `cor_Z₂` and the failure of `bad_Z₆Z₉` can be immediately seen through the fault matrix,
+as the wrong "correction" does not result in the same logical flips ad the error:
+
+```jldoctest faults_matrix
+julia> O = faults_matrix(Shor9())
+2×18 BitMatrix:
+ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1
+ 1  0  0  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0
+
+julia> O * stab_to_gf2(err_Z₁)
+2-element Vector{Int64}:
+ 0
+ 0
+
+julia> O * stab_to_gf2(cor_Z₂)
+2-element Vector{Int64}:
+ 0
+ 0
+
+julia> O * stab_to_gf2(bad_Z₆Z₉)
+2-element Vector{Int64}:
+ 1
+ 0
+```
+
+While its use in this situation is rather contrived, the fault matrix is incredibly useful
+when running large scale simulations in which we want a separate fast error sampling process,
+(e.g. with Pauli frames) and a syndrome decoding process, without coupling between them.
+We just gather all our syndrome measurement **and logical observables** from the Pauli frame simulations,
+and then use them with the fault matrix in the syndrome decoding simulation.
+"""
+function faults_matrix(c::AbstractECC)
+    n = code_n(c)
+    k = code_k(c)
+    O = falses(2k, 2n)
+    md = MixedDestabilizer(c)
+    logviews = [logicalxview(md); logicalzview(md)]
+    errors = [one(Stabilizer,n; basis=:X);one(Stabilizer,n)]
+    for i in 1:2k
+        O[i, :] = comm(logviews[i], errors)
+    end
+    return O
 end
 
 # TODO implement isdegenerate

test/test_doctests.jl

@@ -1,6 +1,8 @@
 using Documenter
 using QuantumClifford
 
+ENV["LINES"] = 80    # for forcing `displaysize(io)` to be big enough
+ENV["COLUMNS"] = 80
 @testset "Doctests" begin
     DocMeta.setdocmeta!(QuantumClifford, :DocTestSetup, :(using QuantumClifford); recursive=true)
     doctest(QuantumClifford)