Use FastInterpolations for splining

Project.toml

@@ -8,6 +8,7 @@ Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 ForwardDiffChainRules = "c9556dd2-1aed-4cfe-8560-1557cf593001"
 GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
@@ -49,6 +50,7 @@ Bessels = "0.2.8"
 ChainRulesCore = "1.25.1"
 CommonSolve = "0.2.4"
 DataInterpolations = "8.0.0"
+FastInterpolations = "0.4.8"
 ForwardDiff = "0.10, 1.0.0"
 ForwardDiffChainRules = "0.3.0"
 GraphRecipes = "0.5.13"

src/utils.jl

@@ -1,4 +1,4 @@
-import DataInterpolations: CubicSpline, CubicHermiteSpline
+using FastInterpolations
 import Symbolics: taylor, operation, sorted_arguments, unwrap
 import Base: identity, replace
 using QuadGK
@@ -105,27 +105,26 @@ function spline(y, x)
     dx[end] == 0 && return spline(y[begin:end-1], x[begin:end-1]) # endpoints are duplicated when ODE solver ends with callback; in that case remove it
     all(diff(x) .< 0) && return spline(reverse(y), reverse(x)) # reverse if monotonically decreasing
     all(diff(x) .> 0) || error("x is not monotonically increasing")
-    return CubicSpline(y, x; extrapolation = ExtrapolationType.Linear)
+    return cubic_interp(x, y)
 end
 
 function spline(y, ẏ, x)
-    dx = diff(x)
     all(diff(x) .< 0) && return spline(reverse(y), reverse(ẏ), reverse(x)) # reverse if monotonically decreasing
     all(diff(x) .> 0) || error("x is not monotonically increasing")
-    return CubicHermiteSpline(ẏ, y, x) # TODO: use PCHIP instead? https://docs.sciml.ai/DataInterpolations/stable/methods/#PCHIP-Interpolation
+    return hermite_interp(x, y, ẏ) # TODO: use PCHIP instead? https://docs.sciml.ai/DataInterpolations/stable/methods/#PCHIP-Interpolation
 end
 
 function spline(sol::ODESolution)
     ts = sol.t
     Nu, _ = size(sol)
     us = map(u -> SVector{Nu}(u), sol(ts, Val{0}))
     dus = map(u -> SVector{Nu}(u), sol(ts, Val{1}))
-    return CubicHermiteSpline(dus, us, ts; extrapolation = ExtrapolationType.Extension, cache_parameters = true) # TODO: use PCHIP instead? https://docs.sciml.ai/DataInterpolations/stable/methods/#PCHIP-Interpolation
+    return hermite_interp(ts, us, dus) # TODO: use PCHIP instead? https://docs.sciml.ai/DataInterpolations/stable/methods/#PCHIP-Interpolation
 end
 
 function dummyspline(N)
     nans = SVector{N}(zeros(N))
-    return CubicHermiteSpline([nans, nans], [nans, nans], [0.0, 1.0])
+    return hermite_interp([0.0, 1.0], [nans, nans], [nans, nans])
 end
 
 # https://github.com/SciML/ModelingToolkit.jl/issues/4202#issuecomment-3799583124
@@ -183,7 +182,7 @@ function mtkcompile_spline(sys::System, vars)
     # Build mapping from variables to spline parameters
     splname = :bgspline
     spldummy = dummyspline(length(vars))
-    uprototype = spldummy.u[begin]
+    uprototype = spldummy.y[begin]
     spl, = @parameters $splname::Any
 
     iv = ModelingToolkit.get_iv(sys)

test/runtests.jl

@@ -208,6 +208,15 @@ end
     @test τspans[3] == (sol.bg.t[begin], sol.bg.t[end]) # very high k; should start at same time as background
 end
 
+@testset "Background splining" begin
+    bgsol = solve(prob).bg
+    spl = SymBoltz.spline(bgsol)
+    t = range(bgsol.t[begin], bgsol.t[end], length=500)
+    u1 = stack(bgsol(t))
+    u2 = stack(spl.(t))
+    @test isapprox(u1, u2)
+end
+
 @testset "Automatic background/thermodynamics splining" begin
     sol = solve(prob, 1.0) # solve with one perturbation mode to activate splining
     τs = SymBoltz.timeseries.(sol, log10(M.g.a), range(-8, 0, length=100)) # TODO a => as syntax