From a0aefe6671258b3ba9bd14b102eb56ceea2ee86c Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 30 Jan 2025 14:05:08 -0500
Subject: [PATCH 01/12] Allow custom axes when constructing operators, ITensor
 extension

---
 Project.toml                                  |  4 +-
 ...uantumOperatorDefinitionsITensorBaseExt.jl | 10 +++-
 src/op.jl                                     | 48 +++++++++++++------
 src/sitetype.jl                               |  9 +++-
 src/state.jl                                  |  3 +-
 5 files changed, 55 insertions(+), 19 deletions(-)

diff --git a/Project.toml b/Project.toml
index 85d6394..f2f65f2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,12 +9,14 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [weakdeps]
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 
 [extensions]
-QuantumOperatorDefinitionsITensorBaseExt = "ITensorBase"
+QuantumOperatorDefinitionsITensorBaseExt = ["ITensorBase", "NamedDimsArrays"]
 
 [compat]
 ITensorBase = "0.1.10"
 LinearAlgebra = "1.10"
+NamedDimsArrays = "0.4.0"
 Random = "1.10"
 julia = "1.10"
diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index 280b531..df5f1c3 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -1,11 +1,18 @@
 module QuantumOperatorDefinitionsITensorBaseExt
 
 using ITensorBase: ITensor, Index, dag, gettag, prime
+using NamedDimsArrays: dename
 using QuantumOperatorDefinitions:
   QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
 
 function QuantumOperatorDefinitions.SiteType(r::Index)
-  return SiteType(gettag(r, "sitetype", "Qudit"); dim=Int(length(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 (and 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 QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index)
@@ -14,6 +21,7 @@ end
 
 function Base.AbstractArray(n::OpName, r::Index)
   # TODO: Define this with mapped dimnames.
+  # Generalize beyond prime levels with codomain and domain indices.
   return ITensor(AbstractArray(n, SiteType(r)), (prime(r), dag(r)))
 end
 
diff --git a/src/op.jl b/src/op.jl
index c7fbf43..2e2d2e8 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -8,8 +8,8 @@ struct OpName{Name,Params}
 end
 name(::OpName{Name}) where {Name} = Name
 params(n::OpName) = getfield(n, :params)
-
 Base.getproperty(n::OpName, name::Symbol) = getfield(params(n), name)
+Base.get(t::OpName, name::Symbol, default) = get(params(t), name, default)
 
 OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...))
 
@@ -54,7 +54,9 @@ end
 # Generic to `StateName` or `OpName`.
 const StateOrOpName = Union{StateName,OpName}
 alias(n::StateOrOpName) = n
-function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Integer...)
+function (arrtype::Type{<:AbstractArray})(
+  n::StateOrOpName, domain::Union{Integer,AbstractUnitRange}...
+)
   return arrtype(n, domain)
 end
 (arrtype::Type{<:AbstractArray})(n::StateOrOpName, ts::SiteType...) = arrtype(n, ts)
@@ -87,32 +89,50 @@ function nsites(n::StateOrOpName)
   return nsites(n′)
 end
 
+function array(a::AbstractArray, ax::Tuple{Vararg{AbstractUnitRange}})
+  return a[ax...]
+end
+
 function op_convert(
   arrtype::Type{<:AbstractArray{<:Any,N}},
-  domain::Tuple{Vararg{Integer}},
+  domain::Tuple{Vararg{AbstractUnitRange}},
   a::AbstractArray{<:Any,N},
 ) where {N}
-  # TODO: Check the dimensions.
-  return convert(arrtype, a)
+  ax = (domain..., domain...)
+  a′ = array(a, ax)
+  return convert(arrtype, a′)
 end
 function op_convert(
-  arrtype::Type{<:AbstractArray}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
+  arrtype::Type{<:AbstractArray}, domain::Tuple{Vararg{AbstractUnitRange}}, a::AbstractArray
 )
-  # TODO: Check the dimensions.
-  return convert(arrtype, a)
+  ax = (domain..., domain...)
+  a′ = array(a, ax)
+  return convert(arrtype, a′)
 end
 function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}}, domain::Tuple{Vararg{Integer}}, a::AbstractArray
+  arrtype::Type{<:AbstractArray{<:Any,N}},
+  domain::Tuple{Vararg{AbstractUnitRange}},
+  a::AbstractArray,
 ) where {N}
-  size = (domain..., domain...)
-  @assert length(size) == N
-  return convert(arrtype, reshape(a, size))
+  ax = (domain..., domain...)
+  @assert length(ax) == N
+  a′ = reshape(a, length.(ax))
+  a′′ = array(a′, ax)
+  return convert(arrtype, a′′)
 end
 function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{SiteType}})
-  return op_convert(arrtype, length.(domain), n(domain...))
+  domain′ = AbstractUnitRange.(domain)
+  return op_convert(arrtype, domain′, n(domain...))
+end
+
+function (arrtype::Type{<:AbstractArray})(
+  n::OpName, domain::Tuple{Vararg{AbstractUnitRange}}
+)
+  # TODO: Make `(::OpName)(domain...)` constructor process more general inputs.
+  return op_convert(arrtype, domain, n(Int.(length.(domain))...))
 end
 function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{Integer}})
-  return op_convert(arrtype, domain, n(Int.(domain)...))
+  return arrtype(n, Base.oneto.(domain))
 end
 
 function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 0172f61..24ec663 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -4,14 +4,15 @@ struct SiteType{T,Params}
     return new{N,typeof(params)}(params)
   end
 end
+value(::SiteType{T}) where {T} = T
 params(t::SiteType) = getfield(t, :params)
 Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name)
+Base.get(t::SiteType, name::Symbol, default) = get(params(t), name, default)
 
 SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...))
 
 SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...)
 SiteType(i::Integer; kwargs...) = SiteType{Symbol(i)}(; kwargs...)
-value(::SiteType{T}) where {T} = T
 macro SiteType_str(s)
   return SiteType{Symbol(s)}
 end
@@ -26,7 +27,11 @@ function Base.length(t::SiteType)
   end
   return length(t′)
 end
-Base.AbstractUnitRange(t::SiteType) = Base.OneTo(length(t))
+function Base.AbstractUnitRange(t::SiteType)
+  # This logic allows specifying a range with extra properties,
+  # like ones with symmetry sectors.
+  return get(t, :range, Base.OneTo(length(t)))
+end
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]
 Base.axes(t::SiteType) = (AbstractUnitRange(t),)
