From 4ea5e403879d1ed22506125bfbcef45c7f43a803 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Fri, 23 Feb 2024 17:27:07 +1100
Subject: [PATCH 01/27] Extend Base.stack to DimArrays

---
 src/array/methods.jl | 50 ++++++++++++++++++++++++++++++++++++++------
 src/tables.jl        |  2 +-
 2 files changed, 45 insertions(+), 7 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index ad9d16ea7..e2bd271cc 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -131,7 +131,7 @@ end
     """
     function Base.eachslice(A::AbstractDimArray; dims)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ()) 
+        if !(dimtuple == ())
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple)
@@ -139,7 +139,7 @@ end
 else
     @inline function Base.eachslice(A::AbstractDimArray; dims, drop=true)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ()) 
+        if !(dimtuple == ())
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple, drop)
@@ -195,7 +195,7 @@ function Base.sortslices(A::AbstractDimArray; dims, kw...)
     return rebuild(A, newdata, newdims)
 end
 
-                
+
 Base.cumsum(A::AbstractDimVector) = rebuild(A, Base.cumsum(parent(A)))
 Base.cumsum(A::AbstractDimArray; dims) = rebuild(A, cumsum(parent(A); dims=dimnum(A, dims)))
 Base.cumsum!(B::AbstractArray, A::AbstractDimVector) = cumsum!(B, parent(A))
@@ -293,7 +293,7 @@ function _cat(catdims::Tuple, A1::AbstractDimArray, As::AbstractDimArray...)
             else
                 # vcat the index for the catdim in each of Xin
                 joindims = map(A -> dims(A, catdim), Xin)
-                if !check_cat_lookups(joindims...) 
+                if !check_cat_lookups(joindims...)
                     return rebuild(catdim, NoLookup())
                 end
                 _vcat_dims(joindims...)
@@ -367,7 +367,7 @@ function Base.vcat(As::Union{AbstractDimVector,AbstractDimMatrix}...)
         (catdim,)
     else
         # Make sure this is the same dimension for all arrays
-        if !comparedims(Bool, map(x -> dims(x, 2), As)...; 
+        if !comparedims(Bool, map(x -> dims(x, 2), As)...;
             val=true, warn = " Can't `vcat` AbstractDimArray, applying to `parent` object."
         )
             return Base.vcat(map(parent, As)...)
@@ -542,11 +542,49 @@ end
 _span_string(D, S, span) = _cat_warn_string(D, "not all lookups have `$S` spans. Found $(basetypeof(span))")
 _cat_warn_string(D, message) = """
 `cat` cannot concatenate `Dimension`s, falling back to `parent` type:
-$message on dimension $D. 
+$message on dimension $D.
 
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
 
+function Base._typed_stack(::Colon, ::Type{T}, ::Type{S}, A, Aax=_iterator_axes(A)) where {T,S<:AbstractDimArray}
+    origdims = dims.(A)
+    _A = parent.(A)
+    t = eltype(_A)
+    _A = Base._typed_stack(:, T, t, A)
+
+    if !comparedims(Bool, origdims...;
+        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
+    )
+        return _A
+    else
+        DimArray(_A, (first(origdims)..., AnonDim()))
+    end
+end
+
+function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:AbstractDimArray}
+    origdims = dims.(A)
+    _A = parent.(A)
+    t = eltype(_A)
+    _A = Base._dim_stack(newdim, T, t, A)
+
+    if !comparedims(Bool, origdims...;
+        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
+    )
+        return _A
+    end
+
+    newdims = first(origdims)
+    newdims = ntuple(d -> d == newdim ? AnonDim() : newdims[d-(d>newdim)], length(newdims) + 1)
+    DimArray(_A, newdims)
+end
+
+function Base.stack(dim::Dimension, A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
+    B = Base.stack(A; dims, kwargs...)
+    newdims = ntuple(d -> d == dims ? dim : DimensionalData.dims(B, d), ndims(B))
+    rebuild(B; dims=newdims)
+end
+
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))
     newdims = reverse(dims(A))
diff --git a/src/tables.jl b/src/tables.jl
index 5470139a2..26c72463f 100644
--- a/src/tables.jl
+++ b/src/tables.jl
@@ -149,7 +149,7 @@ Tables.columnaccess(::Type{<:DimTable}) = true
 Tables.columns(t::DimTable) = t
 Tables.columnnames(c::DimTable) = colnames(c)
 
-function Tables.schema(t::DimTable) 
+function Tables.schema(t::DimTable)
     types = vcat([map(eltype, dimcolumns(t))...], [map(eltype, dimarraycolumns(t))...])
     Tables.Schema(colnames(t), types)
 end

From ec3b2cfbd0e2d12467106f65bf23628e09f9b2a1 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Fri, 23 Feb 2024 17:28:55 +1100
Subject: [PATCH 02/27] Initial Base.stack tests

---
 test/methods.jl | 84 ++++++++++++++++++++++++++++---------------------
 1 file changed, 49 insertions(+), 35 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index 01bb29fd5..b7b499510 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -135,7 +135,7 @@ end
     @test dropdims(da[X(1:1)]; dims=X) == [1, 2, 3]
     @test dropdims(da[2:2, 1:1]; dims=(X(), Y()))[] == 4
     @test typeof(dropdims(da[2:2, 1:1]; dims=(X(), Y()))) <: DimArray{Int,0,Tuple{}}
-    @test refdims(dropdims(da[X(1:1)]; dims=X)) == 
+    @test refdims(dropdims(da[X(1:1)]; dims=X)) ==
         (X(Sampled(143:2:143, ForwardOrdered(), Regular(2), Points(), NoMetadata())),)
     dropped = dropdims(da[X(1:1)]; dims=X)
     @test dropped[1:2] == [1, 2]
@@ -240,7 +240,7 @@ end
     @test tda == transpose(parent(da))
     resultdims = (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
                   Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
-    @test typeof(dims(tda)) == typeof(resultdims) 
+    @test typeof(dims(tda)) == typeof(resultdims)
     @test dims(tda) == resultdims
     @test size(tda) == (4, 5)
 
@@ -319,14 +319,14 @@ end
     for dims in xs
         cvda = cov(da; dims=X)
         @test cvda == cov(a; dims=2)
-        @test DimensionalData.dims(cvda) == 
+        @test DimensionalData.dims(cvda) ==
             (Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())),
              Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
     end
     for dims in ys
         crda = cor(da; dims)
         @test crda == cor(a; dims=1)
-        @test DimensionalData.dims(crda) == 
+        @test DimensionalData.dims(crda) ==
             (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
              X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())))
     end
@@ -336,7 +336,7 @@ end
     a = [1 2 3 4
          3 4 5 6
          5 6 7 8]
-    y = Y(Sampled(10:10:30; sampling=Intervals())) 
+    y = Y(Sampled(10:10:30; sampling=Intervals()))
     ti = Ti(Sampled(1:4; sampling=Intervals()))
     da = DimArray(a, (y, ti))
     ys = (1, Y, Y(), :Y, y)
@@ -344,7 +344,7 @@ end
     for dims in ys
         ms = mapslices(sum, da; dims)
         @test ms == [9 12 15 18]
-        @test DimensionalData.dims(ms) == 
+        @test DimensionalData.dims(ms) ==
             (Y(NoLookup(Base.OneTo(1))), Ti(Sampled(1:4, ForwardOrdered(), Regular(1), Intervals(Start()), NoMetadata())))
         @test refdims(ms) == ()
     end
@@ -390,13 +390,13 @@ end
     @test dims(sort(v)) == (X(NoLookup(Base.OneTo(10))),)
     A = rand((X(5:-1:1), Y(11:15)))
     @test sort(A; dims=X) == sort(parent(A); dims=1)
-    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y)) 
+    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y))
 end
 
 @testset "sortslices" begin
     M = rand((X(5:-1:1), Y(11:15)))
     @test sortslices(M; dims=X) == sortslices(parent(M); dims=1)
-    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y)) 
+    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y))
     M = rand((X(5:-1:1), Y(11:15), Ti(3:10)))
     @test sortslices(M; dims=(X, Y)) == sortslices(parent(M); dims=(1, 2))
     @test dims(sortslices(M; dims=(X, Y))) == (X(NoLookup(Base.OneTo(5))), Y(NoLookup(Base.OneTo(5))), dims(M, Ti))
@@ -419,7 +419,7 @@ end
         @test_throws DimensionMismatch vcat(dims(da, 1), dims(de, 1))
         testdims = (X(Sampled([4.0, 5.0, 6.0, 7.0], ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
                     Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
-        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X()) 
+        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X())
         @test cat(da, db; dims=X) == cat(da, db; dims=(X,)) == cat(da, db; dims=1) == cat(da, db; dims=(1,))
         @test dims(cat(da, db; dims=X)) == testdims
         @test val(cat(da, db; dims=X)) == val(testdims)
@@ -460,7 +460,7 @@ end
 
     @testset "Irregular Sampled" begin
         @testset "Intervals" begin
-            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata())) 
+            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata()))
             d2 = X(Sampled([7, 8], ForwardOrdered(), Irregular(7, 9), Intervals(), NoMetadata()))
             iri_dim = vcat(d1, d2)
             @test span(iri_dim) == Irregular(1, 9)
@@ -488,7 +488,7 @@ end
         d2 = X(Sampled([7.5, 9], ForwardOrdered(), Explicit([7 8; 8 10]), Intervals(Center()), NoMetadata()))
         ed = vcat(d1, d2)
         @test span(ed) == Explicit([1 3 4 7 8; 3 4 7 8 10])
