Skip to content

WIP: review before version 3 #119

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
Dec 9, 2020
Merged

WIP: review before version 3 #119

Changes from 1 commit
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
2c58a40
review part 1: some spaces and linebreaks, in/outtype -> In/Out for b…
Jutho Dec 5, 2020
b9036aa
review part 2
Jutho Dec 5, 2020
04b923d
review part 2b: new compositemul attempt
Jutho Dec 6, 2020
59312cd
typo (don't push without testing locally :D )
Jutho Dec 6, 2020
e7bb9f5
review part 2c: add comment
Jutho Dec 6, 2020
b1208fb
review part 3: blockmap
Jutho Dec 6, 2020
f000f83
review part 3b: change comments
Jutho Dec 7, 2020
ba51f26
review part 4:kronecker
Jutho Dec 7, 2020
7d3300b
review part 4b: recursive kron
Jutho Dec 8, 2020
122baed
review part 5: conversion spaces and questions
Jutho Dec 8, 2020
cb63bbc
review part 6: show, linebreak for 92-character rule
Jutho Dec 8, 2020
025cd53
review part 6b: remove comments
Jutho Dec 9, 2020
94bab04
clean up final comments and suggestions
Jutho Dec 9, 2020
10b2531
remove parent, remove a few comments
dkarrasch Dec 9, 2020
85e350e
improve coverage, minor kronecker fix
dkarrasch Dec 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 38 additions & 16 deletions src/conversion.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,10 @@ SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A)
# special cases

# ScaledMap
Base.Matrix{T}(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) where {T} = convert(Matrix{T}, A.λ*A.lmap.lmap)
SparseArrays.sparse(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) = convert(SparseMatrixCSC, A.λ*A.lmap.lmap)
Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) where {T} =
convert(Matrix{T}, A.λ * A.lmap.lmap)
SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) =
convert(SparseMatrixCSC, A.λ * A.lmap.lmap)

# UniformScalingMap
Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J))
Expand All @@ -62,25 +64,29 @@ SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)

# TransposeMap & AdjointMap
for (TT, T) in ((AdjointMap, adjoint), (TransposeMap, transpose))
@eval Base.convert(::Type{AbstractMatrix}, A::$TT) = $T(convert(AbstractMatrix, A.lmap))
@eval SparseArrays.sparse(A::$TT) = $T(convert(SparseMatrixCSC, A.lmap))
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 Base.Matrix{T}(ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}}) where {T}
function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) where {T}
maps = ΣA.maps
mats = map(A->getfield(A, :lmap), maps)
return Matrix{T}(sum(mats))
end
function SparseArrays.sparse(ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}})
function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}})
maps = ΣA.maps
mats = map(A->getfield(A, :lmap), maps)
return convert(SparseMatrixCSC, sum(mats))
end

# CompositeMap
function Base.Matrix{T}(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,LinearMap}}) where {T}
# Are all these methods below really useful?
# The conversion to matrices is not something that belongs to ordinary use of LinearMaps.
# It's rather for debugging purposes, where efficiency is not the primary concern.
# Also, I don't think they are really that much more efficient
Copy link
Member

Choose a reason for hiding this comment

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

Generally, the Matrix conversion is used in the Kronecker code in the vec trick, which is why I started introducing these methods (at least with the intention to write code that is faster than the columnwise mat-vec-muls):

https://github.com/Jutho/LinearMaps.jl/blob/b38c1ac1900d4751237bbef62927477017f4f7d0/src/kronecker.jl#L121-L131

X there is promoted to a MatrixMap, that's why we have all those MatrixMap conversion methods here.

function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T}
B, A = AB.maps
require_one_based_indexing(B)
Y = Matrix{eltype(AB)}(undef, size(AB))
Expand All @@ -92,39 +98,48 @@ end
for ((TA, fieldA), (TB, fieldB)) in (((MatrixMap, :lmap), (MatrixMap, :lmap)),
((MatrixMap, :lmap), (UniformScalingMap, :λ)),
((UniformScalingMap, :λ), (MatrixMap, :lmap)))
@eval function Base.convert(::Type{AbstractMatrix}, AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}})
@eval function Base.convert(::Type{AbstractMatrix},
AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}})
B, A = AB.maps
return A.$fieldA*B.$fieldB
end
end
function Base.Matrix{T}(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}}) where {T}
function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}}) where {T}
B, A = AB.maps
return convert(Matrix{T}, A.lmap*B.lmap)
end
function SparseArrays.sparse(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}})
function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}})
B, A = AB.maps
return convert(SparseMatrixCSC, A.lmap*B.lmap)
end
function Base.Matrix{T}(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}}) where {T}
function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) where {T}
A, J = λA.maps
return convert(Matrix{T}, J.λ*A.lmap)
# if you want less allocations: lmul!(J.λ, convert(Matrix{T}, A.lmap))
end
function SparseArrays.sparse(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}})
A, J = λA.maps
return convert(SparseMatrixCSC, J.λ*A.lmap)
# probably also works with sparse: lmul!(J.λ, convert(SparseMatrixCSC, A.lmap))

end
function Base.Matrix{T}(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}}) where {T}
function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T}
J, A = Aλ.maps
return convert(Matrix{T}, A.lmap*J.λ)
# if you want less allocations: rmul!(convert(Matrix{T}, A.lmap), J.λ)
end
function SparseArrays.sparse(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}})
J, A = Aλ.maps
return convert(SparseMatrixCSC, A.lmap*J.λ)
# probably also works with sparse: rmul!(convert(SparseMatrixCSC, A.lmap), J.λ)
end

# BlockMap & BlockDiagonalMap
Base.Matrix{T}(A::BlockMap) where {T} = hvcat(A.rows, convert.(Matrix{T}, A.maps)...)
Base.convert(::Type{AbstractMatrix}, A::BlockMap) = hvcat(A.rows, convert.(AbstractMatrix, A.maps)...)
# seems like code duplication:
# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A)

function SparseArrays.sparse(A::BlockMap)
return hvcat(
A.rows,
Expand All @@ -140,13 +155,17 @@ 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)...)
Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) =
kron(convert.(AbstractMatrix, A.maps)...)
# seems like code duplication:
# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A)
function SparseArrays.sparse(A::KroneckerMap)
return kron(
convert(SparseMatrixCSC, first(A.maps)),
convert.(AbstractMatrix, Base.tail(A.maps))...
)
end

function Base.Matrix{T}(L::KroneckerSumMap) where {T}
A, B = L.maps
IA = Diagonal(ones(Bool, size(A, 1)))
Expand All @@ -159,6 +178,9 @@ 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
# seems like code duplication:
# convert(AbstractMatrix, L) will anyway fall back to Matrix{eltype(L)}(L)

function SparseArrays.sparse(L::KroneckerSumMap)
A, B = L.maps
IA = sparse(Diagonal(ones(Bool, size(A, 1))))
Expand Down