From 6c54a66711f885efb932cb7c999be5ecb6452632 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Thu, 8 Oct 2026 23:52:42 -0400 Subject: [PATCH 1/2] Make sparse `map`, `broadcast` and reductions agree with dense on non-finite values, order, `missing` and failed writes Reductions: - `prod` no longer returns zero when a stored entry is `Inf` or `NaN`. - `mapreduce(f, +, A; dims=2)` does not add `f(0)` to a row that stores every entry. - `mapreduce(pred, |, A; dims=1, init)` and `&` fold into `init`. - `mapreduce` with an `op` not known to commute folds the elements in their order, for a matrix, a vector, a view of columns, along a dimension and with `sparse = true`. - `findmin` and `findmax` of a matrix accept `missing`. - `any` and `all` reject a predicate that is not Boolean. `map` and `broadcast`: - The kernels for a function that is not zero at zero compare with the fill value by identity, so `missing` and a value equal to it in another type are stored. - A function whose value has no zero (`String`, `Tuple`, `nothing`, `Pair`) is accepted. - `broadcast!` of another shape into a fixed destination computes the result aside and checks it against the pattern before writing. - `x .* y'` falls back to the general kernel when a stored entry is not finite. - A dense result is not built through an index type that cannot address every entry. - `map!` and `broadcast!` leave the destination empty, and valid, when `f` or a conversion throws part-way; densifying checks the index type before writing. - A `Diagonal` argument keeps the sign of a zero on its diagonal. Co-Authored-By: Claude Opus 5.5 --- src/higherorderfns.jl | 135 ++++++++++++++++++++------ src/reductions.jl | 211 ++++++++++++++++++++++++++++++++++------- src/sparsevector.jl | 6 +- test/higherorderfns.jl | 74 +++++++++++++++ test/reductions.jl | 54 +++++++++++ 5 files changed, 417 insertions(+), 63 deletions(-) diff --git a/src/higherorderfns.jl b/src/higherorderfns.jl index 95620629..30a8e714 100644 --- a/src/higherorderfns.jl +++ b/src/higherorderfns.jl @@ -192,13 +192,35 @@ end # Avoid calculating f(zero) unless necessary as it may fail. let r = if _haszeros(A) && all(_haszeros, Bs) fofzeros = f(_zeros_eltypes(A, Bs...)...) - _iszero(fofzeros) ? _map_zeropres!(f, C, A, Bs...) : + _iszeroresult(fofzeros) ? _guardpattern(_map_zeropres!, f, C, A, Bs...) : _map_notzeropres!(f, fofzeros, C, A, Bs...) else - _map_zeropres!(f, C, A, Bs...) + _guardpattern(_map_zeropres!, f, C, A, Bs...) end return @if_move_fixed A Bs... r end +# A zero-preserving kernel rebuilds the pattern of C as it goes. When `f` or the conversion +# of its result throws part-way, C is left structurally empty rather than half-written; a +# fixed C keeps its pattern throughout. +function _guardpattern(kernel!::K, f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {K,Tf,N} + try + return kernel!(f, C, As...) + catch + _is_fixed(C) || (trimstorage!(C, 0); _finishempty!(C)) + rethrow() + end +end +# Whether `f` at the zeros is a zero. A value of a type that has none, such as a `String`, +# a `Tuple` or `nothing`, is not one, and fills the result like any other nonzero. +_iszeroresult(x::Union{Number,AbstractArray,Missing}) = _iszero(x) +function _iszeroresult(x) + try + return _iszero(x) + catch e + e isa MethodError || rethrow() + return false + end +end function _noshapecheck_map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N} # Avoid calculating f(zero) unless necessary as it may fail. @@ -206,7 +228,7 @@ function _noshapecheck_map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N indextypeC = _promote_indtype(A, Bs...) r = if _haszeros(A) && all(_haszeros, Bs) fofzeros = f(_zeros_eltypes(A, Bs...)...) - fpreszeros = _iszero(fofzeros) + fpreszeros = _iszeroresult(fofzeros) maxnnzC = Int(fpreszeros ? min(widelength(A), _sumnnzs(A, Bs...)) : widelength(A)) C = _allocres(size(A), indextypeC, entrytypeC, maxnnzC) fpreszeros ? _map_zeropres!(f, C, A, Bs...) : @@ -227,7 +249,7 @@ copy(bc::SpBroadcasted1) = _asksthezeros(bc) && _fillszeros(bc.f, bc.args[1]) ? isempty(C) && return _finishempty!(C) f = bc.f fofnoargs = f() - if _iszero(fofnoargs) # f() is zero, so empty C + if _iszeroresult(fofnoargs) # f() is zero, so empty C trimstorage!(C, 0) _finishempty!(C) else # f() is nonzero, so densify C and fill with independent calls to f() @@ -249,7 +271,7 @@ function _diffshape_broadcast(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMa # Avoid calculating f(zero) unless necessary as it may fail. r = if _haszeros(A) && all(_haszeros, Bs) fofzeros = f(_zeros_eltypes(A, Bs...)...) - _iszero(fofzeros) ? + _iszeroresult(fofzeros) ? _broadcast_zeropres!(f, _allocres(shapeC, indextypeC, entrytypeC, maxnnzC), A, Bs...) : _broadcast_notzeropres!(f, fofzeros, _allocres(shapeC, indextypeC, entrytypeC, _densennz(shapeC)), A, Bs...) else @@ -323,12 +345,15 @@ Densifies `C`, storing `fillvalue` in place of each unstored entry in `A` and function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat) where Tf # Build dense matrix structure in C, expanding storage if necessary _densestructure!(C) + # Below, and in the other `notzeropres` kernels, a value is left as filled only when it is + # identical to `fillvalue`: one that compares equal may be of another type or sign, and + # the comparison itself may be `missing`. # Populate values fill!(storedvals(C), fillvalue) @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C)) for Ak in colrange(A, j) Cx = f(storedvals(A)[Ak]) - Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx) end end # NOTE: Combining the fill! above into the loop above to avoid multiple sweeps over / @@ -414,7 +439,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMa Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi Bk += oneunit(Bk); Bi = Bk < stopBk ? W(storedinds(B)[Bk]) : rowsentinelB end - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end return _checkbuffers(C) @@ -470,7 +495,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{Spars while activerow < rowsentinel vals, ks, rows = _fusedupdate_all(rowsentinel, activerow, rows, ks, stopks, As) Cx = f(vals...) - Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + activerow] = Cx) activerow = min(rows...) end end @@ -611,7 +636,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV bccolrangejA = numcols(A) == 1 ? colrange(A, 1) : colrange(A, j) for Ak in bccolrangejA Cx, Ci = f(storedvals(A)[Ak]), storedinds(A)[Ak] - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end # Cases with vertical expansion @@ -621,7 +646,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j)) Ax = Ak < stopAk ? svA[Ak] : zero(eltype(A)) fofAx = f(Ax) - if fofAx != fillvalue + if fofAx !== fillvalue for i in (jo + 1):(jo + numrows(C)) svC[i] = fofAx end @@ -839,7 +864,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV Ak += oneunit(Ak); Ai = Ak < stopAk ? W(storedinds(A)[Ak]) : rowsentinelA Bk += oneunit(Bk); Bi = Bk < stopBk ? W(storedinds(B)[Bk]) : rowsentinelB end - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end # Cases with vertical expansion @@ -850,7 +875,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A)) Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B)) Cx = f(Ax, Bx) - if Cx != fillvalue + if Cx !== fillvalue for Ck::Int in (jo + 1):(jo + numrows(C)) storedvals(C)[Ck] = Cx end @@ -862,10 +887,10 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j)) Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A)) fvAzB = f(Ax, zero(eltype(B))) - if fvAzB == fillvalue + if fvAzB === fillvalue while Bk < stopBk Cx = f(Ax, storedvals(B)[Bk]) - Cx != fillvalue && (storedvals(C)[jo + storedinds(B)[Bk]] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + storedinds(B)[Bk]] = Cx) Bk += oneunit(Bk) end else @@ -877,7 +902,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV else Cx = fvAzB end - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end end @@ -887,10 +912,10 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j)) Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B)) fzAvB = f(zero(eltype(A)), Bx) - if fzAvB == fillvalue + if fzAvB === fillvalue while Ak < stopAk Cx = f(storedvals(A)[Ak], Bx) - Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx) Ak += oneunit(Ak) end else @@ -902,7 +927,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV else Cx = fzAvB end - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end end @@ -952,7 +977,14 @@ _finishempty!(C::AbstractCompressedVector) = C _finishempty!(C::AbstractSparseMatrixCSC) = (fill!(getcolptr(C), 1); C) # special case - vector outer product -_copy(f::typeof(*), x::SparseVectorOrView, y::AdjOrTransSparseVectorOrView) = _outer(x, y) +_copy(f::typeof(*), x::SparseVectorOrView, y::AdjOrTransSparseVectorOrView) = + _outerkeepszeros(x, y) ? _outer(x, y) : _copy(f, _sparsifystructured(x), _sparsifystructured(y)) +# `_outer` multiplies stored entries with stored entries alone, which gives the product when +# every stored entry times a zero is a zero. An `Inf` or a `NaN` times a zero is not. +function _outerkeepszeros(x, y) + zx, zy = zero(eltype(x)), zero(eltype(y)) + return all(v -> _iszero(v * zy), nonzeros(x)) && all(v -> _iszero(zx * v), nonzeros(parent(y))) +end @inline _outer(x::SparseVectorOrView, y::Adjoint) = return _outer(conj, x, parent(y)) @inline _outer(x::SparseVectorOrView, y::Transpose) = return _outer(identity, x, parent(y)) function _outer(trans::Tf, x, y) where Tf @@ -1077,11 +1109,11 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg rows = _initrowforcol_all(j, rowsentinel, isemptys, expandsverts, ks, As) defaultCx = f(defargs...) activerow = min(rows...) - if defaultCx == fillvalue # fillvalue-preserving column scan + if defaultCx === fillvalue # fillvalue-preserving column scan while activerow < rowsentinel args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As) Cx = f(args...) - Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + activerow] = Cx) activerow = min(rows...) end else # fillvalue-non-preserving column scan @@ -1093,7 +1125,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg else Cx = defaultCx end - Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) + Cx !== fillvalue && (storedvals(C)[jo + Ci] = Cx) end end end @@ -1188,15 +1220,52 @@ end @inline function _copyto!(f, dest, As::SparseVecOrMat...) _aresameshape(dest, As...) && return _noshapecheck_map!(f, dest, As...) Base.Broadcast.check_broadcast_axes(axes(dest), As...) + # The kernels rebuild the pattern of their destination, which a fixed one does not allow: + # its result is computed aside and stored once it is known to fit. + _is_fixed(dest) && return _storeinfixed!(dest, + _diffshape_broadcast!(f, _allocres(size(dest), indtype(dest), eltype(dest), 0), As...)) + return _diffshape_broadcast!(f, dest, As...) +end +function _diffshape_broadcast!(f::Tf, dest::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N} # Avoid calculating f(zero) unless necessary as it may fail. - all(_haszeros, As) || return _broadcast_zeropres!(f, dest, As...) + all(_haszeros, As) || return _guardpattern(_broadcast_zeropres!, f, dest, As...) fofzeros = f(_zeros_eltypes(As...)...) - if _iszero(fofzeros) - return _broadcast_zeropres!(f, dest, As...) + if _iszeroresult(fofzeros) + return _guardpattern(_broadcast_zeropres!, f, dest, As...) else return _broadcast_notzeropres!(f, fofzeros, dest, As...) end end +# Every entry of `S` that is not a zero must lie in the pattern of the fixed `C`, which is +# checked in full before anything is written. The other entries of `C` become zero. +function _storeinfixed!(C::SparseVecOrMat, S::SparseVecOrMat) + for write in (false, true) + @inbounds for j in columns(C) + Ck, stopCk = colstartind(C, j), colboundind(C, j) + for Sk in colstartind(S, j):(colboundind(S, j) - 1) + Si = storedinds(S)[Sk] + while Ck < stopCk && storedinds(C)[Ck] < Si + write && (storedvals(C)[Ck] = zero(eltype(C))) + Ck += oneunit(Ck) + end + if Ck < stopCk && storedinds(C)[Ck] == Si + write && (storedvals(C)[Ck] = storedvals(S)[Sk]) + Ck += oneunit(Ck) + else + _iszero(storedvals(S)[Sk]) || _throwfixedresult(C) + end + end + if write + for k in Ck:(stopCk - oneunit(stopCk)) + storedvals(C)[k] = zero(eltype(C)) + end + end + end + end + return C +end +@noinline _throwfixedresult(C) = + throw(ArgumentError("the result has a nonzero entry outside the sparsity pattern of the $(nameof(typeof(C))) destination, whose pattern is read-only")) @inline function _copyto!(f, dest, args...) # args contains nothing but SparseVecOrMat and scalars @@ -1322,11 +1391,13 @@ _fillszeros(f, args...) = false # A function that throws at the zeros, or an element type without a zero, is left to the # sparse kernels, which evaluate the function there only if some entry is not stored. function _preserveszeros(f::F, args::Vararg{Any,N}) where {F,N} + local fofzeros try - return _iszero(f(_zeros_eltypes(args...)...)) + fofzeros = f(_zeros_eltypes(args...)...) catch return true end + return _iszeroresult(fofzeros) end # The sparse kernels evaluate the function once at the zeros and then at the stored entries @@ -1336,8 +1407,12 @@ end function _densenonzero(bc::Broadcasted, f::F, args::Vararg{Any,N}) where {F,N} T = Base.promote_typejoin_union(Base.promote_op(f, map(eltype, args)...)) (isconcretetype(T) && T !== Bool) || return _densebroadcast(bc) - return _storedasdense(_copy(f, map(_sparsifystructured, args)...)) + As = map(_sparsifystructured, args) + # the kernels store every entry, which the index type of the arguments may not address + _addresses(_promote_indtype(As...), map(length, axes(bc))) || return _densebroadcast(bc) + return _storedasdense(_copy(f, As...)) end +_addresses(::Type{Ti}, shape) where {Ti} = !(Ti <: Base.BitInteger) || prod(Int64, shape) < typemax(Ti) _storedasdense(C::Union{SparseVector,SparseMatrixCSC}) = length(getnzval(C)) == nnz(C) == length(C) ? reshape(getnzval(C), size(C)) : Array(C) _storedasdense(C) = Array(C) @@ -1374,6 +1449,12 @@ _densifysparse(x) = x _sparsifystructured(M::AbstractMatrix) = _isdenselike(M) ? _fullystored(M) : SparseMatrixCSC(M) _sparsifystructured(V::AbstractVector) = _isdenselike(V) ? _fullystored(V) : SparseVector(V) _sparsifystructured(S::SparseVecOrMat) = S +# stores the whole diagonal, as for the other banded types: a kernel that read a `-0.0` there +# as a structural zero would lose its sign +function _sparsifystructured(D::Diagonal) + n = size(D, 1) + return SparseMatrixCSC(n, n, collect(1:(n + 1)), collect(1:n), copyto!(Vector{eltype(D)}(undef, n), D.diag)) +end _sparsifystructured(S::SparseViewOfColumns) = copy(S) _sparsifystructured(x) = x diff --git a/src/reductions.jl b/src/reductions.jl index 3930d243..ad60d9b5 100644 --- a/src/reductions.jl +++ b/src/reductions.jl @@ -61,17 +61,84 @@ end function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} z = nnz(A) n = widelength(A) - if z == 0 - if n == 0 - Base.mapreduce_empty(f, op, T) - else - _mapreducezeros(f, op, T, n-z-1, f(zero(T))) - end + if n == 0 + Base.mapreduce_empty(f, op, T) + elseif !_commutes(f, op, T) + # every entry stored: they are folded in their order, without `f(0)`, which might throw + z == n ? Base._mapreduce(f, op, nzvalview(A)) : _mapreduce_ordered(f, op, A) + elseif z == 0 + _mapreducezeros(f, op, T, n-z-1, f(zero(T))) else _mapreducezeros(f, op, T, n-z, Base._mapreduce(f, op, nzvalview(A))) end end +# Whether `op` may fold the mapped values out of their order. It is taken to be associative, +# and to be commutative only where it is known to be: a product commutes when its factors are +# real or complex numbers, which is asked of inference. +_commutes(f, ::Union{typeof(+),typeof(Base.add_sum),typeof(max),typeof(min),typeof(|),typeof(&), + typeof(Base.or_any),typeof(Base.and_all),typeof(Base._extrema_rf)}, ::Type) = true +_commutes(f::F, ::Union{typeof(*),typeof(Base.mul_prod)}, ::Type{T}) where {F,T} = + Base.promote_op(f, T) <: Union{Real,Complex} +_commutes(f, op, ::Type) = false + +# An `op` not known to commute is folded over the positions in their order. `A` has an +# unstored entry, so `f(0)` is taken once. A run of unstored entries is folded into the stored +# entry after it, which an associative `op` allows and which keeps the accumulator one call +# of `op` away from its next value. +function _mapreduce_ordered(f::F, op::G, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} + rf = Base.BottomRF(op) + rowval, nzval = getrowval(A), getnzval(A) + m = size(A, 1) + z, idem = _unstoredvalue(f, op, T) + v = Base._InitialValue() + unstored = Int64(0) # unstored positions since the last stored one + @inbounds for col in axes(A, 2) + prev = 0 + for k in getnzrange(A, col) + row = Int(rowval[k]) + v = rf(v, _afterunstored(op, z, idem, unstored + (row - prev - 1), f(nzval[k]))) + unstored, prev = Int64(0), row + end + unstored += m - prev + end + return unstored == 0 ? v : rf(v, _foldunstored(op, z, idem, unstored)) +end +function _mapreduce_ordered(f::F, op::G, A::SparseVectorOrView{T}) where {F,G,T} + rf = Base.BottomRF(op) + nzind, nzval = nonzeroinds(A), nonzeros(A) + z, idem = _unstoredvalue(f, op, T) + v = Base._InitialValue() + prev = 0 + @inbounds for k in eachindex(nzind, nzval) + i = Int(nzind[k]) + v = rf(v, _afterunstored(op, z, idem, i - prev - 1, f(nzval[k]))) + prev = i + end + n = length(A) - prev + return n == 0 ? v : rf(v, _foldunstored(op, z, idem, n)) +end +# `f(0)` as `op` first sees it, and whether folding it with itself leaves it as it is, which +# makes a run of any length that one value +function _unstoredvalue(f::F, op::G, ::Type{T}) where {F,G,T} + z = Base.reduce_first(op, f(zero(T))) + return z, isequal(op(z, z), z) +end +# `z` folded over `n > 0` consecutive positions, in O(log(n)) calls of `op` +@inline function _foldunstored(op::G, z, idem::Bool, n) where {G} + (idem || n == 1) && return z + p = r = z + n -= 1 + while n > 0 + isodd(n) && (r = op(r, p)) + n >>= 1 + n > 0 && (p = op(p, p)) + end + return r +end +# the stored value `x` with the `n` unstored positions before it folded in +@inline _afterunstored(op::G, z, idem::Bool, n, x) where {G} = n == 0 ? x : op(_foldunstored(op, z, idem, n), x) + # Specialized mapreduce for +/*/min/max/_extrema_rf. The entries a sparse array does not # store are folded in at once, as `reduce_first` of `f(0)` (an `Int` for `sum` and `prod` of # small integers, as for dense input) times or to the power of their count, in the type the @@ -93,10 +160,11 @@ _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nze nzeros == 0 ? v0 : op(v0, f(zero(T))) # Specialized mapreduce for any and all +# The assertion rejects a predicate that is not Boolean, as the dense reductions do. Base._any(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) = - iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A) + iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)::Union{Bool,Missing} Base._all(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) = - iszero(widelength(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A) + iszero(widelength(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A)::Union{Bool,Missing} function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,T} nnzA = nnz(A) @@ -104,10 +172,13 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas if nzeros == 0 # No zeros, so don't compute f(0) since it might throw Base._mapreduce(f, op, nzvalview(A)) + elseif !_commutes(f, op, T) + _mapreduce_ordered(f, op, A) else v = Base.reduce_first(op, f(zero(T)))^nzeros - # Bail out early if initial reduction value is zero or if there are no stored elements - (_iszero(v) || nnzA == 0) ? v : v*Base._mapreduce(f, op, nzvalview(A)) + # A zero ends the product only among integers: a stored `Inf` or `NaN` makes it `NaN` + (nnzA == 0 || (_iszero(v) && Base.promote_op(f, T) <: Integer)) ? v : + v*Base._mapreduce(f, op, nzvalview(A)) end end @@ -339,6 +410,8 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,T getcolptr(R)[2] = 2 end return R + elseif (rm == 1 || rn == 1) && !_inorder(f, op, A) + return _mapreduce_dim_sparse_ordered(f, op, init, A, dims, Tr, rm == 1) elseif rm == 1 return _mapreducerows_sparse!(f, op, init, spzeros(Tr, Ti, 1, n), A) elseif rn == 1 @@ -349,6 +422,37 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,T end end +# The two kernels below fold each slice's stored entries first. For an `op` not known to +# commute, the dense reduction, which keeps the order of the elements and takes time +# proportional to nnz(A) + length(result) as well, is given the pattern they produce. +function _mapreduce_dim_sparse_ordered(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,Ti}, dims, ::Type{Tr}, + alongrows::Bool) where {T,Ti,Tr} + m, n = size(A) + D = Base._mapreduce_dim(f, op, init, A, dims) + slicelen, nslices = alongrows ? (m, n) : (n, m) + z = zero(Tr) + zunstored = nnz(A) == m*n && slicelen > 0 ? z : _reduce_unstored(f, op, init, T, slicelen) + stores = fill(!isequal(zunstored, z), nslices) + if alongrows + @inbounds for col in 1:n + stores[col] |= !isempty(getnzrange(A, col)) + end + else + @inbounds for i in view(getrowval(A), _storedinds(A)) + stores[i] = true + end + end + inds = convert(Vector{Ti}, findall(stores)) + vals = Tr[D[i] for i in inds] + alongrows || return SparseMatrixCSC(m, 1, Ti[1, length(inds) + 1], inds, vals) + colptr = Vector{Ti}(undef, n + 1) + colptr[1] = 1 + @inbounds for col in 1:n + colptr[col + 1] = colptr[col] + stores[col] + end + return SparseMatrixCSC(1, n, colptr, fill!(inds, 1), vals) +end + # `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCOrColumnSubset{T}) where T nzval = getnzval(A) @@ -475,6 +579,44 @@ function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCOrColumnSubs R end +# The kernels above fold each slice's stored entries first and its unstored ones last, which +# is the order of the elements when `op` commutes or every entry is stored. Otherwise the +# unstored entries are folded where they lie, with `f(0)` taken once. +_inorder(f, op, A::AbstractArray{T}) where {T} = _commutes(f, op, T) || nnz(A) == length(A) +function _mapreducerows_ordered!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} + require_one_based_indexing(A, R) + rowval, nzval = getrowval(A), getnzval(A) + m = size(A, 1) + z, idem = _unstoredvalue(f, op, T) + @inbounds for col in axes(A, 2) + r = R[1, col] + prev = 0 + for j in getnzrange(A, col) + row = Int(rowval[j]) + r = op(r, _afterunstored(op, z, idem, row - prev - 1, f(nzval[j]))) + prev = row + end + R[1, col] = prev == m ? r : op(r, _foldunstored(op, z, idem, m - prev)) + end + R +end +function _mapreducecols_ordered!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} + require_one_based_indexing(A, R) + rowval, nzval = getrowval(A), getnzval(A) + m, n = size(A) + z, idem = _unstoredvalue(f, op, T) + lastcol = zeros(Int, m) # the last column folded into each row + @inbounds for col in axes(A, 2), j in getnzrange(A, col) + row = rowval[j] + R[row, 1] = op(R[row, 1], _afterunstored(op, z, idem, col - lastcol[row] - 1, f(nzval[j]))) + lastcol[row] = col + end + @inbounds for i in 1:m + lastcol[i] == n || (R[i, 1] = op(R[i, 1], _foldunstored(op, z, idem, n - lastcol[i]))) + end + R +end + function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} require_one_based_indexing(A, R) lsiz = Base.check_reducedims(R,A) @@ -485,10 +627,10 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOr R[1, 1] = op(R[1, 1], mapreduce(f, op, A)) elseif size(R, 1) == 1 # Reduction along rows - _mapreducerows!(f, op, R, A) + _inorder(f, op, A) ? _mapreducerows!(f, op, R, A) : _mapreducerows_ordered!(f, op, R, A) elseif size(R, 2) == 1 # Reduction along columns - _mapreducecols!(f, op, R, A) + _inorder(f, op, A) ? _mapreducecols!(f, op, R, A) : _mapreducecols_ordered!(f, op, R, A) else # Reduction along a dimension > 2 # Compute op(R, f(A)) @@ -582,7 +724,8 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCO end end for i = 1:m - R[i, 1] += _addzeros(zeroval, rownz[i]) + # a row stored in full has no `f(0)` to add, which need not be finite + rownz[i] == 0 || (R[i, 1] += _addzeros(zeroval, rownz[i])) end end end @@ -590,37 +733,35 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCO end # any(pred, A, dims = 1) => mapreduce(pred, |, A, dims = 1) -# final argument `post` is to allow post-mapping each columnar mapreduce +# `R` holds the initial value of each column, which the column's reduction is folded into. +# With `negate`, `pred` is the negated predicate of an `&` reduction, see below. function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::SparseMatrixCSCOrColumnSubset{Tv}, - post::F = identity) where {P, F, Tv} + negate::Bool = false) where {P, Tv} nzval = getnzval(A) m, n = size(A) @inbounds for ii in axes(A,2) + # an initial value that already decides the column + R[1, ii] == negate || continue rng = getnzrange(A, ii) len = length(rng) - # An empty column is trivial - if len == 0 - R[1, ii] = post(pred(zero(Tv))) - continue - end - # If predicate on zero is true, then sparse column can be short-circuited - if pred(zero(Tv)) && len < m - R[1, ii] = post(true) - continue - end - # Otherwise reduce over the stored values - r = false - for jj in rng - r = pred(nzval[jj]) - r && break + if len < m && pred(zero(Tv)) + # an unstored entry satisfies the predicate + r = true + else + # reduce over the stored values + r = false + for jj in rng + r = pred(nzval[jj]) + r && break + end end - R[1, ii] = post(r) + R[1, ii] = r != negate end return R end # all(pred, A, dims = 1) => mapreduce(pred, &, A, dims = 1) == .!mapreduce(!pred, |, A, dims = 1) _mapreducerows!(pred::P, ::typeof(&), R::AbstractMatrix{Bool}, - A::SparseMatrixCSCOrColumnSubset) where {P} = _mapreducerows!(!pred, |, R, A, !) + A::SparseMatrixCSCOrColumnSubset) where {P} = _mapreducerows!(!pred, |, R, A, true) # findmax/min and argmax/min methods # find the first entry of the sparse matrix that is not stored - return its cartesian index @@ -724,8 +865,10 @@ function _findr(op, A::AbstractSparseMatrixCSC{Tv}, region) where {Tv} end end -_isless_fm(a, b) = b == b && ( a != a || isless(a, b) ) -_isgreater_fm(a, b) = b == b && ( a != a || isless(b, a) ) +# whether `a` replaces the best value so far `b`, by the test of Base's dense `findmin` and +# `findmax`, so that the result agrees with dense on `NaN`, `missing` and signed zeros +_isless_fm(a, b) = Base.isgreater(b, a) +_isgreater_fm(a, b) = isless(b, a) findmin(A::AbstractSparseMatrixCSC{Tv}, region::Union{Integer,Tuple{Integer},NTuple{2,Integer}}) where {Tv} = _findr(_isless_fm, A, region) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index dd87cd92..125d9680 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1463,6 +1463,7 @@ function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) w T = eltype(A) isempty(A) && return Base.mapreduce_empty(f, op, T) z = nnz(A) + (z == length(A) || _commutes(f, op, T)) || return _mapreduce_ordered(f, op, A) rest, ini = if z == 0 length(A)-z-1, f(zero(T)) else @@ -1472,9 +1473,9 @@ function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) w end Base._any(f, A::SparseVectorOrView, ::Colon) = - iszero(length(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A) + iszero(length(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)::Union{Bool,Missing} Base._all(f, A::SparseVectorOrView, ::Colon) = - iszero(length(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A) + iszero(length(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A)::Union{Bool,Missing} function Base.mapreducedim!(f::F, op::G, R::AbstractVector, A::SparseVectorOrView) where {F,G} # dim1 reduction could be safely replaced with a mapreduce @@ -1883,6 +1884,7 @@ function _densestructure!(A::AbstractSparseMatrixCSC) _checkdensifiable(A) _is_fixed(A) && return A m, n = size(A) + convert(indtype(A), m * n + 1) # the last column pointer: throws before anything is written resize!(getrowval(A), m * n) resize!(getnzval(A), m * n) colptr = resize!(getcolptr(A), n + 1) diff --git a/test/higherorderfns.jl b/test/higherorderfns.jl index ebbbf695..7066e7d5 100644 --- a/test/higherorderfns.jl +++ b/test/higherorderfns.jl @@ -1225,4 +1225,78 @@ end end end +@static if COMPREHENSIVE +@testset "map and broadcast of a function whose values are `missing`, have no zero or equal f(0) in another type" begin + A, Ad = fixturepair(Float64, 4, 3) + f = x -> x == 0 ? 1.0 : missing + @test isequal(map(f, A), map(f, Ad)) + @test isequal(A .== missing, Ad .== missing) + g = (x, y) -> x == 0 ? 1.0 : missing + @test isequal(map(g, A, A), map(g, Ad, Ad)) + @test isequal(broadcast(g, A, A[:, 1]), broadcast(g, Ad, Ad[:, 1])) + g3 = (x, y, z) -> x == 0 ? 1.0 : missing + @test isequal(broadcast(g3, A, A, A[:, 1]), broadcast(g3, Ad, Ad, Ad[:, 1])) + h = x -> x == 0 ? 1 : 1.0 + @test typeof.(map(h, A)) == typeof.(map(h, Ad)) + for k in (string, x -> (x,), x -> nothing, x -> x => 1) + @test map(k, A) == map(k, Ad) + end + @test string.(A) == string.(Ad) && string.(A) isa Matrix{String} + @test (A .=> 1) == (Ad .=> 1) +end + +@testset "broadcast! of another shape into a fixed destination checks the result first" begin + F = SparseArrays.fixed(sparse([1.0 0; 1.0 2.0])) + Fd, nz = Array(F), copy(nonzeros(F)) + @test_throws ArgumentError F .= F .+ sparsevec([3.0, 0.0]) # a nonzero at (1, 2) + @test nonzeros(F) == nz + @test_throws ArgumentError F .= F .+ [1.0 2.0] + @test nonzeros(F) == nz + Fi = SparseArrays.fixed(sparse([1 0; 1 2])) + @test_throws InexactError Fi .= Fi .* [0.5, 1.0] + @test nonzeros(Fi) == [1, 1, 2] + # a result that fits the pattern is stored, whatever zeros the kernel kept + @test (F .= F .* sparsevec([3.0, 0.0])) === F + @test F == Fd .* [3.0, 0.0] && rowvals(F) == [1, 2, 2] + copyto!(nonzeros(F), nz) + F .= F .* [3.0, 2.0] + @test F == Fd .* [3.0, 2.0] +end + +@testset "broadcast! and map! leave a valid destination when `f` or a conversion throws" begin + A, _ = fixturepair(Float64, 4, 3) + top = maximum(nonzeros(A)) # not a zero: `f` is evaluated at the zeros before anything is written + @test top > 0 + boom = (x, y) -> x == top ? error("boom") : x + y + for B in (copy(A), A[:, 1]) # the `map!` kernel, and the kernel for another shape + C = copy(A) + @test_throws ErrorException C .= boom.(A, B) + @test nnz(C) == 0 && iszero(Array(C)) + C = copy(A) + @test_throws ErrorException C .= boom.(C, B) + @test nnz(C) == 0 && iszero(Array(C)) + end + Ci = sparse([1 0; 0 2]) + @test_throws InexactError Ci .= Ci .* sparse([0.5 1.0; 1.0 1.0]) + @test nnz(Ci) == 0 && iszero(Array(Ci)) + # every entry does not fit the index type: found before anything is written + S = sparse(Int8[1], Int8[1], [1.0], 12, 12) + S0 = copy(S) + @test_throws InexactError S .= 1.0 + @test_throws InexactError map!(x -> x + 1, S, S0) + @test mismatch(S, Array(S0); Ti=Int8) === nothing + # a dense result does not go through that index type + @test cos.(S) == cos.(Array(S)) && cos.(S) isa Matrix{Float64} +end + +@testset "broadcast keeps a non-finite value times an unstored zero, and the sign of a zero on a `Diagonal`" begin + for (x, y) in ((sparsevec([Inf, 0.0, 2.0]), sparsevec([0.0, 1.0, 3.0])), + (sparsevec([1.0, 0.0, 2.0]), sparsevec([0.0, NaN, 3.0])), + (sparsevec([Inf + 0im, 0, 2im]), sparsevec([0, 1im, 3.0 + 0im]))), w in (adjoint, transpose) + @test isequal(Array(x .* w(y)), Array(x) .* w(Array(y))) + end + @test isequal(Array(sparse(ones(2, 2)) ./ Diagonal([-0.0, 1.0])), ones(2, 2) ./ Diagonal([-0.0, 1.0])) +end +end + end # module diff --git a/test/reductions.jl b/test/reductions.jl index 24f310da..66ce7d8f 100644 --- a/test/reductions.jl +++ b/test/reductions.jl @@ -737,4 +737,58 @@ end end end +@static if COMPREHENSIVE +@testset "reductions agree with dense on non-finite values, `init`, the order of the elements and `missing`" begin + # a stored `Inf` or `NaN` is not absorbed by an unstored zero + for D in ([0.0 Inf; 1.0 2.0], [0.0 NaN; 1.0 2.0]), dims in (:, 1, 2, (1, 2)) + @test isequal(prod(sparse(D); dims), prod(D; dims)) + end + @test prod(sparse([0 3; 1 2])) === 0 + # `f(0)` is not added to a row that stores every entry + f = x -> x == 0 ? Inf : x + D = [1.0 2.0; 0.0 3.0] + @test isequal(mapreduce(f, +, sparse(D); dims=2), mapreduce(f, +, D; dims=2)) + # `|` and `&` along the rows fold into `init` + for (pred, op) in ((>(5), |), (!iszero, &)), init in (true, false), + D in ([1.0 2.0; 0.0 3.0], [1.0 2.0; 3.0 4.0], [0.0 6.0; 0.0 0.0]) + @test mapreduce(pred, op, sparse(D); dims=1, init) == mapreduce(pred, op, D; dims=1, init) + end + # an `op` that is not known to commute folds the elements in their order + takelast = (a, b) -> b + A, Ad = fixturepair(Float64, 4, 3) + x, xd = fixturepair(Float64, 5) + B = sparse([1.0 2.0; 3.0 4.0]) # stores every entry + for (S, D) in ((A, Ad), (x, xd), (view(A, :, [3, 1]), Ad[:, [3, 1]]), (view(A, :, 2), Ad[:, 2]), (B, Array(B))) + @test mapreduce(v -> [v], vcat, S) == mapreduce(v -> [v], vcat, D) + @test mapreduce(string, *, S) == mapreduce(string, *, D) + @test reduce(takelast, S) === reduce(takelast, D) + S isa AbstractMatrix || continue + for dims in (1, 2, (1, 2)) + rd = mapreduce(identity, takelast, D; dims, init=7.0) + @test mapreduce(identity, takelast, S; dims, init=7.0) == rd + @test mismatch(mapreduce(identity, takelast, S; dims, init=7.0, sparse=true), rd) === nothing + @test mapreduce(string, *, S; dims, init="") == mapreduce(string, *, D; dims, init="") + end + end + # the ordered fold takes `f(0)` once and each stored entry once + calls = Ref(0) + mapreduce(v -> (calls[] += 1; v), takelast, A) + @test calls[] == nnz(A) + 1 + # `missing` in `findmin` and `findmax`, which dense ranks above `NaN` + for D in (Union{Missing,Float64}[1.0 missing; 0.0 2.0], Union{Missing,Float64}[NaN missing; 0.0 2.0], + Union{Missing,Float64}[0.0 NaN; -1.0 0.0]), fm in (findmin, findmax) + S = sparse(D) + @test isequal(fm(S), fm(D)) + @test isequal(fm(S; dims=1), fm(D; dims=1)) && isequal(fm(S; dims=2), fm(D; dims=2)) + end + # a predicate that is not Boolean is rejected, as for dense + for S in (sparse([1 0; 0 0]), sparsevec([1, 0])) + @test_throws TypeError any(S) + @test_throws TypeError all(S) + end + M = sparse(Union{Missing,Bool}[missing false; false false]) + @test any(M) === missing && all(M) === false +end +end + end # module From ab82a141d4e4f2950db8d8bbe1a0f69d42c5b068 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Fri, 9 Oct 2026 10:12:38 -0400 Subject: [PATCH 2/2] Check every predicate value in sparse `any` and `all` Asserting the type of the result let `missing | 1` through as `missing`. Co-Authored-By: Claude Opus 5.5 --- src/reductions.jl | 10 +++++++--- src/sparsevector.jl | 4 ++-- test/reductions.jl | 10 ++++++++++ 3 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/reductions.jl b/src/reductions.jl index ad60d9b5..bfabbd43 100644 --- a/src/reductions.jl +++ b/src/reductions.jl @@ -160,11 +160,15 @@ _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nze nzeros == 0 ? v0 : op(v0, f(zero(T))) # Specialized mapreduce for any and all -# The assertion rejects a predicate that is not Boolean, as the dense reductions do. +# Each value of the predicate is checked to be Boolean, as the dense reductions check it: +# `|` and `&` would absorb another value next to a `missing`, or combine two of them. +_boolean(v::Union{Bool,Missing}) = v +_boolean(v) = v::Bool +_booleanpredicate(f::F) where {F} = _boolean ∘ f Base._any(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) = - iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)::Union{Bool,Missing} + iszero(widelength(A)) ? false : Base._mapreduce(_booleanpredicate(f), |, IndexCartesian(), A) Base._all(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) = - iszero(widelength(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A)::Union{Bool,Missing} + iszero(widelength(A)) ? true : Base._mapreduce(_booleanpredicate(f), &, IndexCartesian(), A) function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,T} nnzA = nnz(A) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 125d9680..cabd9910 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1473,9 +1473,9 @@ function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) w end Base._any(f, A::SparseVectorOrView, ::Colon) = - iszero(length(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)::Union{Bool,Missing} + iszero(length(A)) ? false : Base._mapreduce(_booleanpredicate(f), |, IndexCartesian(), A) Base._all(f, A::SparseVectorOrView, ::Colon) = - iszero(length(A)) ? true : Base._mapreduce(f, &, IndexCartesian(), A)::Union{Bool,Missing} + iszero(length(A)) ? true : Base._mapreduce(_booleanpredicate(f), &, IndexCartesian(), A) function Base.mapreducedim!(f::F, op::G, R::AbstractVector, A::SparseVectorOrView) where {F,G} # dim1 reduction could be safely replaced with a mapreduce diff --git a/test/reductions.jl b/test/reductions.jl index 66ce7d8f..4daf310e 100644 --- a/test/reductions.jl +++ b/test/reductions.jl @@ -786,8 +786,18 @@ end @test_throws TypeError any(S) @test_throws TypeError all(S) end + # also next to a `missing`, which `|` and `&` would let absorb the other value + for S in (sparse(reshape(Union{Missing,Int}[missing, 1], 2, 1)), sparse(Union{Missing,Int}[missing, 1]), + sparse(Union{Missing,Int}[missing 0; 1 2])) + @test_throws TypeError any(S) + @test_throws TypeError all(S) + @test_throws TypeError any(Array(S)) + @test_throws TypeError all(Array(S)) + end M = sparse(Union{Missing,Bool}[missing false; false false]) @test any(M) === missing && all(M) === false + @test any(iszero, sparse([1.0 0; 2 3])) === true && all(iszero, sparse([1.0 0; 2 3])) === false + @test any(ismissing, M) === true && all(x -> x === false, sparsevec([false, false])) === true end end