Adjust initialization in maximum and minimum

Fixes #27836

Author: andreasnoack

base/reducedim.jl

@@ -79,25 +79,6 @@ end
 reducedim_initarray(A::AbstractArray, region, v0, ::Type{R}) where {R} = fill!(similar(A,R,reduced_indices(A,region)), v0)
 reducedim_initarray(A::AbstractArray, region, v0::T) where {T} = reducedim_initarray(A, region, v0, 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, region)
-reducedim_initarray0_empty(A::AbstractArray, region,::typeof(identity), ops) = mapslices(ops, A, region)
-
 # TODO: better way to handle reducedim initialization
 #
 # The current scheme is basically following Steven G. Johnson's original implementation
@@ -124,8 +105,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, supval) in ((min, max, Inf), (max, min, -Inf))
+    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 = Base.reduced_indices(A, region)
+
+        # Next, throw if reduction is over a region with length zero
+        any(i -> iszero(size(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)(supval) : 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

@@ -1608,8 +1608,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/reducedim.jl

@@ -18,7 +18,7 @@ safe_maxabs(A::Array{T}, region) where {T} = safe_mapslices(maximum, abs.(A), re
 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")
@@ -318,10 +318,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)