Trying "Base.mapreduce!" on an AxisArray

[Balinus]

Hello, I'm trying the following code (which simply calculate the number of element that are over a threshold 1). Probably a problem on my side (e.g. not using the correct command for what I'm trying to do).

Code
julia> Base.mapreducedim!(t -> t >= 1, +, view(FD, 1:1, :, :), view(A.data, idx, :, :))

Some Outputs

julia> typeof(A.data)
AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:time,Array{Date,1}},AxisArrays.Axis{:lon,Array{Float64,1}},AxisArrays.Axis{:lat,Array{Float64,1}}}}

julia> typeof(FD)
Array{Int64,3}

julia> Base.mapreducedim!(t -> t >= 1, +, view(FD, 1:1, :, :), view(A.data, 1:365, :, :))
------ MethodError --------------------- Stacktrace (most recent call last)

 [1] — mapreducedim!(::Function, ::Function, ::SubArray{Int64,3,Array{Int64,3},Tuple{UnitRange{Int64},Colon,Colon},false}, ::AxisArrays.AxisArray{Float64,3,SubArray{Float64,3,Array{Float64,3},Tuple{UnitRange{Int64},Colon,Colon},false},Tuple{AxisArrays.Axis{:time,Array{Date,1}},AxisArrays.Axis{:lon,Array{Float64,1}},AxisArrays.Axis{:lat,Array{Float64,1}}}}) at reducedim.jl:222

 [2] — _mapreducedim!(::##39#40, ::Base.#+, ::SubArray{Int64,3,Array{Int64,3},Tuple{UnitRange{Int64},Colon,Colon},false}, ::AxisArrays.AxisArray{Float64,3,SubArray{Float64,3,Array{Float64,3},Tuple{UnitRange{Int64},Colon,Colon},false},Tuple{AxisArrays.Axis{:time,Array{Date,1}},AxisArrays.Axis{:lon,Array{Float64,1}},AxisArrays.Axis{:lat,Array{Float64,1}}}}) at reducedim.jl:206

 [3] — macro expansion at simdloop.jl:73 [inlined]

 [4] — macro expansion at reducedim.jl:207 [inlined]

MethodError: no method matching isless(::AxisArrays.AxisArray{Float64,1,Array{Float64,1},Tuple{AxisArrays.Axis{:dim_4,Base.OneTo{Int64}}}}, ::Int64)
Closest candidates are:
  isless(::Char, ::Integer) at deprecated.jl:49
  isless(::AbstractFloat, ::Real) at operators.jl:42
  isless(::Real, ::Real) at operators.jl:75
  ...

Using a regular Array{Float64,3}, the command works:

julia> R = randn(20440,128,64)
julia> Base.mapreducedim!(t -> t >= 1, +, view(FD, 1:1, :, :), view(R, 1:365, :, :))
1×128×64 SubArray{Int64,3,Array{Int64,3},Tuple{UnitRange{Int64},Colon,Colon},false}:
[:, :, 1] =
 58  59  61  67  63  46  64  54  52  54  50  51  59  55  73  53  …  55  56  51  56  64  62  68  46  61  66  57  55  61  56  48

...

[:, :, 64] =
 66  48  55  49  70  67  56  59  53  70  40  52  60  49  52  59  …  45  64  59  56  59  63  71  72  67  58  48  46  65  67  48

julia> typeof(R)
Array{Float64,3}

[Balinus]

I tried creating the output Array as a proper AxisArray but I got another error. Thanks for any help!

julia> G = AxisArray(FD, Axis{:time}(numYears[1]:numYears[end]))
3-dimensional AxisArray{Int64,3,...} with axes:
    :time, 1950:2005
    :col, Base.OneTo(128)
    :page, Base.OneTo(64)
And data, a 56×128×64 Array{Int64,3}:
...

julia> Base.mapreducedim!(t -> t >= 1, +, view(FD, 1:1, :, :), view(A.data, 1:365, :, :))
------ MethodError --------------------- Stacktrace (most recent call last)

 [1] — mapreducedim!(::Function, ::Function, ::SubArray{Int64,3,Array{Int64,3},Tuple{UnitRange{Int64},Colon,Colon},false}, ::AxisArrays.AxisArray{Float64,3,SubArray{Float64,3,Array{Float64,3},Tuple{UnitRange{Int64},Colon,Colon},false},Tuple{AxisArrays.Axis{:time,Array{Date,1}},AxisArrays.Axis{:lon,Array{Float64,1}},AxisArrays.Axis{:lat,Array{Float64,1}}}}) at reducedim.jl:222

 [2] — _mapreducedim!(::##45#46, ::Base.#+, ::SubArray{Int64,3,Array{Int64,3},Tuple{UnitRange{Int64},Colon,Colon},false}, ::AxisArrays.AxisArray{Float64,3,SubArray{Float64,3,Array{Float64,3},Tuple{UnitRange{Int64},Colon,Colon},false},Tuple{AxisArrays.Axis{:time,Array{Date,1}},AxisArrays.Axis{:lon,Array{Float64,1}},AxisArrays.Axis{:lat,Array{Float64,1}}}}) at reducedim.jl:206

 [3] — macro expansion at simdloop.jl:73 [inlined]

 [4] — macro expansion at reducedim.jl:207 [inlined]

MethodError: no method matching isless(::AxisArrays.AxisArray{Float64,1,Array{Float64,1},Tuple{AxisArrays.Axis{:dim_4,Base.OneTo{Int64}}}}, ::Int64)
Closest candidates are:
  isless(::Char, ::Integer) at deprecated.jl:49
  isless(::AbstractFloat, ::Real) at operators.jl:42
  isless(::Real, ::Real) at operators.jl:75
  ...

[timholy]

Simpler repro:

using AxisArrays

using Base: safe_tail, Broadcast, indices1

function simple_reducedim!(f, op, R::AbstractArray, A::AbstractArray)
    indsAt, indsRt = safe_tail(indices(A)), safe_tail(indices(R)) # handle d=1 manually
    keep, Idefault = Broadcast.shapeindexer(indsAt, indsRt)
    i1 = first(indices1(R))
    for IA in CartesianRange(indsAt)
        IR = Broadcast.newindex(IA, keep, Idefault)
        r0 = R[i1,IR]
        r = r0
        for i in indices(A, 1)
            val = A[i, IA]
            r = op(r, f(val))
        end
        R[i1,IR] = r
    end
    R
end

a1 = zeros(Int, 1, 5, 5)
b = rand(5,5,5)
ba = AxisArray(b)

simple_reducedim!(t->t>.5, +, a1, ba)

There are some weird issues in the inference, hence a few of the temporary variables. Still no idea what's causing this, except that there appears to be an unexpected creation of an Axis{:dim_4}.

[timholy]

Curiously, the bug seems to have been introduced in #40. Strange that this hasn't manifested before now. Perhaps it was compensated elsewhere and is now exposed.

[Balinus]

Nice! Will try at work next monday.