SphericalHarmonicArrays

Arrays to store spherical-harmonic coefficients, that may be indexed by modes as well as array indices. The type used for this is named SHArray, and the aliases SHVector and SHMatrix are exported for covenience. SHArray is a wrapper around an underlying parent array, usually dense, that is indexed according to iterators that are specified while constructing the type. The arrays may have mixed axes, where certain axes are indexed using spherical harmonic modes whereas the others are indexed like the parent array.

Getting Started

Installing

julia> ]
pkg> add SphericalHarmonicArrays

julia> using SphericalHarmonicArrays

Usage

This package uses iterators from SphericalHarmonicModes.jl for indexing. All the iterators available in the package may be used as axes. Take a look at that package to understand more about the iterators being used here.

SHVector

An SHVector is a 1D array indexed using spherical harmonic modes.

Creating an SHVector

The simplest constructor assigns the array automatically based on the number of modes specified.

julia> modes=LM(1:1)
Spherical harmonic modes with l increasing faster than m
(l_min = 1, l_max = 1, m_min = -1, m_max = 1)

julia> SHVector{Float64}(modes)
3-element SHArray(::Array{Float64,1}, (LM(1:1, -1:1),)):
 0.0
 0.0
 0.0

The parent array may be preallocated.

julia> v = ones(3);

julia> shv = SHVector(v,modes)
3-element SHArray(::Array{Float64,1}, (LM(1:1, -1:1),)):
 1.0
 1.0
 1.0

julia> v[1] = 4 # update the parent array
4

julia> shv # updated as well
3-element SHArray{Float64,1,Array{Float64,1},Tuple{LM},1}:
 4.0
 1.0
 1.0

Indexing

An SHVector may be indexed either linearly as a normal vector, or using the modes that are stored in the array. Linear indexing is faster as this simply passes the indices to the parent, so this is what should be used if all the indices of the array are being iterated over. Modes need to be specified as a tuple of integers, eg. (l,m), corresponding to the type of axis iterator that was used to create the SHVector.

julia> v = [1,2,3]; shv=SHVector(v, LM(1:1))
3-element SHArray{Int64,1,Array{Int64,1},Tuple{LM},1}:
 1
 2
 3

julia> shv[2]
2

julia> shv[(1,0)] # indexed using (l,m)
2

julia> @btime $shv[2] # Linear indexing
  2.026 ns (0 allocations: 0 bytes)
2

julia> mode=(1,0); @btime $shv[$mode] # Indexing using modes
  9.774 ns (0 allocations: 0 bytes)
2

julia> @. shv = 56 # broadcasting works as expected
3-element Array{Int64,1}:
 56
 56
 56

julia> shv[(1,-1)] = 6 # can set indices using modes, in this case it's specificed as an (l,m) pair
6

julia> shv
3-element SHArray{Int64,1,Array{Int64,1},Tuple{LM},1}:
  6
 56
 56

Note that in this specific case it is more efficient to define the mode range as LM(SingleValuedRange(1)). Using this, we obtain

julia> shv = SHVector(v, LM(SingleValuedRange(1)));

julia> mode=(1,0); @btime $shv[$mode]
  3.734 ns (0 allocations: 0 bytes)
2

Indexing operations are significantly more performant for arrays constructed using these special ranges.

SHMatrix

An SHMatrix is a 2D array with both axes storing spherical harmonic coefficients.

Creating an SHMatrix

The constructors are similar to those of SHVector.

julia> SHMatrix{Float64}(LM(1:2, 1:1),LM(1:1))
2×3 SHArray(::Array{Float64,2}, (LM(1:2, 1:1), LM(1:1, -1:1))):
 0.0  0.0  0.0
 0.0  0.0  0.0

# may combine different iterators as axes
julia> SHMatrix{Float64}(LM(1:2, 1:1),ML(1:1))
2×3 SHArray(::Array{Float64,2}, (LM(1:2, 1:1), ML(1:1, -1:1))):
 0.0  0.0  0.0
 0.0  0.0  0.0

Indexing

The matrix elements may be accessed with a combination of mode indices or the index style of the parent array.

