ContinuousSpace{-1} (#58)

* Support Spline{-1}

* Use ArrowheadMatrix in diff(::ContinuousPloy)

* Update arrowhead.jl

* Update piecewisepolynomial.jl

* tests pass

* start backslash

* expansions work

* differentiation

* Update test_arrowhead.jl

Author: dlfivefifty

src/PiecewiseOrthogonalPolynomials.jl

@@ -7,13 +7,14 @@ import BlockArrays: BlockSlice, block, blockindex, blockvec
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, subblockbandwidths, blockbandwidths, AbstractBandedBlockBandedLayout, layout_replace_in_print_matrix
 import ClassicalOrthogonalPolynomials: grid, plotgrid, ldiv, pad, adaptivetransform_ldiv, Plan, singularities, basis_singularities, singularitiesbroadcast
 import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_transform, plan_grid_transform, grammatrix, weaklaplacian, MulPlan, InvPlan
-import LazyArrays: paddeddata
+import LazyArrays: paddeddata, AbstractLazyLayout
 import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
 import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv, show, summary
 import LinearAlgebra: BlasInt
 import InfiniteArrays: OneToInf
 import FillArrays: AbstractFill
 import MatrixFactorizations: reversecholcopy
+import FillArrays: SquareEye
 
 export PiecewisePolynomial, ContinuousPolynomial, DirichletPolynomial, Derivative, Block, weaklaplacian, grammatrix
 

src/arrowhead.jl

@@ -16,7 +16,8 @@ struct ArrowheadMatrix{T, AA<:AbstractMatrix{T}, BB, CC, DD} <: AbstractBandedBl
     ArrowheadMatrix{T, AA, BB, CC, DD}(A, B, C, D) where {T,AA,BB,CC,DD} = new{T,AA,BB,CC,DD}(A, B, C, D)
 end
 
-ArrowheadMatrix{T}(A, B, C, D) where T = ArrowheadMatrix{T, typeof(A), typeof(B), typeof(C), typeof(D)}(A, B, C, D)
+ArrowheadMatrix{T}(A::AbstractMatrix{T}, B, C, D) where T = ArrowheadMatrix{T, typeof(A), typeof(B), typeof(C), typeof(D)}(A, B, C, D)
+ArrowheadMatrix{T}(A, B, C, D) where T = ArrowheadMatrix(convert(AbstractMatrix{T}, A), B, C, D)
 
 const ArrowheadMatrices = Union{ArrowheadMatrix,Symmetric{<:Any,<:ArrowheadMatrix},Hermitian{<:Any,<:ArrowheadMatrix},
                                 UpperOrLowerTriangular{<:Any,<:ArrowheadMatrix}}
@@ -109,6 +110,7 @@ ArrowheadLayouts = Union{ArrowheadLayout,LazyArrowheadLayout,
                     TriangularLayout{'U', 'U', LazyArrowheadLayout}, TriangularLayout{'L', 'U', LazyArrowheadLayout}}
 arrowheadlayout(_) = ArrowheadLayout()
 arrowheadlayout(::BandedLazyLayouts) = LazyArrowheadLayout()
+arrowheadlayout(::DiagonalLayout{<:AbstractLazyLayout}) = LazyArrowheadLayout()
 symmetriclayout(lay::ArrowheadLayouts) = SymmetricLayout{typeof(lay)}()
 
 MemoryLayout(::Type{<:ArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{D}}}) where D = arrowheadlayout(MemoryLayout(D))

src/continuouspolynomial.jl

@@ -1,3 +1,9 @@
+"""
+   ContinuousPolynomial{0}  is piecewise Legendre
+   ContinuousPolynomial{1}  is hat functions in first block and rest are bubble functions Ultraspherical(-1/2)
+   ContinuousPolynomial{-1} is delta functions (non-boundary) and rest are Ultraspherical(1/2) == diff(Legendre())[:,2:end] in the interior of the elements.
+   Between elements we are consistent with differentiating Legendre: we have delta functions of (-1)^k on the left and -1 on the right.
+"""
 struct ContinuousPolynomial{order,T,P<:AbstractVector} <: AbstractPiecewisePolynomial{order,T,P}
     points::P
 end
@@ -11,8 +17,8 @@ ContinuousPolynomial{o}(P::ContinuousPolynomial) where {o} = ContinuousPolynomia
 PiecewisePolynomial(P::ContinuousPolynomial{o,T}) where {o,T} = PiecewisePolynomial(Legendre{T}(), P.points)
 
 axes(B::ContinuousPolynomial{0}) = axes(PiecewisePolynomial(B))
-axes(B::ContinuousPolynomial{1}) =
-    (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points), Fill(length(B.points) - 1, ∞))))
+axes(B::ContinuousPolynomial{1}) = (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points), Fill(length(B.points) - 1, ∞))))
+axes(B::ContinuousPolynomial{-1}) = (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points)-2, Fill(length(B.points) - 1, ∞))))
 
 show(io::IO, Q::ContinuousPolynomial{λ}) where λ = summary(io, Q)
 summary(io::IO, Q::ContinuousPolynomial{λ}) where λ = print(io, "ContinuousPolynomial{$λ}($(Q.points))")