-        @test index(ed) == [2, 3.5, 5, 7.5, 9] 
+        @test index(ed) == [2, 3.5, 5, 7.5, 9]
         @test lookup(ed) == Sampled([2, 3.5, 5, 7.5, 9], ForwardOrdered(), Explicit([1 3 4 7 8; 3 4 7 8 10]), Intervals(Center()), NoMetadata())
         @test_warn "lookups are mixed `ForwardOrdered` and `ReverseOrdered`" vcat(d1, reverse(d2))
         @test_warn "lookups are misaligned" vcat(d2, d1)
@@ -556,16 +556,16 @@ end
 
         @test vcat(da) isa DimArray{Int,1}
         @test vcat(da) == da
-        @test dims(vcat(da)) == 
+        @test dims(vcat(da)) ==
             dims(cat(da; dims=1)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) == 
+        @test dims(vcat(da, db)) ==
             dims(cat(da, db; dims=1)) ==
             (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "X and Y dims on the same axis" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) == 
+        @test dims(vcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=1)) ==
             (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "do not match" hcat(da, db, dd)
@@ -582,21 +582,21 @@ end
         dd = DimArray(c, (X(8.0:9.0), Z(6.0:8.0)))
 
         @test vcat(da) == da
-        @test dims(vcat(da)) == 
+        @test dims(vcat(da)) ==
             dims(cat(da; dims=1)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) == 
+        @test dims(vcat(da, db)) ==
             dims(cat(da, db; dims=1)) ==
-            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "lookups are misaligned" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) == 
+        @test dims(vcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=1)) ==
-            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "lookups are misaligned" hcat(da, db, dd)
     end
 end
@@ -612,16 +612,16 @@ end
         dd = DimArray(c, X(8.0:9.0))
 
         @test hcat(da) == permutedims([1 2])
-        @test dims(hcat(da)) == 
+        @test dims(hcat(da)) ==
             dims(cat(da; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(1))))
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) == 
+        @test dims(hcat(da, db)) ==
             dims(cat(da, db; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(2))))
         @test_warn "do not match" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
-        @test dims(hcat(da, db, dc)) == 
+        @test dims(hcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(3))))
         @test_warn "do not match" hcat(da, db, dd)
@@ -637,24 +637,38 @@ end
 
         @test hcat(da) isa DimArray{Int,2}
         @test hcat(da) == da
-        @test dims(hcat(da)) == 
+        @test dims(hcat(da)) ==
             dims(cat(da; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) == 
+        @test dims(hcat(da, db)) ==
             dims(cat(da, db; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "do not join with the correct step size" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
         @test dims(hcat(da, db, dc)) == dims(cat(da, db, dc; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "do not match" hcat(da, db, dd)
     end
 end
 
+@testset "Base.stack" begin
+    a = [1 2 3; 4 5 6]
+    da = DimArray(a, (X(4.0:5.0), Y(6.0:8.0)))
+    b = [7 8 9; 10 11 12]
+    ca = DimArray(b, (X(4.0:5.0), Y(6.0:8.0)))
+    db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))
+
+    @test stack([da, db]; dims=3) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=3)
+    @test_warn "Lookup values for X" stack([da, db]; dims=3)
+
+    @test stack([da, ca]; dims=1) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=1)
+    @test_warn "Lookup values for X" stack([da, db]; dims=1)
+end
+
 @testset "unique" begin
     a = [1 1 6; 1 1 6]
     xs = (X, X(), :X)
@@ -684,7 +698,7 @@ end
             @test diff(A; dims) == DimArray([111 93 -169 142; 98 -55 110 -131], (Y(['a', 'b']), ti))
         end
         for dims in tis
-            @test diff(A; dims) == 
+            @test diff(A; dims) ==
                 DimArray([38 156 -125; 20 -106 186; -133 59 -55], (y, Ti(DateTime(2021, 1):Month(1):DateTime(2021, 3))))
         end
         @test_throws ArgumentError diff(A; dims='X')

From 0b762ced69d63fe4a5f8c249b2b5d164867662f7 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Sat, 24 Feb 2024 00:01:03 +1100
Subject: [PATCH 03/27] Fix whitespace changes

---
 src/array/methods.jl | 38 ++++++++++++++++++++++++++++++++++++++
 test/methods.jl      | 14 ++++++++++++++
 2 files changed, 52 insertions(+)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index ad9d16ea7..8ebcabbb4 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -547,6 +547,44 @@ $message on dimension $D.
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
 
+function Base._typed_stack(::Colon, ::Type{T}, ::Type{S}, A, Aax=_iterator_axes(A)) where {T,S<:AbstractDimArray}
+    origdims = dims.(A)
+    _A = parent.(A)
+    t = eltype(_A)
+    _A = Base._typed_stack(:, T, t, A)
+
+    if !comparedims(Bool, origdims...;
+        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
+    )
+        return _A
+    else
+        DimArray(_A, (first(origdims)..., AnonDim()))
+    end
+end
+
+function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:AbstractDimArray}
+    origdims = dims.(A)
+    _A = parent.(A)
+    t = eltype(_A)
+    _A = Base._dim_stack(newdim, T, t, A)
+
+    if !comparedims(Bool, origdims...;
+        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
+    )
+        return _A
+    end
+
+    newdims = first(origdims)
+    newdims = ntuple(d -> d == newdim ? AnonDim() : newdims[d-(d>newdim)], length(newdims) + 1)
+    DimArray(_A, newdims)
+end
+
+function Base.stack(dim::Dimension, A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
+    B = Base.stack(A; dims, kwargs...)
+    newdims = ntuple(d -> d == dims ? dim : DimensionalData.dims(B, d), ndims(B))
+    rebuild(B; dims=newdims)
+end
+
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))
     newdims = reverse(dims(A))
diff --git a/test/methods.jl b/test/methods.jl
index 01bb29fd5..d0aa669f3 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -655,6 +655,20 @@ end
     end
 end
 
+@testset "Base.stack" begin
+    a = [1 2 3; 4 5 6]
+    da = DimArray(a, (X(4.0:5.0), Y(6.0:8.0)))
+    b = [7 8 9; 10 11 12]
+    ca = DimArray(b, (X(4.0:5.0), Y(6.0:8.0)))
+    db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))
+
+    @test stack([da, db]; dims=3) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=3)
+    @test_warn "Lookup values for X" stack([da, db]; dims=3)
+
+    @test stack([da, ca]; dims=1) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=1)
+    @test_warn "Lookup values for X" stack([da, db]; dims=1)
+end
+
 @testset "unique" begin
     a = [1 1 6; 1 1 6]
     xs = (X, X(), :X)

From ac3affcbdd7249d414d93e019bf3edff7e807116 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Tue, 27 Feb 2024 15:12:59 +1100
Subject: [PATCH 04/27] Clean up Base.stack, add docstrings, add tests

