From a2cedf568deb3891b82b34b9171aec7390183291 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Fri, 12 Jun 2026 17:28:52 -0400 Subject: [PATCH] Fix QuadraticSpline edge-case crashes Fixes three construction crashes reported in #542: - n = 2 threw a BoundsError: the knot vector for two data points has boundary multiplicity 2 < degree + 1, so the midpoint B-spline basis evaluation in quadratic_spline_parameters under-indexed. The spline now degenerates to the linear interpolant (alpha = 0, beta = u2 - u1), which evaluation, derivative, and integral all share. - Rational t with t[1] == 0 threw a DivideError: the coefficient eltype was computed as typeof(t[1] / t[1]), which evaluates 0//0. Use typeof(one(eltype(t)) / one(eltype(t))) instead. - Duplicate time points threw a raw SingularException from the tridiagonal collocation solve. Construction now throws an informative ArgumentError. Co-Authored-By: Chris Rackauckas --- src/interpolation_caches.jl | 4 ++-- src/interpolation_utils.jl | 10 +++++++++- src/parameter_caches.jl | 18 +++++++++++++----- test/interpolation_tests.jl | 27 +++++++++++++++++++++++++++ 4 files changed, 51 insertions(+), 8 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 806fb0ea..0d12c3a6 100644 --- a/src/interpolation_caches.jl +++ b/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) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index cce10b3a..5e12edcf 100644 --- a/src/interpolation_utils.jl +++ b/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) diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 891827eb..d30930a4 100644 --- a/src/parameter_caches.jl +++ b/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]) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 11c04650..ce5ad0e6 100644 --- a/test/interpolation_tests.jl +++ b/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