Add stack(array_of_arrays)

[mcabbott]

This adds a simple version of the function stack described in #21672. It generalises reduce(hcat, vector_of_vectors) to handle more dimensions, and to handle iterators efficiently.

For iterators, like #31644 this handles only cases for which repeated append! is the fallback, except that it will also widen the type if necessary. You would need quite a different implementation to handle things like reduce(vcat, ([1 i i^2] for i in 1:10 if i>pi)), probably collecting first to know where to write the first slice. Because of that, it does not take a dims keyword. Thus it can be the inverse of eachcol(::Matrix) or eachslice(::Array{T,3}; dims=3) but not of eachrow.

Edit: The current implementation introduced in this message below takes a dims keyword. With this, stack(eachslice(A; dims); dims) == A. Without the dims keyword, the dimensions of the slice always come first, thus it should obey vec(stack(xs)) == collect(Iterators.flatten(xs)). Iterators of unknown length are collected first, as this seems quicker than falling back to append!-and-reshape.

Handling iterators also makes it trivial to allow stack(f, xs) = stack(f(x) for x in xs), as suggested below. This seems nice to have with a do block.

It does not currently handle the empty case , as I'm not sure what it should do. Ideally stack(empty!([1:5])) would make a 5×0 Matrix, but that's impossible. At present it makes a 1×0 Matrix. The empty handling might be more complicated than ideal, see below. Reverted to an error because handling weird array types is hard, and my attempt produced type instabilities.

Should work fine with OffsetArrays, it preserves all axes. Should mostly work with CUDA.

The name unfortunately clashes with DataFrames.stack. Other names suggested in #21672 include pack or nested; perhaps seeing it as the inverse of eachslice might suggest others?

[antoine-levitt]

As I stated in #21672, a concern is that stack([1,2],[3,4]) could reasonably mean [1,2,3,4] as well as [1 3;2 4]. pack actually suggests the first interpretation to me; same for flatten. nested feels like the opposite operation to stack. unnest is ugly but might convey the intention better? Or matrixify? Matrix(::AbstractVector) is free, but it would feel weird if the result was not a Matrix? Not saying that it should block the PR (it's a useful functionality that deserves a name, even if that name is not unambiguous), but worth thinking over in case there's a better alternative.

[mcabbott]

Re names, first English: If you want to "stack" some papers, or some hay, then you don't put the pieces together on their long dimensions as that would fall over. Maybe it becomes less clear if you stack textbooks. Whereas "concatenate" seems to come from liturgy, to mean what it does for strings, or architecture, like hcat of matrices. To me what "pack" suggest is that it fills the output completely, unlike e.g. cat(1,2,3; dims=(1,2)); maybe it leans towards new-axes?

Python's numpy.stack has this meaning, and explicitly distinguishes this "new axes" vs. cat's "same axes". In Mathematica there is no such operation, multi-dimensional arrays are lists of lists. It doesn't look like Matlab has such a function (separate from cat). I guess hcat([1,2], [3,4]) makes most sense if you think of these as really being 2×1 matrices.

In Julia, if this is the inverse of eachslice, it could be unslice? I think this perspective might make what it does clearer: Nobody expects eachcol to cut its input mid-way along a column, the dimensions are partitioned into those of the elements and that of the container. JuliennedArrays calls this Align, but perhaps that makes more sense for a view. There is also Iterators.flatten which always collects to a vector.

My objection to matrixify is that it ties this to 2 dimensions, which is where you need it least. It can't be a method of Array since [1:2, 3:4] isa Array, and it may make a CuArray. But something like dense (it makes a DenseArray, probably), or "solid" or "continguous"?

[simeonschaub]

Could we also add a two-arg version with a function argument? (Could just be implemented as stack(Iterators.map(f, x)).) I find myself using the mapreduce(f, vcat, x) pattern quite often, but unfortunately that's currently really slow.

[jonas-schulze]

Maybe concat (or even concatenate) would be in good analogy to e.g. max and maximum, as maximum(x) = max(x...). Here, concat(v, dims=dims) = cat(v..., dims=dims), and dims=ndims(first(v))+1 could be the default if no dims is present.

[mcabbott]

I'm sure I've suggested this max/maximum analogy too somewhere. Arguments against it, and in favour of a new verb, are:

