From 2dec43617a66f80240a04fec1f5145122c3306ce Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 10:01:35 -0400 Subject: [PATCH 01/14] Add `map` interface for carefully checking zeros --- src/SparseArraysBase.jl | 1 + src/map.jl | 122 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 src/map.jl diff --git a/src/SparseArraysBase.jl b/src/SparseArraysBase.jl index 0e0db90..c51b99b 100644 --- a/src/SparseArraysBase.jl +++ b/src/SparseArraysBase.jl @@ -20,6 +20,7 @@ export SparseArrayDOK, include("abstractsparsearrayinterface.jl") include("sparsearrayinterface.jl") include("indexing.jl") +include("map.jl") include("wrappers.jl") include("abstractsparsearray.jl") include("sparsearraydok.jl") diff --git a/src/map.jl b/src/map.jl new file mode 100644 index 0000000..2e12439 --- /dev/null +++ b/src/map.jl @@ -0,0 +1,122 @@ +# zero-preserving Traits +# ---------------------- +""" + abstract type ZeroPreserving end + +Holy Trait to indicate how a function interacts with abstract zero values: + +- `StrongPreserving` : output is guaranteed to be zero if **any** input is. +- `WeakPreserving` : output is guaranteed to be zero if **all** inputs are. +- `NonPreserving` : no guarantees on output. + +To attempt to automatically determine this, either `ZeroPreserving(f, A::AbstractArray...)` or +`ZeroPreserving(f, T::Type...)` can be used/overloaded. + +!!! warning + incorrectly registering a function to be zero-preserving will lead to silently wrong results. +""" +abstract type ZeroPreserving end +struct StrongPreserving <: ZeroPreserving end +struct WeakPreserving <: ZeroPreserving end +struct NonPreserving <: ZeroPreserving end + +# warning: cannot automatically detect WeakPreserving since this would mean checking all values +function ZeroPreserving(f, A::AbstractArray, Bs::AbstractArray...) + return ZeroPreserving(f, eltype(A), eltype.(Bs)...) +end +# TODO: the following might not properly specialize on the types +# TODO: non-concrete element types +function ZeroPreserving(f, T::Type, Ts::Type...) + return iszero(f(zero(T), zero.(Ts)...)) ? WeakPreserving() : NonPreserving() +end + +const _WEAK_FUNCTIONS = (:+, :-) +for f in _WEAK_FUNCTIONS + @eval begin + ZeroPreserving(::typeof($f), ::AbstractArray, ::AbstractArray...) = WeakPreserving() + ZeroPreserving(::typeof($f), ::Type, ::Type...) = WeakPreserving() + end +end + +const _STRONG_FUNCTIONS = (:*,) +for f in _STRONG_FUNCTIONS + @eval begin + ZeroPreserving(::typeof($f), ::AbstractArray, ::AbstractArray...) = StrongPreserving() + ZeroPreserving(::typeof($f), ::Type, ::Type...) = StrongPreserving() + end +end + +# map(!) +# ------ +@interface I::AbstractSparseArrayInterface function Base.map( + f, A::AbstractArray, Bs::AbstractArray... +) + T = Base.Broadcast.combine_eltypes(f, (A, Bs...)) + C = similar(I, T, size(A)) + return @interface I map!(f, C, A, Bs...) +end + +@interface ::AbstractSparseArrayInterface function Base.map!( + f, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... +) + return _map!(f, ZeroPreserving(f, A, Bs...), C, A, Bs...) +end + +function _map!( + f, ::StrongPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... +) + checkshape(C, A, Bs...) + style = IndexStyle(C, A, Bs...) + unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) + zero!(C) + for I in intersect(eachstoredindex.(Ref(style), unaliased)...) + @inbounds C[I] = f(ith_all(I, unaliased)...) + end + return C +end +function _map!( + f, ::WeakPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... +) + checkshape(C, A, Bs...) + style = IndexStyle(C, A, Bs...) + unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) + zero!(C) + for I in union(eachstoredindex.(Ref(style), unaliased)...) + @inbounds C[I] = f(ith_all(I, unaliased)...) + end + return C +end +function _map!(f, ::NonPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray...) + checkshape(C, A, Bs...) + unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) + for I in eachindex(C, A, Bs...) + @inbounds C[I] = f(ith_all(I, unaliased)...) + end + return C +end + +# Derived functions +# ----------------- +@interface I::AbstractSparseArrayInterface Base.copyto!( + C::AbstractArray, A::AbstractArray +) = @interface I map!(identity, C, A) + +# Utility functions +# ----------------- +# shape check similar to checkbounds +checkshape(::Type{Bool}, A::AbstractArray) = true +checkshape(::Type{Bool}, A::AbstractArray, B::AbstractArray) = size(A) == size(B) +function checkshape(::Type{Bool}, A::AbstractArray, Bs::AbstractArray...) + return allequal(size, (A, Bs...)) +end + +function checkshape(A::AbstractArray, Bs::AbstractArray...) + return checkshape(Bool, A, Bs...) || + throw(DimensionMismatch("argument shapes must match")) +end + +@inline ith_all(i, ::Tuple{}) = () +function ith_all(i, as) + @_propagate_inbounds_meta + return (as[1][i], ith_all(i, Base.tail(as))...) +end From b4b2328fdac1fc300aa191c145b56e4dcfd8165b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 10:40:49 -0400 Subject: [PATCH 02/14] disable previous implementation --- src/abstractsparsearrayinterface.jl | 32 ++++++++++++++--------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/abstractsparsearrayinterface.jl b/src/abstractsparsearrayinterface.jl index 6af80e6..b502b9a 100644 --- a/src/abstractsparsearrayinterface.jl +++ b/src/abstractsparsearrayinterface.jl @@ -153,22 +153,22 @@ function preserves_unstored(f, a_dest::AbstractArray, as::AbstractArray...) return iszero(f(map(a -> getunstoredindex(a, I), as)...)) end -@interface interface::AbstractSparseArrayInterface function Base.map!( - f, a_dest::AbstractArray, as::AbstractArray... -) - isempty(a_dest) && return a_dest # special case to avoid trying to access empty array - indices = if !preserves_unstored(f, a_dest, as...) - eachindex(a_dest) - elseif any(a -> a_dest !== a, as) - as = map(a -> Base.unalias(a_dest, a), as) - @interface interface zero!(a_dest) - eachstoredindex(as...) - else - eachstoredindex(a_dest) - end - @interface interface map_indices!(indices, f, a_dest, as...) - return a_dest -end +# @interface interface::AbstractSparseArrayInterface function Base.map!( +# f, a_dest::AbstractArray, as::AbstractArray... +# ) +# isempty(a_dest) && return a_dest # special case to avoid trying to access empty array +# indices = if !preserves_unstored(f, a_dest, as...) +# eachindex(a_dest) +# elseif any(a -> a_dest !== a, as) +# as = map(a -> Base.unalias(a_dest, a), as) +# @interface interface zero!(a_dest) +# eachstoredindex(as...) +# else +# eachstoredindex(a_dest) +# end +# @interface interface map_indices!(indices, f, a_dest, as...) +# return a_dest +# end # `f::typeof(norm)`, `op::typeof(max)` used by `norm`. function reduce_init(f, op, as...) From bd917909f4940074225b8f962acf12ab83daaae8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 10:48:32 -0400 Subject: [PATCH 03/14] Update indexstyle for eachstoredindex --- src/indexing.jl | 16 ++++++++++++---- src/oneelementarray.jl | 2 +- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/indexing.jl b/src/indexing.jl index 3cd089d..cbd13f7 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -308,10 +308,18 @@ end end end -# required: -@interface ::AbstractSparseArrayInterface eachstoredindex(style::IndexStyle, A::AbstractArray) = throw( - MethodError(eachstoredindex, Tuple{typeof(style),typeof(A)}) -) +# required: one implementation for canonical index style +@interface ::AbstractSparseArrayInterface function eachstoredindex(style::IndexStyle, A::AbstractArray) + if style == IndexStyle(A) + throw(MethodError(eachstoredindex, Tuple{typeof(style),typeof(A)})) + elseif style == IndexCartesian() + return map(Base.Fix1(Base.getindex, CartesianIndices(A)), eachindex(A)) + elseif style == IndexLinear() + return map(Base.Fix1(Base.getindex, LinearIndices(A)), eachindex(A)) + else + throw(ArgumentError(lazy"unknown index style $style")) + end +end # derived but may be specialized: @interface ::AbstractSparseArrayInterface function eachstoredindex( diff --git a/src/oneelementarray.jl b/src/oneelementarray.jl index 7936c6c..251acfa 100644 --- a/src/oneelementarray.jl +++ b/src/oneelementarray.jl @@ -287,7 +287,7 @@ storedindex(a::OneElementArray) = getfield(a, :index) function isstored(a::OneElementArray, I::Int...) return I == storedindex(a) end -function eachstoredindex(a::OneElementArray) +function eachstoredindex(::IndexCartesian, a::OneElementArray) return Fill(CartesianIndex(storedindex(a)), 1) end From 0966414de1cf50135d986674d0b4a87d84996d06 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 14:21:21 -0400 Subject: [PATCH 04/14] `haszero` --- src/map.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/map.jl b/src/map.jl index 2e12439..859d274 100644 --- a/src/map.jl +++ b/src/map.jl @@ -20,6 +20,14 @@ struct StrongPreserving <: ZeroPreserving end struct WeakPreserving <: ZeroPreserving end struct NonPreserving <: ZeroPreserving end +# Backport: remove in 1.12 +@static if !isdefined(Base, :haszero) + _haszero(T::Type) = false + _haszero(::Type{<:Number}) = true +else + _haszero = Base.haszero +end + # warning: cannot automatically detect WeakPreserving since this would mean checking all values function ZeroPreserving(f, A::AbstractArray, Bs::AbstractArray...) return ZeroPreserving(f, eltype(A), eltype.(Bs)...) @@ -27,7 +35,11 @@ end # TODO: the following might not properly specialize on the types # TODO: non-concrete element types function ZeroPreserving(f, T::Type, Ts::Type...) - return iszero(f(zero(T), zero.(Ts)...)) ? WeakPreserving() : NonPreserving() + if all(_haszero, (T, Ts...)) + return iszero(f(zero(T), zero.(Ts)...)) ? WeakPreserving() : NonPreserving() + else + return NonPreserving() + end end const _WEAK_FUNCTIONS = (:+, :-) From 95a636ee908d8427faa2344cba8ca09aa4ed8981 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 14:22:17 -0400 Subject: [PATCH 05/14] Formatting --- src/indexing.jl | 4 +++- src/map.jl | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/indexing.jl b/src/indexing.jl index cbd13f7..7bd5419 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -309,7 +309,9 @@ end end # required: one implementation for canonical index style -@interface ::AbstractSparseArrayInterface function eachstoredindex(style::IndexStyle, A::AbstractArray) +@interface ::AbstractSparseArrayInterface function eachstoredindex( + style::IndexStyle, A::AbstractArray +) if style == IndexStyle(A) throw(MethodError(eachstoredindex, Tuple{typeof(style),typeof(A)})) elseif style == IndexCartesian() diff --git a/src/map.jl b/src/map.jl index 859d274..01e6d02 100644 --- a/src/map.jl +++ b/src/map.jl @@ -109,9 +109,9 @@ end # Derived functions # ----------------- -@interface I::AbstractSparseArrayInterface Base.copyto!( - C::AbstractArray, A::AbstractArray -) = @interface I map!(identity, C, A) +@interface I::AbstractSparseArrayInterface Base.copyto!(C::AbstractArray, A::AbstractArray) = @interface I map!( + identity, C, A +) # Utility functions # ----------------- From 9961f4c1d3dfd3f432e81bf740b2e77d3bec6fe5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 14:39:46 -0400 Subject: [PATCH 06/14] More careful with preserving functions --- src/map.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/map.jl b/src/map.jl index 01e6d02..fa31eb5 100644 --- a/src/map.jl +++ b/src/map.jl @@ -45,16 +45,14 @@ end const _WEAK_FUNCTIONS = (:+, :-) for f in _WEAK_FUNCTIONS @eval begin - ZeroPreserving(::typeof($f), ::AbstractArray, ::AbstractArray...) = WeakPreserving() - ZeroPreserving(::typeof($f), ::Type, ::Type...) = WeakPreserving() + ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = WeakPreserving() end end const _STRONG_FUNCTIONS = (:*,) for f in _STRONG_FUNCTIONS @eval begin - ZeroPreserving(::typeof($f), ::AbstractArray, ::AbstractArray...) = StrongPreserving() - ZeroPreserving(::typeof($f), ::Type, ::Type...) = StrongPreserving() + ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = StrongPreserving() end end From af5fd7a437f3a122865e9462d49fe7f90a5c458c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 15:27:02 -0400 Subject: [PATCH 07/14] Bump v0.6.0 --- Project.toml | 2 +- docs/Project.toml | 2 +- examples/Project.toml | 2 +- test/Project.toml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 7f585fd..ef6906d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseArraysBase" uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208" authors = ["ITensor developers and contributors"] -version = "0.5.11" +version = "0.6.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/docs/Project.toml b/docs/Project.toml index 898915b..ad4f567 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,4 +8,4 @@ SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" Dictionaries = "0.4.4" Documenter = "1.8.1" Literate = "2.20.1" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" diff --git a/examples/Project.toml b/examples/Project.toml index 01ec629..5c0e154 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Dictionaries = "0.4.4" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" Test = "<0.0.1, 1" diff --git a/test/Project.toml b/test/Project.toml index 9176f07..d7c7bb6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -21,7 +21,7 @@ JLArrays = "0.2.0" LinearAlgebra = "<0.0.1, 1" Random = "<0.0.1, 1" SafeTestsets = "0.1.0" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" StableRNGs = "1.0.2" Suppressor = "0.2.8" Test = "<0.0.1, 1" From 0e163b433bd762853fc8688cf099071187774269 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 16:52:03 -0400 Subject: [PATCH 08/14] Revert "Bump v0.6.0" This reverts commit af5fd7a437f3a122865e9462d49fe7f90a5c458c. --- Project.toml | 2 +- docs/Project.toml | 2 +- examples/Project.toml | 2 +- test/Project.toml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index ef6906d..7f585fd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseArraysBase" uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208" authors = ["ITensor developers and contributors"] -version = "0.6.0" +version = "0.5.11" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/docs/Project.toml b/docs/Project.toml index ad4f567..898915b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,4 +8,4 @@ SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" Dictionaries = "0.4.4" Documenter = "1.8.1" Literate = "2.20.1" -SparseArraysBase = "0.6.0" +SparseArraysBase = "0.5.0" diff --git a/examples/Project.toml b/examples/Project.toml index 5c0e154..01ec629 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Dictionaries = "0.4.4" -SparseArraysBase = "0.6.0" +SparseArraysBase = "0.5.0" Test = "<0.0.1, 1" diff --git a/test/Project.toml b/test/Project.toml index d7c7bb6..9176f07 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -21,7 +21,7 @@ JLArrays = "0.2.0" LinearAlgebra = "<0.0.1, 1" Random = "<0.0.1, 1" SafeTestsets = "0.1.0" -SparseArraysBase = "0.6.0" +SparseArraysBase = "0.5.0" StableRNGs = "1.0.2" Suppressor = "0.2.8" Test = "<0.0.1, 1" From 80f6805d5c6267689d9eae448811894a2212635a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 16:52:33 -0400 Subject: [PATCH 09/14] Bump v0.5.12 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7f585fd..06169e6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseArraysBase" uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208" authors = ["ITensor developers and contributors"] -version = "0.5.11" +version = "0.5.12" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From 4632408e9293e0e691a7fd32a5c2df32df36d40c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 16 Jul 2025 17:15:00 -0400 Subject: [PATCH 10/14] Refactor `map` to allow specifying kind of preserving --- src/map.jl | 84 ++++++++++++++++++++++++++++++------------------------ 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/src/map.jl b/src/map.jl index fa31eb5..437e9c9 100644 --- a/src/map.jl +++ b/src/map.jl @@ -1,7 +1,7 @@ # zero-preserving Traits # ---------------------- """ - abstract type ZeroPreserving end + abstract type ZeroPreserving <: Function end Holy Trait to indicate how a function interacts with abstract zero values: @@ -15,10 +15,17 @@ To attempt to automatically determine this, either `ZeroPreserving(f, A::Abstrac !!! warning incorrectly registering a function to be zero-preserving will lead to silently wrong results. """ -abstract type ZeroPreserving end -struct StrongPreserving <: ZeroPreserving end -struct WeakPreserving <: ZeroPreserving end -struct NonPreserving <: ZeroPreserving end +abstract type ZeroPreserving <: Function end + +struct StrongPreserving{F} <: ZeroPreserving + f::F +end +struct WeakPreserving{F} <: ZeroPreserving + f::F +end +struct NonPreserving{F} <: ZeroPreserving + f::F +end # Backport: remove in 1.12 @static if !isdefined(Base, :haszero) @@ -36,23 +43,25 @@ end # TODO: non-concrete element types function ZeroPreserving(f, T::Type, Ts::Type...) if all(_haszero, (T, Ts...)) - return iszero(f(zero(T), zero.(Ts)...)) ? WeakPreserving() : NonPreserving() + return iszero(f(zero(T), zero.(Ts)...)) ? WeakPreserving(f) : NonPreserving(f) else - return NonPreserving() + return NonPreserving(f) end end const _WEAK_FUNCTIONS = (:+, :-) for f in _WEAK_FUNCTIONS @eval begin - ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = WeakPreserving() + ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = WeakPreserving($f) end end const _STRONG_FUNCTIONS = (:*,) for f in _STRONG_FUNCTIONS @eval begin - ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = StrongPreserving() + ZeroPreserving(::typeof($f), ::Type{<:Number}, ::Type{<:Number}...) = StrongPreserving( + $f + ) end end @@ -61,47 +70,48 @@ end @interface I::AbstractSparseArrayInterface function Base.map( f, A::AbstractArray, Bs::AbstractArray... ) - T = Base.Broadcast.combine_eltypes(f, (A, Bs...)) + f_pres = ZeroPreserving(f, A, Bs...) + return @interface I map(f_pres, A, Bs...) +end +@interface I::AbstractSparseArrayInterface function Base.map( + f::ZeroPreserving, A::AbstractArray, Bs::AbstractArray... +) + T = Base.Broadcast.combine_eltypes(f.f, (A, Bs...)) C = similar(I, T, size(A)) return @interface I map!(f, C, A, Bs...) end -@interface ::AbstractSparseArrayInterface function Base.map!( +@interface I::AbstractSparseArrayInterface function Base.map!( f, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... ) - return _map!(f, ZeroPreserving(f, A, Bs...), C, A, Bs...) + f_pres = ZeroPreserving(f, A, Bs...) + return @interface I map!(f_pres, C, A, Bs...) end -function _map!( - f, ::StrongPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... -) - checkshape(C, A, Bs...) - style = IndexStyle(C, A, Bs...) - unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) - zero!(C) - for I in intersect(eachstoredindex.(Ref(style), unaliased)...) - @inbounds C[I] = f(ith_all(I, unaliased)...) - end - return C -end -function _map!( - f, ::WeakPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... +@interface ::AbstractSparseArrayInterface function Base.map!( + f::ZeroPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray... ) checkshape(C, A, Bs...) - style = IndexStyle(C, A, Bs...) unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) - zero!(C) - for I in union(eachstoredindex.(Ref(style), unaliased)...) - @inbounds C[I] = f(ith_all(I, unaliased)...) + + if f isa StrongPreserving + style = IndexStyle(C, unaliased...) + inds = intersect(eachstoredindex.(Ref(style), unaliased)...) + zero!(C) + elseif f isa WeakPreserving + style = IndexStyle(C, unaliased...) + inds = union(eachstoredindex.(Ref(style), unaliased)...) + zero!(C) + elseif f isa NonPreserving + inds = eachindex(C, unaliased...) + else + error(lazy"unknown zero-preserving type $(typeof(f))") end - return C -end -function _map!(f, ::NonPreserving, C::AbstractArray, A::AbstractArray, Bs::AbstractArray...) - checkshape(C, A, Bs...) - unaliased = map(Base.Fix1(Base.unalias, C), (A, Bs...)) - for I in eachindex(C, A, Bs...) - @inbounds C[I] = f(ith_all(I, unaliased)...) + + @inbounds for I in inds + C[I] = f.f(ith_all(I, unaliased)...) end + return C end From 1aa1964a175c505ca76c5d3bd3f74a52993663a9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Jul 2025 10:42:31 -0400 Subject: [PATCH 11/14] Refactor derived map functions --- src/abstractsparsearrayinterface.jl | 63 ----------------------------- src/map.jl | 20 +++++++++ 2 files changed, 20 insertions(+), 63 deletions(-) diff --git a/src/abstractsparsearrayinterface.jl b/src/abstractsparsearrayinterface.jl index b502b9a..699b976 100644 --- a/src/abstractsparsearrayinterface.jl +++ b/src/abstractsparsearrayinterface.jl @@ -92,39 +92,6 @@ end # TODO: Define `default_similartype` or something like that? return SparseArrayDOK{T}(undef, size) end - -# map over a specified subset of indices of the inputs. -function map_indices! end - -@interface interface::AbstractArrayInterface function map_indices!( - indices, f, a_dest::AbstractArray, as::AbstractArray... -) - for I in indices - a_dest[I] = f(map(a -> a[I], as)...) - end - return a_dest -end - -# Only map the stored values of the inputs. -function map_stored! end - -@interface interface::AbstractArrayInterface function map_stored!( - f, a_dest::AbstractArray, as::AbstractArray... -) - @interface interface map_indices!(eachstoredindex(as...), f, a_dest, as...) - return a_dest -end - -# Only map all values, not just the stored ones. -function map_all! end - -@interface interface::AbstractArrayInterface function map_all!( - f, a_dest::AbstractArray, as::AbstractArray... -) - @interface interface map_indices!(eachindex(as...), f, a_dest, as...) - return a_dest -end - using DerivableInterfaces: DerivableInterfaces, zero! # `zero!` isn't defined in `Base`, but it is defined in `ArrayLayouts` @@ -140,36 +107,6 @@ using DerivableInterfaces: DerivableInterfaces, zero! return @interface interface map_stored!(f, a, a) end -# Determines if a function preserves the stored values -# of the destination sparse array. -# The current code may be inefficient since it actually -# accesses an unstored element, which in the case of a -# sparse array of arrays can allocate an array. -# Sparse arrays could be expected to define a cheap -# unstored element allocator, for example -# `get_prototypical_unstored(a::AbstractArray)`. -function preserves_unstored(f, a_dest::AbstractArray, as::AbstractArray...) - I = first(eachindex(as...)) - return iszero(f(map(a -> getunstoredindex(a, I), as)...)) -end - -# @interface interface::AbstractSparseArrayInterface function Base.map!( -# f, a_dest::AbstractArray, as::AbstractArray... -# ) -# isempty(a_dest) && return a_dest # special case to avoid trying to access empty array -# indices = if !preserves_unstored(f, a_dest, as...) -# eachindex(a_dest) -# elseif any(a -> a_dest !== a, as) -# as = map(a -> Base.unalias(a_dest, a), as) -# @interface interface zero!(a_dest) -# eachstoredindex(as...) -# else -# eachstoredindex(a_dest) -# end -# @interface interface map_indices!(indices, f, a_dest, as...) -# return a_dest -# end - # `f::typeof(norm)`, `op::typeof(max)` used by `norm`. function reduce_init(f, op, as...) # TODO: Generalize this. diff --git a/src/map.jl b/src/map.jl index 437e9c9..f35426e 100644 --- a/src/map.jl +++ b/src/map.jl @@ -121,6 +121,26 @@ end identity, C, A ) +# Only map the stored values of the inputs. +function map_stored! end + +@interface interface::AbstractArrayInterface function map_stored!( + f, a_dest::AbstractArray, as::AbstractArray... +) + @interface interface map!(WeakPreserving(f), a_dest, as...) + return a_dest +end + +# Only map all values, not just the stored ones. +function map_all! end + +@interface interface::AbstractArrayInterface function map_all!( + f, a_dest::AbstractArray, as::AbstractArray... +) + @interface interface map!(NonPreserving(f), a_dest, as...) + return a_dest +end + # Utility functions # ----------------- # shape check similar to checkbounds From f6148a041c3bb50c5b78519333532771271397a6 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Jul 2025 11:27:27 -0400 Subject: [PATCH 12/14] Avoid infinite loop --- src/abstractsparsearrayinterface.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/abstractsparsearrayinterface.jl b/src/abstractsparsearrayinterface.jl index 699b976..719c94e 100644 --- a/src/abstractsparsearrayinterface.jl +++ b/src/abstractsparsearrayinterface.jl @@ -104,7 +104,10 @@ using DerivableInterfaces: DerivableInterfaces, zero! # More generally, this codepath could be taking if `zero(eltype(a))` # is defined and the elements are immutable. f = eltype(a) <: Number ? Returns(zero(eltype(a))) : zero! - return @interface interface map_stored!(f, a, a) + @inbounds for I in eachstoredindex(a) + a[I] = f(a[I]) + end + return a end # `f::typeof(norm)`, `op::typeof(max)` used by `norm`. From d316b9f19862a3652003ded82bc1f560c97c4126 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Jul 2025 11:28:15 -0400 Subject: [PATCH 13/14] Revert Revert Bump v0.6.0 --- Project.toml | 2 +- docs/Project.toml | 2 +- examples/Project.toml | 2 +- test/Project.toml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 06169e6..ef6906d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseArraysBase" uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208" authors = ["ITensor developers and contributors"] -version = "0.5.12" +version = "0.6.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/docs/Project.toml b/docs/Project.toml index 898915b..ad4f567 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,4 +8,4 @@ SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" Dictionaries = "0.4.4" Documenter = "1.8.1" Literate = "2.20.1" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" diff --git a/examples/Project.toml b/examples/Project.toml index 01ec629..5c0e154 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Dictionaries = "0.4.4" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" Test = "<0.0.1, 1" diff --git a/test/Project.toml b/test/Project.toml index 9176f07..d7c7bb6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -21,7 +21,7 @@ JLArrays = "0.2.0" LinearAlgebra = "<0.0.1, 1" Random = "<0.0.1, 1" SafeTestsets = "0.1.0" -SparseArraysBase = "0.5.0" +SparseArraysBase = "0.6.0" StableRNGs = "1.0.2" Suppressor = "0.2.8" Test = "<0.0.1, 1" From 8214b8f48576a8318740c0f2c8d4ba465faccb07 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Jul 2025 11:41:34 -0400 Subject: [PATCH 14/14] Refactor eachstoredindex implementation --- src/indexing.jl | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/indexing.jl b/src/indexing.jl index 7bd5419..2c2d6ca 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -308,18 +308,25 @@ end end end +@noinline function error_if_canonical_eachstoredindex(style::IndexStyle, A::AbstractArray) + style === IndexStyle(A) && throw(Base.CanonicalIndexError("eachstoredindex", typeof(A))) + return nothing +end + # required: one implementation for canonical index style @interface ::AbstractSparseArrayInterface function eachstoredindex( style::IndexStyle, A::AbstractArray ) - if style == IndexStyle(A) - throw(MethodError(eachstoredindex, Tuple{typeof(style),typeof(A)})) - elseif style == IndexCartesian() - return map(Base.Fix1(Base.getindex, CartesianIndices(A)), eachindex(A)) - elseif style == IndexLinear() - return map(Base.Fix1(Base.getindex, LinearIndices(A)), eachindex(A)) + error_if_canonical_eachstoredindex(style, A) + inds = eachstoredindex(A) + if style === IndexCartesian() + eltype(inds) === CartesianIndex{ndims(A)} && return inds + return map(Base.Fix1(Base.getindex, CartesianIndices(A)), inds) + elseif style === IndexLinear() + eltype(inds) === Int && return inds + return map(Base.Fix1(Base.getindex, LinearIndices(A)), inds) else - throw(ArgumentError(lazy"unknown index style $style")) + error(lazy"unkown index style $style") end end