DimensionalData.jl icon indicating copy to clipboard operation
DimensionalData.jl copied to clipboard

Weighted mean from `StatsBase` does not work on DimArray

Open bgroenks96 opened this issue 1 year ago • 6 comments

The following does not work with DimArray:

using DimensionalData, StatsBase

data = ones(X(1:10), Y(1:3))
mean(data, weights([0.3,0.3,0.4]), 2)

This produces an error:

ERROR: ArgumentError: No method is implemented for reducing index range of type DimensionalData.Dimensions.DimUnitRange{Int64, Base.OneTo{Int64}, Y{DimensionalData.Dimensions.Lookups.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}. Please implement
reduced_index for this index type or report this as an issue.

Stacktrace:
 [1] reduced_index(i::DimensionalData.Dimensions.DimUnitRange{Int64, Base.OneTo{Int64}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}})
   @ Base ./reducedim.jl:8
 [2] reduced_indices(inds::Tuple{DimensionalData.Dimensions.DimUnitRange{…}, DimensionalData.Dimensions.DimUnitRange{…}}, d::Int64)
   @ Base ./reducedim.jl:25
 [3] _mean
   @ ~/.julia/packages/StatsBase/ebrT3/src/weights.jl:684 [inlined]
 [4] mean
   @ ~/.julia/packages/StatsBase/ebrT3/src/weights.jl:680 [inlined]
 [5] mean(A::DimMatrix{…}, w::Weights{…}, dims::Int64)
   @ StatsBase ./deprecated.jl:105
 [6] top-level scope
   @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.

On a related note, maybe DimensionalData could add dispatches for these functions that takes a Dim instead of an integer?

bgroenks96 avatar Jul 11 '24 23:07 bgroenks96

We got a bit tricky with those unit ranges, guess we missed a method.

It would be good to allow a dim there too, do we need a StatsBase extension for dispatch?

rafaqz avatar Jul 12 '24 04:07 rafaqz

Yes because the second argument has to be typed as AbstractWeights which comes from StatsBase.

bgroenks96 avatar Jul 12 '24 09:07 bgroenks96

Ok we should add an extension for StatsBase. We need it to make sample dimensional too.

If you have time feel free to PR, otherwise it may take a while to get to it.

rafaqz avatar Jul 12 '24 11:07 rafaqz

So turns out this can be fixed by defining Base.reduced_index

using DimensionalData, StatsBase

function Base.reduced_index(dur::DimensionalData.Dimensions.DimUnitRange) 
    r1 = Base.reduced_index(parent(dur))
    d = dims(dur)
    # Type changes are not allowed here so we just 
    # take the "first" value
    d1 = if isreverse(d)
        d[end:end]
    else
        d[begin:begin]
    end
    DimensionalData.Dimensions.DimUnitRange(r1, d1)
end

And weighted mean just works.

julia> mean(data, weights([0.3,0.3,0.4]), 2)
╭──────────────────────────╮
│ 10×1 DimArray{Float64,2} │
├──────────────────────────┴─────────────────────── dims ┐
  ↓ X Sampled{Int64} 1:10 ForwardOrdered Regular Points,
  → Y Sampled{Int64} 1:1 ForwardOrdered Regular Points
└────────────────────────────────────────────────────────┘
  ↓ →  1
  1    1.0
  2    1.0
  3    1.0
  4    1.0
  5    1.0
  6    1.0
  7    1.0
  8    1.0
  9    1.0
 10    1.0

The annoying part is we cant use reducelookup here because it often changes the lookup type in some way and that errors. So we are just reducing the lookup to the first value in whatever order its in. Not quite satisfying, as reducelookup normally keeps the bounds the same after the reduction.

rafaqz avatar Aug 17 '24 20:08 rafaqz

I am not sure I see the problem here from the example. But I think the suboptimal behavior of the lookup reduction is probably acceptable in exchange for having the basic functionality of reduced_index. The user can always manually redefine the lookup as a workaround. We can put the fix for reducelookup in a new issue.

bgroenks96 avatar Aug 19 '24 13:08 bgroenks96

Yeah let's just add it as is

rafaqz avatar Aug 20 '24 10:08 rafaqz

Fixed

rafaqz avatar Nov 02 '24 12:11 rafaqz