make mean work stably

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.5.2"
+version = "3.6.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/history.md

@@ -1,5 +1,15 @@
 # Version history
 
+## What's new in v3.6
+
+* Support for Julia versions below v1.6 has been dropped.
+* `Block[Diagonal]Map`, `CompositeMap`, `KroneckerMap` and `LinearCombination` type objects
+  can now be backed by a `Vector` of `LinearMap`-type elements. This can be beneficial in
+  cases where these higher-order `LinearMap`s are constructed from many maps where a tuple
+  backend may get inefficient or impose hard work for the compiler at construction.
+  The default behavior, however, does not change, and construction of vector-based
+  `LinearMap`s requires usage of the unexported constructors ("expert usage").
+
 ## What's new in v3.5
 
 * `WrappedMap`, `ScaledMap`, `LinearCombination`, `AdjointMap`, `TransposeMap` and
@@ -24,7 +34,7 @@
 ## What's new in v3.3
 
 * `AbstractVector`s can now be wrapped by a `LinearMap` just like `AbstractMatrix``
-  typed objects. Upon wrapping, there are not implicitly reshaped to matrices. This
+  typed objects. Upon wrapping, they are not implicitly reshaped to matrices. This
   feature might be helpful, for instance, in the lazy representation of rank-1
   operators `kron(LinearMap(u), v') == ⊗(u, v') == u ⊗ v'` for vectors `u` and `v`.
   The action on vectors,`(u⊗v')*x`, is implemented optimally via `u*(v'x)`.

src/linearcombination.jl

@@ -14,10 +14,18 @@ end
 
 LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps)
 
+# this method avoids the afoldl-mechanism even for LinearMapTuple
 Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::LinearMapTupleOrVector) =
     LinearCombination{promote_type(map(eltype, maps)...)}(maps)
+# this method is required for type stability in the mixed-map-equal-eltype case
 Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::AbstractVector{<:LinearMap{T}}) where {T} =
     LinearCombination{T}(maps)
+# the following two methods are needed to make e.g. `mean` work,
+# for which `f` is some sort of promotion function
+Base.mapreduce(f::F, ::typeof(Base.add_sum), maps::LinearMapTupleOrVector) where {F} =
+    LinearCombination{promote_type(map(eltype, maps)...)}(f.(maps))
+Base.mapreduce(f::F, ::typeof(Base.add_sum), maps::AbstractVector{<:LinearMap{T}}) where {F,T} =
+    LinearCombination{T}(f.(maps))
 
 MulStyle(A::LinearCombination) = MulStyle(A.maps...)
 

test/Project.toml

@@ -6,6 +6,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]

test/linearcombination.jl

@@ -1,4 +1,4 @@
-using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
+using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Statistics
 using LinearMaps: FiveArg, LinearMapTuple, LinearMapVector
 
 @testset "linear combinations" begin
@@ -14,10 +14,15 @@ using LinearMaps: FiveArg, LinearMapTuple, LinearMapVector
     @test (@inferred sum(L.maps::LinearMapTuple)) == L
     Lv = @inferred LinearMaps.LinearCombination{ComplexF64}(fill(CS!, n))
     @test (@inferred sum(Lv.maps::LinearMapVector)) == Lv
+    @test isa((@inferred mean(Lv.maps)),
+        LinearMaps.ScaledMap{ComplexF64,Float64,<:LinearMaps.LinearCombination{ComplexF64,<:LinearMapVector}})
+    @test (@inferred mean(x -> x*x, Lv.maps)) == (@inferred sum(x -> x*x, Lv.maps)/n)
     @test L == Lv
     @test isa((@inferred sum([CS!, LinearMap(randn(eltype(CS!), size(CS!)))])),
         LinearMaps.LinearCombination{<:ComplexF64,<:LinearMapVector})
-    @test isa(sum([CS!, LinearMap(randn(real(eltype(CS!)), size(CS!)))]),
+    A = randn(eltype(CS!), size(CS!))
+    @test (@inferred mean([CS!, LinearMap(A)])) == (@inferred sum([CS!, LinearMap(A)])/2)
+    @test isa(sum([CS!, LinearMap(real(A))]),
         LinearMaps.LinearCombination{<:ComplexF64,<:LinearMapVector})
     for sum1 in (CS!, L, Lv), sum2 in (CS!, L, Lv)
         m1 = sum1 == CS! ? 1 : 10