Add Derivative for spline

Author: dlfivefifty

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.7
 IntervalSets
 LazyArrays
+BandedMatrices

src/ContinuumArrays.jl

@@ -1,14 +1,16 @@
 module ContinuumArrays
-using IntervalSets, LinearAlgebra, LazyArrays
+using IntervalSets, LinearAlgebra, LazyArrays, BandedMatrices
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear
 import Base.Broadcast: materialize
+import LazyArrays: Mul2
+import BandedMatrices: AbstractBandedLayout
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
 import .QuasiArrays: _length, checkindex, Adjoint, Transpose
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative
 
 ####
 # Interval indexing support

src/QuasiArrays/QuasiArrays.jl

@@ -2,7 +2,7 @@ module QuasiArrays
 using Base, LinearAlgebra, LazyArrays
 import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                 Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, zero,
-                map, eachindex, eltype
+                map, eachindex, eltype, first, last
 import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
                 DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type, _default_type,
                 _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
@@ -14,7 +14,7 @@ import Base.Broadcast: materialize
 
 import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose
 
-import LazyArrays: MemoryLayout, UnknownLayout
+import LazyArrays: MemoryLayout, UnknownLayout, Mul2
 
 export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize
 

src/QuasiArrays/matmul.jl

@@ -1,28 +1,44 @@
 
-@inline MemoryLayout(A::AbstractQuasiArray{T}) where T = UnknownLayout()
-
-const QuasiArrayMulArray{TV, styleA, styleB, p, q, T, V} =
-    Mul{TV, styleA, styleB, <:AbstractQuasiArray{T,p}, <:AbstractArray{V,q}}
+const QuasiArrayMulArray{styleA, styleB, p, q, T, V} =
+    Mul2{styleA, styleB, <:AbstractQuasiArray{T,p}, <:AbstractArray{V,q}}
 
+const QuasiArrayMulQuasiArray{styleA, styleB, p, q, T, V} =
+    Mul2{styleA, styleB, <:AbstractQuasiArray{T,p}, <:AbstractQuasiArray{V,q}}
 ####
 # Matrix * Vector
 ####
-let (p,q) = (2,1)
-    global const QuasiMatMulVec{TV, styleA, styleB, T, V} = QuasiArrayMulArray{TV, styleA, styleB, p, q, T, V}
-end
+const QuasiMatMulVec{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 1, T, V}
 
-axes(M::QuasiMatMulVec) = (axes(M.A,1),)
 
-function getindex(M::QuasiMatMulVec{T}, k::Real) where T
-    ret = zero(T)
-    for j = 1:size(M.A,2)
-        ret += M.A[k,j] * M.B[j]
+function getindex(M::QuasiMatMulVec, k::Real)
+    A,B = M.factors
+    ret = zero(eltype(M))
+    @inbounds for j = 1:size(A,2)
+        ret += A[k,j] * B[j]
     end
     ret
 end
 
