From ba5a543ff982082a119a53a68f4e619de97dca60 Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Thu, 4 Dec 2025 14:08:25 +0000 Subject: [PATCH 1/6] initial commit --- ext/FillArraysSparseArraysExt.jl | 14 +++++- src/fillalgebra.jl | 23 +++++++++ src/oneelement.jl | 81 ++++++++++++++++++++++++++++++++ 3 files changed, 117 insertions(+), 1 deletion(-) diff --git a/ext/FillArraysSparseArraysExt.jl b/ext/FillArraysSparseArraysExt.jl index 4cb29dcd..59e22e48 100644 --- a/ext/FillArraysSparseArraysExt.jl +++ b/ext/FillArraysSparseArraysExt.jl @@ -4,7 +4,7 @@ using SparseArrays using SparseArrays: SparseVectorUnion import Base: convert, kron using FillArrays -using FillArrays: RectDiagonalFill, RectOrDiagonalFill, ZerosVector, ZerosMatrix, getindex_value, AbstractFillVector, _fill_dot +using FillArrays: RectDiagonalFill, RectOrDiagonalFill, ZerosVector, ZerosMatrix, getindex_value, AbstractFillVector, _fill_dot, OneElementVector # Specifying the full namespace is necessary because of https://github.com/JuliaLang/julia/issues/48533 # See https://github.com/JuliaStats/LogExpFunctions.jl/pull/63 using FillArrays.LinearAlgebra @@ -63,4 +63,16 @@ end # Ambiguity. see #178 dot(x::AbstractFillVector, y::SparseVectorUnion) = _fill_dot(x, y) +# OneElements with different indices should return +# SparseArrays under addition +function FillArrays.oneelement_addsub(a::OneElementVector, b::OneElementVector, aval, bval) + return sparsevec([a.ind[1], b.ind[1]], [aval, bval], length(a)) +end + +function FillArrays.oneelement_addsub(a::OneElement, b::OneElement, aval, bval) + nzval = [aval, bval] + nzind = ntuple(i -> [a.ind[i], b.ind[i]], ndims(a)) + return sparse(nzind..., nzval, size(a)...) +end + end # module diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl index deb7e74d..85738e7c 100644 --- a/src/fillalgebra.jl +++ b/src/fillalgebra.jl @@ -36,6 +36,29 @@ Base.@propagate_inbounds function reverse(A::AbstractFill, start::Integer, stop: end reverse(A::AbstractFill; dims=:) = A +## addition/subtraction with Zeros +# Zeros should function as an identity + +function add_zeros(z::AbstractZeros, v::AbstractArray) + axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) + return v +end + +function sub_zeros(v::AbstractArray, z::AbstractZeros) + axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) + return v +end + +function sub_zeros(z::AbstractZeros, v::AbstractArray) + axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) + return -v +end + ++(a::AbstractZeros, b::AbstractArray) = add_zeros(a, b) ++(a::AbstractArray, b::AbstractZeros) = add_zeros(b, a) +-(a::AbstractArray, b::AbstractZeros) = sub_zeros(a, b) +-(a::AbstractZeros, b::AbstractArray) = sub_zeros(a, b) + ## Algebraic identities # Default outputs, can overload to customize diff --git a/src/oneelement.jl b/src/oneelement.jl index f2db8565..b4b88c77 100644 --- a/src/oneelement.jl +++ b/src/oneelement.jl @@ -148,6 +148,87 @@ end /(x::OneElement, b::Number) = OneElement(x.val / b, x.ind, x.axes) \(b::Number, x::OneElement) = OneElement(b \ x.val, x.ind, x.axes) +# Addition + +# O(1) addition with arbitrary array types +function add_one_elem(a::OneElement, b::AbstractArray) + axes(a) == axes(b) || throw(DimensionMismatch(LazyString("A has dimensions ", size(a), " but B has dimensions ", size(b)))) + + ret = copy(b) + try + ret[a.ind...] += getindex_value(a) + catch + print("Exception") + # Fallback to materialising dense array if setindex! + # goes wrong (e.g on a Diagonal) + ret = Array(ret) + ret[a.ind...] += getindex_value(a) + end + return ret +end + +function sub_one_elem(a::AbstractArray, b::OneElement) + axes(a) == axes(b) || throw(DimensionMismatch(LazyString("A has dimensions ", size(a), " but B has dimensions ", size(b)))) + ret = copy(a) + try + ret[b.ind...] -= getindex_value(b) + catch + # Fallback to materialising dense array if setindex! + # goes wrong (e.g on a Diagonal) + ret = Array(ret) + ret[b.ind...] -= getindex_value(b) + end + return ret +end + +function sub_one_elem(a::OneElement, b::AbstractArray) + axes(a) == axes(b) || throw(DimensionMismatch(LazyString("A has dimensions ", size(a), " but B has dimensions ", size(b)))) + ret = copy(-b) + try + ret[a.ind...] += getindex_value(a) + catch + # Fallback to materialising dense array if setindex! + # goes wrong (e.g on a Diagonal) + ret = Array(ret) + ret[a.ind...] += getindex_value(a) + end + return ret +end + ++(a::OneElement, b::AbstractArray) = add_one_elem(a, b) ++(a::AbstractArray, b::OneElement) = add_one_elem(b, a) +-(a::AbstractArray, b::OneElement) = sub_one_elem(a, b) +-(a::OneElement, b::AbstractArray) = sub_one_elem(a, b) +# disambiguity ++(a::AbstractZeros, b::OneElement) = add_zeros(a, b) ++(a::OneElement, b::AbstractZeros) = add_zeros(b, a) + +# Adding/subtracting OneElements +# (Without SparseArrays) materialise dense vector if indices are different + +# This is overridden by the SparseArrays extension to return sparse arrays +function oneelement_addsub(a::OneElement, b::OneElement, aval, bval) + ret = similar(a) + fill!(ret, zero(eltype(ret))) + ret[a.ind...] = aval + ret[b.ind...] = bval + return ret +end + +for (op, bop) in (:+ => :(getindex_value(b)), + :- => :(-getindex_value(b))) + @eval begin + function $op(a::OneElement, b::OneElement) + axes(a) == axes(b) || throw(DimensionMismatch(LazyString("A has dimensions ", size(a), " but B has dimensions ", size(b)))) + if a.ind == b.ind + return OneElement($op(getindex_value(a), getindex_value(b)), a.ind, axes(a)) + else + return oneelement_addsub(a, b, getindex_value(a), $bop) + end + end + end +end + # matrix-vector and matrix-matrix multiplication # Fill and OneElement From 90e4cfa87216e74750e7e5a2c4da2457702a4c86 Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Thu, 4 Dec 2025 16:15:43 +0000 Subject: [PATCH 2/6] fix bugs, add staticarrays extension --- Project.toml | 3 ++ ext/FillArraysStaticArraysExt.jl | 28 ++++++++++++++++++ src/fillalgebra.jl | 49 ++++++++++++++------------------ src/oneelement.jl | 7 ++--- 4 files changed, 56 insertions(+), 31 deletions(-) create mode 100644 ext/FillArraysStaticArraysExt.jl diff --git a/Project.toml b/Project.toml index 17c28a38..9dba54ba 100644 --- a/Project.toml +++ b/Project.toml @@ -5,14 +5,17 @@ version = "1.15.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + [weakdeps] PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] FillArraysPDMatsExt = "PDMats" FillArraysSparseArraysExt = "SparseArrays" +FillArraysStaticArraysExt = "StaticArrays" FillArraysStatisticsExt = "Statistics" [compat] diff --git a/ext/FillArraysStaticArraysExt.jl b/ext/FillArraysStaticArraysExt.jl new file mode 100644 index 00000000..f120be1a --- /dev/null +++ b/ext/FillArraysStaticArraysExt.jl @@ -0,0 +1,28 @@ +module FillArraysStaticArraysExt + +using FillArrays +using StaticArrays + +import Base: promote_op +import FillArrays: elconvert + +# Disambiguity methods for StaticArrays + +function Base.:+(a::FillArrays.Zeros, b::StaticArray) + promote_shape(a,b) + return elconvert(promote_op(+,eltype(a),eltype(b)),b) +end +function Base.:+(a::StaticArray, b::FillArrays.Zeros) + promote_shape(a,b) + return elconvert(promote_op(+,eltype(a),eltype(b)),a) +end +function Base.:-(a::StaticArray, b::FillArrays.Zeros) + promote_shape(a,b) + return elconvert(promote_op(-,eltype(a),eltype(b)),a) +end +function Base.:-(a::FillArrays.Zeros, b::StaticArray) + promote_shape(a,b) + return elconvert(promote_op(-,eltype(a),eltype(b)),-b) +end + +end # module diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl index 85738e7c..47f470fb 100644 --- a/src/fillalgebra.jl +++ b/src/fillalgebra.jl @@ -36,29 +36,6 @@ Base.@propagate_inbounds function reverse(A::AbstractFill, start::Integer, stop: end reverse(A::AbstractFill; dims=:) = A -## addition/subtraction with Zeros -# Zeros should function as an identity - -function add_zeros(z::AbstractZeros, v::AbstractArray) - axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) - return v -end - -function sub_zeros(v::AbstractArray, z::AbstractZeros) - axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) - return v -end - -function sub_zeros(z::AbstractZeros, v::AbstractArray) - axes(z) == axes(v) || throw(DimensionMismatch(LazyString("A has dimensions ", size(z), " but B has dimensions ", size(v)))) - return -v -end - -+(a::AbstractZeros, b::AbstractArray) = add_zeros(a, b) -+(a::AbstractArray, b::AbstractZeros) = add_zeros(b, a) --(a::AbstractArray, b::AbstractZeros) = sub_zeros(a, b) --(a::AbstractZeros, b::AbstractArray) = sub_zeros(a, b) - ## Algebraic identities # Default outputs, can overload to customize @@ -459,14 +436,32 @@ function +(a::AbstractZeros{T}, b::AbstractZeros{V}) where {T, V} # for disambig promote_shape(a,b) return elconvert(promote_op(+,T,V),a) end -# no AbstractArray. Otherwise incompatible with StaticArrays.jl +# concrete Zeros used to stop a variety of ambiguity issues # AbstractFill for disambiguity -for TYPE in (:Array, :AbstractFill, :AbstractRange, :Diagonal) - @eval function +(a::$TYPE{T}, b::AbstractZeros{V}) where {T, V} +for TYPE in (:Array, :AbstractFill, :AbstractRange, :AbstractArray) + @eval function +(a::$TYPE{T}, b::Zeros{V}) where {T, V} promote_shape(a,b) return elconvert(promote_op(+,T,V),a) end - @eval +(a::AbstractZeros, b::$TYPE) = b + a + @eval function -(a::$TYPE{T}, b::Zeros{V}) where {T, V} + promote_shape(a,b) + return elconvert(promote_op(-,T,V),a) + end + @eval function -(a::Zeros{T}, b::$TYPE{V}) where {T, V} + promote_shape(a,b) + return elconvert(promote_op(-,T,V),-b) + end + @eval +(a::Zeros, b::$TYPE) = b + a +end + +function +(a::Zeros, b::Zeros) + promote_shape(a,b) + return elconvert(promote_type(eltype(a), eltype(b)),a) +end + +function -(a::Zeros, b::Zeros) + promote_shape(a,b) + return elconvert(promote_type(eltype(a), eltype(b)),-b) end # for VERSION other than 1.6, could use ZerosMatrix only diff --git a/src/oneelement.jl b/src/oneelement.jl index b4b88c77..a5f83352 100644 --- a/src/oneelement.jl +++ b/src/oneelement.jl @@ -158,11 +158,10 @@ function add_one_elem(a::OneElement, b::AbstractArray) try ret[a.ind...] += getindex_value(a) catch - print("Exception") # Fallback to materialising dense array if setindex! # goes wrong (e.g on a Diagonal) ret = Array(ret) - ret[a.ind...] += getindex_value(a) + checkbounds(Bool, ret, a.ind...) && (ret[a.ind...] += getindex_value(a)) end return ret end @@ -176,7 +175,7 @@ function sub_one_elem(a::AbstractArray, b::OneElement) # Fallback to materialising dense array if setindex! # goes wrong (e.g on a Diagonal) ret = Array(ret) - ret[b.ind...] -= getindex_value(b) + checkbounds(Bool, ret, b.ind...) && (ret[b.ind...] -= getindex_value(b)) end return ret end @@ -190,7 +189,7 @@ function sub_one_elem(a::OneElement, b::AbstractArray) # Fallback to materialising dense array if setindex! # goes wrong (e.g on a Diagonal) ret = Array(ret) - ret[a.ind...] += getindex_value(a) + checkbounds(Bool, ret, a.ind...) && (ret[a.ind...] += getindex_value(a)) end return ret end From a99e4aab31a96325f0ee6decf133cd164d85091a Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Thu, 4 Dec 2025 16:22:48 +0000 Subject: [PATCH 3/6] revert to AbstractZeros --- src/fillalgebra.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl index 47f470fb..62e390ef 100644 --- a/src/fillalgebra.jl +++ b/src/fillalgebra.jl @@ -436,30 +436,25 @@ function +(a::AbstractZeros{T}, b::AbstractZeros{V}) where {T, V} # for disambig promote_shape(a,b) return elconvert(promote_op(+,T,V),a) end -# concrete Zeros used to stop a variety of ambiguity issues -# AbstractFill for disambiguity + +# AbstractFill and Array for disambiguity for TYPE in (:Array, :AbstractFill, :AbstractRange, :AbstractArray) - @eval function +(a::$TYPE{T}, b::Zeros{V}) where {T, V} + @eval function +(a::$TYPE{T}, b::AbstractZeros{V}) where {T, V} promote_shape(a,b) return elconvert(promote_op(+,T,V),a) end - @eval function -(a::$TYPE{T}, b::Zeros{V}) where {T, V} + @eval function -(a::$TYPE{T}, b::AbstractZeros{V}) where {T, V} promote_shape(a,b) return elconvert(promote_op(-,T,V),a) end - @eval function -(a::Zeros{T}, b::$TYPE{V}) where {T, V} + @eval function -(a::AbstractZeros{T}, b::$TYPE{V}) where {T, V} promote_shape(a,b) return elconvert(promote_op(-,T,V),-b) end - @eval +(a::Zeros, b::$TYPE) = b + a -end - -function +(a::Zeros, b::Zeros) - promote_shape(a,b) - return elconvert(promote_type(eltype(a), eltype(b)),a) + @eval +(a::AbstractZeros, b::$TYPE) = b + a end -function -(a::Zeros, b::Zeros) +function -(a::AbstractZeros, b::AbstractZeros) promote_shape(a,b) return elconvert(promote_type(eltype(a), eltype(b)),-b) end From 159119fa0dbcf9529e4db1da325bb666f3b725c6 Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Fri, 5 Dec 2025 11:35:28 +0000 Subject: [PATCH 4/6] add unit tests --- test/runtests.jl | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 174075cc..7ab9585a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -505,6 +505,14 @@ end for A in As, Z in (TZ -> Zeros{TZ}(3)).((Int, Float64, Int8, ComplexF64)) test_addition_and_subtraction_dim_mismatch(A, Z) end + + # Zeros should be additive identity for matrices + D = Diagonal([1, 1]) + Z = Zeros(2, 2) + @test D + Z isa Diagonal + @test D + Z == D + @test D - Z == D + @test Z - D == -D end end @@ -588,6 +596,12 @@ end testsparsediag(E) end end + + # Adding one elements with different indices should return a sparse array. + A = OneElement(2, 5) + B = OneElement(3, 5) + @test A + B isa SparseVector + @test A + B == [0, 1, 1, 0, 0] end @testset "==" begin @@ -2374,6 +2388,19 @@ end @test_throws ArgumentError isassigned(f, true) end + @testset "Addition/Subtraction" begin + A = OneElement(2, 5) + B = OneElement(3, 5) + + @test A + A isa OneElement + @test A + A == OneElement(2, 2, 5) + @test A + B == [0, 1, 1, 0, 0] + + @test B - B isa OneElement + @test B - B == OneElement(0, 2, 5) + @test B - A == [0, -1, 1, 0, 0] + end + @testset "matmul" begin A = reshape(Float64[1:9;], 3, 3) v = reshape(Float64[1:3;], 3) From 6d37f34743cf2b9ca0c9b0c9aea2e387e7c6cd02 Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Fri, 5 Dec 2025 14:00:09 +0000 Subject: [PATCH 5/6] fix redundant method --- src/oneelement.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/oneelement.jl b/src/oneelement.jl index a5f83352..d40d82d3 100644 --- a/src/oneelement.jl +++ b/src/oneelement.jl @@ -199,8 +199,12 @@ end -(a::AbstractArray, b::OneElement) = sub_one_elem(a, b) -(a::OneElement, b::AbstractArray) = sub_one_elem(a, b) # disambiguity -+(a::AbstractZeros, b::OneElement) = add_zeros(a, b) -+(a::OneElement, b::AbstractZeros) = add_zeros(b, a) +function +(a::AbstractZeros, b::OneElement) + promote_shape(a,b) + return elconvert(promote_op(+,eltype(a),eltype(b)),b) +end + ++(a::OneElement, b::AbstractZeros) = b + a # Adding/subtracting OneElements # (Without SparseArrays) materialise dense vector if indices are different From 33466e1c757b396cea96cd8f8abc970e11948fc7 Mon Sep 17 00:00:00 2001 From: Max Vassiliev Date: Fri, 5 Dec 2025 14:42:43 +0000 Subject: [PATCH 6/6] fix sparse arrays extension --- Project.toml | 1 - ext/FillArraysSparseArraysExt.jl | 4 ++-- src/oneelement.jl | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 9dba54ba..0a5c562e 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "1.15.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - [weakdeps] PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/ext/FillArraysSparseArraysExt.jl b/ext/FillArraysSparseArraysExt.jl index 59e22e48..68bfc5dd 100644 --- a/ext/FillArraysSparseArraysExt.jl +++ b/ext/FillArraysSparseArraysExt.jl @@ -4,7 +4,7 @@ using SparseArrays using SparseArrays: SparseVectorUnion import Base: convert, kron using FillArrays -using FillArrays: RectDiagonalFill, RectOrDiagonalFill, ZerosVector, ZerosMatrix, getindex_value, AbstractFillVector, _fill_dot, OneElementVector +using FillArrays: RectDiagonalFill, RectOrDiagonalFill, ZerosVector, ZerosMatrix, getindex_value, AbstractFillVector, _fill_dot, OneElementVector, OneElementMatrix # Specifying the full namespace is necessary because of https://github.com/JuliaLang/julia/issues/48533 # See https://github.com/JuliaStats/LogExpFunctions.jl/pull/63 using FillArrays.LinearAlgebra @@ -69,7 +69,7 @@ function FillArrays.oneelement_addsub(a::OneElementVector, b::OneElementVector, return sparsevec([a.ind[1], b.ind[1]], [aval, bval], length(a)) end -function FillArrays.oneelement_addsub(a::OneElement, b::OneElement, aval, bval) +function FillArrays.oneelement_addsub(a::OneElementMatrix, b::OneElementMatrix, aval, bval) nzval = [aval, bval] nzind = ntuple(i -> [a.ind[i], b.ind[i]], ndims(a)) return sparse(nzind..., nzval, size(a)...) diff --git a/src/oneelement.jl b/src/oneelement.jl index d40d82d3..dabfbf06 100644 --- a/src/oneelement.jl +++ b/src/oneelement.jl @@ -209,7 +209,7 @@ end # Adding/subtracting OneElements # (Without SparseArrays) materialise dense vector if indices are different -# This is overridden by the SparseArrays extension to return sparse arrays +# Sparse Arrays extension overrides this for OneElementVector and OneElementMatrix function oneelement_addsub(a::OneElement, b::OneElement, aval, bval) ret = similar(a) fill!(ret, zero(eltype(ret)))