ComponentVectors are not stack-able

[gdalle]

Expected behavior:

julia> x = [1, 2];

julia> y = [3, 4];

julia> hcat(x, y)
2×2 Matrix{Int64}:
 1  3
 2  4

julia> stack([x, y])
2×2 Matrix{Int64}:
 1  3
 2  4

Observed behavior:

julia> x = ComponentVector(a=[1, 2]);

julia> y = ComponentVector(a=[3, 4]);

julia> hcat(x, y)
2×2 ComponentMatrix{Int64} with axes Axis(a = 1:2,) × FlatAxis()
 1  3
 2  4

julia> reduce(hcat, [x,y])  # different type for some reason
2×2 Matrix{Int64}:
 1  3
 2  4

julia> stack([x, y])
ERROR: BoundsError: attempt to access 2×1 ComponentMatrix{Int64, Matrix{Int64}, Tuple{Axis{(a = 1:2,)}, FlatAxis}} with indices 1:1:2×1:1:1 at index [3:4]
Stacktrace:
  [1] copyto!(dest::ComponentMatrix{…}, dstart::Int64, src::ComponentVector{…}, sstart::Int64, n::Int64)
    @ Base ./abstractarray.jl:1134
  [2] copyto!
    @ ./abstractarray.jl:1118 [inlined]
  [3] _typed_stack(::Colon, ::Type{Int64}, ::Type{ComponentVector{…}}, A::Vector{ComponentVector{…}}, Aax::Tuple{Base.OneTo{…}})
    @ Base ./abstractarray.jl:2827
  [4] _typed_stack
    @ ./abstractarray.jl:2817 [inlined]
  [5] _stack
    @ ./abstractarray.jl:2807 [inlined]
  [6] _stack
    @ ./abstractarray.jl:2799 [inlined]
  [7] stack(iter::Vector{ComponentVector{Int64, Vector{Int64}, Tuple{Axis{(a = 1:2,)}}}})
    @ Base ./abstractarray.jl:2767
  [8] top-level scope

[jonniedie]

Oh interesting. I'll try to take a look at it soon. Unfortunately it probably won't be this week as I'm a little swamped with work stuff. In the meantime, I think the old reduce(hcat, [x, y]) will work, it's probably just not as efficient. Thanks for finding this!

[gdalle]

Thanks for your answer @jonniedie! I tried to help by digging a little more, and I think the problem lies with the definition of similar for ComponentArrays. Here's what happens at the beginning of the problematic function _typed_stack in the above MWE:

julia> ax1 = Base._iterator_axes(x)
(1:1:2,)

julia> Aax = Base._iterator_axes([x, y])
(Base.OneTo(2),)

julia> B = similar(x, ax1..., Aax...)  # should be 2x2
2×1 ComponentMatrix{Int64} with axes Axis(a = 1:2,) × FlatAxis()
 138725991131840
 138725991140464

In other words, the created matrix has the wrong number of columns. Note that if instead we do

julia> B = similar(x, Aax..., Aax...)
2×2 Matrix{Int64}:
 1  138723500513936
 2  138723500514000

then the construction is correct. Unfortunately I am struggling to go beyond this diagnosis, but hopefully this is useful for you to narrow down the issue.

[dingraha]

The MWE @gdalle posted here works now on this PR: https://github.com/jonniedie/ComponentArrays.jl/pull/249, but the behavior is such that it ignores the names of the components in the second array (y in the MWE):

julia> x = ComponentVector(a=[1, 2])
ComponentVector{Int64}(a = [1, 2])

julia> y = ComponentVector(b=[3, 4])
ComponentVector{Int64}(b = [3, 4])

julia> xy = stack([x, y])
2×2 ComponentMatrix{Int64} with axes Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),) × FlatAxis()
 1  3
 2  4

julia> xy[:a, 1]
2-element Vector{Int64}:
 1
 2

julia> xy[:a, 2]
2-element Vector{Int64}:
 3
 4

julia> 

Similar for 3 vectors, only the first argument's axis matters:

julia> z = ComponentVector(c = [5, 6])
ComponentVector{Int64}(c = [5, 6])

julia> xyz = stack([x, y, z])
2×3 ComponentMatrix{Int64} with axes Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),) × FlatAxis()
 1  3  5
 2  4  6

julia> xyz[:a, 1]
2-element Vector{Int64}:
 1
 2

julia> xyz[:a, 2]
2-element Vector{Int64}:
 3
 4

julia> xyz[:a, 3]
2-element Vector{Int64}:
 5
 6

julia> 

[gdalle]

To me this is much better than erroring, and arguably you wanna stack ComponentVectors with the same structure anyway. Thank you! Hope your PR gets merged, if there's anything I can do to help either you or @jonniedie let me know

