init (migrating and renaming)

This commit is contained in:
qwjyh 2024-05-07 11:11:07 +09:00
commit f772e27182
19 changed files with 2262 additions and 0 deletions

204
src/PauliOps/PauliOps.jl Normal file
View file

@ -0,0 +1,204 @@
"""
Pauli operator
"""
module PauliOps
using StaticArrays
using IterTools
export SinglePauliOp, PauliOp
export single_pauliop, pauliop
export weight, xweight, zweight
export GeneratedPauliGroup
"""
@enum SinglePauliOp begin
I
X
Y
Z
end
Pauli Operator on a single qubit.
"""
@enum SinglePauliOp begin
I
X
Y
Z
end
"""
PauliOp{N}
Pauli operator on multiple qubits.
"""
const PauliOp{N} = SVector{N, SinglePauliOp}
"""
single_pauliop(char::Char)::SinglePauliOp
Convert `char` to `SinglePauliOp`.
"""
function single_pauliop(char::Char)::SinglePauliOp
char == 'I' && return I
char == 'X' && return X
char == 'Y' && return Y
char == 'Z' && return Z
throw(ArgumentError("invalid char for pauli operator (must be I, X, Y, or Z)"))
end
function Base.String(p::SinglePauliOp)
p == I && return "I"
p == X && return "X"
p == Y && return "Y"
p == Z && return "Z"
end
@inline function single_pauli_product(x::T, y::T)::T where {T <: SinglePauliOp}
if x == I
return y
elseif y == I
return x
elseif x == y
return I
else
x == X && y == Y && return Z
x == X && y == Z && return Y
x == Y && y == Z && return X
return single_pauli_product(y, x)
end
end
Base.:(*)(x::T, y::T) where {T <: SinglePauliOp} = single_pauli_product(x, y)
"""
pauliop(str::AbstractString)::PauliOp
Convert `str` to `PauliOp`.
"""
function pauliop(str::AbstractString)::PauliOp
SVector(single_pauliop.(collect(str))...)
end
Base.String(p::PauliOp) = join(String.(p))
Base.show(io::IO, p::PauliOp) = print(io, "pauliop(\"$(String(p))\")")
Base.summary(io::IO, p::PauliOp) = print(io, "$(length(p))-element PauliOp")
# function Base.show(io::IO, ::MIME"text/plain", p::PauliOp)
# if get(io, :compact, false)
# print(io, String(p))
# else
# summary(io, p)
# print(io, ": ", String(p))
# end
# end
"""
weight(p::PauliOp, [init = 1])
Weight of the operator `p`, i.e. non \$I\$ operator.
"""
function weight(p::PauliOp, init::Integer = 1)
# length(filter(!=(PauliOps.I), p))
count = 0
for i in eachindex(p)
i < init && continue
if p[i] != PauliOps.I
count += 1
end
end
count
end
"""
xweight(p::PauliOp, [init = 1])
Number of \$X, Y\$ in `p`.
"""
function xweight(p::PauliOp, init::Integer = 1)
count = 0
for i in eachindex(p)
i < init && continue
if p[i] == PauliOps.X || p[i] == PauliOps.Y
count += 1
end
end
count
end
"""
zweight(p::PauliOp, [init = 1])
Number of \$Z, Y\$ in `p`.
a"""
function zweight(p::PauliOp, init::Integer = 1)
# length(filter(!=(PauliOps.I), p))
count = 0
for i in eachindex(p)
i < init && continue
if p[i] == PauliOps.Z || p[i] == PauliOps.Y
count += 1
end
end
count
end
# should have compatibility with other packages like AbstractAlgebra?
"""
struct GeneratedPauliGroup
Iterator for group generated from `gens`.
GeneratedPauliGroup(gens::AbstractVector{T}) where {T <: PauliOp}
# Examples
```jldoctest
julia> gens = pauliop.(["IIXX", "IZZI"])
2-element Vector{StaticArraysCore.SVector{4, QuantumLegos.PauliOps.SinglePauliOp}}:
pauliop("IIXX")
pauliop("IZZI")
julia> g = PauliOps.GeneratedPauliGroup(gens)
GeneratedPauliGroup{4}(StaticArraysCore.SVector{4, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXX"), pauliop("IZZI")], IterTools.Subsets{Vector{StaticArraysCore.SVector{4, QuantumLegos.PauliOps.SinglePauliOp}}}(StaticArraysCore.SVector{4, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXX"), pauliop("IZZI")]))
julia> collect(g)
4-element Vector{StaticArraysCore.SVector{4, QuantumLegos.PauliOps.SinglePauliOp}}:
pauliop("IIII")
pauliop("IIXX")
pauliop("IZZI")
pauliop("IZYX")
```
"""
struct GeneratedPauliGroup{N}
gens::AbstractVector{PauliOp{N}}
subsets::IterTools.Subsets
function GeneratedPauliGroup(gens::AbstractVector{T}) where {T <: PauliOp}
length.(gens) |> unique |> length |> ==(1) ||
throw(ArgumentError("All generators must have the same length."))
subsets = IterTools.subsets(gens)
N = length(gens[1])
new{N}(gens, subsets)
end
end
function Base.iterate(g::GeneratedPauliGroup)
next = iterate(g.subsets)
isnothing(next) && return nothing
subset, state = next
init::PauliOp{length(g.gens[1])} = fill(I, length(g.gens[1]))
return (init, state)
end
function Base.iterate(g::GeneratedPauliGroup, state)
ret = iterate(g.subsets, state)
isnothing(ret) && return nothing
subset, state = ret
return (reduce(.*, subset), state)
end
Base.length(g::GeneratedPauliGroup) = length(g.subsets)
Base.eltype(::Type{GeneratedPauliGroup{N}}) where {N} = PauliOp{N}
end # module PauliOps