diff --git a/src/state.jl b/src/state.jl
index d9904a5..65b6591 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -6,14 +6,15 @@ struct StateName{Name,Params}
     return new{N,typeof(params)}(params)
   end
 end
+name(::StateName{N}) where {N} = N
 params(n::StateName) = getfield(n, :params)
 Base.getproperty(n::StateName, name::Symbol) = getfield(params(n), name)
+Base.get(t::StateName, name::Symbol, default) = get(params(t), name, default)
 
 StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...))
 
 StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...)
 StateName(s::Symbol; kwargs...) = StateName{s}(; kwargs...)
-name(::StateName{N}) where {N} = N
 macro StateName_str(s)
   return StateName{Symbol(s)}
 end

From 08db5010dd172a06b5a3433bee9339dae719be99 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 30 Jan 2025 17:28:30 -0500
Subject: [PATCH 02/12] Output higher dimensional arrays for multi-site
 operators

---
 README.md           |   3 --
 examples/README.jl  |   3 --
 src/op.jl           |  60 +++++++++++--------------
 src/state.jl        |  11 +----
 test/test_basics.jl | 104 +++++++++++++++++++++++++-------------------
 5 files changed, 86 insertions(+), 95 deletions(-)

diff --git a/README.md b/README.md
index 978d575..480173d 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,6 @@ julia> Pkg.add("QuantumOperatorDefinitions")
 
 ````julia
 using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state
-using LinearAlgebra: Diagonal
 using SparseArrays: SparseMatrixCSC, SparseVector
 using Test: @test
 
@@ -64,8 +63,6 @@ using Test: @test
 @test op("Y") == [0 -im; im 0]
 @test op("Z") == [1 0; 0 -1]
 
-@test op("Z") isa Diagonal
-
 @test op(Float32, "X") == [0 1; 1 0]
 @test eltype(op(Float32, "X")) === Float32
 @test op(SparseMatrixCSC, "X") == [0 1; 1 0]
diff --git a/examples/README.jl b/examples/README.jl
index 5b8b5b1..f4f246d 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -38,7 +38,6 @@ julia> Pkg.add("QuantumOperatorDefinitions")
 # ## Examples
 
 using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state
-using LinearAlgebra: Diagonal
 using SparseArrays: SparseMatrixCSC, SparseVector
 using Test: @test
 
@@ -69,8 +68,6 @@ using Test: @test
 @test op("Y") == [0 -im; im 0]
 @test op("Z") == [1 0; 0 -1]
 
-@test op("Z") isa Diagonal
-
 @test op(Float32, "X") == [0 1; 1 0]
 @test eltype(op(Float32, "X")) === Float32
 @test op(SparseMatrixCSC, "X") == [0 1; 1 0]
diff --git a/src/op.jl b/src/op.jl
index 2e2d2e8..7d6d571 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -59,6 +59,9 @@ function (arrtype::Type{<:AbstractArray})(
 )
   return arrtype(n, domain)
 end
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{Integer}})
+  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?
@@ -89,50 +92,37 @@ function nsites(n::StateOrOpName)
   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...]
 end
 
-function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}},
-  domain::Tuple{Vararg{AbstractUnitRange}},
-  a::AbstractArray{<:Any,N},
-) where {N}
-  ax = (domain..., domain...)
-  a′ = array(a, ax)
-  return convert(arrtype, a′)
-end
-function op_convert(
-  arrtype::Type{<:AbstractArray}, domain::Tuple{Vararg{AbstractUnitRange}}, a::AbstractArray
-)
-  ax = (domain..., domain...)
-  a′ = array(a, ax)
-  return convert(arrtype, a′)
+function state_or_op_axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}})
+  return (domain..., domain...)
 end
-function op_convert(
-  arrtype::Type{<:AbstractArray{<:Any,N}},
+
+function state_or_op_convert(
+  n::StateOrOpName,
+  arrtype::Type{<:AbstractArray},
   domain::Tuple{Vararg{AbstractUnitRange}},
   a::AbstractArray,
-) where {N}
-  ax = (domain..., domain...)
-  @assert length(ax) == N
+)
+  ax = state_or_op_axes(n, domain)
   a′ = reshape(a, length.(ax))
   a′′ = array(a′, ax)
   return convert(arrtype, a′′)
 end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{SiteType}})
+
+function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{SiteType}})
   domain′ = AbstractUnitRange.(domain)
-  return op_convert(arrtype, domain′, n(domain...))
+  return state_or_op_convert(n, arrtype, domain′, n(domain...))
 end
-
 function (arrtype::Type{<:AbstractArray})(
-  n::OpName, domain::Tuple{Vararg{AbstractUnitRange}}
+  n::StateOrOpName, domain::Tuple{Vararg{AbstractUnitRange}}
 )
   # TODO: Make `(::OpName)(domain...)` constructor process more general inputs.
-  return op_convert(arrtype, domain, n(Int.(length.(domain))...))
-end
-function (arrtype::Type{<:AbstractArray})(n::OpName, domain::Tuple{Vararg{Integer}})
-  return arrtype(n, Base.oneto.(domain))
+  return state_or_op_convert(n, arrtype, domain, n(Int.(length.(domain))...))
 end
 
 function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
@@ -498,10 +488,10 @@ function (n::OpName"Controlled")(domain...)
   d_control = prod(to_dim.(domain[1:nc]))
   return cat(I(d_control), n.arg(domain[(nc + 1):end]...); dims=(1, 2))
 end
-@op_alias "CNOT" "Controlled" op = OpName"X"()
-@op_alias "CX" "Controlled" op = OpName"X"()
-@op_alias "CY" "Controlled" op = OpName"Y"()
-@op_alias "CZ" "Controlled" op = OpName"Z"()
+@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)...))
 end
@@ -524,17 +514,17 @@ function alias(::OpName"CRn")
 end
 @op_alias "CRn̂" "CRn"
 
-@op_alias "CCNOT" "Controlled" ncontrol = 2 op = OpName"X"()
+@op_alias "CCNOT" "Controlled" ncontrol = 2 arg = OpName"X"()
 @op_alias "Toffoli" "CCNOT"
 @op_alias "CCX" "CCNOT"
 @op_alias "TOFF" "CCNOT"
 
-@op_alias "CSWAP" "Controlled" ncontrol = 2 op = OpName"SWAP"()
+@op_alias "CSWAP" "Controlled" ncontrol = 2 arg = OpName"SWAP"()
 @op_alias "Fredkin" "CSWAP"
 @op_alias "CSwap" "CSWAP"
 @op_alias "CS" "CSWAP"
 
