From b447596fed7d0645b7c471e1c3bcc98ebba3071a Mon Sep 17 00:00:00 2001 From: Adrian Hill Date: Fri, 26 Jul 2024 17:41:50 +0200 Subject: [PATCH] Fix sparse array construction (#142) * Fix sparse array construction on Julia >=1.9 * Stop supporting sparse array construction on Julia <1.9 --- src/overloads/arrays.jl | 49 ++++-------------- test/test_arrays.jl | 112 ++++++++++++++++++++++------------------ 2 files changed, 72 insertions(+), 89 deletions(-) diff --git a/src/overloads/arrays.jl b/src/overloads/arrays.jl index 1b3a6b1..87d2c53 100644 --- a/src/overloads/arrays.jl +++ b/src/overloads/arrays.jl @@ -52,11 +52,6 @@ function split_dual_array(A::AbstractArray{D}) where {D<:Dual} tracers = getproperty.(A, :tracer) return primals, tracers end -function split_dual_array(A::SparseArrays.SparseMatrixCSC{D}) where {D<:Dual} - primals = getproperty.(A, :primal) - tracers = getproperty.(A, :tracer) - return primals, tracers -end #==================# # LinearAlgebra.jl # @@ -81,13 +76,6 @@ function LinearAlgebra.norm(A::AbstractArray{T}, p::Real=2) where {T<:AbstractTr return second_order_or(A) end end -function LinearAlgebra.opnorm(A::AbstractArray{T}, p::Real=2) where {T<:AbstractTracer} - if isone(p) || isinf(p) - return first_order_or(A) - else - return second_order_or(A) - end -end function LinearAlgebra.opnorm(A::AbstractMatrix{T}, p::Real=2) where {T<:AbstractTracer} if isone(p) || isinf(p) return first_order_or(A) @@ -196,31 +184,14 @@ end # SparseArrays # #==============# -# Conversion of matrices of tracers to SparseMatrixCSC has to be rewritten -# due to use of `count(_isnotzero, M)` in SparseArrays.jl -# -# Code modified from MIT licensed SparseArrays.jl source: -# https://github.com/JuliaSparse/SparseArrays.jl/blob/45dfe459ede2fa1419e7068d4bda92d9d22bd44d/src/sparsematrix.jl#L901-L920 -# Copyright (c) 2009-2024: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors: https://github.com/JuliaLang/julia/contributors -function SparseArrays.SparseMatrixCSC{Tv,Ti}( - M::StridedMatrix{Tv} -) where {Tv<:AbstractTracer,Ti} - nz = count(!isemptytracer, M) - colptr = zeros(Ti, size(M, 2) + 1) - nzval = Vector{Tv}(undef, nz) - rowval = Vector{Ti}(undef, nz) - colptr[1] = 1 - cnt = 1 - @inbounds for j in 1:size(M, 2) - for i in 1:size(M, 1) - v = M[i, j] - if !isemptytracer(v) - rowval[cnt] = i - nzval[cnt] = v - cnt += 1 - end - end - colptr[j + 1] = cnt - end - return SparseArrays.SparseMatrixCSC(size(M, 1), size(M, 2), colptr, rowval, nzval) +# Helper function needed in SparseArrays's sparsematrix, sparsevector and higherorderfns. +# On Tracers, `iszero` and `!iszero` don't return a boolean, +# but we need a function that does to handle the structure of the array. + +if VERSION >= v"1.9" # _iszero was added in JuliaSparse/SparseArrays.jl#177 + SparseArrays._iszero(t::AbstractTracer) = isemptytracer(t) + SparseArrays._iszero(d::Dual) = isemptytracer(tracer(d)) && iszero(primal(d)) + + SparseArrays._isnotzero(t::AbstractTracer) = !isemptytracer(t) + SparseArrays._isnotzero(d::Dual) = !isemptytracer(tracer(d)) || !iszero(primal(d)) end diff --git a/test/test_arrays.jl b/test/test_arrays.jl index a1871c2..c814d0f 100644 --- a/test/test_arrays.jl +++ b/test/test_arrays.jl @@ -66,25 +66,30 @@ end test_patterns(logabsdet_first, A) test_patterns(logabsdet_last, A; jac=iszero, hes=iszero) end - @testset "`SparseMatrixCSC` (3×3)" begin - # TODO: this is a temporary solution until sparse matrix inputs are supported (#28) - test_patterns(A -> det(sparse(A)), rand(3, 3)) - test_patterns(A -> logdet(sparse(A)), rand(3, 3)) - test_patterns(A -> norm(sparse(A)), rand(3, 3)) - test_patterns(A -> eigmax(sparse(A)), rand(3, 3)) - test_patterns(A -> eigmin(sparse(A)), rand(3, 3)) - test_patterns(A -> opnorm1(sparse(A)), rand(3, 3); hes=iszero) - test_patterns(A -> logabsdet_first(sparse(A)), rand(3, 3)) - test_patterns(A -> logabsdet_last(sparse(A)), rand(3, 3); jac=iszero, hes=iszero) - - test_patterns(v -> det(spdiagm(v)), rand(3)) - test_patterns(v -> logdet(spdiagm(v)), rand(3)) - test_patterns(v -> norm(spdiagm(v)), rand(3)) - test_patterns(v -> eigmax(spdiagm(v)), rand(3)) - test_patterns(v -> eigmin(spdiagm(v)), rand(3)) - test_patterns(v -> opnorm1(spdiagm(v)), rand(3); hes=iszero) - test_patterns(v -> logabsdet_first(spdiagm(v)), rand(3)) - test_patterns(v -> logabsdet_last(spdiagm(v)), rand(3); jac=iszero, hes=iszero) + + if VERSION >= v"1.9" + @testset "`SparseMatrixCSC` (3×3)" begin + # TODO: this is a temporary solution until sparse matrix inputs are supported (#28) + test_patterns(A -> det(sparse(A)), rand(3, 3)) + test_patterns(A -> logdet(sparse(A)), rand(3, 3)) + test_patterns(A -> norm(sparse(A)), rand(3, 3)) + test_patterns(A -> eigmax(sparse(A)), rand(3, 3)) + test_patterns(A -> eigmin(sparse(A)), rand(3, 3)) + test_patterns(A -> opnorm1(sparse(A)), rand(3, 3); hes=iszero) + test_patterns(A -> logabsdet_first(sparse(A)), rand(3, 3)) + test_patterns( + A -> logabsdet_last(sparse(A)), rand(3, 3); jac=iszero, hes=iszero + ) + + test_patterns(v -> det(spdiagm(v)), rand(3)) + test_patterns(v -> logdet(spdiagm(v)), rand(3)) + test_patterns(v -> norm(spdiagm(v)), rand(3)) + test_patterns(v -> eigmax(spdiagm(v)), rand(3)) + test_patterns(v -> eigmin(spdiagm(v)), rand(3)) + test_patterns(v -> opnorm1(spdiagm(v)), rand(3); hes=iszero) + test_patterns(v -> logabsdet_first(spdiagm(v)), rand(3)) + test_patterns(v -> logabsdet_last(spdiagm(v)), rand(3); jac=iszero, hes=iszero) + end end end @@ -99,33 +104,36 @@ end test_patterns(pow0, A; outsum=true, con=iszero, jac=iszero, hes=iszero) test_patterns(pow3, A; outsum=true) end - @testset "`SparseMatrixCSC` (3×3)" begin - # TODO: this is a temporary solution until sparse matrix inputs are supported (#28) - - test_patterns(A -> exp(sparse(A)), rand(3, 3); outsum=true) - test_patterns( - A -> pow0(sparse(A)), - rand(3, 3); - outsum=true, - con=iszero, - jac=iszero, - hes=iszero, - ) - test_patterns(A -> pow3(sparse(A)), rand(3, 3); outsum=true) - - test_patterns(v -> exp(spdiagm(v)), rand(3); outsum=true) - - if VERSION >= v"1.10" - # issue with custom _mapreducezeros in SparseArrays on Julia 1.6 + + if VERSION >= v"1.9" + @testset "`SparseMatrixCSC` (3×3)" begin + # TODO: this is a temporary solution until sparse matrix inputs are supported (#28) + + test_patterns(A -> exp(sparse(A)), rand(3, 3); outsum=true) test_patterns( - v -> pow0(spdiagm(v)), - rand(3); + A -> pow0(sparse(A)), + rand(3, 3); outsum=true, con=iszero, jac=iszero, hes=iszero, ) - test_patterns(v -> pow3(spdiagm(v)), rand(3); outsum=true) + test_patterns(A -> pow3(sparse(A)), rand(3, 3); outsum=true) + + test_patterns(v -> exp(spdiagm(v)), rand(3); outsum=true) + + if VERSION >= v"1.10" + # issue with custom _mapreducezeros in SparseArrays on Julia 1.6 + test_patterns( + v -> pow0(spdiagm(v)), + rand(3); + outsum=true, + con=iszero, + jac=iszero, + hes=iszero, + ) + test_patterns(v -> pow3(spdiagm(v)), rand(3); outsum=true) + end end end @@ -133,8 +141,10 @@ end @testset "$name" for (name, A) in TEST_MATRICES test_patterns(pinv, A; outsum=true) end - @testset "`SparseMatrixCSC` (3×4)" begin - test_patterns(A -> pinv(sparse(A)), rand(3, 4); outsum=true) + if VERSION >= v"1.9" + @testset "`SparseMatrixCSC` (3×4)" begin + test_patterns(A -> pinv(sparse(A)), rand(3, 4); outsum=true) + end end end @@ -165,13 +175,15 @@ end @test all(t -> SCT.gradient(t) == s_out, vectors) end -@testset "SparseMatrixCSC construction" begin - t1 = TG(P(S(1))) - t2 = TG(P(S(2))) - t3 = TG(P(S(3))) - SA = sparse([t1 t2; t3 0]) - @test length(SA.nzval) == 3 +if VERSION >= v"1.9" + @testset "SparseMatrixCSC construction" begin + t1 = TG(P(S(1))) + t2 = TG(P(S(2))) + t3 = TG(P(S(3))) + SA = sparse([t1 t2; t3 0]) + @test length(SA.nzval) == 3 - res = opnorm(SA, 1) - @test SCT.gradient(res) == S([1, 2, 3]) + res = opnorm(SA, 1) + @test SCT.gradient(res) == S([1, 2, 3]) + end end