18
src/QuantumLegos.jl Normal file
View file

@ -0,0 +1,18 @@
module QuantumLegos
export PauliOps
export pauliop, GeneratedPauliGroup, weight
export Lego, LegoLeg, edge, CheckMatrix, checkmatrix, generators, State, add_lego!, add_edge!, distance
using StaticArrays
using IterTools
# PauliOps submodule
include("PauliOps/PauliOps.jl")
using .PauliOps
include("checkmatrix.jl")
include("game.jl")
end

565
src/checkmatrix.jl Normal file
View file

@ -0,0 +1,565 @@
"""
CheckMatrix
# Fields
- `cmat`: matrix itself (matrix of `Bool` since we consider only pauli ops.)
- `nlegs`: number of columns ÷ 2
- `ngens`: number of rows
# Constructor
CheckMatrix(cmat::AbstractMatrix{Bool})
See also [`checkmatrix`](@ref).
"""
mutable struct CheckMatrix
# change this to AbstractMatrix{Bool} to accept BitMatrix?
cmat::Matrix{Bool}
nlegs::Int64
ngens::Int64
function CheckMatrix(cmat::Matrix{Bool}, nlegs::Integer, ngens::Integer)
size(cmat)[1] == ngens ||
throw(ArgumentError("`ngens` must equal to number of row of `cmat`"))
size(cmat)[2] == nlegs * 2 ||
throw(ArgumentError("`nlegs * 2` must equal to number of row of `cmat`"))
new(cmat, nlegs, ngens)
end
end
function CheckMatrix(cmat::AbstractMatrix{Bool})
s = size(cmat)
s[2] % 2 == 0 || throw(ArgumentError("invalid `cmat` size"))
CheckMatrix(cmat, s[2] ÷ 2, s[1])
end
Base.:(==)(x::T, y::T) where {T <: CheckMatrix} = (==)(x.cmat, y.cmat)
Base.copy(x::CheckMatrix) = CheckMatrix(copy(x.cmat))
function Base.show(io::IO, ::MIME"text/plain", cmat::CheckMatrix)
print(io, "CheckMatrix with $(cmat.ngens) generators, $(cmat.nlegs) legs:\n")
display_height, display_width = displaysize(io)
cur_vpos = 1
for row in eachrow(cmat.cmat)
cur_hpos = 1
for (i, elem) in enumerate(row)
sep = i == cmat.nlegs + 1 ? " | " : " "
elem_width = length(sep) + 1
print(io, sep)
if elem
printstyled(io, 1, color = :default)
cur_hpos += elem_width
else
printstyled(io, 0, color = :light_black)
cur_hpos += elem_width
end
if cur_hpos > display_width - 5
print(io, "")
break
end
end
print(io, "\n")
cur_vpos += 1
if cur_vpos > display_height - 6
cur_hpos = 1
for i in 1:(cmat.nlegs * 2)
sep = i == cmat.nlegs + 1 ? " | " : " "
elem_width = length(sep) + 1
print(io, sep)
print(io, "")
cur_hpos += elem_width
if cur_hpos > display_width - 5
print(io, "\n")
break
end
end
return nothing
end
end
end
"""
checkmatrix(stbgens::AbstractVector{T})::CheckMatrix where {T <: PauliOp}
Get [`CheckMatrix`](@ref) from stabilizer generator.
"""
function checkmatrix(stbgens::AbstractVector{T})::CheckMatrix where {T <: PauliOp}
# no stbgens
if length(stbgens) == 0
throw(ArgumentError("No stabilizer is provided. Need at least one."))
end
# check whether all stbgens have the same length
stbgens_length = length.(stbgens)
all(==(stbgens_length[1]), stbgens_length[2:end]) || begin
@error "lengths of stbgens: " length.(stbgens)
throw(ArgumentError("All stabilizer must have the same length"))
end
ngens = length(stbgens)
nlegs = length(stbgens[1])
cmat = zeros(Bool, (ngens, nlegs * 2))
for (i, gen) in enumerate(stbgens)
# get X Z component (can be more optimized)
xs = gen .== PauliOps.X
zs = gen .== PauliOps.Z
ys = gen .== PauliOps.Y
xs = xs .| ys
zs = zs .| ys
cmat[i, :] .= vcat(xs, zs)
end
return CheckMatrix(cmat)
end
"""
xpart(cmat::CheckMatrix)
Get X part (left half) of CheckMatrix.
"""
function xpart(cmat::CheckMatrix)
cmat.cmat[:, 1:(cmat.nlegs)]
end
"""
zpart(cmat::CheckMatrix)
Get Z part (right half) of CheckMatrix.
"""
function zpart(cmat::CheckMatrix)
cmat.cmat[:, (cmat.nlegs + 1):end]
end
function xzparts(cmat::CheckMatrix)
return xpart(cmat), zpart(cmat)
end
"""
generators(cmat::CheckMatrix)::Vector{PauliOp}
Get generators from [`CheckMatrix`](@ref).
"""
function generators(cmat::CheckMatrix)::Vector{PauliOp}
gens = PauliOp{cmat.nlegs}[]
for row in eachrow(cmat.cmat)
gen = map(zip(row[1:(cmat.nlegs)], row[(cmat.nlegs + 1):end])) do (x, z)
g = PauliOps.I
x && (g = g * PauliOps.X)
z && (g = g * PauliOps.Z)
g
end
push!(gens, gen)
end
gens
end
@doc raw"""
direct_sum(cmat_1::T, cmat_2::T)::T where {T <: CheckMatrix}
Returns block diagonal `CheckMatrix` consists of `cmat_1` and `cmat_2`.
# CheckMatrix transform
When adding a lego $l$ with check matrix $H_l$ to a state with check matrix $H_s$,
the resulting check matrix of the state will be
```math
\begin{pmatrix}
H_s & O \\
O & H_l \\
\end{pmatrix}
```
"""
function direct_sum(cmat_1::T, cmat_2::T)::T where {T <: CheckMatrix}
# expand size and fill?
old_x, old_z = xzparts(cmat_1)
add_x, add_z = xzparts(cmat_2)
new_x = cat(old_x, add_x, dims = (1, 2))
new_z = cat(old_z, add_z, dims = (1, 2))
CheckMatrix(hcat(new_x, new_z))
end
"""
eliminate_column!(
cmat::CheckMatrix,
col::Integer,
avoid_row::AbstractVector{T},
)::Union{Nothing, Int64} where {T <: Integer}
Perform Gauss Elimination on `col`.
Keep only one `1` on `col` and remove from other rows by performing multiplication.
Choose row not in `avoid_row`.
If all rows with 1 on `col` are in `avoid_row`, then the last row of them is chose to keep 1.
If all rows on `col` is 0, return nothing.
# Return
Row index which have 1 on `col`.
If all row on `col` is 0, return nothing.
# Example
```jldoctest
julia> ex_cmat = CheckMatrix(Bool[1 0 1 0; 1 1 1 1])
CheckMatrix with 2 generators, 2 legs:
1 0 | 1 0
1 1 | 1 1
julia> QuantumLegos.eliminate_column!(ex_cmat, 1, Int64[])
1
julia> ex_cmat
CheckMatrix with 2 generators, 2 legs:
1 0 | 1 0
0 1 | 0 1
```
"""
function eliminate_column!(
cmat::CheckMatrix,
col::Integer,
avoid_row::AbstractVector{T},
)::Union{Nothing, Int64} where {T <: Integer}
all(avoid_row .≤ cmat.ngens) || throw(ArgumentError("avoid_row is out of index"))
@assert col cmat.nlegs * 2
row_ids = findall(cmat.cmat[:, col]) # row with 1 on `col`
@debug "ids of rows with 1 on col $(col) = $(row_ids)"
if isempty(row_ids)
@debug "No 1 on column $(col)."
return nothing
end
if length(row_ids) == 1
# already have only 1 1
@debug "Already have only 1 1. Nothing to do."
return row_ids[1]
else
@debug "More than 1 1. Need some operation."
row_id_id = findfirst(!(avoid_row), row_ids)
row_id = if isnothing(row_id_id)
@debug "Rows with 1 are duplicated with avoid_row."
# select last row_ids
pop!(row_ids)
else
popat!(row_ids, row_id_id)
end
@debug "Selected row:$(row_id) to keep 1 on col:$(col)."
for row in row_ids
cmat.cmat[row, :] .⊻= cmat.cmat[row_id, :]
end
@debug "Finished. now only row:$(row_id) has 1 on col:$(col)"
return row_id
end
error("Unreachable")
end
"""
eliminate_column!(cmat::CheckMatrix, col::Integer,
avoid_row::AbstractVector{T}) where {T <: Union{Nothing, Integer}}
`avoid_row` can include `Nothing`, which is ignored in actual evaluation.
"""
function eliminate_column!(
cmat::CheckMatrix,
col::Integer,
avoid_row::AbstractVector{T},
) where {T <: Union{Nothing, Integer}}
eliminate_column!(cmat, col, Vector{Integer}(filter(!isnothing, avoid_row)))
end
@inline function swap_row!(m::AbstractMatrix, i::Integer, j::Integer)
i == j && return nothing
@inbounds m[i, :], m[j, :] = m[j, :], m[i, :]
end
"""
align_row!(m::AbstractMatrix, row::Integer, occupied::Vector{Union{Nothing, Integer}}) -> Integer
align_row!(m::AbstractMatrix, row::Nothing, occupied::Vector{Union{Nothing, Integer}}) -> Nothing
Swap row at `row` in `m` and row at next to the maximum in `occupied`.
`occupied` is supposed to be a list of returns from [`eliminate_column!`](@ref).
If `row` is in `occupied`, do nothing and returns `row`.
If `row` is `nothing`, return `nothing`.
# Arguments
- `m::AbstractMatrix`: mutated
- `row::Union{Nothing, Integer}`: row to be aligned
- `occupied::Vector{Union{Nothing, Integer}}`: indices of already occupied rows. `row` will be next to these rows.
# Return
- Row index where `row` is moved.
"""
function align_row! end
function align_row!(
m::AbstractMatrix,
row::Integer,
occupied::AbstractVector{T},
) where {T <: Integer}
if row occupied
return row
end
target = if isempty(occupied)
1 # begin
else
maximum(occupied) + 1 # TODO: might cause out of index
end
swap_row!(m, row, target)
target
end
function align_row!(
_::AbstractMatrix,
row::Nothing,
_::AbstractVector{T},
) where {T <: Union{Nothing, Integer}}
nothing
end
function align_row!(
m::AbstractMatrix,
row::Integer,
occupied::AbstractVector{T},
) where {T <: Union{Nothing, Integer}}
filter!(!isnothing, occupied)
occupied = Vector{Integer}(occupied)
align_row!(m, row, occupied)
end
# TODO: improve perf (see README.md)
"""
ref!(cmat::CheckMatrix) -> Int
Convert `cmat` to row echelon form.
Returns rank of check matrix.
# Examples
```jldoctest
julia> cmat = CheckMatrix(Bool[
1 0 1 0 1 1 0 1
0 1 0 0 0 1 0 0
1 1 1 0 1 0 1 1
0 1 0 0 0 1 0 0
])
CheckMatrix with 4 generators, 4 legs:
1 0 1 0 | 1 1 0 1
0 1 0 0 | 0 1 0 0
1 1 1 0 | 1 0 1 1
0 1 0 0 | 0 1 0 0
julia> QuantumLegos.ref!(cmat)
3
julia> cmat
CheckMatrix with 4 generators, 4 legs:
1 0 1 0 | 1 1 0 1
0 1 0 0 | 0 1 0 0
0 0 0 0 | 0 0 1 0
0 0 0 0 | 0 0 0 0
```
"""
function ref!(cmat::CheckMatrix)
# TODO: remain cmat shape when rank is max
# For manual debugging stabilizer generators
# Requires LinearAlgebra
# if cmat.ngens == rank(cmat.cmat)
# return cmat.ngens
# end
r = 0 # row
# for col in eachcol(cmat.cmat)
for i in 1:(cmat.nlegs * 2)
isall0 = true
rid = 0
for k in (r + 1):(cmat.ngens)
if cmat.cmat[k, i]
isall0 = false
rid = k
break
end
end
if isall0
continue
end
r += 1
swap_row!(cmat.cmat, r, rid)
for k in (rid + 1):(cmat.ngens)
if cmat.cmat[k, i]
@. cmat.cmat[k, :] = cmat.cmat[r, :] cmat.cmat[k, :]
end
end
#
# rids = Int[]
# for j in (r + 1):(cmat.ngens)
# cmat.cmat[j, i] && push!(rids, j)
# end
# # rids = findall(col) |> filter(>(r))
# if isempty(rids) # all 0
# continue
# end
# r += 1
# rid = popfirst!(rids)
# swap_row!(cmat.cmat, r, rid)
# if isempty(rids) # only 1 1
# continue
# end
# for i in rids
# @. cmat.cmat[i, :] = cmat.cmat[r, :] ⊻ cmat.cmat[i, :]
# end
if r == cmat.ngens
break
end
end
r # rank
end
"""
eliminate_dependent_row!(cmat::CheckMatrix) -> CheckMatrix
Remove dependent rows to keep only independent generators.
"""
function eliminate_dependent_row!(cmat::CheckMatrix)
# convert to reduced row echelon form: ref!
# remove all-zero rows
r = ref!(cmat)
num_zero = cmat.ngens - r
if num_zero > 0
cmat.cmat = cmat.cmat[1:r, :]
cmat.ngens = r
end
return cmat
end
"""
self_trace!(cmat::CheckMatrix, col_1::Integer, col_2::Integer)
Take a self-trace of checkmatrix `cmat` with column `col_1` and `col_2`.
# Example
TODO
"""
function self_trace!(cmat::CheckMatrix, col_1::Integer, col_2::Integer)#::CheckMatrix
if !(col_1 cmat.nlegs && col_2 cmat.nlegs)
throw(
ArgumentError(
"Invalid column index(specified index is too large: $(max(col_1, col_2)) > $(cmat.nlegs))",
),
)
end
if col_1 == col_2
throw(ArgumentError("Can't trace two same legs"))
end
# sort to col_1 < col_2
if col_1 > col_2
col_1, col_2 = col_2, col_1
end
@debug "Initial cmat" cmat
# TODO: cmat.nlegs ≤ 2 ?
# organize cmat to where col_1 and col_2 have less than 3 true
## do for X on col_1
col_1_x_id = eliminate_column!(cmat, col_1, Int64[])
col_1_x_id = align_row!(cmat.cmat, col_1_x_id, Int64[])
@debug "Finished eliminating of X" cmat
## do for Z on col_1
col_1_z_id = eliminate_column!(cmat, col_1 + cmat.nlegs, [col_1_x_id])
col_1_z_id = align_row!(cmat.cmat, col_1_z_id, [col_1_x_id])
@debug "Finished elimination of Z" cmat
@debug "col_1 X:" findall(cmat.cmat[:, col_1])
@debug "col_1 Z:" findall(cmat.cmat[:, col_1 + cmat.nlegs])
## repeat for col_2
## for X on col_2
col_2_x_id = eliminate_column!(cmat, col_2, [col_1_x_id, col_1_z_id])
col_2_x_id = align_row!(cmat.cmat, col_2_x_id, [col_1_x_id, col_1_z_id])
@debug "Finished eliminating of X" cmat
## do for Z on col_1
col_2_z_id =
eliminate_column!(cmat, col_2 + cmat.nlegs, [col_1_x_id, col_1_z_id, col_2_x_id])
col_2_z_id = align_row!(cmat.cmat, col_2_z_id, [col_1_x_id, col_1_z_id, col_2_x_id])
@debug "Finished elimination of Z" cmat
@debug "col_2 X:" findall(cmat.cmat[:, col_2])
@debug "col_2 Z:" findall(cmat.cmat[:, col_2 + cmat.nlegs])
@debug "Finished Gauss Elimination" cmat
# @info "col x/z s $([col_1_x_id, col_1_z_id, col_2_x_id, col_2_z_id])"
# then adding them
# select rows which have same value on col_1 and col_2
# trace (remove col_1, col_2) and remove duplicate row(or dependent row: i.e. row-elim.)
@debug "Tracing"
# assert that col_1, col_2 have 1s only on row 1, 2, 3, 4
@assert findall(cmat.cmat[:, col_1]) .|> (1) |> all
@assert findall(cmat.cmat[:, col_1 + cmat.nlegs]) .|> (2) |> all
@assert findall(cmat.cmat[:, col_2]) .|> (3) |> all
@assert findall(cmat.cmat[:, col_2 + cmat.nlegs]) .|> (4) |> all
col_12_xz_ids = [col_1_z_id, col_1_x_id, col_2_x_id, col_2_z_id] |> filter(!isnothing)
# note that col_1 < col_2
remaining_x_cols =
[1:(col_1 - 1)..., (col_1 + 1):(col_2 - 1)..., (col_2 + 1):(cmat.nlegs)...]
remaining_cols = [remaining_x_cols..., (remaining_x_cols .+ cmat.nlegs)...]
if cmat.ngens 3
# TODO:
# @warn "cmat has ≤ 3 generators" cmat
subsets_row = eachrow(cmat.cmat) |> subsets |> collect
filter!(!isempty, subsets_row)
generated_rows = map(subsets_row) do vs
reduce(.⊻, vs)
end
# matching on col
filter!(generated_rows) do v
v[[col_1, col_2]] == v[[col_1 + cmat.nlegs, col_2 + cmat.nlegs]]
end
new_ngens = length(generated_rows)
new_cmat = zeros(Bool, (new_ngens, 2 * (cmat.nlegs - 2)))
for (i, v) in enumerate(generated_rows)
new_cmat[i, :] .= v[remaining_cols]
end
cmat.cmat = new_cmat
cmat.nlegs -= 2
@assert size(cmat.cmat)[2] == 2 * cmat.nlegs "cmat size mismatched"
cmat.ngens = size(cmat.cmat)[1]
@assert cmat.ngens == new_ngens
elseif cmat.cmat[1, col_1] &&
cmat.cmat[2, col_1 + cmat.nlegs] &&
cmat.cmat[3, col_2] &&
cmat.cmat[4, col_2 + cmat.nlegs]
@assert Set([col_1_x_id, col_1_z_id, col_2_x_id, col_2_z_id]) == Set([1, 2, 3, 4])
# All errors correctable (D.9) (row 2, 3 swapped)
cmat.cmat[1, :] .⊻= cmat.cmat[3, :]
cmat.cmat[2, :] .⊻= cmat.cmat[4, :]
@assert cmat.cmat[1, col_1] && cmat.cmat[3, col_2] "Test whether the form is D.10"
@assert cmat.cmat[2, col_1 + cmat.nlegs] && cmat.cmat[4, col_2 + cmat.nlegs] "Test whether the form is D.10"
new_cmat = zeros(Bool, (cmat.ngens - 2, 2 * (cmat.nlegs - 2)))
new_cmat[1, :] .= cmat.cmat[1, remaining_cols]
new_cmat[2, :] .= cmat.cmat[2, remaining_cols]
new_cmat[3:end, :] .= cmat.cmat[5:end, remaining_cols]
cmat.cmat = new_cmat
cmat.nlegs -= 2
cmat.ngens -= 2
else#if maximum(col_12_xz_ids) == 3 # TODO: can be split
# TODO need some cases
# @warn "Implementing..."
# @info cmat cmat
row_123 = eachrow(cmat.cmat)[1:3]
subsets_from_123 = collect(subsets(row_123))
filter!(!isempty, subsets_from_123)
generated_from_123 = map(subsets_from_123) do vs
reduce(.⊻, vs)
end
filter!(generated_from_123) do v
# matched
v[[col_1, col_2]] == v[[col_1 + cmat.nlegs, col_2 + cmat.nlegs]]
end
n_generated_from_123 = length(generated_from_123)
new_ngens = n_generated_from_123 + cmat.ngens - 3
new_cmat = zeros(Bool, (new_ngens, 2 * (cmat.nlegs - 2)))
for (i, v) in enumerate(generated_from_123)
new_cmat[i, :] .= v[remaining_cols]
end
new_cmat[(n_generated_from_123 + 1):end, :] .= cmat.cmat[4:end, remaining_cols]
cmat.cmat = new_cmat
cmat.nlegs -= 2
@assert cmat.nlegs == size(new_cmat)[2] ÷ 2
cmat.ngens = size(new_cmat)[1]
end
# ngens is updated below
eliminate_dependent_row!(cmat)
return cmat
end

