Merge pull request #616 from maltezfaria/cubic-free-bc

Fix free boundary condition on `Cubic`

Author: mkitti

src/b-splines/cubic.jl

@@ -232,16 +232,33 @@ end
 continuous derivative at the second-to-last cell boundary; this means
 `y_1'''(2) = y_2'''(2)`, yielding
 
-    1 cm -3 c + 3 cp -1 cpp = 0
-"""
-function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{<:Free}) where {T,TC}
-    dl, d, du = inner_system_diags(T,n,degree)
 
-    specs = WoodburyMatrices.sparse_factors(T, n,
-                                  (1, n, du[1]),
-                                  (n, 1, dl[end])
-                                  )
+    -1 cm + 4 c - 6 cp + 4 cpp - 1 cppp = 0
 
-    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
+This enforces continuity of the third derivative, which is appropriate for cubic splines
+(in contrast to the quadratic case which enforces continuity of the second derivative).
+
+"""
+function prefiltering_system(
+    ::Type{T},
+    ::Type{TC},
+    n::Int,
+    degree::Cubic{<:Free},
+) where {T,TC}
+    dl, d, du = inner_system_diags(T, n, degree)
+    d[1] = d[end] = -1
+    du[1] = dl[end] = 4
+
+    specs = WoodburyMatrices.sparse_factors(
+        T,
+        n,
+        (1, 3, -6),
+        (1, 4, 4),
+        (1, 5, -1),
+        (n, n - 2, -6),
+        (n, n - 3, 4),
+        (n, n - 4, -1),
+    )
+
+    return Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
 end

test/b-splines/cubic.jl

@@ -90,4 +90,21 @@
     itp = interpolate(rand(5, 9), BSpline(Cubic(Flat(OnGrid()))))
     @test itp == Interpolations.BSplineInterpolation(itp.coefs, itp.it, itp.parentaxes)
     @test_throws ArgumentError Interpolations.BSplineInterpolation(itp.coefs, itp.it, (2:3, 2:10))
+
+    # Issue 594
+    @testset "Free BC" begin
+        # should be exact for polynomials of degree <= 3 (up to rounding error)
+        ix = 1:15
+        f = x -> x^3 - 2x^2 + 3x - 4
+        A = map(f, ix)
+        itp = interpolate(A, BSpline(Cubic(Free(OnGrid()))))
+        @test all(x -> itp(x) ≈ f(x), 1:0.1:15)
+
+        # also in 2d
+        f2(x, y) = x^3 - 2x^2 + 3x - 4 + 2y^3 - 2y^2 + 3y - 4
+        iy = 1:15
+        A = map(I -> f2(I[1], I[2]), Iterators.product(ix, iy))
+        itp = interpolate(A, BSpline(Cubic(Free(OnGrid()))))
+        @test all(x -> itp(x...) ≈ f2(x...), Iterators.product(1:0.1:15, 1:0.1:15))
+    end
 end