Fix linear spline interpolation of closed parametric curves (#96)

* Fix linear spline interpolation of closed curves

* Add tests

Author: jipolanco

src/Collocation/Collocation.jl

@@ -12,13 +12,16 @@ using ..DifferentialOps
 using ..Recombinations: num_constraints
 
 using BandedMatrices
+using LinearAlgebra: LinearAlgebra
 using SparseArrays
 using StaticArrays: Size, SVector, MVector, SMatrix, MMatrix
 
 # For Documenter only:
 using ..Recombinations: RecombinedBSplineBasis
 using ..BoundaryConditions: Natural
 
+const IdentityMatrix = typeof(LinearAlgebra.I)
+
 include("matrix.jl")
 include("points.jl")
 include("cyclic.jl")
@@ -140,6 +143,19 @@ _default_matrix_type(::BSplineOrder, ::Type{<:PeriodicBSplineBasis}, ::Type{T})
 _default_matrix_type(::BSplineOrder{4}, ::Type{<:PeriodicBSplineBasis}, ::Type{T}) where {T} =
     CyclicTridiagonalMatrix{T}
 
+# In this case the collocation matrix is simply the identity matrix.
+# We assume knots are located at collocation points (otherwise this is not true).
+_default_matrix_type(::BSplineOrder{2}, ::Type{<:PeriodicBSplineBasis}, ::Type{T}) where {T} =
+    IdentityMatrix
+
+allocate_collocation_matrix(::Type{IdentityMatrix}, x, B) = LinearAlgebra.I
+
+function collocation_matrix!(C::IdentityMatrix, B::PeriodicBSplineBasis, x, ::Derivative{0}; kwargs...)
+    ts = @view parent(knots(B))[1:end-1]  # knot locations in [0, L)
+    ts == x || throw(ArgumentError("expected collocation points to match B-spline knots"))
+    C
+end
+
 allocate_collocation_matrix(::Type{M}, x, B) where {M <: AbstractMatrix} =
     M(undef, length(x), length(B))
 

src/SplineInterpolations/SplineInterpolations.jl

@@ -2,7 +2,7 @@ module SplineInterpolations
 
 using ..BSplines
 using ..Collocation
-using ..Collocation: CyclicTridiagonalMatrix
+using ..Collocation: CyclicTridiagonalMatrix, IdentityMatrix
 using ..Splines
 
 using BandedMatrices
@@ -50,7 +50,7 @@ locations.
 """
 struct SplineInterpolation{
         S <: Spline,
-        F <: Factorization,
+        F <: Union{Factorization, IdentityMatrix},
         Points <: AbstractVector,
         ColMatCopy <: Union{Nothing, AbstractMatrix},
     } <: SplineWrapper{S}
@@ -66,12 +66,12 @@ struct SplineInterpolation{
     C_copy :: ColMatCopy
 
     function SplineInterpolation(
-            B::AbstractBSplineBasis, C::Factorization, x::AbstractVector,
+            B::AbstractBSplineBasis, C::Union{Factorization, IdentityMatrix}, x::AbstractVector,
             ::Type{T};
             C_copy = nothing,
         ) where {T}
         N = length(B)
-        size(C) == (N, N) ||
+        C isa IdentityMatrix || size(C) == (N, N) ||
             throw(DimensionMismatch("collocation matrix has wrong dimensions"))
         length(x) == N ||
             throw(DimensionMismatch("wrong number of collocation points"))
@@ -100,10 +100,12 @@ end
 
 # Factorise in-place when possible.
 _factorise!(C::AbstractMatrix) = lu!(C)
+_factorise!(I::IdentityMatrix) = I  # identity matrix (for periodic linear splines)
 _factorise!(C::SparseMatrixCSC) = lu(C)  # SparseMatrixCSC doesn't support `lu!`
 
 # Create copy if required by the matrix type.
 _maybe_copy_matrix(::AbstractMatrix) = nothing
+_maybe_copy_matrix(::IdentityMatrix) = nothing   # identity matrix (for periodic linear splines)
 _maybe_copy_matrix(C::CyclicTridiagonalMatrix) = copy(C)  # periodic cubic splines
 
 interpolation_points(S::SplineInterpolation) = S.x

test/periodic.jl

@@ -4,23 +4,21 @@ using LinearAlgebra
 using Random
 using Test
 
-# For now this only works for cubic splines.
-function test_parametric_curve(ord::BSplineOrder{4})
+# For now this only works for linear (k = 2) and cubic (k = 4) splines.
+function test_parametric_curve(ord::BSplineOrder{k}) where {k}
     L = 2π
     θs = range(0, L; length = 6)[1:5]  # input angles in [0, 2π)
     points = [SVector(cos(θ), sin(θ)) for θ ∈ θs]
-    S = interpolate(θs, copy(points), ord, Periodic(L))
+    S = @inferred interpolate(θs, copy(points), ord, Periodic(L))
     a, b = boundaries(basis(S))
     @test a == 0
     @test b == L
     # Verify continuity of the curve
-    S′ = Derivative(1) * S
-    S″ = Derivative(2) * S
-    S‴ = Derivative(3) * S
     @test S(a) == S(b)
-    @test S′(a) == S′(b)
-    @test S″(a) == S″(b)
-    @test S‴(a) == S‴(b)
+    for n ∈ 1:(k - 1)
+        dS = Derivative(n) * S
+        @test dS(a) == dS(b)
+    end
     nothing
 end
 
@@ -117,7 +115,8 @@ function test_periodic_splines(ord::BSplineOrder)
         end
 
         @testset "Approximations" begin
-            @test isapprox(app_in(x), ftest(x); rtol=0.01)
+            rtol = k ≤ 2 ? 0.02 : 0.01
+            @test isapprox(app_in(x), ftest(x); rtol)  # requires slightly larger tolerance for k = 2
             @test isapprox(app_L2(x), ftest(x); rtol=0.01)
         end
     end
@@ -148,7 +147,9 @@ function test_periodic_splines(ord::BSplineOrder)
             coefficients(Sp) .= @view coefs[(δ+1):(δ+length(Sp))]
             @test Sn(x) == Sp(x)
             @test (Derivative(1) * Sn)(x) == (Derivative(1) * Sp)(x)
-            @test (Derivative(2) * Sn)(x) == (Derivative(2) * Sp)(x)
+            if k > 2
+                @test (Derivative(2) * Sn)(x) == (Derivative(2) * Sp)(x)
+            end
             ∫n = integral(Sn)
             ∫p = integral(Sp)
             @test ∫n(x′) - ∫n(x) ≈ ∫p(x′) - ∫p(x)
@@ -159,8 +160,10 @@ function test_periodic_splines(ord::BSplineOrder)
         xs = @inferred collocation_points(B)
         @test iseven(k) == (xs == breaks[begin:end-1])  # this is the default for even k
         C = @inferred collocation_matrix(B, xs)
-        @test cond(Array(C)) < 1e3
-        @test all(≈(1), sum(C; dims=2))  # partition of unity
+        if C !== LinearAlgebra.I  # case of linear splines (k = 2)
+            @test cond(Array(C)) < 1e3
+            @test all(≈(1), sum(C; dims=2))  # partition of unity
+        end
 
         # Interpolate manually
         ys = ftest.(xs)
@@ -192,7 +195,7 @@ function test_periodic_splines(ord::BSplineOrder)
         end
     end
 
-    if k == 4
+    if k ∈ (2, 4)
         @testset "Parametric closed curve" test_parametric_curve(ord)
     end