From 5e7d88ec271df4bd12accc16eb56e7a9e14043fb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 17 Jan 2025 04:03:18 -0600 Subject: [PATCH] Coercing constructors & convert (#212) This adds a "constructor cascade" that allows one to coerce the inputs to the chosen type. This is frequently useful in generic programming where you want an object to be of the same type as something else, for example when appending to a list of objects. Co-authored-by: David Widmann --- Project.toml | 2 +- src/chol.jl | 2 +- src/pdmat.jl | 65 +++++++++++++++++++++++++++++++------------ test/pdmtypes.jl | 25 ++++++++++++++++- test/specialarrays.jl | 8 ++++++ 5 files changed, 81 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 811c890..11e2f07 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.31" +version = "0.11.32" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/chol.jl b/src/chol.jl index ec3ec0c..92c679e 100644 --- a/src/chol.jl +++ b/src/chol.jl @@ -73,7 +73,7 @@ function invquad(A::Cholesky, x::AbstractVector) @check_argdims size(A, 1) == size(x, 1) return sum(abs2, chol_lower(A) \ x) end -function invquad(A::Cholesky, X::AbstractMatrix) +function invquad(A::Cholesky, X::AbstractMatrix) @check_argdims size(A, 1) == size(X, 1) Z = chol_lower(A) \ X return vec(sum(abs2, Z; dims=1)) diff --git a/src/pdmat.jl b/src/pdmat.jl index 99c2deb..9c42844 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -1,23 +1,54 @@ """ Full positive definite matrix together with a Cholesky factorization object. """ -struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T} +struct PDMat{T<:Real,S<:AbstractMatrix{T}} <: AbstractPDMat{T} mat::S chol::Cholesky{T,S} - PDMat{T,S}(m::AbstractMatrix{T},c::Cholesky{T,S}) where {T,S} = new{T,S}(m,c) -end - -function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S} - d = LinearAlgebra.checksquare(mat) - if size(chol, 1) != d - throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) + function PDMat{T,S}(m::AbstractMatrix, c::Cholesky) where {T<:Real,S<:AbstractMatrix{T}} + d = LinearAlgebra.checksquare(m) + if size(c, 1) != d + throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) + end + # in principle we might want to check that `c` is a Cholesky factorization of `m`, + # but that's slow + return new{T,S}(m,c) end - PDMat{T,S}(convert(S, mat), chol) end +function PDMat{T}(m::AbstractMatrix, c::Cholesky) where T<:Real + c = convert(Cholesky{T}, c) + return PDMat{T,typeof(c.factors)}(m, c) +end +PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(mat, chol) +# Construction from another PDMat +PDMat{T,S}(pdm::PDMat{T,S}) where {T<:Real,S<:AbstractMatrix{T}} = pdm # since PDMat doesn't support `setindex!` it's not mutable (xref https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#Mutable-collections) +PDMat{T,S}(pdm::PDMat) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(pdm.mat, pdm.chol) +PDMat{T}(pdm::PDMat{T}) where T<:Real = pdm +PDMat{T}(pdm::PDMat) where T<:Real = PDMat{T}(pdm.mat, pdm.chol) +PDMat(pdm::PDMat) = pdm + +# Construction from an AbstractMatrix +function PDMat{T,S}(mat::AbstractMatrix) where {T<:Real,S<:AbstractMatrix{T}} + mat = convert(S, mat) + return PDMat{T,S}(mat, cholesky(mat)) +end +function PDMat{T}(mat::AbstractMatrix) where T<:Real + mat = convert(AbstractMatrix{T}, mat) + return PDMat{T}(mat, cholesky(mat)) +end PDMat(mat::AbstractMatrix) = PDMat(mat, cholesky(mat)) -PDMat(fac::Cholesky) = PDMat(AbstractMatrix(fac), fac) + +# Construction from a Cholesky factorization +function PDMat{T,S}(c::Cholesky) where {T<:Real,S<:AbstractMatrix{T}} + c = convert(Cholesky{T,S}, c) + return PDMat{T,S}(AbstractMatrix(c), c) +end +function PDMat{T}(c::Cholesky) where T<:Real + c = convert(Cholesky{T}, c) + return PDMat{T}(AbstractMatrix(c), c) +end +PDMat(c::Cholesky) = PDMat(AbstractMatrix(c), c) function Base.getproperty(a::PDMat, s::Symbol) if s === :dim @@ -30,13 +61,11 @@ Base.propertynames(::PDMat) = (:mat, :chol, :dim) AbstractPDMat(A::Cholesky) = PDMat(A) ### Conversion -Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a -function Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} - chol = convert(Cholesky{T}, a.chol) - S = typeof(chol.factors) - mat = convert(S, a.mat) - return PDMat{T,S}(mat, chol) -end +# This next method isn't needed because PDMat{T}(a) returns `a` directly +# Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a +Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} = PDMat{T}(a) +Base.convert(::Type{PDMat{T,S}}, a::PDMat) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(a) + Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat{T}, a) ### Basics @@ -55,7 +84,7 @@ Base.IndexStyle(::Type{PDMat{T,S}}) where {T,S} = Base.IndexStyle(S) # Linear Indexing Base.@propagate_inbounds Base.getindex(a::PDMat, i::Int) = getindex(a.mat, i) # Cartesian Indexing -Base.@propagate_inbounds Base.getindex(a::PDMat, I::Vararg{Int, 2}) = getindex(a.mat, I...) +Base.@propagate_inbounds Base.getindex(a::PDMat, i::Int, j::Int) = getindex(a.mat, i, j) ### Arithmetics diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 7be3898..8d1a78d 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -18,13 +18,36 @@ using Test M = convert(Array{T,2}, [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]) V = convert(Array{T,1}, [1.5, 2.5, 2.0]) X = convert(T,2.0) + f64M = Float64.(M) @testset "PDMat from Matrix" begin + pdf64M = PDMat(f64M) test_pdmat(PDMat(M), M, cmat_eq=true, verbose=1) + test_pdmat(PDMat{Float64}(M), f64M, cmat_eq=true, verbose=1) + test_pdmat(PDMat{Float64,Matrix{Float64}}(M), f64M, cmat_eq=true, verbose=1) + @test_throws TypeError PDMat{Float32,Matrix{Float64}}(M) + end + @testset "PDMat from PDMat" begin + pdM = PDMat(M) + pdf64M = PDMat(f64M) + @test pdM === PDMat(pdM) + @test pdf64M === PDMat{Float64}(pdf64M) === PDMat{Float64,Matrix{Float64}}(pdf64M) + test_pdmat(PDMat(pdM), M, cmat_eq=true, verbose=1) + test_pdmat(PDMat{Float64}(pdf64M), f64M, cmat_eq=true, verbose=1) + test_pdmat(PDMat{Float64,Matrix{Float64}}(pdf64M), f64M, cmat_eq=true, verbose=1) + if Base.VERSION >= v"1.12.0-DEV.1654" # julia #56562 + @test isa(convert(typeof(pdf64M), pdM), typeof(pdf64M)) + end + @test_throws TypeError PDMat{Float32,Matrix{Float64}}(pdM) end @testset "PDMat from Cholesky" begin cholL = Cholesky(Matrix(transpose(cholesky(M).factors)), 'L', 0) + cholLf64 = Cholesky(Matrix(transpose(cholesky(f64M).factors)), 'L', 0) test_pdmat(PDMat(cholL), M, cmat_eq=true, verbose=1) + test_pdmat(PDMat{Float64}(cholLf64), f64M, cmat_eq=true, verbose=1) + if Base.VERSION >= v"1.12.0-DEV.1654" # julia #56562 + test_pdmat(PDMat{Float64,Matrix{Float64}}(cholLf64), f64M, cmat_eq=true, verbose=1) + end end @testset "PDiagMat" begin test_pdmat(PDiagMat(V), Matrix(Diagonal(V)), cmat_eq=true, verbose=1) @@ -142,7 +165,7 @@ using Test # right division not defined for CHOLMOD: # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined if !HAVE_CHOLMOD - z = x / PDSparseMat(sparse(first(A), 1, 1)) + z = x / PDSparseMat(sparse(first(A), 1, 1)) @test typeof(z) === typeof(y) @test size(z) == size(y) @test z ≈ y diff --git a/test/specialarrays.jl b/test/specialarrays.jl index 8b4aaad..8342df2 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -15,6 +15,14 @@ using StaticArrays @test PDMat(S, C) === PDS @test @allocated(PDMat(S)) == @allocated(PDMat(C)) == @allocated(PDMat(S, C)) + if Base.VERSION >= v"1.12.0-DEV.1654" # julia #56562 + A = PDMat(Matrix{Float64}(I, 2, 2)) + B = PDMat(SMatrix{2,2,Float64}(I)) + @test !isa(A.mat, typeof(B.mat)) + S = convert(typeof(B), A) + @test isa(S.mat, typeof(B.mat)) + end + # Diagonal matrix D = PDiagMat(@SVector(rand(4))) @test D isa PDiagMat{Float64, <:SVector{4, Float64}}