Additional group theory functions for error correction

src/QuantumClifford.jl

@@ -76,6 +76,7 @@ export
     graphstate, graphstate!, graph_gatesequence, graph_gate,
     # Group theory tools
     groupify, minimal_generating_set, pauligroup, normalizer, centralizer, contractor, delete_columns,
+    logical_operator_canonicalize, commutavise, embed, SubsystemCodeTableau,
     # Clipped Gauge
     canonicalize_clip!, bigram, entanglement_entropy,
     # mctrajectories

src/grouptableaux.jl

@@ -1,6 +1,92 @@
 using Graphs
 using LinearAlgebra
 
+"""
+A tableau representation of the logical operator canonical form of a set of Paulis.
+The index of the first operator that commutes with all others is tracked so that the
+stabilizer, destabilizer, logical X, and logical Z aspects of the tableau can be 
+represented.
+"""
+mutable struct SubsystemCodeTableau <: AbstractStabilizer
+    tab::Tableau
+    index::Int
+    r::Int
+    m::Int
+    k::Int
+end
+function SubsystemCodeTableau(t::Tableau)
+    # identity = zero(PauliOperator, nqubits(t))
+    # num = 0
+    # for p in t 
+    #     if p == identity num+=1 end
+    # end
+    index = 1
+    for i in range(1, stop=length(t), step=2)
+        if i + 1 > length(t) 
+            break
+        end
+        if comm(t[i], t[i+1]) == 0x01 
+            index = i+2 # index to split loc into non-commuting pairs and commuting operators
+        end
+    end
+    s = Stabilizer(t[index:length(t)])
+    ind = 1
+    if length(s)>nqubits(s)#if stabilizer is overdetermined, Destabilizer constructor throws error 
+        m = length(s)
+        tab = zero(Tableau, length(t)+length(s), nqubits(t))
+        for i in s 
+            tab[ind] = zero(PauliOperator, nqubits(s))
+            ind+=1
+        end
+    else
+        d = Destabilizer(s)
+        m = length(d)
+        tab = zero(Tableau, length(t)+length(d), nqubits(t))
+        for p in destabilizerview(d) 
+            tab[ind] = p
+            ind+=1
+        end
+    end
+    for i in range(1, stop=index-1, step=2)
+        tab[ind] = t[i]
+        ind+=1
+    end
+    for p in s
+
+            tab[ind] = p
+            ind+=1
+
+    end
+    for i in range(2, stop=index, step=2)
+        tab[ind] = t[i]
+        ind+=1
+    end
+    return SubsystemCodeTableau(tab, index, length(s), m, (index-1)/2)
+
+end
+
+Base.length(t::SubsystemCodeTableau) = length(t.tab)
+
+Base.copy(t::SubsystemCodeTableau) = SubsystemCodeTableau(copy(t.tab))
+
+"""A view of the subtableau corresponding to the stabilizer. See also [`tab`](@ref), [`destabilizerview`](@ref), [`logicalxview`](@ref), [`logicalzview`](@ref)"""
+function stabilizerview(s::SubsystemCodeTableau)
+    return Stabilizer(@view tab(s)[s.m+s.k+1:s.m+s.k+s.r])
+end
+"""A view of the subtableau corresponding to the destabilizer. See also [`tab`](@ref), [`stabilizerview`](@ref), [`logicalxview`](@ref), [`logicalzview`](@ref)"""
+function destabilizerview(s::SubsystemCodeTableau)
+    return Stabilizer(@view tab(s)[1:s.m])
+end
+"""A view of the subtableau corresponding to the logical X operators. See also [`tab`](@ref), [`stabilizerview`](@ref), [`destabilizerview`](@ref), [`logicalzview`](@ref)"""
+function logicalxview(s::SubsystemCodeTableau)
+    return Stabilizer(tab(s)[s.m+1:s.m+s.k])
+end
+"""A view of the subtableau corresponding to the logical Z operators. See also [`tab`](@ref), [`stabilizerview`](@ref), [`destabilizerview`](@ref), [`logicalzview`](@ref)"""
+function logicalzview(s::SubsystemCodeTableau)
+    return Stabilizer(tab(s)[s.m+1+s.k+s.r:length(s)])
+end
+
+
 """
 Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).
 
@@ -18,16 +104,17 @@ 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 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))
+    gen_set = minimal_generating_set(s)
+    n = length(gen_set)::Int
+    group = zero(Tableau, 2^n, nqubits(gen_set))
     for i in 0:2^n-1
         for (digit_order, j) in enumerate(digits(i, base=2, pad=n))
             if j == 1
-                group[i+1] *= s[digit_order]
+                group[i+1] *= gen_set[digit_order]
             end
         end
     end
-    return group 
+    return group
 end
 
 
@@ -59,9 +146,143 @@ function minimal_generating_set(s::Stabilizer)
     end
 end
 
+"""
+For a not-necessarily commutative set of Paulis, returning a generating set of the form 
+<A₁, B₁, A₂, B₂, ... Aₖ, Bₖ, Aₖ₊₁, ... Aₘ> where two operators anticommute if and only if they 
+are of the form Aₖ, Bₖ and commute otherwise.
+
+Returns the generating set as a data structure of type [`SubsystemCodeTableau`](@ref).
+
+```jldoctest
+julia> tab(logical_operator_canonicalize(QuantumClifford.Tableau([P"XX", P"XZ", P"XY"])))
++ Z_
++ XX
+-iX_
++ XZ
+```
+"""
+function logical_operator_canonicalize(t:: QuantumClifford.Tableau)
+    loc = zero(QuantumClifford.Tableau, length(t), nqubits(t))
+    index = 1
+    for i in eachindex(t)
+        for j in eachindex(t)
+            if comm(t[i], t[j]) == 0x01
+                for k in eachindex(t)
+                    if k !=i && k != j
+                        if comm(t[k], t[i]) == 0x01
+                            t[k] = t[j] *t[k]
+                        end
+                        if comm(t[k], t[j]) == 0x01
+                            t[k] = t[i] *t[k]
+                        end
+                    end
+                end
+                if !(t[i] in loc || -1 *t[i] in loc || 1im *t[i] in loc || -1im * t[i] in loc)
+                    loc[index]= t[i]
+                    index+=1
+                end
+                if !(t[j] in loc || -1 *t[j] in loc || 1im *t[j] in loc || -1im * t[j] in loc)
+                    loc[index]= t[j]
+                    index+=1
+                end
+            end
+        end
+    end
+    for i in eachindex(t)
+        if !(t[i] in loc || -1 *t[i] in loc || 1im *t[i] in loc || -1im * t[i] in loc)
+            loc[index]= t[i]
+            index+=1
+        end
+    end
+    while length(loc) > 1 && loc[length(loc)-1] == zero(PauliOperator, nqubits(loc)) 
+        loc = loc[1:(length(loc)-1)]
+    end
+    return SubsystemCodeTableau(loc)
+end
+
+"""
+For a not-necessarily commutative set of Paulis, return the logical operator canonical
+form of that set, and for each pair Aₖ, Bₖ of anticommutative Paulis, add an X to Aₖ,
+a Z to Bₖ, and an I to each other operator, such that the set is now fully commutative 
+as well as return a list of the indices of the added qubits.
+
+Phases are ignored, and returned Stabilizers contain only the + phase
+
+```jldoctest
+julia> commutavise(QuantumClifford.Tableau([P"XX", P"XZ", P"XY"]))[1]
++ XXX
++ XZZ
++ X__
+
+julia> commutavise(QuantumClifford.Tableau([P"XX", P"XZ", P"XY"]))[2]
+1-element Vector{Any}:
+ 3
+```
+"""
+function commutavise(t)
+    loc = QuantumClifford.logical_operator_canonicalize(t)
+    commutative = zero(Stabilizer, 2*loc.k+loc.r, nqubits(loc)+loc.k)
+    ind = 1
+    for i in 1:loc.k
+        dummy = commutative[ind]
+        for j in eachindex(logicalxview(loc)[i]) dummy[j] = logicalxview(loc)[i][j] end
+        dummy[nqubits(loc)+i] = (true, false)
+        commutative[ind] = dummy
+        ind+=1
+    end
+    for i in 1:loc.k
+        dummy = commutative[ind]
+        for j in eachindex(logicalzview(loc)[i]) dummy[j] = logicalzview(loc)[i][j] end
+        dummy[nqubits(loc)+i] = (false, true)
+        commutative[ind] = dummy
+        ind+=1
+    end
+    for i in 1:loc.r
+        dummy = commutative[ind]
+        for j in eachindex(stabilizerview(loc)[i]) dummy[j] = stabilizerview(loc)[i][j] end
+        commutative[ind] = dummy
+        ind+=1   
+    end
+    to_delete = []
+    for i in 1:loc.k
+        push!(to_delete, i+nqubits(tab(loc)))
+    end
+    return commutative, to_delete
+end
+
+"""
+For a given set of Paulis that does not necessarily represent a state, return a set of
+Paulis that contains the commutativised given set and represents a state, as well as deletions 
+needed to return to the original set of Paulis. 
+
+By deleting the qubits in the first returned array, taking the normalizer, then deleting the qubits in the 
+second returned array from the normalizer, the original set is produced. 
+
+```jldoctest
+julia> embed(QuantumClifford.Tableau([P"XX"]))[1]
++ X_X
++ ZZZ
++ XX_
++ ___
+
+julia> embed(QuantumClifford.Tableau([P"XX"]))[2]
+1-element Vector{Any}:
+ 3
+
+julia> embed(QuantumClifford.Tableau([P"XX"]))[3]
+Any[]
+```
+"""
+function embed(t::QuantumClifford.Tableau)
+    com, d1= QuantumClifford.commutavise(t)
+    norm = QuantumClifford.normalizer(com.tab)
+    state, d2 = QuantumClifford.commutavise(norm)
+    return state, d2, d1
+end
+
 """
 Return the full Pauli group of a given length. Phases are ignored by default, 