---
 src/array/methods.jl | 48 +++++++++++++++++++++++++++++++++++++++++---
 test/methods.jl      | 10 +++++++--
 2 files changed, 53 insertions(+), 5 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 8ebcabbb4..ca233f5fe 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -548,7 +548,7 @@ To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(new
 """
 
 function Base._typed_stack(::Colon, ::Type{T}, ::Type{S}, A, Aax=_iterator_axes(A)) where {T,S<:AbstractDimArray}
-    origdims = dims.(A)
+    origdims = map(dims, A)
     _A = parent.(A)
     t = eltype(_A)
     _A = Base._typed_stack(:, T, t, A)
@@ -575,13 +575,55 @@ function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:A
     end
 
     newdims = first(origdims)
-    newdims = ntuple(d -> d == newdim ? AnonDim() : newdims[d-(d>newdim)], length(newdims) + 1)
+    newdims = ntuple(length(newdims) + 1) do d
+        if d == newdim
+            AnonDim()
+        else # Return the old dimension, shifted across once if it comes after the new dim
+            newdims[d-(d>newdim)]
+        end
+    end
     DimArray(_A, newdims)
 end
 
+"""
+    Base.stack([dim::Dimension], A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
+
+Stack arrays along a specified axis `dims`, while preserving the dimensional
+information of other axes.
+If the optional `Dimension` argument `dim` is supplied, it must have
+`length(dim) == length(A)`; the resulting axis `dims` is given the dimension `dim`.
+If `dim` is not supplied, the new dimension will be an `AnonDim()`.
+
+# Examples
+```julia
+julia> da = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(300:-100:100)));
+julia> db = DimArray([6 5 4; 3 2 1], (X(10:10:20), Y(300:-100:100)));
+
+# Stack along a new dimension `Z`
+julia> dc = stack(Z(1:2), [da, db], dims=3)
+╭─────────────────────────╮
+│ 2×3×2 DimArray{Int64,3} │
+├─────────────────────────┴──────────────────────────────── dims ┐
+  ↓ X Sampled{Int64} 10:10:20 ForwardOrdered Regular Points,
+  → Y Sampled{Int64} 300:-100:100 ReverseOrdered Regular Points,
+  ↗ Z 1:2
+└────────────────────────────────────────────────────────────────┘
+
+julia> dims(dc, 3) == Z(1:2)
+true
+julia> parent(dc) == stack(map(parent, [da, db]), dims=3)
+true
+```
+"""
 function Base.stack(dim::Dimension, A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
     B = Base.stack(A; dims, kwargs...)
-    newdims = ntuple(d -> d == dims ? dim : DimensionalData.dims(B, d), ndims(B))
+    newdims = ntuple(ndims(B)) do d
+        if d == dims # Use the new provided dimension
+            dim
+        else
+            DimensionalData.dims(B, d)
+        end
+    end
     rebuild(B; dims=newdims)
 end
 
diff --git a/test/methods.jl b/test/methods.jl
index d0aa669f3..26597a3c7 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -662,11 +662,17 @@ end
     ca = DimArray(b, (X(4.0:5.0), Y(6.0:8.0)))
     db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))
 
-    @test stack([da, db]; dims=3) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=3)
+    @test stack([da, db]; dims=3) == stack([parent(da), parent(db)], dims=3)
     @test_warn "Lookup values for X" stack([da, db]; dims=3)
 
-    @test stack([da, ca]; dims=1) == stack([[1 2 3; 4 5 6], [7 8 9; 10 11 12]], dims=1)
+    @test stack([da, ca]; dims=1) == stack([parent(da), parent(ca)], dims=1)
     @test_warn "Lookup values for X" stack([da, db]; dims=1)
+
+    for d = 1:3
+        dc = stack(Z(1:2), [da, ca], dims=d)
+        @test dims(dc, d) == Z(1:2)
+        @test parent(dc) == stack(map(parent, [da, db]), dims=d)
+    end
 end
 
 @testset "unique" begin

From 83ace6072cd5240875bdfd5e17a31281ce9ceb3e Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Tue, 27 Feb 2024 15:54:19 +1100
Subject: [PATCH 05/27] Fix Base.stack docstring

---
 src/array/methods.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index ca233f5fe..9dc66c881 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -586,16 +586,16 @@ function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:A
 end
 
 """
-    Base.stack([dim::Dimension], A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
+    Base.stack([dim::Dimension], A::AbstractVector{<:AbstractDimArray}; dims=nothing)
 
 Stack arrays along a specified axis `dims`, while preserving the dimensional
 information of other axes.
 If the optional `Dimension` argument `dim` is supplied, it must have
 `length(dim) == length(A)`; the resulting axis `dims` is given the dimension `dim`.
-If `dim` is not supplied, the new dimension will be an `AnonDim()`.
+If `dim` is not supplied, the new dimension will be an `AnonDim`.
 
 # Examples
-```julia
+```julia-repl
 julia> da = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(300:-100:100)));
 julia> db = DimArray([6 5 4; 3 2 1], (X(10:10:20), Y(300:-100:100)));
 

From e6cfb1bad10a2989d0938bcf3eb7ca4d89123935 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Fri, 1 Mar 2024 16:13:26 +1100
Subject: [PATCH 06/27] New Pair syntax for extended Base.stack

---
 src/array/methods.jl | 53 +++++++++++++++++-----------
 test/methods.jl      | 82 +++++++++++++++++++++++++-------------------
 2 files changed, 79 insertions(+), 56 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 9dc66c881..b3ae7db5f 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -131,7 +131,7 @@ end
     """
     function Base.eachslice(A::AbstractDimArray; dims)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ()) 
+        if !(dimtuple == ())
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple)
@@ -139,7 +139,7 @@ end
 else
     @inline function Base.eachslice(A::AbstractDimArray; dims, drop=true)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ()) 
+        if !(dimtuple == ())
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple, drop)
@@ -195,7 +195,7 @@ function Base.sortslices(A::AbstractDimArray; dims, kw...)
     return rebuild(A, newdata, newdims)
 end
 
-                
+
 Base.cumsum(A::AbstractDimVector) = rebuild(A, Base.cumsum(parent(A)))
 Base.cumsum(A::AbstractDimArray; dims) = rebuild(A, cumsum(parent(A); dims=dimnum(A, dims)))
 Base.cumsum!(B::AbstractArray, A::AbstractDimVector) = cumsum!(B, parent(A))
@@ -293,7 +293,7 @@ function _cat(catdims::Tuple, A1::AbstractDimArray, As::AbstractDimArray...)
             else
                 # vcat the index for the catdim in each of Xin
                 joindims = map(A -> dims(A, catdim), Xin)
-                if !check_cat_lookups(joindims...) 
+                if !check_cat_lookups(joindims...)
                     return rebuild(catdim, NoLookup())
                 end
                 _vcat_dims(joindims...)
@@ -367,7 +367,7 @@ function Base.vcat(As::Union{AbstractDimVector,AbstractDimMatrix}...)
         (catdim,)
     else
         # Make sure this is the same dimension for all arrays
-        if !comparedims(Bool, map(x -> dims(x, 2), As)...; 
+        if !comparedims(Bool, map(x -> dims(x, 2), As)...;
             val=true, warn = " Can't `vcat` AbstractDimArray, applying to `parent` object."
         )
             return Base.vcat(map(parent, As)...)
@@ -542,7 +542,7 @@ end
 _span_string(D, S, span) = _cat_warn_string(D, "not all lookups have `$S` spans. Found $(basetypeof(span))")
 _cat_warn_string(D, message) = """
 `cat` cannot concatenate `Dimension`s, falling back to `parent` type:
-$message on dimension $D. 
+$message on dimension $D.
 
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
@@ -586,13 +586,17 @@ function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:A
 end
 
 """
-    Base.stack([dim::Dimension], A::AbstractVector{<:AbstractDimArray}; dims=nothing)
+    Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1]), AnonDim()))
+
+Stack arrays along a new axis while preserving the dimensional information of other axes.
+
+The optional keyword argument `dims` has the following behavior:
+- `dims isa Integer`: The dimension of the new axis is an `AnonDim` at position `dims`
+- `dims isa Dimension`: The new axis is at `ndims(A[1])+1` and has a dimension of `dims`.
+- `dims isa Pair{Integer, Dimension}`: The new axis is at `first(dims)` and has a dimension
+  of `last(dims)`.
 
-Stack arrays along a specified axis `dims`, while preserving the dimensional
-information of other axes.
-If the optional `Dimension` argument `dim` is supplied, it must have
-`length(dim) == length(A)`; the resulting axis `dims` is given the dimension `dim`.
-If `dim` is not supplied, the new dimension will be an `AnonDim`.
+If `dims` contains a `Dimension`, that `Dimension` must have the same length as A.
 
 # Examples
 ```julia-repl
@@ -615,16 +619,25 @@ julia> parent(dc) == stack(map(parent, [da, db]), dims=3)
 true
 ```
 """
-function Base.stack(dim::Dimension, A::AbstractVector{<:AbstractDimArray}; dims=nothing, kwargs...)
-    B = Base.stack(A; dims, kwargs...)
-    newdims = ntuple(ndims(B)) do d
-        if d == dims # Use the new provided dimension
-            dim
-        else
-            DimensionalData.dims(B, d)
+function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1]), AnonDim()), kwargs...)
+    dims isa Integer && (dims = dims => AnonDim())
+    dims isa Dimension && (dims = ndims(A[1])+1 => dims)
+
+    B = Base._stack(first(dims), A)
+
+    if B isa AbstractDimArray
+        newdims = ntuple(ndims(B)) do d
+            if d == first(dims) # Use the new provided dimension
+                last(dims)
+            else
+                DimensionalData.dims(B, d)
+            end
         end
+        newdims = format(newdims, B)
+
+        B = rebuild(B; dims=newdims)
     end
-    rebuild(B; dims=newdims)
+    return B
 end
 
 function Base.inv(A::AbstractDimArray{T,2}) where T
diff --git a/test/methods.jl b/test/methods.jl
index 26597a3c7..51c993cf6 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -135,7 +135,7 @@ end
     @test dropdims(da[X(1:1)]; dims=X) == [1, 2, 3]
     @test dropdims(da[2:2, 1:1]; dims=(X(), Y()))[] == 4
     @test typeof(dropdims(da[2:2, 1:1]; dims=(X(), Y()))) <: DimArray{Int,0,Tuple{}}
-    @test refdims(dropdims(da[X(1:1)]; dims=X)) == 
+    @test refdims(dropdims(da[X(1:1)]; dims=X)) ==
         (X(Sampled(143:2:143, ForwardOrdered(), Regular(2), Points(), NoMetadata())),)
     dropped = dropdims(da[X(1:1)]; dims=X)
     @test dropped[1:2] == [1, 2]
@@ -240,7 +240,7 @@ end
     @test tda == transpose(parent(da))
     resultdims = (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
                   Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
-    @test typeof(dims(tda)) == typeof(resultdims) 
+    @test typeof(dims(tda)) == typeof(resultdims)
     @test dims(tda) == resultdims
     @test size(tda) == (4, 5)
 
@@ -319,14 +319,14 @@ end
     for dims in xs
         cvda = cov(da; dims=X)
         @test cvda == cov(a; dims=2)
-        @test DimensionalData.dims(cvda) == 
+        @test DimensionalData.dims(cvda) ==
             (Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())),
              Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
     end
     for dims in ys
         crda = cor(da; dims)
         @test crda == cor(a; dims=1)
-        @test DimensionalData.dims(crda) == 
+        @test DimensionalData.dims(crda) ==
             (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
              X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())))
     end
@@ -336,7 +336,7 @@ end
     a = [1 2 3 4
          3 4 5 6
          5 6 7 8]
-    y = Y(Sampled(10:10:30; sampling=Intervals())) 
+    y = Y(Sampled(10:10:30; sampling=Intervals()))
     ti = Ti(Sampled(1:4; sampling=Intervals()))
     da = DimArray(a, (y, ti))
     ys = (1, Y, Y(), :Y, y)
@@ -344,7 +344,7 @@ end
     for dims in ys
         ms = mapslices(sum, da; dims)
         @test ms == [9 12 15 18]
-        @test DimensionalData.dims(ms) == 
+        @test DimensionalData.dims(ms) ==
             (Y(NoLookup(Base.OneTo(1))), Ti(Sampled(1:4, ForwardOrdered(), Regular(1), Intervals(Start()), NoMetadata())))
         @test refdims(ms) == ()
     end
@@ -390,13 +390,13 @@ end
     @test dims(sort(v)) == (X(NoLookup(Base.OneTo(10))),)
     A = rand((X(5:-1:1), Y(11:15)))
     @test sort(A; dims=X) == sort(parent(A); dims=1)