• One, if M is a matrix, this is the inverse of eachcol in that stack(eachslice(M, dims=2)) == M. If a future version took It takes a dims keyword, then a plausible meaning would be which allows stack(eachrow(M); dims=1) == M. That can be useful (it's what Flux.stack does). But it's not the meaning of vcat, which would put the vectors of eachrow end-to-end.

• Two, the column-major append! behaviour very naturally handles more dimensions, including container dimensions. When add Slices array type for eachslice/eachrow/eachcol #32310 is eventually merged Now, eachslice(A4; dims=(3,4)) will be is a matrix of matrices, which stack puts back together. This also doesn't fit, cat(mats...; dims=(3,4)) instead some kind of block-diagonal arrangement.

• [Edit] Three, this implementation traverses non-array elements, stack(((1,2), (x=3, y=4))) is quite unlike hcat((1,2), (x=3, y=4)) or any other cat function. I think that a function with a name closely connected to cat should not differ on this. And that this is a nice bonus, rather than an essential feature.

• [Edit] Related to one, I think there's clarity in the idea that cat always extends existing axes (possibly trivial, size([1,2,3],4) == 1) while stack always places arrays along new axes. This is the same distinction made in Python, mentioned above. It is also why stack has obvious rules for OffsetArrays.

[Edit] Since this thread is now very long, I've enlarged this slightly. There are examples below where the point was raised again.

[mcabbott]

I added the stack(f, iter) method suggested.

Notice that this allows f to return tuples, and treats them just like vectors. This is not how cat behaves, but if it's a different verb then perhaps we have more freedom. The reason to want this is that, if your function always returns (say) 3-vectors, then this might let you avoid allocating many small vectors.

You could go further and make stack([(1,2), (3,4)]) into reinterpret(reshape, Int, [(1,2), (3,4)]). But that's perhaps a surprising return type.

Edit: The analogy to Iterators.flatten is another reason to think tuples should be allowed. Having vec(stack(iter)) == collect(Iterators.flatten(iter)) seems like a simple and consistent rule.

[mcabbott]

However, this has type-stability problems right now? No, I was testing wrong. Corrected version:

julia> x = randn(3, 1000);

julia> @btime stack(eachcol($x));
  min 4.774 μs, mean 7.023 μs (4 allocations, 23.58 KiB. GC mean 12.95%)

julia> ft(x) = @inbounds (x[1], x[2]/2, x[3]^3);
julia> fv(x) = @inbounds [x[1], x[2]/2, x[3]^3];

julia> @btime stack(fv, eachcol($x));
  min 21.125 μs, mean 25.090 μs (1004 allocations, 101.70 KiB)

julia> @btime stack(ft, eachcol($x));  # making tuples not vectors
  min 4.768 μs, mean 6.436 μs (4 allocations, 23.58 KiB)
  
julia> @btime mapslices(fv, $x, dims=1);  # alternative (although this can be improved)
  min 222.333 μs, mean 237.788 μs (7502 allocations, 266.73 KiB)

[longemen3000]

There is also PackedVectorsOfVectors.pack

[mschauer]

Ref #16708

[mcabbott]

Thanks, I hadn't seen that iterator approach. Perhaps something like that could now be written as mixture of Iterators.flatten (same iteration) and Iterators.product (same shape, given the first element).

But can it be fast? It seems that flatten is a lot slower than reduce(vcat, xs):

julia> xs = [rand(128) for _ in 1:100];

julia> v = @btime collect(Iterators.flatten($xs));
  min 55.750 μs, mean 64.555 μs (9 allocations, 326.55 KiB)

julia> v ≈ @btime collect(reduce(vcat, $xs))
  min 4.542 μs, mean 13.543 μs (5 allocations, 200.97 KiB)
true

[mschauer]

Yes, it is a tricky thing, that's perhaps why it didn't went anywhere. On the other hand

julia> using StaticArrays

julia> xs = [@SVector rand(3) for _ in 1:100];

julia> v = @btime collect(Iterators.flatten($xs));
  592.251 ns (1 allocation: 2.50 KiB)

julia> v ≈ @btime collect(reduce(vcat, $xs))
  861.942 ns (3 allocations: 5.88 KiB)

so reduce(hcat,...) isn't the last word either...

[mcabbott]

Indeed, that's pretty fast. (Notice we're penalising reduce(hcat, xs) by collecting twice, my mistake.)

The analogy to Iterators.flatten is perhaps a useful way to think about this. If it collects an iterator of iterators, then the slices need not be arrays. Which means that something like this can avoid allocating them, for an efficient small-slice mapslices:

