Support algebra operations on expansions

Author: dlfivefifty

src/ContinuumArrays.jl

@@ -7,7 +7,7 @@ import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert,
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport,
                         adjointlayout, LdivApplyStyle, arguments, _arguments, call, broadcastlayout, layout_getindex,
-                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles
+                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 import FillArrays: AbstractFill, getindex_value, SquareEye

src/bases/bases.jl

@@ -149,6 +149,53 @@ function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:Abstrac
     T \ L.B[p]
 end
 
+
+##
+# Algebra
+##
+
+# struct ExpansionLayout <: MemoryLayout end
+# applylayout(::Type{typeof(*)}, ::BasisLayout, _) = ExpansionLayout()
+
+const Expansion{T,Space<:Basis,Coeffs<:AbstractVector} = ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}
+
+for op in (:*, :\)
+    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Number, f::Expansion)
+        S,c = arguments(f)
+        S * broadcast($op, x, c)
+    end
+end
+for op in (:*, :/)
+    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), f::Expansion, x::Number)
+        S,c = arguments(f)
+        S * broadcast($op, c, x)
+    end
+end
+
+
+function broadcastbasis(::typeof(+), a, b)
+    a ≠ b && error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
+    a
+end
+
+broadcastbasis(::typeof(-), a, b) = broadcastbasis(+, a, b)
+
+for op in (:+, :-)
+    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), f::Expansion, g::Expansion)
+        S,c = arguments(f)
+        T,d = arguments(g)
+        ST = broadcastbasis($op, S, T)
+        ST * $op((ST \ S) * c , (ST \ T) * d)
+    end
+end
+
+@eval function ==(f::Expansion, g::Expansion)
+    S,c = arguments(f)
+    T,d = arguments(g)
+    ST = broadcastbasis(+, S, T)
+    (ST \ S) * c == (ST \ T) * d
+end
+
 ## materialize views
 
 # materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =

src/operators.jl

@@ -158,6 +158,7 @@ function diff(d::AbstractQuasiVector)
     Derivative(x)*d
 end
 
+^(D::Derivative, k::Integer) = ApplyQuasiArray(^, D, k)
 
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
 #     f::F

test/runtests.jl

@@ -1,6 +1,6 @@
 using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, FastTransforms, Test
 import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, 
-                        MappedBasisLayout, igetindex, TransformFactorization, Weight, WeightedBasisLayout
+                        MappedBasisLayout, igetindex, TransformFactorization, Weight, WeightedBasisLayout, Expansion
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivApplyStyle
 
@@ -124,10 +124,28 @@ end
     end
 end
 
-@testset "Derivative" begin
+@testset "Algebra" begin
     L = LinearSpline([1,2,3])
     f = L*[1,2,4]
+    g = L*[5,6,7]
+    
+    @test f isa Expansion
+    @test 2f isa Expansion
+    @test f*2 isa Expansion
+    @test 2\f isa Expansion
+    @test f/2 isa Expansion
+    @test f+g isa Expansion
+    @test f-g isa Expansion
     @test f[1.2] == 1.2
+    @test (2f)[1.2] == (f*2)[1.2] == 2.4
+    @test (2\f)[1.2] == (f/2)[1.2] == 0.6
+    @test (f+g)[1.2] ≈ f[1.2] + g[1.2]
+    @test (f-g)[1.2] ≈ f[1.2] - g[1.2]
+end
+
+@testset "Derivative" begin
+    L = LinearSpline([1,2,3])
+    f = L*[1,2,4]
 
     D = Derivative(axes(L,1))
     @test ApplyStyle(*,typeof(D),typeof(L)) isa SimplifyStyle
@@ -151,11 +169,13 @@ end
     @test fp[1.1] ≈ 1
     @test fp[2.2] ≈ 2
 
-
     fp = D*f
     @test length(fp.args) == 2
     @test fp[1.1] ≈ 1
     @test fp[2.2] ≈ 2
+
+    @test D^2 isa ApplyQuasiMatrix{Float64,typeof(*)}
+    
 end
 
 @testset "Weak Laplacian" begin