-    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y)) 
+    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y))
 end
 
 @testset "sortslices" begin
     M = rand((X(5:-1:1), Y(11:15)))
     @test sortslices(M; dims=X) == sortslices(parent(M); dims=1)
-    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y)) 
+    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y))
     M = rand((X(5:-1:1), Y(11:15), Ti(3:10)))
     @test sortslices(M; dims=(X, Y)) == sortslices(parent(M); dims=(1, 2))
     @test dims(sortslices(M; dims=(X, Y))) == (X(NoLookup(Base.OneTo(5))), Y(NoLookup(Base.OneTo(5))), dims(M, Ti))
@@ -419,7 +419,7 @@ end
         @test_throws DimensionMismatch vcat(dims(da, 1), dims(de, 1))
         testdims = (X(Sampled([4.0, 5.0, 6.0, 7.0], ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
                     Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
-        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X()) 
+        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X())
         @test cat(da, db; dims=X) == cat(da, db; dims=(X,)) == cat(da, db; dims=1) == cat(da, db; dims=(1,))
         @test dims(cat(da, db; dims=X)) == testdims
         @test val(cat(da, db; dims=X)) == val(testdims)
@@ -460,7 +460,7 @@ end
 
     @testset "Irregular Sampled" begin
         @testset "Intervals" begin
-            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata())) 
+            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata()))
             d2 = X(Sampled([7, 8], ForwardOrdered(), Irregular(7, 9), Intervals(), NoMetadata()))
             iri_dim = vcat(d1, d2)
             @test span(iri_dim) == Irregular(1, 9)
@@ -488,7 +488,7 @@ end
         d2 = X(Sampled([7.5, 9], ForwardOrdered(), Explicit([7 8; 8 10]), Intervals(Center()), NoMetadata()))
         ed = vcat(d1, d2)
         @test span(ed) == Explicit([1 3 4 7 8; 3 4 7 8 10])
-        @test index(ed) == [2, 3.5, 5, 7.5, 9] 
+        @test index(ed) == [2, 3.5, 5, 7.5, 9]
         @test lookup(ed) == Sampled([2, 3.5, 5, 7.5, 9], ForwardOrdered(), Explicit([1 3 4 7 8; 3 4 7 8 10]), Intervals(Center()), NoMetadata())
         @test_warn "lookups are mixed `ForwardOrdered` and `ReverseOrdered`" vcat(d1, reverse(d2))
         @test_warn "lookups are misaligned" vcat(d2, d1)
@@ -556,16 +556,16 @@ end
 
         @test vcat(da) isa DimArray{Int,1}
         @test vcat(da) == da
-        @test dims(vcat(da)) == 
+        @test dims(vcat(da)) ==
             dims(cat(da; dims=1)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) == 
+        @test dims(vcat(da, db)) ==
             dims(cat(da, db; dims=1)) ==
             (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "X and Y dims on the same axis" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) == 
+        @test dims(vcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=1)) ==
             (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "do not match" hcat(da, db, dd)
@@ -582,21 +582,21 @@ end
         dd = DimArray(c, (X(8.0:9.0), Z(6.0:8.0)))
 
         @test vcat(da) == da
-        @test dims(vcat(da)) == 
+        @test dims(vcat(da)) ==
             dims(cat(da; dims=1)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) == 
+        @test dims(vcat(da, db)) ==
             dims(cat(da, db; dims=1)) ==
-            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "lookups are misaligned" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) == 
+        @test dims(vcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=1)) ==
-            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "lookups are misaligned" hcat(da, db, dd)
     end
 end
@@ -612,16 +612,16 @@ end
         dd = DimArray(c, X(8.0:9.0))
 
         @test hcat(da) == permutedims([1 2])
-        @test dims(hcat(da)) == 
+        @test dims(hcat(da)) ==
             dims(cat(da; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(1))))
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) == 
+        @test dims(hcat(da, db)) ==
             dims(cat(da, db; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(2))))
         @test_warn "do not match" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
-        @test dims(hcat(da, db, dc)) == 
+        @test dims(hcat(da, db, dc)) ==
             dims(cat(da, db, dc; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(3))))
         @test_warn "do not match" hcat(da, db, dd)
@@ -637,20 +637,20 @@ end
 
         @test hcat(da) isa DimArray{Int,2}
         @test hcat(da) == da
-        @test dims(hcat(da)) == 
+        @test dims(hcat(da)) ==
             dims(cat(da; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) == 
+        @test dims(hcat(da, db)) ==
             dims(cat(da, db; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "do not join with the correct step size" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
         @test dims(hcat(da, db, dc)) == dims(cat(da, db, dc; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
-             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
+             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
         @test_warn "do not match" hcat(da, db, dd)
     end
 end
@@ -668,8 +668,18 @@ end
     @test stack([da, ca]; dims=1) == stack([parent(da), parent(ca)], dims=1)
     @test_warn "Lookup values for X" stack([da, db]; dims=1)
 
+    dc = stack([da, ca], dims=Z(1:2))
+    @test dims(dc, ndims(da)+1) == Z(1:2)
+    @test parent(dc) == stack(map(parent, [da, db]))
+
+    for d = 1:3
+        dc = stack([da, ca], dims=d)
+        @test dims(dc, d) isa AnonDim
+        @test parent(dc) == stack(map(parent, [da, db]), dims=d)
+    end
+
     for d = 1:3
-        dc = stack(Z(1:2), [da, ca], dims=d)
+        dc = stack([da, ca], dims=d=>Z(1:2))
         @test dims(dc, d) == Z(1:2)
         @test parent(dc) == stack(map(parent, [da, db]), dims=d)
     end
@@ -704,7 +714,7 @@ end
             @test diff(A; dims) == DimArray([111 93 -169 142; 98 -55 110 -131], (Y(['a', 'b']), ti))
         end
         for dims in tis
-            @test diff(A; dims) == 
+            @test diff(A; dims) ==
                 DimArray([38 156 -125; 20 -106 186; -133 59 -55], (y, Ti(DateTime(2021, 1):Month(1):DateTime(2021, 3))))
         end
         @test_throws ArgumentError diff(A; dims='X')

From ce16068889487bd96ca7b2fc8bf5029b609db9b0 Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Thu, 7 Mar 2024 14:31:08 +1100
Subject: [PATCH 07/27] Better Base.stack formatting

---
 src/array/methods.jl | 17 +++++++++--------
 test/methods.jl      |  1 +
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index b3ae7db5f..9d0cc82d6 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -586,7 +586,7 @@ function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:A
 end
 
 """
-    Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1]), AnonDim()))
+    Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()))
 
 Stack arrays along a new axis while preserving the dimensional information of other axes.
 
@@ -604,7 +604,7 @@ julia> da = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(300:-100:100)));
 julia> db = DimArray([6 5 4; 3 2 1], (X(10:10:20), Y(300:-100:100)));
 
 # Stack along a new dimension `Z`
-julia> dc = stack(Z(1:2), [da, db], dims=3)
+julia> dc = stack([da, db], dims=3=>Z(1:2))
 ╭─────────────────────────╮
 │ 2×3×2 DimArray{Int64,3} │
 ├─────────────────────────┴──────────────────────────────── dims ┐
@@ -619,9 +619,12 @@ julia> parent(dc) == stack(map(parent, [da, db]), dims=3)
 true
 ```
 """
-function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1]), AnonDim()), kwargs...)
-    dims isa Integer && (dims = dims => AnonDim())
-    dims isa Dimension && (dims = ndims(A[1])+1 => dims)
+function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()), kwargs...)
+    if dims isa Integer
+        dims = dims => AnonDim()
+    elseif dims isa Dimension
+        dims = ndims(A[1])+1 => dims
+    end
 
     B = Base._stack(first(dims), A)
 
@@ -633,9 +636,7 @@ function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])
                 DimensionalData.dims(B, d)
             end
         end
-        newdims = format(newdims, B)
-
-        B = rebuild(B; dims=newdims)
+        B = rebuild(B; dims=format(newdims, B))
     end
     return B
 end
diff --git a/test/methods.jl b/test/methods.jl
index 51c993cf6..f70f650bb 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -663,6 +663,7 @@ end
     db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))
 
     @test stack([da, db]; dims=3) == stack([parent(da), parent(db)], dims=3)
+    @test stack([da, db]; dims=3) == stack([da, db]) # Test default dims
     @test_warn "Lookup values for X" stack([da, db]; dims=3)
 
     @test stack([da, ca]; dims=1) == stack([parent(da), parent(ca)], dims=1)

From bdf151f27215bd5bee2465832655f500424eea8f Mon Sep 17 00:00:00 2001
From: brendanjohnharris <bhar9988@uni.sydney.edu.au>
Date: Thu, 7 Mar 2024 15:31:07 +1100
Subject: [PATCH 08/27] Tweak docstring for Base.stack, remove kwargs