-@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"()
+@op_alias "CCCNOT" "Controlled" ncontrol = 3 arg = OpName"X"()
 
 ## # 1-qudit rotation around generic axis n̂.
 ## # exp(-im * α / 2 * n̂ ⋅ σ⃗)
diff --git a/src/state.jl b/src/state.jl
index 65b6591..55c5155 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -32,15 +32,8 @@ macro state_alias(name1, name2, params...)
   return state_alias_expr(name1, name2)
 end
 
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain::Tuple{Vararg{SiteType}})
-  # TODO: Define `state_convert` to handle reshaping multisite states
-  # to higher order arrays.
-  return convert(arrtype, n(domain...))
-end
-function (arrtype::Type{<:AbstractArray})(n::StateName, domain::Tuple{Vararg{Integer}})
-  # TODO: Define `state_convert` to handle reshaping multisite states
-  # to higher order arrays.
-  return convert(arrtype, n(Int.(domain)...))
+function state_or_op_axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
+  return domain
 end
 
 # This compiles operator expressions, such as:
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 53f2595..31d0c96 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -33,22 +33,30 @@ const elts = (real_elts..., complex_elts...)
         [0 -im; im 0],
         [1 0; 0 -1],
         [0 0; 0 1],
-        [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1],
-        [1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1],
-        (_, θ) -> [
-          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)
-        ],
-        (_, θ) -> [
-          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)
-        ],
-        (_, θ) ->
+        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],
@@ -60,31 +68,37 @@ const elts = (real_elts..., complex_elts...)
         √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],
-        [
-          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
-        ],
-        [
-          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
-        ],
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
-        (O, θ) -> exp(-im * (θ / 2) * kron(O, O)),
+        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],
@@ -114,9 +128,9 @@ const elts = (real_elts..., complex_elts...)
         (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("√SWAP"), 2, complex_elts, √SWAPmat),
         (OpName("iSWAP"), 2, complex_elts, iSWAPmat),
-        (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)),
@@ -128,7 +142,7 @@ const elts = (real_elts..., complex_elts...)
         (OpName("StandardBasis"; index=(1, 2)), 1, elts, StandardBasis12mat),
       )
         @test nsites(o) == nbits
-        for arraytype in (AbstractArray, AbstractMatrix, Array, Matrix)
+        for arraytype in (AbstractArray, Array)
           for elt in elts
             ts = ntuple(Returns(t), nbits)
             lens = ntuple(Returns(len), nbits)
@@ -149,7 +163,7 @@ const elts = (real_elts..., complex_elts...)
     @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))") ==
-      exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z"))))
+      reshape(exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))), (2, 2, 2, 2))
     @test op("Ry{θ=π/2}") == op("Ry"; θ=π / 2)
     # Awkward parsing corner cases.
     @test op("S+") == Matrix(OpName("S+"))
@@ -187,7 +201,7 @@ const elts = (real_elts..., complex_elts...)
     @test state("2", 3) == [0, 0, 1]
 
     @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+")
-    @test state("|0⟩ ⊗ |+⟩") == kron(state("0"), state("+"))
+    @test state("|0⟩ ⊗ |+⟩") == reshape(kron(state("0"), state("+")), (2, 2))
   end
   @testset "Electron/tJ" begin
     for (ns, x) in (

From 833cc47842d741a8fcffe1d9da1485cf466f2b7f Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 30 Jan 2025 21:38:27 -0500
Subject: [PATCH 03/12] Support for specifying symmetries of axes

---
 Project.toml                                  |  9 ++
 ...uantumOperatorDefinitionsITensorBaseExt.jl |  6 +-
 ...umOperatorDefinitionsSymmetrySectorsExt.jl | 86 +++++++++++++++++++
 src/sitetype.jl                               | 56 ++++++++++--
 4 files changed, 150 insertions(+), 7 deletions(-)
 create mode 100644 ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl

diff --git a/Project.toml b/Project.toml
index f2f65f2..f1caba2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,15 +8,24 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [weakdeps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
+LabelledNumbers = "f856a3a6-4152-4ec4-b2a7-02c1a55d7993"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
+SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e"
 
 [extensions]
 QuantumOperatorDefinitionsITensorBaseExt = ["ITensorBase", "NamedDimsArrays"]
+QuantumOperatorDefinitionsSymmetrySectorsExt = ["BlockArrays", "GradedUnitRanges", "LabelledNumbers", "SymmetrySectors"]
 
 [compat]
+BlockArrays = "1.3.0"
+GradedUnitRanges = "0.1.2"
 ITensorBase = "0.1.10"
+LabelledNumbers = "0.1.0"
 LinearAlgebra = "1.10"
 NamedDimsArrays = "0.4.0"
 Random = "1.10"
+SymmetrySectors = "0.1.3"
 julia = "1.10"
diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index df5f1c3..067b87e 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -1,6 +1,6 @@
 module QuantumOperatorDefinitionsITensorBaseExt
 
-using ITensorBase: ITensor, Index, dag, gettag, prime
+using ITensorBase: ITensorBase, ITensor, Index, dag, gettag, prime
 using NamedDimsArrays: dename
 using QuantumOperatorDefinitions:
   QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
@@ -15,6 +15,10 @@ function QuantumOperatorDefinitions.SiteType(r::Index)
   )
 end
 
+function ITensorBase.Index(t::SiteType; kwargs...)
+  return Index(AbstractUnitRange(t); kwargs...)
+end
+
 function QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index)
   return has_fermion_string(OpName(n), SiteType(r))
 end
