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

Show file tree
Hide file tree
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
2 changes: 2 additions & 0 deletions src/LinearMaps.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -235,6 +235,8 @@ function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α, β
return _generic_mapmat_mul!(y, A, x, α, β)
end

const LinearMapTuple = Tuple{Vararg{LinearMap}}

include("left.jl") # left multiplication by a transpose or adjoint vector
include("transpose.jl") # transposing linear maps
include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
Expand Down
49 changes: 30 additions & 19 deletions src/composition.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,18 +1,18 @@
struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
maps::As # stored in order of application to vector
function CompositeMap{T, As}(maps::As) where {T, As}
N = length(maps)
for n in 2:N
check_dim_mul(maps[n], maps[n-1])
end
for n in eachindex(maps)
A = maps[n]
@assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in CompositeMap constructor"
for TA in Base.Generator(eltype, maps)
promote_type(T, TA) == T ||
error("eltype $TA cannot be promoted to $T in CompositeMap constructor")
end
new{T, As}(maps)
end
end
CompositeMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = CompositeMap{T, As}(maps)
CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(maps)

# basic methods
Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
Expand All @@ -27,7 +27,8 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
$_f(maps::Tuple{}) = true
$_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
$_f(maps::Tuple{Vararg{<:LinearMap}}) = maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
$_f(maps::Tuple{Vararg{<:LinearMap}}) =
maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
# since the introduction of ScaledMap, the following cases cannot occur
# function $_f(maps::Tuple{Vararg{<:LinearMap}}) # length(maps) >= 2
# if maps[1] isa UniformScalingMap{<:RealOrComplex}
Expand Down Expand Up @@ -117,30 +118,40 @@ function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
end

# special transposition behavior
LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(transpose, reverse(A.maps)))
LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjoint, reverse(A.maps)))
LinearAlgebra.transpose(A::CompositeMap{T}) where {T} =
CompositeMap{T}(map(transpose, reverse(A.maps)))
LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} =
CompositeMap{T}(map(adjoint, reverse(A.maps)))

# comparison of CompositeMap objects
Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

# multiplication with vectors
function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, x::AbstractVector)
return _unsafe_mul!(y, A.maps[1], x)
end
function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector)
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)))
end
function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector)
_compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1)))
end
const CompositeMap1 = CompositeMap{<:Any,<:Tuple{LinearMap}}
const CompositeMap2 = CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}
_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap1, x::AbstractVector) =
_unsafe_mul!(y, A.maps[1], x)
# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector)
# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)))
# end
# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector)
# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1)))
# end
_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) =
_compositemul!(y, A, x)

function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector)
function _compositemul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector)
# was there any advantage in having this z an argument, this is not a public method?
z = similar(y, size(A.maps[1], 1))
# no size checking, will be done by individual maps
_unsafe_mul!(z, A.maps[1], x)
_unsafe_mul!(y, A.maps[2], z)
return y
end
function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector)
function _compositemul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector)
# was there any advantage in having these arguments, this is not a public method?
source = similar(y, size(A.maps[1], 1))
dest = similar(y, size(A.maps[2], 1))
# no size checking, will be done by individual maps
N = length(A.maps)
_unsafe_mul!(source, A.maps[1], x)
Expand Down
22 changes: 14 additions & 8 deletions src/functionmap.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,26 @@ function FunctionMap{T}(f::F1, fc::F2, M::Int, N::Int;
end

# additional constructors
FunctionMap{T}(f, M::Int; kwargs...) where {T} = FunctionMap{T}(f, nothing, M, M; kwargs...)
FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} = FunctionMap{T}(f, nothing, M, N; kwargs...)
FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} = FunctionMap{T}(f, fc, M, M; kwargs...)
FunctionMap{T}(f, M::Int; kwargs...) where {T} =
FunctionMap{T}(f, nothing, M, M; kwargs...)
FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} =
FunctionMap{T}(f, nothing, M, N; kwargs...)
FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} =
FunctionMap{T}(f, fc, M, M; kwargs...)

# properties
Base.size(A::FunctionMap) = (A.M, A.N)
Base.parent(A::FunctionMap) = (A.f, A.fc)
Base.parent(A::FunctionMap) = (A.f, A.fc) # is this meanigful/useful?
Copy link
Member