+QuasiMatMulMat{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 2, T, V}
+QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styleB, 2, 2, T, V}
+
+
 *(A::AbstractQuasiArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 *(A::AbstractQuasiArray, B::AbstractArray) = materialize(Mul(A,B))
 *(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 inv(A::AbstractQuasiArray) = materialize(Inv(A))
-*(A::Inv{<:Any,<:Any,<:AbstractQuasiArray}, B::AbstractQuasiArray) = materialize(Mul(A,B))
+*(A::Inv{<:Any,<:AbstractQuasiArray}, B::AbstractQuasiArray) = materialize(Mul(A,B))
+
+
+_Mul(A::Mul, B::Mul) = Mul(A.factors..., B.factors...)
+_Mul(A::Mul, B) = Mul(A.factors..., B)
+_Mul(A, B::Mul) = Mul(A, B.factors...)
+_Mul(A, B) = Mul(A, B)
+_lsimplify2(A, B...) = _Mul(A, _lsimplify(B...))
+_lsimplify2(A::Mul, B...) = _lsimplify2(A.factors..., B...)
+_lsimplify(A) = materialize(A)
+_lsimplify(A, B) = materialize(Mul(A,B))
+_lsimplify(A, B, C, D...) = _lsimplify2(materialize(Mul(A,B)), C, D...)
+lsimplify(M::Mul) = _lsimplify(M.factors...)
+
+*(A::AbstractQuasiArray, B::Mul) = lsimplify(_Mul(A, B))

src/bases/splines.jl

@@ -35,28 +35,68 @@ function getindex(B::HeavisideSpline{T}, x::Real, k::Int) where T
 end
 
 
-function convert(::Type{SymTridiagonal}, AB::Mul{T,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) where T
-    Ac,B = AB.A, AB.B
+function copyto!(dest::SymTridiagonal, AB::Mul2{<:Any,<:Any,<:Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) where T
+    Ac,B = AB.factors
     A = parent(Ac)
-    @assert A.points == B.points
+    A.points == B.points || throw(ArgumentError())
+    dv,ev = dest.dv,dest.ev
     x = A.points; n = length(x)
-    dv = Vector{T}(undef, n)
+    length(dv) == n || throw(DimensionMismatch())
+
     dv[1] = (x[2]-x[1])/3
-    for k = 2:n-1
+    @inbounds for k = 2:n-1
         dv[k] = (x[k+1]-x[k-1])/3
     end
     dv[n] = (x[n] - x[n-1])/3
 
-    SymTridiagonal(dv, diff(x)./6)
+    @inbounds for k = 1:n-1
+        ev[k] = (x[k+1]-x[k])/6
+    end
+
+    dest
+end
+
+## Mass matrix
+function similar(AB::Mul2{<:Any,<:Any,<:Adjoint{<:Any,<:LinearSpline},<:LinearSpline}, ::Type{T}) where T
+    n = size(AB,1)
+    SymTridiagonal(Vector{T}(undef, n), Vector{T}(undef, n-1))
 end
 #
-materialize(M::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
-    convert(SymTridiagonal, M)
+materialize(M::Mul2{<:Any,<:Any,<:Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
+    copyto!(similar(M, eltype(M)), M)
 
-function materialize(M::Mul{T,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:HeavisideSpline},<:HeavisideSpline}) where T
-    Ac, B = M.A, M.B
+function materialize(M::Mul2{<:Any,<:Any,<:Adjoint{<:Any,<:HeavisideSpline},<:HeavisideSpline})
+    Ac, B = M.factors
     axes(Ac,2) == axes(B,1) || throw(DimensionMismatch("axes must be same"))
     A = parent(Ac)
     A.points == B.points || throw(ArgumentError("Cannot multiply incompatible splines"))
     Diagonal(diff(A.points))
 end
+
+## Derivative
+function copyto!(dest::Mul2{<:Any,<:Any,<:HeavisideSpline},
+                 M::Mul2{<:Any,<:Any,<:Derivative,<:LinearSpline})
+    D, L = M.factors
+    H, A = dest.factors
+    x = H.points
+
+    axes(dest) == axes(M) || throw(DimensionMismatch("axes must be same"))
+    x == L.points || throw(ArgumentError("Cannot multiply incompatible splines"))
+    bandwidths(A) == (0,1) || throw(ArgumentError("Not implemented"))
+
+    d = diff(x)
+    A[band(0)] .= (-).(d)
+    A[band(1)] .= d
+
+    dest
+end
+
+function similar(M::Mul2{<:Any,<:Any,<:Derivative,<:LinearSpline}, ::Type{T}) where T
+    D, B = M.factors
+    n = size(B,2)
+    Mul(HeavisideSpline{T}(B.points),
+        BandedMatrix{T}(undef, (n-1,n), (0,1)))
+end
+
+materialize(M::Mul2{<:Any,<:Any,<:Derivative,<:LinearSpline}) =
+    copyto!(similar(M, eltype(M)), M)

src/operators.jl

@@ -17,14 +17,23 @@ function getindex(δ::DiracDelta{T}, x::Real) where T
 end
 
 
-function materialize(M::Mul{<:Any,<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelta},<:AbstractQuasiVector})
+function materialize(M::Mul2{<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelta},<:AbstractQuasiVector})
     A, B = M.A, M.B
     axes(A,2) == axes(A,1) || throw(DimensionMismatch())
     B[parent(A).x]
 end
 
-function materialize(M::Mul{<:Any,<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelta},<:AbstractQuasiMatrix})
-    A, B = M.A, M.B
+function materialize(M::Mul2{<:Any,<:Any,<:QuasiArrays.Adjoint{<:Any,<:DiracDelta},<:AbstractQuasiMatrix})
+    A, B = M.factors
     axes(A,2) == axes(B,1) || throw(DimensionMismatch())
     B[parent(A).x,:]
 end
+
+struct Derivative{T,A} <: AbstractQuasiVector{T}
+    axis::A
+end
+
+Derivative{T}(axis::A) where {T,A} = Derivative{T,A}(axis)
+Derivative(axis) = Derivative{Float64}(axis)
+
+axes(D::Derivative) = (D.axis, D.axis)

test/runtests.jl

@@ -10,8 +10,6 @@ using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, Test
     @test Base.IndexStyle(δ) === Base.IndexLinear()
 end
 
-
-
 @testset "HeavisideSpline" begin
     H = HeavisideSpline([1,2,3])
     @test axes(H) === (axes(H,1),axes(H,2)) === (1.0..3.0, Base.OneTo(2))
@@ -69,9 +67,15 @@ end
     L = LinearSpline([1,2,3])
     @test δ'L ≈ [0.8, 0.2, 0.0]
 
+    @test L'L == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])
 end
 
+@testset "Derivative" begin
+    L = LinearSpline([1,2,3])
+    f = L*[1,2,4]
+    D = Derivative(axes(L,1))
+    fp = D*f
 
-
-
-@test L'L == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])
+    @test fp[1.1] ≈ 1
+    @test fp[2.2] ≈ 2
+end