Expand functions in P'(x .^2 .* P)

.github/workflows/ci.yml

@@ -20,7 +20,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - 'lts'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -39,7 +39,7 @@ Infinities = "0.1"
 IntervalSets = "0.7"
 LazyArrays = "2"
 Makie = "0.20, 0.21"
-QuasiArrays = "0.11.5"
+QuasiArrays = "0.11.8"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
 julia = "1.10"

src/bases/bases.jl

@@ -145,15 +145,54 @@
     a * (A \ b)
 end
 
-_broadcast_mul_ldiv(::Tuple{ScalarLayout,AbstractBasisLayout}, A, B) =
-    _broadcast_mul_ldiv((ScalarLayout(),UnknownLayout()), A, B)
+_broadcast_mul_ldiv(::Tuple{ScalarLayout,AbstractBasisLayout}, A, B) = _broadcast_mul_ldiv((ScalarLayout(),UnknownLayout()), A, B)
 _broadcast_mul_ldiv(_, A, B) = copy(Ldiv{typeof(MemoryLayout(A)),UnknownLayout}(A,B))
 
 copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
 copy(L::Ldiv{<:MappedBasisLayouts,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
 
 
 
+# multiplication operators, reexpand in basis A
+@inline function _broadcast_mul_adj(::Tuple{Any,AbstractBasisLayout}, Ac, B)
+    a,b = arguments(B)
+    @assert a isa AbstractQuasiVector # Only works for vec .* mat
+    A = Ac'
+    ab = (A * (A \ a)) .* b # broadcasted should be overloaded
+    MemoryLayout(ab) isa BroadcastLayout && return Ac*transform_ldiv(A, ab)
+    Ac*ab
+end
+
+@inline function _broadcast_mul_adj(::Tuple{Any,ApplyLayout{typeof(*)}}, Ac, B)
+    a,b = arguments(B)
+    @assert a isa AbstractQuasiVector # Only works for vec .* mat
+    args = arguments(*, b)
+    *(Ac*(a .* first(args)), tail(args)...)
+end
+
+
+function _broadcast_mul_adj(::Tuple{ScalarLayout,Any}, Ac, B)
+    a,b = arguments(B)
+    a * (Ac*b)
+end
+
+function _broadcast_mul_adj(::Tuple{ScalarLayout,ApplyLayout{typeof(*)}}, Ac, B)
+    a,b = arguments(B)
+    a * (Ac*b)
+end
+
+_broadcast_mul_adj(::Tuple{ScalarLayout,AbstractBasisLayout}, A, B) = _broadcast_mul_adj((ScalarLayout(),UnknownLayout()), A, B)
+_broadcast_mul_adj(_, A, B) = copy(Mul{typeof(MemoryLayout(A)),UnknownLayout}(A,B))
+
+_broadcast_mul_adj_simplifiable(_, ::AbstractBasisLayout) = Val(true)
+_broadcast_mul_adj_simplifiable(_, ::ApplyLayout{typeof(*)}) = Val(true)
+_broadcast_mul_adj_simplifiable(::ScalarLayout, _) = Val(true)
+_broadcast_mul_adj_simplifiable(::ScalarLayout, ::ApplyLayout{typeof(*)}) = Val(true)
+_broadcast_mul_adj_simplifiable(::ScalarLayout, ::AbstractBasisLayout) = Val(true)
+_broadcast_mul_adj_simplifiable(_, _) = Val(false)
+
+simplifiable(L::Mul{<:AdjointBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul_adj_simplifiable(map(MemoryLayout,arguments(L.B))...)
+copy(L::Mul{<:AdjointBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul_adj(map(MemoryLayout,arguments(L.B)), L.A, L.B)
 
 
 """
@@ -651,6 +690,7 @@
 ####
 
 simplifiable(::Mul{<:AdjointBasisLayout, <:AbstractBasisLayout}) = Val(true)
+@inline simplifiable(L::Mul{<:AdjointBasisLayout,ApplyLayout{typeof(*)}}) = simplifiable(*, L.A, first(arguments(*, L.B)))
 function copy(M::Mul{<:AdjointBasisLayout, <:AbstractBasisLayout})
     A = (M.A)'
     A == M.B && return grammatrix(A)

test/test_chebyshev.jl

@@ -1,4 +1,5 @@
-using ContinuumArrays, LinearAlgebra, FastTransforms, QuasiArrays, ArrayLayouts, Base64, LazyArrays, Test
+using ContinuumArrays, LinearAlgebra, QuasiArrays, ArrayLayouts, Base64, LazyArrays, Test
+using FastTransforms
 import ContinuumArrays: Basis, Weight, Map, LazyQuasiArrayStyle, TransformFactorization,
                         ExpansionLayout, checkpoints, MappedBasisLayout, MappedWeightedBasisLayout,
                         SubWeightedBasisLayout, WeightedBasisLayout, WeightLayout, basis, grammatrix
@@ -159,6 +160,25 @@ Base.:(==)(::FooBasis, ::FooBasis) = true
 
         ã = T * (T \ a)
         @test T \ (ã .* ã) ≈ [1.5,1,0.5,0,0]
+
+        @test T'*(a .* T) isa Matrix
+        @test T'*(a .* (T * (T \ a))) isa Vector
+        @test_broken T'f isa Vector
+        @test T'ã isa Vector
+        @test T'*(ã .* ã) isa Vector
+        @test (2T)'*(a .* T) isa Matrix
+        @test T'*(2T) isa Matrix
+        @test T'*(2T*randn(5)) isa Vector
+        @test (2T)'*(T*(1:5)) ≈ T'*(2T*(1:5)) ≈ T'BroadcastQuasiMatrix(*, 2, T*(1:5))
+        @test T' * (a .* (T * (1:5))) ≈ T' * ((a .* T) * (1:5))
+        @test T'BroadcastQuasiMatrix(*, 2, 2T) == 4*(T'T)
+
+        @test LazyArrays.simplifiable(*, T', T*(1:5)) == Val(true)
+        @test LazyArrays.simplifiable(*, T', (a .* (T * (1:5)))) == Val(true)
+        @test LazyArrays.simplifiable(*, T', a .* T) == Val(true)
+        @test LazyArrays.simplifiable(*, T', 2T) == Val(true)
+        @test LazyArrays.simplifiable(*, T', BroadcastQuasiMatrix(*, 2, T*(1:5))) == Val(true)
+        @test LazyArrays.simplifiable(*, T', BroadcastQuasiMatrix(*, 2, 2T)) == Val(true)
     end
 
     @testset "sum/dot/diff" begin

test/test_splines.jl

@@ -362,7 +362,7 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
         @test Δ == -(*(B',D',D,B))
         @test Δ == -(B'D'D*B)
         @test Δ == -((B'D')*(D*B))
-        @test_broken Δ == -B'*(D'D)*B
+        @test Δ == -B'*(D'D)*B
         @test Δ == -(B'*(D'D)*B)
 
         f = L*exp.(L.points) # project exp(x)