Add Mul, mass matrices for splines

Author: dlfivefifty

REQUIRE

@@ -1,2 +1,3 @@
 julia 0.7
 IntervalSets
+LazyArrays

examples/BSplines.jl

@@ -1,5 +1,5 @@
-using ContinuumArrays, LazyArrays, IntervalSets
-import ContinuumArrays: AbstractAxisMatrix, ℵ₀
+using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, Test
+import ContinuumArrays: AbstractAxisMatrix, ℵ₀, materialize
 import Base: axes, getindex, convert
 
 struct Spline{order,T} <: AbstractAxisMatrix{T}
@@ -11,7 +11,7 @@ const HeavisideSpline{T} = Spline{0,T}
 
 Spline{o}(pts::AbstractVector{T}) where {o,T} = Spline{o,float(T)}(pts)
 
-axes(B::Spline{o}) where o = (first(B.points)..last(B.points), Base.OneTo(length(B.points)-o-1))
+axes(B::Spline{o}) where o = (first(B.points)..last(B.points), Base.OneTo(length(B.points)+o-1))
 
 function getindex(B::LinearSpline{T}, x::Real, k::Int) where T
     x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())
@@ -37,34 +37,83 @@ function getindex(B::HeavisideSpline{T}, x::Real, k::Int) where T
     return one(T)
 end
 
-# getindex(B::LinearSpline, ::Colon, k::Int) = Mul(B, [Zeros{Int}(k-1); 1; Zeros{Int}(size(B,2)-k)])
 
-# function convert(::Type{SymTridiagonal}, AB::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline})
-#     Ac,B = AB.A, AB.B
-#     A = parent(Ac)
-#     @assert A.points == B.points
-#     x = A.points
-#     SymTridiagonal(x, x/2) # TODO fix
-# end
+function convert(::Type{SymTridiagonal}, AB::Mul{T,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) where T
+    Ac,B = AB.A, AB.B
+    A = parent(Ac)
+    @assert A.points == B.points
+    x = A.points; n = length(x)
+    dv = Vector{T}(undef, n)
+    dv[1] = (x[2]-x[1])/3
+    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)
+end
 #
-# materialize(M::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
-#     convert(SymTridiagonal, M)
+materialize(M::Mul{<:Any,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:LinearSpline},<:LinearSpline}) =
+    convert(SymTridiagonal, M)
+
+function materialize(M::Mul{T,<:Any,<:Any,<:ContinuumArrays.Adjoint{<:Any,<:HeavisideSpline},<:HeavisideSpline}) where T
+    Ac, B = M.A, M.B
+    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
+
 
 ## tests
 
-B = HeavisideSpline([1,2,3])
-@test size(B) == (ℵ₀, 2)
+H = HeavisideSpline([1,2,3])
+@test size(H) == (ℵ₀, 2)
+
+@test_throws BoundsError H[0.1, 1]
+@test H[1.1,1] === H'[1,1.1] === transpose(H)[1,1.1] === 1.0
+@test H[2.1,1] === H'[1,2.1] === transpose(H)[1,2.1] === 0.0
+@test H[1.1,2] === 0.0
+@test H[2.1,2] === 1.0
+@test_throws BoundsError H[2.1,3]
+@test_throws BoundsError H'[3,2.1]
+@test_throws BoundsError transpose(H)[3,2.1]
+@test_throws BoundsError H[3.1,2]
+
+@test all(H[[1.1,2.1], 1] .=== H'[1,[1.1,2.1]] .=== transpose(H)[1,[1.1,2.1]] .=== [1.0,0.0])
+@test all(H[1.1,1:2] .=== [1.0,0.0])
+@test all(H[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])
+
+@test_throws BoundsError H[[0.1,2.1], 1]
+
+
+L = LinearSpline([1,2,3])
+@test size(L) == (ℵ₀, 3)
+
+@test_throws BoundsError L[0.1, 1]
+@test L[1.1,1] == L'[1,1.1] == transpose(L)[1,1.1] ≈ 0.9
+@test L[2.1,1] === L'[1,2.1] === transpose(L)[1,2.1] === 0.0
+@test L[1.1,2] ≈ 0.1
+@test L[2.1,2] ≈ 0.9
+@test L[2.1,3] == L'[3,2.1] == transpose(L)[3,2.1] ≈ 0.1
+@test_throws BoundsError L[3.1,2]
+L[[1.1,2.1], 1]
+@test L[[1.1,2.1], 1] == L'[1,[1.1,2.1]] == transpose(L)[1,[1.1,2.1]] ≈ [0.9,0.0]
+@test L[1.1,1:2] ≈ [0.9,0.1]
+@test L[[1.1,2.1], 1:2] ≈ [0.9 0.1; 0.0 0.9]
+
+@test_throws BoundsError L[[0.1,2.1], 1]
+
 
-@test_throws BoundsError B[0.1, 1]
-@test B[1.1,1] === 1.0
-@test B[2.1,1] === 0.0
-@test B[1.1,2] === 0.0
-@test B[2.1,2] === 1.0
-@test_throws BoundsError B[2.1,3]
-@test_throws BoundsError B[3.1,2]
+f = H*[1,2]
+@test axes(f) == (1.0..3.0,)
+@test f[1.1] ≈ 1
+@test f[2.1] ≈ 2
 
-@test all(B[[1.1,2.1], 1] .=== [1.0,0.0])
-@test all(B[1.1,1:2] .=== [1.0,0.0])
-@test all(B[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])
+f = L*[1,2,4]
+@test axes(f) == (1.0..3.0,)
+@test f[1.1] ≈ 1.1
+@test f[2.1] ≈ 2.2
 
-@test_throws BoundsError B[[0.1,2.1], 1]
+@test H'H == Eye(2)
+@test L'L == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])

src/ContinuumArrays.jl

@@ -1,30 +1,27 @@
 module ContinuumArrays
-using Base, LinearAlgebra, IntervalSets
-import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
-                Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, zero,
-                map, eachindex
-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,
-                index_shape, to_shape, unsafe_length, @nloops, @ncall
+include("axisarrays/AbstractAxisArrays.jl")
+using .AbstractAxisArrays
 
 
-import LinearAlgebra: transpose, adjoint
+####
+# Interval indexing support
+####
 
-abstract type AbstractAxisArray{T,N} end
-AbstractAxisVector{T} = AbstractAxisArray{T,1}
-AbstractAxisMatrix{T} = AbstractAxisArray{T,2}
-AbstractAxisVecOrMat{T} = Union{AbstractAxisVector{T}, AbstractAxisMatrix{T}}
+using IntervalSets
+import .AbstractAxisArrays: _length, checkindex, Adjoint, Transpose
+import Base: @_inline_meta
 
 struct ℵ₀ <: Number end
 _length(::AbstractInterval) = ℵ₀
-_length(d) = length(d)
 
-size(A::AbstractAxisArray) = _length.(axes(A))
-axes(A::AbstractAxisArray) = error("Override axes for $(typeof(A))")
+checkindex(::Type{Bool}, inds::AbstractInterval, i::Real) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
+function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
+    @_inline_meta
+    b = true
+    for i in I
+        b &= checkindex(Bool, inds, i)
+    end
+    b
+end
 
-include("indices.jl")
-include("abstractaxisarray.jl")
-include("adjtrans.jl")
-include("multidimensional.jl")
 end