Adjust initialization in maximum and minimum (version 2) (#28116)

* Adjust initialization in maximum and minimum

Fixes #27836

* Reduce to first element in index range for OffsetArrays

* Adjust mapreduce to use first element in offset range as well.

Adjust tests

* Avoid defining global Areduc variable in reducedim tests

* Add comment in mapslices to explain the wrapping of scalars

* Change reduced_index for Slices and add test for previously failing
maximum of OffsetArray example

* Add checks that reductions along dimensions are inferred

base/abstractarray.jl

@@ -1926,7 +1926,12 @@ function mapslices(f, A::AbstractArray; dims)
     Rsize = copy(dimsA)
     # TODO: maybe support removing dimensions
     if !isa(r1, AbstractArray) || ndims(r1) == 0
-        r1 = [r1]
+        # If the result of f on a single slice is a scalar then we add singleton
+        # dimensions. When adding the dimensions, we have to respect the
+        # index type of the input array (e.g. in the case of OffsetArrays)
+        tmp = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
+        tmp[firstindex(tmp)] = r1
+        r1 = tmp
     end
     nextra = max(0, length(dims)-ndims(r1))
     if eltype(Rsize) == Int

base/reducedim.jl

@@ -3,28 +3,41 @@
 ## Functions to compute the reduced shape
 
 # for reductions that expand 0 dims to 1
+reduced_index(i::OneTo) = OneTo(1)
+reduced_index(i::Slice) = first(i):first(i)
+reduced_index(i::AbstractUnitRange) =
+    throw(ArgumentError(
+"""
+No method is implemented for reducing index range of type $typeof(i). Please implement
+reduced_index for this index type or report this as an issue.
+"""
+    ))
 reduced_indices(a::AbstractArray, region) = reduced_indices(axes(a), region)
 
 # for reductions that keep 0 dims as 0
 reduced_indices0(a::AbstractArray, region) = reduced_indices0(axes(a), region)
 
-function reduced_indices(inds::Indices{N}, d::Int, rd::AbstractUnitRange) where N
+function reduced_indices(inds::Indices{N}, d::Int) where N
     d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
     if d == 1
-        return (oftype(inds[1], rd), tail(inds)...)
+        return (reduced_index(inds[1]), tail(inds)...)
     elseif 1 < d <= N
-        return tuple(inds[1:d-1]..., oftype(inds[d], rd), inds[d+1:N]...)::typeof(inds)
+        return tuple(inds[1:d-1]..., oftype(inds[d], reduced_index(inds[d])), inds[d+1:N]...)::typeof(inds)
     else
         return inds
     end
 end
-reduced_indices(inds::Indices, d::Int) = reduced_indices(inds, d, OneTo(1))
 
 function reduced_indices0(inds::Indices{N}, d::Int) where N
     d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
     if d <= N
         ind = inds[d]
-        return reduced_indices(inds, d, (isempty(ind) ? ind : OneTo(1)))
+        rd = isempty(ind) ? ind : reduced_index(inds[d])
+        if d == 1
+            return (rd, tail(inds)...)
+        else
+            return tuple(inds[1:d-1]..., oftype(inds[d], rd), inds[d+1:N]...)::typeof(inds)
+        end
     else
         return inds
     end
@@ -38,7 +51,7 @@ function reduced_indices(inds::Indices{N}, region) where N
         if d < 1
             throw(ArgumentError("region dimension(s) must be ≥ 1, got $d"))
         elseif d <= N
-            rinds[d] = oftype(rinds[d], OneTo(1))
+            rinds[d] = reduced_index(rinds[d])
         end
     end
     tuple(rinds...)::typeof(inds)
@@ -53,7 +66,7 @@ function reduced_indices0(inds::Indices{N}, region) where N
             throw(ArgumentError("region dimension(s) must be ≥ 1, got $d"))
         elseif d <= N
             rind = rinds[d]
-            rinds[d] = oftype(rind, (isempty(rind) ? rind : OneTo(1)))
+            rinds[d] = isempty(rind) ? rind : reduced_index(rind)
         end
     end
     tuple(rinds...)::typeof(inds)
@@ -79,25 +92,6 @@ end
 reducedim_initarray(A::AbstractArray, region, init, ::Type{R}) where {R} = fill!(similar(A,R,reduced_indices(A,region)), init)
 reducedim_initarray(A::AbstractArray, region, init::T) where {T} = reducedim_initarray(A, region, init, T)
 
-function reducedim_initarray0(A::AbstractArray{T}, region, f, ops) where T
-    ri = reduced_indices0(A, region)
-    if isempty(A)
-        if prod(length, reduced_indices(A, region)) != 0
-            reducedim_initarray0_empty(A, region, f, ops) # ops over empty slice of A
-        else
-            R = f == identity ? T : Core.Compiler.return_type(f, (T,))
-            similar(A, R, ri)
-        end
-    else
-        R = f == identity ? T : typeof(f(first(A)))
-        si = similar(A, R, ri)
-        mapfirst!(f, si, A)
-    end
-end
-
-reducedim_initarray0_empty(A::AbstractArray, region, f, ops) = mapslices(x->ops(f.(x)), A, dims = region)
-reducedim_initarray0_empty(A::AbstractArray, region,::typeof(identity), ops) = mapslices(ops, A, dims = region)
-
 # TODO: better way to handle reducedim initialization
 #
 # The current scheme is basically following Steven G. Johnson's original implementation
@@ -124,8 +118,33 @@ function _reducedim_init(f, op, fv, fop, A, region)
     return reducedim_initarray(A, region, z, Tr)
 end
 
-reducedim_init(f, op::typeof(max), A::AbstractArray{T}, region) where {T} = reducedim_initarray0(A, region, f, maximum)
-reducedim_init(f, op::typeof(min), A::AbstractArray{T}, region) where {T} = reducedim_initarray0(A, region, f, minimum)
+# initialization when computing minima and maxima requires a little care
+for (f1, f2, initval) in ((:min, :max, :Inf), (:max, :min, :(-Inf)))
+    @eval function reducedim_init(f, op::typeof($f1), A::AbstractArray, region)
+        # First compute the reduce indices. This will throw an ArgumentError
+        # if any region is invalid
+        ri = reduced_indices(A, region)
+
+        # Next, throw if reduction is over a region with length zero
+        any(i -> isempty(axes(A, i)), region) && _empty_reduce_error()
+
+        # Make a view of the first slice of the region
+        A1 = view(A, ri...)
+
+        if isempty(A1)
+            # If the slice is empty just return non-view version as the initial array
+            return copy(A1)
+        else
+            # otherwise use the min/max of the first slice as initial value
+            v0 = mapreduce(f, $f2, A1)
+
+            # but NaNs need to be avoided as intial values
+            v0 = v0 != v0 ? typeof(v0)($initval) : v0
+
+            return reducedim_initarray(A, region, v0)
+        end
+    end
+end
 reducedim_init(f::Union{typeof(abs),typeof(abs2)}, op::typeof(max), A::AbstractArray{T}, region) where {T} =
     reducedim_initarray(A, region, zero(f(zero(T))))
 

stdlib/SparseArrays/src/sparsematrix.jl

@@ -1616,8 +1616,6 @@ end
 # non-trivial, so use Arrays for output
 Base.reducedim_initarray(A::SparseMatrixCSC, region, v0, ::Type{R}) where {R} =
     fill(v0, Base.reduced_indices(A,region))
-Base.reducedim_initarray0(A::SparseMatrixCSC, region, v0, ::Type{R}) where {R} =
-    fill(v0, Base.reduced_indices0(A,region))
 
 # General mapreduce
 function _mapreducezeros(f, op, ::Type{T}, nzeros::Int, v0) where T

stdlib/SparseArrays/test/higherorderfns.jl

@@ -628,5 +628,8 @@ end
     @test Diagonal(spzeros(5)) \ view(rand(10), 1:5) == [Inf,Inf,Inf,Inf,Inf]
 end
 
+@testset "Issue #27836" begin
+    @test minimum(sparse([1, 2], [1, 2], ones(Int32, 2)), dims = 1) isa Matrix
+end
 
 end # module

test/offsetarray.jl

@@ -213,7 +213,7 @@ targets2 = ["(1.0, 1.0)",
             "([1.0], [1.0])",
             "([1.0], [1.0])",
             "([1.0], [1.0])"]
-for n = 0:4
+@testset "printing of OffsetArray with n=$n" for n = 0:4
     a = OffsetArray(fill(1.,ntuple(d->1,n)), ntuple(identity,n))
     show(IOContext(io, :limit => true), MIME("text/plain"), a)
     @test String(take!(io)) == targets1[n+1]
@@ -362,9 +362,9 @@ A = OffsetArray(rand(4,4), (-3,5))
 @test maximum(A) == maximum(parent(A))
 @test minimum(A) == minimum(parent(A))
 @test extrema(A) == extrema(parent(A))
-@test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), (0,A.offsets[2]))
-@test maximum(A, dims=2) == OffsetArray(maximum(parent(A), dims=2), (A.offsets[1],0))
-@test maximum(A, dims=1:2) == maximum(parent(A), dims=1:2)
+@test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), A.offsets)
+@test maximum(A, dims=2) == OffsetArray(maximum(parent(A), dims=2), A.offsets)
+@test maximum(A, dims=1:2) == OffsetArray(maximum(parent(A), dims=1:2), A.offsets)
 C = similar(A)
 cumsum!(C, A, dims=1)
 @test parent(C) == cumsum(parent(A), dims=1)
