Can't iterate over an AbstractArray with custom indexing?

General Usage
1

So (like many other physicists trying out Julia), I went and tried to build a spherical tensor type where the size L is a type parameter and indexing runs from -L: L. I followed the manual and defined an abstract unit range type and corresponding methods for Base.axes, similar, getindex, and setindex!.

struct MRange <: AbstractUnitRange{Int}
    start::Int
    stop::Int
end

function MRange(L)
    MRange(-L,L)
end

# definition
struct Spherical{L} <: AbstractArray{Complex{Float64},1}
    x::Array{Complex{Float64},1}
end

# standard constructor
function Spherical(x)
    Spherical{Int((length(x)-1)/2)}(x)
end

# initialization: unit vector in m = 0
function Spherical{L}() where L
    Spherical(Complex{Float64}.( (-L:L) .== 0))
end

Base.axes(y::Spherical{L}) where L = (MRange(L),)
Base.size(y::Spherical{L}) where L = (2*L+1,)
function Base.getindex(y::Spherical{L}, m::Int) where L
    y.x[L+m+1]
end
function Base.setindex!(y::Spherical{L}, v ,m::Int) where L
    y.x[L+m+1] = v
end
function Base.similar(A::AbstractArray, T::Type, shape::Tuple{MRange,Vararg{MRange}})
    Spherical(A)
end
function Base.similar(f::Any, shape::Tuple{MRange,Vararg{MRange}}) 
    #something didn't work here from the recipe book, might be a bug
    Spherical(first(shape))
end

The dispatcher could not find the f::Union{Function,DataType} signature recommended in the manual, but I loosened it a little and suddenly most of Base was talking to my Spherical’s like they were there all along! Some beautiful code like this clebsch-gordon product worked instantly <3

# clebsch product into L3 irrep 
function cleb(y1::Spherical{L1},y2::Spherical{L2},L3::Int) where {L1, L2}
    if abs(L1 - L2) <= L3 <= L1 + L2
        return Spherical{L3}(
            [ sum( clebschgordan(L1, m1, L2, m3 - m1, L3, m3) * y1[m1] * y2[m3 - m1] 
                for m1 in -L1:L1 if abs(m3 - m1) <= L2)
            for m3 in -L3:L3 ]
         )
    else
        nothing
    end
end

Getting the indices right the normal way would have been a nightmare.

But, of course, not everything works. Arrays of Spherical don’t print well (possibly this?)

julia> [Spherical{2}(),Spherical{2}()]
2-element Array{Spherical{2},1}:
 [0.0 + 0.0im, 0.0 + 0.0im, #undef, #undef, #undef]
 [0.0 + 0.0im, 0.0 + 0.0im, #undef, #undef, #undef]

but most disturbingly and surprisingly I can’t iterate over Spherical’s safely as AbstractArrays!


 julia> [ii for ii in Spherical{2}()]
ERROR: BoundsError: attempt to access 5-element Array{Complex{Float64},1} at index [6]
Stacktrace:
 [1] setindex! at ./array.jl:826 [inlined]
 [2] setindex! at /Users/lamont/Desktop/WrightFisherSim/SphericalTensor.jl:61 [inlined]
 [3] collect_to! at ./array.jl:715 [inlined]
 [4] collect_to_with_first!(::Spherical{2}, ::Complex{Float64}, ::Base.Generator{Spherical{2},typeof(identity)}, ::Tuple{MRange,Int64}) at ./array.jl:689
 [5] collect(::Base.Generator{Spherical{2},typeof(identity)}) at ./array.jl:670
 [6] top-level scope at REPL[47]:1

Do I need more methods? How do I figure out how to fix this? How broken are things going to be down the line? Also, would there be any advantage to using AbstractVector instead?

1 Like
2

Not sure I can easily explain why, but notice this important difference from OffsetArrays:

julia> Spherical(ones(3))
3-element Spherical{1} with indices -1:1:1:
 1.0 + 0.0im
 1.0 + 0.0im
 1.0 + 0.0im

julia> axes(ans,1)
-1:1:1

julia> axes(ans,1)
Base.OneTo(3)

julia> using OffsetArrays

julia> OffsetArray(ones(3).+0im, -2)
3-element OffsetArray(::Array{Complex{Float64},1}, -1:1) with eltype Complex{Float64} with indices -1:1:
 1.0 + 0.0im
 1.0 + 0.0im
 1.0 + 0.0im

julia> axes(ans,1)
-1:1

julia> axes(ans,1) # this is different
-1:1 

julia> [Spherical(ones(3))]
1-element Array{Spherical{1},1}:
 [1.0 + 0.0im, #undef, #undef]

julia> [OffsetArray(ones(3), -2)]
1-element Array{OffsetArray{Float64,1,Array{Float64,1}},1}:
 [1.0, 1.0, 1.0]
3

Interesting clue