stack(xs, eachcol(Y)) do x, y
    x/y[1], y[2]/y[3], y[4]/x
end

Allowing this seems like an advantage of not calling this *cat, since that has other rules: hcat([1,2], [3,4]) != hcat((1,2), (3,4)).

I think I have a better implementation now, which doesn't make an iterator but does accept them. It deletes all the code dealing with append! and type-widening, instead requiring that it knows the size and type up front, if neccessary by collecting. On this example:

julia> xs = [@SVector rand(3) for _ in 1:100];

julia> m = @btime reduce(hcat, $xs);
  min 539.243 ns, mean 843.973 ns (1 allocation, 2.50 KiB)

julia> m ≈ @btime reshape(collect(Iterators.flatten($xs)), 3, :)
  min 344.329 ns, mean 637.694 ns (3 allocations, 2.59 KiB)
true

julia> m ≈ @btime reshape(collect(Iterators.flatten(Tuple(x) for x in xs)), 3, :)
  min 1.492 μs, mean 2.170 μs (10 allocations, 7.73 KiB)
true

julia> m ≈ @btime stack(Tuple(x) for x in $xs)
  min 123.719 ns, mean 417.366 ns (1 allocation, 2.50 KiB)
true

Still a bit WIP, am not entirely sure whether to push to this PR or make a separate one.
It also wants to allow dims, but what's still slow is stackT(Tuple(x) for x in xs; dims=1).

[JeffBezanson]

Triage likes this. The name situation with DataFrames is tricky, but we do think stack is the canonical name so it would be good to go with that.

[MasonProtter]

FWIW, DataFrames.stack only has methods for AbstractDataFrame, so it could be done to just share the name with Base.

julia> using DataFrames

julia> methods(stack)
# 3 methods for generic function "stack":
[1] stack(df::AbstractDataFrame) in DataFrames at /Users/mason/.julia/packages/DataFrames/ORSVA/src/abstractdataframe/reshape.jl:134
[2] stack(df::AbstractDataFrame, measure_vars) in DataFrames at /Users/mason/.julia/packages/DataFrames/ORSVA/src/abstractdataframe/reshape.jl:134
[3] stack(df::AbstractDataFrame, measure_vars, id_vars; variable_name, value_name, view, variable_eltype) in DataFrames at /Users/mason/.julia/packages/DataFrames/ORSVA/src/abstractdataframe/reshape.jl:134

[mcabbott]

OK, I've replaced the implementation here with one which is I think better in every way, except brevity. It will never call append!, as for iterators of unknown length producing arrays, it seems almost always faster to collect, and always knowing the size is simpler.

Its goals are:

• Accept an iterator of iterators, not necessarily an array of arrays.
  Partly because stack(f, A) with a do block sounds like a useful thing.
• Be as fast as reduce(hcat, A) on a vector of vectors, where they agree.
  This seems to be possible without indexing A, only iterating.
• Take a dims keyword, such that stack(eachslice(A; dims); dims) == A.
  Also lines up nicely with selectdim. This must be an integer (in all 3 functions).
• Make OffsetArrays work automatically, by preserving axes.
  Unlike cat, this never needs to put existing axes end-to-end.

I've attached some comments on the implementation below.

[mcabbott]

This implementation accepts iterators, but demands to know the size and element type before proceeding. If these are unknown, or abstract, it will collect.

If the elements (the slices being stacked) themselves have abstract eltype, then so will the answer. Which matches this:

julia> Iterators.flatten(Any[[1,2], [3,4]]) |> collect
4-element Vector{Int64}:
 1
 2
 3
 4

julia> Iterators.flatten([Any[1,2], Any[3,4]]) |> collect
4-element Vector{Any}:
 1
 2
 3
 4

I tried some versions which instead promoted the types within individual slices to make the final output concrete. But these were always much slower. Perhaps getting back a Matrix{Real} is better, a visible sign of a problem you should fix, rather than an invisible slow path.

I thought cat always promoted, but in fact it varies:

julia> vcat(Any[1,2], Any[3,4])
4-element Vector{Any}:
 1
 2
 3
 4

julia> hcat(Any[1,2], Any[3,4])
2×2 Matrix{Int64}:
 1  3
 2  4

julia> cat(Any[1,2], Any[3,4]; dims=2)
2×2 Matrix{Any}:
 1  3
 2  4

[mcabbott]

Here _axes(x) generalises axes(x) to regard HasLength as being like HasShape{1}. And I've got what's essentially the same thing for ndims.