diff --git a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
new file mode 100644
index 0000000..4d2b9e5
--- /dev/null
+++ b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
@@ -0,0 +1,86 @@
+module QuantumOperatorDefinitionsSymmetrySectorsExt
+
+using BlockArrays: blocklasts, blocklengths
+using GradedUnitRanges: GradedOneTo, gradedrange
+using LabelledNumbers: label, labelled, unlabel
+using QuantumOperatorDefinitions:
+  QuantumOperatorDefinitions, @SiteType_str, @SymmetryType_str, SiteType, SymmetryType, name
+using SymmetrySectors: ×, SectorProduct, U1, Z
+
+sortedunion(a, b) = sort(union(a, b))
+function QuantumOperatorDefinitions.combine_axes(a1::GradedOneTo, a2::GradedOneTo)
+  return gradedrange(
+    map(blocklengths(a1), blocklengths(a2)) do s1, s2
+      l1 = unlabel(s1)
+      l2 = unlabel(s2)
+      @assert l1 == l2
+      labelled(l1, label(s1) × label(s2))
+    end,
+  )
+end
+QuantumOperatorDefinitions.combine_axes(a::GradedOneTo, b::Base.OneTo) = a
+QuantumOperatorDefinitions.combine_axes(a::Base.OneTo, b::GradedOneTo) = b
+
+function Base.AbstractUnitRange(::SymmetryType"N", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; N=U1(i - 1))) => 1, 1:length(t)))
+end
+function Base.AbstractUnitRange(::SymmetryType"Sz", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; Sz=U1(i - 1))) => 1, 1:length(t)))
+end
+function Base.AbstractUnitRange(::SymmetryType"Sz↑", t::SiteType)
+  return AbstractUnitRange(SymmetryType"Sz"(), t)
+end
+function Base.AbstractUnitRange(::SymmetryType"Sz↓", t::SiteType)
+  return gradedrange(map(i -> SectorProduct((; Sz=U1(-(i - 1)))) => 1, 1:length(t)))
+end
+
+function sector(symmetrytype::SymmetryType, sec)
+  sectorname = Symbol(get(symmetrytype, :name, name(symmetrytype)))
+  return SectorProduct(NamedTuple{(sectorname,)}((sec,)))
+end
+
+function Base.AbstractUnitRange(s::SymmetryType"Nf", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+# TODO: Write in terms of `SymmetryType"Nf"` definition.
+function Base.AbstractUnitRange(s::SymmetryType"NfParity", t::SiteType"Fermion")
+  return gradedrange([sector(s, Z{2}(0)) => 1, sector(s, Z{2}(1)) => 1])
+end
+function Base.AbstractUnitRange(s::SymmetryType"Sz", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+function Base.AbstractUnitRange(s::SymmetryType"Sz↑", t::SiteType"Fermion")
+  return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
+end
+function Base.AbstractUnitRange(s::SymmetryType"Sz↓", t::SiteType"Fermion")
+  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::SymmetryType"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,
+  ])
+end
+# TODO: Write in terms of `SymmetryType"Nf"` definition.
+function Base.AbstractUnitRange(s::SymmetryType"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,
+  ])
+end
+function Base.AbstractUnitRange(s::SymmetryType"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,
+  ])
+end
+
+end
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 24ec663..39042c9 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -4,33 +4,77 @@ struct SiteType{T,Params}
     return new{N,typeof(params)}(params)
   end
 end
-value(::SiteType{T}) where {T} = T
+name(::SiteType{T}) where {T} = T
 params(t::SiteType) = getfield(t, :params)
 Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name)
 Base.get(t::SiteType, name::Symbol, default) = get(params(t), name, default)
-
+Base.haskey(t::SiteType, name::Symbol) = haskey(params(t), name)
 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{Symbol(s)}
+  return :(SiteType{$(Expr(:quote, Symbol(s)))})
 end
 
 alias(t::SiteType) = t
 alias(i::Integer) = i
 
+# Like `Base.Broadcast.axistype` (https://github.com/JuliaLang/julia/blob/v1.11.3/base/broadcast.jl#L536-L538)
+# and `BlockArrays.combine_blockaxes` (https://github.com/JuliaArrays/BlockArrays.jl/blob/v1.3.0/src/blockbroadcast.jl#L37-L38).
+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)
+end
+combine_axes(a) = a
+combine_axes(a, b, rest...) = combine_axes(combine_axes(a, b), rest...)
+
+struct SymmetryType{T,Params}
+  params::Params
+  function SymmetryType{N}(params::NamedTuple) where {N}
+    return new{N,typeof(params)}(params)
+  end
+end
+name(::SymmetryType{T}) where {T} = T
+params(t::SymmetryType) = getfield(t, :params)
+Base.getproperty(t::SymmetryType, name::Symbol) = getfield(params(t), name)
+Base.get(t::SymmetryType, name::Symbol, default) = get(params(t), name, default)
+Base.haskey(t::SymmetryType, name::Symbol) = haskey(params(t), name)
+SymmetryType{N}(; kwargs...) where {N} = SymmetryType{N}((; kwargs...))
+SymmetryType(s::AbstractString; kwargs...) = SymmetryType{Symbol(s)}(; kwargs...)
+function SymmetryType(s::Pair{<:AbstractString,<:AbstractString}; kwargs...)
+  return SymmetryType(first(s); kwargs..., name=last(s))
+end
+function SymmetryType(s::Pair{<:AbstractString,<:NamedTuple}; kwargs...)
+  return SymmetryType(first(s); kwargs..., last(s)...)
+end
+macro SymmetryType_str(s)
+  return :(SymmetryType{$(Expr(:quote, Symbol(s)))})
+end
+
+function Base.AbstractUnitRange(symmetry::SymmetryType, t::SiteType)
+  return error("Not implemented.")
+end
+function Base.AbstractUnitRange(symmetry::SymmetryType"Trivial", t::SiteType)
+  return Base.OneTo(length(t))
+end
+
 function Base.length(t::SiteType)
   t′ = alias(t)
   if t == t′
-    return t.length
+    return t.dim
   end
   return length(t′)
 end
 function Base.AbstractUnitRange(t::SiteType)
   # This logic allows specifying a range with extra properties,
   # like ones with symmetry sectors.
-  return get(t, :range, Base.OneTo(length(t)))
+  haskey(t, :range) && return t.range
+  if haskey(t, :symmetries)
+    rs = map(symmetry -> AbstractUnitRange(SymmetryType(symmetry), t), t.symmetries)
+    return combine_axes(Base.OneTo(length(t)), rs...)
+  end
+  return Base.OneTo(length(t))
 end
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]

From e5bddb7c25c17c4e23a8fba9bdbf20248f732abf Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 30 Jan 2025 21:39:09 -0500
Subject: [PATCH 04/12] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f1caba2..5f28d05 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumOperatorDefinitions"
 uuid = "826dd319-6fd5-459a-a990-3a4f214664bf"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 44105b0755be1848a80ff0bc65f3cdd4909eef64 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 08:22:16 -0500
Subject: [PATCH 05/12] Change terminology from symmetry to grading

---
 ...umOperatorDefinitionsSymmetrySectorsExt.jl | 36 +++++++++---------
 src/sitetype.jl                               | 38 +++++++++----------
 2 files changed, 37 insertions(+), 37 deletions(-)

diff --git a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
index 4d2b9e5..9814efc 100644
--- a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
+++ b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
@@ -4,7 +4,7 @@ using BlockArrays: blocklasts, blocklengths
 using GradedUnitRanges: GradedOneTo, gradedrange
 using LabelledNumbers: label, labelled, unlabel
 using QuantumOperatorDefinitions:
