Support algebra operations on expansions (#54)

* Support algebra operations on expansions

* Start new simplify

* Tests pass!

* v0.3

* Remove \ simplify

* * for AlephInfinity

* Update Project.toml

* Update runtests.jl

* Use new simplifiable from LazyArrays.jl

* Update ContinuumArrays.jl

* Update Project.toml

* Update bases.jl

* Update Project.toml

* Update bases.jl

* require Julia v1.5

* Update bases.jl

* Increase coverage

Author: dlfivefifty

.travis.yml

@@ -5,8 +5,7 @@ os:
   - osx
   - windows
 julia:
-  - 1.3
-  - 1.4
+  - 1.5
   - nightly
 matrix:
   allow_failures:

Project.toml

@@ -1,24 +1,28 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.2.5"
+version = "0.3.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
-ArrayLayouts = "0.2.4, 0.3"
-BandedMatrices = "0.15"
-FillArrays = "0.8.2, 0.9"
+ArrayLayouts = "0.4.3"
+BandedMatrices = "0.15.17"
+FillArrays = "0.9.3"
+InfiniteArrays = "0.8"
 IntervalSets = "0.3.2, 0.4, 0.5"
-LazyArrays = "0.16"
-QuasiArrays = "0.2.2"
-julia = "1.3"
+LazyArrays = "0.17.1"
+QuasiArrays = "0.3"
+julia = "1.5"
 
 [extras]
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"

src/ContinuumArrays.jl

@@ -1,41 +1,75 @@
 module ContinuumArrays
-using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays
+using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, InfiniteArrays
 import Base: @_inline_meta, @_propagate_inbounds_meta, axes, getindex, convert, prod, *, /, \, +, -, ==,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, diff,
                 first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
-                getproperty
+                getproperty, isone, iszero, zero, abs, <, ≤, >, ≥, string
 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
+import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
+                        adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex,
+                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
+                        simplifiable, _simplify
 import LinearAlgebra: pinv
 import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 import FillArrays: AbstractFill, getindex_value, SquareEye
-
+import ArrayLayouts: mul
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
-                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, quasimulapplystyle,
-                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
-                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, quasildivapplystyle, _factorize
+                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
+                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
+                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize
+import InfiniteArrays: Infinity
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, fullmaterialize, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, WeightedBasis, grid, transform, affine
 
 ####
 # Interval indexing support
 ####
 struct AlephInfinity{N} <: Integer end
 
+isone(::AlephInfinity) = false
+iszero(::AlephInfinity) = false
+
 ==(::AlephInfinity, ::Int) = false
 ==(::Int, ::AlephInfinity) = false
 
 *(::AlephInfinity{N}, ::AlephInfinity{N}) where N = AlephInfinity{N}()
+*(::AlephInfinity{N}, ::Infinity) where N = AlephInfinity{N}()
+*(::Infinity, ::AlephInfinity{N}) where N = AlephInfinity{N}()
+function *(a::Integer, b::AlephInfinity)
+    a > 0 || throw(ArgumentError("$a is negative"))
+    b
+end
+
+*(a::AlephInfinity, b::Integer) = b*a
+
+
+abs(a::AlephInfinity) = a
+zero(::AlephInfinity) = 0
+
+for OP in (:<, :≤)
+    @eval begin
+        $OP(::Real, ::AlephInfinity) = true
+        $OP(::AlephInfinity, ::Real) = false
+    end
+end
+
+for OP in (:>, :≥)
+    @eval begin
+        $OP(::Real, ::AlephInfinity) = false
+        $OP(::AlephInfinity, ::Real) = true
+    end
+end
+
 
 const ℵ₁ = AlephInfinity{1}()
 
+string(::AlephInfinity{1}) = "ℵ₁"
+
 show(io::IO, F::AlephInfinity{1}) where N =
     print(io, "ℵ₁")
 
 
-const QMul2{A,B} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B}}
+const QMul2{A,B} = Mul{<:Any, <:Any, <:A,<:B}
 const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}
 
 cardinality(::AbstractInterval) = ℵ₁