Are these safe? Can a HasLength iterator in fact collect to have multiple dimensions?

And if they are safe, can they just be made methods of axes and ndims globally?

[mcabbott]

You could do something like this globally:

Base.size(iter) = _size(IteratorSize(iter), iter)
_size(::HasLength, iter) = (length(iter),)
_size(::HasShape, iter) = map(length, axes(iter))

One reason not to is that types which are iterable but scalar under broadcast have two notions of size. Choosing this one breaks the link between axes and broadcast:

julia> axes("str")  # previously an error
(Base.OneTo(3),)

julia> axes(collect("str"))  # matches this
(Base.OneTo(3),)

julia> "str" .* '+'  # scalar under broadcasting
"str+"

julia> collect("str") .* '+' 
3-element Vector{String}:
 "s+"
 "t+"
 "r+"

[mcabbott]

This forbids e.g. stack([1:3, 2:4]; dims=3), which you might hope would make a 3×1×2 Array like cat([1:3, 2:4]...; dims=3). Doing so lets it be type-stable.

[mcabbott]

The output array is always B = similar(_prototype(x1, A), T, axes...), with the aim of making an array of CuArrays just work. Perhaps not worth the complexity and better to handle there?

Even without that, using say similar(1:0, T, axes...) instead of Array{T}(undef, ...) lets this work with OffsetArrays automatically. Preserving axes not size seems like the right thing to do.

[mcabbott]

Some empty behaviours, most to least specific:

julia> stack(Tuple{Int,Int,Int}[])
3×0 Matrix{Int64}

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

julia> stack(Matrix{Int}[])
1×1×0 Array{Int64, 3}

julia> stack([[1 2; 3 4]])
2×2×1 Array{Int64, 3}:
[:, :, 1] =
 1  2
 3  4

julia> stack(Matrix[])  # eltype is abstract, perhaps not ideal?
Any[]

julia> stack([])
Any[]

julia> stack([1,2,3])  # because numbers are iterable
3-element Vector{Int64}:
 1
 2
 3

My thinking re the Tuple handling is partly that SVectors could share these paths. It's possible I haven't thought enough about SMatrices.

[mcabbott]

Some comments here below about empty behaviour inducing type-instabilities. Even if the collection type is known, without an instance unfortunately you can't get similar to make one for you, e.g.:

julia> similar([1,2,3], (4,5)) |> summary
"4×5 Matrix{Int64}"

julia> similar(typeof([1,2,3]), (4,5)) |> summary
ERROR: MethodError: no method matching Vector{Int64}(::UndefInitializer, ::Tuple{Int64, Int64})

It's possible that empty cases should all (or mostly?) be made errors. This is what reduce(hcat, Vector{Int}[]) does.

[mcabbott]

One-month bump?

[mcabbott]

Should I read anything into the test failures? They appear the same on re-running, and don't seem to appear on other recent PRs. But they also look unrelated to me, e.g.:

:mac: test x86_64-apple-darwin

stdlib/v1.9/InteractiveUtils/test/runtests.jl:252
--
  | Test threw exception
  | Expression: !(occursin("Environment:", read(setenv(`$exename -e 'using InteractiveUtils; versioninfo()'`, String[]), String)))

:linux: test powerpc64le-linux-gnu

