Skip to content

Commit

Permalink
Fix sparse array construction (#142)
Browse files Browse the repository at this point in the history 
* Fix sparse array construction on Julia >=1.9

* Stop supporting sparse array construction on Julia <1.9
  • Loading branch information
@adrhill
adrhill authored Jul 26, 2024
1 parent ac94586 commit b447596
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 89 deletions.
49 changes: 10 additions & 39 deletions src/overloads/arrays.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 #
Expand All @@ -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)
Expand Down Expand Up @@ -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
112 changes: 62 additions & 50 deletions test/test_arrays.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -99,42 +104,47 @@ 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

# Functions that work on all matrices
@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

Expand Down Expand Up @@ -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

0 comments on commit b447596

Please sign in to comment.