Skip to content

Commit

Permalink
Handle dependencies via Pkg extensions (#208)
Browse files Browse the repository at this point in the history
  • Loading branch information
dkarrasch authored Jun 27, 2023
1 parent f522a8f commit b6b1044
Show file tree
Hide file tree
Showing 8 changed files with 140 additions and 108 deletions.
6 changes: 5 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "LinearMaps"
uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
version = "3.10.1"
version = "3.10.2"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand All @@ -10,9 +10,13 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[weakdeps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[extensions]
LinearMapsChainRulesCoreExt = "ChainRulesCore"
LinearMapsSparseArraysExt = "SparseArrays"
LinearMapsStatisticsExt = "Statistics"

[compat]
ChainRulesCore = "1"
Expand Down
117 changes: 117 additions & 0 deletions ext/LinearMapsSparseArraysExt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
module LinearMapsSparseArraysExt

import SparseArrays: sparse, blockdiag, SparseMatrixCSC
using SparseArrays: AbstractSparseMatrix

using LinearAlgebra, LinearMaps
import LinearMaps: _issymmetric, _ishermitian
using LinearMaps: WrappedMap, CompositeMap, LinearCombination, ScaledMap, UniformScalingMap,
AdjointMap, TransposeMap, BlockMap, BlockDiagonalMap, KroneckerMap, KroneckerSumMap,
VecOrMatMap, AbstractVecOrMatOrQ, MapOrVecOrMat
using LinearMaps: convert_to_lmaps, _tail, _unsafe_mul!

_issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
_ishermitian(A::AbstractSparseMatrix) = ishermitian(A)

# blockdiagonal concatenation via extension of blockdiag

"""
blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap
Construct a (lazy) representation of the diagonal concatenation of the arguments.
To avoid fallback to the generic `blockdiag`, there must be a `LinearMap`
object among the first 8 arguments.
"""
blockdiag

for k in 1:8 # is 8 sufficient?
Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMatOrQ), Val(k-1))
# yields (:A1, :A2, :A3, ..., :A(k-1))
L = :($(Symbol(:A, k))::LinearMap)
# yields :Ak
mapargs = ntuple(n ->:($(Symbol(:A, n))), Val(k-1))
# yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))

@eval function blockdiag($(Is...), $L, As::MapOrVecOrMat...)
return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
$(Symbol(:A, k)),
convert_to_lmaps(As...)...)
end
end

# conversion to sparse arrays
# sparse: create sparse matrix representation of LinearMap
function sparse(A::LinearMap{T}) where {T}
M, N = size(A)
rowind = Int[]
nzval = T[]
colptr = Vector{Int}(undef, N+1)
v = fill(zero(T), N)
Av = Vector{T}(undef, M)

@inbounds for i in eachindex(v)
v[i] = one(T)
_unsafe_mul!(Av, A, v)
js = findall(!iszero, Av)
colptr[i] = length(nzval) + 1
if length(js) > 0
append!(rowind, js)
append!(nzval, Av[js])
end
v[i] = zero(T)
end
colptr[N+1] = length(nzval) + 1

return SparseMatrixCSC(M, N, colptr, rowind, nzval)
end
Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
SparseMatrixCSC(A::LinearMap) = sparse(A)

sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) =
A.λ * sparse(A.lmap.lmap)
sparse(A::WrappedMap) = sparse(A.lmap)
Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
for (T, t) in ((:AdjointMap, adjoint), (:TransposeMap, transpose))
@eval sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap))
end
function sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}})
mats = map(A->getfield(A, :lmap), ΣA.maps)
return sum(sparse, mats)
end
function sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}})
B, A = AB.maps
return sparse(A.lmap)*sparse(B.lmap)
end
function sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}})
A, J = λA.maps
return J.λ*sparse(A.lmap)
end
function sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}})
J, A =.maps
return sparse(A.lmap)*J.λ
end
function sparse(A::BlockMap)
return hvcat(
A.rows,
convert(SparseMatrixCSC, first(A.maps)),
convert.(AbstractArray, _tail(A.maps))...
)
end
function sparse(A::BlockDiagonalMap)
return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
end
Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
function sparse(A::KroneckerMap)
return kron(
convert(SparseMatrixCSC, first(A.maps)),
convert.(AbstractMatrix, _tail(A.maps))...
)
end
function sparse(L::KroneckerSumMap)
A, B = L.maps
IA = sparse(Diagonal(ones(Bool, size(A, 1))))
IB = sparse(Diagonal(ones(Bool, size(B, 1))))
return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
end