-  QuantumOperatorDefinitions, @SiteType_str, @SymmetryType_str, SiteType, SymmetryType, name
+  QuantumOperatorDefinitions, @SiteType_str, @GradingType_str, SiteType, GradingType, name
 using SymmetrySectors: ×, SectorProduct, U1, Z
 
 sortedunion(a, b) = sort(union(a, b))
@@ -21,43 +21,43 @@ end
 QuantumOperatorDefinitions.combine_axes(a::GradedOneTo, b::Base.OneTo) = a
 QuantumOperatorDefinitions.combine_axes(a::Base.OneTo, b::GradedOneTo) = b
 
-function Base.AbstractUnitRange(::SymmetryType"N", t::SiteType)
+function Base.AbstractUnitRange(::GradingType"N", t::SiteType)
   return gradedrange(map(i -> SectorProduct((; N=U1(i - 1))) => 1, 1:length(t)))
 end
-function Base.AbstractUnitRange(::SymmetryType"Sz", t::SiteType)
+function Base.AbstractUnitRange(::GradingType"Sz", t::SiteType)
   return gradedrange(map(i -> SectorProduct((; Sz=U1(i - 1))) => 1, 1:length(t)))
 end
-function Base.AbstractUnitRange(::SymmetryType"Sz↑", t::SiteType)
-  return AbstractUnitRange(SymmetryType"Sz"(), t)
+function Base.AbstractUnitRange(::GradingType"Sz↑", t::SiteType)
+  return AbstractUnitRange(GradingType"Sz"(), t)
 end
-function Base.AbstractUnitRange(::SymmetryType"Sz↓", t::SiteType)
+function Base.AbstractUnitRange(::GradingType"Sz↓", t::SiteType)
   return gradedrange(map(i -> SectorProduct((; Sz=U1(-(i - 1)))) => 1, 1:length(t)))
 end
 
-function sector(symmetrytype::SymmetryType, sec)
-  sectorname = Symbol(get(symmetrytype, :name, name(symmetrytype)))
+function sector(gradingtype::GradingType, sec)
+  sectorname = Symbol(get(gradingtype, :name, name(gradingtype)))
   return SectorProduct(NamedTuple{(sectorname,)}((sec,)))
 end
 
-function Base.AbstractUnitRange(s::SymmetryType"Nf", t::SiteType"Fermion")
+function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Fermion")
   return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
 end
-# TODO: Write in terms of `SymmetryType"Nf"` definition.
-function Base.AbstractUnitRange(s::SymmetryType"NfParity", t::SiteType"Fermion")
+# 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])
 end
-function Base.AbstractUnitRange(s::SymmetryType"Sz", t::SiteType"Fermion")
+function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Fermion")
   return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
 end
-function Base.AbstractUnitRange(s::SymmetryType"Sz↑", t::SiteType"Fermion")
+function Base.AbstractUnitRange(s::GradingType"Sz↑", t::SiteType"Fermion")
   return gradedrange([sector(s, U1(0)) => 1, sector(s, U1(1)) => 1])
 end
-function Base.AbstractUnitRange(s::SymmetryType"Sz↓", t::SiteType"Fermion")
+function Base.AbstractUnitRange(s::GradingType"Sz↓", t::SiteType"Fermion")
   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::SymmetryType"Nf", t::SiteType"Electron")
+function Base.AbstractUnitRange(s::GradingType"Nf", t::SiteType"Electron")
   return gradedrange([
     sector(s, U1(0)) => 1,
     sector(s, U1(1)) => 1,
@@ -65,8 +65,8 @@ function Base.AbstractUnitRange(s::SymmetryType"Nf", t::SiteType"Electron")
     sector(s, U1(2)) => 1,
   ])
 end
-# TODO: Write in terms of `SymmetryType"Nf"` definition.
-function Base.AbstractUnitRange(s::SymmetryType"NfParity", t::SiteType"Electron")
+# 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,
@@ -74,7 +74,7 @@ function Base.AbstractUnitRange(s::SymmetryType"NfParity", t::SiteType"Electron"
     sector(s, Z{2}(0)) => 1,
   ])
 end
-function Base.AbstractUnitRange(s::SymmetryType"Sz", t::SiteType"Electron")
+function Base.AbstractUnitRange(s::GradingType"Sz", t::SiteType"Electron")
   return gradedrange([
     sector(s, U1(0)) => 1,
     sector(s, U1(1)) => 1,
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 39042c9..c58cb46 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -29,33 +29,33 @@ end
 combine_axes(a) = a
 combine_axes(a, b, rest...) = combine_axes(combine_axes(a, b), rest...)
 
-struct SymmetryType{T,Params}
+struct GradingType{T,Params}
   params::Params
-  function SymmetryType{N}(params::NamedTuple) where {N}
+  function GradingType{N}(params::NamedTuple) where {N}
     return new{N,typeof(params)}(params)
   end
 end
-name(::SymmetryType{T}) where {T} = T
-params(t::SymmetryType) = getfield(t, :params)
-Base.getproperty(t::SymmetryType, name::Symbol) = getfield(params(t), name)
-Base.get(t::SymmetryType, name::Symbol, default) = get(params(t), name, default)
-Base.haskey(t::SymmetryType, name::Symbol) = haskey(params(t), name)
-SymmetryType{N}(; kwargs...) where {N} = SymmetryType{N}((; kwargs...))
-SymmetryType(s::AbstractString; kwargs...) = SymmetryType{Symbol(s)}(; kwargs...)
-function SymmetryType(s::Pair{<:AbstractString,<:AbstractString}; kwargs...)
-  return SymmetryType(first(s); kwargs..., name=last(s))
+name(::GradingType{T}) where {T} = T
+params(t::GradingType) = getfield(t, :params)
+Base.getproperty(t::GradingType, name::Symbol) = getfield(params(t), name)
+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))
 end
-function SymmetryType(s::Pair{<:AbstractString,<:NamedTuple}; kwargs...)
-  return SymmetryType(first(s); kwargs..., last(s)...)
+function GradingType(s::Pair{<:AbstractString,<:NamedTuple}; kwargs...)
+  return GradingType(first(s); kwargs..., last(s)...)
 end
-macro SymmetryType_str(s)
-  return :(SymmetryType{$(Expr(:quote, Symbol(s)))})
+macro GradingType_str(s)
+  return :(GradingType{$(Expr(:quote, Symbol(s)))})
 end
 
-function Base.AbstractUnitRange(symmetry::SymmetryType, t::SiteType)
+function Base.AbstractUnitRange(grading::GradingType, t::SiteType)
   return error("Not implemented.")
 end