src/bases/bases.jl

@@ -4,14 +4,14 @@ abstract type Weight{T} <: LazyQuasiVector{T} end
 
 const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}
 
-struct WeightLayout <: MemoryLayout end
-abstract type AbstractBasisLayout <: MemoryLayout end
+struct WeightLayout <: AbstractQuasiLazyLayout end
+abstract type AbstractBasisLayout <: AbstractQuasiLazyLayout end
 struct BasisLayout <: AbstractBasisLayout end
 struct SubBasisLayout <: AbstractBasisLayout end
 struct MappedBasisLayout <: AbstractBasisLayout end
 struct WeightedBasisLayout <: AbstractBasisLayout end
 
-abstract type AbstractAdjointBasisLayout <: MemoryLayout end
+abstract type AbstractAdjointBasisLayout <: AbstractQuasiLazyLayout end
 struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
 struct AdjointSubBasisLayout <: AbstractAdjointBasisLayout end
 struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end
@@ -25,9 +25,7 @@ adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
 broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = WeightedBasisLayout()
 broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = WeightedBasisLayout()
 
-combine_mul_styles(::AbstractBasisLayout) = LazyQuasiArrayApplyStyle()
-combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()
-
+# Default is lazy
 ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
 pinv(J::Basis) = apply(pinv,J)
 
@@ -37,10 +35,6 @@ function ==(A::Basis, B::Basis)
     false
 end
 
-@inline quasildivapplystyle(::AbstractBasisLayout, ::AbstractBasisLayout) = LdivApplyStyle()
-@inline quasildivapplystyle(::AbstractBasisLayout, _) = LdivApplyStyle()
-@inline quasildivapplystyle(_, ::AbstractBasisLayout) = LdivApplyStyle()
-
 
 @inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
 @inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) = 
@@ -138,7 +132,7 @@ copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiV
     transform_ldiv(L.A, L.B)
 
 function copy(L::Ldiv{ApplyLayout{typeof(*)},<:AbstractBasisLayout})
-    args = arguments(L.A)
+    args = arguments(ApplyLayout{typeof(*)}(), L.A)
     @assert length(args) == 2 # temporary
     apply(\, last(args), apply(\, first(args), L.B))
 end
@@ -149,6 +143,55 @@ 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}}
+
+basis(v::AbstractQuasiVector) = v.args[1]
+
+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}}}) =
@@ -164,9 +207,8 @@ end
 _sub_getindex(A, kr, jr) = A[kr, jr]
 _sub_getindex(A, ::Slice, ::Slice) = A
 
-function copy(M::QMul2{<:QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
-                        <:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
-    Ac, B = M.args
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
+             B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     A = Ac'
     PA,PB = parent(A),parent(B)
     kr,jr = parentindices(B)
@@ -175,14 +217,12 @@ end
 
 
 # Differentiation of sub-arrays
-function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
-    A, B = M.args
+@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})
     P = parent(B)
     (Derivative(axes(P,1))*P)[parentindices(B)...]
 end
 
-function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
-    A, B = M.args
+@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(B)
     kr,jr = parentindices(B)
     (Derivative(axes(P,1))*P*kr.A)[kr,jr]