end # module LinearMapsSparseArraysExt
12 changes: 12 additions & 0 deletions ext/LinearMapsStatisticsExt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
module LinearMapsStatisticsExt

import Statistics: mean

using LinearMaps
using LinearMaps: LinearMapTupleOrVector, LinearCombination

mean(f::F, maps::LinearMapTupleOrVector) where {F} = sum(f, maps) / length(maps)
mean(maps::LinearMapTupleOrVector) = mean(identity, maps)
mean(A::LinearCombination) = mean(A.maps)

end # module ChainRulesCore
5 changes: 2 additions & 3 deletions src/LinearMaps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,6 @@ export ⊗, squarekron, kronsum, ⊕, sumkronsum, khatrirao, facesplitting
using LinearAlgebra
using LinearAlgebra: AbstractQ
import LinearAlgebra: mul!
using SparseArrays

import Statistics: mean

using Base: require_one_based_indexing

Expand Down Expand Up @@ -422,6 +419,8 @@ LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwarg

@static if !isdefined(Base, :get_extension)
include("../ext/LinearMapsChainRulesCoreExt.jl")
include("../ext/LinearMapsSparseArraysExt.jl")
include("../ext/LinearMapsStatisticsExt.jl")
end

end # module
25 changes: 4 additions & 21 deletions src/blockmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -505,34 +505,17 @@ for k in 1:8 # is 8 sufficient?
mapargs = ntuple(n ->:($(Symbol(:A, n))), Val(k-1))
# yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))

@eval begin
function SparseArrays.blockdiag($(Is...), $L, As::MapOrVecOrMat...)
@eval function Base.cat($(Is...), $L, As::MapOrVecOrMat...; dims::Dims{2})
if dims == (1,2)
return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
$(Symbol(:A, k)),
convert_to_lmaps(As...)...)
end

function Base.cat($(Is...), $L, As::MapOrVecOrMat...; dims::Dims{2})
if dims == (1,2)
return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
$(Symbol(:A, k)),
convert_to_lmaps(As...)...)
else
throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
end
else
throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
end
end
end

"""
blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap
Construct a (lazy) representation of the diagonal concatenation of the arguments.
To avoid fallback to the generic `SparseArrays.blockdiag`, there must be a `LinearMap`
object among the first 8 arguments.
"""
SparseArrays.blockdiag