-function Base.AbstractUnitRange(symmetry::SymmetryType"Trivial", t::SiteType)
+function Base.AbstractUnitRange(grading::GradingType"Trivial", t::SiteType)
   return Base.OneTo(length(t))
 end
 
@@ -70,8 +70,8 @@ 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, :symmetries)
-    rs = map(symmetry -> AbstractUnitRange(SymmetryType(symmetry), t), t.symmetries)
+  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))

From 2fa4a5e2f361180053d0c561ea072948e4b59390 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 09:00:12 -0500
Subject: [PATCH 06/12] Simplifications

---
 ...uantumOperatorDefinitionsITensorBaseExt.jl | 19 ++++---------------
 src/sitetype.jl                               |  7 +++++++
 2 files changed, 11 insertions(+), 15 deletions(-)

diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index 067b87e..bea962e 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -8,29 +8,18 @@ using QuantumOperatorDefinitions:
 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 (and it may not even
+  # 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 ITensorBase.Index(t::SiteType; kwargs...)
-  return Index(AbstractUnitRange(t); kwargs...)
-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))
 end
 
-function Base.AbstractArray(n::OpName, r::Index)
-  # TODO: Define this with mapped dimnames.
-  # Generalize beyond prime levels with codomain and domain indices.
-  return ITensor(AbstractArray(n, SiteType(r)), (prime(r), dag(r)))
-end
-
-function Base.AbstractArray(n::StateName, r::Index)
-  return ITensor(AbstractArray(n, SiteType(r)), (r,))
-end
-
 end
diff --git a/src/sitetype.jl b/src/sitetype.jl
index c58cb46..28a72cc 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -66,6 +66,9 @@ function Base.length(t::SiteType)
   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.
@@ -76,6 +79,10 @@ function Base.AbstractUnitRange(t::SiteType)
   end
   return Base.OneTo(length(t))
 end
+# kwargs are passed for fancier constructors, like `ITensors.Index`.
+function (rangetype::Type{<:AbstractUnitRange})(t::SiteType; kwargs...)
+  return rangetype(AbstractUnitRange(t); kwargs...)
+end
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]
 Base.axes(t::SiteType) = (AbstractUnitRange(t),)

From c6c5ba0df6642c28601b8873676b75105599355d Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 09:10:42 -0500
Subject: [PATCH 07/12] Add back missing prime

---
 .../QuantumOperatorDefinitionsITensorBaseExt.jl               | 4 ++++
 src/op.jl                                                     | 4 ++--
 src/state.jl                                                  | 2 +-
 3 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index bea962e..ec49ac7 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -22,4 +22,8 @@ function QuantumOperatorDefinitions.has_fermion_string(n::String, r::Index)
   return has_fermion_string(OpName(n), SiteType(r))
 end
 
+function Base.axes(::OpName, domain::Tuple{Vararg{Index}})
+  return (prime.(domain)..., dag.(domain)...)
+end
+
 end
diff --git a/src/op.jl b/src/op.jl
index 7d6d571..c772b38 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -98,7 +98,7 @@ function array(a::AbstractArray, ax::Tuple{Vararg{AbstractUnitRange}})
   return a[ax...]
 end
 
-function state_or_op_axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}})
+function Base.axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}})
   return (domain..., domain...)
 end
 
@@ -108,7 +108,7 @@ function state_or_op_convert(
   domain::Tuple{Vararg{AbstractUnitRange}},
   a::AbstractArray,
 )
-  ax = state_or_op_axes(n, domain)
+  ax = axes(n, domain)
   a′ = reshape(a, length.(ax))
   a′′ = array(a′, ax)
   return convert(arrtype, a′′)
diff --git a/src/state.jl b/src/state.jl
index 55c5155..25b87f2 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -32,7 +32,7 @@ macro state_alias(name1, name2, params...)
   return state_alias_expr(name1, name2)
 end
 
-function state_or_op_axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
+function Base.axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
   return domain
 end
 

From 3a8da93c03715c8ef5fc10c168f2278eaa315542 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 12:37:33 -0500
Subject: [PATCH 08/12] Change memory ordering of operators and states when
 converting to arrays

---
 ...uantumOperatorDefinitionsITensorBaseExt.jl |  5 +++-
 src/op.jl                                     | 27 ++++++++++++++++---
 src/state.jl                                  | 13 ++++++---
 3 files changed, 37 insertions(+), 8 deletions(-)

diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index ec49ac7..86ad084 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -3,7 +3,7 @@ module QuantumOperatorDefinitionsITensorBaseExt
 using ITensorBase: ITensorBase, ITensor, Index, dag, gettag, prime
 using NamedDimsArrays: dename
 using QuantumOperatorDefinitions:
-  QuantumOperatorDefinitions, OpName, SiteType, StateName, has_fermion_string
+  QuantumOperatorDefinitions, @OpName_str, OpName, SiteType, StateName, has_fermion_string
 
 function QuantumOperatorDefinitions.SiteType(r::Index)
   # We pass the axis of the (unnamed) Index because
@@ -25,5 +25,8 @@ end
 function Base.axes(::OpName, domain::Tuple{Vararg{Index}})
   return (prime.(domain)..., dag.(domain)...)
 end
+## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{Index}})
+##   return (prime.(reverse(domain))..., dag.(domain)...)
+## end
 
 end
diff --git a/src/op.jl b/src/op.jl
index c772b38..4c4b514 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -101,6 +101,27 @@ end
 function Base.axes(::OpName, domain::Tuple{Vararg{AbstractUnitRange}})
   return (domain..., domain...)
 end
