Fix QuadraticSpline edge-case crashes: n=2 BoundsError, Rational t starting at 0, duplicate knots

src/interpolation_caches.jl

@@ -641,7 +641,7 @@ function QuadraticSpline(
     u, t = munge_data(u, t)
 
     n = length(t)
-    dtype_sc = typeof(t[1] / t[1])
+    dtype_sc = typeof(one(eltype(t)) / one(eltype(t)))
     sc = zeros(dtype_sc, n)
     k, A = quadratic_spline_params(t, sc)
     c = A \ u
@@ -671,7 +671,7 @@ function QuadraticSpline(
     u, t = munge_data(u, t)
 
     n = length(t)
-    dtype_sc = typeof(t[1] / t[1])
+    dtype_sc = typeof(one(eltype(t)) / one(eltype(t)))
     sc = zeros(dtype_sc, n)
     k, A = quadratic_spline_params(t, sc)
 

src/interpolation_utils.jl

@@ -60,6 +60,14 @@ function spline_coefficients!(N, d, k, u::AbstractVector)
 end
 
 function quadratic_spline_params(t::AbstractVector, sc::AbstractVector)
+    # Duplicate time points make the collocation system singular
+    if any(i -> t[i] == t[i + 1], 1:(length(t) - 1))
+        throw(
+            ArgumentError(
+                "The time points `t` must be unique for `QuadraticSpline`, but duplicate values were found."
+            )
+        )
+    end
 
     # Create knot vector
     # Don't use x[end-1] as knot to match number of degrees of freedom with data
@@ -72,7 +80,7 @@ function quadratic_spline_params(t::AbstractVector, sc::AbstractVector)
     # - A consists of basis function evaluations in t
     # - c are 1D control points
     n = length(t)
-    dtype_sc = typeof(t[1] / t[1])
+    dtype_sc = typeof(one(eltype(t)) / one(eltype(t)))
 
     diag = Vector{dtype_sc}(undef, n)
     diag_hi = Vector{dtype_sc}(undef, n - 1)

src/parameter_caches.jl

@@ -146,11 +146,19 @@ function QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
 end
 
 function quadratic_spline_parameters(u, t, k, c, sc, idx)
-    tᵢ₊ = (t[idx] + t[idx + 1]) / 2
-    nonzero_coefficient_idxs = spline_coefficients!(sc, 2, k, tᵢ₊)
-    uᵢ₊ = zero(first(u))
-    for j in nonzero_coefficient_idxs
-        uᵢ₊ += sc[j] * c[j]
+    uᵢ₊ = if length(t) == 2
+        # For 2 data points the knot vector has boundary multiplicity 2 < degree + 1,
+        # so the B-spline basis cannot be evaluated at interior points; the spline
+        # degenerates to the linear interpolant (α = 0, β = u₂ - u₁).
+        (u[1] + u[2]) / 2
+    else
+        tᵢ₊ = (t[idx] + t[idx + 1]) / 2
+        nonzero_coefficient_idxs = spline_coefficients!(sc, 2, k, tᵢ₊)
+        uᵢ₊ = zero(first(u))
+        for j in nonzero_coefficient_idxs
+            uᵢ₊ += sc[j] * c[j]
+        end
+        uᵢ₊
     end
     α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊
     β = 4 * (uᵢ₊ - u[idx]) - (u[idx + 1] - u[idx])

test/interpolation_tests.jl

@@ -1029,6 +1029,33 @@ end
     A = @inferred(QuadraticSpline(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-2.0)
     @test_throws DataInterpolations.RightExtrapolationError A(2.0)
+
+    # Two data points degenerate to the linear interpolant
+    u = [1.0, 3.0]
+    t = [0.0, 1.0]
+    A = @inferred(QuadraticSpline(u, t))
+    @test A(0.0) == 1.0
+    @test A(1.0) == 3.0
+    @test A(0.25) == 1.5
+    @test DataInterpolations.derivative(A, 0.5) == 2.0
+    @test DataInterpolations.integral(A, 0.0, 1.0) == 2.0
+    A = QuadraticSpline(u, t; cache_parameters = true)
+    @test A(0.25) == 1.5
+    A = QuadraticSpline([[1.0, 2.0], [3.0, 6.0]], t)
+    @test A(0.5) == [2.0, 4.0]
+
+    # Rational data with t[1] == 0; u = (t + 1)^2 is reproduced exactly
+    u = [1 // 1, 4 // 1, 9 // 1]
+    t = [0 // 1, 1 // 1, 2 // 1]
+    A = QuadraticSpline(u, t)
+    @test A(1 // 2) == 9 // 4
+    @test A(3 // 2) == 25 // 4
+
+    # Duplicate time points throw an informative error
+    @test_throws ArgumentError QuadraticSpline(
+        [1.0, 2.0, 3.0, 4.0, 5.0], [0.0, 1.0, 1.0, 2.0, 3.0]
+    )
+    @test_throws ArgumentError QuadraticSpline([1.0, 2.0, 3.0], [0.0, 0.0, 1.0])
 end
 
 @testset "CubicSpline Interpolation" begin