@@ -396,11 +396,11 @@ I = findall(!iszero, z)
 @test findall(x->x==0, h) == [2]
 @test mean(A_3_3) == median(A_3_3) == 5
 @test mean(x->2x, A_3_3) == 10
-@test mean(A_3_3, dims=1) == median(A_3_3, dims=1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2]))
-@test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0))
+@test mean(A_3_3, dims=1) == median(A_3_3, dims=1) == OffsetArray([2 5 8], A_3_3.offsets)
+@test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), A_3_3.offsets)
 @test var(A_3_3) == 7.5
-@test std(A_3_3, dims=1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2]))
-@test std(A_3_3, dims=2) == OffsetArray(reshape([3,3,3], (3,1)), (A_3_3.offsets[1],0))
+@test std(A_3_3, dims=1) == OffsetArray([1 1 1], A_3_3.offsets)
+@test std(A_3_3, dims=2) == OffsetArray(reshape([3,3,3], (3,1)), A_3_3.offsets)
 @test sum(OffsetArray(fill(1,3000), -1000)) == 3000
 
 @test norm(v) ≈ norm(parent(v))
@@ -484,3 +484,8 @@ module SimilarUR
     @test_throws MethodError similar(a, Float64, (ur,))
     @test_throws MethodError similar(a, (2.0,3.0))
 end