---
 src/array/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 9d0cc82d6..84c51f5d7 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -619,7 +619,7 @@ julia> parent(dc) == stack(map(parent, [da, db]), dims=3)
 true
 ```
 """
-function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()), kwargs...)
+function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()))
     if dims isa Integer
         dims = dims => AnonDim()
     elseif dims isa Dimension

From 331485210cf413a40c01cbc0548f90237b0274d7 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Thu, 7 Mar 2024 16:10:40 +1100
Subject: [PATCH 09/27] Update whitespace

---
 src/array/methods.jl | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 23ec67cbb..50af61962 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -109,7 +109,7 @@ end
     ds = DD.dims(A, _astuple(dims))
     # Run one slice with dimensions to get the transformed dim
     d_inds = map(d -> rebuild(d, 1), otherdims(A, ds))
-    example_dims = length(d_inds) > 0 ? DD.dims(f(view(A, d_inds...))) : ()
+    example_dims = DD.dims(f(view(A, d_inds...)))
     replacement_dims = if isnothing(example_dims) || length(example_dims) != length(ds)
         map(d -> rebuild(d, NoLookup()), ds)
     else
@@ -131,7 +131,7 @@ end
     """
     function Base.eachslice(A::AbstractDimArray; dims)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ())
+        if !(dimtuple == ()) 
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple)
@@ -139,7 +139,7 @@ end
 else
     @inline function Base.eachslice(A::AbstractDimArray; dims, drop=true)
         dimtuple = _astuple(dims)
-        if !(dimtuple == ())
+        if !(dimtuple == ()) 
             all(hasdim(A, dimtuple...)) || throw(DimensionMismatch("A doesn't have all dimensions $dims"))
         end
         _eachslice(A, dimtuple, drop)
@@ -195,7 +195,7 @@ function Base.sortslices(A::AbstractDimArray; dims, kw...)
     return rebuild(A, newdata, newdims)
 end
 
-
+                
 Base.cumsum(A::AbstractDimVector) = rebuild(A, Base.cumsum(parent(A)))
 Base.cumsum(A::AbstractDimArray; dims) = rebuild(A, cumsum(parent(A); dims=dimnum(A, dims)))
 Base.cumsum!(B::AbstractArray, A::AbstractDimVector) = cumsum!(B, parent(A))
@@ -293,7 +293,7 @@ function _cat(catdims::Tuple, A1::AbstractDimArray, As::AbstractDimArray...)
             else
                 # vcat the index for the catdim in each of Xin
                 joindims = map(A -> dims(A, catdim), Xin)
-                if !check_cat_lookups(joindims...)
+                if !check_cat_lookups(joindims...) 
                     return rebuild(catdim, NoLookup())
                 end
                 _vcat_dims(joindims...)
@@ -367,7 +367,7 @@ function Base.vcat(As::Union{AbstractDimVector,AbstractDimMatrix}...)
         (catdim,)
     else
         # Make sure this is the same dimension for all arrays
-        if !comparedims(Bool, map(x -> dims(x, 2), As)...;
+        if !comparedims(Bool, map(x -> dims(x, 2), As)...; 
             val=true, warn = " Can't `vcat` AbstractDimArray, applying to `parent` object."
         )
             return Base.vcat(map(parent, As)...)
@@ -390,8 +390,8 @@ end
 check_cat_lookups(dims::Dimension...) =
     _check_cat_lookups(basetypeof(first(dims)), lookup(dims)...)
 
-# Lookups may need adjustment for `cat`
-_check_cat_lookups(D, lookups::Lookup...) = _check_cat_lookup_order(D, lookups...)
+# LookupArrays may need adjustment for `cat`
+_check_cat_lookups(D, lookups::LookupArray...) = _check_cat_lookup_order(D, lookups...)
 _check_cat_lookups(D, l1::NoLookup, lookups::NoLookup...) = true
 function _check_cat_lookups(D, l1::AbstractSampled, lookups::AbstractSampled...)
     length(lookups) > 0 || return true
@@ -430,7 +430,7 @@ function _check_cat_lookups(D, ::Irregular, lookups...)
     end |> all
 end
 
-function _check_cat_lookup_order(D, lookups::Lookup...)
+function _check_cat_lookup_order(D, lookups::LookupArray...)
     l1 = first(lookups)
     length(l1) == 0 && return _check_cat_lookup_order(D, Base.tail(lookups)...)
     L = basetypeof(l1)
@@ -481,8 +481,8 @@ function _vcat_dims(d1::Dimension, ds::Dimension...)
     return rebuild(d1, newlookup)
 end
 
-# Lookups may need adjustment for `cat`
-function _vcat_lookups(lookups::Lookup...)
+# LookupArrays may need adjustment for `cat`
+function _vcat_lookups(lookups::LookupArray...)
     newindex = _vcat_index(lookups...)
     return rebuild(lookups[1]; data=newindex)
 end
@@ -520,10 +520,10 @@ end
 _vcat_index(A1::NoLookup, A::NoLookup...) = OneTo(mapreduce(length, +, (A1, A...)))
 # TODO: handle vcat OffsetArrays?
 # Otherwise just vcat. TODO: handle order breaking vcat?
-# function _vcat_index(lookup::Lookup, lookups...)
+# function _vcat_index(lookup::LookupArray, lookups...)
     # _vcat_index(span(lookup), lookup, lookups...)
 # end
-function _vcat_index(lookup1::Lookup, lookups::Lookup...)
+function _vcat_index(lookup1::LookupArray, lookups::LookupArray...)
     shifted = map((lookup1, lookups...)) do l
         parent(maybeshiftlocus(locus(lookup1), l))
     end
@@ -542,7 +542,7 @@ end
 _span_string(D, S, span) = _cat_warn_string(D, "not all lookups have `$S` spans. Found $(basetypeof(span))")
 _cat_warn_string(D, message) = """
 `cat` cannot concatenate `Dimension`s, falling back to `parent` type:
-$message on dimension $D.
+$message on dimension $D. 
 
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """

From 0ac2eb057b1644db26c2fa612bae46359395f686 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Thu, 7 Mar 2024 16:16:54 +1100
Subject: [PATCH 10/27] Fix whitespace

---
 src/array/methods.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 50af61962..e5ca00349 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -109,7 +109,7 @@ end
     ds = DD.dims(A, _astuple(dims))
     # Run one slice with dimensions to get the transformed dim
     d_inds = map(d -> rebuild(d, 1), otherdims(A, ds))
-    example_dims = DD.dims(f(view(A, d_inds...)))
+    example_dims = length(d_inds) > 0 ? DD.dims(f(view(A, d_inds...))) : ()
     replacement_dims = if isnothing(example_dims) || length(example_dims) != length(ds)
         map(d -> rebuild(d, NoLookup()), ds)
     else
@@ -390,8 +390,8 @@ end
 check_cat_lookups(dims::Dimension...) =
     _check_cat_lookups(basetypeof(first(dims)), lookup(dims)...)
 
-# LookupArrays may need adjustment for `cat`
-_check_cat_lookups(D, lookups::LookupArray...) = _check_cat_lookup_order(D, lookups...)
+# Lookups may need adjustment for `cat`
+_check_cat_lookups(D, lookups::Lookup...) = _check_cat_lookup_order(D, lookups...)
 _check_cat_lookups(D, l1::NoLookup, lookups::NoLookup...) = true
 function _check_cat_lookups(D, l1::AbstractSampled, lookups::AbstractSampled...)
     length(lookups) > 0 || return true
@@ -430,7 +430,7 @@ function _check_cat_lookups(D, ::Irregular, lookups...)
     end |> all
 end
 
-function _check_cat_lookup_order(D, lookups::LookupArray...)
+function _check_cat_lookup_order(D, lookups::Lookup...)
     l1 = first(lookups)
     length(l1) == 0 && return _check_cat_lookup_order(D, Base.tail(lookups)...)
     L = basetypeof(l1)
@@ -481,8 +481,8 @@ function _vcat_dims(d1::Dimension, ds::Dimension...)
     return rebuild(d1, newlookup)
 end
 
-# LookupArrays may need adjustment for `cat`
-function _vcat_lookups(lookups::LookupArray...)
+# Lookups may need adjustment for `cat`
+function _vcat_lookups(lookups::Lookup...)
     newindex = _vcat_index(lookups...)
     return rebuild(lookups[1]; data=newindex)
 end
@@ -520,10 +520,10 @@ end
 _vcat_index(A1::NoLookup, A::NoLookup...) = OneTo(mapreduce(length, +, (A1, A...)))
 # TODO: handle vcat OffsetArrays?
 # Otherwise just vcat. TODO: handle order breaking vcat?
-# function _vcat_index(lookup::LookupArray, lookups...)
+# function _vcat_index(lookup::Lookup, lookups...)
     # _vcat_index(span(lookup), lookup, lookups...)
 # end
-function _vcat_index(lookup1::LookupArray, lookups::LookupArray...)
+function _vcat_index(lookup1::Lookup, lookups::Lookup...)
     shifted = map((lookup1, lookups...)) do l
         parent(maybeshiftlocus(locus(lookup1), l))
     end

From 884ab687a47e8bc25ccc1198ed008623e511ec11 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Thu, 7 Mar 2024 16:19:13 +1100
Subject: [PATCH 11/27] Fix whitespace

---
 test/methods.jl | 70 ++++++++++++++++++++++++-------------------------
 1 file changed, 35 insertions(+), 35 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index f616d9da6..fcacf2de9 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -135,7 +135,7 @@ end
     @test dropdims(da[X(1:1)]; dims=X) == [1, 2, 3]
     @test dropdims(da[2:2, 1:1]; dims=(X(), Y()))[] == 4
     @test typeof(dropdims(da[2:2, 1:1]; dims=(X(), Y()))) <: DimArray{Int,0,Tuple{}}
-    @test refdims(dropdims(da[X(1:1)]; dims=X)) ==
+    @test refdims(dropdims(da[X(1:1)]; dims=X)) == 
         (X(Sampled(143:2:143, ForwardOrdered(), Regular(2), Points(), NoMetadata())),)
     dropped = dropdims(da[X(1:1)]; dims=X)
     @test dropped[1:2] == [1, 2]
@@ -240,7 +240,7 @@ end
     @test tda == transpose(parent(da))
     resultdims = (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
                   Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
-    @test typeof(dims(tda)) == typeof(resultdims)
+    @test typeof(dims(tda)) == typeof(resultdims) 
     @test dims(tda) == resultdims
     @test size(tda) == (4, 5)
 
@@ -319,14 +319,14 @@ end
     for dims in xs
         cvda = cov(da; dims=X)
         @test cvda == cov(a; dims=2)
-        @test DimensionalData.dims(cvda) ==
+        @test DimensionalData.dims(cvda) == 
             (Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())),
              Y(Sampled(LinRange(10.0, 20.0, 5), ForwardOrdered(), Regular(2.5), Points(), NoMetadata())))
     end
     for dims in ys
         crda = cor(da; dims)
         @test crda == cor(a; dims=1)
-        @test DimensionalData.dims(crda) ==
+        @test DimensionalData.dims(crda) == 
             (X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
              X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())))
     end
@@ -336,7 +336,7 @@ end
     a = [1 2 3 4
          3 4 5 6
          5 6 7 8]
-    y = Y(Sampled(10:10:30; sampling=Intervals()))
+    y = Y(Sampled(10:10:30; sampling=Intervals())) 
     ti = Ti(Sampled(1:4; sampling=Intervals()))
     da = DimArray(a, (y, ti))
     ys = (1, Y, Y(), :Y, y)
@@ -344,7 +344,7 @@ end
     for dims in ys
         ms = mapslices(sum, da; dims)
         @test ms == [9 12 15 18]
-        @test DimensionalData.dims(ms) ==
+        @test DimensionalData.dims(ms) == 
             (Y(NoLookup(Base.OneTo(1))), Ti(Sampled(1:4, ForwardOrdered(), Regular(1), Intervals(Start()), NoMetadata())))
         @test refdims(ms) == ()
     end
@@ -396,13 +396,13 @@ end
     @test dims(sort(v)) == (X(NoLookup(Base.OneTo(10))),)
     A = rand((X(5:-1:1), Y(11:15)))
     @test sort(A; dims=X) == sort(parent(A); dims=1)
-    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y))
+    @test dims(sort(A; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(A, Y)) 
 end
 
 @testset "sortslices" begin
     M = rand((X(5:-1:1), Y(11:15)))
     @test sortslices(M; dims=X) == sortslices(parent(M); dims=1)
-    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y))
+    @test dims(sort(M; dims=X)) == (X(NoLookup(Base.OneTo(5))), dims(M, Y)) 
     M = rand((X(5:-1:1), Y(11:15), Ti(3:10)))
     @test sortslices(M; dims=(X, Y)) == sortslices(parent(M); dims=(1, 2))
     @test dims(sortslices(M; dims=(X, Y))) == (X(NoLookup(Base.OneTo(5))), Y(NoLookup(Base.OneTo(5))), dims(M, Ti))
@@ -425,7 +425,7 @@ end
         @test_throws DimensionMismatch vcat(dims(da, 1), dims(de, 1))
         testdims = (X(Sampled([4.0, 5.0, 6.0, 7.0], ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
                     Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
-        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X())
+        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X()) 
         @test cat(da, db; dims=X) == cat(da, db; dims=(X,)) == cat(da, db; dims=1) == cat(da, db; dims=(1,))
         @test dims(cat(da, db; dims=X)) == testdims
         @test val(cat(da, db; dims=X)) == val(testdims)
@@ -466,7 +466,7 @@ end
 
     @testset "Irregular Sampled" begin
         @testset "Intervals" begin
-            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata()))
+            d1 = X(Sampled([1, 3, 4], ForwardOrdered(), Irregular(1, 5), Intervals(), NoMetadata())) 
             d2 = X(Sampled([7, 8], ForwardOrdered(), Irregular(7, 9), Intervals(), NoMetadata()))
             iri_dim = vcat(d1, d2)
             @test span(iri_dim) == Irregular(1, 9)
@@ -494,7 +494,7 @@ end
         d2 = X(Sampled([7.5, 9], ForwardOrdered(), Explicit([7 8; 8 10]), Intervals(Center()), NoMetadata()))
         ed = vcat(d1, d2)
         @test span(ed) == Explicit([1 3 4 7 8; 3 4 7 8 10])
-        @test index(ed) == [2, 3.5, 5, 7.5, 9]
+        @test index(ed) == [2, 3.5, 5, 7.5, 9] 
         @test lookup(ed) == Sampled([2, 3.5, 5, 7.5, 9], ForwardOrdered(), Explicit([1 3 4 7 8; 3 4 7 8 10]), Intervals(Center()), NoMetadata())
         @test_warn "lookups are mixed `ForwardOrdered` and `ReverseOrdered`" vcat(d1, reverse(d2))
         @test_warn "lookups are misaligned" vcat(d2, d1)
@@ -562,16 +562,16 @@ end
 
         @test vcat(da) isa DimArray{Int,1}
         @test vcat(da) == da
-        @test dims(vcat(da)) ==
+        @test dims(vcat(da)) == 
             dims(cat(da; dims=1)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) ==
+        @test dims(vcat(da, db)) == 
             dims(cat(da, db; dims=1)) ==
             (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "X and Y dims on the same axis" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) ==
+        @test dims(vcat(da, db, dc)) == 
             dims(cat(da, db, dc; dims=1)) ==
             (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),)
         @test_warn "do not match" hcat(da, db, dd)
@@ -588,21 +588,21 @@ end
         dd = DimArray(c, (X(8.0:9.0), Z(6.0:8.0)))
 
         @test vcat(da) == da
-        @test dims(vcat(da)) ==
+        @test dims(vcat(da)) == 
             dims(cat(da; dims=1)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test vcat(da, db) == cat(da, db; dims=1)
-        @test dims(vcat(da, db)) ==
+        @test dims(vcat(da, db)) == 
             dims(cat(da, db; dims=1)) ==
-            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:7.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test_warn "lookups are misaligned" hcat(da, dd)
         @test vcat(da, db, dc) == cat(da, db, dc; dims=1)
-        @test dims(vcat(da, db, dc)) ==
+        @test dims(vcat(da, db, dc)) == 
             dims(cat(da, db, dc; dims=1)) ==
-            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:9.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test_warn "lookups are misaligned" hcat(da, db, dd)
     end
 end
@@ -618,16 +618,16 @@ end
         dd = DimArray(c, X(8.0:9.0))
 
         @test hcat(da) == permutedims([1 2])
-        @test dims(hcat(da)) ==
+        @test dims(hcat(da)) == 
             dims(cat(da; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(1))))
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) ==
+        @test dims(hcat(da, db)) == 
             dims(cat(da, db; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(2))))
         @test_warn "do not match" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
-        @test dims(hcat(da, db, dc)) ==
+        @test dims(hcat(da, db, dc)) == 
             dims(cat(da, db, dc; dims=2)) ==
             (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), AnonDim(NoLookup(Base.OneTo(3))))
         @test_warn "do not match" hcat(da, db, dd)
@@ -643,20 +643,20 @@ end
 
         @test hcat(da) isa DimArray{Int,2}
         @test hcat(da) == da
-        @test dims(hcat(da)) ==
+        @test dims(hcat(da)) == 
             dims(cat(da; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:8.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test hcat(da, db) == cat(da, db; dims=2)
-        @test dims(hcat(da, db)) ==
+        @test dims(hcat(da, db)) == 
             dims(cat(da, db; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:11.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test_warn "do not join with the correct step size" hcat(da, dd)
         @test hcat(da, db, dc) == cat(da, db, dc; dims=2)
         @test dims(hcat(da, db, dc)) == dims(cat(da, db, dc; dims=2)) ==
-            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())),
-             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())))
+            (X(Sampled(4.0:5.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata())), 
+             Y(Sampled(6.0:14.0, ForwardOrdered(), Regular(1.0), Points(), NoMetadata()))) 
         @test_warn "do not match" hcat(da, db, dd)
     end
 end
@@ -721,7 +721,7 @@ end
             @test diff(A; dims) == DimArray([111 93 -169 142; 98 -55 110 -131], (Y(['a', 'b']), ti))
         end
         for dims in tis
-            @test diff(A; dims) ==
+            @test diff(A; dims) == 
                 DimArray([38 156 -125; 20 -106 186; -133 59 -55], (y, Ti(DateTime(2021, 1):Month(1):DateTime(2021, 3))))
         end
         @test_throws ArgumentError diff(A; dims='X')

From c369d068bfd76204fa4e261154086df13299b689 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 19 May 2024 09:58:13 +0200
Subject: [PATCH 12/27] use rebuild

---
 src/array/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index e5ca00349..b3153c3eb 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -558,7 +558,7 @@ function Base._typed_stack(::Colon, ::Type{T}, ::Type{S}, A, Aax=_iterator_axes(
     )
         return _A
     else
-        DimArray(_A, (first(origdims)..., AnonDim()))
+        rebuild(A, _A, format((first(origdims)..., AnonDim()), _A))
     end
 end
 

From 2c7b8b9000631fdad4118c33da6ccadfdbaed652 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 19 May 2024 10:00:18 +0200
Subject: [PATCH 13/27] rebuild

---
 src/array/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index b3153c3eb..dd8d61df5 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -582,7 +582,7 @@ function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:A
             newdims[d-(d>newdim)]
         end
     end
-    DimArray(_A, newdims)
+    rebuild(A, _A, format(newdims, _A))
 end
 
 """

From 0a278a6d3273430b78776eae773295b321f2ff7f Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Fri, 15 Nov 2024 17:30:31 +1100
Subject: [PATCH 14/27] Update Base.stack tests

---
 test/methods.jl | 30 ++++++++++++++----------------
 1 file changed, 14 insertions(+), 16 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index 99f2829b9..0ca7e5201 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -809,27 +809,25 @@ end
     ca = DimArray(b, (X(4.0:5.0), Y(6.0:8.0)))
     db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))
 
-    @test stack([da, db]; dims=3) == stack([parent(da), parent(db)], dims=3)
-    @test stack([da, db]; dims=3) == stack([da, db]) # Test default dims
-    @test_warn "Lookup values for X" stack([da, db]; dims=3)
+    x = DimArray([da, db], (Z(4.0:5.0)))
+    @test_warn "Lookup values" stack(x)
 
-    @test stack([da, ca]; dims=1) == stack([parent(da), parent(ca)], dims=1)
-    @test_warn "Lookup values for X" stack([da, db]; dims=1)
+    x = DimArray(fill(da, 2, 2), (Dim{:a}(4.0:5.0), Dim{:b}(6.0:7.0)))
+    sx = stack(x)
+    @test sx isa AbstractDimArray
+    @test dims(sx) == (dims(first(x))..., dims(x)...)
 
-    dc = stack([da, ca], dims=Z(1:2))
-    @test dims(dc, ndims(da)+1) == Z(1:2)
-    @test parent(dc) == stack(map(parent, [da, db]))
+    x = DimArray([da, ca], (Dim{:a}(1:2),))
+    sx = stack(x; dims=1)
+    @test sx == stack([parent(da), parent(ca)], dims=1)
+    @test sx isa AbstractDimArray
+    @test dims(sx) == (dims(x)..., dims(first(x))...)
 
     for d = 1:3
-        dc = stack([da, ca], dims=d)
+        x = DimArray([da, ca], (AnonDim(),))
+        dc = stack(x, dims=d)
         @test dims(dc, d) isa AnonDim
-        @test parent(dc) == stack(map(parent, [da, db]), dims=d)
-    end
-
-    for d = 1:3
-        dc = stack([da, ca], dims=d=>Z(1:2))
-        @test dims(dc, d) == Z(1:2)
-        @test parent(dc) == stack(map(parent, [da, db]), dims=d)
+        @test parent(dc) == stack([da, db], dims=d)
     end
 end
 

From f41f424bee247a3d272ad649ef7fa3d25bf73ec7 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Fri, 15 Nov 2024 17:33:06 +1100
Subject: [PATCH 15/27] Rework Base.stack

---
 src/array/methods.jl | 98 +++-----------------------------------------
 1 file changed, 6 insertions(+), 92 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 926008983..5ff3276ac 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -529,99 +529,13 @@ $message on dimension $D.
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
 
-function Base._typed_stack(::Colon, ::Type{T}, ::Type{S}, A, Aax=_iterator_axes(A)) where {T,S<:AbstractDimArray}
-    origdims = map(dims, A)
-    _A = parent.(A)
-    t = eltype(_A)
-    _A = Base._typed_stack(:, T, t, A)
-
-    if !comparedims(Bool, origdims...;
-        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
-    )
-        return _A
-    else
-        rebuild(A, _A, format((first(origdims)..., AnonDim()), _A))
-    end
-end
-
-function Base._dim_stack(newdim::Integer, ::Type{T}, ::Type{S}, A) where {T,S<:AbstractDimArray}
-    origdims = dims.(A)
-    _A = parent.(A)
-    t = eltype(_A)
-    _A = Base._dim_stack(newdim, T, t, A)
-
-    if !comparedims(Bool, origdims...;
-        order=true, val=true, warn=" Can't `stack` AbstractDimArray, applying to `parent` object."
-    )
-        return _A
-    end
-
-    newdims = first(origdims)
-    newdims = ntuple(length(newdims) + 1) do d
-        if d == newdim
-            AnonDim()
-        else # Return the old dimension, shifted across once if it comes after the new dim
-            newdims[d-(d>newdim)]
-        end
-    end
-    rebuild(A, _A, format(newdims, _A))
-end
-
-"""
-    Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()))
-
-Stack arrays along a new axis while preserving the dimensional information of other axes.
-
-The optional keyword argument `dims` has the following behavior:
-- `dims isa Integer`: The dimension of the new axis is an `AnonDim` at position `dims`
-- `dims isa Dimension`: The new axis is at `ndims(A[1])+1` and has a dimension of `dims`.
-- `dims isa Pair{Integer, Dimension}`: The new axis is at `first(dims)` and has a dimension
-  of `last(dims)`.
-
-If `dims` contains a `Dimension`, that `Dimension` must have the same length as A.
-
-# Examples
-```julia-repl
-julia> da = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(300:-100:100)));
-julia> db = DimArray([6 5 4; 3 2 1], (X(10:10:20), Y(300:-100:100)));
-
-# Stack along a new dimension `Z`
-julia> dc = stack([da, db], dims=3=>Z(1:2))
-╭─────────────────────────╮
-│ 2×3×2 DimArray{Int64,3} │
-├─────────────────────────┴──────────────────────────────── dims ┐
-  ↓ X Sampled{Int64} 10:10:20 ForwardOrdered Regular Points,
-  → Y Sampled{Int64} 300:-100:100 ReverseOrdered Regular Points,
-  ↗ Z 1:2
-└────────────────────────────────────────────────────────────────┘
-
-julia> dims(dc, 3) == Z(1:2)
-true
-julia> parent(dc) == stack(map(parent, [da, db]), dims=3)
-true
-```
-"""
-function Base.stack(A::AbstractVector{<:AbstractDimArray}; dims=Pair(ndims(A[1])+1, AnonDim()))
-    if dims isa Integer
-        dims = dims => AnonDim()
-    elseif dims isa Dimension
-        dims = ndims(A[1])+1 => dims
-    end
-
-    B = Base._stack(first(dims), A)
-
-    if B isa AbstractDimArray
-        newdims = ntuple(ndims(B)) do d
-            if d == first(dims) # Use the new provided dimension
-                last(dims)
-            else
-                DimensionalData.dims(B, d)
-            end
-        end
-        B = rebuild(B; dims=format(newdims, B))
-    end
-    return B
+function check_stack_dims(iter)
+    comparedims(Bool, dims.(iter)...; order=true, val=true, msg=Dimensions.Warn(" Can't `stack` AbstractDimArray, applying to `parent` object."))
+    iter
 end
+Base.stack(iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, check_stack_dims(iter))
+Base.stack(f, iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, f(x) for x in check_stack_dims(iter))
+Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; dims=:) = _stack(dims, f(xy...) |> check_stack_dims for xy in zip(xs, yzs...))
 
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))

From d9ee5e7064c97a02cea282d3c0249b7b13786c5a Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 12:11:19 +1100
Subject: [PATCH 16/27] Update Base.stack to accept dims<:Dimension

---
 src/array/methods.jl | 29 ++++++++++++++++++++++++-----
 1 file changed, 24 insertions(+), 5 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 5ff3276ac..fc680ef43 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -529,13 +529,32 @@ $message on dimension $D.
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
 
-function check_stack_dims(iter)
-    comparedims(Bool, dims.(iter)...; order=true, val=true, msg=Dimensions.Warn(" Can't `stack` AbstractDimArray, applying to `parent` object."))
-    iter
-end
+"""
+    stack(x::AbstractDimArray; [dims])
+
+Combine a nested `AbstractDimArray` into a single array by arranging the inner arrays
+along one or more new dimensions taken from the outer array.
+The inner `AbstractDimArray`s must all have the same dimensions.
+
+The keyword `dims` specifies where the new dimensions will be placed in the resulting array.
+The index of all existing dimensions equal to or greater than `dims` will be increased to
+accomodate the new dimensions.
+The default `dims` places the new dimensions at the end of the resulting array.
+
+# Examples
+```julia
+using DimensionalData
+a = DimArray([1 2 3; 4 5 6], (X(4.0:5.0), Y(6.0:8.0)))
+b = DimArray([7 8 9; 10 11 12], (X(4.0:5.0), Y(6.0:8.0)))
+x = DimArray([a, b], Z(4.0:5.0)) # Construct a nested DimArray
+y = stack(x; dims=2) # Resulting array has dims (X, *Z*, Y)
+```
+"""
 Base.stack(iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, check_stack_dims(iter))
 Base.stack(f, iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, f(x) for x in check_stack_dims(iter))
-Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; dims=:) = _stack(dims, f(xy...) |> check_stack_dims for xy in zip(xs, yzs...))
+Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; dims=:) = Base._stack(dims, check_stack_dims(f(xy...) for xy in zip(xs, yzs...)))
+Base._stack(dims::Dimension, iter) = Base._stack(typeof(dims), iter)
+Base._stack(dims::Type{<:Dimension}, iter) = Base._stack(dimnum(first(iter), dims), iter)
 
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))

From 9b3b3b7dd46b6c7f2cca88bb362aa8b2e3465b28 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 12:12:27 +1100
Subject: [PATCH 17/27] Improve Base.stack tests

---
 test/methods.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/test/methods.jl b/test/methods.jl
index 0ca7e5201..8a7839365 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -819,6 +819,9 @@ end
 
     x = DimArray([da, ca], (Dim{:a}(1:2),))
     sx = stack(x; dims=1)
+    sy = @test_nowarn stack(x; dims=:)
+    sz = @test_nowarn stack(x; dims=X)
+    @test sx == sz
     @test sx == stack([parent(da), parent(ca)], dims=1)
     @test sx isa AbstractDimArray
     @test dims(sx) == (dims(x)..., dims(first(x))...)
@@ -828,6 +831,9 @@ end
         dc = stack(x, dims=d)
         @test dims(dc, d) isa AnonDim
         @test parent(dc) == stack([da, db], dims=d)
+
+        @test stack(x -> x .^ 2, x) == stack(parent(x)) .^ 2
+        @test stack(+, x, x, x) == stack(x + x + x)
     end
 end
 

From 25dc8a25d98a3f8852eef852cd150371f215262d Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 15:46:40 +1100
Subject: [PATCH 18/27] Add test for Base.stack with symbol dims

---
 test/methods.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/methods.jl b/test/methods.jl
index 8a7839365..8ffdac077 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -821,7 +821,9 @@ end
     sx = stack(x; dims=1)
     sy = @test_nowarn stack(x; dims=:)
     sz = @test_nowarn stack(x; dims=X)
+    sw = @test_nowarn stack(X, dims=:X)
     @test sx == sz
+    @test sz == sw
     @test sx == stack([parent(da), parent(ca)], dims=1)
     @test sx isa AbstractDimArray
     @test dims(sx) == (dims(x)..., dims(first(x))...)

From e2aeb0e6c40d2c30c336f5af7a55a4062a060aec Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 15:53:45 +1100
Subject: [PATCH 19/27] Replace accidentally deleted line for Base.stack

---
 src/array/methods.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index fc680ef43..65e6dd4a8 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -529,6 +529,11 @@ $message on dimension $D.
 To fix for `AbstractDimArray`, pass new lookup values as `cat(As...; dims=$D(newlookupvals))` keyword or `dims=$D()` for empty `NoLookup`.
 """
 
+function check_stack_dims(iter)
+    comparedims(Bool, dims.(iter)...; order=true, val=true, msg=Dimensions.Warn(" Can't `stack` AbstractDimArray, applying to `parent` object."))
+    iter
+end
+
 """
     stack(x::AbstractDimArray; [dims])
 

From 39c78e884ee32ffb6c18c2cc22eed90f58dd3109 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 16:10:34 +1100
Subject: [PATCH 20/27] Update Base.stack for symbol dims

---
 src/array/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 65e6dd4a8..2b56a5f65 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -559,7 +559,7 @@ Base.stack(iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims,
 Base.stack(f, iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, f(x) for x in check_stack_dims(iter))
 Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; dims=:) = Base._stack(dims, check_stack_dims(f(xy...) for xy in zip(xs, yzs...)))
 Base._stack(dims::Dimension, iter) = Base._stack(typeof(dims), iter)
-Base._stack(dims::Type{<:Dimension}, iter) = Base._stack(dimnum(first(iter), dims), iter)
+Base._stack(dims::Union{Symbol,Type{<:Dimension}}, iter) = Base._stack(dimnum(first(iter), dims), iter)
 
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))

From 16c205a3abfc798eeb311f3b3d0a49f9ba8dbb6f Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Sun, 17 Nov 2024 16:37:44 +1100
Subject: [PATCH 21/27] Fix typo

---
 test/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/methods.jl b/test/methods.jl
index 8ffdac077..e2b7d93f3 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -821,7 +821,7 @@ end
     sx = stack(x; dims=1)
     sy = @test_nowarn stack(x; dims=:)
     sz = @test_nowarn stack(x; dims=X)
-    sw = @test_nowarn stack(X, dims=:X)
+    sw = @test_nowarn stack(x; dims=:X)
     @test sx == sz
     @test sz == sw
     @test sx == stack([parent(da), parent(ca)], dims=1)

From 15cc9fd0285c3e8d500a08c3d6f6c6cfa3f6d4e8 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Tue, 19 Nov 2024 13:50:08 +1100
Subject: [PATCH 22/27] Avoid dispatching on Base._stack

---
 src/array/methods.jl | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 2b56a5f65..1bcf72889 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -534,6 +534,10 @@ function check_stack_dims(iter)
     iter
 end
 
+_maybe_dimnum(x, dim) = hasdim(x, dim) ? dimnum(x, dim) : ndims(x) + 1
+_maybe_dimnum(_, ::Colon) = Colon()
+_maybe_dimnum(x, dims::Tuple) = _maybe_dimnum.([x], dims)
+
 """
     stack(x::AbstractDimArray; [dims])
 
@@ -555,11 +559,15 @@ x = DimArray([a, b], Z(4.0:5.0)) # Construct a nested DimArray
 y = stack(x; dims=2) # Resulting array has dims (X, *Z*, Y)
 ```
 """
-Base.stack(iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, check_stack_dims(iter))
-Base.stack(f, iter::AbstractArray{<:AbstractDimArray}; dims=:) = Base._stack(dims, f(x) for x in check_stack_dims(iter))
-Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; dims=:) = Base._stack(dims, check_stack_dims(f(xy...) for xy in zip(xs, yzs...)))
-Base._stack(dims::Dimension, iter) = Base._stack(typeof(dims), iter)
-Base._stack(dims::Union{Symbol,Type{<:Dimension}}, iter) = Base._stack(dimnum(first(iter), dims), iter)
+function Base.stack(iter::AbstractArray{<:AbstractDimArray}; dims=:)
+    Base._stack(_maybe_dimnum(first(iter), dims), check_stack_dims(iter))
+end
+function Base.stack(f, iter::AbstractArray{<:AbstractDimArray}; kwargs...)
+    Base.stack(f(x) for x in check_stack_dims(iter); kwargs...)
+end
+function Base.stack(f, xs::AbstractArray{<:AbstractDimArray}, yzs...; kwargs...)
+    Base.stack(f(xy...) for xy in zip(xs, yzs...); kwargs...)
+end
 
 function Base.inv(A::AbstractDimArray{T,2}) where T
     newdata = inv(parent(A))

From 52ce34a06d72183590eee138a32cbf67c4ce7d0c Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Tue, 19 Nov 2024 13:51:18 +1100
Subject: [PATCH 23/27] Update tests

---
 test/methods.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index e2b7d93f3..29b660aa2 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -821,9 +821,7 @@ end
     sx = stack(x; dims=1)
     sy = @test_nowarn stack(x; dims=:)
     sz = @test_nowarn stack(x; dims=X)
-    sw = @test_nowarn stack(x; dims=:X)
     @test sx == sz
-    @test sz == sw
     @test sx == stack([parent(da), parent(ca)], dims=1)
     @test sx isa AbstractDimArray
     @test dims(sx) == (dims(x)..., dims(first(x))...)
@@ -834,8 +832,12 @@ end
         @test dims(dc, d) isa AnonDim
         @test parent(dc) == stack([da, db], dims=d)
 
-        @test stack(x -> x .^ 2, x) == stack(parent(x)) .^ 2
-        @test stack(+, x, x, x) == stack(x + x + x)
+        dc = stack(x -> x .^ 2, x; dims=d)
+        @test dims(dc, d) isa AnonDim
+        @test dc == stack(parent(x); dims=d) .^ 2
+        dc = stack(+, x, x, x; dims=d)
+        @test dims(dc, d) isa AnonDim
+        @test dc == stack(x + x + x; dims=d)
     end
 end
 

From f868ddf325f3d9a62e62c66afcb13494334c7a9e Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Tue, 19 Nov 2024 14:00:07 +1100
Subject: [PATCH 24/27] Update test with error on out-of-range dims

---
 test/methods.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/methods.jl b/test/methods.jl
index 29b660aa2..7bc943d4a 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -817,6 +817,8 @@ end
     @test sx isa AbstractDimArray
     @test dims(sx) == (dims(first(x))..., dims(x)...)
 
+    @test_throws "ArgumentError" stack(x; dims=4)
+
     x = DimArray([da, ca], (Dim{:a}(1:2),))
     sx = stack(x; dims=1)
     sy = @test_nowarn stack(x; dims=:)

From 911acdaa1e841b708080e64c1639a1f6c075bdcd Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Tue, 19 Nov 2024 14:03:08 +1100
Subject: [PATCH 25/27] Update _maybe_dimnum to error if specified Integer dims
 are too large

---
 src/array/methods.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 1bcf72889..12b5ab1fc 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -535,6 +535,13 @@ function check_stack_dims(iter)
 end
 
 _maybe_dimnum(x, dim) = hasdim(x, dim) ? dimnum(x, dim) : ndims(x) + 1
+function _maybe_dimnum(x, dim::Integer)
+    if dim < ndims(x) + 2
+        return dim
+    else
+        throw(ArgumentError(LazyString("cannot stack slices ndims(x) = ", ndims(x) + 1, " along dims = ", dim)))
+    end
+end
 _maybe_dimnum(_, ::Colon) = Colon()
 _maybe_dimnum(x, dims::Tuple) = _maybe_dimnum.([x], dims)
 

From 9517d16af1171bdc39ad3c6ce755d06fdb4cc7fc Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Thu, 21 Nov 2024 13:57:02 +1100
Subject: [PATCH 26/27] Use map instead of dot broadcast

Co-authored-by: Rafael Schouten <rafaelschouten@gmail.com>
---
 src/array/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 12b5ab1fc..c8a2fc7ab 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -543,7 +543,7 @@ function _maybe_dimnum(x, dim::Integer)
     end
 end
 _maybe_dimnum(_, ::Colon) = Colon()
-_maybe_dimnum(x, dims::Tuple) = _maybe_dimnum.([x], dims)
+_maybe_dimnum(x, dims::Tuple) = map(d -> _maybe_dimnum(x, d), dims)
 
 """
     stack(x::AbstractDimArray; [dims])

From b32951988186a5de86cb15b84bfd6346c00e0e07 Mon Sep 17 00:00:00 2001
From: Brendan Harris <brendanjohnharris@gmail.com>
Date: Thu, 21 Nov 2024 14:24:28 +1100
Subject: [PATCH 27/27] Update reference.md

---
 docs/src/api/reference.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/api/reference.md b/docs/src/api/reference.md
index 3665d92af..ae4c77302 100644
--- a/docs/src/api/reference.md
+++ b/docs/src/api/reference.md
@@ -105,6 +105,7 @@ Base methods
 
 ```@docs
 Base.cat
+Base.stack
 Base.copy!
 Base.eachslice
 ```