-but can be included by setting `phases=true`.
+but can be included by setting `phases = true`.
 
 ```jldoctest
 julia> pauligroup(1)
@@ -130,27 +351,37 @@ julia> normalizer(T"X")
 + X
 ```
 """
-function normalizer(t::Tableau)
+function normalizer(t::Tableau; phases=false)
     # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each 
     # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector 
     # converted to `Tableau`.
     n = nqubits(t)
-    pgroup = pauligroup(n, phases=false)
-    ptype = typeof(t[1])
-    normalizer = ptype[]
-    for p in pgroup
+    ptype =typeof(P"I")
+    norm = ptype[]
+
+    p = zero(PauliOperator, n)
+    paulis = ((false, false), (true, false), (false, true), (true, true))
+    for i in Iterators.product(Iterators.repeated(paulis, n)...)
+        QuantumClifford.zero!(p)
+        for (j, k) in enumerate(i)
+            p[j] = k
+        end
         commutes = true
         for q in t
             if comm(p, q) == 0x01
                 commutes = false
             end
         end
         if commutes
-            push!(normalizer, p)
+            push!(norm, copy(p))
+        end
+        if phases
+            for phase in [-1, 1im, -1im]   
+                push!(norm, phase *p)
+            end
         end
     end
-
-    return Tableau(normalizer)
+    return QuantumClifford.Tableau(norm)
 end
 
 """