Choose a reason for hiding this comment

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

IIRC the reason for parent was for manual inspection of wrapped maps, because show is now rather sparse. It may be used in the package only rarely, though.

Copy link
Member

Choose a reason for hiding this comment

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

perhaps any parent-like method would be most useful if it returned all the arguments needed to recreate the object. at a minimum that would include size too. maybe a named tuple would be more useful too.
not sure i would use it much though.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

From Base:

help?> Base.parent
  parent(A)

  Return the underlying "parent array”. ...

So I see this making sense for WrappedMap, but other than that I am not entirely sure. Where is this used?

Copy link
Member

Choose a reason for hiding this comment

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

I think we use parent only somewhere in kronecker.jl, and then potentially in an inconsistent way. It's not obvious what parent should return given the wide range of LinearMap subtypes. I'd be fine with actually removing it. I introduced this when I cut down the show method to show much less details, and then I used it in one or two places as it was handy. But this is not really crucial.

Copy link
Member

Choose a reason for hiding this comment

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

Another decision to be made: what shall we do with parent? Just remove? I think it's used internally only in one spot in the kronecker code, and could be worked around with a ternary if-clause, as I had it originally.

LinearAlgebra.issymmetric(A::FunctionMap) = A._issymmetric
LinearAlgebra.ishermitian(A::FunctionMap) = A._ishermitian
LinearAlgebra.isposdef(A::FunctionMap) = A._isposdef
ismutating(A::FunctionMap) = A._ismutating
_ismutating(f) = first(methods(f)).nargs == 3

# multiplication with vector
const TransposeFunctionMap = TransposeMap{<:Any, <:FunctionMap}
const AdjointFunctionMap = AdjointMap{<:Any, <:FunctionMap}

function Base.:(*)(A::FunctionMap, x::AbstractVector)
length(x) == A.N || throw(DimensionMismatch())
if ismutating(A)
Expand All @@ -41,7 +47,7 @@ function Base.:(*)(A::FunctionMap, x::AbstractVector)
end
return y
end
function Base.:(*)(A::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector)
function Base.:(*)(A::AdjointFunctionMap, x::AbstractVector)
Afun = A.lmap
ishermitian(Afun) && return Afun*x
length(x) == size(A, 2) || throw(DimensionMismatch())
Expand All @@ -67,7 +73,7 @@ function Base.:(*)(A::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector)
error("adjoint not implemented for $(A.lmap)")
end
end
function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
function Base.:(*)(A::TransposeFunctionMap, x::AbstractVector)
Afun = A.lmap
(issymmetric(Afun) || (isreal(A) && ishermitian(Afun))) && return Afun*x
length(x) == size(A, 2) || throw(DimensionMismatch())
Expand Down Expand Up @@ -105,7 +111,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, A::FunctionMap, x::AbstractVector)
return y
end

function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeFunctionMap, x::AbstractVector)
A = At.lmap
(issymmetric(A) || (isreal(A) && ishermitian(A))) && return _unsafe_mul!(y, A, x)
if A.fc !== nothing
Expand All @@ -126,7 +132,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeMap{<:Any,<:FunctionMap}
end
end

function _unsafe_mul!(y::AbstractVecOrMat, Ac::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector)
function _unsafe_mul!(y::AbstractVecOrMat, Ac::AdjointFunctionMap, x::AbstractVector)
A = Ac.lmap
ishermitian(A) && return _unsafe_mul!(y, A, x)
if A.fc !== nothing
Expand Down
29 changes: 15 additions & 14 deletions src/linearcombination.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
struct LinearCombination{T, As<:LinearMapTuple} <: LinearMap{T}
maps::As
function LinearCombination{T, As}(maps::As) where {T, As}
N = length(maps)
Expand All @@ -19,9 +19,10 @@ MulStyle(A::LinearCombination) = MulStyle(A.maps...)
# basic methods
Base.size(A::LinearCombination) = size(A.maps[1])
Base.parent(A::LinearCombination) = A.maps
LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) # sufficient but not necessary
LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) # sufficient but not necessary
LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficient but not necessary
# following conditions are sufficient but not necessary
LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps)
LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps)

