Add a groupby function (#591)

* add groupby

* ambiguity fixes

* update test

* fix groupby more

Author: rafaqz

Project.toml

@@ -7,6 +7,7 @@ version = "0.25.8"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 Interfaces = "85a1e053-f937-4924-92a5-1367d23b7b87"
@@ -38,6 +39,7 @@ ColorTypes = "0.11"
 Combinatorics = "1"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6"
+DataAPI = "1"
 DataFrames = "1"
 Dates = "1"
 Distributions = "0.25"

docs/Project.toml

@@ -6,3 +6,5 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

docs/crash/course/groupby.jl

@@ -0,0 +1,115 @@
+using DimensionalData
+using Dates
+using Statistics
+const DD = DimensionalData
+
+# # Basics: DateTime operations we can use for grouping
+
+# First lets look at the kind of functions that can be used to group `DateTime`.
+# Other types will follow the same principles, but are usually simpler.
+
+# Create a demo `DateTime` range
+tempo = range(DateTime(2000), step=Hour(1), length=365*24*2)
+
+# Now lets see how some common functions work.
+
+# The `hour` function will transform values to hour of the day - the integers `0:23`
+
+# ### hour
+hour.(tempo)
+
+# These do similar things with other time periods
+
+# ### dayofweek
+dayofweek.(tempo)
+
+# ### month
+month.(tempo)
+
+# ### dayofyear
+dayofyear.(tempo)
+
+# ## Tuple grouping
+
+# Some functions return a tuple - we can also use tuples for grouping.
+# They are sorted by the left to right values.
+
+# ### yearmonth
+yearmonth.(tempo)
+
+# We can creat our own anonymous function that return tuples
+yearday(x) = year(x), dayofyear(x)
+
+yearhour(x) = year(x), hour(x)
+
+# And you can probably guess what they do:
+yearday.(tempo)
+
+# All of these functions can be used in `groupby` on `DateTime` objects.
+
+# # Practical example: grouping by season
+
+# ### TODOS: We will need several functions.
+
+# # groupby operations
+# Here we use the same time functions from above
+
+ds = rand(X(1:0.01:2), Ti(tempo))
+
+# ## select by month, days, years and seasons
+# ### TODO, how do we select month 1 or 2, and even a group of them, i.e. [1,3,5]? Same for days, years and seasons.
+mean.(groupby(ds, Ti=>Bins(month, 1:2))) 
+mean.(groupby(ds, Ti=>Bins(month, [1, 3, 5]))) 
+mean.(groupby(ds, Ti => season(; start=December))) 
+mean.(groupby(ds, Ti => Bins(dayofyear, intervals(1:8:370))))
+mean.(groupby(ds, Ti => Bins(yearday, [[1,2,3], [4,5,6]], labels=x -> join(string.(x), ','))))
+mean.(groupby(ds, Ti => week))
+mean.(groupby(ds, Ti => hours(12; start=6, labels=x -> 6 in x ? :night : :day)))
+mean.(groupby(ds, Ti => dims(ds, Ti)))
+
+# ### TODO, we need a new function that can return DJF (Dec-Jan-Feb), MAM (Mar-Apr-May)... etc.
+# THIS IS HARD. We need a succinct way to select around the end-start of the year. 
+
+# is combining month from different years
+mean.(groupby(ds, Ti=>month)) 
+
+# Use three-month bins. The 13 is the open side of the last interval.
+mean.(groupby(ds, Ti=>Bins(yearmonth, intervals(1:3:12))))
+
+mean.(groupby(ds, Ti=>Bins(month, 4))) # is combining month from different years
+
+#
+mean.(groupby(ds, Ti=>year))
+
+#
+mean.(groupby(ds, Ti=>yearmonth))
+
+#
+mean.(groupby(ds, Ti=>hour))
+
+#
+mean.(groupby(ds, Ti=>Dates.hour12))
+
+# ### TODO. How do could we incorporate resample? Let's say if we have hour resolution I want to resample every 3,6,12.. hours?
+mean.(groupby(ds, Ti=>Bins(yearhour, intervals(1:3:24)))) # it will combine the same day from different year.
+
+mean.(groupby(ds, Ti=>dayofyear)) # it will combine the same day from different year.
+
+#
+mean.(groupby(ds, Ti=>yearmonthday)) # this does the a daily mean aggregation.
+
+#
+mean.(groupby(ds, Ti=>yearmonth)) # this does a monthly mean aggregation
+
+# 
+mean.(groupby(ds, Ti=>yearday)) # this does a daily mean aggregation
+
+mean.(groupby(ds, Ti=>Bins(yearhour, 12))) # this does a daily mean aggregation
+
+# ### TODO. Similar to the hourly resample, how do we do it for more than 1 day, let's say 8daily?
+mean.(groupby(ds, Ti=>Bins(dayofyear, map(x -> x:x+7, 1:8:370))))
+
+# ### TODO: Group by Dims. This should include the rasters input sampling.
+mean.(groupby(ds, dims(ds, Ti)))
+
+# ## Apply custom function (i.e. normalization) to grouped output.

src/DimensionalData.jl

@@ -17,14 +17,16 @@ using Dates,
 using Base.Broadcast: Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle,
       Unknown
 
-using Base: tail, OneTo, @propagate_inbounds
+using Base: tail, OneTo, Callable, @propagate_inbounds
 
 # Ecosystem
 import Adapt, 
        ArrayInterface,
        ConstructionBase, 
+       DataAPI,
        Extents,
        Interfaces,
+       IntervalSets,
        InvertedIndices,
        IteratorInterfaceExtensions,
        RecipesBase,
@@ -34,23 +36,27 @@ import Adapt,
 
 using RecipesBase: @recipe
 
+# using IntervalSets: .., Interval
+
 include("Dimensions/Dimensions.jl")
 
 using .Dimensions
 using .Dimensions.LookupArrays
 using .Dimensions: StandardIndices, DimOrDimType, DimTuple, DimTupleOrEmpty, DimType, AllDims
 import .LookupArrays: metadata, set, _set, rebuild, basetypeof, 
-    order, span, sampling, locus, val, index, bounds, intervalbounds, 
+    order, span, sampling, locus, val, index, bounds, intervalbounds,
     hasselection, units, SelectorOrInterval
 import .Dimensions: dims, refdims, name, lookup, dimstride, kwdims, hasdim, label, _astuple
 
+import DataAPI.groupby
+
 export LookupArrays, Dimensions
 
 # Dimension
 export X, Y, Z, Ti, Dim, Coord
 
 # Selector
-export At, Between, Touches, Contains, Near, Where, All, .., Not
+export At, Between, Touches, Contains, Near, Where, All, .., Not, Bins
 
 export AbstractDimArray, DimArray
 
@@ -71,6 +77,8 @@ export dimnum, hasdim, hasselection, otherdims
 # utils
 export set, rebuild, reorder, modify, broadcast_dims, broadcast_dims!, mergedims, unmergedims
 
+export groupby, season, months, hours, yeardays, monthdays, intervals, ranges
+
 const DD = DimensionalData
 
 # Common
@@ -96,6 +104,7 @@ include("tables.jl")
 include("plotrecipes.jl")
 include("utils.jl")
 include("set.jl")
+include("groupby.jl")
 include("precompile.jl")
 include("interface_tests.jl")
 

src/Dimensions/Dimensions.jl

@@ -13,7 +13,7 @@ using DimensionalData.Dimensions
 """
 module Dimensions
 
-import Adapt, ConstructionBase, Extents
+import Adapt, ConstructionBase, Extents, IntervalSets
 using Dates 
 
 include("../LookupArrays/LookupArrays.jl")

src/Dimensions/format.jl

@@ -26,11 +26,12 @@ function format(dims::Tuple{<:Pair,Vararg{Pair}}, A::AbstractArray)
     end
     return format(dims, A)
 end
-format(dims::Tuple{Vararg{Any,N}}, A::AbstractArray{<:Any,N}) where N =
-    format(dims, axes(A))
+format(dims::Tuple{Vararg{Any,N}}, A::AbstractArray{<:Any,N}) where N = format(dims, axes(A))
 @noinline format(dims::Tuple{Vararg{Any,M}}, A::AbstractArray{<:Any,N}) where {N,M} =
     throw(DimensionMismatch("Array A has $N axes, while the number of dims is $M: $(map(basetypeof, dims))"))
 format(dims::Tuple{Vararg{Any,N}}, axes::Tuple{Vararg{Any,N}}) where N = map(_format, dims, axes)
+format(d::Dimension{<:AbstractArray}) = _format(d, axes(val(d), 1))
+format(d::Dimension, axis::AbstractRange) = _format(d, axis)
 
 _format(dimname::Symbol, axis::AbstractRange) = Dim{dimname}(NoLookup(axes(axis, 1)))
 _format(::Type{D}, axis::AbstractRange) where D<:Dimension = D(NoLookup(axes(axis, 1)))
@@ -110,8 +111,13 @@ _format(span::Irregular{<:Tuple}, D, index) = span
 _format(span::Explicit, D, index) = span
 # Sampling
 _format(sampling::AutoSampling, span::Span, D::Type, index) = Points()
+_format(::AutoSampling, ::Span, D::Type, ::AbstractArray{<:IntervalSets.Interval}) =
+    Intervals(Start())
 _format(sampling::AutoSampling, span::Explicit, D::Type, index) =
     Intervals(_format(locus(sampling), D, index))
+# For ambiguity, not likely to happen in practice
+_format(::AutoSampling, ::Explicit, D::Type, ::AbstractArray{<:IntervalSets.Interval}) =
+    Intervals(_format(locus(sampling), D, index))
 _format(sampling::Points, span::Span, D::Type, index) = sampling
 _format(sampling::Points, span::Explicit, D::Type, index) = _explicitpoints_error() 
 _format(sampling::Intervals, span::Span, D::Type, index) =

src/Dimensions/primitives.jl

@@ -203,7 +203,7 @@ julia> dimnum(A, Y)
 ```
 """
 @inline function dimnum(x, q1, query...)
-    all(hasdim(x, q1, query...)) || _extradimserror()
+    all(hasdim(x, q1, query...)) || _extradimserror(otherdims(x, (q1, query)))
     _call_primitive(_dimnum, MaybeFirst(), x, q1, query...)
 end
 @inline dimnum(x, query::Function) =
@@ -757,5 +757,5 @@ _typemsg(a, b) = "Lookups do not all have the same type: $(order(a)), $(order(b)
 @noinline _valerror(a, b) = throw(DimensionMismatch(_valmsg(a, b)))
 @noinline _ordererror(a, b) = throw(DimensionMismatch(_ordermsg(a, b)))
 @noinline _metadataerror(a, b) = throw(DimensionMismatch(_metadatamsg(a, b)))
-@noinline _extradimserror(args...) = throw(ArgumentError(_extradimsmsg(args)))
+@noinline _extradimserror(args) = throw(ArgumentError(_extradimsmsg(args)))
 @noinline _dimsnotdefinederror() = throw(ArgumentError("Object does not define a `dims` method"))

src/Dimensions/show.jl

@@ -13,25 +13,27 @@ function dimcolors(i)
     colors[c]
 end
 
-function Base.show(io::IO, mime::MIME"text/plain", dims::DimTuple)
+function show_dims(io::IO, mime::MIME"text/plain", dims::DimTuple;
+    colors=map(x -> get(io, :dimcolo, dimcolors(x)), ntuple(identity, length(dims)))
+)
     ctx = IOContext(io, :compact => true)
     inset = get(io, :inset, "")
     print(io, inset)
     if all(map(d -> !(parent(d) isa AbstractArray) || (parent(d) isa NoLookup), dims))
-        dc = get(ctx, :dimcolor, dimcolors(1))
+        dc = colors[1]
         printstyled(io, dimsymbols(1), ' '; color=dc)
         show(IOContext(ctx, :dimcolor => dc, :dimname_len => 0), mime, first(dims))
         foreach(enumerate(Base.tail(dims))) do (i, d)
             n = i + 1
             print(io, ", ")
-            dc = get(ctx, :dimcolor, dimcolors(n))
+            dc = colors[n]
             printstyled(io, dimsymbols(n), ' '; color=dc)
             show(IOContext(ctx, :dimcolor => dc, :dimname_len => 0), mime, d)
         end
         return 0
     else # Dims get a line each
         lines = 3
-        dc = get(ctx, :dimcolor, dimcolors(1))
+        dc = colors[1]
         printstyled(io, dimsymbols(1), ' '; color=dc)
         maxname = maximum(length ∘ string ∘ dim2key, dims) 
         dim_ctx = IOContext(ctx, :dimcolor => dc, :dimname_len=> maxname)
@@ -41,14 +43,17 @@ function Base.show(io::IO, mime::MIME"text/plain", dims::DimTuple)
             lines += 1
             s = dimsymbols(n)
             print(io, ",\n", inset)
-            dc = get(ctx, :dimcolor, dimcolors(n))
+            dc = colors[n]
             printstyled(io, s, ' '; color=dc)
             dim_ctx = IOContext(ctx, :dimcolor => dc, :dimname_len => maxname)
             show(dim_ctx, mime, d)
         end
         return lines
     end
 end
+
+Base.show(io::IO, mime::MIME"text/plain", dims::DimTuple) =
+    show_dims(io, mime, dims)
 function Base.show(io::IO, mime::MIME"text/plain", dim::Dimension)
     get(io, :compact, false) && return show_compact(io, mime, dim)
     print_dimname(io, dim)
@@ -77,14 +82,14 @@ end
 dimcolor(io) = get(io, :dimcolor, dimcolors(1))
 
 # print dims with description string and inset
-function print_dims(io::IO, mime, dims::Tuple{})
+function print_dims(io::IO, mime, dims::Tuple{}; kw...)
     @nospecialize io mime dims
     print(io, ": ")
     return 0
 end
-function print_dims(io::IO, mime, dims::Tuple)
+function print_dims(io::IO, mime, dims::Tuple; kw...)
     ctx = IOContext(io, :inset => "  ")
-    return show(ctx, mime, dims)
+    return show_dims(ctx, mime, dims; kw...)
 end
 # print refdims with description string and inset
 function print_refdims(io::IO, mime, refdims::Tuple)

src/LookupArrays/lookup_arrays.jl

@@ -46,11 +46,11 @@ ordered_lastindex(o::ForwardOrdered, l::LookupArray) = lastindex(parent(l))
 ordered_lastindex(o::ReverseOrdered, l::LookupArray) = firstindex(parent(l))
 ordered_lastindex(o::Unordered, l::LookupArray) = lastindex(parent(l))
 
-function Base.searchsortedfirst(lookup::LookupArray, val; lt=<)
-    searchsortedfirst(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt)
+function Base.searchsortedfirst(lookup::LookupArray, val; lt=<, kw...)
+    searchsortedfirst(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt, kw...)
 end
-function Base.searchsortedlast(lookup::LookupArray, val; lt=<)
-    searchsortedlast(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt)
+function Base.searchsortedlast(lookup::LookupArray, val; lt=<, kw...)
+    searchsortedlast(parent(lookup), unwrap(val); order=ordering(order(lookup)), lt=lt, kw...)
 end
 
 function Adapt.adapt_structure(to, l::LookupArray)
@@ -593,7 +593,8 @@ Base.:(==)(l1::Transformed, l2::Transformed) = typeof(l1) == typeof(l2) && f(l1)
 
 intervalbounds(l::LookupArray, args...) = _intervalbounds_no_interval_error()
 intervalbounds(l::AbstractSampled, args...) = intervalbounds(span(l), sampling(l), l, args...)
-intervalbounds(span::Span, ::Points, l::LookupArray, args...) = _intervalbounds_no_interval_error()
+intervalbounds(span::Span, ::Points, ls::LookupArray) = map(l -> (l, l), ls) 
+intervalbounds(span::Span, ::Points, ls::LookupArray, i::Int) = ls[i], ls[i]
 intervalbounds(span::Span, sampling::Intervals, l::LookupArray, i::Int) =
     intervalbounds(order(l), locus(sampling), span, l, i)
 function intervalbounds(order::ForwardOrdered, locus::Start, span::Span, l::LookupArray, i::Int)

src/LookupArrays/metadata.jl

@@ -13,7 +13,7 @@ objects to [`set`](@ref) without ambiguity about where to put them.
 """
 abstract type AbstractMetadata{X,T} end
 
-const _MetadataContents =Union{AbstractDict,NamedTuple}
+const _MetadataContents = Union{AbstractDict,NamedTuple}
 const AllMetadata = Union{AbstractMetadata,AbstractDict}
 
 Base.get(m::AbstractMetadata, args...) = get(val(m), args...)