Error when mapslices changes the array size

[sethaxen]

julia> using DimensionalData

julia> x = DimArray(randn(4, 100, 2), (:chain, :draw, :x_dim_1));

julia> y = mapslices(vec, x; dims=(:chain, :draw));

julia> size(y)
(400, 1, 2)

julia> DimensionalData.dims(y)
Dim{:chain} , Dim{:draw} , Dim{:x_dim_1} 

julia> y
400×1×2 DimArray{Float64,3} with dimensions: Dim{:chain} , Dim{:draw} , Dim{:x_dim_1} 
[:, :, 1]
Error showing value of type DimArray{Float64, 3, Tuple{Dim{:chain, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:draw, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:x_dim_1, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}:
ERROR: BoundsError: attempt to access 4-element Base.OneTo{Int64} at index [1:8]
Stacktrace:
  [1] throw_boundserror(A::Base.OneTo{Int64}, I::Tuple{UnitRange{Int64}})
    @ Base ./abstractarray.jl:691
  [2] checkbounds
    @ ./abstractarray.jl:656 [inlined]
  [3] getindex
    @ ./range.jl:918 [inlined]
  [4] getindex
    @ ~/.julia/packages/DimensionalData/1iOoH/src/LookupArrays/indexing.jl:8 [inlined]
  [5] getindex
    @ ~/.julia/packages/DimensionalData/1iOoH/src/Dimensions/indexing.jl:11 [inlined]
  [6] _slicedims
    @ ~/.julia/packages/DimensionalData/1iOoH/src/Dimensions/primitives.jl:418 [inlined]
  [7] _slicedims
    @ ~/.julia/packages/DimensionalData/1iOoH/src/Dimensions/primitives.jl:410 [inlined]
  [8] _slicedims
    @ ~/.julia/packages/DimensionalData/1iOoH/src/Dimensions/primitives.jl:402 [inlined]
  [9] slicedims
    @ ~/.julia/packages/DimensionalData/1iOoH/src/Dimensions/primitives.jl:382 [inlined]
 [10] rebuildsliced (repeats 2 times)
    @ ~/.julia/packages/DimensionalData/1iOoH/src/array/array.jl:56 [inlined]
 [11] getindex
    @ ~/.julia/packages/DimensionalData/1iOoH/src/array/indexing.jl:49 [inlined]
 [12] print_matrix(io::IOContext{Base.TTY}, A::DimArray{Float64, 2, Tuple{Dim{:chain, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:draw, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{Dim{:x_dim_1, DimensionalData.Dimensions.LookupArrays.NoLookup{UnitRange{Int64}}}}, SubArray{Float64, 2, Array{Float64, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata})
    @ DimensionalData ~/.julia/packages/DimensionalData/1iOoH/src/array/show.jl:87
 [13] print_array(io::IOContext{Base.TTY}, mime::MIME{Symbol("text/plain")}, A::DimArray{Float64, 3, Tuple{Dim{:chain, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:draw, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:x_dim_1, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata})
    @ DimensionalData ~/.julia/packages/DimensionalData/1iOoH/src/array/show.jl:51
 [14] show_after
    @ ~/.julia/packages/DimensionalData/1iOoH/src/array/show.jl:33 [inlined]
 [15] show(io::IOContext{Base.TTY}, mime::MIME{Symbol("text/plain")}, A::DimArray{Float64, 3, Tuple{Dim{:chain, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:draw, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Dim{:x_dim_1, DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata})
    @ DimensionalData ~/.julia/packages/DimensionalData/1iOoH/src/array/show.jl:25
 [16] (::REPL.var"#43#44"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:266
 [17] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:510
 [18] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:259
 [19] display(d::REPL.REPLDisplay, x::Any)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:271
 [20] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:328
 [21] #invokelatest#2
    @ ./essentials.jl:716 [inlined]
 [22] invokelatest
    @ ./essentials.jl:714 [inlined]
 [23] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:293
 [24] (::REPL.var"#45#46"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:277
 [25] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:510
 [26] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:275
 [27] (::REPL.var"#do_respond#66"{Bool, Bool, REPL.var"#77#87"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:846
 [28] (::VSCodeServer.var"#98#101"{REPL.var"#do_respond#66"{Bool, Bool, REPL.var"#77#87"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.6.30/scripts/packages/VSCodeServer/src/repl.jl:122
 [29] #invokelatest#2
    @ ./essentials.jl:716 [inlined]
 [30] invokelatest
    @ ./essentials.jl:714 [inlined]
 [31] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/LineEdit.jl:2493
 [32] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/.julia/juliaup/julia-1.7.3+0.x64/share/julia/stdlib/v1.7/REPL/src/REPL.jl:1232
 [33] (::REPL.var"#49#54"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:429

[rafaqz]

This isn't actually a problem with show, just where it surfaces. The issue is your mapslices resizes the axes of the array, but we assume the dimensions have the same axes as the array.

With your named NoLookup dimensions it shouldn't matter, but they do have a stored axis so other operations work, so in practice the vec resize breaks the relationship between the dimensions and the array.

To fix we could add a resize check to mapslices to allow this for NoLookup, and convert to NoLookup for other Lookup types, like Sampled, where the dimension no longer matches the array axis. But does it even make sense to have the same dimension names after you vec the slices?

[felixcremer]

I have a related problem, when I do a mean along one dimension I get a Raster with the same number of dimensions but the time dimension is reduced to a singleton dimension. I think in that case it make sense to keep the name of the dimension, but the lookup type could be changed and the length of the time dimension should also be reduced to one.
I always assumed, that the length of the dimensions and the size(raster, dim) need to be the same.

I can write down an explicit example of my analysis if this is helpful.

[rafaqz]

This seems like a solution. If the dims have changed we strip the lookups but keep the names, if not keep everything:

@inline function Base.mapslices(f, A::AbstractDimArray; dims=1, kw...)
    # Run `mapslices` on the parent array
    dimnums = dimnum(A, _astuple(dims))
    newdata = mapslices(f, parent(A); dims=dimnums, kw...)
    # Run one slice with dimensions to get the transformed dim
    d_inds = map(d -> rebuild(d, 1), DD.dims(A, _astuple(dims)))
    example_dims = DD.dims(f(view(A, d_inds...)))
    newdims = if isnothing(example_dims) || size(example_dims) != size(dims(A, d_inds))
        replacement_dims = map(d -> rebuild(d, NoLookup()), d_inds)
        newdims = format(setdims(DD.dims(A), replacement_dims), newdata)
    else
        example_dims
    end

    return rebuild(A, newdata, newdims)
end

Theres a bit of type instability but it works fine:

julia> y = mapslices(vec, x; dims=(:chain, :draw))
400×1×2 DimArray{Float64,3} with dimensions: Dim{:chain}, Dim{:draw}, Dim{:x_dim_1}
[:, :, 1]
 -0.82393
  0.0772902
 -0.143477
 -0.820947
 -0.995829
  0.0400354
  0.520451
  0.696105
  0.322051
 -1.12305
  0.777396
  ⋮
 -0.073463
  2.46297
 -1.40587
 -0.936701
  0.691183
  0.991734
 -0.757005
  1.49445
  0.857546
 -1.05538
 -0.0514016
[and 19 more slices...]