+
+@testset "Issue 28101" begin
+    A = OffsetArray(reshape(16:-1:1, (4, 4)), (-3,5))
+    @test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), A.offsets)
+end

test/reducedim.jl

@@ -17,11 +17,10 @@ safe_sumabs2(A::Array{T}, region) where {T} = safe_mapslices(sum, abs2.(A), regi
 safe_maxabs(A::Array{T}, region) where {T} = safe_mapslices(maximum, abs.(A), region)
 safe_minabs(A::Array{T}, region) where {T} = safe_mapslices(minimum, abs.(A), region)
 
-Areduc = rand(3, 4, 5, 6)
-for region in Any[
+@testset "test reductions over region: $region" for region in Any[
     1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4),
     (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)]
-    # println("region = $region")
+    Areduc = rand(3, 4, 5, 6)
     r = fill(NaN, map(length, Base.reduced_indices(axes(Areduc), region)))
     @test sum!(r, Areduc) ≈ safe_sum(Areduc, region)
     @test prod!(r, Areduc) ≈ safe_prod(Areduc, region)
@@ -51,14 +50,14 @@ for region in Any[
     fill!(r, -1.5)
     @test minimum!(abs, r, Areduc, init=false) ≈ fill!(r2, -1.5)
 
-    @test sum(Areduc, dims=region) ≈ safe_sum(Areduc, region)
-    @test prod(Areduc, dims=region) ≈ safe_prod(Areduc, region)
-    @test maximum(Areduc, dims=region) ≈ safe_maximum(Areduc, region)
-    @test minimum(Areduc, dims=region) ≈ safe_minimum(Areduc, region)
-    @test sum(abs, Areduc, dims=region) ≈ safe_sumabs(Areduc, region)
-    @test sum(abs2, Areduc, dims=region) ≈ safe_sumabs2(Areduc, region)
-    @test maximum(abs, Areduc, dims=region) ≈ safe_maxabs(Areduc, region)
-    @test minimum(abs, Areduc, dims=region) ≈ safe_minabs(Areduc, region)
+    @test @inferred(sum(Areduc, dims=region)) ≈ safe_sum(Areduc, region)
+    @test @inferred(prod(Areduc, dims=region)) ≈ safe_prod(Areduc, region)
+    @test @inferred(maximum(Areduc, dims=region)) ≈ safe_maximum(Areduc, region)
+    @test @inferred(minimum(Areduc, dims=region)) ≈ safe_minimum(Areduc, region)
+    @test @inferred(sum(abs, Areduc, dims=region)) ≈ safe_sumabs(Areduc, region)
+    @test @inferred(sum(abs2, Areduc, dims=region)) ≈ safe_sumabs2(Areduc, region)
+    @test @inferred(maximum(abs, Areduc, dims=region)) ≈ safe_maxabs(Areduc, region)
+    @test @inferred(minimum(abs, Areduc, dims=region)) ≈ safe_minabs(Areduc, region)
 end
 
 # Test reduction along first dimension; this is special-cased for
@@ -330,10 +329,10 @@ end
 @test sum(Any[1 2;3 4], dims=1) == [4 6]
 @test sum(Vector{Int}[[1,2],[4,3]], dims=1)[1] == [5,5]
 
-# issue #10461
-Areduc = rand(3, 4, 5, 6)
-for region in Any[-1, 0, (-1, 2), [0, 1], (1,-2,3), [0 1;
+@testset "Issue #10461. region=$region" for region in Any[-1, 0, (-1, 2), [0, 1], (1,-2,3), [0 1;
                                                      2 3], "hello"]
+    Areduc = rand(3, 4, 5, 6)
+
     @test_throws ArgumentError sum(Areduc, dims=region)
     @test_throws ArgumentError prod(Areduc, dims=region)
     @test_throws ArgumentError maximum(Areduc, dims=region)