[mcabbott]

Similar for 3 vectors, only the first argument's axis matters:

FWIW, stack demands that axes match, for all the inner objects. This property seems somewhat surprising to me:

julia> axes(x) == axes(z)
true

but if you have reasons to have it that way, then stack((x,y,z)) should probably also work without error. Presumably any of axes(x) or axes(z) would be acceptable?

Alternatively you could decide that, in mismatched cases, it's better not to give an unpredictable fancy axis... follow this instead:

julia> hcat(x,x)
2×2 ComponentMatrix{Int64} with axes Axis(a = 1:2,) × FlatAxis()
 1  1
 2  2

julia> hcat(x,z)
2×2 Matrix{Int64}:
 1  5
 2  6

For testing things, note also that stack accepts iterators of iterators, and hence these should probably all work:

stack((x,y,x); dims=1)  # tuple of arrays
stack(ComponentArray(a=(1,2,3), b=(4,5,6)))  # array of tuples
stack(ComponentArray(z=[x,x]) for _ in 1:4)  # generator of arrays
stack(x -> ComponentArray(a=x, b=[x,x]), [5 6; 7 8]; dims=2)  # map then stack

[gdalle]

@dingraha do you think your changes are compatible with what @mcabbott mentioned?

[dingraha]

I'm actually not sure what the result of those are supposed to be 😬. But I tried each of them. Here's what I get:

First one:

julia> x = ComponentVector(a=[1, 2])
ComponentVector{Int64}(a = [1, 2])

julia> y = ComponentVector(b=[3, 4])
ComponentVector{Int64}(b = [3, 4])

julia> x_noca = getdata(x)
2-element Vector{Int64}:
 1
 2

julia> y_noca = getdata(y)
2-element Vector{Int64}:
 3
 4

julia> stack((x, y, x); dims=1)
3×2 ComponentMatrix{Int64} with axes FlatAxis() × Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),)
 1  2
 3  4
 1  2

julia> stack((x_noca, y_noca, x_noca); dims=1)
3×2 Matrix{Int64}:
 1  2
 3  4
 1  2

julia> 

I think that's what we'd want.

Second one:

julia> stack(ComponentArray(a=(1,2,3), b=(4,5,6)))
3×2 ComponentMatrix{Int64} with axes FlatAxis() × Axis(a = 1, b = 2)
 1  4
 2  5
 3  6

julia> 

I think that does what's intended. It's the same as stacking a tuple of two vectors:

julia> stack((x, y))
2×2 ComponentMatrix{Int64} with axes Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),) × FlatAxis()
 1  3
 2  4

julia> 

Third one:

julia> x
ComponentVector{Int64}(a = [1, 2])

julia> stack(ComponentArray(z=[x,x]) for _ in 1:4)
4×4 ComponentMatrix{Int64} with axes Axis(z = ViewAxis(1:4, PartitionedAxis(2, Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),))),) × FlatAx
is()
 1  1  1  1
 2  2  2  2
 1  1  1  1
 2  2  2  2

julia>

Not sure about that one.

Fourth one:

julia> x
ComponentVector{Int64}(a = [1, 2])

julia> stack(x -> ComponentArray(a=x, b=[x,x]), [5 6; 7 8]; dims=2)
3×4 ComponentMatrix{Int64} with axes Axis(a = 1, b = ViewAxis(2:3, Shaped1DAxis((2,)))) × FlatAxis()
 5  7  6  8
 5  7  6  8
 5  7  6  8

julia> 

Also not sure if that's what's intended.

[mcabbott]

These look fine -- or at least they look like they contain the right Matrix, ignoring special axis types.

In the 4th, x is a local variable (sorry), hence I think it ought to match this:

julia> stack(x -> [x,x,x], [5 6; 7 8]; dims=2)
3×4 Matrix{Int64}:
 5  7  6  8
 5  7  6  8
 5  7  6  8

Perhaps without dims would be a better test, as no earlier case returns an array with ndims > 2:

 julia> stack(x -> [x,x,x], [5 6; 7 8])
3×2×2 Array{Int64, 3}:
[:, :, 1] =
 5  7
 5  7
 5  7

[:, :, 2] =
 6  8
 6  8
 6  8

[dingraha]

@mcabbott Excellent, thanks for the help. Here's what I'm getting:

julia> x2
ERROR: UndefVarError: `x2` not defined

julia> stack(x2 -> [x2,x2,x2], [5 6; 7 8]; dims=2)
3×4 Matrix{Int64}:
 5  7  6  8
 5  7  6  8
 5  7  6  8

