# Can't iterate over an AbstractArray with custom indexing?

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 
Stacktrace:
 setindex! at ./array.jl:826 [inlined]
 setindex! at /Users/lamont/Desktop/WrightFisherSim/SphericalTensor.jl:61 [inlined]
 collect_to! at ./array.jl:715 [inlined]
 collect_to_with_first!(::Spherical{2}, ::Complex{Float64}, ::Base.Generator{Spherical{2},typeof(identity)}, ::Tuple{MRange,Int64}) at ./array.jl:689
 collect(::Base.Generator{Spherical{2},typeof(identity)}) at ./array.jl:670
 top-level scope at REPL: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

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]
``````

Interesting clue