From 99584739f849382dedecead557ed95ec0a760052 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Wed, 15 Jan 2025 00:13:35 +0100
Subject: [PATCH 1/3] add combine method

---
 src/DimensionalData.jl |  2 +-
 src/groupby.jl         | 32 +++++++++++++++++++++++++++++++-
 src/stack/indexing.jl  | 26 ++++++++++++--------------
 src/stack/stack.jl     | 28 +++++++++++++++++++++++++++-
 test/groupby.jl        |  5 +++++
 5 files changed, 76 insertions(+), 17 deletions(-)

diff --git a/src/DimensionalData.jl b/src/DimensionalData.jl
index 7fe8c7cca..44c82fa89 100644
--- a/src/DimensionalData.jl
+++ b/src/DimensionalData.jl
@@ -78,7 +78,7 @@ export dimnum, hasdim, hasselection, otherdims
 export set, rebuild, reorder, modify, broadcast_dims, broadcast_dims!,
     mergedims, unmergedims, maplayers
 
-export groupby, seasons, months, hours, intervals, ranges
+export groupby, combine, seasons, months, hours, intervals, ranges
 
 
 export @d
diff --git a/src/groupby.jl b/src/groupby.jl
index 70f7add8b..6faa1b864 100644
--- a/src/groupby.jl
+++ b/src/groupby.jl
@@ -447,7 +447,8 @@ end
 
 Generate a `Vector` of `UnitRange` with length `step(A)`
 """
-intervals(rng::AbstractRange) = IntervalSets.Interval{:closed,:open}.(rng, rng .+ step(rng))
+intervals(rng::AbstractRange) =
+    IntervalSets.Interval{:closed,:open}.(rng, rng .+ step(rng))
 
 """
     ranges(A::AbstractRange{<:Integer})