julia> stack(x2 -> ComponentArray(a=x2, b=[x2,x2]), [5 6; 7 8]; dims=2)  # map then stack
3×4 ComponentMatrix{Int64} with axes Axis(a = 1, b = ViewAxis(2:3, Shaped1DAxis((2,)))) × FlatAxis()
 5  7  6  8
 5  7  6  8
 5  7  6  8

julia> stack(x2 -> [x2,x2,x2], [5 6; 7 8])
3×2×2 Array{Int64, 3}:
[:, :, 1] =
 5  7
 5  7
 5  7

[:, :, 2] =
 6  8
 6  8
 6  8

julia> stack(x2 -> ComponentArray(a=x2, b=[x2,x2]), [5 6; 7 8])
3×2×2 ComponentArray{Int64, 3, Array{Int64, 3}, Tuple{Axis{(a = 1, b = ViewAxis(2:3, Shaped1DAxis((2,))))}, FlatAxis, FlatAxis}} with in
dices 1:1:3×1:1:2×1:1:2:
[:, :, 1] =
 5  7
 5  7
 5  7

[:, :, 2] =
 6  8
 6  8
 6  8

julia> 

[dingraha]

I should have actually tried to index the component arrays with components. Here are the four examples, but this time I try indexing them with symbols. Sadly the third one doesn't work. :-(

julia> Xstack1 = stack((x, y, x); dims=1)
3×2 ComponentMatrix{Int64} with axes FlatAxis() × Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),)
 1  2
 3  4
 1  2

julia> Xstack1[:, :a]
3×2 Matrix{Int64}:
 1  2
 3  4
 1  2

julia> Xstack2 = stack(ComponentArray(a=(1,2,3), b=(4,5,6)))
3×2 ComponentMatrix{Int64} with axes FlatAxis() × Axis(a = 1, b = 2)
 1  4
 2  5
 3  6

julia> Xstack2[:, :a]
3-element Vector{Int64}:
 1
 2
 3

julia> Xstack2[:, :b]
3-element Vector{Int64}:
 4
 5
 6

julia> Xstack3 = stack(ComponentArray(z=[x,x]) for _ in 1:4)
4×4 ComponentMatrix{Int64} with axes Axis(z = ViewAxis(1:4, PartitionedAxis(2, Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),))),) × FlatAx
is()
 1  1  1  1
 2  2  2  2
 1  1  1  1
 2  2  2  2

julia> getaxes(Xstack3)
(Axis(z = ViewAxis(1:4, PartitionedAxis(2, Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),))),), FlatAxis())

julia> Xstack3[:z, :]
Error showing value of type ComponentArrays.LazyArray{Any, 1, Base.Generator{Base.Generator{StepRange{Int64, Int64}, ComponentArrays.var
"#1#2"{Matrix{Int64}, Int64}}, ComponentArrays.var"#18#19"{Tuple{Axis{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}}}}:
ERROR: type NamedTuple has no field a
Stacktrace:
  [1] getproperty
    @ ./Base.jl:37 [inlined]
  [2] getindex(#unused#::FlatAxis, s::Symbol)
    @ ComponentArrays ~/projects/componentarrays_bugs/dev/ComponentArrays/src/axis.jl:169
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [5] (::Base.Broadcast.var"#31#32"{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(getindex), Tuple{Tuple{Axis{
(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}, Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(ComponentA
rrays.getval), Tuple{Tuple{DataType}}}}}})(k::Int64)
    @ Base.Broadcast ./broadcast.jl:1088
  [6] ntuple
    @ ./ntuple.jl:49 [inlined]
  [7] copy
    @ ./broadcast.jl:1088 [inlined]
  [8] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(getindex), Tuple{Tuple{Axis{(a = ViewAxis(
1:2, Shaped1DAxis((2,))),)}, FlatAxis}, Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, typeof(ComponentArrays.getval),
 Tuple{Tuple{DataType}}}}})
    @ Base.Broadcast ./broadcast.jl:873
  [9] #s40#57
    @ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:120 [inlined]
 [10] var"#s40#57"(::Any, index_fun::Any, x::Any, idx::Any)
    @ ComponentArrays ./none:0
 [11] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any})
    @ Core ./boot.jl:602
 [12] getindex
    @ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:103 [inlined]
 [13] getindex
    @ ~/projects/componentarrays_bugs/dev/ComponentArrays/src/array_interface.jl:102 [inlined]
 [14] show(io::IOContext{IOBuffer}, x::ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}, Tuple{Ax
is{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}})
    @ ComponentArrays ~/projects/componentarrays_bugs/dev/ComponentArrays/src/show.jl:52
 [15] sprint(f::Function, args::ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}, Tuple{Axis{(a =
 ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}}; context::IOContext{IOBuffer}, sizehint::Int64)
    @ Base ./strings/io.jl:112
 [16] sprint
    @ ./strings/io.jl:107 [inlined]
 [17] alignment_from_show
    @ ./show.jl:2817 [inlined]
 [18] alignment(io::IOContext{IOBuffer}, x::ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}, Tup
