diff --git a/Project.toml b/Project.toml index 0ebecfc7..74c645b4 100644 --- a/Project.toml +++ b/Project.toml @@ -49,7 +49,7 @@ QuadGK = "2.9.1" RecipesBase = "1.3.4" Reexport = "1" SafeTestsets = "0.1" -SciMLTesting = "1" +SciMLTesting = "2.1" SparseConnectivityTracer = "1" StableRNGs = "1" StaticArrays = "1.9.8" diff --git a/docs/src/arclength_interpolation.md b/docs/src/arclength_interpolation.md index ffc88936..18a5121c 100644 --- a/docs/src/arclength_interpolation.md +++ b/docs/src/arclength_interpolation.md @@ -53,7 +53,7 @@ SmoothArcLengthInterpolation(::AbstractMatrix, ::AbstractMatrix) ## Method derivation -Say we have an ordered set of points $u_1, \ldots, u_n \in \mathbb{R}^N$ and we want to make a lightweight $C^1$ smooth interpolation by arc-length $\tilde{\gamma}: [0,T] \rightarrow \mathbb{R}^N$ through these points. The first part is easy, just pick your favorite established interpolation method that achieves $C^1$ smoothness. The arc-length part however turns out to be quite [nasty](https://ijpam.eu/contents/2006-31-3/10/10.pdf). Here I propose a method that is quite general and cheap to compute. +Say we have an ordered set of points $u_1, \ldots, u_n \in \mathbb{R}^N$ and we want to make a lightweight $C^1$ smooth interpolation by arc-length $\tilde{\gamma}: [0,T] \rightarrow \mathbb{R}^N$ through these points. The first part is easy, just pick your favorite established interpolation method that achieves $C^1$ smoothness. The arc-length part however turns out to be quite nasty; see Gil, "On the Arc Length Parametrization Problem", International Journal of Pure and Applied Mathematics 31(3), 2006. Here I propose a method that is quite general and cheap to compute. ### The 2-dimensional case diff --git a/docs/src/extrapolation_methods.md b/docs/src/extrapolation_methods.md index d7475b2b..603a9b57 100644 --- a/docs/src/extrapolation_methods.md +++ b/docs/src/extrapolation_methods.md @@ -15,6 +15,10 @@ plot(A) Extrapolation behavior can be set left and right of the data simultaneously with the `extrapolation` keyword, or left and right separately with the `extrapolation_left` and `extrapolation_right` keywords respectively. +```@docs +ExtrapolationType +``` + ## `ExtrapolationType.None` This extrapolation type will throw an error when the input `t` is beyond the data in the specified direction. diff --git a/docs/src/methods.md b/docs/src/methods.md index 65bbac8d..68489420 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -228,3 +228,7 @@ And the parameters show the issue: ```@example tutorial A.pmin ``` + +```@docs +Curvefit +``` diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 73f0c34b..b3db98b8 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -9,10 +9,38 @@ using RecipesBase: RecipesBase, @recipe, @series using PrettyTables: PrettyTables, pretty_table using ForwardDiff: ForwardDiff using EnumX: EnumX, @enumx -using StaticArrays: MVector, SVector +using StaticArrays: SVector import FindFirstFunctions import FindFirstFunctions: Guesser +""" + ExtrapolationType + +Enumeration of extrapolation modes used by interpolation constructors through +the `extrapolation`, `extrapolation_left`, and `extrapolation_right` keyword +arguments. + +## Values + + - `ExtrapolationType.None`: throw an error outside the interpolation interval. + - `ExtrapolationType.Constant`: use the nearest endpoint value. + - `ExtrapolationType.Linear`: extend the interpolation linearly from the nearest interval. + - `ExtrapolationType.Extension`: evaluate the interpolation formula outside the data range. + - `ExtrapolationType.Periodic`: wrap the query point periodically onto the data range. + - `ExtrapolationType.Reflective`: reflect the query point back onto the data range. + +## Examples + +```julia +using DataInterpolations + +t = [0.0, 1.0] +u = [1.0, 3.0] +A = LinearInterpolation(u, t; extrapolation = ExtrapolationType.Extension) + +A(1.5) +``` +""" @enumx ExtrapolationType None Constant Linear Extension Periodic Reflective include("parameter_caches.jl") @@ -168,7 +196,13 @@ struct CurvefitCache{ end end -# Define an empty function, so that it can be extended via `DataInterpolationsOptimExt` +""" + Curvefit(u, t, m, p0, alg; extrapolate = false, ub = nothing, lb = nothing) + +Construct an interpolation by fitting the model `m(t, p)` to data values `u` at +timepoints `t`, using `p0` as the initial parameter guess and `alg` as the +optimization algorithm. +""" function Curvefit() error("CurveFit requires loading Optim and ForwardDiff, e.g. `using Optim, ForwardDiff`") end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 72bf7c6f..0f1c7135 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -70,41 +70,54 @@ end # A degree-`d` B-spline has only `d + 1` nonzero basis functions at any point, so # evaluation needs a scratch buffer of just `d + 1` entries — not one sized to the # full knot/control-point vector. `bspline_nonzero_coefficients` computes those -# values into a stack-allocated `MVector` (no heap allocation, hence reentrant / -# thread-safe, #532) for degrees up to `BSPLINE_STACK_MAXLEN - 1`; callers fall back -# to the heap `spline_coefficients!` for larger degrees, which essentially never occur. +# values into an immutable static vector for degrees up to `BSPLINE_STACK_MAXLEN - 1`; +# callers fall back to the heap `spline_coefficients!` for larger degrees, which +# essentially never occur. const BSPLINE_STACK_MAXLEN = 16 +@inline function _static_setindex(N::SVector{L, T}, x, idx) where {L, T} + return SVector{L, T}(ntuple(j -> ifelse(j == idx, convert(T, x), N[j]), Val(L))) +end + # Returns `(vals, offset, m)`: `vals[l]` is the basis weight of control point # `offset + l` for `l in 1:m`. `ncp` is the number of control points (needed only to # place the single nonzero weight at the right boundary). Mirrors the locator and # recurrence of `spline_coefficients!`, but writes into a local window indexed `1:m` # instead of absolute knot indices. @inline function bspline_nonzero_coefficients(d::Integer, k, u::T, ncp::Integer) where {T} - N = zero(MVector{BSPLINE_STACK_MAXLEN, T}) + N = zero(SVector{BSPLINE_STACK_MAXLEN, T}) if u == k[1] - @inbounds N[1] = one(u) - return SVector(N), 0, 1 + N = _static_setindex(N, one(u), 1) + return N, 0, 1 elseif u == k[end] - @inbounds N[1] = one(u) - return SVector(N), ncp - 1, 1 + N = _static_setindex(N, one(u), 1) + return N, ncp - 1, 1 end i = searchsortedlast(k, u) # Local index of global knot index `g` is `g + off` (so `i - d → 1`, `i → d + 1`). off = d + 1 - i @inbounds begin - N[i + off] = one(u) + N = _static_setindex(N, one(u), i + off) for deg in 1:d - N[i - deg + off] = (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) * - N[i - deg + 1 + off] + N = _static_setindex( + N, + (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) * + N[i - deg + 1 + off], + i - deg + off, + ) for j in (i - deg + 1):(i - 1) - N[j + off] = (u - k[j]) / (k[j + deg] - k[j]) * N[j + off] + - (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * N[j + 1 + off] + N = _static_setindex( + N, + (u - k[j]) / (k[j + deg] - k[j]) * N[j + off] + + (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * + N[j + 1 + off], + j + off, + ) end - N[i + off] = (u - k[i]) / (k[i + deg] - k[i]) * N[i + off] + N = _static_setindex(N, (u - k[i]) / (k[i + deg] - k[i]) * N[i + off], i + off) end end - return SVector(N), i - d - 1, d + 1 + return N, i - d - 1, d + 1 end function quadratic_spline_params(t::AbstractVector, sc::AbstractVector) diff --git a/test/Extensions/Project.toml b/test/Extensions/Project.toml index 19742179..ca6ff891 100644 --- a/test/Extensions/Project.toml +++ b/test/Extensions/Project.toml @@ -9,16 +9,13 @@ SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" -[sources] -DataInterpolations = {path = "../.."} - [compat] DataInterpolations = "9" ForwardDiff = "1" LinearAlgebra = "1.10" Mooncake = "0.4.175, 0.5" SafeTestsets = "0.1" -SciMLTesting = "1.7" +SciMLTesting = "2.1" SparseConnectivityTracer = "1" Test = "1.10" Zygote = "0.6.77, 0.7" diff --git a/test/qa/Project.toml b/test/qa/Project.toml index 51784413..a8882212 100644 --- a/test/qa/Project.toml +++ b/test/qa/Project.toml @@ -7,13 +7,10 @@ SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[sources] -DataInterpolations = {path = "../.."} - [compat] AllocCheck = "0.2" Aqua = "0.8" SafeTestsets = "0.0.1, 0.1" -SciMLTesting = "1.7" +SciMLTesting = "2.1" StaticArrays = "1.9.8" Test = "1" diff --git a/test/qa/alloc_tests.jl b/test/qa/alloc_tests.jl index a5e187b4..017e0faa 100644 --- a/test/qa/alloc_tests.jl +++ b/test/qa/alloc_tests.jl @@ -1,10 +1,9 @@ using DataInterpolations -using DataInterpolations: derivative using AllocCheck: @check_allocs using StaticArrays: SVector @check_allocs(test_allocs(itp, x) = itp(x)) # Reuse function definition to save on compilation time -@check_allocs(test_derivative_allocs(itp, x) = derivative(itp, x)) +@check_allocs(test_derivative_allocs(itp, x) = DataInterpolations.derivative(itp, x)) @testset "Allocation-free interpolation tests" begin @testset "LinearInterpolation" begin diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 10ac9ce4..8d975b93 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -14,4 +14,5 @@ run_qa( # declares them public. all_qualified_accesses_are_public = (; ignore = (:Dual, :derivative, :value)), ), + api_docs_kwargs = (; rendered = true), )