+function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{Integer}})
+  return axes(n, Base.OneTo.(domain))
+end
+function Base.axes(n::StateOrOpName, domain::Tuple{Vararg{SiteType}})
+  return axes(n, AbstractUnitRange.(domain))
+end
+
+## function Base.axes(::OpName"SWAP", domain::Tuple{Vararg{AbstractUnitRange}})
+##   return (reverse(domain)..., domain...)
+## end
+
+function reversed_sites(n::StateOrOpName, 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)
+end
 
 function state_or_op_convert(
   n::StateOrOpName,
@@ -116,13 +137,13 @@ end
 
 function (arrtype::Type{<:AbstractArray})(n::StateOrOpName, domain::Tuple{Vararg{SiteType}})
   domain′ = AbstractUnitRange.(domain)
-  return state_or_op_convert(n, arrtype, domain′, n(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, n(Int.(length.(domain))...))
+  return state_or_op_convert(n, arrtype, domain, reversed_sites(n, Int.(length.(domain))))
 end
 
 function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...)
@@ -485,7 +506,7 @@ function (n::OpName"Controlled")(domain...)
   # Number of control sites.
   nc = get(params(n), :ncontrol, length(domain) - nt)
   @assert length(domain) == nc + nt
-  d_control = prod(to_dim.(domain[1:nc]))
+  d_control = prod(to_dim.(domain)) - prod(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"()
diff --git a/src/state.jl b/src/state.jl
index 25b87f2..56262c3 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -32,10 +32,6 @@ macro state_alias(name1, name2, params...)
   return state_alias_expr(name1, name2)
 end
 
-function Base.axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
-  return domain
-end
-
 # This compiles operator expressions, such as:
 # ```julia
 # stateexpr("0 + 1") == StateName("0") + StateName("1")
@@ -120,6 +116,15 @@ function state_or_op_expr(ntype::Type, ex::Expr, depth::Int)
   return error("Can't parse expression $ex.")
 end
 
+function Base.axes(::StateName, domain::Tuple{Vararg{AbstractUnitRange}})
+  return domain
+end
+
+function reverse_sites(n::StateName, a::AbstractArray)
+  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...)
 end

From 4540d8c7457b19de6b29e017e225f49c3786375c Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 17:45:28 -0500
Subject: [PATCH 09/12] Add more tests for range construction and ITensor

---
 ...uantumOperatorDefinitionsITensorBaseExt.jl | 30 ++++++++-
 src/op.jl                                     |  2 +-
 src/sitetype.jl                               |  4 +-
 test/test_basics.jl                           | 20 +++++-
 test/test_itensorbaseext.jl                   | 65 +++++++++++++++++--
 5 files changed, 106 insertions(+), 15 deletions(-)

diff --git a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
index 86ad084..84b8b4b 100644
--- a/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
+++ b/ext/QuantumOperatorDefinitionsITensorBaseExt/QuantumOperatorDefinitionsITensorBaseExt.jl
@@ -1,9 +1,15 @@
 module QuantumOperatorDefinitionsITensorBaseExt
 
-using ITensorBase: ITensorBase, ITensor, Index, dag, gettag, prime
+using ITensorBase: ITensorBase, ITensor, Index, dag, gettag, prime, settag
 using NamedDimsArrays: dename
 using QuantumOperatorDefinitions:
-  QuantumOperatorDefinitions, @OpName_str, OpName, SiteType, StateName, has_fermion_string
+  QuantumOperatorDefinitions,
+  @OpName_str,
+  OpName,
+  SiteType,
+  StateName,
+  has_fermion_string,
+  name
 
 function QuantumOperatorDefinitions.SiteType(r::Index)
   # We pass the axis of the (unnamed) Index because
@@ -16,6 +22,10 @@ function QuantumOperatorDefinitions.SiteType(r::Index)
   )
 end
 
+function (rangetype::Type{<:Index})(t::SiteType)
+  return settag(rangetype(AbstractUnitRange(t)), "sitetype", String(name(t)))
+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)
@@ -29,4 +39,20 @@ end
 ##   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)
+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))
+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))
+end
+
 end
diff --git a/src/op.jl b/src/op.jl
index 4c4b514..7aee6a0 100644
--- a/src/op.jl
+++ b/src/op.jl
@@ -506,7 +506,7 @@ function (n::OpName"Controlled")(domain...)
   # Number of control sites.
   nc = get(params(n), :ncontrol, length(domain) - nt)
   @assert length(domain) == nc + nt
-  d_control = prod(to_dim.(domain)) - prod(domain[(nc + 1):end])
+  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"()
diff --git a/src/sitetype.jl b/src/sitetype.jl
index 28a72cc..67936d9 100644
--- a/src/sitetype.jl
+++ b/src/sitetype.jl
@@ -80,8 +80,8 @@ function Base.AbstractUnitRange(t::SiteType)
   return Base.OneTo(length(t))
 end
 # kwargs are passed for fancier constructors, like `ITensors.Index`.
-function (rangetype::Type{<:AbstractUnitRange})(t::SiteType; kwargs...)
-  return rangetype(AbstractUnitRange(t); kwargs...)
+function (rangetype::Type{<:AbstractUnitRange})(t::SiteType)
+  return rangetype(AbstractUnitRange(t))
 end
 Base.size(t::SiteType) = (length(t),)
 Base.size(t::SiteType, dim::Integer) = size(t)[dim]
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 31d0c96..dfafddc 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -162,8 +162,10 @@ const elts = (real_elts..., complex_elts...)
     @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))") ==
-      reshape(exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))), (2, 2, 2, 2))
+    @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+"))
@@ -201,7 +203,19 @@ const elts = (real_elts..., complex_elts...)
     @test state("2", 3) == [0, 0, 1]
 
     @test state("|0⟩ + 2|+⟩") == state("0") + 2 * state("+")