@@ -455,3 +456,32 @@ intervals(rng::AbstractRange) = IntervalSets.Interval{:closed,:open}.(rng, rng .
 Generate a `Vector` of `UnitRange` with length `step(A)`
 """
 ranges(rng::AbstractRange{<:Integer}) = map(x -> x:x+step(rng)-1, rng)
+
+
+"""
+    combine(f::Function, gb::DimGroupByArray; dims=:)
+
+Combine the `DimGroupByArray` using funciton `f` over the group dimensions.
+
+If `dims` is given, combine only the dimensions in `dims`. The reducing function
+`f` must accept a `dims` keyword.
+"""
+function combine(f::Function, gb::DimGroupByArray{G}; dims=:) where G
+    # This works for both arrays and stacks
+    # Combine the remaining dimensions after reduction and the group dimensions
+    destdims = (otherdims(DD.dims(first(gb)), dims)..., DD.dims(gb)...)
+    # Get the output eltype 
+    T = Base.promote_op(f, G)
+    # Create a output array with the combined dimensions
+    dest = similar(first(gb), T, destdims)
+    for D in DimIndices(gb)
+        if dims isa Colon
+            # Assigned reduced scalar to dest
+            dest[D...] = f(gb[D])
+        else
+            # Broadcast the reduced array to dest
+            dest[D...] .= f(gb[D]; dims)
+        end
+    end
+    return dest
+end
\ No newline at end of file
diff --git a/src/stack/indexing.jl b/src/stack/indexing.jl
index 19752d97b..f73b9ca11 100644
--- a/src/stack/indexing.jl
+++ b/src/stack/indexing.jl
@@ -147,6 +147,9 @@ for f in (:getindex, :view, :dotview)
     end
 end
 
+@generated function _any_dimarray(v::Union{NamedTuple,Tuple})
+    any(T -> T <: AbstractDimArray, v.types)
+end
 
 #### setindex ####
 @propagate_inbounds Base.setindex!(s::AbstractDimStack, xs, I...; kw...) =
@@ -157,22 +160,17 @@ end
     hassamedims(s) ? _map_setindex!(s, xs, i; kw...) : _setindex_mixed!(s, xs, i; kw...)
 @propagate_inbounds Base.setindex!(s::AbstractDimStack, xs::NamedTuple, i::AbstractArray; kw...) =
     hassamedims(s) ? _map_setindex!(s, xs, i; kw...) : _setindex_mixed!(s, xs, i; kw...)
+@propagate_inbounds Base.setindex!(s::AbstractDimStack, xs::NamedTuple, i::DimensionIndsArrays; kw...) =
+    _map_setindex!(s, xs, i; kw...)
+@propagate_inbounds Base.setindex!(s::AbstractDimStack, xs::NamedTuple, I...; kw...) =
+    _map_setindex!(s, xs, I...; kw...)
 
-@propagate_inbounds function Base.setindex!(
-    s::AbstractDimStack, xs::NamedTuple, I...; kw...
-)
-    map((A, x) -> setindex!(A, x, I...; kw...), layers(s), xs)
-end
-
-_map_setindex!(s, xs, i; kw...) = map((A, x) -> setindex!(A, x, i...; kw...), layers(s), xs)
+_map_setindex!(s, xs, i...; kw...) = map((A, x) -> setindex!(A, x, i...; kw...), layers(s), xs)
 
-_setindex_mixed!(s::AbstractDimStack, x, i::AbstractArray) =
-    map(A -> setindex!(A, x, DimIndices(dims(s))[i]), layers(s))
-_setindex_mixed!(s::AbstractDimStack, i::Integer) =
-    map(A -> setindex!(A, x, DimIndices(dims(s))[i]), layers(s))
-function _setindex_mixed!(s::AbstractDimStack, x, i::Colon)
-    map(DimIndices(dims(s))) do D
-        map(A -> setindex!(A, D), x, layers(s))
+function _setindex_mixed!(s::AbstractDimStack, xs::NamedTuple, i)
+    D = DimIndices(dims(s))[i]
+    map(layers(s), xs) do A, x
+        A[D] = x
     end
 end
 
diff --git a/src/stack/stack.jl b/src/stack/stack.jl
index 597297627..750655798 100644
--- a/src/stack/stack.jl
+++ b/src/stack/stack.jl
@@ -153,7 +153,6 @@ Base.length(s::AbstractDimStack) = prod(size(s))
 Base.axes(s::AbstractDimStack) = map(first ∘ axes, dims(s))
 Base.axes(s::AbstractDimStack, dims::DimOrDimType) = axes(s, dimnum(s, dims))
 Base.axes(s::AbstractDimStack, dims::Integer) = axes(s)[dims]
-Base.similar(s::AbstractDimStack, args...) = maplayers(A -> similar(A, args...), s)
 Base.eltype(::AbstractDimStack{<:Any,T}) where T = T
 Base.ndims(::AbstractDimStack{<:Any,<:Any,N}) where N = N
 Base.CartesianIndices(s::AbstractDimStack) = CartesianIndices(dims(s))
@@ -197,6 +196,33 @@ Base.get(f::Base.Callable, st::AbstractDimStack, k::Symbol) =
 @propagate_inbounds Base.iterate(st::AbstractDimStack, i) =
     i > length(st) ? nothing : (st[DimIndices(st)[i]], i + 1)
 
+Base.similar(s::AbstractDimStack) = similar(s, eltype(s))
+Base.similar(s::AbstractDimStack, dims::Tuple{Vararg{Dimension}}) = 
+    similar(s, eltype(s), dims)
+Base.similar(s::AbstractDimStack, ::Type{T}) where T = 
+    similar(s, T, dims(s))
+function Base.similar(s::AbstractDimStack, ::Type{T}, dims::Tuple) where T
+    # Any dims not in the stack are added to all layers
+    ods = otherdims(s, dims)
+    maplayers(s) do A
+        # Original layer dims are maintained, other dims are added
+        D = DD.commondims(dims, (dims(A)..., ods))
+        similar(A, T, D)
+    end
+end
+function Base.similar(s::AbstractDimStack, ::Type{T}, dims::Tuple) where T<:NamedTuple
+    ods = otherdims(s, dims)
+    maplayers(s, _nt_types(T)) do A, Tx 
+        D = DD.commondims(dims, (DD.dims(A)..., ods))
+        similar(A, Tx, D)
+    end
+end
+
+@generated function _nt_types(::Type{NamedTuple{K,T}}) where {K,T}
+    expr = Expr(:tuple, T.parameters...)
+    return :(NamedTuple{K}($expr))
+end
+
 # `merge` for AbstractDimStack and NamedTuple.
 # One of the first three arguments must be an AbstractDimStack for dispatch to work.
 Base.merge(s::AbstractDimStack) = s
diff --git a/test/groupby.jl b/test/groupby.jl
index 57e1103ed..9611d6c3d 100644
--- a/test/groupby.jl
+++ b/test/groupby.jl
@@ -19,6 +19,10 @@ st = DimStack((a=A, b=A, c=A[X=1]))
         mean(st[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1])
     end
     @test mean.(groupby(st, Ti=>month)) == manualmeans_st
+    combined_st = combine(mean, groupby(st, Ti=>month))
+    @test combined_st isa DimStack{(:a, :b, :c), @NamedTuple{a::Float64, b::Float64, c::Float64}}
+    @test collect(combined_st) == manualmeans_st
+    st[1] = (a= 1, b=2, c=3)
 
     manualsums = mapreduce(hcat, months) do m
         vcat(sum(A[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1, X=1 .. 1.5]), 
@@ -52,6 +56,7 @@ end
     @test mean.(groupby(A, Ti=>Bins(month, ranges(1:3:12)))) == manualmeans
     @test mean.(groupby(A, Ti=>Bins(month, intervals(1:3:12)))) == manualmeans
     @test mean.(groupby(A, Ti=>Bins(month, 4))) == manualmeans
+    @test DimensionalData.combine(mean, groupby(A, Ti=>Bins(month, ranges(1:3:12)))) == manualmeans
 end
 
 @testset "dimension matching groupby" begin

From bc646f597b858a1b6ac8eafdefbf3dcea6e401b9 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Wed, 15 Jan 2025 23:10:32 +0100
Subject: [PATCH 2/3] test groupby and similar

---
 src/array/methods.jl |  2 +-
 src/groupby.jl       | 11 ++++++++---
 src/stack/stack.jl   | 11 +++++++----
 src/utils.jl         |  5 +++++
 test/groupby.jl      | 33 +++++++++++++++++++++++++--------
 test/stack.jl        | 16 ++++++++++++++--
 6 files changed, 60 insertions(+), 18 deletions(-)

diff --git a/src/array/methods.jl b/src/array/methods.jl
index 0088751b9..7ad5a2da6 100644
--- a/src/array/methods.jl
+++ b/src/array/methods.jl
@@ -392,7 +392,7 @@ function _check_cat_lookups(D, ::Regular, lookups...)
             @warn _cat_warn_string(D, "step sizes $(step(span(l))) and $s do not match")
             return false
         end
-        if !(lastval + s ≈ first(l))
+        if !(s isa Dates.AbstractTime) && !(lastval + s ≈ first(l))
             @warn _cat_warn_string(D, "`Regular` lookups do not join with the correct step size: $(lastval) + $s ≈ $(first(l)) should hold")
             return false
         end
diff --git a/src/groupby.jl b/src/groupby.jl
index 6faa1b864..882d2efac 100644
--- a/src/groupby.jl
+++ b/src/groupby.jl
@@ -356,6 +356,7 @@ end
 function _group_indices(dim::Dimension, f::Base.Callable; labels=nothing)
     orig_lookup = lookup(dim)
     k1 = f(first(orig_lookup))
+    # TODO: using a Dict here is a bit slow
     indices_dict = Dict{typeof(k1),Vector{Int}}()
     for (i, x) in enumerate(orig_lookup)
          k = f(x)
@@ -457,7 +458,6 @@ Generate a `Vector` of `UnitRange` with length `step(A)`
 """
 ranges(rng::AbstractRange{<:Integer}) = map(x -> x:x+step(rng)-1, rng)
 
-
 """
     combine(f::Function, gb::DimGroupByArray; dims=:)
 
@@ -467,6 +467,9 @@ If `dims` is given, combine only the dimensions in `dims`. The reducing function
 `f` must accept a `dims` keyword.
 """
 function combine(f::Function, gb::DimGroupByArray{G}; dims=:) where G
+    targetdims = DD.commondims(first(gb), dims)
+    all(hasdim(first(gb), targetdims)) || throw(ArgumentError("dims must be a subset of the groupby dimensions"))
+    all(hasdim(targetdims, DD.dims(gb))) || throw(ArgumentError("grouped dimensions $(DD.basedims(gb)) must be included in dims"))
     # This works for both arrays and stacks
     # Combine the remaining dimensions after reduction and the group dimensions
     destdims = (otherdims(DD.dims(first(gb)), dims)..., DD.dims(gb)...)
@@ -475,12 +478,14 @@ function combine(f::Function, gb::DimGroupByArray{G}; dims=:) where G
     # Create a output array with the combined dimensions
     dest = similar(first(gb), T, destdims)
     for D in DimIndices(gb)
-        if dims isa Colon
+        if all(hasdim(targetdims, DD.dims(first(gb))))
             # Assigned reduced scalar to dest
             dest[D...] = f(gb[D])
         else
+            # Reduce with `f` and drop length 1 dimensions
+            xs = dropdims(f(gb[D]; dims); dims)
             # Broadcast the reduced array to dest
-            dest[D...] .= f(gb[D]; dims)
+            broadcast_dims!(identity, view(dest, D...), xs)
         end
     end
     return dest
diff --git a/src/stack/stack.jl b/src/stack/stack.jl
index 750655798..e5b4e5fe8 100644
--- a/src/stack/stack.jl
+++ b/src/stack/stack.jl
@@ -197,23 +197,26 @@ Base.get(f::Base.Callable, st::AbstractDimStack, k::Symbol) =
     i > length(st) ? nothing : (st[DimIndices(st)[i]], i + 1)
 
 Base.similar(s::AbstractDimStack) = similar(s, eltype(s))
+Base.similar(s::AbstractDimStack, dims::Dimension...) = similar(s, dims)
+Base.similar(s::AbstractDimStack, ::Type{T},dims::Dimension...) where T =
+    similar(s, T, dims)
 Base.similar(s::AbstractDimStack, dims::Tuple{Vararg{Dimension}}) = 
     similar(s, eltype(s), dims)
 Base.similar(s::AbstractDimStack, ::Type{T}) where T = 
     similar(s, T, dims(s))
 function Base.similar(s::AbstractDimStack, ::Type{T}, dims::Tuple) where T
     # Any dims not in the stack are added to all layers
-    ods = otherdims(s, dims)
+    ods = otherdims(dims, DD.dims(s))
     maplayers(s) do A
         # Original layer dims are maintained, other dims are added
-        D = DD.commondims(dims, (dims(A)..., ods))
+        D = DD.commondims(dims, (DD.dims(A)..., ods...))
         similar(A, T, D)
     end
 end
 function Base.similar(s::AbstractDimStack, ::Type{T}, dims::Tuple) where T<:NamedTuple
-    ods = otherdims(s, dims)
+    ods = otherdims(dims, DD.dims(s))
     maplayers(s, _nt_types(T)) do A, Tx 
-        D = DD.commondims(dims, (DD.dims(A)..., ods))
+        D = DD.commondims(dims, (DD.dims(A)..., ods...))
         similar(A, Tx, D)
     end
 end
diff --git a/src/utils.jl b/src/utils.jl
index 8c9f5cb43..f4b23b7db 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -160,6 +160,11 @@ function broadcast_dims!(f, dest::AbstractDimArray{<:Any,N}, As::AbstractBasicDi
     od = map(A -> otherdims(dest, dims(A)), As)
     return _broadcast_dims_inner!(f, dest, As, od)
 end
+function broadcast_dims!(f, dest::AbstractDimStack, stacks::AbstractDimStack...)
+    maplayers(dest, stacks...) do d, layers...
+        broadcast_dims!(f, d, layers...)
+    end
+end
 
 # Function barrier
 function _broadcast_dims_inner!(f, dest, As, od)
diff --git a/test/groupby.jl b/test/groupby.jl
index 9611d6c3d..34ef61f35 100644
--- a/test/groupby.jl
+++ b/test/groupby.jl
@@ -8,21 +8,23 @@ days = DateTime(2000):Day(1):DateTime(2000, 12, 31)
 A = DimArray((1:6) * (1:366)', (X(1:0.2:2), Ti(days)))
 st = DimStack((a=A, b=A, c=A[X=1]))
 
-@testset "manual groupby comparisons" begin
+@testset "manual and groupby comparisons" begin
     # Group by month and even/odd Y axis values
     months = DateTime(2000):Month(1):DateTime(2000, 12, 31)
     manualmeans = map(months) do m
         mean(A[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1])
     end
     @test mean.(groupby(A, Ti=>month)) == manualmeans
+    combinedmeans = combine(mean, groupby(A, Ti=>month))
+    @test combinedmeans isa DimArray
+    @test combinedmeans == manualmeans
     manualmeans_st = map(months) do m
         mean(st[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1])
     end
     @test mean.(groupby(st, Ti=>month)) == manualmeans_st
-    combined_st = combine(mean, groupby(st, Ti=>month))
-    @test combined_st isa DimStack{(:a, :b, :c), @NamedTuple{a::Float64, b::Float64, c::Float64}}
-    @test collect(combined_st) == manualmeans_st
-    st[1] = (a= 1, b=2, c=3)
+    combinedmeans_st = combine(mean, groupby(st, Ti=>month))
+    @test combinedmeans_st isa DimStack{(:a, :b, :c), @NamedTuple{a::Float64, b::Float64, c::Float64}}
+    @test collect(combinedmeans_st) == manualmeans_st
 
     manualsums = mapreduce(hcat, months) do m
         vcat(sum(A[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1, X=1 .. 1.5]), 
@@ -33,6 +35,8 @@ st = DimStack((a=A, b=A, c=A[X=1]))
     @test dims(gb_sum, Ti) == Ti(Sampled([1:12...], ForwardOrdered(), Irregular((nothing, nothing)), Points(), NoMetadata()))
     @test typeof(dims(gb_sum, X)) == typeof(X(Sampled(BitVector([false, true]), ForwardOrdered(), Irregular((nothing, nothing)), Points(), NoMetadata())))
     @test gb_sum == manualsums
+    combined_sum = combine(sum, groupby(A, Ti=>month, X => >(1.5)))
+    @test collect(combined_sum) == manualsums
 
     manualsums_st = mapreduce(hcat, months) do m
         vcat(sum(st[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1, X=1 .. 1.5]), 
@@ -43,10 +47,22 @@ st = DimStack((a=A, b=A, c=A[X=1]))
     @test dims(gb_sum_st, Ti) == Ti(Sampled([1:12...], ForwardOrdered(), Irregular((nothing, nothing)), Points(), NoMetadata()))
     @test typeof(dims(gb_sum_st, X)) == typeof(X(Sampled(BitVector([false, true]), ForwardOrdered(), Irregular((nothing, nothing)), Points(), NoMetadata())))
     @test gb_sum_st == manualsums_st
+    combined_sum_st = combine(sum, groupby(st, Ti=>month, X => >(1.5)))
+    @test collect(combined_sum_st) == manualsums_st
 
     @test_throws ArgumentError groupby(st, Ti=>month, Y=>isodd)
 end
 
+@testset "partial reductions in combine" begin
+    months = DateTime(2000):Month(1):DateTime(2000, 12, 31)
+    using BenchmarkTools
+    manualmeans = cat(map(months) do m
+        mean(A[Ti=dayofyear(m):dayofyear(m)+daysinmonth(m)-1]; dims=Ti)
+    end...; dims=Ti(collect(1:12)))
+    combinedmeans = combine(mean, groupby(A, Ti()=>month); dims=Ti())
+    @test combinedmeans == manualmeans
+end
+
 @testset "bins" begin
     seasons = DateTime(2000):Month(3):DateTime(2000, 12, 31)
     manualmeans = map(seasons) do s
@@ -56,7 +72,7 @@ end
     @test mean.(groupby(A, Ti=>Bins(month, ranges(1:3:12)))) == manualmeans
     @test mean.(groupby(A, Ti=>Bins(month, intervals(1:3:12)))) == manualmeans
     @test mean.(groupby(A, Ti=>Bins(month, 4))) == manualmeans
-    @test DimensionalData.combine(mean, groupby(A, Ti=>Bins(month, ranges(1:3:12)))) == manualmeans
+    @test combine(mean, groupby(A, Ti=>Bins(month, ranges(1:3:12)))) == manualmeans
 end
 
 @testset "dimension matching groupby" begin
@@ -73,9 +89,10 @@ end
     end
     @test all(collect(mean.(gb)) .=== manualmeans)
     @test all(mean.(gb) .=== manualmeans)
+    @test all(combine(mean, gb) .=== manualmeans)
 end
 
-@testset "broadcastdims runs after groupby" begin
+@testset "broadcast_dims runs after groupby" begin
     dimlist = (
         Ti(Date("2021-12-01"):Day(1):Date("2022-12-31")),
         X(range(1, 10, length=10)),
@@ -85,7 +102,7 @@ end
     data = rand(396, 10, 15, 2)
     A = DimArray(data, dimlist)
     month_length = DimArray(daysinmonth, dims(A, Ti))
-    g_tempo = DimensionalData.groupby(month_length, Ti=>seasons(; start=December))
+    g_tempo = DimensionalData.groupby(month_length, Ti => seasons(; start=December))
     sum_days = sum.(g_tempo, dims=Ti)
     @test sum_days isa DimArray
     weights = map(./, g_tempo, sum_days)
diff --git a/test/stack.jl b/test/stack.jl
index cfde24449..ee6d2eac0 100644
--- a/test/stack.jl
+++ b/test/stack.jl
@@ -3,7 +3,7 @@ using DimensionalData, Test, LinearAlgebra, Statistics, ConstructionBase, Random
 using DimensionalData: data
 using DimensionalData: Sampled, Categorical, AutoLookup, NoLookup, Transformed,
     Regular, Irregular, Points, Intervals, Start, Center, End,
-    Metadata, NoMetadata, ForwardOrdered, ReverseOrdered, Unordered, layers, basedims
+    Metadata, NoMetadata, ForwardOrdered, ReverseOrdered, Unordered, layers, basedims, layerdims
 
 A = [1.0 2.0 3.0;
      4.0 5.0 6.0]
@@ -94,11 +94,23 @@ end
     @test all(maplayers(similar(mixed), mixed) do s, m
         dims(s) == dims(m) && dims(s) === dims(m) && eltype(s) === eltype(m)
     end)
-    @test eltype(similar(s, Int)) === @NamedTuple{one::Int, two::Int, three::Int}
+    @test eltype(similar(s, Int)) === 
+    @NamedTuple{one::Int, two::Int, three::Int}
+    @test eltype(similar(s, @NamedTuple{one::Int, two::Float32, three::Bool})) === 
+        @NamedTuple{one::Int, two::Float32, three::Bool}
     st2 = similar(mixed, Bool, x, y)
     @test dims(st2) === (x, y)
     @test dims(st2[:one]) === (x, y)
     @test eltype(st2) === @NamedTuple{one::Bool, two::Bool, extradim::Bool}
+    @test eltype(similar(mixed)) == eltype(mixed)
+    @test size(similar(mixed)) == size(mixed)
+    @test keys(similar(mixed)) == keys(mixed)
+    @test layerdims(similar(mixed)) == layerdims(mixed)
+    xy = (X(), Y()) 
+    @test layerdims(similar(mixed, dims(mixed, (X, Y)))) == (one=xy, two=xy, extradim=xy)
+    st3 = similar(mixed, @NamedTuple{one::Int, two::Float32, extradim::Bool}, (Z([:a, :b, :c]), Ti(1:12), X(1:3)))
+    @test layerdims(st3) == (one=(Ti(), X()), two=(Ti(), X()), extradim=(Z(), Ti(), X()))
+    @test eltype(st3) == @NamedTuple{one::Int, two::Float32, extradim::Bool}
 end
 
 @testset "merge" begin

From d07239767451fef48811723dba8d046802bd5a61 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 19 Jan 2025 17:55:04 +0100
Subject: [PATCH 3/3] docs entry

---
 docs/src/api/reference.md |  1 +
 src/groupby.jl            | 18 ++++++++++++++----
 2 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/docs/src/api/reference.md b/docs/src/api/reference.md
index 3665d92af..2427fed92 100644
--- a/docs/src/api/reference.md
+++ b/docs/src/api/reference.md
@@ -58,6 +58,7 @@ For transforming DimensionalData objects:
 
 ```@docs
 groupby
+combine
 DimensionalData.DimGroupByArray
 Bins
 ranges
diff --git a/src/groupby.jl b/src/groupby.jl
index 882d2efac..1cdd9e509 100644
--- a/src/groupby.jl
+++ b/src/groupby.jl
@@ -249,7 +249,6 @@ Group some data along the time dimension:
 
 ```jldoctest groupby; setup = :(using Random; Random.seed!(123))
 julia> using DimensionalData, Dates
-
 julia> A = rand(X(1:0.1:20), Y(1:20), Ti(DateTime(2000):Day(3):DateTime(2003)));
 
 julia> groups = groupby(A, Ti => month) # Group by month
@@ -461,10 +460,21 @@ ranges(rng::AbstractRange{<:Integer}) = map(x -> x:x+step(rng)-1, rng)
 """
     combine(f::Function, gb::DimGroupByArray; dims=:)
 
-Combine the `DimGroupByArray` using funciton `f` over the group dimensions.
+Combine the `DimGroupByArray` using function `f` over the group dimensions.
+Unlike broadcasting a reducing function over a `DimGroupByArray`, this function
+always returns a new flattened `AbstractDimArray` even where not all dimensions 
+are reduced. It will also work over grouped `AbstractDimStack`.
+
+If `dims` is given, it will combine only the dimensions in `dims`, the 
+others will be present in the final array. Note that all grouped dimensions
+must be reduced and included in `dims`.
+
+The reducing function `f` must also accept a `dims` keyword.
 
-If `dims` is given, combine only the dimensions in `dims`. The reducing function
-`f` must accept a `dims` keyword.
+# Example
+
+```jldoctest groupby
+````
 """
 function combine(f::Function, gb::DimGroupByArray{G}; dims=:) where G
     targetdims = DD.commondims(first(gb), dims)