julia> shm = SHMatrix{Float64}(LM(1:2,1:1), LM(1:1))
2×3 SHArray(::Array{Float64,2}, (LM(1:2, 1:1), LM(1:1, -1:1))):
 0.0  0.0  0.0
 0.0  0.0  0.0

# Linear indexing works if the parent array supports it
julia> for i in eachindex(shm)
       shm[i]=i
       end

julia> shm
2×3 SHArray(::Array{Float64,2}, (LM(1:2, 1:1), LM(1:1, -1:1))):
 1.0  3.0  5.0
 2.0  4.0  6.0

# Both axes may be indexed using modes
julia> shm[(1,1),(1,0)]
3.0

# Any axis may be indexed using the corresponding mode
julia> shm[(1,1), 2]
3.0 + 0.0im

# May use any combination of Cartesian and fancy mode indexing
julia> shm[1,2] == shm[(1,1),2] == shm[1,(1,0)] ==shm[(1,1),(1,0)]
true

julia> mode1=(1,1);mode2=(1,0); @btime $shm[$mode1,$mode2] # twice as expensive as SHVector
  19.026 ns (0 allocations: 0 bytes)
3.0 + 0.0im

SHArray

This is the most general type of an array with any axis possibly being indexed using a collection of modes. This differs from the previously described aliases as some (or all) axes need not be indexed using modes.

Creating an SHArray

julia> sha = SHArray(zeros(1),(LM(1:1,0:0),))
1-element SHArray(::Array{Float64,1}, (LM(1:1, 0:0),)):
 0.0

# SHVector is an alias for a 1D SHArray that is indexed with modes
julia> sha isa SHVector
true

julia> SHArray{Float64}((LM(1:1,0:0),1:2)) # supports OffsetArrays
1×2 SHArray(OffsetArray(::Array{Float64,2}, 1:1, 1:2), (LM(1:1, 0:0), 1:2)) with indices 1:1×1:2:
 0.0  0.0

The arrays may have mixed axes, where some store spherical harmonic modes and some don't.

julia> sha = SHArray{Float64}((LM(1:1,0:0),-1:1,ML(0:1,0:0)))
1×3×2 SHArray(OffsetArray(::Array{Float64,3}, 1:1, -1:1, 1:2), (LM(1:1, 0:0), -1:1, ML(0:1, 0:0))) with indices 1:1×-1:1×1:2:
[:, :, 1] =
 0.0  0.0  0.0

[:, :, 2] =
 0.0  0.0  0.0

Indexing

Indexing is similar to SHVector and SHMatrix.

julia> SHArray{Float64}((1:1, LM(1:1,0:1)))
1×2 SHArray(OffsetArray(::Array{Float64,2}, 1:1, 1:2), (1:1, LM(1:1, 0:1))) with indices 1:1×1:2:
 0.0  0.0

julia> sha[1,(1,0)] = 4 # first index
4

julia> sha[1,2] = 5 # second index
5

julia> sha
1×2 SHArray(OffsetArray(::Array{Float64,2}, 1:1, 1:2), (1:1, LM(1:1, 0:1))) with indices 1:1×1:2:
 4.0  5.0

Broadcasting

SHArrays retain information about their modes upon broadcasting. If multiple SHArrays are involved in a broadcast operation, the result has the same axes as the one with the most dimensions. The dimensions being broadcasted over, if indexed with modes, have to exactly match for all the SHArrays involved in the operation.

julia> s = SHMatrix{Float64}(LM(1:1,0:0),LM(1:1,-1:0)); s .= 4
1×2 SHArray(::Array{Float64,2}, (LM(1:1, 0:0), LM(1:1, -1:0))):
 4.0  4.0

julia> s + s
1×2 SHArray(::Array{Float64,2}, (LM(1:1, 0:0), LM(1:1, -1:0))):
 8.0  8.0

julia> s .* s
1×2 SHArray(::Array{Float64,2}, (LM(1:1, 0:0), LM(1:1, -1:0))):
 16.0  16.0

julia> sv = SHVector{Float64}(first(SphericalHarmonicArrays.modes(s))); sv .= 6;

julia> s .* sv # Leading dimensions of s and sv are the same
1×2 SHArray(::Array{Float64,2}, (LM(1:1, 0:0), LM(1:1, -1:0))):
 24.0  24.0

Broadcasting operations might be slow, so watch out for performance drops.