-    @test state("|0⟩ ⊗ |+⟩") == reshape(kron(state("0"), state("+")), (2, 2))
+    @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}
+
+    @test op("X", 2) == op("X", Base.OneTo(2))
+    @test op("X", 3) == op("X", Base.OneTo(3))
   end
   @testset "Electron/tJ" begin
     for (ns, x) in (
diff --git a/test/test_itensorbaseext.jl b/test/test_itensorbaseext.jl
index e71fe46..1df61b0 100644
--- a/test/test_itensorbaseext.jl
+++ b/test/test_itensorbaseext.jl
@@ -1,13 +1,64 @@
-using ITensorBase: ITensor, Index, prime
-using QuantumOperatorDefinitions: op, state
-using Test: @test, @testset
+using ITensorBase: ITensor, Index, gettag, prime, settag
+using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state
+using Test: @test, @test_broken, @testset
+
+@testset "ITensorBaseExt" begin
+  i = Index(SiteType("S=1/2"))
+  @test gettag(i, "sitetype") == "S=1/2"
 
-@testset "QuantumOperatorDefinitionsITensorBaseExt" begin
   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 = Index(2)
   a = op("X", i)
-  @test a == ITensor([0 1; 1 0], (prime(i), 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(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)
 
-  a = state(1, i)
-  @test a == ITensor([1, 0], (i,))
+  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)
 end

From 4fc466de9eebe7cf8e8933ee9fbc86f9a6ae8d14 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 19:01:35 -0500
Subject: [PATCH 10/12] Test symmetries

---
 ...umOperatorDefinitionsSymmetrySectorsExt.jl |  16 ++-
 test/Project.toml                             |   4 +
 test/test_itensorbaseext.jl                   |   2 +-
 test/test_symmetrysectorsext.jl               | 113 ++++++++++++++++++
 4 files changed, 131 insertions(+), 4 deletions(-)
 create mode 100644 test/test_symmetrysectorsext.jl

diff --git a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
index 9814efc..e21f6c1 100644
--- a/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
+++ b/ext/QuantumOperatorDefinitionsSymmetrySectorsExt/QuantumOperatorDefinitionsSymmetrySectorsExt.jl
@@ -1,11 +1,21 @@
 module QuantumOperatorDefinitionsSymmetrySectorsExt
 
 using BlockArrays: blocklasts, blocklengths
-using GradedUnitRanges: GradedOneTo, gradedrange
+using GradedUnitRanges: AbstractGradedUnitRange, GradedOneTo, gradedrange
 using LabelledNumbers: label, labelled, unlabel
 using QuantumOperatorDefinitions:
-  QuantumOperatorDefinitions, @SiteType_str, @GradingType_str, SiteType, GradingType, name
-using SymmetrySectors: ×, SectorProduct, U1, Z
+  QuantumOperatorDefinitions,
+  @SiteType_str,
+  @GradingType_str,
+  SiteType,
+  GradingType,
+  OpName,
+  name
+using SymmetrySectors: ×, dual, SectorProduct, U1, Z
+
+function Base.axes(::OpName, domain::Tuple{Vararg{AbstractGradedUnitRange}})
+  return (domain..., dual.(domain)...)
+end
 
 sortedunion(a, b) = sort(union(a, b))
 function QuantumOperatorDefinitions.combine_axes(a1::GradedOneTo, a2::GradedOneTo)
diff --git a/test/Project.toml b/test/Project.toml
index 32b3272..43654f1 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,11 +1,15 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
diff --git a/test/test_itensorbaseext.jl b/test/test_itensorbaseext.jl
index 1df61b0..cb274b5 100644
--- a/test/test_itensorbaseext.jl
+++ b/test/test_itensorbaseext.jl
@@ -1,6 +1,6 @@
 using ITensorBase: ITensor, Index, gettag, prime, settag
 using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state
-using Test: @test, @test_broken, @testset
+using Test: @test, @testset
 
 @testset "ITensorBaseExt" begin
   i = Index(SiteType("S=1/2"))
diff --git a/test/test_symmetrysectorsext.jl b/test/test_symmetrysectorsext.jl
new file mode 100644
index 0000000..af65135
--- /dev/null
+++ b/test/test_symmetrysectorsext.jl
@@ -0,0 +1,113 @@
+using BlockArrays: AbstractBlockArray, blocklengths
+using BlockSparseArrays: BlockSparseArray
+using GradedUnitRanges: blocklabels
+using ITensorBase: ITensor, Index, gettag, prime, settag
+using QuantumOperatorDefinitions: OpName, SiteType, StateName, op, state
+using SymmetrySectors: SectorProduct, U1, Z
+using NamedDimsArrays: dename
+using Test: @test, @test_broken, @testset
+
+@testset "SymmetrySectorsExt" begin
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  r = AbstractUnitRange(t)
+  @test r == 1:2
+  @test blocklabels(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 blocklabels(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 blocklabels(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 blocklabels(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 blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  @test blocklabels(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+
+  # TODO: There is a bug slicing `BitVector` by `GradedOneTo`, investigate.
+  # See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
+  t = SiteType("S=1/2"; gradings=("Sz",))
+  @test_broken state("0", t)
+
+  # 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 AbstractBlockArray
+  # TODO: Currently slicing a dense array by graded ranges outputs a `BlockedArray`
+  # rather than a `BlockSparseArray`.
+  # See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
+  @test_broken a isa BlockSparseArray
+  (r1,) = axes(a)
+  @test blocklabels(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 AbstractBlockArray
+  @test_broken a isa BlockSparseArray
+  (r1, r2) = axes(a)
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  # TODO: This is a bug in indexing with GradedUnitRangeDual, fix this.
+  @test_broken blocklabels(r2) ==
+    [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+end
+
+@testset "SymmetrySectorsExt + 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 blocklabels(dename(i)) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(dename(i)) == [1, 1]
+
+  i′ = prime(i)
+  a = op("σ⁺", i)
+  # TODO: The indices should be `(i′, dual(i))`.
+  @test a == ITensor([0 2; 0 0], (i′, i))
+  a′ = dename(a)
+  @test a′ isa AbstractBlockArray
+  @test_broken a′ isa BlockSparseArray
+  # TODO: Test these without denaming `a`.
+  (r1, r2) = axes(a′)
+  @test blocklabels(r1) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(1)))]
+  @test blocklengths(r1) == [1, 1]
+  # TODO: This is a bug in indexing with GradedUnitRangeDual, fix this.
+  @test_broken blocklabels(r2) ==
+    [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
+  @test blocklengths(r2) == [1, 1]
+end

From 3a5a6d7863d427545abdc923b6ab3215ea3a6b56 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 19:09:50 -0500
Subject: [PATCH 11/12] Add missing test dependency

---
 test/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/Project.toml b/test/Project.toml
index 43654f1..69e5ba0 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"

From b36566fe908bec02b210b4b3037250ab00a91172 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 31 Jan 2025 20:50:16 -0500
Subject: [PATCH 12/12] Fix tests in Julia 1.10

---
 test/test_symmetrysectorsext.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/test_symmetrysectorsext.jl b/test/test_symmetrysectorsext.jl
index af65135..e7da36d 100644
--- a/test/test_symmetrysectorsext.jl
+++ b/test/test_symmetrysectorsext.jl
@@ -54,10 +54,10 @@ using Test: @test, @test_broken, @testset
   @test blocklabels(r2) == [SectorProduct((; Sz=U1(0))), SectorProduct((; Sz=U1(-1)))]
   @test blocklengths(r2) == [1, 1]
 
-  # TODO: There is a bug slicing `BitVector` by `GradedOneTo`, investigate.
-  # See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
+  # TODO: There is a bug slicing `BitVector` by `GradedOneTo` in Julia 1.11,
+  # investigate. See: https://github.com/ITensor/GradedUnitRanges.jl/issues/9
   t = SiteType("S=1/2"; gradings=("Sz",))
-  @test_broken state("0", t)
+  @test state("0", t) == [1, 0] broken = VERSION ≥ v"1.11"
 
   # Force conversion to `Vector{Float64}` before conversion,
   # since there is a bug slicing `BitVector` by `GradedOneTo`.