Adjust initialization in maximum and minimum (version 2) (JuliaLang#2…

…8116)

* Adjust initialization in maximum and minimum

Fixes JuliaLang#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)