281
src/game.jl Normal file
View file

@ -0,0 +1,281 @@
macro todo()
error("TODO: Unimplemented")
end
macro todo(msg)
error("TODO: Unimplemented: $msg")
end
"""
Quantum lego with `N` legs.
# Fields
- `nlegs::Int64`: number of legs, equals `N`
- `stabgens::SVector{N, PauliOp{N}}`: stabilizer generators. vector of [`PauliOp`](@ref)
# Constructor
Lego([nlegs::Integer], stabgens::AbstractVector{PauliOp{N}})
Constructor for [`Lego`](@ref).
`nlegs` is optional (default is length of the first stabilizer generator).
# Example
```jldoctest
julia> stabgens = pauliop.(["II", "XX"])
2-element Vector{StaticArraysCore.SVector{2, QuantumLegos.PauliOps.SinglePauliOp}}:
pauliop("II")
pauliop("XX")
julia> Lego(stabgens)
Lego{2}(2, StaticArraysCore.SVector{2, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("II"), pauliop("XX")])
```
"""
struct Lego{N}
nlegs::Int64
stabgens::SVector{N, PauliOp{N}} # There two Ns must be the same?
function Lego(nlegs::Integer, stabgens::AbstractVector{PauliOp{N}}) where {N}
all(length(stabgens[1]) .== length.(stabgens)) ||
throw(ArgumentError("All stabgens have the same length"))
nlegs == length(stabgens[1]) ||
throw(ArgumentError("`nlegs` must equal to length of stabilizer"))
new{N}(nlegs, stabgens)
end
end
Lego(stabgens::AbstractVector{PauliOp{N}}) where {N} = Lego(length(stabgens[1]), stabgens)
function nlegs(lego::Lego)::Integer
lego.nlegs
end
"""
mutable struct State
To be used in [`State`](@ref).
# Fields
- `lego_id::Int64`: index in `legos` in `State`
- `edge_id::Int64`: index in `Lego` in `legos` in `State`. No validation check included in `LegoLeg`.
# Example
```jldoctest
julia> x = LegoLeg.([(2, 1), (1, 1), (1, 0)])
3-element Vector{LegoLeg}:
LegoLeg(2, 1)
LegoLeg(1, 1)
LegoLeg(1, 0)
julia> sort(x)
3-element Vector{LegoLeg}:
LegoLeg(1, 0)
LegoLeg(1, 1)
LegoLeg(2, 1)
```
"""
struct LegoLeg
lego_id::Int64
leg_id::Int64
end
LegoLeg(t::Tuple{Integer, Integer}) = LegoLeg(t...)
"""
Helper function to create `Tuple{LegoLeg, LegoLeg}` to represent edge.
"""
function edge end
"""
edge(t::Tuple{T, T}) where {T <: Tuple{Integer, Integer}}
"""
function edge(t::Tuple{T, T}) where {T <: Tuple{Integer, Integer}}
(LegoLeg(t[1]), LegoLeg(t[2]))
end
"""
edge(t::Tuple{T, T, T, T}) where {T <: Integer}
"""
function edge(t::Tuple{T, T, T, T}) where {T <: Integer}
edge(((t[1], t[2]), (t[3], t[4])))
end
"""
edge(x::T, y::T, z::T, w::T) where {T <: Integer}
"""
function edge(x::T, y::T, z::T, w::T) where {T <: Integer}
edge(((x, y), (z, w)))
end
function Base.isless(x::T, y::T) where {T <: LegoLeg}
if x.lego_id != y.lego_id
return x.lego_id < y.lego_id
else
return x.leg_id < y.leg_id
end
end
"""
mutable struct State
State (in p.4)
# Fields
- `legos`: `Vector{Lego}`
- `edges`: Vector of ((`lego_i, leg_n`), (`lego_j, leg_m`)).
Each element is sorted (i.e. `lego_i < lego_j` or `lego_i == lego_j && leg_n < leg_m`).
This feature is used in [`is_connected_to_firstlego`](@ref).
- `cmat::CheckMatrix`: CheckMatrix
# Constructor
State(legos::Vector{Lego{N}}, edges::Vector{Tuple{LegoLeg, LegoLeg}})
# Methods with
- [`add_lego!`](@ref)
- [`add_edge!`](@ref)
# Example
TODO
"""
mutable struct State
legos::Vector{Lego}
edges::Vector{Tuple{LegoLeg, LegoLeg}}
cmat::CheckMatrix
function State(legos::Vector{Lego{N}}, edges::Vector{Tuple{LegoLeg, LegoLeg}}) where {N}
if length(edges) == 0
if length(legos) == 0
throw(ArgumentError("Need at least one lego"))
elseif length(legos) == 1
cmat = checkmatrix(legos[1].stabgens)
return new(legos, edges, cmat)
else
# just iterate instead of recursion?
return add_lego!(State(legos[1:(end - 1)], edges), legos[end])
end
else
new_edge = pop!(edges)
if new_edge[1] == new_edge[2]
throw(ArgumentError("Can't make edges between the same leg."))
end
# sort
if new_edge[1] > new_edge[2]
new_edge = (new_edge[2], new_edge[1])
end
state = State(legos, edges)
return add_edge!(state, new_edge)
end
end
end
function Base.:(==)(x::T, y::T) where {T <: State}
x.legos == y.legos && x.edges == y.edges && x.cmat == y.cmat
end
"""
add_lego!(state::State, lego::Lego) -> State
Add a new lego, updating `state`.
"""
function add_lego!(state::State, lego::Lego)::State
push!(state.legos, lego)
# direct sum
state.cmat = direct_sum(state.cmat, checkmatrix(lego.stabgens))
return state
end
"""
cmat_index(state::State, leg::LegoLeg)::Int64
Get column index corresponds to `leg` in check matrix of `state`.
If given `leg` is already connected, it throws `ArgumentError`.
If given `lego_id` of `leg` is out of `state.legos`, throws `ArgumentError`.
"""
function cmat_index(state::State, leg::LegoLeg)::Int64
connected_legs = if isempty(state.edges)
LegoLeg[]
else
mapreduce(x -> [x...], vcat, state.edges)
end
if leg in connected_legs
throw(ArgumentError("The specified leg:$(leg) is already connected."))
end
if leg.lego_id > length(state.legos)
throw(ArgumentError("state doesn't have lego:$(leg.lego_id)"))
end
sort!(connected_legs)
filter!(<(leg), connected_legs)
n_connected_edges = length(connected_legs)
n_lego_legs = state.legos[1:(leg.lego_id - 1)] .|> nlegs |> sum
return n_lego_legs + leg.leg_id - n_connected_edges
end
"""
add_edge!(state::State, leg_1::LegoLeg, leg_2::LegoLeg)::State
Add a new edge between `leg_1` and `leg_2`, updating `state`.
"""
function add_edge!(state::State, leg_1::LegoLeg, leg_2::LegoLeg)::State
# sort
if leg_1 > leg_2
leg_1, leg_2 = leg_2, leg_1
end
col_1 = cmat_index(state, leg_1)
col_2 = cmat_index(state, leg_2)
self_trace!(state.cmat, col_1, col_2) # mutates cmat
push!(state.edges, (leg_1, leg_2))
return state
end
add_edge!(state::State, edge::Tuple{LegoLeg, LegoLeg}) = add_edge!(state, edge...)
"""
is_connected_to_firstlego(state::State)::BitVector
Returns vector which stores whether each lego is connected to the first lego.
"""
function is_connected_to_firstlego(state::State)::BitVector
connected_to_first = falses(state.legos |> length)
connected_to_first[1] = true
sort!(state.edges)
for edge in state.edges
is_edge_1_connected = connected_to_first[edge[1].lego_id]
is_edge_2_connected = connected_to_first[edge[2].lego_id]
if is_edge_1_connected && is_edge_2_connected
continue
elseif is_edge_1_connected
connected_to_first[edge[2].lego_id] = true
elseif is_edge_2_connected
connected_to_first[edge[1].lego_id] = true
else
continue
end
end
return connected_to_first
end
"""
distance(state::State) -> Int
Calculate code distance when the first leg of `state` is assigned as logical.
"""
distance(state::State) = _naive_distance(state)
# TODO
"""
Calculate distance.
Use minimum distance of all generated normalizers.
"""
function _naive_distance(state::State)
# not optimized(can use less alloc)?
_naive_distance(state.cmat)
end
function _naive_distance(cmat::CheckMatrix)
# not optimized(can use less alloc)?
if cmat.ngens != cmat.nlegs
return 0
end
normalizers = cmat |> generators |> GeneratedPauliGroup |> collect
filter!(x -> x[1] != PauliOps.I, normalizers)
isempty(normalizers) && @warn "No stabilizer in the code" state
distance::Integer = minimum(weight, normalizers) - 1 # substitute 1 because first op is always ≠I
@assert distance >= 0
return distance
end