From 6342d2f9aca17a8219f0e5ddda7aa48ec001216a Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Tue, 21 Sep 2021 11:56:23 -0400 Subject: [PATCH 01/11] bump patch --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6c4f7fb0..71f60af6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearMaps" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "3.4.0" +version = "3.4.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From bd3dad50539122e04b08eeb2d3245726ed808b71 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 1 Dec 2022 18:11:02 -0500 Subject: [PATCH 02/11] Towards handling Unitful linear maps --- src/composition.jl | 78 +++++++++++++++++++++++++++++++-------------- test/composition.jl | 2 +- 2 files changed, 55 insertions(+), 25 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index afb20148..67088d86 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -5,15 +5,25 @@ struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T} for n in 2:N check_dim_mul(maps[n], maps[n-1]) end - for TA in Base.Iterators.map(eltype, maps) - promote_type(T, TA) == T || - error("eltype $TA cannot be promoted to $T in CompositeMap constructor") - end + Tprod = eltype(prod(m -> oneunit(eltype(m)), maps)) +# todo: why force caller to provide Tprod? just figure out Tprod in inner constructor? + T == Tprod || @warn("eltype $Tprod != $T in CompositeMap constructor") +# for TA in Base.Iterators.map(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<:LinearMapTupleOrVector} = CompositeMap{T, As}(maps) +# constructor with eltype inferred from the product +function CompositeMap(maps::As) where {As <: LinearMapTupleOrVector} + T = eltype(prod(m -> oneunit(eltype(m)), maps)) + return CompositeMap{T, As}(maps) +end + Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTupleOrVector) = CompositeMap{promote_type(map(eltype, maps)...)}(reverse(maps)) Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector{<:LinearMap{T}}) where {T} = @@ -78,39 +88,51 @@ end # scalar multiplication and division (non-commutative case) function Base.:(*)(α::Number, A::LinearMap) - T = promote_type(typeof(α), eltype(A)) - return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) +# T = promote_type(typeof(α), eltype(A)) +# T = eltype(oneunit(α) * oneunit(eltype(A))) +# return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) + return CompositeMap(_combine(A, UniformScalingMap(α, size(A, 1)))) end function Base.:(*)(α::Number, A::CompositeMap) - T = promote_type(typeof(α), eltype(A)) +# T = promote_type(typeof(α), eltype(A)) +# T = eltype(oneunit(α) * oneunit(eltype(A))) Alast = last(A.maps) if Alast isa UniformScalingMap - return CompositeMap{T}(_combine(_front(A.maps), α * Alast)) +# return CompositeMap{T}(_combine(_front(A.maps), α * Alast)) + return CompositeMap(_combine(_front(A.maps), α * Alast)) else - return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) +# return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) + return CompositeMap(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) end end # needed for disambiguation function Base.:(*)(α::RealOrComplex, A::CompositeMap{<:RealOrComplex}) - T = Base.promote_op(*, typeof(α), eltype(A)) +# T = Base.promote_op(*, typeof(α), eltype(A)) + T = eltype(oneunit(α) * oneunit(eltype(A))) return ScaledMap{T}(α, A) end function Base.:(*)(A::LinearMap, α::Number) - T = promote_type(typeof(α), eltype(A)) - return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A)) +# T = promote_type(typeof(α), eltype(A)) +# T = eltype(oneunit(eltype(A)) * oneunit(α)) +# return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A)) + return CompositeMap(_combine(UniformScalingMap(α, size(A, 2)), A)) end function Base.:(*)(A::CompositeMap, α::Number) - T = promote_type(typeof(α), eltype(A)) +# T = promote_type(typeof(α), eltype(A)) +# T = eltype(oneunit(eltype(A)) * oneunit(α)) Afirst = first(A.maps) if Afirst isa UniformScalingMap - return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps))) +# return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps))) + return CompositeMap(_combine(Afirst * α, _tail(A.maps))) else - return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) +# return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) + return CompositeMap(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) end end # needed for disambiguation function Base.:(*)(A::CompositeMap{<:RealOrComplex}, α::RealOrComplex) - T = Base.promote_op(*, typeof(α), eltype(A)) +# T = Base.promote_op(*, typeof(α), eltype(A)) + T = eltype(oneunit(eltype(A)) * oneunit(α)) return ScaledMap{T}(α, A) end @@ -135,20 +157,28 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3); ``` """ function Base.:(*)(A₁::LinearMap, A₂::LinearMap) - T = promote_type(eltype(A₁), eltype(A₂)) - return CompositeMap{T}(_combine(A₂, A₁)) +# T = promote_type(eltype(A₁), eltype(A₂)) +# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) +# return CompositeMap{T}(_combine(A₂, A₁)) + return CompositeMap(_combine(A₂, A₁)) end function Base.:(*)(A₁::LinearMap, A₂::CompositeMap) - T = promote_type(eltype(A₁), eltype(A₂)) - return CompositeMap{T}(_combine(A₂.maps, A₁)) +# T = promote_type(eltype(A₁), eltype(A₂)) +# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) +# return CompositeMap{T}(_combine(A₂.maps, A₁)) + return CompositeMap(_combine(A₂.maps, A₁)) end function Base.:(*)(A₁::CompositeMap, A₂::LinearMap) - T = promote_type(eltype(A₁), eltype(A₂)) - return CompositeMap{T}(_combine(A₂, A₁.maps)) +# T = promote_type(eltype(A₁), eltype(A₂)) +# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) +# return CompositeMap{T}(_combine(A₂, A₁.maps)) + return CompositeMap(_combine(A₂, A₁.maps)) end function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap) - T = promote_type(eltype(A₁), eltype(A₂)) - return CompositeMap{T}(_combine(A₂.maps, A₁.maps)) +# T = promote_type(eltype(A₁), eltype(A₂)) +# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) +# return CompositeMap{T}(_combine(A₂.maps, A₁.maps)) + return CompositeMap(_combine(A₂.maps, A₁.maps)) end # needed for disambiguation Base.:(*)(A₁::ScaledMap, A₂::CompositeMap) = A₁.λ * (A₁.lmap * A₂) diff --git a/test/composition.jl b/test/composition.jl index a6076165..c45f6259 100644 --- a/test/composition.jl +++ b/test/composition.jl @@ -7,7 +7,7 @@ using LinearMaps: LinearMapVector, LinearMapTuple FCM = LinearMaps.CompositeMap{ComplexF64}((FC,)) L = LowerTriangular(ones(10,10)) @test_throws DimensionMismatch F * LinearMap(rand(2,2)) - @test_throws ErrorException 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)) From 840558ddfe84c6c822f6ae6a193addf063ea9fe8 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 1 Dec 2022 18:29:47 -0500 Subject: [PATCH 03/11] Add .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..ba39cc53 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +Manifest.toml From 4392a92ec4ba27c0ebcf23fb4889df00749f60e0 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 1 Dec 2022 18:30:11 -0500 Subject: [PATCH 04/11] Revert throws --- src/composition.jl | 10 +++++----- test/composition.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index 67088d86..21a76388 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -5,10 +5,10 @@ struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T} for n in 2:N check_dim_mul(maps[n], maps[n-1]) end - Tprod = eltype(prod(m -> oneunit(eltype(m)), maps)) -# todo: why force caller to provide Tprod? just figure out Tprod in inner constructor? - T == Tprod || @warn("eltype $Tprod != $T in CompositeMap constructor") -# for TA in Base.Iterators.map(eltype, maps) + Tprod = eltype(prod(m -> oneunit(eltype(m)), maps)) # handles units + promote_type(T, Tprod) == T || + error("eltype $Tprod and $T incompatible in CompositeMap constructor") +# for TA in Base.Iterators.map(eltype, maps) # todo: cut # promote_type(T, TA) == T || # error("eltype $TA cannot be promoted to $T in CompositeMap constructor") # end @@ -20,7 +20,7 @@ CompositeMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} = CompositeMap{T # constructor with eltype inferred from the product function CompositeMap(maps::As) where {As <: LinearMapTupleOrVector} - T = eltype(prod(m -> oneunit(eltype(m)), maps)) + T = eltype(prod(m -> oneunit(eltype(m)), maps)) # todo: can the compiler infer? return CompositeMap{T, As}(maps) end diff --git a/test/composition.jl b/test/composition.jl index c45f6259..a6076165 100644 --- a/test/composition.jl +++ b/test/composition.jl @@ -7,7 +7,7 @@ using LinearMaps: LinearMapVector, LinearMapTuple FCM = LinearMaps.CompositeMap{ComplexF64}((FC,)) L = LowerTriangular(ones(10,10)) @test_throws DimensionMismatch F * LinearMap(rand(2,2)) -# @test_throws ErrorException 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)) From 84aef7741727497a0f2f6d9ca63941e4bc4a5335 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 09:00:12 -0500 Subject: [PATCH 05/11] Apply suggestions from code review `typeof` Co-authored-by: Daniel Karrasch --- src/composition.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index 21a76388..663ccffb 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -5,7 +5,7 @@ struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T} for n in 2:N check_dim_mul(maps[n], maps[n-1]) end - Tprod = eltype(prod(m -> oneunit(eltype(m)), maps)) # handles units + Tprod = typeof(*(map(oneunit∘eltype, maps)...)) # handles units promote_type(T, Tprod) == T || error("eltype $Tprod and $T incompatible in CompositeMap constructor") # for TA in Base.Iterators.map(eltype, maps) # todo: cut @@ -88,10 +88,8 @@ end # scalar multiplication and division (non-commutative case) function Base.:(*)(α::Number, A::LinearMap) -# T = promote_type(typeof(α), eltype(A)) -# T = eltype(oneunit(α) * oneunit(eltype(A))) -# return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) - return CompositeMap(_combine(A, UniformScalingMap(α, size(A, 1)))) + T = typeof(oneunit(α) * oneunit(eltype(A))) + return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) end function Base.:(*)(α::Number, A::CompositeMap) # T = promote_type(typeof(α), eltype(A)) From 9c9c3d106478a7140ab6ada7d22c827464fe076e Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 09:36:29 -0500 Subject: [PATCH 06/11] Use typeof --- src/composition.jl | 69 +++++++++++++++++----------------------------- 1 file changed, 26 insertions(+), 43 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index 663ccffb..c15be32a 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -8,24 +8,25 @@ struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T} Tprod = typeof(*(map(oneunit∘eltype, maps)...)) # handles units promote_type(T, Tprod) == T || error("eltype $Tprod and $T incompatible in CompositeMap constructor") -# for TA in Base.Iterators.map(eltype, maps) # todo: cut -# 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<:LinearMapTupleOrVector} = CompositeMap{T, As}(maps) +#= # constructor with eltype inferred from the product +# todo - cut? function CompositeMap(maps::As) where {As <: LinearMapTupleOrVector} - T = eltype(prod(m -> oneunit(eltype(m)), maps)) # todo: can the compiler infer? + T = typeof(*(map(oneunit∘eltype, maps)...)) # todo: can the compiler infer? return CompositeMap{T, As}(maps) end +=# -Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTupleOrVector) = - CompositeMap{promote_type(map(eltype, maps)...)}(reverse(maps)) +function Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTupleOrVector) + Tprod = typeof(*(map(oneunit∘eltype, maps)...)) # handles units + return CompositeMap{Tprod}(reverse(maps)) +end Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector{<:LinearMap{T}}) where {T} = CompositeMap{T}(reverse(maps)) @@ -92,45 +93,35 @@ function Base.:(*)(α::Number, A::LinearMap) return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) end function Base.:(*)(α::Number, A::CompositeMap) -# T = promote_type(typeof(α), eltype(A)) -# T = eltype(oneunit(α) * oneunit(eltype(A))) + T = typeof(oneunit(α) * oneunit(eltype(A))) Alast = last(A.maps) if Alast isa UniformScalingMap -# return CompositeMap{T}(_combine(_front(A.maps), α * Alast)) - return CompositeMap(_combine(_front(A.maps), α * Alast)) + return CompositeMap{T}(_combine(_front(A.maps), α * Alast)) else -# return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) - return CompositeMap(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) + return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1)))) end end # needed for disambiguation function Base.:(*)(α::RealOrComplex, A::CompositeMap{<:RealOrComplex}) -# T = Base.promote_op(*, typeof(α), eltype(A)) - T = eltype(oneunit(α) * oneunit(eltype(A))) + T = typeof(oneunit(α) * oneunit(eltype(A))) return ScaledMap{T}(α, A) end function Base.:(*)(A::LinearMap, α::Number) -# T = promote_type(typeof(α), eltype(A)) -# T = eltype(oneunit(eltype(A)) * oneunit(α)) -# return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A)) - return CompositeMap(_combine(UniformScalingMap(α, size(A, 2)), A)) + T = typeof(oneunit(eltype(A)) * oneunit(α)) + return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A)) end function Base.:(*)(A::CompositeMap, α::Number) -# T = promote_type(typeof(α), eltype(A)) -# T = eltype(oneunit(eltype(A)) * oneunit(α)) + T = typeof(oneunit(eltype(A)) * oneunit(α)) Afirst = first(A.maps) if Afirst isa UniformScalingMap -# return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps))) - return CompositeMap(_combine(Afirst * α, _tail(A.maps))) + return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps))) else -# return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) - return CompositeMap(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) + return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps)) end end # needed for disambiguation function Base.:(*)(A::CompositeMap{<:RealOrComplex}, α::RealOrComplex) -# T = Base.promote_op(*, typeof(α), eltype(A)) - T = eltype(oneunit(eltype(A)) * oneunit(α)) + T = typeof(oneunit(eltype(A)) * oneunit(α)) return ScaledMap{T}(α, A) end @@ -155,28 +146,20 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3); ``` """ function Base.:(*)(A₁::LinearMap, A₂::LinearMap) -# T = promote_type(eltype(A₁), eltype(A₂)) -# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) -# return CompositeMap{T}(_combine(A₂, A₁)) - return CompositeMap(_combine(A₂, A₁)) + T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + return CompositeMap{T}(_combine(A₂, A₁)) end function Base.:(*)(A₁::LinearMap, A₂::CompositeMap) -# T = promote_type(eltype(A₁), eltype(A₂)) -# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) -# return CompositeMap{T}(_combine(A₂.maps, A₁)) - return CompositeMap(_combine(A₂.maps, A₁)) + T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + return CompositeMap{T}(_combine(A₂.maps, A₁)) end function Base.:(*)(A₁::CompositeMap, A₂::LinearMap) -# T = promote_type(eltype(A₁), eltype(A₂)) -# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) -# return CompositeMap{T}(_combine(A₂, A₁.maps)) - return CompositeMap(_combine(A₂, A₁.maps)) + T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + return CompositeMap{T}(_combine(A₂, A₁.maps)) end function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap) -# T = promote_type(eltype(A₁), eltype(A₂)) -# T = eltype(prod(m -> oneunit(eltype(m)), (A₁, A₂))) -# return CompositeMap{T}(_combine(A₂.maps, A₁.maps)) - return CompositeMap(_combine(A₂.maps, A₁.maps)) + T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + return CompositeMap{T}(_combine(A₂.maps, A₁.maps)) end # needed for disambiguation Base.:(*)(A₁::ScaledMap, A₂::CompositeMap) = A₁.λ * (A₁.lmap * A₂) From 2547151d32880bf298b0809094ce0ac8a661901b Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 09:38:09 -0500 Subject: [PATCH 07/11] Add tests --- test/runtests.jl | 2 ++ test/units.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 test/units.jl diff --git a/test/runtests.jl b/test/runtests.jl index 9ad9bb06..87d5d171 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,3 +41,5 @@ include("getindex.jl") include("inversemap.jl") include("rrules.jl") + +include("units.jl") diff --git a/test/units.jl b/test/units.jl new file mode 100644 index 00000000..e4db202e --- /dev/null +++ b/test/units.jl @@ -0,0 +1,26 @@ +# test/units + +using Test: @test, @testset, @inferred, @test_throws +using LinearMaps: LinearMap +using Unitful: m, s, g + +@testset "units" begin + A = rand(4,3) * 1m + B = rand(3,2) * 1s + C = A * B + D = 1f0g * C + + Ma = @inferred LinearMap(A) + Mb = @inferred LinearMap(B) + Mc = Ma * Mb + Md = 1f0g * Mc # todo - UniformScalingMap + + @test Matrix(Ma) == A + @test Matrix(Mc) == C + @test_throws ArgumentError Matrix(Md) == D + + x = randn(2) + @test B * x == Mb * x + @test_broken C * x == Mc * x + @test_broken D * x == Md * x +end From 2e9c4040e0b44d77a23b4d55cc69e1614a255aff Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 09:39:46 -0500 Subject: [PATCH 08/11] Cut unneeded constructor --- src/composition.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index c15be32a..5ff36947 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -14,15 +14,6 @@ end CompositeMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} = CompositeMap{T, As}(maps) -#= -# constructor with eltype inferred from the product -# todo - cut? -function CompositeMap(maps::As) where {As <: LinearMapTupleOrVector} - T = typeof(*(map(oneunit∘eltype, maps)...)) # todo: can the compiler infer? - return CompositeMap{T, As}(maps) -end -=# - function Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTupleOrVector) Tprod = typeof(*(map(oneunit∘eltype, maps)...)) # handles units return CompositeMap{Tprod}(reverse(maps)) From f34841eeb455b66df39200dae7bb44667f92aaf1 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 15:23:09 -0500 Subject: [PATCH 09/11] ScaledMap constructor --- src/scaledmap.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index 549bf7a8..d6f82eb9 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -7,14 +7,16 @@ struct ScaledMap{T, S<:RealOrComplex, L<:LinearMap} <: LinearMap{T} λ::S lmap::L function ScaledMap{T}(λ::S, A::L) where {T, S <: RealOrComplex, L <: LinearMap{<:RealOrComplex}} - @assert Base.promote_op(*, S, eltype(A)) == T "target type $T cannot hold products of $S and $(eltype(A)) objects" + Tprod = typeof(oneunit(S) * oneunit(eltype(A))) + promote_type(T, Tprod) == T || + error("target type $T vs product of $S and $(eltype(A))") new{T,S,L}(λ, A) end end # constructor ScaledMap(λ::RealOrComplex, lmap::LinearMap{<:RealOrComplex}) = - ScaledMap{Base.promote_op(*, typeof(λ), eltype(lmap))}(λ, lmap) + ScaledMap{typeof(oneunit(λ) * oneunit(eltype(lmap)))}(λ, lmap) # basic methods Base.size(A::ScaledMap) = size(A.lmap) From 42a8100d8add9176b80460573ec959dd5473db1d Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 15:23:50 -0500 Subject: [PATCH 10/11] Multiplication --- src/LinearMaps.jl | 16 ++++++------- src/composition.jl | 58 ++++++++++++++++++++++++++++++++-------------- test/scaledmap.jl | 2 +- test/units.jl | 8 +++---- 4 files changed, 52 insertions(+), 32 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 80d09f8a..2923ea1f 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -133,7 +133,7 @@ julia> A(x) """ function Base.:(*)(A::LinearMap, x::AbstractVector) check_dim_mul(A, x) - T = promote_type(eltype(A), eltype(x)) + T = typeof(oneunit(eltype(A)) * oneunit(eltype(x))) y = similar(x, T, axes(A)[1]) return mul!(y, A, x) end @@ -311,25 +311,23 @@ function _generic_map_mul!(Y, A, X::AbstractMatrix, α, β) end return Y end -function _generic_map_mul!(Y, A, s::Number) - T = promote_type(eltype(A), typeof(s)) +function _generic_map_mul!(Y, A, s::S) where {S <: Number} ax2 = axes(A)[2] - xi = zeros(T, ax2) + xi = zeros(S, ax2) @inbounds for (i, Yi) in zip(ax2, eachcol(Y)) xi[i] = s mul!(Yi, A, xi) - xi[i] = zero(T) + xi[i] = zero(S) end return Y end -function _generic_map_mul!(Y, A, s::Number, α, β) - T = promote_type(eltype(A), typeof(s)) +function _generic_map_mul!(Y, A, s::S, α, β) where {S <: Number} ax2 = axes(A)[2] - xi = zeros(T, ax2) + xi = zeros(S, ax2) @inbounds for (i, Yi) in zip(ax2, eachcol(Y)) xi[i] = s mul!(Yi, A, xi, α, β) - xi[i] = zero(T) + xi[i] = zero(S) end return Y end diff --git a/src/composition.jl b/src/composition.jl index 5ff36947..5f30822b 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -1,3 +1,9 @@ +# appropriate type for product of two (possibly Unitful) quantities: +_multype(a, b) = typeof(oneunit(eltype(a)) * oneunit(eltype(b))) +# this variant is needed because reinterpret changes length of complex vectors +#_multype(a, b, iscomplex::Bool) = +# typeof((1+0im) * oneunit(eltype(a)) * oneunit(eltype(b))) + struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T} maps::As # stored in order of application to vector function CompositeMap{T, As}(maps::As) where {T, As} @@ -80,11 +86,11 @@ end # scalar multiplication and division (non-commutative case) function Base.:(*)(α::Number, A::LinearMap) - T = typeof(oneunit(α) * oneunit(eltype(A))) + T = _multype(α, A) return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1)))) end function Base.:(*)(α::Number, A::CompositeMap) - T = typeof(oneunit(α) * oneunit(eltype(A))) + T = _multype(α, A) Alast = last(A.maps) if Alast isa UniformScalingMap return CompositeMap{T}(_combine(_front(A.maps), α * Alast)) @@ -94,15 +100,15 @@ function Base.:(*)(α::Number, A::CompositeMap) end # needed for disambiguation function Base.:(*)(α::RealOrComplex, A::CompositeMap{<:RealOrComplex}) - T = typeof(oneunit(α) * oneunit(eltype(A))) + T = _multype(α, A) return ScaledMap{T}(α, A) end function Base.:(*)(A::LinearMap, α::Number) - T = typeof(oneunit(eltype(A)) * oneunit(α)) + T = _multype(A, α) return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A)) end function Base.:(*)(A::CompositeMap, α::Number) - T = typeof(oneunit(eltype(A)) * oneunit(α)) + T = _multype(A, α) Afirst = first(A.maps) if Afirst isa UniformScalingMap return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps))) @@ -112,7 +118,7 @@ function Base.:(*)(A::CompositeMap, α::Number) end # needed for disambiguation function Base.:(*)(A::CompositeMap{<:RealOrComplex}, α::RealOrComplex) - T = typeof(oneunit(eltype(A)) * oneunit(α)) + T = _multype(A, α) return ScaledMap{T}(α, A) end @@ -137,19 +143,19 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3); ``` """ function Base.:(*)(A₁::LinearMap, A₂::LinearMap) - T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + T = _multype(A₁, A₂) return CompositeMap{T}(_combine(A₂, A₁)) end function Base.:(*)(A₁::LinearMap, A₂::CompositeMap) - T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + T = _multype(A₁, A₂) return CompositeMap{T}(_combine(A₂.maps, A₁)) end function Base.:(*)(A₁::CompositeMap, A₂::LinearMap) - T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + T = _multype(A₁, A₂) return CompositeMap{T}(_combine(A₂, A₁.maps)) end function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap) - T = typeof(oneunit(eltype(A₁)) * oneunit(eltype(A₂))) + T = _multype(A₁, A₂) return CompositeMap{T}(_combine(A₂.maps, A₁.maps)) end # needed for disambiguation @@ -175,9 +181,13 @@ function _compositemul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, x, dest = nothing) return _unsafe_mul!(y, A.maps[1], x) end -function _compositemul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x, - source = similar(y, (size(A.maps[1],1), size(x)[2:end]...)), - dest = nothing) +function _compositemul!( + y, + A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, + x, + source = similar(y, _multype(A.maps[1], x), (size(A.maps[1],1), size(x)[2:end]...)), + dest = nothing, +) _unsafe_mul!(source, A.maps[1], x) _unsafe_mul!(y, A.maps[2], source) return y @@ -200,15 +210,27 @@ function _resize(dest::AbstractMatrix, sz::Tuple{<:Integer,<:Integer}) similar(dest, sz) end -function _compositemul!(y, A::CompositeMap{<:Any,<:LinearMapTuple}, x, - source = similar(y, (size(A.maps[1],1), size(x)[2:end]...)), - dest = similar(y, (size(A.maps[2],1), size(x)[2:end]...))) +function _compositemul!( + y, + A::CompositeMap{<:Any,<:LinearMapTuple}, + x, + source = similar([], _multype(A.maps[1], x), (size(A.maps[1],1), size(x)[2:end]...)), + dest = nothing, # similar(y, (size(A.maps[2],1), size(x)[2:end]...)), +) N = length(A.maps) + # caution: be careful if y is complex but intermediate products are not +# source = reinterpret(_multype(A.maps[1], x, !isreal(y)), source) # for units +# resize!(source, size(A.maps[1],1)) # trick due to complex case +# todo: build reinterpret into _unsafe_mul! instead? +# only necessary if either source or map has Number type instead of Real|Complex _unsafe_mul!(source, A.maps[1], x) for n in 2:N-1 - dest = _resize(dest, (size(A.maps[n],1), size(x)[2:end]...)) +# dest = _resize(dest, (size(A.maps[n],1), size(x)[2:end]...)) +# dest = reinterpret(_multype(A.maps[n], source), dest) + dest = similar([], _multype(A.maps[n], source), (size(A.maps[n],1), size(source)[2:end]...)) _unsafe_mul!(dest, A.maps[n], source) - dest, source = source, dest # alternate dest and source +# dest, source = source, dest # alternate dest and source + dest, source = [], dest # reuse source as next dest, but abandon dest end _unsafe_mul!(y, A.maps[N], source) return y diff --git a/test/scaledmap.jl b/test/scaledmap.jl index d611dd0d..330371e1 100644 --- a/test/scaledmap.jl +++ b/test/scaledmap.jl @@ -40,7 +40,7 @@ using Test, LinearMaps, LinearAlgebra # complex case β = π + 2π * im C = @inferred β * A - @test_throws AssertionError LinearMaps.ScaledMap{Float64}(β, A) + @test_throws ErrorException LinearMaps.ScaledMap{Float64}(β, A) @inferred conj(β) * A' # needed in left-mul T = ComplexF64 xc = rand(T, N) diff --git a/test/units.jl b/test/units.jl index e4db202e..f2c8119f 100644 --- a/test/units.jl +++ b/test/units.jl @@ -13,14 +13,14 @@ using Unitful: m, s, g Ma = @inferred LinearMap(A) Mb = @inferred LinearMap(B) Mc = Ma * Mb - Md = 1f0g * Mc # todo - UniformScalingMap + Md = 1f0g * Mc @test Matrix(Ma) == A @test Matrix(Mc) == C - @test_throws ArgumentError Matrix(Md) == D + @test Matrix(Md) == D x = randn(2) @test B * x == Mb * x - @test_broken C * x == Mc * x - @test_broken D * x == Md * x + @test C * x ≈ Mc * x + @test D * x ≈ Md * x end From c97726803a9084f604e0a8cffdcb6ca48200a13a Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 2 Dec 2022 15:32:00 -0500 Subject: [PATCH 11/11] Add Unitful to test/ --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 020e0bbc..3453a9b1 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,6 +10,7 @@ Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] Aqua = "0.5"