commit f772e2718285d0c83c8957fb8746f626ab9cf540 Author: qwjyh Date: Tue May 7 11:11:07 2024 +0900 init (migrating and renaming) diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..95731a5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.jl.*.cov +*.jl.cov +*.jl.mem +/Manifest.toml +/docs/Manifest.toml +/docs/build/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..badeb32 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 qwjyh and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..48713a7 --- /dev/null +++ b/Project.toml @@ -0,0 +1,31 @@ +name = "QuantumLegos" +uuid = "6c892f99-6e80-4382-8dc3-97545b5cb80e" +authors = ["qwjyh "] +version = "0.1.0-DEV" + +[deps] +IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[compat] +AbstractAlgebra = "0.34" +Aqua = "0.8" +Documenter = "1" +IterTools = "1.8" +JET = "0.8" +LinearAlgebra = "1" +StaticArrays = "1" +Test = "1" +julia = "1.9" + +[extras] +AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["AbstractAlgebra", "Aqua", "DataStructures", "Documenter", "JET", "LinearAlgebra", "Test"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..78e2447 --- /dev/null +++ b/README.md @@ -0,0 +1,64 @@ +# Legos + +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://qwjyh.github.io/QuantumLegos.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://qwjyh.github.io/QuantumLegos.jl/dev/) +[![Build Status](https://github.com/qwjyh/QuantumLegos.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/qwjyh/QuantumLegos.jl/actions/workflows/CI.yml?query=branch%3Amain) + +[![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) + +## How-Tos for Julia beginners +### How to get Julia +Download binary from https://julialang.org/downloads/. +Or use [juliaup](https://github.com/JuliaLang/juliaup)(recommended). + +### How to run (for dev.)(init) +1. clone this repo +2. cd to here +3. `julia --project` +4. type `]` to enter Pkg mode in REPL. +5. `instantiate` to download all dependencies +6. press backspace to exit Pkg mode +7. `using Legos` to use (normally starts from here) + +Recommended tools for REPL: +- OhMyREPL: for syntax highlight and interactive history search in REPL +- Revise: for auto-loading changes to the package + +### How to build and browse docs locally(init) +1. cd to `docs/` +2. `julia --project` +3. `]` to enter Pkg mode +5. `Pkg> dev ..` to add `Legos` to available package +4. `Pkg> instantiate` +6. Exit Pkg mode(backspace) and `include("make.jl")` to build. Now you can exit julia REPL by `exit()` or `Ctrl-D`. +7. Alternatively, run `julia --project make.jl` from any shell. +8. To browse file on browser, I recommend either to use LiveServer.jl (from Pkg mode REPL, `add LiveServer`) and `julia -e 'using LiveServer; serve(dir="docs/build")'` or cd to `docs/build` and `python -m http.server --bind localhost`. Details on [Note on Documenter.jl docs](https://documenter.juliadocs.org/stable/man/guide/#Building-an-Empty-Document) + +- You can stop the local server by just stopping LiveServer or python http.server (just `Ctrl-C`) +- You can update the docs by `julia --project make.jl` on `docs/` and reload on your browser. + +## TODO +- [x] write test for CheckMatrix constructor + - [x] add `Base.==` method for CheckMatrix +- [x] implement checkmatrix constructor from stabilizer generator +- [x] implement constructor for State + - [x] implement initial constructor for State with only 1 lego and 0 edge + - [x] implement action function on State + - [x] implement tracing function for CheckMatrix + - [ ] test for self-tracing + - [x] implement map function from LegoLeg to checkmatrix index +- [x] implement functions to glean stabilizer generator from CheckMatrix +- [ ] test functions with examples on the paper +- [x] improve perf of ref! +- [ ] implement function to calculate enumerator polynomial + +## ref! optimization +before + +![Benchmark of ref! vs AbstractAlgebra.rref](./docs/src/img/bench_ref_vs_aarref.jpg) + +(so many allocs, take profile) + +after + +![Optimized](./docs/src/img/bench_ref_optimized.jpg) diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..a5252db --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +Legos = "6c892f99-6e80-4382-8dc3-97545b5cb80e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..3599446 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,23 @@ +#! format: off +using Legos +using Documenter + +DocMeta.setdocmeta!(Legos, :DocTestSetup, :(using Legos); recursive=true) + +makedocs(; + modules=[Legos], + authors="", + # repo="", + sitename="Legos.jl", + format=Documenter.HTML(; + prettyurls=get(ENV, "CI", "false") == "true", + edit_link="main", + assets=String[], + ), + pages=[ + "Home" => "index.md", + "PauliOps" => "pauliops.md", + "checkmatrix.md", + "distance.md", + ], +) diff --git a/docs/src/checkmatrix.md b/docs/src/checkmatrix.md new file mode 100644 index 0000000..4d6fa4e --- /dev/null +++ b/docs/src/checkmatrix.md @@ -0,0 +1,82 @@ +# Details on check matrix operations + +## Basic features of check matrix + +- All action on state can be represented as operations on check matrix. +- QuantumLegos which don't connect each other are represented by block diagonal check matrix. + +Thus, one needs only one check matrix to represent state. + +- Generated group from generators represented by check matrix is invariant under row-swap and row-multiplying transforms on the check matrix. + +Thus, one can freely perform these row operations on check matrices. + +## Construction of check matrices and retrieving features + +Use [`CheckMatrix`](@ref) from `Matrix{Bool}` and [`checkmatrix`](@ref) from `Vector{PauliOp}`. + +```@docs +CheckMatrix +checkmatrix +``` + +Get xpart and zpart of check matrix. + +```@docs +QuantumLegos.xpart +QuantumLegos.zpart +``` + +Get generators from check matrix with [`generators`](@ref). + +```@docs +generators +``` + +See also, [`GeneratedPauliGroup`](@ref) and [`Base.Set`](https://docs.julialang.org/en/v1/base/collections/#Base.Set) to get group generated by generators. + +## Direct sum + +Used when adding lego without edge. + +```@docs +QuantumLegos.direct_sum +``` + +## Self-tracing + +### self-tracing + +```@docs +QuantumLegos.self_trace! +``` + +During self-tracing, pre-formatting and post-formatting described below is performed. + +### pre-formatting + +```@docs +QuantumLegos.eliminate_column! +``` + +```@docs +QuantumLegos.align_row! +``` + +[`QuantumLegos.align_row!`](@ref) + +By this process, all columns to be traced have `1`s on only 1st to 4th rows. +So during [`QuantumLegos.self_trace!`](@ref), all stabilizers generated from top three rows are calculated and perform operator matching. + +### post-formatting + +Remove rows which are linear combinations of other rows. + +```@docs +QuantumLegos.ref! +``` + +```@docs +QuantumLegos.eliminate_dependent_row! +``` + diff --git a/docs/src/distance.md b/docs/src/distance.md new file mode 100644 index 0000000..2df98ce --- /dev/null +++ b/docs/src/distance.md @@ -0,0 +1,236 @@ +# How to calculate code distance from the state. + +!!! warning "WIP" + This document is not fully completed. + +## Definition of code distance. + +Let's consider encoding circuit with $1$ logical bit and $k$ physical qubits[^1]. +Then this encoding has two physical basis, $\ket{0}_L$ and $\ket{1}_L$. +The **distance** of this encoding is the minimum bit flip required to convert between $\ket{0}_L$ and $\ket{1}_L$. + +[^1]: Not all state can be formalized like this. TODO + +## Classification of the stabilizers. + +When treating `State`, the logical leg is not assigned and one can treat all stabilizers equally. +However, if logical leg is assigned to the state, these stabilizers can be classified to $4$ groups. + +1. stabilizers on physical qubits +2. ​$\bar{X}$, which corresponds to logical $X$ +3. ​$\bar{Z}$, which corresponds to logical $Z$ +4. ​$\bar{Y}$, which corresponds to logical $Y$ + +Let $\ket{V}$ is the dual state of the channel or encoding map $[[n, 1, d]]$, + +```math +distance = \min_{S \in stabilizers} \#\left\{ i \mid \bar{Z}_i ≠ S_i \right\} +``` + +## Calculating code distance from the check matrix. + +TODO: nor required if the performance doesn't matter. + +## Examples + +### $[[5, 1, 3]]$ code + +​$[[5, 1, 3]]$ code has $4$ stabilizers generators, $XZZXI, IXZZX, XIXZZ, ZXIXZ$ and $2$ logical operators, $\bar{X} = XXXXX$ and $\bar{Z} = ZZZZZ$. +Therefore, stabilizer generators for the corresponding state $[[6, 0]]$ is $IXZZXI, IIXZZX, IXIXZZ, IZXIXZ, XXXXXX, ZZZZZZ$. + +Let's construct $[[6, 0]]$ state on QuantumLegos.jl. +```jldoctest 1 +julia> using QuantumLegos + +julia> stab_513 = pauliop.(["IXZZXI", "IIXZZX", "IXIXZZ", "IZXIXZ", "XXXXXX", "ZZZZZZ"]) +6-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IXZZXI") + pauliop("IIXZZX") + pauliop("IXIXZZ") + pauliop("IZXIXZ") + pauliop("XXXXXX") + pauliop("ZZZZZZ") + +julia> lego_513 = Lego(stab_513) +Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IXZZXI"), pauliop("IIXZZX"), pauliop("IXIXZZ"), pauliop("IZXIXZ"), pauliop("XXXXXX"), pauliop("ZZZZZZ")]) + +julia> state_513 = State([lego_513], edge.([])) +State(Lego[Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IXZZXI"), pauliop("IIXZZX"), pauliop("IXIXZZ"), pauliop("IZXIXZ"), pauliop("XXXXXX"), pauliop("ZZZZZZ")])], Tuple{LegoLeg, LegoLeg}[], CheckMatrix(Bool[0 1 … 0 0; 0 0 … 1 0; … ; 1 1 … 0 0; 0 0 … 1 1], 6, 6)) + +``` + +Then collect generators of the state. +```jldoctest 1 +julia> normalizers = state_513.cmat |> generators |> GeneratedPauliGroup |> collect +64-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IIIIII") + pauliop("IXZZXI") + pauliop("IIXZZX") + pauliop("IXYIYX") + pauliop("IXIXZZ") + pauliop("IIZYYZ") + pauliop("IXXYIY") + pauliop("IIYXXY") + pauliop("IZXIXZ") + pauliop("IYYZIZ") + ⋮ + pauliop("YYIZZI") + pauliop("YXZYZX") + pauliop("YIIXYX") + pauliop("YXYXII") + pauliop("YIXYXI") + pauliop("YIZZIY") + pauliop("YXIIXY") + pauliop("YIYIZZ") + pauliop("YXXZYZ") + +``` + +Get stabilizer and normalizers of the $[[5, 1, 3]]$ code by assigning the first leg as logical. +```jldoctest 1 +julia> stabs = filter(x -> x[1] == PauliOps.I, normalizers) +16-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IIIIII") + pauliop("IXZZXI") + pauliop("IIXZZX") + pauliop("IXYIYX") + pauliop("IXIXZZ") + pauliop("IIZYYZ") + pauliop("IXXYIY") + pauliop("IIYXXY") + pauliop("IZXIXZ") + pauliop("IYYZIZ") + pauliop("IZIZYY") + pauliop("IYZIZY") + pauliop("IYXXYI") + pauliop("IZYYZI") + pauliop("IYIYXX") + pauliop("IZZXIX") + +julia> norm_x = filter(x -> x[1] == PauliOps.X, normalizers) +16-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("XXXXXX") + pauliop("XIYYIX") + pauliop("XXIYYI") + pauliop("XIZXZI") + pauliop("XIXIYY") + pauliop("XXYZZY") + pauliop("XIIZXZ") + pauliop("XXZIIZ") + pauliop("XYIXIY") + pauliop("XZZYXY") + pauliop("XYXYZZ") + pauliop("XZYXYZ") + pauliop("XZIIZX") + pauliop("XYZZYX") + pauliop("XZXZII") + pauliop("XYYIXI") + +julia> norm_y = filter(x -> x[1] == PauliOps.Y, normalizers) +16-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("YYYYYY") + pauliop("YZXXZY") + pauliop("YYZXXZ") + pauliop("YZIYIZ") + pauliop("YZYZXX") + pauliop("YYXIIX") + pauliop("YZZIYI") + pauliop("YYIZZI") + pauliop("YXZYZX") + pauliop("YIIXYX") + pauliop("YXYXII") + pauliop("YIXYXI") + pauliop("YIZZIY") + pauliop("YXIIXY") + pauliop("YIYIZZ") + pauliop("YXXZYZ") + +julia> norm_z = filter(x -> x[1] == PauliOps.Z, normalizers) +16-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("ZZZZZZ") + pauliop("ZYIIYZ") + pauliop("ZZYIIY") + pauliop("ZYXZXY") + pauliop("ZYZYII") + pauliop("ZZIXXI") + pauliop("ZYYXZX") + pauliop("ZZXYYX") + pauliop("ZIYZYI") + pauliop("ZXXIZI") + pauliop("ZIZIXX") + pauliop("ZXIZIX") + pauliop("ZXYYXZ") + pauliop("ZIXXIZ") + pauliop("ZXZXYY") + pauliop("ZIIYZY") + +``` + +These normalizers are generated from one logical operator and stabilizers. +```jldoctest 1 +julia> map(x -> x .* pauliop("XXXXXX"), stabs) |> Set == Set(norm_x) +true + +julia> map(x -> x .* pauliop("ZZZZZZ"), stabs) |> Set == Set(norm_x) +false + +julia> map(x -> x .* pauliop("ZZZZZZ"), stabs) |> Set == Set(norm_z) +true + +julia> map(x -> x .* pauliop("YYYYYY"), stabs) |> Set == Set(norm_y) +true + +julia> using IterTools + +julia> groupby(x -> x[1], normalizers) .|> Set == Set.([stabs, norm_x, norm_z, norm_y]) +true + +``` + +Define a function to get weight of the operator. +```jldoctest 1 +julia> function weight(x, i = 1) + count(x[i:end] .!= PauliOps.I) + end +weight (generic function with 2 methods) + +julia> weight(pauliop("XIXIXI")) +3 + +julia> weight(pauliop("XIXIXI"), 2) +2 + +``` + +Calculate coefficients of enumerator polynomial. +```jldoctest 1 +julia> using DataStructures + +julia> stabs .|> weight |> counter +Accumulator{Int64, Int64} with 2 entries: + 0 => 1 + 4 => 15 + +julia> function weight(i::Integer) + Base.Fix2(weight, i) + end +weight (generic function with 3 methods) + +julia> normalizers .|> weight(2) |> counter +Accumulator{Int64, Int64} with 4 entries: + 0 => 1 + 4 => 15 + 5 => 18 + 3 => 30 + +julia> [norm_x..., norm_y..., norm_z...] .|> weight(2) |> counter +Accumulator{Int64, Int64} with 2 entries: + 5 => 18 + 3 => 30 + +julia> [norm_x..., norm_y..., norm_z...] .|> weight(2) |> counter |> keys |> minimum +3 + +``` +Code distance of the encoding is the minimum degree of the non-zero term in the normalizer's polynomial($B$) and not in the stabilizer's polynomial($A$). +So the code distance of this encoding is $3$. diff --git a/docs/src/img/bench_ref_optimized.jpg b/docs/src/img/bench_ref_optimized.jpg new file mode 100644 index 0000000..227215a Binary files /dev/null and b/docs/src/img/bench_ref_optimized.jpg differ diff --git a/docs/src/img/bench_ref_vs_aarref.jpg b/docs/src/img/bench_ref_vs_aarref.jpg new file mode 100644 index 0000000..b64926b Binary files /dev/null and b/docs/src/img/bench_ref_vs_aarref.jpg differ diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..b32b942 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,201 @@ +```@meta +CurrentModule = QuantumLegos +``` + +# QuantumLegos + +Documentation for QuantumLegos. + +All contents: + +```@contents +``` + +# Example + +## CheckMatrix and defining Lego +```jldoctest +julia> using QuantumLegos + +julia> stabilizers = pauliop.(["IIXXXX", "IIZZZZ", "ZIZZII", "IZZIZI", "IXXXII", "XIXIXI"]) +6-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IIXXXX") + pauliop("IIZZZZ") + pauliop("ZIZZII") + pauliop("IZZIZI") + pauliop("IXXXII") + pauliop("XIXIXI") + +julia> cmat = checkmatrix(stabilizers) +CheckMatrix with 6 generators, 6 legs: + 0 0 1 1 1 1 | 0 0 0 0 0 0 + 0 0 0 0 0 0 | 0 0 1 1 1 1 + 0 0 0 0 0 0 | 1 0 1 1 0 0 + 0 0 0 0 0 0 | 0 1 1 0 1 0 + 0 1 1 1 0 0 | 0 0 0 0 0 0 + 1 0 1 0 1 0 | 0 0 0 0 0 0 + + +julia> cmat.nlegs +6 + +julia> cmat.ngens +6 + +julia> cmat.cmat +6×12 Matrix{Bool}: + 0 0 1 1 1 1 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 1 1 1 1 + 0 0 0 0 0 0 1 0 1 1 0 0 + 0 0 0 0 0 0 0 1 1 0 1 0 + 0 1 1 1 0 0 0 0 0 0 0 0 + 1 0 1 0 1 0 0 0 0 0 0 0 + +julia> # define lego + +julia> lego = Lego(stabilizers) +Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]) + +julia> lego.stabgens |> checkmatrix +CheckMatrix with 6 generators, 6 legs: + 0 0 1 1 1 1 | 0 0 0 0 0 0 + 0 0 0 0 0 0 | 0 0 1 1 1 1 + 0 0 0 0 0 0 | 1 0 1 1 0 0 + 0 0 0 0 0 0 | 0 1 1 0 1 0 + 0 1 1 1 0 0 | 0 0 0 0 0 0 + 1 0 1 0 1 0 | 0 0 0 0 0 0 + +``` + +- [`pauliop`](@ref) +- [`checkmatrix`](@ref) and [`CheckMatrix`](@ref) +- [`Lego`](@ref) + +## Defining and Updating State +```jldoctest +julia> using QuantumLegos + +julia> stabilizers = pauliop.(["IIXXXX", "IIZZZZ", "ZIZZII", "IZZIZI", "IXXXII", "XIXIXI"]) +6-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IIXXXX") + pauliop("IIZZZZ") + pauliop("ZIZZII") + pauliop("IZZIZI") + pauliop("IXXXII") + pauliop("XIXIXI") + +julia> lego = Lego(stabilizers) +Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]) + +julia> # state with 1 lego, 0 leg + +julia> st = State([lego, ], Tuple{LegoLeg, LegoLeg}[]) +State(Lego[Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")])], Tuple{LegoLeg, LegoLeg}[], CheckMatrix(Bool[0 0 … 0 0; 0 0 … 1 1; … ; 0 1 … 0 0; 1 0 … 0 0], 6, 6)) + +julia> st.cmat.cmat +6×12 Matrix{Bool}: + 0 0 1 1 1 1 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 1 1 1 1 + 0 0 0 0 0 0 1 0 1 1 0 0 + 0 0 0 0 0 0 0 1 1 0 1 0 + 0 1 1 1 0 0 0 0 0 0 0 0 + 1 0 1 0 1 0 0 0 0 0 0 0 + +julia> add_lego!(st, lego) +State(Lego[Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]), Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")])], Tuple{LegoLeg, LegoLeg}[], CheckMatrix(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0], 12, 12)) + +julia> st.cmat.cmat +12×24 Matrix{Bool}: + 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 + 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 + 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 + +julia> # state with 2 legos, 0 leg + +julia> st2 = State([lego, lego], Tuple{LegoLeg, LegoLeg}[]) +State(Lego[Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]), Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")])], Tuple{LegoLeg, LegoLeg}[], CheckMatrix(Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0], 12, 12)) + +julia> st == st2 +true +``` + +## 2 Lego 1 edge state +```jldoctest +julia> using QuantumLegos + +julia> stabilizers = pauliop.(["IIXXXX", "IIZZZZ", "ZIZZII", "IZZIZI", "IXXXII", "XIXIXI"]) +6-element Vector{StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}}: + pauliop("IIXXXX") + pauliop("IIZZZZ") + pauliop("ZIZZII") + pauliop("IZZIZI") + pauliop("IXXXII") + pauliop("XIXIXI") + +julia> lego = Lego(stabilizers) +Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]) + +julia> state = State([lego, lego], edge.([((1, 3), (2, 3))])) +State(Lego[Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")]), Lego{6}(6, StaticArraysCore.SVector{6, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("IIXXXX"), pauliop("IIZZZZ"), pauliop("ZIZZII"), pauliop("IZZIZI"), pauliop("IXXXII"), pauliop("XIXIXI")])], Tuple{LegoLeg, LegoLeg}[(LegoLeg(1, 3), LegoLeg(2, 3))], CheckMatrix(Bool[1 0 … 0 0; 0 1 … 0 0; … ; 0 0 … 1 1; 0 0 … 0 1], 10, 10)) + +julia> state.cmat +CheckMatrix with 10 generators, 10 legs: + 1 0 1 0 1 0 0 0 0 0 | 0 0 0 0 0 0 0 0 0 0 + 0 1 0 1 1 0 0 0 0 0 | 0 0 0 0 0 0 0 0 0 0 + 0 0 1 1 1 0 0 1 1 1 | 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 1 0 1 0 1 | 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 1 1 | 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 | 1 0 0 1 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 | 0 1 1 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 | 0 0 1 1 1 0 0 1 1 1 + 0 0 0 0 0 0 0 0 0 0 | 0 0 0 0 0 1 0 0 1 1 + 0 0 0 0 0 0 0 0 0 0 | 0 0 0 0 0 0 1 1 0 1 + + +julia> pg = state.cmat |> generators |> GeneratedPauliGroup +GeneratedPauliGroup{10}(StaticArraysCore.SVector{10, QuantumLegos.PauliOps.SinglePauliOp}[pauliop("XIXIXIIIII"), pauliop("IXIXXIIIII"), pauliop("IIXXXIIXXX"), pauliop("IIIIIXIXIX"), pauliop("IIIIIIXIXX"), pauliop("ZIIZZIIIII"), pauliop("IZZIZIIIII"), pauliop("IIZZZIIZZZ"), pauliop("IIIIIZIIZZ"), pauliop("IIIIIIZZIZ")], IterTools.Subsets{Vector{StaticArraysCore.SVector{N, QuantumLegos.PauliOps.SinglePauliOp} where N}}(StaticArraysCore.SVector{N, QuantumLegos.PauliOps.SinglePauliOp} where N[pauliop("XIXIXIIIII"), pauliop("IXIXXIIIII"), pauliop("IIXXXIIXXX"), pauliop("IIIIIXIXIX"), pauliop("IIIIIIXIXX"), pauliop("ZIIZZIIIII"), pauliop("IZZIZIIIII"), pauliop("IIZZZIIZZZ"), pauliop("IIIIIZIIZZ"), pauliop("IIIIIIZZIZ")])) + +julia> pauliop("XIIXIXIIXI") in pg +true + +``` + +# Internal(how it works) + +## Notes on Overall flow +Details on [^1] +- state is translated to a single check matrix + - the size is ≤ $N \times 2N$ where $N$ is maximum number of lego logs. +- any contraction can be performed on this single check matrix +- if the check matrix can be represented as direct sum of matrices with $k N$ columns where $k ∈ ℕ$, then they are not contracted + +### Construction of State +Construction of `State` is completed by calling `State` constructor recursively. + +1. Construct `State` without edge. Just adding legos. Checkmatrix is just a direct sum of each lego's checkmatrix +2. Concatenate each edges. During this operation, self tracing of checkmatrix is evaluated. + +Each constructor calls action function (which is a map from `State` to `State`). +Therefore, action functions can be used both for direct construction of `State` and action application to `State` during the game. + +# API + +```@index +``` + +```@autodocs +Modules = [QuantumLegos] +Pages = ["game.jl"] +``` + +[^1]: [C. Cao and B. Lackey, ‘Approximate Bacon-Shor code and holography’, J. High Energ. Phys., vol. 2021, no. 5, p. 127, May 2021, doi: 10.1007/JHEP05(2021)127.](https://doi.org/10.1007/JHEP05(2021)127) + diff --git a/docs/src/pauliops.md b/docs/src/pauliops.md new file mode 100644 index 0000000..fe6e58d --- /dev/null +++ b/docs/src/pauliops.md @@ -0,0 +1,46 @@ +# PauliOps + +submodule + +## Example +```jldoctest +julia> using QuantumLegos + +julia> p = PauliOps.single_pauliop('I') +I::SinglePauliOp = 0 + +julia> typeof(p) +Enum QuantumLegos.PauliOps.SinglePauliOp: +I = 0 +X = 1 +Y = 2 +Z = 3 + +julia> pauliop("IXYZ") +4-element PauliOp: + I::SinglePauliOp = 0 + X::SinglePauliOp = 1 + Y::SinglePauliOp = 2 + Z::SinglePauliOp = 3 + +julia> typeof(ans) +SVector{4, SinglePauliOp} (alias for StaticArraysCore.SArray{Tuple{4}, QuantumLegos.PauliOps.SinglePauliOp, 1, 4}) + +julia> PauliOps.I * PauliOps.X +X::SinglePauliOp = 1 + +julia> PauliOps.X * PauliOps.Z +Y::SinglePauliOp = 2 + +julia> pauliop("IIX") .* pauliop("XIY") +3-element PauliOp: + X::SinglePauliOp = 1 + I::SinglePauliOp = 0 + Z::SinglePauliOp = 3 +``` + +## API +```@autodocs +Modules = [PauliOps] +``` + diff --git a/examples/t6_2legos_notebook.jl b/examples/t6_2legos_notebook.jl new file mode 100644 index 0000000..936b62a --- /dev/null +++ b/examples/t6_2legos_notebook.jl @@ -0,0 +1,47 @@ +### A Pluto.jl notebook ### +# v0.17.7 + +using Markdown +using InteractiveUtils + +# ╔═╡ 3c6bf5ec-909f-11ee-06bd-83347655e198 +begin + import Pkg + Pkg.develop(path = "..") + using QuantumLegos +end + +# ╔═╡ 9088d661-7bf1-43fb-ab88-77fa325a5cf3 +stabilizers = pauliop.(["IIXXXX", "IIZZZZ", "ZIZZII", "IZZIZI", "IXXXII", "XIXIXI"]) + +# ╔═╡ f806287c-592d-476b-a912-205d2031fd93 +lego = Lego(stabilizers) + +# ╔═╡ 93251f25-5829-45a7-8aed-f76c834050a9 +state = State([lego, lego], Tuple{LegoLeg, LegoLeg}[]) + +# ╔═╡ 924588fb-0020-47e6-a918-98084d1fabad +state.cmat + +# ╔═╡ 99f153a1-da44-499a-b8af-e5c484b70597 +QuantumLegos.self_trace!(state.cmat, 3, 9) + +# ╔═╡ 726061b5-0d3a-4bf4-aebd-81a2c0fe7ea1 +state.cmat |> generators + +# ╔═╡ 69a71bfd-81d3-4961-9051-5f19be20f286 +pg = state.cmat |> generators |> GeneratedPauliGroup |> collect + +# ╔═╡ 656d8d7a-0ede-4621-99f0-9f83619c6a73 +pauliop("XIIXIXIIXI") in pg # example on Fig.6 + +# ╔═╡ Cell order: +# ╠═3c6bf5ec-909f-11ee-06bd-83347655e198 +# ╠═9088d661-7bf1-43fb-ab88-77fa325a5cf3 +# ╠═f806287c-592d-476b-a912-205d2031fd93 +# ╠═93251f25-5829-45a7-8aed-f76c834050a9 +# ╠═924588fb-0020-47e6-a918-98084d1fabad +# ╠═99f153a1-da44-499a-b8af-e5c484b70597 +# ╠═726061b5-0d3a-4bf4-aebd-81a2c0fe7ea1 +# ╠═69a71bfd-81d3-4961-9051-5f19be20f286 +# ╠═656d8d7a-0ede-4621-99f0-9f83619c6a73 diff --git a/src/PauliOps/PauliOps.jl b/src/PauliOps/PauliOps.jl new file mode 100644 index 0000000..dce3ceb --- /dev/null +++ b/src/PauliOps/PauliOps.jl @@ -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 diff --git a/src/QuantumLegos.jl b/src/QuantumLegos.jl new file mode 100644 index 0000000..c1ecac0 --- /dev/null +++ b/src/QuantumLegos.jl @@ -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 diff --git a/src/checkmatrix.jl b/src/checkmatrix.jl new file mode 100644 index 0000000..ec5d827 --- /dev/null +++ b/src/checkmatrix.jl @@ -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 diff --git a/src/game.jl b/src/game.jl new file mode 100644 index 0000000..b680193 --- /dev/null +++ b/src/game.jl @@ -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 diff --git a/test/pauliops.jl b/test/pauliops.jl new file mode 100644 index 0000000..455a356 --- /dev/null +++ b/test/pauliops.jl @@ -0,0 +1,26 @@ +@testset "SinglePauliOp" begin + @test PauliOps.single_pauliop('I') == PauliOps.I + @test_throws ArgumentError PauliOps.single_pauliop('a') +end + +@testset "SinglePauliOp Product" begin + using QuantumLegos.PauliOps: I, X, Y, Z + @test I * X == X + @test Z * I == Z + @test X * Y == Z + @test Z * X == Y + @test X * X == I +end + +@testset "PauliOp" begin + @test pauliop("IXYZ") == [PauliOps.I, PauliOps.X, PauliOps.Y, PauliOps.Z] + @test pauliop("IIXXZZ") == + [PauliOps.I, PauliOps.I, PauliOps.X, PauliOps.X, PauliOps.Z, PauliOps.Z] +end + +@testset "weight" begin + p = pauliop("IXYZIXYZ") + @test weight(p) == 6 + @test QuantumLegos.xweight(p) == 4 + @test QuantumLegos.zweight(p) == 4 +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..415ca53 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,405 @@ +# for language server completion +# must be removed for test +true || include("../src/QuantumLegos.jl") + +using QuantumLegos +using Test +using Documenter +using Aqua +using JET + +@testset "QuantumLegos.jl" begin + @testset "PauliOps" begin + include("pauliops.jl") + end + + @testset "Lego" begin + stabilizers = pauliop.(["II", "XX"]) + wrong_stabilizers = pauliop.(["I", "XX"]) # not SVector + @info stabilizers + @test_throws ArgumentError QuantumLegos.Lego(3, stabilizers) + @test_throws MethodError QuantumLegos.Lego(2, wrong_stabilizers) + end + + @testset "CheckMatrix" begin + @test_throws ArgumentError QuantumLegos.CheckMatrix([true false; false true], 2, 2) + @test_throws ArgumentError QuantumLegos.CheckMatrix([true false; false true], 1, 1) + @test QuantumLegos.CheckMatrix([true false; false true]) == + QuantumLegos.CheckMatrix(Bool[1 0; 0 1], 1, 2) + @test_throws ArgumentError QuantumLegos.CheckMatrix([true false true; false true true]) + + @test_throws MethodError QuantumLegos.checkmatrix(pauliop.(["I", "IX"])) + @test_throws ArgumentError QuantumLegos.checkmatrix(PauliOps.PauliOp[]) + @test QuantumLegos.checkmatrix(pauliop.(["IIXX", "IZZI"])) == + QuantumLegos.CheckMatrix(Bool[0 0 1 1 0 0 0 0; 0 0 0 0 0 1 1 0]) + @test QuantumLegos.checkmatrix(pauliop.(["YIXX", "IZYI"])) == + QuantumLegos.CheckMatrix(Bool[1 0 1 1 1 0 0 0; 0 0 1 0 0 1 1 0]) + + let + gens = pauliop.(["IIXX", "IZZI"]) + cmat_1 = QuantumLegos.checkmatrix(gens) + @test QuantumLegos.xpart(cmat_1) == Bool[0 0 1 1; 0 0 0 0] + @test QuantumLegos.zpart(cmat_1) == Bool[0 0 0 0; 0 1 1 0] + @test generators(cmat_1) == gens + end + let + gens = pauliop.(["YIXX", "IZYI"]) + cmat_2 = QuantumLegos.checkmatrix(gens) + @test QuantumLegos.xpart(cmat_2) == Bool[1 0 1 1; 0 0 1 0] + @test QuantumLegos.zpart(cmat_2) == Bool[1 0 0 0; 0 1 1 0] + @test generators(cmat_2) == gens + end + + @testset "eliminate_column!" begin + let + cmat = QuantumLegos.CheckMatrix(Bool[ + 1 0 0 1 + 1 0 1 1 + ]) + @test QuantumLegos.eliminate_column!(cmat, 1, Int64[]) == 1 + @test cmat.cmat == Bool[ + 1 0 0 1 + 0 0 1 0 + ] + end + let + cmat = QuantumLegos.CheckMatrix(Bool[ + 1 0 0 1 + 1 0 1 1 + ]) + @test QuantumLegos.eliminate_column!(cmat, 1, Int64[1]) == 2 + @test cmat.cmat == Bool[ + 0 0 1 0 + 1 0 1 1 + ] + end + let + cmat = QuantumLegos.CheckMatrix( + Bool[ + 1 1 0 1 0 0 1 1 + 0 0 1 1 0 1 0 1 + 0 1 1 0 1 0 1 1 + 0 1 0 1 0 1 0 1 + ], + ) + @test QuantumLegos.eliminate_column!(cmat, 2, [1]) == 3 + @test cmat.cmat == Bool[ + 1 0 1 1 1 0 0 0 + 0 0 1 1 0 1 0 1 + 0 1 1 0 1 0 1 1 + 0 0 1 1 1 1 1 0 + ] + end + let + cmat = QuantumLegos.CheckMatrix( + Bool[ + 1 1 0 1 0 0 1 1 + 0 0 1 1 0 1 0 1 + 0 1 1 0 1 0 1 1 + 0 1 0 1 0 1 0 1 + ], + ) + @test QuantumLegos.eliminate_column!(cmat, 2, [1, nothing]) == 3 + @test cmat.cmat == Bool[ + 1 0 1 1 1 0 0 0 + 0 0 1 1 0 1 0 1 + 0 1 1 0 1 0 1 1 + 0 0 1 1 1 1 1 0 + ] + end + let + cmat = QuantumLegos.CheckMatrix( + Bool[ + 1 0 0 0 1 1 0 1 + 0 1 0 1 0 1 1 0 + 0 1 0 0 0 1 0 0 + 1 1 0 1 0 0 0 1 + ], + ) + @test QuantumLegos.eliminate_column!(cmat, 3, Int[]) === nothing + end + let + # random test (robustness) + function random_test()::Bool + test_mat = rand(Bool, (8, 16)) + cmat = QuantumLegos.CheckMatrix(copy(test_mat)) + keep_indices = (1:8)[rand(1:8, rand(0:8))] + keep_index = rand(1:16) + try + QuantumLegos.eliminate_column!(cmat, keep_index, keep_indices) + return true + catch e + # unexpected error + @error "Error on eliminate_column! with keep_index: $(keep_index), avoid: $(keep_indices), matrix:" cmat + @error "error: " e + return false + end + end + + for i in 1:100 + # @info i + @test random_test() + end + end + @testset "ArgumentError" begin + cmat = CheckMatrix(rand(Bool, (6, 12))) + @test_throws ArgumentError QuantumLegos.eliminate_column!(cmat, 3, [1, 7]) + end + end + @testset "swap_row!" begin + let + m = rand(Bool, (10, 10)) + m_copy = copy(m) + QuantumLegos.swap_row!(m, 3, 7) + QuantumLegos.swap_row!(m, 3, 7) + @test m == m_copy + end + let + m = rand(Bool, (10, 10)) + m = BitMatrix(m) + m_copy = copy(m) + QuantumLegos.swap_row!(m, 4, 2) + QuantumLegos.swap_row!(m, 4, 2) + @test m == m_copy + end + end + @testset "align_row!" begin + mat = Bool[ + 1 0 1 1 0 0 1 0 + 0 1 0 1 0 1 0 0 + 1 1 1 0 0 1 0 0 + 1 0 0 1 0 1 0 0 + ] + mat_1 = copy(mat) + @test 1 == QuantumLegos.align_row!(mat_1, 3, Int64[]) + @test mat_1 == Bool[ + 1 1 1 0 0 1 0 0 + 0 1 0 1 0 1 0 0 + 1 0 1 1 0 0 1 0 + 1 0 0 1 0 1 0 0 + ] + mat_1 = copy(mat) + @test 1 == QuantumLegos.align_row!(mat_1, 3, [nothing]) + @test mat_1 == Bool[ + 1 1 1 0 0 1 0 0 + 0 1 0 1 0 1 0 0 + 1 0 1 1 0 0 1 0 + 1 0 0 1 0 1 0 0 + ] + mat_2 = copy(mat) + @test 2 == QuantumLegos.align_row!(mat_2, 3, [1]) + @test mat_2 == Bool[ + 1 0 1 1 0 0 1 0 + 1 1 1 0 0 1 0 0 + 0 1 0 1 0 1 0 0 + 1 0 0 1 0 1 0 0 + ] + mat_3 = copy(mat) + @test 3 == QuantumLegos.align_row!(mat_3, 3, [1, nothing, 3]) + @test mat_3 == Bool[ + 1 0 1 1 0 0 1 0 + 0 1 0 1 0 1 0 0 + 1 1 1 0 0 1 0 0 # 3 + 1 0 0 1 0 1 0 0 # 4 + ] + mat_3_2 = copy(mat_3) + @test 4 == QuantumLegos.align_row!(mat_3, 4, [1, nothing, 3]) + @test mat_3 == mat_3_2 + @test nothing === QuantumLegos.align_row!(mat_3, nothing, [1, 2]) + @test nothing === QuantumLegos.align_row!(mat_3, nothing, [1, 2, nothing]) + @test nothing === QuantumLegos.align_row!(mat_3, nothing, [nothing]) + end + @testset "ref!" begin + @testset "compare with AbstractAlgebra" begin + # test using existing package + using AbstractAlgebra + F₂ = GF(2) # finite field + + """ + Compare self-implemented `ref!` with AbstractAlgebra's `rref!` or `rank`. + """ + function random_test(size::Tuple{T, T}) where {T <: Integer} + mat = rand(Bool, size) + cmat = CheckMatrix(mat) + S = matrix_space(F₂, size...) + cmat_aa = S(F₂.(mat)) + # r_aa, A_aa = AbstractAlgebra.rref(cmat_aa) + # cmat_aa = A_aa.entries .|> ==(1) |> Matrix{Bool} |> CheckMatrix + r_aa = rank(cmat_aa) + r = QuantumLegos.ref!(cmat) + # @info "compare cmat" cmat cmat_aa + @test r == r_aa + end + for _ in 1:10 + random_test((2, 4)) + random_test((4, 8)) + random_test((100, 200)) + end + end + @testset "manual sample" begin + import LinearAlgebra + mat = + Bool[1 0 0 1 1 0 0 0; 0 1 0 0 1 1 0 1; 0 0 1 1 0 1 0 1; 0 1 0 0 1 0 0 1] + cmat = CheckMatrix(mat) + @test QuantumLegos.ref!(cmat) == LinearAlgebra.rank(mat) + @test cmat == CheckMatrix( + Bool[ + 1 0 0 1 1 0 0 0 + 0 1 0 0 1 1 0 1 + 0 0 1 1 0 1 0 1 + 0 0 0 0 0 1 0 0 + ], + 4, + 4, + ) + let + mat = copy(mat) + dependent_row = reduce(.⊻, eachrow(mat)[[1, 2]]) + mat[4, :] = dependent_row + cmat = CheckMatrix(mat) + @test QuantumLegos.ref!(cmat) == LinearAlgebra.rank(mat) == 3 + end + let + mat = copy(mat) + dependent_row = reduce(.⊻, eachrow(mat)[[1, 2]]) + mat = vcat(mat, dependent_row') + dependent_row = reduce(.⊻, eachrow(mat)[[1, 2, 3]]) + mat = vcat(mat, dependent_row') + dependent_row = reduce(.⊻, eachrow(mat)[[2, 4]]) + mat = vcat(mat, dependent_row') + cmat = CheckMatrix(mat) + @test QuantumLegos.ref!(cmat) == rank(mat) == 3 + end + end + @testset "generated group is invariant under ref!" begin + for _ in 1:10 + cmat = CheckMatrix(rand(Bool, (8, 12))) + before = cmat |> generators |> GeneratedPauliGroup |> Set + QuantumLegos.ref!(cmat) + @test cmat |> generators |> GeneratedPauliGroup |> Set |> ==(before) + end + end + end + @testset "eliminate_dependent_row!" begin + @testset "trivial" begin + cmat = CheckMatrix( + Bool[ + 1 0 0 1 0 1 0 1 + 0 1 1 0 0 1 1 0 + 0 0 0 1 1 0 1 0 + 0 0 0 0 0 0 0 0 + ], + ) + @test QuantumLegos.eliminate_dependent_row!(cmat) == CheckMatrix( + Bool[1 0 0 1 0 1 0 1; 0 1 1 0 0 1 1 0; 0 0 0 1 1 0 1 0], + 4, + 3, + ) + end + @testset "less trivial" begin + cmat = CheckMatrix( + Bool[ + 1 0 0 1 0 1 0 1 + 0 1 1 0 0 1 1 0 + 0 0 0 1 1 0 1 0 + 1 1 1 1 0 0 1 1 + ], + ) + @test QuantumLegos.eliminate_dependent_row!(cmat) == CheckMatrix( + Bool[1 0 0 1 0 1 0 1; 0 1 1 0 0 1 1 0; 0 0 0 1 1 0 1 0], + 4, + 3, + ) + end + # TODO: add more? + end + @testset "self_trace!" begin + cmat = CheckMatrix(rand(Bool, (8, 16))) + @test_throws ArgumentError QuantumLegos.self_trace!(cmat, 16, 17) + @test_throws ArgumentError QuantumLegos.self_trace!(cmat, 18, 17) + end + end + + @testset "State" begin + @testset "LegoLeg" begin + @test LegoLeg(0, 1) < LegoLeg(1, 0) + @test LegoLeg(2, 1) > LegoLeg(1, 0) + @test LegoLeg(1, 0) < LegoLeg(1, 1) + @test LegoLeg(1, 1) == LegoLeg(1, 1) + + @test sort(LegoLeg.([(1, 2), (2, 3), (0, 2), (0, 1)])) == + LegoLeg[LegoLeg(0, 1), LegoLeg(0, 2), LegoLeg(1, 2), LegoLeg(2, 3)] + end + @testset "edge" begin + @test edge(1, 2, 3, 4) == (LegoLeg(1, 2), LegoLeg(3, 4)) + @test edge((1, 2, 3, 4)) == (LegoLeg(1, 2), LegoLeg(3, 4)) + @test edge(((1, 2), (3, 4))) == (LegoLeg(1, 2), LegoLeg(3, 4)) + end + @testset "0 lego, 0 leg" begin + @test_throws ArgumentError QuantumLegos.State(Lego{6}[], Tuple{LegoLeg, LegoLeg}[]) + end + @testset "1 lego, 0 leg" begin + stabgens = pauliop.(["IIXX", "XXII", "IZZI", "ZIIZ"]) + lego = QuantumLegos.Lego(stabgens) + state = QuantumLegos.State([lego], Tuple{LegoLeg, LegoLeg}[]) + cmat = QuantumLegos.checkmatrix(stabgens) + @test state == QuantumLegos.State([lego], Tuple{LegoLeg, LegoLeg}[]) + @test state.legos == [lego] + @test state.edges == Tuple{LegoLeg, LegoLeg}[] + @test state.cmat == cmat + end + @testset "2+ legos, 0 leg" begin + stabgens = pauliop.(["IIXX", "XXII", "IZZI", "ZIIZ"]) + lego = QuantumLegos.Lego(stabgens) + state_1 = QuantumLegos.State([lego], Tuple{LegoLeg, LegoLeg}[]) + add_lego!(state_1, lego) + @test state_1.legos == [lego, lego] + @test state_1.edges == Tuple{LegoLeg, LegoLeg}[] + @test state_1.cmat.ngens == 8 + @test state_1.cmat.nlegs == 8 + state_2 = QuantumLegos.State([lego, lego], Tuple{LegoLeg, LegoLeg}[]) + @test state_1 == state_2 + @test all(state_2.cmat.cmat[5:8, 1:4] .== false) # block diagonal + end + @testset "2+ legos, 1+ legs" begin + stabgens = pauliop.(["IIXX", "XXII", "IZZI", "ZIIZ"]) + lego = QuantumLegos.Lego(stabgens) + @test_throws ArgumentError State([lego, lego], edge.([(1, 1, 1, 1)])) + end + @testset "is_connected_to_firstlego" begin + stabilizers = pauliop.(["IIXXXX", "IIZZZZ", "ZIZIZI", "IZIZIZ", "IXIIXX", "XIXXII"]) + lego = Lego(stabilizers) + let + state = + State(fill(lego, 6), edge.([(3, 2, 1, 5), (5, 2, 3, 3), (4, 1, 5, 1)])) + @test QuantumLegos.is_connected_to_firstlego(state) == + BitVector([1, 0, 1, 1, 1, 0]) + end + let + state = State( + fill(lego, 7), + edge.([(3, 2, 1, 5), (5, 2, 3, 3), (4, 1, 5, 1), (2, 2, 3, 1)]), + ) + @test QuantumLegos.is_connected_to_firstlego(state) == + BitVector([1, 1, 1, 1, 1, 0, 0]) + end + end + end + + @testset "Doctest" begin + DocMeta.setdocmeta!( + QuantumLegos, + :DocTestSetup, + :(using QuantumLegos; ENV["JULIA_DEBUG"] = ""); + recursive = true, + ) + doctest(QuantumLegos) + end + + # @testset "Code quality (Aqua.jl)" begin + # Aqua.test_all(QuantumLegos) + # end + # @testset "Code linting (JET.jl)" begin + # JET.test_package(QuantumLegos; target_defined_modules = true) + # end +end