HalfIntegerArrays

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.

Prerequisites

The package is to be used alongside HalfIntegers.jl. This is installed automatically with this package, and may be imported as shown below.

Installation

julia> ]
pkg> add https://github.com/jishnub/HalfIntegerArrays.jl

julia> using HalfIntegerArrays
julia> using HalfIntegerArrays.HalfIntegers

Usage

There 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  ⋅
 ⋅  2

and 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  missing

In 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-310im

A SpinMatrix requires the parent matrix to have 1-based indexing and have a size of (2j+1, 2j+1).

Indexing

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)]
3

An 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]
9

Linear and Cartesian Indexing

Julia'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]]
2

Indexing with CartesianIndices works as well for arrays with integer indices.

Broadcasting

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  8

Comparison 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