@@ -45,6 +51,21 @@ function getindex(P::ContinuousPolynomial{1,T}, x::Number, Kk::BlockIndex{1}) wh
     end
 end
 
+function getindex(P::ContinuousPolynomial{-1,T}, x::Number, Kk::BlockIndex{1}) where {T}
+    K, k = block(Kk), blockindex(Kk)
+    if K == Block(1)
+        Spline{-1,T}(P.points)[x, k]
+    else
+        b = searchsortedlast(P.points, x)
+        if b == k
+            α, β = convert(T, P.points[b]), convert(T, P.points[b+1])
+            Ultraspherical{T}(3one(real(T))/2)[affine(α.. β, ChebyshevInterval{real(T)}())[x], Int(K)-1]
+        else
+            zero(T)
+        end
+    end
+end
+
 
 
 factorize(V::SubQuasiArray{T,N,<:ContinuousPolynomial{0},<:Tuple{Inclusion,BlockSlice}}, dims...) where {T,N} =
@@ -214,6 +235,22 @@ function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
 end
 
 
+@simplify function \(D::ContinuousPolynomial{-1}, P::ContinuousPolynomial{0})
+    T = promote_type(eltype(D), eltype(P))
+    @assert D.points == P.points
+    N = length(P.points)
+    R = Ultraspherical{T}(3one(T)/2)\Legendre{T}()
+
+    # The basis for D is defined in each element a..b as diff(legendre(a..b)) with delta functions. But note that
+    # the delta functions don't change! We need to kill off the deltas in conversion.
+    ArrowheadMatrix(_BandedMatrix(Vcat(Fill(-one(T),1,N-1), Fill(one(T),1,N-1)), N-2, 0, 1),
+        (_BandedMatrix(Ones{T}(2,N-1), N-2, 0, 1),),
+        (SquareEye{T}(N-1),),
+        Fill(R[:,2:end], N-1))
+end
+
+
+
 
 ######
 # Gram matrix
@@ -279,7 +316,7 @@ end
 #####
 
 function diff(C::ContinuousPolynomial{1,T}; dims=1) where T
-    # Legendre() \ (D*Weighted(Jacobi(1,1)))
+    # Legendre() \ (D*Ultraspherical(3/2)))
     r = C.points
     N = length(r)
     s = one(T) ./ (r[2:end]-r[1:end-1])
@@ -296,6 +333,18 @@ function diff(C::ContinuousPolynomial{1,T,<:AbstractRange}; dims=1) where T
                                                  Fill(-2s*Eye{T}(∞), N-1))
 end
 
+function diff(P::ContinuousPolynomial{0,T,<:AbstractRange}; dims=1) where T
+    r = P.points
+    N = length(r)
+    s = step(r)
+    
+    ContinuousPolynomial{-1}(r) * ArrowheadMatrix(_BandedMatrix(Vcat(Fill(one(T),1,N-1), Fill(-one(T),1,N-1)), N-2, 0, 1),
+        (),
+        (),
+        Fill(Eye{T}(∞) * (2/s), N-1))
+end
+
+
 function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
     r = C.points
     N = length(r)

test/test_arrowhead.jl

@@ -6,6 +6,14 @@ import Base: oneto, OneTo
 
 
 @testset "ArrowheadMatrix" begin
+    @testset "Constructor" begin
+        n = 4; p = 5;
+        @test ArrowheadMatrix{Float64}(BandedMatrix(0 => 1:n, 1 => 1:n-1, -1 => 1:n-1),
+                                ntuple(_ -> BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)), 2),
+                                ntuple(_ -> BandedMatrix((0 => randn(n), 1 => randn(n-1)), (n-1,n)), 3),
+                            fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2), -1=> randn(p-1)), (p, p)), n-1)) isa ArrowheadMatrix{Float64}
+    end
+
     @testset "Algebra" begin
         n = 4; p = 5;
         A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)),

test/test_continuouspolynomial.jl

@@ -209,4 +209,18 @@ using ContinuumArrays: plan_grid_transform
 
         @test F_xy ≈ f_xy.(x, reshape(y,1,1,size(y)...))
     end
+
+    @testset "Dirac" begin
+        r = range(-1,1; length=4)
+        C = ContinuousPolynomial{1}(r)
+        P = ContinuousPolynomial{0}(r)
+        H = ContinuousPolynomial{-1}(r)
+
+        f = expand(P, exp)
+        @test (H/H\f)[0.1] ≈ exp(0.1)
+        @test diff(f)[0.1] ≈ exp(0.1)
+
+        f = expand(C,exp)
+        @test diff(diff(f))[0.1] ≈ exp(0.1)
+    end
 end