Arrays that may have half-integer indices, commonly encountered while working with rotations and spin.
This package is very much a WIP, and bugs are expected. Please open an issue if you encounter a bug.
The package is to be used alongside HalfIntegers.jl. This is installed automatically with this package, and may be imported as shown below.
julia> ]
pkg> add https://github.com/jishnub/HalfIntegerArrays.jl
julia> using HalfIntegerArrays
julia> using HalfIntegerArrays.HalfIntegersThere are two types exported: HalfIntArray and SpinMatrix. The former represents an arbitrary array with possibly half-integral axes, whereas the latter represents a square matrix corresponding to an angular momentum j. In the second case the array is guaranteed to have a size of (2j+1, 2j+1), and axes (-j:j, -j:j).
HalfIntArrays are wrappers around AbstractArrays, and may be constructed in two ways: firstly by providing the axes for the parent array, eg:
julia> h = HalfIntArray(ones(Int,3,3), -1:1, -1:1)
3×3 HalfIntArray(::Array{Int64,2}, -1:1, -1:1) with eltype Int64 with indices -1:1×-1:1:
1 1 1
1 1 1
1 1 1
# We may wrap structured arrays
julia> import LinearAlgebra: Diagonal
julia> h = HalfIntArray(Diagonal([1,2]), -1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(::Diagonal{Int64,Array{Int64,1}}, -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
1 ⋅
⋅ 2and secondly by using the typed constructor
julia> h = HalfIntArray{Float64}(undef,-1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(::Array{Float64,2}, -1/2:1/2, -1/2:1/2) with eltype Float64 with indices -1/2:1/2×-1/2:1/2:
0.0 0.0
6.89924e-310 0.0
julia> h = HalfIntArray{Union{Float64,Missing}}(missing,-1:1, -1:1)
3×3 HalfIntArray(::Array{Union{Missing, Float64},2}, -1:1, -1:1) with eltype Union{Missing, Float64} with indices -1:1×-1:1:
missing missing missing
missing missing missing
missing missing missingIn the second case an underlying array of an appropriate size is initialized that has the specified element type.
A SpinMatrix may be constructed as
julia> s = SpinMatrix(zeros(3,3)) # the angular momentum is inferred
3×3 SpinMatrix(::Array{Float64,2}, 1) with eltype Float64 with indices -1:1×-1:1:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
# Specifying the angular momentum will allocate an appropriate array
julia> s = SpinMatrix{ComplexF64}(undef, 1//2)
2×2 SpinMatrix(::Array{Complex{Float64},2}, 1/2) with eltype Complex{Float64} with indices -1/2:1/2×-1/2:1/2:
6.91635e-310+6.91635e-310im 6.91637e-310+6.91637e-310im
6.91635e-310+6.91637e-310im 6.91637e-310+6.91637e-310imA SpinMatrix requires the parent matrix to have 1-based indexing and have a size of (2j+1, 2j+1).
The arrays may be indexed with integral or half-integral values. For optimal performance it's preferable to use integral, floating-point or HalfInt types as indices. A HalfInt may be constructed using the function half, such that half(n) == n/2. Alternately they may also be constructed as HalfInt(n) which is numerically equivalent to n. Rational numbers may be used as well, but these might not be as performant.
An example with a half-integral spin:
julia> s = SpinMatrix(reshape(1:4,2,2), 1//2)
2×2 SpinMatrix(reshape(::UnitRange{Int64}, 2, 2), 1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
1 3
2 4
julia> s[1//2,1//2]
4
julia> import HalfIntegerArrays: half
julia> s[-half(1),half(1)]
3An example with integral spin:
julia> s = SpinMatrix(reshape(1:9,3,3))
3×3 SpinMatrix(reshape(::UnitRange{Int64}, 3, 3), 1) with eltype Int64 with indices -1:1×-1:1:
1 4 7
2 5 8
3 6 9
julia> s[1,1]
9Julia's default CartesianIndex, CartesianIndices and LinearIndices types require integer indices, therefore these are not compatible with HalfIntArrays. This package exports the equivalent types CartesianIndexHalfInt, CartesianIndicesHalfInt and LinearIndicesHalfInt that support both integer and half-integer indices. These are therefore the safe choices for indexing into HalfIntArrays.
julia> h = HalfIntArray(reshape(1:4,2,2), 0:1, 0:1)
2×2 HalfIntArray(reshape(::UnitRange{Int64}, 2, 2), 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
1 3
2 4
julia> cinds = eachindex(IndexCartesian(), h)
2×2 CartesianIndicesHalfInt{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} with indices 0:1×0:1:
CartesianIndexHalfInt(0, 0) CartesianIndexHalfInt(0, 1)
CartesianIndexHalfInt(1, 0) CartesianIndexHalfInt(1, 1)
julia> h[cinds[1,0]]
2Indexing with CartesianIndices works as well for arrays with integer indices.
Broadcasting works, but is a bit slow at the moment. For optimal performance it's better to broadcast on the parent array.
julia> h = HalfIntArray(reshape(1:4,2,2), -1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(reshape(::UnitRange{Int64}, 2, 2), -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
1 3
2 4
julia> h .+ h
2×2 HalfIntArray(::Array{Int64,2}, -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
2 6
4 8
julia> @btime $h .+ $h;
193.130 ns (2 allocations: 192 bytes)
julia> @btime parent($h) .+ parent($h);
94.443 ns (1 allocation: 160 bytes)Upon broadcasting with SpinMatrix types, the result will be a HalfIntArray. This behaviour might change in the future.
julia> s = SpinMatrix(reshape(1:4,2,2), half(1))
2×2 SpinMatrix(reshape(::UnitRange{Int64}, 2, 2), 1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
1 3
2 4
julia> s .+ s
2×2 HalfIntArray(::Array{Int64,2}, -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
2 6
4 8
# It's possible to recreate the wrapper
julia> s .+ s |> SpinMatrix
2×2 SpinMatrix(::Array{Int64,2}, 1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
2 6
4 8Comparison with OffsetArrays
A HalfIntArray with integer axes is equivalent to an OffsetArray, except that a HalfIntArray is somewhat less performant when it comes to indexing.
julia> using OffsetArrays
julia> oa = OffsetArray(reshape(1:9,3,3), -1:1, -1:1)
3×3 OffsetArray(reshape(::UnitRange{Int64}, 3, 3), -1:1, -1:1) with eltype Int64 with indices -1:1×-1:1:
1 4 7
2 5 8
3 6 9
julia> h = HalfIntArray(parent(oa), axes(oa)...)
3×3 HalfIntArray(reshape(::UnitRange{Int64}, 3, 3), -1:1, -1:1) with eltype Int64 with indices -1:1×-1:1:
1 4 7
2 5 8
3 6 9
julia> h == oa
true