le{Axis{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}})
    @ Base ./show.jl:2836
 [19] alignment(io::IOContext{IOBuffer}, X::AbstractVecOrMat, rows::Vector{Int64}, cols::Vector{Int64}, cols_if_complete::Int64, cols_ot
herwise::Int64, sep::Int64, ncols::Int64)
    @ Base ./arrayshow.jl:69
 [20] _print_matrix(io::IOContext{IOBuffer}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, 
ddots::String, hmod::Int64, vmod::Int64, rowsA::UnitRange{Int64}, colsA::UnitRange{Int64})
    @ Base ./arrayshow.jl:207
 [21] print_matrix(io::IOContext{IOBuffer}, X::Vector{ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, 
true}, Tuple{Axis{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}}}, pre::String, sep::String, post::String, hdots::String, vdots::
String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base ./arrayshow.jl:171
 [22] print_matrix
    @ ./arrayshow.jl:171 [inlined]
 [23] print_array
    @ ./arrayshow.jl:358 [inlined]
 [24] show(io::IOBuffer, #unused#::MIME{Symbol("text/plain")}, X::Vector{ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{
UnitRange{Int64}}, true}, Tuple{Axis{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}}})
    @ Base ./arrayshow.jl:399
 [25] __binrepr(m::MIME{Symbol("text/plain")}, x::Vector{ComponentVector{Int64, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}
}, true}, Tuple{Axis{(a = ViewAxis(1:2, Shaped1DAxis((2,))),)}, FlatAxis}}}, context::Nothing)
    @ Base.Multimedia ./multimedia.jl:171
 [26] _binrepr
    @ ./multimedia.jl:0 [inlined]
 [27] #repr#1
    @ ./multimedia.jl:159 [inlined]
 [28] repr
    @ ./multimedia.jl:159 [inlined]
 [29] show(io::IOContext{Base.TTY}, mime::MIME{Symbol("text/plain")}, a::ComponentArrays.LazyArray{Any, 1, Base.Generator{Base.Generator
{StepRange{Int64, Int64}, ComponentArrays.var"#1#2"{Matrix{Int64}, Int64}}, ComponentArrays.var"#18#19"{Tuple{Axis{(a = ViewAxis(1:2, Sh
aped1DAxis((2,))),)}, FlatAxis}}}})
    @ ComponentArrays ~/projects/componentarrays_bugs/dev/ComponentArrays/src/lazyarray.jl:43
 [30] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:276
 [31] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:557
 [32] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:262
 [33] display
    @ ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:281 [inlined]
 [34] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:340
 [35] #invokelatest#2
    @ ./essentials.jl:819 [inlined]
 [36] invokelatest
    @ ./essentials.jl:816 [inlined]
 [37] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:305
 [38] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:287
 [39] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:557
 [40] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:285
 [41] (::REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.Li
neEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:899
 [42] #invokelatest#2
    @ ./essentials.jl:819 [inlined]
 [43] invokelatest
    @ ./essentials.jl:816 [inlined]
 [44] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/LineEdit.jl:2647
 [45] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/local/julia/1.9.4/share/julia/stdlib/v1.9/REPL/src/REPL.jl:1300
 [46] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:514

julia> Xstack3[:, :]
4×4 ComponentMatrix{Int64} with axes Axis(z = ViewAxis(1:4, PartitionedAxis(2, Axis(a = ViewAxis(1:2, Shaped1DAxis((2,))),))),) × FlatAx
is()
 1  1  1  1
 2  2  2  2
 1  1  1  1
 2  2  2  2

julia> Xstack4 = stack(x -> ComponentArray(a=x, b=[x,x]), [5 6; 7 8])
3×2×2 ComponentArray{Int64, 3, Array{Int64, 3}, Tuple{Axis{(a = 1, b = ViewAxis(2:3, Shaped1DAxis((2,))))}, FlatAxis, FlatAxis}} with in
dices 1:1:3×1:1:2×1:1:2:
[:, :, 1] =
 5  7
 5  7
 5  7

[:, :, 2] =
 6  8
 6  8
 6  8

julia> Xstack4[:a, :, :]
2×2 Matrix{Int64}:
 5  6
 7  8

julia> Xstack4[:b, :, :]
2×2×2 Array{Int64, 3}:
[:, :, 1] =
 5  7
 5  7

[:, :, 2] =
 6  8
 6  8

julia>