@@ -220,15 +260,18 @@ end
 # SubLayout behaves like ApplyLayout{typeof(*)}
 
 combine_mul_styles(::SubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
-_arguments(::SubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
+_mul_arguments(::SubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
+arguments(::SubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
 call(::SubBasisLayout, ::SubQuasiArray) = *
 
 combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
-_arguments(::AdjointSubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
+_mul_arguments(::AdjointSubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
 arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
 call(::AdjointSubBasisLayout, ::SubQuasiArray) = *
 
-function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
+copy(M::Mul{AdjointSubBasisLayout,SubBasisLayout}) = apply(*, arguments(M.A)..., arguments(M.B)...)
+
+function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
     A = parent(V)
     _,jr = parentindices(V)
     first(jr) ≥ 1 || throw(BoundsError())
@@ -249,8 +292,8 @@ function __sum(::SubBasisLayout, Vm, dims)
     @assert dims == 1
     sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
 end
-function __sum(::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
-    a = arguments(V)
+function __sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
+    a = arguments(LAY, V)
     first(apply(*, sum(a[1]; dims=1), tail(a)...))
 end
 

src/bases/splines.jl

@@ -48,22 +48,19 @@ end
 
 
 ## Mass matrix
-
-ApplyStyle(::typeof(*), ::Type{<:QuasiAdjoint{<:Any,<:LinearSpline}}, ::Type{<:LinearSpline}) = 
-    SimplifyStyle()    
-
-
 function similar(AB::QMul2{<:QuasiAdjoint{<:Any,<:LinearSpline},<:LinearSpline}, ::Type{T}) where T
     n = size(AB,1)
     SymTridiagonal(Vector{T}(undef, n), Vector{T}(undef, n-1))
 end
 #
-copy(M::QMul2{<:QuasiAdjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:LinearSpline}, B::LinearSpline) 
+    M = Mul(Ac, B)
     copyto!(similar(M, eltype(M)), M)
+end
 
 function copyto!(dest::SymTridiagonal,
                  AB::QMul2{<:QuasiAdjoint{<:Any,<:LinearSpline},<:LinearSpline}) where T
-    Ac,B = AB.args
+    Ac,B = AB.A,AB.B
     A = parent(Ac)
     A.points == B.points || throw(ArgumentError())
     dv,ev = dest.dv,dest.ev
@@ -92,13 +89,9 @@ end
 
 
 ## Derivative
-ApplyStyle(::typeof(*), ::Type{<:Derivative}, ::Type{<:LinearSpline}) = SimplifyStyle()    
-
-
-
 function copyto!(dest::MulQuasiMatrix{<:Any,<:Tuple{<:HeavisideSpline,<:Any}},
                  M::QMul2{<:Derivative,<:LinearSpline})
-    D, L = M.args
+    D, L = M.A, M.B
     H, A = dest.args
     x = H.points
 
@@ -114,20 +107,15 @@ function copyto!(dest::MulQuasiMatrix{<:Any,<:Tuple{<:HeavisideSpline,<:Any}},
 end
 
 function similar(M::QMul2{<:Derivative,<:LinearSpline}, ::Type{T}) where T
-    D, B = M.args
+    D, B = M.A, M.B
     n = size(B,2)
     ApplyQuasiMatrix(*, HeavisideSpline{T}(B.points),
         BandedMatrix{T}(undef, (n-1,n), (0,1)))
 end
 
-copy(M::QMul2{<:Derivative,<:LinearSpline}) =
+@simplify function *(D::Derivative, L::LinearSpline)
+    M = Mul(D, L)
     copyto!(similar(M, eltype(M)), M)
-
-ApplyStyle(::typeof(*), ::Type{<:QuasiAdjoint{<:Any,<:LinearSpline}}, ::Type{<:QuasiAdjoint{<:Any,<:Derivative}}) = SimplifyStyle()        
-
-function copy(M::QMul2{<:QuasiAdjoint{<:Any,<:LinearSpline},<:QuasiAdjoint{<:Any,<:Derivative}})
-    Bc,Ac = M.args
-    apply(*,Ac',Bc')'
 end
 
 
@@ -138,4 +126,15 @@ end
 function _sum(A::HeavisideSpline, dims)
     @assert dims == 1
     permutedims(diff(A.points))
+end
+
+function _sum(P::LinearSpline, dims)
+    d = diff(P.points)
+    ret = Array{float(eltype(d))}(undef, length(d)+1)
+    ret[1] = d[1]/2
+    for k = 2:length(d)
+        ret[k] = (d[k-1] + d[k])/2
+    end
+    ret[end] = d[end]/2
+    permutedims(ret)
 end