diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml deleted file mode 100644 index 4c49a86..0000000 --- a/.JuliaFormatter.toml +++ /dev/null @@ -1,3 +0,0 @@ -# See https://domluna.github.io/JuliaFormatter.jl/stable/ for a list of options -style = "blue" -indent = 2 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 456fa05..0614de9 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -2,7 +2,7 @@ name: "CompatHelper" on: schedule: - - cron: 0 0 * * * + - cron: '0 0 * * *' workflow_dispatch: permissions: contents: write diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 3f78afc..1525861 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,11 +1,14 @@ name: "Format Check" on: - push: - branches: - - 'main' - tags: '*' - pull_request: + pull_request_target: + paths: ['**/*.jl'] + types: [opened, synchronize, reopened, ready_for_review] + +permissions: + contents: read + actions: write + pull-requests: write jobs: format-check: diff --git a/.gitignore b/.gitignore index 10593a9..7085ca8 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,10 @@ .vscode/ Manifest.toml benchmark/*.json +dev/ +docs/LocalPreferences.toml docs/Manifest.toml docs/build/ docs/src/index.md +examples/LocalPreferences.toml +test/LocalPreferences.toml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 88bc8b4..3fc4743 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,5 +1,5 @@ ci: - skip: [julia-formatter] + skip: [runic] repos: - repo: https://github.com/pre-commit/pre-commit-hooks @@ -11,7 +11,7 @@ repos: - id: end-of-file-fixer exclude_types: [markdown] # incompatible with Literate.jl -- repo: "https://github.com/domluna/JuliaFormatter.jl" - rev: v2.1.6 +- repo: https://github.com/fredrikekre/runic-pre-commit + rev: v2.0.1 hooks: - - id: "julia-formatter" + - id: runic diff --git a/Project.toml b/Project.toml index fb9f230..2ae50cc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumOperatorDefinitions" uuid = "826dd319-6fd5-459a-a990-3a4f214664bf" authors = ["ITensor developers and contributors"] -version = "0.2.8" +version = "0.2.9" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index fc418cf..5321428 100644 --- a/README.md +++ b/README.md @@ -82,7 +82,7 @@ using Test: @test @test Matrix(OpName("Y")) == [0 -im; im 0] @test Matrix(OpName("Z")) == [1 0; 0 -1] -@test Matrix(OpName("Rx"; θ=π / 3)) ≈ [sin(π / 3) -cos(π / 3)*im; -cos(π / 3)*im sin(π / 3)] +@test Matrix(OpName("Rx"; θ = π / 3)) ≈ [sin(π / 3) -cos(π / 3) * im; -cos(π / 3) * im sin(π / 3)] ```` --- diff --git a/docs/make.jl b/docs/make.jl index bbcaa3e..1e16d9a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,28 +2,28 @@ using QuantumOperatorDefinitions: QuantumOperatorDefinitions using Documenter: Documenter, DocMeta, deploydocs, makedocs DocMeta.setdocmeta!( - QuantumOperatorDefinitions, - :DocTestSetup, - :(using QuantumOperatorDefinitions); - recursive=true, + QuantumOperatorDefinitions, + :DocTestSetup, + :(using QuantumOperatorDefinitions); + recursive = true, ) include("make_index.jl") makedocs(; - modules=[QuantumOperatorDefinitions], - authors="ITensor developers and contributors", - sitename="QuantumOperatorDefinitions.jl", - format=Documenter.HTML(; - canonical="https://itensor.github.io/QuantumOperatorDefinitions.jl", - edit_link="main", - assets=["assets/favicon.ico", "assets/extras.css"], - ), - pages=["Home" => "index.md", "Reference" => "reference.md"], + modules = [QuantumOperatorDefinitions], + authors = "ITensor developers and contributors", + sitename = "QuantumOperatorDefinitions.jl", + format = Documenter.HTML(; + canonical = "https://itensor.github.io/QuantumOperatorDefinitions.jl", + edit_link = "main", + assets = ["assets/favicon.ico", "assets/extras.css"], + ), + pages = ["Home" => "index.md", "Reference" => "reference.md"], ) deploydocs(; - repo="github.com/ITensor/QuantumOperatorDefinitions.jl", - devbranch="main", - push_preview=true, + repo = "github.com/ITensor/QuantumOperatorDefinitions.jl", + devbranch = "main", + push_preview = true, ) diff --git a/docs/make_index.jl b/docs/make_index.jl index 1db1dd1..02d04f2 100644 --- a/docs/make_index.jl +++ b/docs/make_index.jl @@ -2,20 +2,20 @@ using Literate: Literate using QuantumOperatorDefinitions: QuantumOperatorDefinitions function ccq_logo(content) - include_ccq_logo = """ + include_ccq_logo = """ ```@raw html Flatiron Center for Computational Quantum Physics logo. Flatiron Center for Computational Quantum Physics logo. ``` """ - content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) - return content + content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) + return content end Literate.markdown( - joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), - joinpath(pkgdir(QuantumOperatorDefinitions), "docs", "src"); - flavor=Literate.DocumenterFlavor(), - name="index", - postprocess=ccq_logo, + joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), + joinpath(pkgdir(QuantumOperatorDefinitions), "docs", "src"); + flavor = Literate.DocumenterFlavor(), + name = "index", + postprocess = ccq_logo, ) diff --git a/docs/make_readme.jl b/docs/make_readme.jl index 62c6441..d15914f 100644 --- a/docs/make_readme.jl +++ b/docs/make_readme.jl @@ -2,20 +2,20 @@ using Literate: Literate using QuantumOperatorDefinitions: QuantumOperatorDefinitions function ccq_logo(content) - include_ccq_logo = """ + include_ccq_logo = """ Flatiron Center for Computational Quantum Physics logo. """ - content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) - return content + content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) + return content end Literate.markdown( - joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), - joinpath(pkgdir(QuantumOperatorDefinitions)); - flavor=Literate.CommonMarkFlavor(), - name="README", - postprocess=ccq_logo, + joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), + joinpath(pkgdir(QuantumOperatorDefinitions)); + flavor = Literate.CommonMarkFlavor(), + name = "README", + postprocess = ccq_logo, ) diff --git a/examples/README.jl b/examples/README.jl index 1399492..f49ca79 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -1,5 +1,5 @@ # # QuantumOperatorDefinitions.jl -# +# # [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/QuantumOperatorDefinitions.jl/stable/) # [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/QuantumOperatorDefinitions.jl/dev/) # [![Build Status](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml?query=branch%3Amain) @@ -83,4 +83,4 @@ using Test: @test @test Matrix(OpName("Y")) == [0 -im; im 0] @test Matrix(OpName("Z")) == [1 0; 0 -1] -@test Matrix(OpName("Rx"; θ=π / 3)) ≈ [sin(π / 3) -cos(π / 3)*im; -cos(π / 3)*im sin(π / 3)] +@test Matrix(OpName("Rx"; θ = π / 3)) ≈ [sin(π / 3) -cos(π / 3) * im; -cos(π / 3) * im sin(π / 3)] diff --git a/ext/QuantumOperatorDefinitionsGradedArraysExt/QuantumOperatorDefinitionsGradedArraysExt.jl b/ext/QuantumOperatorDefinitionsGradedArraysExt/QuantumOperatorDefinitionsGradedArraysExt.jl index 28f7927..c913b71 100644 --- a/ext/QuantumOperatorDefinitionsGradedArraysExt/QuantumOperatorDefinitionsGradedArraysExt.jl +++ b/ext/QuantumOperatorDefinitionsGradedArraysExt/QuantumOperatorDefinitionsGradedArraysExt.jl @@ -2,97 +2,103 @@ module QuantumOperatorDefinitionsGradedArraysExt using BlockArrays: blocklasts, blocklength, blocklengths using GradedArrays: - AbstractGradedUnitRange, GradedOneTo, SectorProduct, U1, Z, ×, dual, gradedrange, sectors + AbstractGradedUnitRange, GradedOneTo, SectorProduct, U1, Z, ×, dual, gradedrange, sectors using QuantumOperatorDefinitions: - QuantumOperatorDefinitions, - @GradingType_str, - @SiteType_str, - GradingType, - OpName, - SiteType, - name + QuantumOperatorDefinitions, + @GradingType_str, + @SiteType_str, + GradingType, + OpName, + SiteType, + name function Base.axes(::OpName, domain::Tuple{Vararg{AbstractGradedUnitRange}}) - return (domain..., dual.(domain)...) + return (domain..., dual.(domain)...) end sortedunion(a, b) = sort(union(a, b)) function QuantumOperatorDefinitions.combine_axes(a1::GradedOneTo, a2::GradedOneTo) - blocklength(a1) == blocklength(a2) || - throw(ArgumentError("Axes must have the same number of blocks.")) - nblocks = blocklength(a1) - return gradedrange( - map(Base.OneTo(nblocks)) do i - l1 = blocklengths(a1)[i] - l2 = blocklengths(a2)[i] - l1 == l2 || throw(ArgumentError("Blocks must have the same length.")) - return sectors(a1)[i] × sectors(a2)[i] => l1 - end, - ) + blocklength(a1) == blocklength(a2) || + throw(ArgumentError("Axes must have the same number of blocks.")) + nblocks = blocklength(a1) + return gradedrange( + map(Base.OneTo(nblocks)) do i + l1 = blocklengths(a1)[i] + l2 = blocklengths(a2)[i] + l1 == l2 || throw(ArgumentError("Blocks must have the same length.")) + return sectors(a1)[i] × sectors(a2)[i] => l1 + end, + ) end QuantumOperatorDefinitions.combine_axes(a::GradedOneTo, b::Base.OneTo) = a QuantumOperatorDefinitions.combine_axes(a::Base.OneTo, b::GradedOneTo) = b function Base.AbstractUnitRange(::GradingType"N", t::SiteType) - return gradedrange(map(i -> SectorProduct((; N=U1(i - 1))) => 1, 1:length(t))) + return gradedrange(map(i -> SectorProduct((; N = U1(i - 1))) => 1, 1:length(t))) end function Base.AbstractUnitRange(::GradingType"Sz", t::SiteType) - return gradedrange(map(i -> SectorProduct((; Sz=U1(i - 1))) => 1, 1:length(t))) + return gradedrange(map(i -> SectorProduct((; Sz = U1(i - 1))) => 1, 1:length(t))) end function Base.AbstractUnitRange(::GradingType"Sz↑", t::SiteType) - return AbstractUnitRange(GradingType"Sz"(), t) + return AbstractUnitRange(GradingType"Sz"(), t) end function Base.AbstractUnitRange(::GradingType"Sz↓", t::SiteType) - return gradedrange(map(i -> SectorProduct((; Sz=U1(-(i - 1)))) => 1, 1:length(t))) + return gradedrange(map(i -> SectorProduct((; Sz = U1(-(i - 1)))) => 1, 1:length(t))) end function sector(gradingtype::GradingType, sec) - sectorname = Symbol(get(gradingtype, :name, name(gradingtype))) - return SectorProduct(NamedTuple{(sectorname,)}((sec,))) + sectorname = Symbol(get(gradingtype, :name, name(gradingtype))) + return SectorProduct(NamedTuple{(sectorname,)}((sec,))) end function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Fermion") - return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) + return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) end # TODO: Write in terms of `GradingType"Nf"` definition. function Base.AbstractUnitRange(s::GradingType"NfParity", t::SiteType"Fermion") - return gradedrange([sector(s, Z{2}(0)) => 1, sector(s, Z{2}(1)) => 1]) + return gradedrange([sector(s, Z{2}(0)) => 1, sector(s, Z{2}(1)) => 1]) end function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Fermion") - return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) + return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) end function Base.AbstractUnitRange(s::GradingType"Sz↑", t::SiteType"Fermion") - return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) + return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1]) end function Base.AbstractUnitRange(s::GradingType"Sz↓", t::SiteType"Fermion") - return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(-1)) => 1]) + return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(-1)) => 1]) end # TODO: Write in terms of `SiteType"Fermion"` definitions. function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Electron") - return gradedrange([ - sector(s, U1(0)) => 1, - sector(s, U1(1)) => 1, - sector(s, U1(1)) => 1, - sector(s, U1(2)) => 1, - ]) + return gradedrange( + [ + sector(s, U1(0)) => 1, + sector(s, U1(1)) => 1, + sector(s, U1(1)) => 1, + sector(s, U1(2)) => 1, + ] + ) end # TODO: Write in terms of `GradingType"Nf"` definition. function Base.AbstractUnitRange(s::GradingType"NfParity", t::SiteType"Electron") - return gradedrange([ - sector(s, Z{2}(0)) => 1, - sector(s, Z{2}(1)) => 1, - sector(s, Z{2}(1)) => 1, - sector(s, Z{2}(0)) => 1, - ]) + return gradedrange( + [ + sector(s, Z{2}(0)) => 1, + sector(s, Z{2}(1)) => 1, + sector(s, Z{2}(1)) => 1, + sector(s, Z{2}(0)) => 1, + ] + ) end function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Electron") - return gradedrange([ - sector(s, U1(0)) => 1, - sector(s, U1(1)) => 1, - sector(s, U1(-1)) => 1, - sector(s, U1(0)) => 1, - ]) + return gradedrange( + [ + sector(s, U1(0)) => 1, + sector(s, U1(1)) => 1, + sector(s, U1(-1)) => 1, + sector(s, U1(0)) => 1, + ] + ) end end diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl index 59c573b..2f79398 100644 --- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl +++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl @@ -4,61 +4,61 @@ using ITensorBase: ITensorBase, ITensor, Index, gettag, prime, settag using GradedArrays: dual using NamedDimsArrays: dename using QuantumOperatorDefinitions: - QuantumOperatorDefinitions, - @OpName_str, - OpName, - SiteType, - StateName, - has_fermion_string, - name + QuantumOperatorDefinitions, + @OpName_str, + OpName, + SiteType, + StateName, + has_fermion_string, + name function QuantumOperatorDefinitions.SiteType(r::Index) - # We pass the axis of the (unnamed) Index because - # the Index may have originated from a slice, in which - # case the start may not be 1 (for NonContiguousIndex, - # which we need to add support for, it may not even - # be a unit range). - return SiteType( - gettag(r, "sitetype", "Qudit"); dim=Int.(length(r)), range=only(axes(dename(r))) - ) + # We pass the axis of the (unnamed) Index because + # the Index may have originated from a slice, in which + # case the start may not be 1 (for NonContiguousIndex, + # which we need to add support for, it may not even + # be a unit range). + return SiteType( + gettag(r, "sitetype", "Qudit"); dim = Int.(length(r)), range = only(axes(dename(r))) + ) end function (rangetype::Type{<:Index})(t::SiteType) - i = rangetype(AbstractUnitRange(t)) - i = settag(i, "sitetype", String(name(t))) - if haskey(t, :site) - i = settag(i, "site", string(t.site)) - end - return i + i = rangetype(AbstractUnitRange(t)) + i = settag(i, "sitetype", String(name(t))) + if haskey(t, :site) + i = settag(i, "site", string(t.site)) + end + return i end # TODO: Define in terms of `OpName` directly, and define a generic # forwarding method `has_fermion_string(n::String, t) = has_fermion_string(OpName(n), t)`. function QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index) - return has_fermion_string(OpName(n), SiteType(r)) + return has_fermion_string(OpName(n), SiteType(r)) end function Base.axes(::OpName, domain::Tuple{Vararg{Index}}) - return (prime.(domain)..., dual.(domain)...) + return (prime.(domain)..., dual.(domain)...) end ## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{Index}}) ## return (prime.(reverse(domain))..., dag.(domain)...) ## end # Fix ambiguity error with generic `AbstractArray` version. -function ITensorBase.ITensor(n::Union{OpName,StateName}, domain::Index...) - return ITensor(n, domain) +function ITensorBase.ITensor(n::Union{OpName, StateName}, domain::Index...) + return ITensor(n, domain) end # Fix ambiguity error with generic `AbstractArray` version. -function ITensorBase.ITensor(n::Union{OpName,StateName}, domain::Tuple{Vararg{Index}}) - return ITensor(AbstractArray(n, domain), axes(n, domain)) +function ITensorBase.ITensor(n::Union{OpName, StateName}, domain::Tuple{Vararg{Index}}) + return ITensor(AbstractArray(n, domain), axes(n, domain)) end function (arrtype::Type{<:AbstractArray})( - n::Union{OpName,StateName}, domain::Tuple{Vararg{Index}} -) - # Convert to `SiteType` in case the Index specifies a `"sitetype"` tag. - # TODO: Try to build this into the generic codepath. - return ITensor(arrtype(n, SiteType.(domain)), axes(n, domain)) + n::Union{OpName, StateName}, domain::Tuple{Vararg{Index}} + ) + # Convert to `SiteType` in case the Index specifies a `"sitetype"` tag. + # TODO: Try to build this into the generic codepath. + return ITensor(arrtype(n, SiteType.(domain)), axes(n, domain)) end end diff --git a/src/has_fermion_string.jl b/src/has_fermion_string.jl index 614f940..d86af97 100644 --- a/src/has_fermion_string.jl +++ b/src/has_fermion_string.jl @@ -1,7 +1,7 @@ function has_fermion_string(n::OpName, t::SiteType) - n′ = alias(n) - if n == n′ - return false - end - return has_fermion_string(n′, t) + n′ = alias(n) + if n == n′ + return false + end + return has_fermion_string(n′, t) end diff --git a/src/op.jl b/src/op.jl index ea0257d..9f92476 100644 --- a/src/op.jl +++ b/src/op.jl @@ -1,10 +1,10 @@ using LinearAlgebra: diag, qr -struct OpName{Name,Params} - params::Params - function OpName{N}(params::NamedTuple) where {N} - return new{N,typeof(params)}(params) - end +struct OpName{Name, Params} + params::Params + function OpName{N}(params::NamedTuple) where {N} + return new{N, typeof(params)}(params) + end end name(::OpName{Name}) where {Name} = Name params(n::OpName) = getfield(n, :params) @@ -19,7 +19,7 @@ OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...)) # opexpr("Ry{θ=π/2}") == OpName("Ry"; θ=π/2) # ``` function opexpr(n::String; kwargs...) - return state_or_op_expr(OpName, n; kwargs...) + return state_or_op_expr(OpName, n; kwargs...) end # TODO: Should this parse the string? @@ -27,7 +27,7 @@ OpName(s::AbstractString; kwargs...) = OpName{Symbol(s)}(; kwargs...) OpName(s::Symbol; kwargs...) = OpName{s}(; kwargs...) # TODO: Should this parse the string? macro OpName_str(s) - return :(OpName{$(Expr(:quote, Symbol(s)))}) + return :(OpName{$(Expr(:quote, Symbol(s)))}) end # This version parses. Disabled for now until @@ -43,69 +43,71 @@ end # end function op_alias_expr(name1, name2, pars...) - return :(function alias(n::OpName{Symbol($name1)}) - return OpName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) - end) + return :( + function alias(n::OpName{Symbol($name1)}) + return OpName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) + end + ) end macro op_alias(name1, name2, pars...) - return op_alias_expr(name1, name2, pars...) + return op_alias_expr(name1, name2, pars...) end # Generic to `StateName` or `OpName`. -const StateOrOpName = Union{StateName,OpName} +const StateOrOpName = Union{StateName, OpName} alias(n::StateOrOpName) = n function (arrtype::Type{<:AbstractArray})( - n::StateOrOpName, domain::Union{Integer,AbstractUnitRange}... -) - return arrtype(n, domain) + n::StateOrOpName, domain::Union{Integer, AbstractUnitRange}... + ) + return arrtype(n, domain) end function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{Integer}}) - return arrtype(n, Base.oneto.(domain)) + return arrtype(n, Base.oneto.(domain)) end (arrtype::Type{<:AbstractArray})(n::StateOrOpName, ts::SiteType...) = arrtype(n, ts) function (n::StateOrOpName)(domain...) - # TODO: Try one alias at a time? - # TODO: First call `alias(n, domain...)` - # to allow for aliases specific to certain - # SiteTypes? - n′ = alias(n) - domain′ = alias.(domain) - if n′ == n && domain′ == domain - error("Not implemented.") - end - return n′(domain′...) + # TODO: Try one alias at a time? + # TODO: First call `alias(n, domain...)` + # to allow for aliases specific to certain + # SiteTypes? + n′ = alias(n) + domain′ = alias.(domain) + if n′ == n && domain′ == domain + error("Not implemented.") + end + return n′(domain′...) end # TODO: Decide on this. function (n::StateOrOpName)() - return n(ntuple(Returns(default_sitetype()), nsites(n))...) + return n(ntuple(Returns(default_sitetype()), nsites(n))...) end function (arrtype::Type{<:AbstractArray})(n::StateOrOpName) - return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n))) + return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n))) end # `Int(ndims // 2)`, i.e. `ndims_domain`/`ndims_codomain`. function nsites(n::StateOrOpName) - n′ = alias(n) - if n′ == n - # Default value, assume 1-site operator/state. - return 1 - end - return nsites(n′) + n′ = alias(n) + if n′ == n + # Default value, assume 1-site operator/state. + return 1 + end + return nsites(n′) end # TODO: This does some unwanted conversions, like turning # `Diagonal` dense. function array(a::AbstractArray, ax::Tuple{Vararg{AbstractUnitRange}}) - return a[ax...] + return a[ax...] end function Base.axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}}) - return (domain..., domain...) + return (domain..., domain...) end function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{Integer}}) - return axes(n, Base.OneTo.(domain)) + return axes(n, Base.OneTo.(domain)) end function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{SiteType}}) - return axes(n, AbstractUnitRange.(domain)) + return axes(n, AbstractUnitRange.(domain)) end ## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{AbstractUnitRange}}) @@ -113,181 +115,181 @@ end ## end function reversed_sites(n::StateOrOpName, domain) - return reverse_sites(n, reshape(n(domain...), length.(axes(n, reverse(domain))))) + return reverse_sites(n, reshape(n(domain...), length.(axes(n, reverse(domain))))) end function reverse_sites(n::OpName, a::AbstractArray) - ndomain = Int(ndims(a)//2) - perm1 = reverse(ntuple(identity, ndomain)) - perm2 = perm1 .+ ndomain - perm = (perm1..., perm2...) - return permutedims(a, perm) + ndomain = Int(ndims(a) // 2) + perm1 = reverse(ntuple(identity, ndomain)) + perm2 = perm1 .+ ndomain + perm = (perm1..., perm2...) + return permutedims(a, perm) end function state_or_op_convert( - n::StateOrOpName, - arrtype::Type{<:AbstractArray}, - domain::Tuple{Vararg{AbstractUnitRange}}, - a::AbstractArray, -) - ax = axes(n, domain) - a′ = reshape(a, length.(ax)) - a′′ = array(a′, ax) - return convert(arrtype, a′′) + n::StateOrOpName, + arrtype::Type{<:AbstractArray}, + domain::Tuple{Vararg{AbstractUnitRange}}, + a::AbstractArray, + ) + ax = axes(n, domain) + a′ = reshape(a, length.(ax)) + a′′ = array(a′, ax) + return convert(arrtype, a′′) end function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{SiteType}}) - domain′ = AbstractUnitRange.(domain) - return state_or_op_convert(n, arrtype, domain′, reversed_sites(n, domain)) + domain′ = AbstractUnitRange.(domain) + return state_or_op_convert(n, arrtype, domain′, reversed_sites(n, domain)) end function (arrtype::Type{<:AbstractArray})( - n::StateOrOpName, domain::Tuple{Vararg{AbstractUnitRange}} -) - # TODO: Make `(::OpName)(domain...)` constructor process more general inputs. - return state_or_op_convert(n, arrtype, domain, reversed_sites(n, Int.(length.(domain)))) + n::StateOrOpName, domain::Tuple{Vararg{AbstractUnitRange}} + ) + # TODO: Make `(::OpName)(domain...)` constructor process more general inputs. + return state_or_op_convert(n, arrtype, domain, reversed_sites(n, Int.(length.(domain)))) end function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...) - return arrtype(opexpr(n; kwargs...), domain...) + return arrtype(opexpr(n; kwargs...), domain...) end function op(elt::Type{<:Number}, n::String, domain...; kwargs...) - return op(AbstractArray{elt}, n, domain...; kwargs...) + return op(AbstractArray{elt}, n, domain...; kwargs...) end function op(n::String, domain...; kwargs...) - return op(AbstractArray, n, domain...; kwargs...) + return op(AbstractArray, n, domain...; kwargs...) end # Unary operations for nametype in (:StateName, :OpName) - applied = :($(nametype){:applied}) - @eval begin - nsites(n::$(applied)) = nsites(n.arg) - function (n::$(applied))(domain...) - return n.f(n.arg(domain...)) - end - end - for f in ( - :(Base.real), - :(Base.imag), - :(Base.complex), - # :(Base.adjoint), # Decide on this, since this becomes a matrix. - :(Base.:+), - :(Base.:-), - ) + applied = :($(nametype){:applied}) @eval begin - $f(n::$(nametype)) = $(applied)(; f=($f), arg=n) + nsites(n::$(applied)) = nsites(n.arg) + function (n::$(applied))(domain...) + return n.f(n.arg(domain...)) + end + end + for f in ( + :(Base.real), + :(Base.imag), + :(Base.complex), + # :(Base.adjoint), # Decide on this, since this becomes a matrix. + :(Base.:+), + :(Base.:-), + ) + @eval begin + $f(n::$(nametype)) = $(applied)(; f = ($f), arg = n) + end end - end end for nametype in (:StateName, :OpName) - kronned = :($(nametype){:kronned}) - @eval begin - nsites(n::$(kronned)) = sum(nsites, n.args) - function (n::$(kronned))(domain...) - # TODO: Generalize beyond two arguments. - # Use `cumsum(nsites.(n.args))`. - stops = cumsum(nsites.(n.args)) - starts = [1, stops[1:(end - 1)] .+ 1...] - domains = map((start, stop) -> domain[start:stop], starts, stops) - return kron(map((arg, domain) -> arg(domain...), n.args, domains)...) + kronned = :($(nametype){:kronned}) + @eval begin + nsites(n::$(kronned)) = sum(nsites, n.args) + function (n::$(kronned))(domain...) + # TODO: Generalize beyond two arguments. + # Use `cumsum(nsites.(n.args))`. + stops = cumsum(nsites.(n.args)) + starts = [1, stops[1:(end - 1)] .+ 1...] + domains = map((start, stop) -> domain[start:stop], starts, stops) + return kron(map((arg, domain) -> arg(domain...), n.args, domains)...) + end + Base.kron(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) = $(kronned)(; + args = (n1, n2, n_rest...) + ) + ⊗(n1::$(nametype), n2::$(nametype)) = kron(n1, n2) + ⊗(n1::$(kronned), n2::$(kronned)) = kron(n1.args..., n2.args...) + ⊗(n1::$(nametype), n2::$(kronned)) = kron(n1, n2.args...) + ⊗(n1::$(kronned), n2::$(nametype)) = kron(n1.args..., n2) end - Base.kron(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) = $(kronned)(; - args=(n1, n2, n_rest...) - ) - ⊗(n1::$(nametype), n2::$(nametype)) = kron(n1, n2) - ⊗(n1::$(kronned), n2::$(kronned)) = kron(n1.args..., n2.args...) - ⊗(n1::$(nametype), n2::$(kronned)) = kron(n1, n2.args...) - ⊗(n1::$(kronned), n2::$(nametype)) = kron(n1.args..., n2) - end end for nametype in (:StateName, :OpName) - summed = :($(nametype){:summed}) - @eval begin - function nsites(n::$(summed)) - # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support. - @assert allequal(nsites.(n.args)) - return nsites(first(n.args)) - end - function (n::$(summed))(domain...) - return mapreduce(a -> a(domain...), +, n.args) + summed = :($(nametype){:summed}) + @eval begin + function nsites(n::$(summed)) + # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support. + @assert allequal(nsites.(n.args)) + return nsites(first(n.args)) + end + function (n::$(summed))(domain...) + return mapreduce(a -> a(domain...), +, n.args) + end + Base.:+(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) = $(summed)(; + args = (n1, n2, n_rest...) + ) + Base.:+(n1::$(summed), n2::$(summed)) = +(n1.args..., n2.args...) + Base.:+(n1::$(summed), n2::$(nametype)) = +(n1.args..., n2) + Base.:+(n1::$(nametype), n2::$(summed)) = +(n1, n2.args...) + Base.:-(n1::$(nametype), n2::$(nametype)) = n1 + (-n2) end - Base.:+(n1::$(nametype), n2::$(nametype), n_rest::$(nametype)...) = $(summed)(; - args=(n1, n2, n_rest...) - ) - Base.:+(n1::$(summed), n2::$(summed)) = +(n1.args..., n2.args...) - Base.:+(n1::$(summed), n2::$(nametype)) = +(n1.args..., n2) - Base.:+(n1::$(nametype), n2::$(summed)) = +(n1, n2.args...) - Base.:-(n1::$(nametype), n2::$(nametype)) = n1 + (-n2) - end end for nametype in (:StateName, :OpName) - scaled = :($(nametype){:scaled}) - @eval begin - nsites(n::$(scaled)) = nsites(n.arg) - function (n::$(scaled))(domain...) - return n.arg(domain...) * n.c - end - function Base.:*(c::Number, n::$(nametype)) - return $(scaled)(; arg=n, c) - end - function Base.:*(n::$(nametype), c::Number) - return $(scaled)(; arg=n, c) - end - function Base.:/(n::$(nametype), c::Number) - return $(scaled)(; arg=n, c=inv(c)) - end - - function Base.:*(c::Number, n::$(scaled)) - return $(scaled)(; arg=n.arg, c=(c * n.c)) - end - function Base.:*(n::$(scaled), c::Number) - return $(scaled)(; arg=n.arg, c=(n.c * c)) - end - function Base.:/(n::$(scaled), c::Number) - return $(scaled)(; arg=n.arg, c=(n.c / c)) + scaled = :($(nametype){:scaled}) + @eval begin + nsites(n::$(scaled)) = nsites(n.arg) + function (n::$(scaled))(domain...) + return n.arg(domain...) * n.c + end + function Base.:*(c::Number, n::$(nametype)) + return $(scaled)(; arg = n, c) + end + function Base.:*(n::$(nametype), c::Number) + return $(scaled)(; arg = n, c) + end + function Base.:/(n::$(nametype), c::Number) + return $(scaled)(; arg = n, c = inv(c)) + end + + function Base.:*(c::Number, n::$(scaled)) + return $(scaled)(; arg = n.arg, c = (c * n.c)) + end + function Base.:*(n::$(scaled), c::Number) + return $(scaled)(; arg = n.arg, c = (n.c * c)) + end + function Base.:/(n::$(scaled), c::Number) + return $(scaled)(; arg = n.arg, c = (n.c / c)) + end end - end end # Unary operations unique to operators. for f in (:(Base.sqrt), :(Base.exp), :(Base.cis), :(Base.cos), :(Base.sin), :(Base.adjoint)) - @eval begin - $f(n::OpName) = OpName"applied"(; f=($f), arg=n) - end + @eval begin + $f(n::OpName) = OpName"applied"(; f = ($f), arg = n) + end end nsites(n::OpName"exponentiated") = nsites(n.arg) function (n::OpName"exponentiated")(domain...) - return n.arg(domain...)^n.exponent + return n.arg(domain...)^n.exponent end -Base.:^(n::OpName, exponent) = OpName"exponentiated"(; arg=n, exponent) +Base.:^(n::OpName, exponent) = OpName"exponentiated"(; arg = n, exponent) function nsites(n::OpName"producted") - # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support. - @assert allequal(nsites.(n.args)) - return nsites(first(n.args)) + # TODO: Use `allequal(nsites, n.args)` once we drop Julia 1.10 support. + @assert allequal(nsites.(n.args)) + return nsites(first(n.args)) end function (n::OpName"producted")(domain...) - return mapreduce(a -> a(domain...), *, n.args) + return mapreduce(a -> a(domain...), *, n.args) end function Base.:*(n1::OpName, n2::OpName, n_rest::OpName...) - return OpName"producted"(; args=(n1, n2, n_rest...)) + return OpName"producted"(; args = (n1, n2, n_rest...)) end Base.:*(n1::StateName"producted", n2::StateName"producted") = *(n1.args..., n2.args...) Base.:*(n1::StateName, n2::StateName"producted") = *(n1, n2.args...) Base.:*(n1::StateName"producted", n2::StateName) = *(n1.args..., n2) -controlled(n::OpName; ncontrol=1) = OpName"Controlled"(; ncontrol, arg=n) +controlled(n::OpName; ncontrol = 1) = OpName"Controlled"(; ncontrol, arg = n) using LinearAlgebra: Diagonal function (::OpName"Id")(domain) - return Diagonal(trues(to_dim(domain))) + return Diagonal(trues(to_dim(domain))) end function (n::OpName"Id")(domain1, domain_rest...) - domain = (domain1, domain_rest) - return kron(ntuple(Returns(n), length(domain))...)(domain...) + domain = (domain1, domain_rest) + return kron(ntuple(Returns(n), length(domain))...)(domain...) end @op_alias "I" "Id" @op_alias "σ0" "Id" @@ -297,24 +299,24 @@ end @op_alias "F" "Id" function (n::OpName"StandardBasis")(domain) - d = to_dim(domain) - a = falses(d, d) - a[n.index...] = one(Bool) - return a + d = to_dim(domain) + a = falses(d, d) + a[n.index...] = one(Bool) + return a end function alias(n::OpName"Proj") - return OpName"StandardBasis"(; index=(n.index, n.index)) + return OpName"StandardBasis"(; index = (n.index, n.index)) end @op_alias "proj" "Proj" function (n::OpName"a†")(domain) - d = to_dim(domain) - a = zeros(d, d) - for k in 1:(d - 1) - a[k + 1, k] = √k - end - return a + d = to_dim(domain) + a = zeros(d, d) + for k in 1:(d - 1) + a[k + 1, k] = √k + end + return a end @op_alias "Adag" "a†" @op_alias "adag" "a†" @@ -346,9 +348,9 @@ alias(::OpName"a†a†") = OpName("a†") ⊗ OpName("a†") # https://en.wikipedia.org/wiki/Generalized_Clifford_algebra # https://github.com/QuantumKitHub/MPSKitModels.jl/blob/v0.4.0/src/operators/spinoperators.jl function (n::OpName"σ⁺")(domain) - d = to_dim(domain) - s = (d - 1) / 2 - return [2 * δ(i + 1, j) * √((s + 1) * (i + j - 1) - i * j) for i in 1:d, j in 1:d] + d = to_dim(domain) + s = (d - 1) / 2 + return [2 * δ(i + 1, j) * √((s + 1) * (i + j - 1) - i * j) for i in 1:d, j in 1:d] end alias(::OpName"S+") = OpName("σ⁺") / 2 @op_alias "S⁺" "S+" @@ -383,7 +385,7 @@ alias(n::OpName"R") = cis(-(n.θ / 2) * n.arg) # Rotation around X-axis # exp(-im * θ / 2 * X) -alias(n::OpName"Rx") = OpName"R"(; params(n)..., arg=OpName"X"()) +alias(n::OpName"Rx") = OpName"R"(; params(n)..., arg = OpName"X"()) alias(::OpName"Y") = -im * (OpName"σ⁺"() - OpName"σ⁻"()) / 2 # TODO: No subsript `\_y` available @@ -408,21 +410,21 @@ alias(::OpName"Sy2") = -OpName"iSy"()^2 # Rotation around Y-axis # exp(-im * θ / 2 * Y) = exp(-θ / 2 * iY) -alias(n::OpName"Ry") = OpName"R"(; params(n)..., arg=OpName"Y"()) +alias(n::OpName"Ry") = OpName"R"(; params(n)..., arg = OpName"Y"()) # Ising (XX) coupling gate # exp(-im * θ/2 * X ⊗ X) -alias(n::OpName"Rxx") = OpName"R"(; params(n)..., arg=OpName"X"() ⊗ OpName"X"()) +alias(n::OpName"Rxx") = OpName"R"(; params(n)..., arg = OpName"X"() ⊗ OpName"X"()) @op_alias "RXX" "Rxx" # Ising (YY) coupling gate # exp(-im * θ/2 * Y ⊗ Y) -alias(n::OpName"Ryy") = OpName"R"(; params(n)..., arg=OpName"Y"() ⊗ OpName"Y"()) +alias(n::OpName"Ryy") = OpName"R"(; params(n)..., arg = OpName"Y"() ⊗ OpName"Y"()) @op_alias "RYY" "Ryy" # Ising (ZZ) coupling gate # exp(-im * θ/2 * Z ⊗ Z) -alias(n::OpName"Rzz") = OpName"R"(; params(n)..., arg=OpName"Z"() ⊗ OpName"Z"()) +alias(n::OpName"Rzz") = OpName"R"(; params(n)..., arg = OpName"Z"() ⊗ OpName"Z"()) @op_alias "RZZ" "Rzz" ## TODO: Check this definition and see if it is worth defining this. @@ -432,9 +434,9 @@ alias(n::OpName"Rzz") = OpName"R"(; params(n)..., arg=OpName"Z"() ⊗ OpName"Z"( ## @op_alias "RXY" "Rxy" function (n::OpName"σᶻ")(domain) - d = to_dim(domain) - s = (d - 1) / 2 - return Diagonal([2 * (s + 1 - i) for i in 1:d]) + d = to_dim(domain) + s = (d - 1) / 2 + return Diagonal([2 * (s + 1 - i) for i in 1:d]) end # TODO: No subsript `\_z` available # in unicode. @@ -460,24 +462,24 @@ alias(::OpName"Sz2") = OpName"Sz"()^2 # Rotation around Z-axis # exp(-im * θ / 2 * Z) -alias(n::OpName"Rz") = OpName"R"(; params(n)..., arg=OpName"Z"()) +alias(n::OpName"Rz") = OpName"R"(; params(n)..., arg = OpName"Z"()) using LinearAlgebra: eigen function (n::OpName"H")(domain) - Λ, H = eigen(OpName("X")(domain)) - p = sortperm(Λ; rev=true) - return H[:, p] + Λ, H = eigen(OpName("X")(domain)) + p = sortperm(Λ; rev = true) + return H[:, p] end using LinearAlgebra: Diagonal nsites(::OpName"SWAP") = 2 function (::OpName"SWAP")(domain1, domain2) - domain = to_dim.((domain1, domain2)) - I_matrix = Diagonal(trues(prod(domain))) - I_array = reshape(I_matrix, (domain..., domain...)) - SWAP_array = permutedims(I_array, (2, 1, 3, 4)) - SWAP_matrix = reshape(SWAP_array, (prod(domain), prod(domain))) - return SWAP_matrix + domain = to_dim.((domain1, domain2)) + I_matrix = Diagonal(trues(prod(domain))) + I_array = reshape(I_matrix, (domain..., domain...)) + SWAP_array = permutedims(I_array, (2, 1, 3, 4)) + SWAP_matrix = reshape(SWAP_array, (prod(domain), prod(domain))) + return SWAP_matrix end @op_alias "Swap" "SWAP" alias(::OpName"√SWAP") = √(OpName"SWAP"()) @@ -486,11 +488,11 @@ alias(::OpName"√SWAP") = √(OpName"SWAP"()) using LinearAlgebra: diagind nsites(::OpName"iSWAP") = 2 function (::OpName"iSWAP")(domain1, domain2) - domain = (domain1, domain2) - swap = OpName"SWAP"()(domain...) - iswap = im * swap - iswap[diagind(iswap)] .*= -im - return iswap + domain = (domain1, domain2) + swap = OpName"SWAP"()(domain...) + iswap = im * swap + iswap[diagind(iswap)] .*= -im + return iswap end @op_alias "iSwap" "iSWAP" alias(::OpName"√iSWAP") = √(OpName"iSWAP"()) @@ -503,37 +505,37 @@ alias(::OpName"√iSWAP") = √(OpName"iSWAP"()) # https://arxiv.org/abs/2410.05122 nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.arg) function (n::OpName"Controlled")(domain...) - # Number of target sites. - nt = nsites(n.arg) - # Number of control sites. - nc = get(params(n), :ncontrol, length(domain) - nt) - @assert length(domain) == nc + nt - d_control = prod(to_dim.(domain)) - prod(to_dim.(domain[(nc + 1):end])) - return cat(I(d_control), n.arg(domain[(nc + 1):end]...); dims=(1, 2)) + # Number of target sites. + nt = nsites(n.arg) + # Number of control sites. + nc = get(params(n), :ncontrol, length(domain) - nt) + @assert length(domain) == nc + nt + d_control = prod(to_dim.(domain)) - prod(to_dim.(domain[(nc + 1):end])) + return cat(I(d_control), n.arg(domain[(nc + 1):end]...); dims = (1, 2)) end @op_alias "CNOT" "Controlled" arg = OpName"X"() @op_alias "CX" "Controlled" arg = OpName"X"() @op_alias "CY" "Controlled" arg = OpName"Y"() @op_alias "CZ" "Controlled" arg = OpName"Z"() function alias(n::OpName"CPhase") - return controlled(OpName"Phase"(; params(n)...)) + return controlled(OpName"Phase"(; params(n)...)) end @op_alias "CPHASE" "CPhase" @op_alias "Cphase" "CPhase" function alias(n::OpName"CRx") - return controlled(OpName"Rx"(; params(n)...)) + return controlled(OpName"Rx"(; params(n)...)) end @op_alias "CRX" "CRx" function alias(::OpName"CRy") - return controlled(OpName"Ry"(; params(n)...)) + return controlled(OpName"Ry"(; params(n)...)) end @op_alias "CRY" "CRy" function alias(::OpName"CRz") - return controlled(OpName"Rz"(; params(n)...)) + return controlled(OpName"Rz"(; params(n)...)) end @op_alias "CRZ" "CRz" function alias(::OpName"CRn") - return controlled(; arg=OpName"Rn"(; params(n)...)) + return controlled(; arg = OpName"Rn"(; params(n)...)) end @op_alias "CRn̂" "CRn" @@ -563,36 +565,36 @@ end # Version of `sign` that returns one # if `x == 0`. function nonzero_sign(x) - iszero(x) && return one(x) - return sign(x) + iszero(x) && return one(x) + return sign(x) end function qr_positive(M::AbstractMatrix) - Q, R = qr(M) - Q′ = typeof(R)(Q) - signs = nonzero_sign.(diag(R)) - Q′ = Q′ * Diagonal(signs) - R = Diagonal(conj.(signs)) * R - return Q′, R + Q, R = qr(M) + Q′ = typeof(R)(Q) + signs = nonzero_sign.(diag(R)) + Q′ = Q′ * Diagonal(signs) + R = Diagonal(conj.(signs)) * R + return Q′, R end # TODO: Store a random matrix or seed as a parameter # of the `OpName`? function (n::OpName"RandomUnitary")(domain...) - d = prod(to_dim.(domain)) - elt = get(params(n), :eltype, Complex{Float64}) - Q, _ = qr_positive(randn(elt, (d, d))) - return Q + d = prod(to_dim.(domain)) + elt = get(params(n), :eltype, Complex{Float64}) + Q, _ = qr_positive(randn(elt, (d, d))) + return Q end @op_alias "randU" "RandomUnitary" # Expand the operator in a new basis. using LinearAlgebra: ⋅ function expand(v::AbstractArray, basis) - gramian = [basisᵢ ⋅ basisⱼ for basisᵢ in basis, basisⱼ in basis] - vbasis = [basisᵢ ⋅ v for basisᵢ in basis] - return gramian \ vbasis + gramian = [basisᵢ ⋅ basisⱼ for basisᵢ in basis, basisⱼ in basis] + vbasis = [basisᵢ ⋅ v for basisᵢ in basis] + return gramian \ vbasis end function expand(n::OpName, basis, ts...) - return expand(AbstractArray(n, ts...), map(basisᵢ -> AbstractArray(basisᵢ, ts...), basis)) + return expand(AbstractArray(n, ts...), map(basisᵢ -> AbstractArray(basisᵢ, ts...), basis)) end diff --git a/src/sitetype.jl b/src/sitetype.jl index 81c69b1..2f95a8b 100644 --- a/src/sitetype.jl +++ b/src/sitetype.jl @@ -1,8 +1,8 @@ -struct SiteType{T,Params} - params::Params - function SiteType{N}(params::NamedTuple) where {N} - return new{N,typeof(params)}(params) - end +struct SiteType{T, Params} + params::Params + function SiteType{N}(params::NamedTuple) where {N} + return new{N, typeof(params)}(params) + end end name(::SiteType{T}) where {T} = T params(t::SiteType) = getfield(t, :params) @@ -13,7 +13,7 @@ SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...)) SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...) SiteType(i::Integer; kwargs...) = SiteType{Symbol(i)}(; kwargs...) macro SiteType_str(s) - return :(SiteType{$(Expr(:quote, Symbol(s)))}) + return :(SiteType{$(Expr(:quote, Symbol(s)))}) end alias(t::SiteType) = t @@ -24,16 +24,16 @@ alias(i::Integer) = i combine_axes(a::T, b::T) where {T} = a combine_axes(a::Base.OneTo, b::Base.OneTo) = Base.OneTo{Int}(a) function combine_axes(a, b) - return UnitRange{Int}(a) + return UnitRange{Int}(a) end combine_axes(a) = a combine_axes(a, b, rest...) = combine_axes(combine_axes(a, b), rest...) -struct GradingType{T,Params} - params::Params - function GradingType{N}(params::NamedTuple) where {N} - return new{N,typeof(params)}(params) - end +struct GradingType{T, Params} + params::Params + function GradingType{N}(params::NamedTuple) where {N} + return new{N, typeof(params)}(params) + end end name(::GradingType{T}) where {T} = T params(t::GradingType) = getfield(t, :params) @@ -42,45 +42,45 @@ Base.get(t::GradingType, name::Symbol, default) = get(params(t), name, default) Base.haskey(t::GradingType, name::Symbol) = haskey(params(t), name) GradingType{N}(; kwargs...) where {N} = GradingType{N}((; kwargs...)) GradingType(s::AbstractString; kwargs...) = GradingType{Symbol(s)}(; kwargs...) -function GradingType(s::Pair{<:AbstractString,<:AbstractString}; kwargs...) - return GradingType(first(s); kwargs..., name=last(s)) +function GradingType(s::Pair{<:AbstractString, <:AbstractString}; kwargs...) + return GradingType(first(s); kwargs..., name = last(s)) end -function GradingType(s::Pair{<:AbstractString,<:NamedTuple}; kwargs...) - return GradingType(first(s); kwargs..., last(s)...) +function GradingType(s::Pair{<:AbstractString, <:NamedTuple}; kwargs...) + return GradingType(first(s); kwargs..., last(s)...) end macro GradingType_str(s) - return :(GradingType{$(Expr(:quote, Symbol(s)))}) + return :(GradingType{$(Expr(:quote, Symbol(s)))}) end function Base.AbstractUnitRange(grading::GradingType, t::SiteType) - return error("Not implemented.") + return error("Not implemented.") end function Base.AbstractUnitRange(grading::GradingType"Trivial", t::SiteType) - return Base.OneTo(length(t)) + return Base.OneTo(length(t)) end function Base.length(t::SiteType) - t′ = alias(t) - if t == t′ - return t.dim - end - return length(t′) + t′ = alias(t) + if t == t′ + return t.dim + end + return length(t′) end # TODO: Use a shorthand `(t::SiteType)() = AbstractUnitRange(t)`, # i.e. make `SiteType` callable like `OpName` and `StateName` # are right now. function Base.AbstractUnitRange(t::SiteType) - # This logic allows specifying a range with extra properties, - # like ones with symmetry sectors. - haskey(t, :range) && return t.range - if haskey(t, :gradings) - rs = map(grading -> AbstractUnitRange(GradingType(grading), t), t.gradings) - return combine_axes(Base.OneTo(length(t)), rs...) - end - return Base.OneTo(length(t)) + # This logic allows specifying a range with extra properties, + # like ones with symmetry sectors. + haskey(t, :range) && return t.range + if haskey(t, :gradings) + rs = map(grading -> AbstractUnitRange(GradingType(grading), t), t.gradings) + return combine_axes(Base.OneTo(length(t)), rs...) + end + return Base.OneTo(length(t)) end function (rangetype::Type{<:AbstractUnitRange})(t::SiteType) - return rangetype(AbstractUnitRange(t)) + return rangetype(AbstractUnitRange(t)) end Base.size(t::SiteType) = (length(t),) Base.size(t::SiteType, dim::Integer) = size(t)[dim] @@ -99,21 +99,21 @@ default_sitetype() = SiteType"Qubit"() # (t::SiteType)() = AbstractUnitRange(t) function site(rangetype::Type{<:AbstractUnitRange}, name::String; kwargs...) - return rangetype(SiteType(name; kwargs...)) + return rangetype(SiteType(name; kwargs...)) end site(name::String; kwargs...) = site(AbstractUnitRange, name; kwargs...) function sites(rangetype::Type{<:AbstractUnitRange}, name::String, positions; kwargs...) - return map(position -> site(rangetype, name; site=position, kwargs...), positions) + return map(position -> site(rangetype, name; site = position, kwargs...), positions) end function sites( - rangetype::Type{<:AbstractUnitRange}, name::String, npositions::Integer; kwargs... -) - return sites(rangetype, name, Base.OneTo(npositions); kwargs...) + rangetype::Type{<:AbstractUnitRange}, name::String, npositions::Integer; kwargs... + ) + return sites(rangetype, name, Base.OneTo(npositions); kwargs...) end function sites(name::String, positions; kwargs...) - return sites(AbstractUnitRange, name, positions; kwargs...) + return sites(AbstractUnitRange, name, positions; kwargs...) end function sites(name::String, npositions::Integer; kwargs...) - return sites(AbstractUnitRange, name, Base.OneTo(npositions); kwargs...) + return sites(AbstractUnitRange, name, Base.OneTo(npositions); kwargs...) end diff --git a/src/sitetypes/electron.jl b/src/sitetypes/electron.jl index 0de8438..cb33df7 100644 --- a/src/sitetypes/electron.jl +++ b/src/sitetypes/electron.jl @@ -46,7 +46,7 @@ alias(::OpName"ntot") = OpName"n↑"() + OpName"n↓"() @op_alias "Adn" "a↓" # c ⊗ F function (::OpName"c↓")(::SiteType"Electron") - return (OpName"C"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) + return (OpName"C"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) end @op_alias "Cdn" "c↓" @@ -55,41 +55,41 @@ end @op_alias "Adagdn" "a†↓" # c† ⊗ F function (::OpName"c†↓")(::SiteType"Electron") - return (OpName"Cdag"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) + return (OpName"Cdag"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) end @op_alias "Cdagdn" "c†↓" # F ⊗ F function (::OpName"F")(::SiteType"Electron") - return (OpName"F"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) + return (OpName"F"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) end # I ⊗ F function (::OpName"F↑")(::SiteType"Electron") - return (OpName"I"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) + return (OpName"I"() ⊗ OpName"F"())(SiteType"Fermion"(), SiteType"Fermion"()) end @op_alias "Fup" "F↑" # F ⊗ I function (::OpName"F↓")(::SiteType"Electron") - return (OpName"F"() ⊗ OpName"I"())(SiteType"Fermion"(), SiteType"Fermion"()) + return (OpName"F"() ⊗ OpName"I"())(SiteType"Fermion"(), SiteType"Fermion"()) end @op_alias "Fdn" "F↓" function (n::OpName"SpinOp")(::SiteType"Electron") - return cat(falses(1, 1), n.arg(2), falses(1, 1); dims=(1, 2)) + return cat(falses(1, 1), n.arg(2), falses(1, 1); dims = (1, 2)) end # These implicitly define other spin operators. # TODO: Maybe require calling it as `OpName("SpinOp"; arg=OpName("Sz"))`? function (n::OpName"σ⁺")(domain::SiteType"Electron") - return OpName"SpinOp"(; arg=n)(domain) + return OpName"SpinOp"(; arg = n)(domain) end function (n::OpName"σᶻ")(domain::SiteType"Electron") - return OpName"SpinOp"(; arg=n)(domain) + return OpName"SpinOp"(; arg = n)(domain) end function (n::OpName"R")(domain::SiteType"Electron") - return OpName"SpinOp"(; arg=n)(domain) + return OpName"SpinOp"(; arg = n)(domain) end has_fermion_string(::OpName"c↑", ::SiteType"Electron") = true diff --git a/src/sitetypes/qubit.jl b/src/sitetypes/qubit.jl index a322671..3b097ba 100644 --- a/src/sitetypes/qubit.jl +++ b/src/sitetypes/qubit.jl @@ -42,24 +42,24 @@ Base.length(::SiteType"Qubit") = 2 # SIC-POVMs (::StateName"Tetra0")(::SiteType"Qubit") = [ - 1 - 0 + 1 + 0 ] (::StateName"Tetra2")(::SiteType"Qubit") = [ - 1 / √3 - √2 / √3 + 1 / √3 + √2 / √3 ] (::StateName"Tetra3")(::SiteType"Qubit") = [ - 1 / √3 - √2 / √3 * exp(im * 2π / 3) + 1 / √3 + √2 / √3 * exp(im * 2π / 3) ] (::StateName"Tetra4")(::SiteType"Qubit") = [ - 1 / √3 - √2 / √3 * exp(im * 4π / 3) + 1 / √3 + √2 / √3 * exp(im * 4π / 3) ] # TODO: Define as `(I + σᶻ) / 2`? -(::OpName"ProjUp")(::SiteType"Qubit") = OpName"Proj"(; index=1)(2) +(::OpName"ProjUp")(::SiteType"Qubit") = OpName"Proj"(; index = 1)(2) (::OpName"projUp")(domain::SiteType"Qubit") = OpName"ProjUp"()(domain) (::OpName"Proj↑")(domain::SiteType"Qubit") = OpName"ProjUp"()(domain) (::OpName"proj↑")(domain::SiteType"Qubit") = OpName"ProjUp"()(domain) @@ -68,7 +68,7 @@ Base.length(::SiteType"Qubit") = 2 # TODO: Define as `σ⁺ * σ⁻`? # TODO: Define as `(I - σᶻ) / 2`? -(::OpName"ProjDn")(::SiteType"Qubit") = OpName"Proj"(; index=2)(2) +(::OpName"ProjDn")(::SiteType"Qubit") = OpName"Proj"(; index = 2)(2) (::OpName"projDn")(domain::SiteType"Qubit") = OpName"ProjDn"()(domain) (::OpName"Proj↓")(domain::SiteType"Qubit") = OpName"ProjDn"()(domain) (::OpName"proj↓")(domain::SiteType"Qubit") = OpName"ProjDn"()(domain) diff --git a/src/sitetypes/qudit.jl b/src/sitetypes/qudit.jl index 7205c98..a5ce3c2 100644 --- a/src/sitetypes/qudit.jl +++ b/src/sitetypes/qudit.jl @@ -1,3 +1,3 @@ Base.length(t::SiteType"Qudit") = t.dim -alias(t::SiteType"Boson") = SiteType"Qudit"(; dim=t.dim) -alias(t::SiteType"S") = SiteType"Qudit"(; dim=Int(2t.spin + 1)) +alias(t::SiteType"Boson") = SiteType"Qudit"(; dim = t.dim) +alias(t::SiteType"S") = SiteType"Qudit"(; dim = Int(2t.spin + 1)) diff --git a/src/sitetypes/tj.jl b/src/sitetypes/tj.jl index 4efa7ef..9f7e969 100644 --- a/src/sitetypes/tj.jl +++ b/src/sitetypes/tj.jl @@ -1,11 +1,11 @@ Base.length(::SiteType"tJ") = 3 function (n::StateName)(domain::SiteType"tJ") - return n(SiteType"Electron"())[1:length(domain)] + return n(SiteType"Electron"())[1:length(domain)] end function (n::OpName)(domain::SiteType"tJ") - return n(SiteType"Electron"())[1:length(domain), 1:length(domain)] + return n(SiteType"Electron"())[1:length(domain), 1:length(domain)] end has_fermion_string(::OpName"c↑", ::SiteType"tJ") = true diff --git a/src/state.jl b/src/state.jl index 56262c3..3dc1d17 100644 --- a/src/state.jl +++ b/src/state.jl @@ -1,10 +1,10 @@ using Random: randstring -struct StateName{Name,Params} - params::Params - function StateName{N}(params::NamedTuple) where {N} - return new{N,typeof(params)}(params) - end +struct StateName{Name, Params} + params::Params + function StateName{N}(params::NamedTuple) where {N} + return new{N, typeof(params)}(params) + end end name(::StateName{N}) where {N} = N params(n::StateName) = getfield(n, :params) @@ -16,20 +16,22 @@ StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...)) StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...) StateName(s::Symbol; kwargs...) = StateName{s}(; kwargs...) macro StateName_str(s) - return StateName{Symbol(s)} + return StateName{Symbol(s)} end # Handles special case `state(1) == [1, 0]`. Note the # one-based indexing, as opposed to `state("0") == [1, 0]`. -StateName(i::Integer; kwargs...) = StateName"StandardBasis"(; index=i, kwargs...) +StateName(i::Integer; kwargs...) = StateName"StandardBasis"(; index = i, kwargs...) function state_alias_expr(name1, name2, pars...) - return :(function alias(n::StateName{Symbol($name1)}) - return StateName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) - end) + return :( + function alias(n::StateName{Symbol($name1)}) + return StateName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) + end + ) end macro state_alias(name1, name2, params...) - return state_alias_expr(name1, name2) + return state_alias_expr(name1, name2) end # This compiles operator expressions, such as: @@ -37,7 +39,7 @@ end # stateexpr("0 + 1") == StateName("0") + StateName("1") # ``` function stateexpr(n::String; kwargs...) - return state_or_op_expr(StateName, n; kwargs...) + return state_or_op_expr(StateName, n; kwargs...) end # Handles special case `state(1) == [1, 0]`. Note the @@ -53,103 +55,103 @@ const MINUS_STRING = randcharstring() const VERTICALBAR_STRING = randcharstring() const RANGLE_STRING = randcharstring() const EXPR_REPLACEMENTS_1 = ( - "†" => DAGGER_STRING, - "↑" => UPARROW_STRING, - "↓" => DOWNARROW_STRING, - # Replace trailing plus and minus characters - # in operators, which don't parse properly. - r"(\S)\+" => SubstitutionString("\\1$(PLUS_STRING)"), - r"(\S)\-" => SubstitutionString("\\1$(MINUS_STRING)"), + "†" => DAGGER_STRING, + "↑" => UPARROW_STRING, + "↓" => DOWNARROW_STRING, + # Replace trailing plus and minus characters + # in operators, which don't parse properly. + r"(\S)\+" => SubstitutionString("\\1$(PLUS_STRING)"), + r"(\S)\-" => SubstitutionString("\\1$(MINUS_STRING)"), ) const EXPR_REPLACEMENTS_2 = ( - r"\|(\S*)⟩" => SubstitutionString("$(VERTICALBAR_STRING)\\1$(RANGLE_STRING)"), + r"\|(\S*)⟩" => SubstitutionString("$(VERTICALBAR_STRING)\\1$(RANGLE_STRING)"), ) const INVERSE_EXPR_REPLACEMENTS = ( - DAGGER_STRING => "†", - UPARROW_STRING => "↑", - DOWNARROW_STRING => "↓", - PLUS_STRING => "+", - MINUS_STRING => "-", - # We remove the bra-ket notation and search for states - # with names stored inside. - VERTICALBAR_STRING => "", - RANGLE_STRING => "", + DAGGER_STRING => "†", + UPARROW_STRING => "↑", + DOWNARROW_STRING => "↓", + PLUS_STRING => "+", + MINUS_STRING => "-", + # We remove the bra-ket notation and search for states + # with names stored inside. + VERTICALBAR_STRING => "", + RANGLE_STRING => "", ) function state_or_op_expr(ntype::Type, n::String; kwargs...) - # Do this in two rounds since for some reason - # one round doesn't work for expressions - # like `"|+⟩"`. - n = replace(n, EXPR_REPLACEMENTS_1...) - n = replace(n, EXPR_REPLACEMENTS_2...) - depth = 1 - return state_or_op_expr(ntype, Meta.parse(n), depth; kwargs...) + # Do this in two rounds since for some reason + # one round doesn't work for expressions + # like `"|+⟩"`. + n = replace(n, EXPR_REPLACEMENTS_1...) + n = replace(n, EXPR_REPLACEMENTS_2...) + depth = 1 + return state_or_op_expr(ntype, Meta.parse(n), depth; kwargs...) end function state_or_op_expr(ntype::Type, n::Number, depth::Int; kwargs...) - if depth == 1 - return ntype{Symbol(n)}(; kwargs...) - end - return n + if depth == 1 + return ntype{Symbol(n)}(; kwargs...) + end + return n end function state_or_op_expr(ntype::Type, n::Symbol, depth::Int; kwargs...) - n === :im && return im - n === :π && return π - n = Symbol(replace(String(n), INVERSE_EXPR_REPLACEMENTS...)) - return ntype{n}(; kwargs...) + n === :im && return im + n === :π && return π + n = Symbol(replace(String(n), INVERSE_EXPR_REPLACEMENTS...)) + return ntype{n}(; kwargs...) end function state_or_op_expr(ntype::Type, ex::Expr, depth::Int) - if Meta.isexpr(ex, :call) - return eval(ex.args[1])(state_or_op_expr.(ntype, ex.args[2:end], depth + 1)...) - end - if Meta.isexpr(ex, :curly) - # Syntax for parametrized gates, i.e. - # `state_or_op_expr("Ry{θ=π/2}")`. - params = ex.args[2:end] - kwargs = Dict( - map(params) do param - @assert Meta.isexpr(param, :(=)) - return param.args[1] => eval(param.args[2]) - end, - ) - return ntype{ex.args[1]}(; kwargs...) - end - return error("Can't parse expression $ex.") + if Meta.isexpr(ex, :call) + return eval(ex.args[1])(state_or_op_expr.(ntype, ex.args[2:end], depth + 1)...) + end + if Meta.isexpr(ex, :curly) + # Syntax for parametrized gates, i.e. + # `state_or_op_expr("Ry{θ=π/2}")`. + params = ex.args[2:end] + kwargs = Dict( + map(params) do param + @assert Meta.isexpr(param, :(=)) + return param.args[1] => eval(param.args[2]) + end, + ) + return ntype{ex.args[1]}(; kwargs...) + end + return error("Can't parse expression $ex.") end function Base.axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}}) - return domain + return domain end function reverse_sites(n::StateName, a::AbstractArray) - perm = reverse(ntuple(identity, ndims(a))) - return permutedims(a, perm) + perm = reverse(ntuple(identity, ndims(a))) + return permutedims(a, perm) end -function state(arrtype::Type{<:AbstractArray}, n::Union{Int,String}, domain...; kwargs...) - return arrtype(stateexpr(n; kwargs...), domain...) +function state(arrtype::Type{<:AbstractArray}, n::Union{Int, String}, domain...; kwargs...) + return arrtype(stateexpr(n; kwargs...), domain...) end -function state(elt::Type{<:Number}, n::Union{Int,String}, domain...; kwargs...) - return state(AbstractArray{elt}, n, domain...; kwargs...) +function state(elt::Type{<:Number}, n::Union{Int, String}, domain...; kwargs...) + return state(AbstractArray{elt}, n, domain...; kwargs...) end -function state(n::Union{Int,String}, domain...; kwargs...) - return state(AbstractArray, n, domain...; kwargs...) +function state(n::Union{Int, String}, domain...; kwargs...) + return state(AbstractArray, n, domain...; kwargs...) end function (n::StateName"StandardBasis")(domain) - a = falses(to_dim(domain)) - a[n.index] = one(Bool) - return a + a = falses(to_dim(domain)) + a[n.index] = one(Bool) + return a end function (n::StateName{N})(domain...) where {N} - # TODO: Try one alias at a time? - # TODO: First call `alias(n, domain...)` - # to allow for aliases specific to certain - # SiteTypes? - n′ = alias(n) - domain′ = alias.(domain) - if n == n′ && domain′ == domain - index = parse(Int, String(N)) + 1 - return StateName"StandardBasis"(; index)(domain...) - end - return n′(domain′...) + # TODO: Try one alias at a time? + # TODO: First call `alias(n, domain...)` + # to allow for aliases specific to certain + # SiteTypes? + n′ = alias(n) + domain′ = alias.(domain) + if n == n′ && domain′ == domain + index = parse(Int, String(N)) + 1 + return StateName"StandardBasis"(; index)(domain...) + end + return n′(domain′...) end diff --git a/test/runtests.jl b/test/runtests.jl index 98b2d2b..0008050 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,60 +6,62 @@ using Suppressor: Suppressor const pat = r"(?:--group=)(\w+)" arg_id = findfirst(contains(pat), ARGS) const GROUP = uppercase( - if isnothing(arg_id) - get(ENV, "GROUP", "ALL") - else - only(match(pat, ARGS[arg_id]).captures) - end, + if isnothing(arg_id) + get(ENV, "GROUP", "ALL") + else + only(match(pat, ARGS[arg_id]).captures) + end, ) "match files of the form `test_*.jl`, but exclude `*setup*.jl`" function istestfile(fn) - return endswith(fn, ".jl") && startswith(basename(fn), "test_") && !contains(fn, "setup") + return endswith(fn, ".jl") && startswith(basename(fn), "test_") && !contains(fn, "setup") end "match files of the form `*.jl`, but exclude `*_notest.jl` and `*setup*.jl`" function isexamplefile(fn) - return endswith(fn, ".jl") && !endswith(fn, "_notest.jl") && !contains(fn, "setup") + return endswith(fn, ".jl") && !endswith(fn, "_notest.jl") && !contains(fn, "setup") end @time begin - # tests in groups based on folder structure - for testgroup in filter(isdir, readdir(@__DIR__)) - if GROUP == "ALL" || GROUP == uppercase(testgroup) - groupdir = joinpath(@__DIR__, testgroup) - for file in filter(istestfile, readdir(groupdir)) - filename = joinpath(groupdir, file) - @eval @safetestset $file begin - include($filename) + # tests in groups based on folder structure + for testgroup in filter(isdir, readdir(@__DIR__)) + if GROUP == "ALL" || GROUP == uppercase(testgroup) + groupdir = joinpath(@__DIR__, testgroup) + for file in filter(istestfile, readdir(groupdir)) + filename = joinpath(groupdir, file) + @eval @safetestset $file begin + include($filename) + end + end end - end end - end - # single files in top folder - for file in filter(istestfile, readdir(@__DIR__)) - (file == basename(@__FILE__)) && continue # exclude this file to avoid infinite recursion - @eval @safetestset $file begin - include($file) + # single files in top folder + for file in filter(istestfile, readdir(@__DIR__)) + (file == basename(@__FILE__)) && continue # exclude this file to avoid infinite recursion + @eval @safetestset $file begin + include($file) + end end - end - # test examples - examplepath = joinpath(@__DIR__, "..", "examples") - for (root, _, files) in walkdir(examplepath) - contains(chopprefix(root, @__DIR__), "setup") && continue - for file in filter(isexamplefile, files) - filename = joinpath(root, file) - @eval begin - @safetestset $file begin - $(Expr( - :macrocall, - GlobalRef(Suppressor, Symbol("@suppress")), - LineNumberNode(@__LINE__, @__FILE__), - :(include($filename)), - )) + # test examples + examplepath = joinpath(@__DIR__, "..", "examples") + for (root, _, files) in walkdir(examplepath) + contains(chopprefix(root, @__DIR__), "setup") && continue + for file in filter(isexamplefile, files) + filename = joinpath(root, file) + @eval begin + @safetestset $file begin + $( + Expr( + :macrocall, + GlobalRef(Suppressor, Symbol("@suppress")), + LineNumberNode(@__LINE__, @__FILE__), + :(include($filename)), + ) + ) + end + end end - end end - end end diff --git a/test/test_aqua.jl b/test/test_aqua.jl index eedff41..702d0f9 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -3,5 +3,5 @@ using Aqua: Aqua using Test: @testset @testset "Code quality (Aqua.jl)" begin - Aqua.test_all(QuantumOperatorDefinitions; ambiguities=false) + Aqua.test_all(QuantumOperatorDefinitions; ambiguities = false) end diff --git a/test/test_basics.jl b/test/test_basics.jl index 0df73cf..a2e1b37 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -1,5 +1,5 @@ using QuantumOperatorDefinitions: - OpName, SiteType, ⊗, expand, nsites, op, opexpr, site, sites, state + OpName, SiteType, ⊗, expand, nsites, op, opexpr, site, sites, state using LinearAlgebra: Diagonal, I using Test: @test, @testset @@ -8,304 +8,304 @@ const complex_elts = map(elt -> Complex{elt}, real_elts) const elts = (real_elts..., complex_elts...) @testset "QuantumOperatorDefinitions" begin - @testset "Qubit/Qudit" begin - # https://en.wikipedia.org/wiki/Pauli_matrices - # https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins - for ( - t, - len, - Xmat, - Ymat, - Zmat, - Nmat, - SWAPmat, - iSWAPmat, - RXXmat, - RYYmat, - RZZmat, - Proj1mat, - Proj2mat, - StandardBasis12mat, - ) in ( - ( - SiteType("Qubit"), - 2, - [0 1; 1 0], - [0 -im; im 0], - [1 0; 0 -1], - [0 0; 0 1], - reshape([1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1], (2, 2, 2, 2)), - reshape([1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1], (2, 2, 2, 2)), - (_, θ) -> reshape( - [ - cos(θ / 2) 0 0 -im*sin(θ / 2) - 0 cos(θ / 2) -im*sin(θ / 2) 0 - 0 -im*sin(θ / 2) cos(θ / 2) 0 - -im*sin(θ / 2) 0 0 cos(θ / 2) - ], - (2, 2, 2, 2), - ), - (_, θ) -> reshape( - [ - cos(θ / 2) 0 0 im*sin(θ / 2) - 0 cos(θ / 2) -im*sin(θ / 2) 0 - 0 -im*sin(θ / 2) cos(θ / 2) 0 - im*sin(θ / 2) 0 0 cos(θ / 2) - ], - (2, 2, 2, 2), - ), - (_, θ) -> reshape( - Diagonal([exp(-im * θ / 2), exp(im * θ / 2), exp(im * θ / 2), exp(-im * θ / 2)]), - (2, 2, 2, 2), - ), - [1 0; 0 0], - [0 0; 0 1], - [0 1; 0 0], - ), - ( - SiteType("Qudit"; dim=3), - 3, - √2 * [0 1 0; 1 0 1; 0 1 0], - √2 * [0 -im 0; im 0 -im; 0 im 0], - 2 * [1 0 0; 0 0 0; 0 0 -1], - [0 0 0; 0 1 0; 0 0 2], - reshape( - [ - 1 0 0 0 0 0 0 0 0 - 0 0 0 1 0 0 0 0 0 - 0 0 0 0 0 0 1 0 0 - 0 1 0 0 0 0 0 0 0 - 0 0 0 0 1 0 0 0 0 - 0 0 0 0 0 0 0 1 0 - 0 0 1 0 0 0 0 0 0 - 0 0 0 0 0 1 0 0 0 - 0 0 0 0 0 0 0 0 1 - ], - (3, 3, 3, 3), - ), - reshape( - [ - 1 0 0 0 0 0 0 0 0 - 0 0 0 im 0 0 0 0 0 - 0 0 0 0 0 0 im 0 0 - 0 im 0 0 0 0 0 0 0 - 0 0 0 0 1 0 0 0 0 - 0 0 0 0 0 0 0 im 0 - 0 0 im 0 0 0 0 0 0 - 0 0 0 0 0 im 0 0 0 - 0 0 0 0 0 0 0 0 1 - ], - (3, 3, 3, 3), - ), - (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), - (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), - (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), - [1 0 0; 0 0 0; 0 0 0], - [0 0 0; 0 1 0; 0 0 0], - [0 1 0; 0 0 0; 0 0 0], - ), - ) - @test length(t) == len - @test size(t) == (len,) - @test size(t, 1) == len - @test axes(t) == (Base.OneTo(len),) - @test axes(t, 1) == Base.OneTo(len) - for (o, nbits, elts, ref) in ( - (OpName("X"), 1, elts, Xmat), - (OpName("σˣ"), 1, elts, Xmat), - (OpName("√X"), 1, complex_elts, √Xmat), - (OpName("iX"), 1, complex_elts, im * Xmat), - (OpName("Y"), 1, complex_elts, Ymat), - (OpName("σʸ"), 1, complex_elts, Ymat), - (OpName("iY"), 1, elts, im * Ymat), - (OpName("Z"), 1, elts, Zmat), - (OpName("σᶻ"), 1, elts, Zmat), - (OpName("iZ"), 1, complex_elts, im * Zmat), - (OpName("N"), 1, elts, Nmat), - (OpName("n"), 1, elts, Nmat), - (OpName("Phase"; θ=π / 3), 1, complex_elts, exp(im * π / 3 * Nmat)), - (OpName("π/8"), 1, complex_elts, exp(im * π / 4 * Nmat)), - (OpName("Rx"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Xmat)), - (OpName("Ry"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Ymat)), - (OpName("Rz"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Zmat)), - (OpName("SWAP"), 2, elts, SWAPmat), - # (OpName("√SWAP"), 2, complex_elts, √SWAPmat), - (OpName("iSWAP"), 2, complex_elts, iSWAPmat), - # (OpName("√iSWAP"), 2, complex_elts, √iSWAPmat), - (OpName("Rxx"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), - (OpName("RXX"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), - (OpName("Ryy"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), - (OpName("RYY"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), - (OpName("Rzz"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), - (OpName("RZZ"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), - (OpName("Proj"; index=1), 1, elts, Proj1mat), - (OpName("Proj"; index=2), 1, elts, Proj2mat), - (OpName("StandardBasis"; index=(1, 2)), 1, elts, StandardBasis12mat), - ) - @test nsites(o) == nbits - for arraytype in (AbstractArray, Array) - for elt in elts - ts = ntuple(Returns(t), nbits) - lens = ntuple(Returns(len), nbits) - for domain in (ts, (ts,), lens, (lens,)) - @test arraytype(o, domain...) ≈ ref - @test arraytype{elt}(o, domain...) ≈ ref - @test eltype(arraytype{elt}(o, domain...)) === elt + @testset "Qubit/Qudit" begin + # https://en.wikipedia.org/wiki/Pauli_matrices + # https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins + for ( + t, + len, + Xmat, + Ymat, + Zmat, + Nmat, + SWAPmat, + iSWAPmat, + RXXmat, + RYYmat, + RZZmat, + Proj1mat, + Proj2mat, + StandardBasis12mat, + ) in ( + ( + SiteType("Qubit"), + 2, + [0 1; 1 0], + [0 -im; im 0], + [1 0; 0 -1], + [0 0; 0 1], + reshape([1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1], (2, 2, 2, 2)), + reshape([1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1], (2, 2, 2, 2)), + (_, θ) -> reshape( + [ + cos(θ / 2) 0 0 -im * sin(θ / 2) + 0 cos(θ / 2) -im * sin(θ / 2) 0 + 0 -im * sin(θ / 2) cos(θ / 2) 0 + -im * sin(θ / 2) 0 0 cos(θ / 2) + ], + (2, 2, 2, 2), + ), + (_, θ) -> reshape( + [ + cos(θ / 2) 0 0 im * sin(θ / 2) + 0 cos(θ / 2) -im * sin(θ / 2) 0 + 0 -im * sin(θ / 2) cos(θ / 2) 0 + im * sin(θ / 2) 0 0 cos(θ / 2) + ], + (2, 2, 2, 2), + ), + (_, θ) -> reshape( + Diagonal([exp(-im * θ / 2), exp(im * θ / 2), exp(im * θ / 2), exp(-im * θ / 2)]), + (2, 2, 2, 2), + ), + [1 0; 0 0], + [0 0; 0 1], + [0 1; 0 0], + ), + ( + SiteType("Qudit"; dim = 3), + 3, + √2 * [0 1 0; 1 0 1; 0 1 0], + √2 * [0 -im 0; im 0 -im; 0 im 0], + 2 * [1 0 0; 0 0 0; 0 0 -1], + [0 0 0; 0 1 0; 0 0 2], + reshape( + [ + 1 0 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 0 0 1 + ], + (3, 3, 3, 3), + ), + reshape( + [ + 1 0 0 0 0 0 0 0 0 + 0 0 0 im 0 0 0 0 0 + 0 0 0 0 0 0 im 0 0 + 0 im 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 im 0 + 0 0 im 0 0 0 0 0 0 + 0 0 0 0 0 im 0 0 0 + 0 0 0 0 0 0 0 0 1 + ], + (3, 3, 3, 3), + ), + (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), + (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), + (O, θ) -> reshape(exp(-im * (θ / 2) * kron(O, O)), (3, 3, 3, 3)), + [1 0 0; 0 0 0; 0 0 0], + [0 0 0; 0 1 0; 0 0 0], + [0 1 0; 0 0 0; 0 0 0], + ), + ) + @test length(t) == len + @test size(t) == (len,) + @test size(t, 1) == len + @test axes(t) == (Base.OneTo(len),) + @test axes(t, 1) == Base.OneTo(len) + for (o, nbits, elts, ref) in ( + (OpName("X"), 1, elts, Xmat), + (OpName("σˣ"), 1, elts, Xmat), + (OpName("√X"), 1, complex_elts, √Xmat), + (OpName("iX"), 1, complex_elts, im * Xmat), + (OpName("Y"), 1, complex_elts, Ymat), + (OpName("σʸ"), 1, complex_elts, Ymat), + (OpName("iY"), 1, elts, im * Ymat), + (OpName("Z"), 1, elts, Zmat), + (OpName("σᶻ"), 1, elts, Zmat), + (OpName("iZ"), 1, complex_elts, im * Zmat), + (OpName("N"), 1, elts, Nmat), + (OpName("n"), 1, elts, Nmat), + (OpName("Phase"; θ = π / 3), 1, complex_elts, exp(im * π / 3 * Nmat)), + (OpName("π/8"), 1, complex_elts, exp(im * π / 4 * Nmat)), + (OpName("Rx"; θ = π / 3), 1, complex_elts, exp(-im * π / 6 * Xmat)), + (OpName("Ry"; θ = π / 3), 1, complex_elts, exp(-im * π / 6 * Ymat)), + (OpName("Rz"; θ = π / 3), 1, complex_elts, exp(-im * π / 6 * Zmat)), + (OpName("SWAP"), 2, elts, SWAPmat), + # (OpName("√SWAP"), 2, complex_elts, √SWAPmat), + (OpName("iSWAP"), 2, complex_elts, iSWAPmat), + # (OpName("√iSWAP"), 2, complex_elts, √iSWAPmat), + (OpName("Rxx"; θ = π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), + (OpName("RXX"; θ = π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), + (OpName("Ryy"; θ = π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), + (OpName("RYY"; θ = π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), + (OpName("Rzz"; θ = π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), + (OpName("RZZ"; θ = π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), + (OpName("Proj"; index = 1), 1, elts, Proj1mat), + (OpName("Proj"; index = 2), 1, elts, Proj2mat), + (OpName("StandardBasis"; index = (1, 2)), 1, elts, StandardBasis12mat), + ) + @test nsites(o) == nbits + for arraytype in (AbstractArray, Array) + for elt in elts + ts = ntuple(Returns(t), nbits) + lens = ntuple(Returns(len), nbits) + for domain in (ts, (ts,), lens, (lens,)) + @test arraytype(o, domain...) ≈ ref + @test arraytype{elt}(o, domain...) ≈ ref + @test eltype(arraytype{elt}(o, domain...)) === elt + end + end + end end - end end - end end - end - @testset "op parsing" begin - @test Matrix(opexpr("X * Y")) == op("X") * op("Y") - @test op("X * Y") == op("X") * op("Y") - @test op("X * Y + Z") == op("X") * op("Y") + op("Z") - @test op("X * Y + 2 * Z") == op("X") * op("Y") + 2 * op("Z") - @test op("exp(im * (X * Y + 2 * Z))") == exp(im * (op("X") * op("Y") + 2 * op("Z"))) - @test op("exp(im * (X ⊗ Y + Z ⊗ Z))") == permutedims( - reshape(exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))), (2, 2, 2, 2)), - (2, 1, 4, 3), - ) - @test op("Ry{θ=π/2}") == op("Ry"; θ=π / 2) - # Awkward parsing corner cases. - @test op("S+") == Matrix(OpName("S+")) - @test op("S-") == Matrix(OpName("S-")) - @test op("S+ + S-") == Matrix(OpName("S+") + OpName("S-")) - @test op("S+ - S-") == Matrix(OpName("S+") - OpName("S-")) - @test op("a†") == Matrix(OpName("a†")) - for name in ("c↑", "c†↑", "c↓", "c†↓") - @test op(name, SiteType("Electron")) == Matrix(OpName(name), SiteType("Electron")) + @testset "op parsing" begin + @test Matrix(opexpr("X * Y")) == op("X") * op("Y") + @test op("X * Y") == op("X") * op("Y") + @test op("X * Y + Z") == op("X") * op("Y") + op("Z") + @test op("X * Y + 2 * Z") == op("X") * op("Y") + 2 * op("Z") + @test op("exp(im * (X * Y + 2 * Z))") == exp(im * (op("X") * op("Y") + 2 * op("Z"))) + @test op("exp(im * (X ⊗ Y + Z ⊗ Z))") == permutedims( + reshape(exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))), (2, 2, 2, 2)), + (2, 1, 4, 3), + ) + @test op("Ry{θ=π/2}") == op("Ry"; θ = π / 2) + # Awkward parsing corner cases. + @test op("S+") == Matrix(OpName("S+")) + @test op("S-") == Matrix(OpName("S-")) + @test op("S+ + S-") == Matrix(OpName("S+") + OpName("S-")) + @test op("S+ - S-") == Matrix(OpName("S+") - OpName("S-")) + @test op("a†") == Matrix(OpName("a†")) + for name in ("c↑", "c†↑", "c↓", "c†↓") + @test op(name, SiteType("Electron")) == Matrix(OpName(name), SiteType("Electron")) + end end - end - @testset "Random unitary" begin - U = op("RandomUnitary") - @test eltype(U) === Complex{Float64} - @test U * U' ≈ I(2) - - U = op("RandomUnitary"; eltype=Float32) - @test eltype(U) === Float32 - @test U * U' ≈ I(2) - - U = op("RandomUnitary", 3) - @test eltype(U) === Complex{Float64} - @test U * U' ≈ I(3) - end - @testset "state" begin - @test state(1) == [1, 0] - @test state("0") == [1, 0] - @test state(2) == [0, 1] - @test state("1") == [0, 1] - @test state(1, 3) == [1, 0, 0] - @test state("0", 3) == [1, 0, 0] - @test state(2, 3) == [0, 1, 0] - @test state("1", 3) == [0, 1, 0] - @test state(3, 3) == [0, 0, 1] - @test state("2", 3) == [0, 0, 1] + @testset "Random unitary" begin + U = op("RandomUnitary") + @test eltype(U) === Complex{Float64} + @test U * U' ≈ I(2) - @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+") - @test state("|0⟩ ⊗ |+⟩") == - permutedims(reshape(kron(state("0"), state("+")), (2, 2)), (2, 1)) - end - @testset "Ranges" begin - @test AbstractUnitRange(SiteType("S=1/2")) == Base.OneTo(2) - @test UnitRange{Int32}(SiteType("S=1/2")) == Base.OneTo(2) - @test UnitRange{Int32}(SiteType("S=1/2")) isa UnitRange{Int32} - @test AbstractUnitRange(SiteType("Qudit"; dim=3)) == Base.OneTo(3) - @test UnitRange{Int32}(SiteType("Qudit"; dim=3)) == Base.OneTo(3) - @test UnitRange{Int32}(SiteType("Qudit"; dim=3)) isa UnitRange{Int32} + U = op("RandomUnitary"; eltype = Float32) + @test eltype(U) === Float32 + @test U * U' ≈ I(2) - @test op("X", 2) == op("X", Base.OneTo(2)) - @test op("X", 3) == op("X", Base.OneTo(3)) - end - @testset "site, sites" begin - s = site("Qudit"; dim=3) - @test s === Base.OneTo(3) + U = op("RandomUnitary", 3) + @test eltype(U) === Complex{Float64} + @test U * U' ≈ I(3) + end + @testset "state" begin + @test state(1) == [1, 0] + @test state("0") == [1, 0] + @test state(2) == [0, 1] + @test state("1") == [0, 1] + @test state(1, 3) == [1, 0, 0] + @test state("0", 3) == [1, 0, 0] + @test state(2, 3) == [0, 1, 0] + @test state("1", 3) == [0, 1, 0] + @test state(3, 3) == [0, 0, 1] + @test state("2", 3) == [0, 0, 1] - s = site(AbstractUnitRange{Int32}, "Qudit"; dim=3) - @test s === Base.OneTo(Int32(3)) + @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+") + @test state("|0⟩ ⊗ |+⟩") == + permutedims(reshape(kron(state("0"), state("+")), (2, 2)), (2, 1)) + end + @testset "Ranges" begin + @test AbstractUnitRange(SiteType("S=1/2")) == Base.OneTo(2) + @test UnitRange{Int32}(SiteType("S=1/2")) == Base.OneTo(2) + @test UnitRange{Int32}(SiteType("S=1/2")) isa UnitRange{Int32} + @test AbstractUnitRange(SiteType("Qudit"; dim = 3)) == Base.OneTo(3) + @test UnitRange{Int32}(SiteType("Qudit"; dim = 3)) == Base.OneTo(3) + @test UnitRange{Int32}(SiteType("Qudit"; dim = 3)) isa UnitRange{Int32} - for ss in (sites("Qudit", 4; dim=3), sites("Qudit", 2:5; dim=3)) - @test length(ss) == 4 - @test ss == map(Returns(Base.OneTo(3)), 1:4) - for s in ss - @test s === Base.OneTo(3) - end + @test op("X", 2) == op("X", Base.OneTo(2)) + @test op("X", 3) == op("X", Base.OneTo(3)) end + @testset "site, sites" begin + s = site("Qudit"; dim = 3) + @test s === Base.OneTo(3) - for ss in ( - sites(AbstractUnitRange{Int32}, "Qudit", 4; dim=3), - sites(AbstractUnitRange{Int32}, "Qudit", 2:5; dim=3), - ) - @test length(ss) == 4 - @test ss == map(Returns(Base.OneTo(3)), 1:4) - for s in ss + s = site(AbstractUnitRange{Int32}, "Qudit"; dim = 3) @test s === Base.OneTo(Int32(3)) - end - end - ss = sites("Qudit", (1, 2, 3); dim=3) - @test ss === (Base.OneTo(3), Base.OneTo(3), Base.OneTo(3)) - end - @testset "Electron/tJ" begin - for (ns, x) in ( - (("0", "Emp"), [1, 0, 0, 0]), - (("↑", "Up"), [0, 1, 0, 0]), - (("↓", "Dn"), [0, 0, 1, 0]), - (("↑↓", "UpDn"), [0, 0, 0, 1]), - ) - for n in ns - @test state(n, SiteType("Electron")) == x - @test state(n, SiteType("tJ")) == x[1:3] - end - end - for (ns, x) in ( - (("n↑", "Nup"), [0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]), - (("n↓", "Ndn"), [0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]), - (("n↑↓", "Nupdn"), [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1]), - (("ntot", "Ntot"), [0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2]), - (("c↑", "Cup"), [0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]), - (("Sx",), [0 0 0 0; 0 0 1/2 0; 0 1/2 0 0; 0 0 0 0]), - # TODO: Add more tests. - ) - for n in ns - @test op(n, SiteType("Electron")) == x - @test op(n, SiteType("tJ")) == x[1:3, 1:3] - end + for ss in (sites("Qudit", 4; dim = 3), sites("Qudit", 2:5; dim = 3)) + @test length(ss) == 4 + @test ss == map(Returns(Base.OneTo(3)), 1:4) + for s in ss + @test s === Base.OneTo(3) + end + end + + for ss in ( + sites(AbstractUnitRange{Int32}, "Qudit", 4; dim = 3), + sites(AbstractUnitRange{Int32}, "Qudit", 2:5; dim = 3), + ) + @test length(ss) == 4 + @test ss == map(Returns(Base.OneTo(3)), 1:4) + for s in ss + @test s === Base.OneTo(Int32(3)) + end + end + + ss = sites("Qudit", (1, 2, 3); dim = 3) + @test ss === (Base.OneTo(3), Base.OneTo(3), Base.OneTo(3)) end - end - @testset "Boson and spin" begin - for t in (SiteType("Qudit"; dim=3), SiteType("Boson"; dim=3), SiteType("S"; spin=1)) - @test length(t) == 3 - @test op("X", t) == op("X", 3) + @testset "Electron/tJ" begin + for (ns, x) in ( + (("0", "Emp"), [1, 0, 0, 0]), + (("↑", "Up"), [0, 1, 0, 0]), + (("↓", "Dn"), [0, 0, 1, 0]), + (("↑↓", "UpDn"), [0, 0, 0, 1]), + ) + for n in ns + @test state(n, SiteType("Electron")) == x + @test state(n, SiteType("tJ")) == x[1:3] + end + end + for (ns, x) in ( + (("n↑", "Nup"), [0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]), + (("n↓", "Ndn"), [0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]), + (("n↑↓", "Nupdn"), [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1]), + (("ntot", "Ntot"), [0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2]), + (("c↑", "Cup"), [0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]), + (("Sx",), [0 0 0 0; 0 0 1 / 2 0; 0 1 / 2 0 0; 0 0 0 0]), + # TODO: Add more tests. + ) + for n in ns + @test op(n, SiteType("Electron")) == x + @test op(n, SiteType("tJ")) == x[1:3, 1:3] + end + end end + @testset "Boson and spin" begin + for t in (SiteType("Qudit"; dim = 3), SiteType("Boson"; dim = 3), SiteType("S"; spin = 1)) + @test length(t) == 3 + @test op("X", t) == op("X", 3) + end - for t in (SiteType("S=1"), SiteType("SpinOne")) - for (ns, x) in ( - (("Z+", "↑", "Up"), [1, 0, 0]), - (("Z0", "0"), [0, 1, 0]), - (("Z-", "↓", "Dn"), [0, 0, 1]), - (("X+",), [1 / 2, 1 / sqrt(2), 1 / 2]), - (("X0",), [-1 / sqrt(2), 0, 1 / sqrt(2)]), - (("X-",), [1 / 2, -1 / sqrt(2), 1 / 2]), - (("Y+",), [-1 / 2, -im / sqrt(2), 1 / 2]), - (("Y0",), [1 / sqrt(2), 0, 1 / sqrt(2)]), - (("Y-",), [-1 / 2, im / sqrt(2), 1 / 2]), - ) - for n in ns - @test state(n, t) ≈ x + for t in (SiteType("S=1"), SiteType("SpinOne")) + for (ns, x) in ( + (("Z+", "↑", "Up"), [1, 0, 0]), + (("Z0", "0"), [0, 1, 0]), + (("Z-", "↓", "Dn"), [0, 0, 1]), + (("X+",), [1 / 2, 1 / sqrt(2), 1 / 2]), + (("X0",), [-1 / sqrt(2), 0, 1 / sqrt(2)]), + (("X-",), [1 / 2, -1 / sqrt(2), 1 / 2]), + (("Y+",), [-1 / 2, -im / sqrt(2), 1 / 2]), + (("Y0",), [1 / sqrt(2), 0, 1 / sqrt(2)]), + (("Y-",), [-1 / 2, im / sqrt(2), 1 / 2]), + ) + for n in ns + @test state(n, t) ≈ x + end + end end - end - end - # Check they are eigenvalues - (λ⁺, λ⁰, λ⁻) = (2, 0, -2) - for t in (SiteType("S=1"), SiteType("SpinOne")) - for n in ("X", "Y", "Z") - @test op(n, t) * state("$(n)+", t) ≈ λ⁺ * state("$(n)+", t) atol = eps() - @test op(n, t) * state("$(n)0", t) ≈ λ⁰ * state("$(n)0", t) atol = eps() - @test op(n, t) * state("$(n)-", t) ≈ λ⁻ * state("$(n)-", t) atol = eps() - end + # Check they are eigenvalues + (λ⁺, λ⁰, λ⁻) = (2, 0, -2) + for t in (SiteType("S=1"), SiteType("SpinOne")) + for n in ("X", "Y", "Z") + @test op(n, t) * state("$(n)+", t) ≈ λ⁺ * state("$(n)+", t) atol = eps() + @test op(n, t) * state("$(n)0", t) ≈ λ⁰ * state("$(n)0", t) atol = eps() + @test op(n, t) * state("$(n)-", t) ≈ λ⁻ * state("$(n)-", t) atol = eps() + end + end end - end end diff --git a/test/test_exports.jl b/test/test_exports.jl index 92a25fe..f296ab6 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -1,17 +1,17 @@ using QuantumOperatorDefinitions: QuantumOperatorDefinitions using Test: @test, @testset @testset "Test exports" begin - exports = [ - :QuantumOperatorDefinitions, - Symbol("@OpName_str"), - Symbol("@SiteType_str"), - Symbol("@StateName_str"), - :OpName, - :SiteType, - :StateName, - :⊗, - :op, - :state, - ] - @test issetequal(names(QuantumOperatorDefinitions), exports) + exports = [ + :QuantumOperatorDefinitions, + Symbol("@OpName_str"), + Symbol("@SiteType_str"), + Symbol("@StateName_str"), + :OpName, + :SiteType, + :StateName, + :⊗, + :op, + :state, + ] + @test issetequal(names(QuantumOperatorDefinitions), exports) end diff --git a/test/test_gradedarraysext.jl b/test/test_gradedarraysext.jl index 20dd60d..f95bbd9 100644 --- a/test/test_gradedarraysext.jl +++ b/test/test_gradedarraysext.jl @@ -7,94 +7,94 @@ using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state using Test: @test, @testset @testset "GradedArraysExt" begin - t = SiteType("S=1/2"; gradings=("Sz",)) - r = AbstractUnitRange(t) - @test r == 1:2 - @test sectors(r) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r) == [1, 1] + t = SiteType("S=1/2"; gradings = ("Sz",)) + r = AbstractUnitRange(t) + @test r == 1:2 + @test sectors(r) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r) == [1, 1] - t = SiteType("Electron"; gradings=("Nf", "Sz")) - r = AbstractUnitRange(t) - @test r == 1:4 - @test sectors(r) == [ - SectorProduct((; Nf=U1(0), Sz=U1(0))), - SectorProduct((; Nf=U1(1), Sz=U1(1))), - SectorProduct((; Nf=U1(1), Sz=U1(-1))), - SectorProduct((; Nf=U1(2), Sz=U1(0))), - ] - @test blocklengths(r) == [1, 1, 1, 1] + t = SiteType("Electron"; gradings = ("Nf", "Sz")) + r = AbstractUnitRange(t) + @test r == 1:4 + @test sectors(r) == [ + SectorProduct((; Nf = U1(0), Sz = U1(0))), + SectorProduct((; Nf = U1(1), Sz = U1(1))), + SectorProduct((; Nf = U1(1), Sz = U1(-1))), + SectorProduct((; Nf = U1(2), Sz = U1(0))), + ] + @test blocklengths(r) == [1, 1, 1, 1] - t = SiteType("Electron"; gradings=("Nf" => "NfA", "Sz" => "SzA")) - r = AbstractUnitRange(t) - @test r == 1:4 - @test sectors(r) == [ - SectorProduct((; NfA=U1(0), SzA=U1(0))), - SectorProduct((; NfA=U1(1), SzA=U1(1))), - SectorProduct((; NfA=U1(1), SzA=U1(-1))), - SectorProduct((; NfA=U1(2), SzA=U1(0))), - ] - @test blocklengths(r) == [1, 1, 1, 1] + t = SiteType("Electron"; gradings = ("Nf" => "NfA", "Sz" => "SzA")) + r = AbstractUnitRange(t) + @test r == 1:4 + @test sectors(r) == [ + SectorProduct((; NfA = U1(0), SzA = U1(0))), + SectorProduct((; NfA = U1(1), SzA = U1(1))), + SectorProduct((; NfA = U1(1), SzA = U1(-1))), + SectorProduct((; NfA = U1(2), SzA = U1(0))), + ] + @test blocklengths(r) == [1, 1, 1, 1] - t = SiteType("Electron"; gradings=("NfParity", "Sz")) - r = AbstractUnitRange(t) - @test r == 1:4 - @test sectors(r) == [ - SectorProduct((; NfParity=Z{2}(0), Sz=U1(0))), - SectorProduct((; NfParity=Z{2}(1), Sz=U1(1))), - SectorProduct((; NfParity=Z{2}(1), Sz=U1(-1))), - SectorProduct((; NfParity=Z{2}(0), Sz=U1(0))), - ] - @test blocklengths(r) == [1, 1, 1, 1] + t = SiteType("Electron"; gradings = ("NfParity", "Sz")) + r = AbstractUnitRange(t) + @test r == 1:4 + @test sectors(r) == [ + SectorProduct((; NfParity = Z{2}(0), Sz = U1(0))), + SectorProduct((; NfParity = Z{2}(1), Sz = U1(1))), + SectorProduct((; NfParity = Z{2}(1), Sz = U1(-1))), + SectorProduct((; NfParity = Z{2}(0), Sz = U1(0))), + ] + @test blocklengths(r) == [1, 1, 1, 1] - t = SiteType("S=1/2"; gradings=("Sz",)) - (r1, r2) = axes(OpName("σ⁺"), (t,)) - @test sectors(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r1) == [1, 1] - @test sectors(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r2) == [1, 1] + t = SiteType("S=1/2"; gradings = ("Sz",)) + (r1, r2) = axes(OpName("σ⁺"), (t,)) + @test sectors(r1) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r1) == [1, 1] + @test sectors(r2) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r2) == [1, 1] - t = SiteType("S=1/2"; gradings=("Sz",)) - @test state("0", t) == [1, 0] + t = SiteType("S=1/2"; gradings = ("Sz",)) + @test state("0", t) == [1, 0] - # Force conversion to `Vector{Float64}` before conversion, - # since there is a bug slicing `BitVector` by `GradedOneTo`. - t = SiteType("S=1/2"; gradings=("Sz",)) - a = AbstractArray(2.0 * StateName("0"), t) - @test a == [2, 0] - @test a isa BlockSparseArray - (r1,) = axes(a) - @test sectors(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r1) == [1, 1] + # Force conversion to `Vector{Float64}` before conversion, + # since there is a bug slicing `BitVector` by `GradedOneTo`. + t = SiteType("S=1/2"; gradings = ("Sz",)) + a = AbstractArray(2.0 * StateName("0"), t) + @test a == [2, 0] + @test a isa BlockSparseArray + (r1,) = axes(a) + @test sectors(r1) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r1) == [1, 1] - t = SiteType("S=1/2"; gradings=("Sz",)) - a = op("σ⁺", t) - @test a == [0 2; 0 0] - @test a isa BlockSparseArray - (r1, r2) = axes(a) - @test sectors(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r1) == [1, 1] - @test sectors(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r2) == [1, 1] + t = SiteType("S=1/2"; gradings = ("Sz",)) + a = op("σ⁺", t) + @test a == [0 2; 0 0] + @test a isa BlockSparseArray + (r1, r2) = axes(a) + @test sectors(r1) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r1) == [1, 1] + @test sectors(r2) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r2) == [1, 1] end @testset "GradedArraysExt + ITensorBaseExt" begin - i = Index(SiteType("S=1/2"; gradings=("Sz",))) - @test gettag(i, "sitetype") == "S=1/2" - # TODO: Test without denaming. - @test dename(i) == 1:2 - @test sectors(dename(i)) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(dename(i)) == [1, 1] + i = Index(SiteType("S=1/2"; gradings = ("Sz",))) + @test gettag(i, "sitetype") == "S=1/2" + # TODO: Test without denaming. + @test dename(i) == 1:2 + @test sectors(dename(i)) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(dename(i)) == [1, 1] - i′ = prime(i) - a = op("σ⁺", i) - @test a == ITensor([0 2; 0 0], (i′, dual(i))) - @test all(isdual.(axes(a)) .== (false, true)) - a′ = dename(a) - @test a′ isa BlockSparseArray - # TODO: Test these without denaming `a`. - (r1, r2) = axes(a′) - @test sectors(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r1) == [1, 1] - @test sectors(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))] - @test blocklengths(r2) == [1, 1] + i′ = prime(i) + a = op("σ⁺", i) + @test a == ITensor([0 2; 0 0], (i′, dual(i))) + @test all(isdual.(axes(a)) .== (false, true)) + a′ = dename(a) + @test a′ isa BlockSparseArray + # TODO: Test these without denaming `a`. + (r1, r2) = axes(a′) + @test sectors(r1) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r1) == [1, 1] + @test sectors(r2) == [SectorProduct((; Sz = U1(0))), SectorProduct((; Sz = U1(1)))] + @test blocklengths(r2) == [1, 1] end diff --git a/test/test_itensorbaseext.jl b/test/test_itensorbaseext.jl index 3407ac3..87ede22 100644 --- a/test/test_itensorbaseext.jl +++ b/test/test_itensorbaseext.jl @@ -5,88 +5,88 @@ using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, site, sites, using Test: @test, @testset @testset "ITensorBaseExt" begin - i = Index(SiteType("S=1/2")) - @test gettag(i, "sitetype") == "S=1/2" + i = Index(SiteType("S=1/2")) + @test gettag(i, "sitetype") == "S=1/2" - i = Index(2) - a = state("0", i) - @test a == ITensor(state("0", 2), i) - @test a == state("0", (i,)) - @test a == ITensor(StateName("0"), i) - @test a == ITensor(StateName("0"), (i,)) + i = Index(2) + a = state("0", i) + @test a == ITensor(state("0", 2), i) + @test a == state("0", (i,)) + @test a == ITensor(StateName("0"), i) + @test a == ITensor(StateName("0"), (i,)) - i = settag(Index(2), "sitetype", "S=1/2") - a = state("X+", i) - @test a == ITensor(state("X+", SiteType("S=1/2")), i) + i = settag(Index(2), "sitetype", "S=1/2") + a = state("X+", i) + @test a == ITensor(state("X+", SiteType("S=1/2")), i) - i = Index(2) - a = op("X", i) - @test a == ITensor(op("X", 2), (prime(i), i)) - @test a == op("X", (i,)) - @test a == ITensor(OpName("X"), i) - @test a == ITensor(OpName("X"), (i,)) + i = Index(2) + a = op("X", i) + @test a == ITensor(op("X", 2), (prime(i), i)) + @test a == op("X", (i,)) + @test a == ITensor(OpName("X"), i) + @test a == ITensor(OpName("X"), (i,)) - i1 = Index(2) - i2 = Index(2) - i1′ = prime(i1) - i2′ = prime(i2) - a = ITensor(OpName("CX"), (i1, i2)) - @test a == ITensor(op("CX", (2, 2)), (i1′, i2′, i1, i2)) - @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) - @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) - @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) - @test a[i1′[2], i2′, i1[2], i2] == op("X", i2) + i1 = Index(2) + i2 = Index(2) + i1′ = prime(i1) + i2′ = prime(i2) + a = ITensor(OpName("CX"), (i1, i2)) + @test a == ITensor(op("CX", (2, 2)), (i1′, i2′, i1, i2)) + @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) + @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) + @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) + @test a[i1′[2], i2′, i1[2], i2] == op("X", i2) - i1 = Index(3) - i2 = Index(2) - i1′ = prime(i1) - i2′ = prime(i2) - a = ITensor(OpName("CX"), (i1, i2)) - @test a == ITensor(op("CX", (3, 2)), (i1′, i2′, i1, i2)) - @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) - @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) - @test a[i1′[3], i2′, i1[1], i2] == zeros(i2′, i2) - @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) - @test a[i1′[2], i2′, i1[2], i2] == op("I", i2) - @test a[i1′[3], i2′, i1[2], i2] == zeros(i2′, i2) - @test a[i1′[1], i2′, i1[3], i2] == zeros(i2′, i2) - @test a[i1′[2], i2′, i1[3], i2] == zeros(i2′, i2) - @test a[i1′[3], i2′, i1[3], i2] == op("X", i2) + i1 = Index(3) + i2 = Index(2) + i1′ = prime(i1) + i2′ = prime(i2) + a = ITensor(OpName("CX"), (i1, i2)) + @test a == ITensor(op("CX", (3, 2)), (i1′, i2′, i1, i2)) + @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) + @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) + @test a[i1′[3], i2′, i1[1], i2] == zeros(i2′, i2) + @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) + @test a[i1′[2], i2′, i1[2], i2] == op("I", i2) + @test a[i1′[3], i2′, i1[2], i2] == zeros(i2′, i2) + @test a[i1′[1], i2′, i1[3], i2] == zeros(i2′, i2) + @test a[i1′[2], i2′, i1[3], i2] == zeros(i2′, i2) + @test a[i1′[3], i2′, i1[3], i2] == op("X", i2) - i1 = Index(2) - i2 = Index(3) - i1′ = prime(i1) - i2′ = prime(i2) - a = ITensor(OpName("CX"), (i1, i2)) - @test a == ITensor(op("CX", (2, 3)), (i1′, i2′, i1, i2)) - @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) - @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) - @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) - @test a[i1′[2], i2′, i1[2], i2] == op("X", i2) + i1 = Index(2) + i2 = Index(3) + i1′ = prime(i1) + i2′ = prime(i2) + a = ITensor(OpName("CX"), (i1, i2)) + @test a == ITensor(op("CX", (2, 3)), (i1′, i2′, i1, i2)) + @test a[i1′[1], i2′, i1[1], i2] == op("I", i2) + @test a[i1′[2], i2′, i1[1], i2] == zeros(i2′, i2) + @test a[i1′[1], i2′, i1[2], i2] == zeros(i2′, i2) + @test a[i1′[2], i2′, i1[2], i2] == op("X", i2) - i = site(Index, "Qudit"; dim=3) - @test dename(i) == Base.OneTo(3) - @test gettag(i, "sitetype") == "Qudit" - @test !hastag(i, "site") + i = site(Index, "Qudit"; dim = 3) + @test dename(i) == Base.OneTo(3) + @test gettag(i, "sitetype") == "Qudit" + @test !hastag(i, "site") - for is in ( - sites(Index, "Qudit", 3; dim=3), - sites(Index, "Qudit", 1:3; dim=3), - sites(Index, "Qudit", (1, 2, 3); dim=3), - ) - @test length(is) == 3 - for (pos, i) in pairs(is) - @test dename(i) == Base.OneTo(3) - @test gettag(i, "sitetype") == "Qudit" - @test gettag(i, "site") == "$pos" + for is in ( + sites(Index, "Qudit", 3; dim = 3), + sites(Index, "Qudit", 1:3; dim = 3), + sites(Index, "Qudit", (1, 2, 3); dim = 3), + ) + @test length(is) == 3 + for (pos, i) in pairs(is) + @test dename(i) == Base.OneTo(3) + @test gettag(i, "sitetype") == "Qudit" + @test gettag(i, "site") == "$pos" + end end - end - is = sites(Index, "Qudit", 2:4; dim=3) - @test dename(is[1]) == Base.OneTo(3) - @test gettag(is[1], "site") == "2" - @test dename(is[2]) == Base.OneTo(3) - @test gettag(is[2], "site") == "3" - @test dename(is[3]) == Base.OneTo(3) - @test gettag(is[3], "site") == "4" + is = sites(Index, "Qudit", 2:4; dim = 3) + @test dename(is[1]) == Base.OneTo(3) + @test gettag(is[1], "site") == "2" + @test dename(is[2]) == Base.OneTo(3) + @test gettag(is[2], "site") == "3" + @test dename(is[3]) == Base.OneTo(3) + @test gettag(is[3], "site") == "4" end