# adding linear maps
"""
Expand Down Expand Up @@ -59,24 +60,24 @@ end
Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂)

# comparison of LinearCombination objects, sufficient but not necessary
Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps)
Base.:(==)(A::LinearCombination, B::LinearCombination) =
(eltype(A) == eltype(B) && A.maps == B.maps)

# special transposition behavior
LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps))
LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A)}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::LinearCombination) =
LinearCombination{eltype(A)}(map(transpose, A.maps))
LinearAlgebra.adjoint(A::LinearCombination) =
LinearCombination{eltype(A)}(map(adjoint, A.maps))

# multiplication with vectors & matrices
for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
@eval begin
function _unsafe_mul!(y::$outtype, A::LinearCombination, x::$intype)

function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In)
_unsafe_mul!(y, first(A.maps), x)
_mul!(MulStyle(A), y, A, x, true)
return y
end

function _unsafe_mul!(y::$outtype, A::LinearCombination, x::$intype,
α::Number, β::Number)
function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In, α::Number, β::Number)
if iszero(α) # trivial cases
iszero(β) && return fill!(y, zero(eltype(y)))
isone(β) && return y
Expand Down Expand Up @@ -109,7 +110,7 @@ function _mul!(::ThreeArg, y, A::LinearCombination, x, α)
return __mul!(y, Base.tail(A.maps), x, α, similar(y))
end

__mul!(y, As::Tuple{Vararg{LinearMap}}, x, α, z) =
__mul!(y, As::LinearMapTuple, x, α, z) =
__mul!(__mul!(y, first(As), x, α, z), Base.tail(As), x, α, z)
__mul!(y, A::Tuple{LinearMap}, x, α, z) = __mul!(y, first(A), x, α, z)
__mul!(y, A::LinearMap, x, α, z) = muladd!(MulStyle(A), y, A, x, α, z)
Expand Down
4 changes: 2 additions & 2 deletions src/scaledmap.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} =
Base.size(A::ScaledMap) = size(A.lmap)
Base.parent(A::ScaledMap) = (A.λ, A.lmap)
Base.isreal(A::ScaledMap) = isreal(A.λ) && isreal(A.lmap)
LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap)
LinearAlgebra.ishermitian(A::ScaledMap) = ishermitian(A.lmap)
LinearAlgebra.issymmetric(A::ScaledMap) = isreal(A.λ) && issymmetric(A.lmap)
LinearAlgebra.ishermitian(A::ScaledMap) = isreal(A.λ) && ishermitian(A.lmap)
LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap)

Base.transpose(A::ScaledMap) = A.λ * transpose(A.lmap)
Expand Down
2 changes: 1 addition & 1 deletion src/show.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ function _show_typeof(A::LinearMap{T}) where {T}
split(string(typeof(A)), '{')[1] * '{' * string(T) * '}'
end

function print_maps(io::IO, maps::Tuple{Vararg{LinearMap}}, k)
function print_maps(io::IO, maps::LinearMapTuple, k)
n = length(maps)
str = ""
if get(io, :limit, true) && n > 10
Expand Down
1 change: 1 addition & 0 deletions src/transpose.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,7 @@ _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β
const ConjugateMap = AdjointMap{<:Any, <:TransposeMap}
# canonical order of adjoint followed by transpose
LinearAlgebra.transpose(A::AdjointMap) = adjoint(transpose(A.lmap))
LinearAlgebra.transpose(A::ConjugateMap) = adjoint(A.lmap.lmap)
for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
@eval begin
function _unsafe_mul!(y::$Out, Ac::ConjugateMap, x::$In)
Expand Down
2 changes: 1 addition & 1 deletion test/composition.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
FCM = LinearMaps.CompositeMap{ComplexF64}((FC,))
L = LowerTriangular(ones(10,10))
@test_throws DimensionMismatch F * LinearMap(rand(2,2))
@test_throws AssertionError LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10))))
@test_throws ErrorException LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10))))
A = 2 * rand(ComplexF64, (10, 10)) .- 1
B = rand(size(A)...)
H = LinearMap(Hermitian(A'A))
Expand Down