"""
cat(As::Union{LinearMap,AbstractVecOrMatOrQ}...; dims=(1,2))::BlockDiagonalMap
Expand Down
77 changes: 0 additions & 77 deletions src/conversion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,39 +19,8 @@ Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A)
Base.convert(::Type{AbstractMatrix}, A::LinearMap) = AbstractMatrix(A)
Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A)

# sparse: create sparse matrix representation of LinearMap
function SparseArrays.sparse(A::LinearMap{T}) where {T}
M, N = size(A)
rowind = Int[]
nzval = T[]
colptr = Vector{Int}(undef, N+1)
v = fill(zero(T), N)
Av = Vector{T}(undef, M)

@inbounds for i in eachindex(v)
v[i] = one(T)
_unsafe_mul!(Av, A, v)
js = findall(!iszero, Av)
colptr[i] = length(nzval) + 1
if length(js) > 0
append!(rowind, js)
append!(nzval, Av[js])
end
v[i] = zero(T)
end
colptr[N+1] = length(nzval) + 1

return SparseMatrixCSC(M, N, colptr, rowind, nzval)
end
Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A)

# special cases

# ScaledMap
SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) =
A.λ * sparse(A.lmap.lmap)

# UniformScalingMap
Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, J.M))

Expand All @@ -61,19 +30,10 @@ Base.convert(::Type{T}, A::WrappedMap) where {T<:Matrix} = convert(T, A.lmap)
Base.Matrix{T}(A::VectorMap) where {T} = copyto!(Matrix{eltype(T)}(undef, size(A)), A.lmap)
Base.convert(::Type{T}, A::VectorMap) where {T<:Matrix} = T(A)
Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)

# TransposeMap & AdjointMap
for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose))
@eval Base.convert(::Type{AbstractMatrix}, A::$T) = $t(convert(AbstractMatrix, A.lmap))
@eval SparseArrays.sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap))
end

# LinearCombination
function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}})
mats = map(A->getfield(A, :lmap), ΣA.maps)
return sum(sparse, mats)
end

# CompositeMap
Expand All @@ -99,50 +59,19 @@ function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap
B, A = AB.maps
return mul!(Matrix{T}(undef, size(AB)), A.lmap, B.lmap)
end
function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}})
B, A = AB.maps
return sparse(A.lmap)*sparse(B.lmap)
end
function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) where {T}
A, J = λA.maps
return mul!(Matrix{T}(undef, size(λA)), J.λ, A.lmap)
end
function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}})
A, J = λA.maps
return J.λ*sparse(A.lmap)
end
function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) where {T}
J, A =.maps
return mul!(Matrix{T}(undef, size(Aλ)), A.lmap, J.λ)
end
function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}})
J, A =.maps
return sparse(A.lmap)*J.λ
end

# BlockMap & BlockDiagonalMap
function SparseArrays.sparse(A::BlockMap)
return hvcat(
A.rows,
convert(SparseMatrixCSC, first(A.maps)),
convert.(AbstractArray, _tail(A.maps))...
)
end
Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
function SparseArrays.sparse(A::BlockDiagonalMap)
return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
end

# KroneckerMap & KroneckerSumMap
Base.Matrix{T}(A::KroneckerMap) where {T} = kron(convert.(Matrix{T}, A.maps)...)
Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) =
kron(convert.(AbstractMatrix, A.maps)...)
function SparseArrays.sparse(A::KroneckerMap)
return kron(
convert(SparseMatrixCSC, first(A.maps)),
convert.(AbstractMatrix, _tail(A.maps))...
)
end

function Base.Matrix{T}(L::KroneckerSumMap) where {T}
A, B = L.maps
Expand All @@ -156,12 +85,6 @@ function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap)
IB = Diagonal(ones(Bool, size(B, 1)))
return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
end
function SparseArrays.sparse(L::KroneckerSumMap)
A, B = L.maps
IA = sparse(Diagonal(ones(Bool, size(A, 1))))
IB = sparse(Diagonal(ones(Bool, size(B, 1))))
return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
end

# FillMap
Base.Matrix{T}(A::FillMap) where {T} = fill(T(A.λ), size(A))
4 changes: 0 additions & 4 deletions src/linearcombination.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,6 @@ Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::LinearMapTupleO
Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::AbstractVector{<:LinearMap{T}}) where {T} =
LinearCombination{T}(maps)

mean(f::F, maps::LinearMapTupleOrVector) where {F} = sum(f, maps) / length(maps)
mean(maps::LinearMapTupleOrVector) = mean(identity, maps)
mean(A::LinearCombination) = mean(A.maps)

MulStyle(A::LinearCombination) = MulStyle(A.maps...)

# basic methods
Expand Down
2 changes: 0 additions & 2 deletions src/wrappedmap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,14 +24,12 @@ WrappedMap(lmap::MapOrVecOrMat{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kw
# cheap property checks (usually by type)
_issymmetric(A::AbstractMatrix) = false
_issymmetric(A::AbstractQ) = false
_issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
_issymmetric(A::LinearMap) = issymmetric(A)
_issymmetric(A::LinearAlgebra.RealHermSymComplexSym) = issymmetric(A)
_issymmetric(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = issymmetric(A)

_ishermitian(A::AbstractMatrix) = false
_ishermitian(A::AbstractQ) = false
_ishermitian(A::AbstractSparseMatrix) = ishermitian(A)
_ishermitian(A::LinearMap) = ishermitian(A)
_ishermitian(A::LinearAlgebra.RealHermSymComplexHerm) = ishermitian(A)
_ishermitian(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = ishermitian(A)
Expand Down

2 comments on commit b6b1044

@dkarrasch
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/86386

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v3.10.2 -m "<description of version>" b6b10448e0ac041c05581b578248a605cc8025dd
git push origin v3.10.2

Please sign in to comment.