stdlib/v1.9/LinearAlgebra/test/qr.jl:426
--
  | Got exception outside of a @test
  | BoundsError: attempt to access 10×4 view(::Matrix{Float16}, 1:10, 2:5) with eltype Float16 at index [2:10, 2]
  | Stacktrace:
  | [1] throw_boundserror(A::SubArray{Float16, 2, Matrix{Float16}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, I::Tuple{UnitRange{Int64}, Int64})

[N5N3]

They should be unrelated, as we are pretty sure that stack is not used in Base and InteractiveUtils.
I think this PR is ready to merge. Not sure will there be other folks want to take a look.

Triage likes this. The name situation with DataFrames is tricky, but we do think stack is the canonical name so it would be good to go with that.

Let's merge this at this weekend if there's no other objections.

[raphaelchinchilla]

I have been waiting for a function like this for years. Thank you @mcabbott and the Julia community at large, you all are amazing

[nlw0]

I wouldn't say I object to this, but I feel we have some conversations going on in #21672 and in #45985, #46003 and #44792. The conversations didn't really reach any form of transcendental synthesis or synergistic consensus, in my opinion. A conclusion was apparently reached, though, since we're going ahead and merging this PR. It feels like we're in a kind of patch-merging race, there are different ideas for how this feature, important one in my opinion, should be brought to Base, but conversations don't really go anywhere and yet code is getting merged. I feel a bit lost as a contributor, what is the way to go?

[mcabbott]

This may not be the end of the story. stack is no help for vcat(vec_of_vecs...) which is a common desire, and I think there's some agreement that the reduce(vcat, vv) pattern isn't great.

One possibility is something like flat(x) == collect(Iterators.flatten(x)) (which always makes a vector). Another would be something like the function you call blockstack (which concatenates inner dims along outer dims -- perhaps seen as generalising hvcat). Making a prototype somewhere (and figuring out how to make it fast, what to do about OffsetArrays, how to clearly describe precisely what it does) isn't a bad idea.

[aplavin]

One possibility is something like flat(x) == collect(Iterators.flatten(x)) (which always makes a vector).

Please don't :) It's clearly feasible to preserve the array type and eltype better than that.

I've put flatten and flatmap that keep types (and infer) as far as possible - generally and performantly - into FlexiMaps.jl: https://gitlab.com/aplavin/FlexiMaps.jl.
And there's no need to create new function names given that flatten/flatmap are already present in Iterators in Julia, and are common in other languages as well.

I'd be glad to have this functionality in base, but won't pursue it myself - a package works just as well in practice, and much easier/faster to get done.

[mcabbott]

Yes I don't take any position on precise naming (sometimes Base & Iterators agree, sometimes they differ, not sure why) nor on return types (vector not Vector).

At present collect(Iterators.flatten(x)) is called by some comprehensions, and is quite slow, so one goal would be to speed that up without breaking anything. Perhaps independently of adding a new verb.

Perhaps someone should make an issue though, since this isn't the right place, stack is done.

[nlw0]

There's an open issue already, as I have just mentioned. It didn't make a difference.

I feel the concerns raised by me and others were brushed away with not enough attention, with little to no participation of core developers. In the meantime I have tried to offer alternative working code, to enrich the conversation and illustrate my points concretely, and I even did a short conference talk that was actually related to the subject. It was all completely ineffective, I guess I didn't do enough.

I feel this function lacks a proper semantic understanding, and I'm not convinced it warrants this special implementation independent of reduce(hcat, ...) and alternatives. I could be wrong, but I feel it's a big deal, and core developers did not pay enough attention to this.

I don't mean to criticize the contribution. I'm very interested in it. Maybe I'm wrong and eventually I'll change my mind, who knows. But right now it's not my feeling at all. This is an important feature that the language lacks. It's not such a big deal, reason why there's only now a PR for it. But I don't feel we have achieved the best design. Again, I don't mean to dunk on your work and I'm happy to see progress being done, but I'm not personally satisfied with the direction this is going. Hope I'm wrong.

[nlw0]

This may not be the end of the story. stack is no help for vcat(vec_of_vecs...) which is a common desire, and I think there's some agreement that the reduce(vcat, vv) pattern isn't great.

One possibility is something like flat(x) == collect(Iterators.flatten(x)) (which always makes a vector). Another would be something like the function you call blockstack (which concatenates inner dims along outer dims -- perhaps seen as generalising hvcat). Making a prototype somewhere (and figuring out how to make it fast, what to do about OffsetArrays, how to clearly describe precisely what it does) isn't a bad idea.

I'm sorry, but doesn't stack(vec_of_vecs, dims=1) implement reduce(vcat, vec_of_vecs)? If it doesn't, why not? And if it does, how come we go on merging without clarity about this? That's precisely my concern. I think this is a very important new feature being introduced to Base, with very little thought about the design. Nevermind the final result, I find the lack of curiosity very disappointing.

[mcabbott]

doesn't stack(vec_of_vecs, dims=1) implement reduce(vcat, vec_of_vecs)?

It does not. stack(vecs; dims=1) is one of the early examples in the docstring. It also matches numpy.stack referenced above.

without clarity about this?

But there is clarity about this. Much thought has gone into the design.

lack of curiosity

And you didn't try it out, by typing using Compat? Read the numpy docs about their variants?

I don't think this is the right level to discuss on this thread, and encourage you to show up on zulip or slack or something less formal.

[nlw0]

OK, well, thanks for pointing out the example to me. I'm sorry I'll probably keep struggling to intuitively understand what the method does.