From 2cca6312629fb9241df1dfbe0762b637a52bfe6d Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 11:06:34 +0100 Subject: [PATCH] Merge remote-tracking branch 'upstream/master' into stable_quadratic_spline --- .github/workflows/CompatHelper.yml | 2 +- .github/workflows/Downgrade.yml | 4 +- .github/workflows/Invalidations.yml | 2 + .github/workflows/TagBot.yml | 2 +- .github/workflows/Tests.yml | 5 + Project.toml | 4 +- README.md | 20 ++ docs/Project.toml | 2 +- docs/src/index.md | 19 ++ ext/DataInterpolationsChainRulesCoreExt.jl | 45 ++- joss/paper.md | 28 +- src/DataInterpolations.jl | 17 +- src/derivatives.jl | 63 +++- src/integral_inverses.jl | 16 +- src/interpolation_caches.jl | 335 +++++++++++++++++---- src/interpolation_methods.jl | 65 +++- src/interpolation_utils.jl | 48 ++- src/parameter_caches.jl | 29 +- test/derivative_tests.jl | 44 ++- test/interpolation_tests.jl | 90 +++++- test/zygote_tests.jl | 33 +- 21 files changed, 724 insertions(+), 149 deletions(-) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 35cc34ba..f95ed99c 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -13,7 +13,7 @@ jobs: CompatHelper: runs-on: ubuntu-latest steps: - - uses: julia-actions/setup-julia@780022b48dfc0c2c6b94cfee6a9284850107d037 + - uses: julia-actions/setup-julia@9b79636afcfb07ab02c256cede01fe2db6ba808c with: version: 1.3 - name: Pkg.add("CompatHelper") diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index 4546ebd0..3030e057 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -21,14 +21,14 @@ jobs: group: - Core version: - - '1' + - '1.10' os: - ubuntu-latest - macos-latest - windows-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2.3.0 + - uses: julia-actions/setup-julia@v2.6.0 with: version: ${{ matrix.version }} - uses: julia-actions/julia-downgrade-compat@v1 diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 34eb7a92..0133ad34 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -13,3 +13,5 @@ jobs: evaluate-invalidations: name: "Evaluate Invalidations" uses: "SciML/.github/.github/workflows/invalidations.yml@v1" + with: + julia-version: "1.10" diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 43c8233a..3d99dc0b 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -25,7 +25,7 @@ jobs: if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' runs-on: ubuntu-latest steps: - - uses: JuliaRegistries/TagBot@aa5545ecce2ae3b2cd7d3a8a0a286ec6bf25838f + - uses: JuliaRegistries/TagBot@29c35fccdd29270e3560ede0c1b77b4b6e12abce with: token: ${{ secrets.GITHUB_TOKEN }} # Edit the following line to reflect the actual name of the GitHub Secret containing your private key diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 4d53b02c..9235c8c0 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -24,11 +24,16 @@ jobs: strategy: fail-fast: false matrix: + version: + - "1" + - "lts" + - "pre" os: - "ubuntu-latest" - "macos-latest" - "windows-latest" uses: "SciML/.github/.github/workflows/tests.yml@v1" with: + julia-version: "${{ matrix.version }}" os: "${{ matrix.os }}" secrets: "inherit" diff --git a/Project.toml b/Project.toml index 59d27ef9..f91cd3bc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DataInterpolations" uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" -version = "6.1.0" +version = "6.5.2" [deps] FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" @@ -39,7 +39,7 @@ Reexport = "1" RegularizationTools = "0.6" SafeTestsets = "0.1" StableRNGs = "1" -Symbolics = "5.29" +Symbolics = "5.29, 6" Test = "1" Zygote = "0.6.70" julia = "1.10" diff --git a/README.md b/README.md index 5348685f..df100c0d 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,7 @@ [![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle) +[![DOI](https://joss.theoj.org/papers/10.21105/joss.06917/status.svg)](https://doi.org/10.21105/joss.06917) DataInterpolations.jl is a library for performing interpolations of one-dimensional data. By "data interpolations" we mean techniques for interpolating possibly noisy data, and thus @@ -93,3 +94,22 @@ The series types defined are: - `:quintic_hermite_spline` By and large, these accept the same keywords as their function counterparts. + +## Citing + +If you use this software in your work, please cite: + +```bib +@article{Bhagavan2024, + doi = {10.21105/joss.06917}, + url = {https://doi.org/10.21105/joss.06917}, + year = {2024}, + publisher = {The Open Journal}, + volume = {9}, + number = {101}, + pages = {6917}, + author = {Sathvik Bhagavan and Bart de Koning and Shubham Maddhashiya and Christopher Rackauckas}, + title = {DataInterpolations.jl: Fast Interpolations of 1D data}, + journal = {Journal of Open Source Software} +} +``` diff --git a/docs/Project.toml b/docs/Project.toml index 540ffc47..6440a7b7 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -20,4 +20,4 @@ OrdinaryDiffEq = "6" Plots = "1" RegularizationTools = "0.6" StableRNGs = "1" -Symbolics = "5.29" +Symbolics = "5.29, 6.0" diff --git a/docs/src/index.md b/docs/src/index.md index f2075f8e..1b19d50f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -65,6 +65,25 @@ The series types defined are: By and large, these accept the same keywords as their function counterparts. +## Citing + +If you use this software in your work, please cite: + +```bib +@article{Bhagavan2024, + doi = {10.21105/joss.06917}, + url = {https://doi.org/10.21105/joss.06917}, + year = {2024}, + publisher = {The Open Journal}, + volume = {9}, + number = {101}, + pages = {6917}, + author = {Sathvik Bhagavan and Bart de Koning and Shubham Maddhashiya and Christopher Rackauckas}, + title = {DataInterpolations.jl: Fast Interpolations of 1D data}, + journal = {Journal of Open Source Software} +} +``` + ## Contributing - Please refer to the diff --git a/ext/DataInterpolationsChainRulesCoreExt.jl b/ext/DataInterpolationsChainRulesCoreExt.jl index 988e519b..c59a7f11 100644 --- a/ext/DataInterpolationsChainRulesCoreExt.jl +++ b/ext/DataInterpolationsChainRulesCoreExt.jl @@ -4,17 +4,27 @@ if isdefined(Base, :get_extension) LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, BSplineInterpolation, BSplineApprox, get_idx, get_parameters, - _quad_interp_indices + _quad_interp_indices, munge_data using ChainRulesCore else using ..DataInterpolations: _interpolate, derivative, AbstractInterpolation, LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation, BSplineInterpolation, BSplineApprox, get_parameters, - _quad_interp_indices + _quad_interp_indices, munge_data using ..ChainRulesCore end +function ChainRulesCore.rrule(::typeof(munge_data), u, t) + u_out, t_out = munge_data(u, t) + + # For now modifications by munge_data not supported + @assert (u == u_out && t == t_out) + + munge_data_pullback = Δ -> (NoTangent(), Δ[1], Δ[2]) + (u_out, t_out), munge_data_pullback +end + function ChainRulesCore.rrule( ::Type{LinearInterpolation}, u, t, I, p, extrapolate, cache_parameters) A = LinearInterpolation(u, t, I, p, extrapolate, cache_parameters) @@ -51,16 +61,21 @@ function ChainRulesCore.rrule( end function u_tangent(A::LinearInterpolation, t, Δ) - out = zero(A.u) + out = zero.(A.u) idx = get_idx(A, t, A.iguesser) t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) - out[idx] = Δ * (one(eltype(out)) - t_factor) - out[idx + 1] = Δ * t_factor + if eltype(out) <: Number + out[idx] = Δ * (one(eltype(out)) - t_factor) + out[idx + 1] = Δ * t_factor + else + @. out[idx] = Δ * (true - t_factor) + @. out[idx + 1] = Δ * t_factor + end out end function u_tangent(A::QuadraticInterpolation, t, Δ) - out = zero(A.u) + out = zero.(A.u) i₀, i₁, i₂ = _quad_interp_indices(A, t, A.iguesser) t₀ = A.t[i₀] t₁ = A.t[i₁] @@ -68,9 +83,15 @@ function u_tangent(A::QuadraticInterpolation, t, Δ) Δt₀ = t₁ - t₀ Δt₁ = t₂ - t₁ Δt₂ = t₂ - t₀ - out[i₀] = Δ * (t - A.t[i₁]) * (t - A.t[i₂]) / (Δt₀ * Δt₂) - out[i₁] = -Δ * (t - A.t[i₀]) * (t - A.t[i₂]) / (Δt₀ * Δt₁) - out[i₂] = Δ * (t - A.t[i₀]) * (t - A.t[i₁]) / (Δt₂ * Δt₁) + if eltype(out) <: Number + out[i₀] = Δ * (t - A.t[i₁]) * (t - A.t[i₂]) / (Δt₀ * Δt₂) + out[i₁] = -Δ * (t - A.t[i₀]) * (t - A.t[i₂]) / (Δt₀ * Δt₁) + out[i₂] = Δ * (t - A.t[i₀]) * (t - A.t[i₁]) / (Δt₂ * Δt₁) + else + @. out[i₀] = Δ * (t - A.t[i₁]) * (t - A.t[i₂]) / (Δt₀ * Δt₂) + @. out[i₁] = -Δ * (t - A.t[i₀]) * (t - A.t[i₂]) / (Δt₀ * Δt₁) + @. out[i₂] = Δ * (t - A.t[i₀]) * (t - A.t[i₁]) / (Δt₂ * Δt₁) + end out end @@ -90,7 +111,7 @@ function ChainRulesCore.rrule(::typeof(_interpolate), t::Number) deriv = derivative(A, t) function interpolate_pullback(Δ) - (NoTangent(), Tangent{typeof(A)}(; u = u_tangent(A, t, Δ)), deriv * Δ) + (NoTangent(), Tangent{typeof(A)}(; u = u_tangent(A, t, Δ)), sum(deriv .* Δ)) end return _interpolate(A, t), interpolate_pullback end @@ -100,4 +121,8 @@ function ChainRulesCore.frule((_, _, Δt), ::typeof(_interpolate), A::AbstractIn return _interpolate(A, t), derivative(A, t) * Δt end +function ChainRulesCore.frule((_, Δt), A::AbstractInterpolation, t::Number) + return A(t), derivative(A, t) * Δt +end + end # module diff --git a/joss/paper.md b/joss/paper.md index 66d15aa0..a2c77cd4 100644 --- a/joss/paper.md +++ b/joss/paper.md @@ -8,22 +8,22 @@ authors: orcid: 0000-0003-0785-3586 corresponding: true affiliation: 1 - - name: Christopher Rackauckas - orcid: 0000-0001-5850-0663 - affiliation: "1, 2, 3" - - name: Shubham Maddhashiya - affiliation: 3 - name: Bart de Koning orcid: 0009-0005-6134-6608 - affiliation: 4 + affiliation: 2 + - name: Shubham Maddhashiya + affiliation: 3 + - name: Christopher Rackauckas + orcid: 0000-0001-5850-0663 + affiliation: "1, 3, 4" affiliations: - name: JuliaHub index: 1 - - name: Massachusetts Institute of Technology + - name: Deltares index: 2 - name: Pumas-AI index: 3 - - name: Deltares + - name: Massachusetts Institute of Technology index: 4 date: 6 June 2024 bibliography: paper.bib @@ -31,12 +31,12 @@ bibliography: paper.bib # Summary -Interpolations are used to estimate values between known data points using an approximate continuous function.DataInterpolations.jl is a Julia [@Bezanson2017] package containing 1D implementations of some of the most commonly used interpolation functions. These include: +Interpolations are used to estimate values between known data points using an approximate continuous function. DataInterpolations.jl is a Julia [@Bezanson2017] package containing 1D implementations of some of the most commonly used interpolation functions. These include: - Constant Interpolation - Linear Interpolation - Quadratic Interpolation - - Lagrange Interpolation [@lagrange] + - Lagrange Interpolation [@lagrange1898lectures] - Quadratic Splines - Cubic Splines [@Schoenberg1988] - Akima Splines [@10.1145/321607.321609] @@ -46,7 +46,7 @@ Interpolations are used to estimate values between known data points using an ap - B-Splines [@Curry1988] [@DEBOOR197250] - Regression based B-Splines -and a continually growing list. Along with these, the package also has methods to fit parameterized curves with the data points and Tikhonov regularization [@Tikhonov1943OnTS] [@amt-14-7909-2021] for obtaining smooth curves. The package also provides functionality to compute integrals and derivatives upto second order for those interpolations methods. It is also automatic differentiation friendly. It can also be used symbolically with Symbolics.jl [@gowda2021high] and plugged into models defined using ModelingToolkit.jl [@ma2021modelingtoolkit]. +and a continually growing list. Along with these, the package also has methods to fit parameterized curves with the data points and Tikhonov regularization [@Tikhonov1943OnTS] [@amt-14-7909-2021] for obtaining smooth curves. The package also provides functionality to compute integrals and derivatives upto second order for those interpolations methods. It is also automatic differentiation friendly. It can also be used symbolically with Symbolics.jl [@10.1145/3511528.3511535] and plugged into models defined using ModelingToolkit.jl [@ma2021modelingtoolkit]. # Statement of need @@ -54,7 +54,11 @@ Interpolations are a very important component of many modeling workflows. Often, # Example -The following tutorials in the documentation [1](https://docs.sciml.ai/DataInterpolations/stable/methods/) provides how to define each of the interpolation methods and compute the value at any point. [2](https://docs.sciml.ai/DataInterpolations/stable/interface/) provides explanation for using the interface and interpolated objects for evaluating at any point, computing the derivative at any point and computing the integral between any two points. [3](https://docs.sciml.ai/DataInterpolations/stable/symbolics/) provides how to use interpolation objects with Symbolics.jl and ModelingToolkit.jl. +The following tutorials are provided in the documentation: + + - [Tutorial 1](https://docs.sciml.ai/DataInterpolations/stable/methods/) provides how to define each of the interpolation methods and compute the value at any point. + - [Tutorial 2](https://docs.sciml.ai/DataInterpolations/stable/interface/) provides explanation for using the interface and interpolated objects for evaluating at any point, computing the derivative at any point and computing the integral between any two points. + - [Tutorial 3](https://docs.sciml.ai/DataInterpolations/stable/symbolics/) provides how to use interpolation objects with Symbolics.jl and ModelingToolkit.jl. A simple demonstration here: diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 5f035e62..393f1160 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -2,7 +2,7 @@ module DataInterpolations ### Interface Functionality -abstract type AbstractInterpolation{T} end +abstract type AbstractInterpolation{T, N} end using LinearAlgebra, RecipesBase using PrettyTables @@ -92,8 +92,8 @@ export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation -struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}} <: - AbstractInterpolation{T} +struct RegularizationSmooth{uType, tType, T, T2, N, ITP <: AbstractInterpolation{T, N}} <: + AbstractInterpolation{T, N} u::uType û::uType t::tType @@ -116,7 +116,8 @@ struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T} alg, Aitp, extrapolate) - new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), eltype(u), typeof(λ), N, typeof(Aitp)}( u, û, t, @@ -143,8 +144,9 @@ struct CurvefitCache{ lbType, algType, pminType, - T -} <: AbstractInterpolation{T} + T, + N +} <: AbstractInterpolation{T, N} u::uType t::tType m::mType # model type @@ -155,9 +157,10 @@ struct CurvefitCache{ pmin::pminType # optimized params extrapolate::Bool function CurvefitCache(u, t, m, p0, ub, lb, alg, pmin, extrapolate) + N = get_output_dim(u) new{typeof(u), typeof(t), typeof(m), typeof(p0), typeof(ub), typeof(lb), - typeof(alg), typeof(pmin), eltype(u)}(u, + typeof(alg), typeof(pmin), eltype(u), N}(u, t, m, p0, diff --git a/src/derivatives.jl b/src/derivatives.jl index c0835060..f7a2b0e1 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -114,12 +114,12 @@ function _derivative(A::ConstantInterpolation, t::Number, iguess) return zero(first(A.u)) end -function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number) +function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) end -function _derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number) +function _derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1] end @@ -153,19 +153,43 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num n = length(A.t) scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) t_ = A.p[idx] + (t - A.t[idx]) * scale - N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N - spline_coefficients!(N, A.d - 1, A.k, t_) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc + spline_coefficients!(sc, A.d - 1, A.k, t_) ducum = zero(eltype(A.u)) if t == A.t[1] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) else for i in 1:(n - 1) - ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) + ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end ducum * A.d * scale end +function _derivative( + A::BSplineInterpolation{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N} + # change t into param [0 1] + ax_u = axes(A.u)[1:(end - 1)] + t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...) + t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...) + idx = get_idx(A, t, iguess) + n = length(A.t) + scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) + t_ = A.p[idx] + (t - A.t[idx]) * scale + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc + spline_coefficients!(sc, A.d - 1, A.k, t_) + ducum = zeros(size(A.u)[1:(end - 1)]...) + if t == A.t[1] + ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2]) + else + for i in 1:(n - 1) + ducum = ducum + + sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / + (A.k[i + A.d + 1] - A.k[i + 1]) + end + end + ducum * A.d * scale +end # BSpline Curve Approx function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess) # change t into param [0 1] @@ -174,19 +198,42 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig idx = get_idx(A, t, iguess) scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) t_ = A.p[idx] + (t - A.t[idx]) * scale - N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N - spline_coefficients!(N, A.d - 1, A.k, t_) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc + spline_coefficients!(sc, A.d - 1, A.k, t_) ducum = zero(eltype(A.u)) if t == A.t[1] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) else for i in 1:(A.h - 1) - ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) + ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end ducum * A.d * scale end +function _derivative( + A::BSplineApprox{<:AbstractArray{<:Number, N}}, t::Number, iguess) where {N} + # change t into param [0 1] + ax_u = axes(A.u)[1:(end - 1)] + t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...) + t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...) + idx = get_idx(A, t, iguess) + scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx]) + t_ = A.p[idx] + (t - A.t[idx]) * scale + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc + spline_coefficients!(sc, A.d - 1, A.k, t_) + ducum = zeros(size(A.u)[1:(end - 1)]...) + if t == A.t[1] + ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2]) + else + for i in 1:(A.h - 1) + ducum = ducum + + sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / + (A.k[i + A.d + 1] - A.k[i + 1]) + end + end + ducum * A.d * scale +end # Cubic Hermite Spline function _derivative( A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 4cddf7d3..33621f1a 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -1,4 +1,4 @@ -abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end +abstract type AbstractIntegralInverseInterpolation{T, N} <: AbstractInterpolation{T, N} end """ invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation @@ -33,15 +33,16 @@ Can be easily constructed with `invert_integral(A::LinearInterpolation{<:Abstrac - `t` : Given by `A.I` (the cumulative integral of `A`) - `A` : The `LinearInterpolation` object """ -struct LinearInterpolationIntInv{uType, tType, itpType, T} <: - AbstractIntegralInverseInterpolation{T} +struct LinearInterpolationIntInv{uType, tType, itpType, T, N} <: + AbstractIntegralInverseInterpolation{T, N} u::uType t::tType extrapolate::Bool iguesser::Guesser{tType} itp::itpType function LinearInterpolationIntInv(u, t, A) - new{typeof(u), typeof(t), typeof(A), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(A), eltype(u), N}( u, t, A.extrapolate, Guesser(t), A) end end @@ -79,15 +80,16 @@ Can be easily constructed with `invert_integral(A::ConstantInterpolation{<:Abstr - `t` : Given by `A.I` (the cumulative integral of `A`) - `A` : The `ConstantInterpolation` object """ -struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: - AbstractIntegralInverseInterpolation{T} +struct ConstantInterpolationIntInv{uType, tType, itpType, T, N} <: + AbstractIntegralInverseInterpolation{T, N} u::uType t::tType extrapolate::Bool iguesser::Guesser{tType} itp::itpType function ConstantInterpolationIntInv(u, t, A) - new{typeof(u), typeof(t), typeof(A), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(A), eltype(u), N}( u, t, A.extrapolate, Guesser(t), A ) end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 309a17a6..f15b4154 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -19,7 +19,7 @@ Extrapolation extends the last linear polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T} +struct LinearInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType @@ -30,7 +30,8 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati linear_lookup::Bool function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u), N}( u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -66,7 +67,8 @@ Extrapolation extends the last quadratic polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T} +struct QuadraticInterpolation{uType, tType, IType, pType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType I::IType @@ -81,7 +83,8 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u), N}( u, t, I, p, mode, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -116,8 +119,8 @@ It is the method of interpolation using Lagrange polynomials of (k-1)th order pa - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. """ -struct LagrangeInterpolation{uType, tType, T, bcacheType} <: - AbstractInterpolation{T} +struct LagrangeInterpolation{uType, tType, T, bcacheType, N} <: + AbstractInterpolation{T, N} u::uType t::tType n::Int @@ -129,7 +132,8 @@ struct LagrangeInterpolation{uType, tType, T, bcacheType} <: bcache = zeros(eltype(u[1]), n + 1) idxs = zeros(Int, n + 1) fill!(bcache, NaN) - new{typeof(u), typeof(t), eltype(u), typeof(bcache)}(u, + N = get_output_dim(u) + new{typeof(u), typeof(t), eltype(u), typeof(bcache), N}(u, t, n, bcache, @@ -169,8 +173,8 @@ Extrapolation extends the last cubic polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: - AbstractInterpolation{T} +struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType I::IType @@ -184,8 +188,9 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <: function AkimaInterpolation( u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) + N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), - typeof(d), eltype(u)}(u, + typeof(d), eltype(u), N}(u, t, I, b, @@ -251,7 +256,7 @@ Extrapolation extends the last constant polynomial at the end points on each sid for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} +struct ConstantInterpolation{uType, tType, IType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType @@ -264,7 +269,8 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T} function ConstantInterpolation( u, t, I, dir, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(I), eltype(u), N}( u, t, I, nothing, dir, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -299,30 +305,31 @@ Extrapolation extends the last quadratic polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct QuadraticSpline{uType, tType, IType, pType, kType, cType, NType, T} <: - AbstractInterpolation{T} +struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType I::IType p::QuadraticSplineParameterCache{pType} k::kType # knot vector c::cType # 1D B-spline control points - N::NType # Spline coefficients (preallocated memory) + sc::scType # Spline coefficients (preallocated memory) extrapolate::Bool iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticSpline( - u, t, I, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) + N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k), - typeof(c), typeof(N), eltype(u)}(u, + typeof(c), typeof(sc), eltype(u), N}(u, t, I, p, k, c, - N, + sc, extrapolate, Guesser(t), cache_parameters, @@ -338,16 +345,16 @@ function QuadraticSpline( u, t = munge_data(u, t) n = length(t) - dtype_N = typeof(t[1] / t[1]) - N = zeros(dtype_N, n) - k, A = quadratic_spline_params(t, N) + dtype_sc = typeof(t[1] / t[1]) + sc = zeros(dtype_sc, n) + k, A = quadratic_spline_params(t, sc) c = A \ u - p = QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) + p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) + u, t, nothing, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) + QuadraticSpline(u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) end function QuadraticSpline( @@ -399,7 +406,8 @@ Second derivative on both ends are zero, which are also called "natural" boundar for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T} +struct CubicSpline{uType, tType, IType, pType, hType, zType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType I::IType @@ -412,7 +420,9 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter linear_lookup::Bool function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), + typeof(h), typeof(z), eltype(u), N}( u, t, I, @@ -457,6 +467,40 @@ function CubicSpline(u::uType, CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) end +function CubicSpline(u::uType, + t; + extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractArray{T, N}} where {T, N} + u, t = munge_data(u, t) + n = length(t) - 1 + h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) + dl = vcat(h[2:n], zero(eltype(h))) + d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) + du = vcat(zero(eltype(h)), h[3:(n + 1)]) + tA = Tridiagonal(dl, d_tmp, du) + + # zero for element type of d, which we don't know yet + ax = axes(u)[1:(end - 1)] + typed_zero = zero(6(u[ax..., begin + 2] - u[ax..., begin + 1]) / h[begin + 2] - + 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) + + h_ = reshape(h, ones(Int64, N - 1)..., :) + ax_h = axes(h_)[1:(end - 1)] + d = 6 * ((u[ax..., 3:(n + 1)] - u[ax..., 2:n]) ./ h_[ax_h..., 3:(n + 1)]) - + 6 * ((u[ax..., 2:n] - u[ax..., 1:(n - 1)]) ./ h_[ax_h..., 2:n]) + d = cat(typed_zero, d, typed_zero; dims = ndims(d)) + d_reshaped = reshape(d, prod(size(d)[1:(end - 1)]), :) + z = (tA \ d_reshaped')' + z = reshape(z, size(u)...) + linear_lookup = seems_linear(assume_linear_t, t) + p = CubicSplineParameterCache(u, h, z, cache_parameters) + A = CubicSpline( + u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + I = cumulative_integral(A, cache_parameters) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) +end + function CubicSpline( u::uType, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: @@ -505,15 +549,15 @@ Extrapolation is a constant polynomial of the end points on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: - AbstractInterpolation{T} +struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType d::Int # degree p::pType # params vector k::kType # knot vector c::cType # control points - N::NType # Spline coefficients (preallocated memory) + sc::scType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol extrapolate::Bool @@ -525,19 +569,21 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: p, k, c, - N, + sc, pVecType, knotVecType, extrapolate, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u, + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}( + u, t, d, p, k, c, - N, + sc, pVecType, knotVecType, extrapolate, @@ -548,7 +594,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <: end function BSplineInterpolation( - u, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) + u::AbstractVector, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) u, t = munge_data(u, t) n = length(t) n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.") @@ -608,12 +654,85 @@ function BSplineInterpolation( end end # control points - N = zeros(eltype(t), n, n) - spline_coefficients!(N, d, k, p) - c = vec(N \ u[:, :]) - N = zeros(eltype(t), n) + sc = zeros(eltype(t), n, n) + spline_coefficients!(sc, d, k, p) + c = vec(sc \ u[:, :]) + sc = zeros(eltype(t), n) BSplineInterpolation( - u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) +end + +function BSplineInterpolation( + u::AbstractArray{T, N}, t, d, pVecType, knotVecType; extrapolate = false, + assume_linear_t = 1e-2) where {T, N} + u, t = munge_data(u, t) + n = length(t) + n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.") + s = zero(eltype(u)) + p = zero(t) + k = zeros(eltype(t), n + d + 1) + l = zeros(eltype(u), n - 1) + p[1] = zero(eltype(t)) + p[end] = one(eltype(t)) + + ax_u = axes(u)[1:(end - 1)] + + for i in 2:n + s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2)) + l[i - 1] = s + end + if pVecType == :Uniform + for i in 2:(n - 1) + p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1) + end + elseif pVecType == :ArcLen + for i in 2:(n - 1) + p[i] = p[1] + l[i - 1] / s * (p[end] - p[1]) + end + end + + lidx = 1 + ridx = length(k) + while lidx <= (d + 1) && ridx >= (length(k) - d) + k[lidx] = p[1] + k[ridx] = p[end] + lidx += 1 + ridx -= 1 + end + + ps = zeros(eltype(t), n - 2) + s = zero(eltype(t)) + for i in 2:(n - 1) + s += p[i] + ps[i - 1] = s + end + + if knotVecType == :Uniform + # uniformly spaced knot vector + # this method is not recommended because, if it is used with the chord length method for global interpolation, + # the system of linear equations would be singular. + for i in (d + 2):n + k[i] = k[1] + (i - d - 1) // (n - d) * (k[end] - k[1]) + end + elseif knotVecType == :Average + # average spaced knot vector + idx = 1 + if d + 2 <= n + k[d + 2] = 1 // d * ps[d] + end + for i in (d + 3):n + k[i] = 1 // d * (ps[idx + d] - ps[idx]) + idx += 1 + end + end + # control points + sc = zeros(eltype(t), n, n) + spline_coefficients!(sc, d, k, p) + c = (sc \ reshape(u, prod(size(u)[1:(end - 1)]), :)')' + c = reshape(c, size(u)...) + sc = zeros(eltype(t), n) + BSplineInterpolation( + u, t, d, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) end """ @@ -640,8 +759,8 @@ Extrapolation is a constant polynomial of the end points on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: - AbstractInterpolation{T} +struct BSplineApprox{uType, tType, pType, kType, cType, scType, T, N} <: + AbstractInterpolation{T, N} u::uType t::tType d::Int # degree @@ -649,7 +768,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: p::pType # params vector k::kType # knot vector c::cType # control points - N::NType # Spline coefficients (preallocated memory) + sc::scType # Spline coefficients (preallocated memory) pVecType::Symbol knotVecType::Symbol extrapolate::Bool @@ -662,21 +781,23 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: p, k, c, - N, + sc, pVecType, knotVecType, extrapolate, assume_linear_t ) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u, + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u), N}( + u, t, d, h, p, k, c, - N, + sc, pVecType, knotVecType, extrapolate, @@ -687,7 +808,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <: end function BSplineApprox( - u, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) + u::AbstractVector, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2) u, t = munge_data(u, t) n = length(t) h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.") @@ -752,30 +873,125 @@ function BSplineApprox( c[1] = u[1] c[end] = u[end] q = zeros(eltype(u), n) - N = zeros(eltype(t), n, h) + sc = zeros(eltype(t), n, h) for i in 1:n - spline_coefficients!(view(N, i, :), d, k, p[i]) + spline_coefficients!(view(sc, i, :), d, k, p[i]) end for k in 2:(n - 1) - q[k] = u[k] - N[k, 1] * u[1] - N[k, h] * u[end] + q[k] = u[k] - sc[k, 1] * u[1] - sc[k, h] * u[end] end Q = Matrix{eltype(u)}(undef, h - 2, 1) for i in 2:(h - 1) s = 0.0 for k in 2:(n - 1) - s += N[k, i] * q[k] + s += sc[k, i] * q[k] end Q[i - 1] = s end - N = N[2:(end - 1), 2:(h - 1)] - M = transpose(N) * N + sc = sc[2:(end - 1), 2:(h - 1)] + M = transpose(sc) * sc P = M \ Q c[2:(end - 1)] .= vec(P) - N = zeros(eltype(t), h) + sc = zeros(eltype(t), h) BSplineApprox( - u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t) + u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) end +function BSplineApprox( + u::AbstractArray{T, N}, t, d, h, pVecType, knotVecType; extrapolate = false, + assume_linear_t = 1e-2) where {T, N} + u, t = munge_data(u, t) + n = length(t) + h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.") + s = zero(eltype(u)) + p = zero(t) + k = zeros(eltype(t), h + d + 1) + l = zeros(eltype(u), n - 1) + p[1] = zero(eltype(t)) + p[end] = one(eltype(t)) + + ax_u = axes(u)[1:(end - 1)] + + for i in 2:n + s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2)) + l[i - 1] = s + end + if pVecType == :Uniform + for i in 2:(n - 1) + p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1) + end + elseif pVecType == :ArcLen + for i in 2:(n - 1) + p[i] = p[1] + l[i - 1] / s * (p[end] - p[1]) + end + end + + lidx = 1 + ridx = length(k) + while lidx <= (d + 1) && ridx >= (length(k) - d) + k[lidx] = p[1] + k[ridx] = p[end] + lidx += 1 + ridx -= 1 + end + + ps = zeros(eltype(t), n - 2) + s = zero(eltype(t)) + for i in 2:(n - 1) + s += p[i] + ps[i - 1] = s + end + + if knotVecType == :Uniform + # uniformly spaced knot vector + # this method is not recommended because, if it is used with the chord length method for global interpolation, + # the system of linear equations would be singular. + for i in (d + 2):h + k[i] = k[1] + (i - d - 1) // (h - d) * (k[end] - k[1]) + end + elseif knotVecType == :Average + # NOTE: verify that average method can be applied when size of k is less than size of p + # average spaced knot vector + idx = 1 + if d + 2 <= h + k[d + 2] = 1 // d * ps[d] + end + for i in (d + 3):h + k[i] = 1 // d * (ps[idx + d] - ps[idx]) + idx += 1 + end + end + # control points + c = zeros(eltype(u), size(u)[1:(end - 1)]..., h) + c[ax_u..., 1] = u[ax_u..., 1] + c[ax_u..., end] = u[ax_u..., end] + q = zeros(eltype(u), size(u)[1:(end - 1)]..., n) + sc = zeros(eltype(t), n, h) + for i in 1:n + spline_coefficients!(view(sc, i, :), d, k, p[i]) + end + for k in 2:(n - 1) + q[ax_u..., k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - + sc[k, h] * u[ax_u..., end] + end + Q = Array{eltype(u), N}(undef, size(u)[1:(end - 1)]..., h - 2) + for i in 2:(h - 1) + s = zeros(eltype(sc), size(u)[1:(end - 1)]...) + for k in 2:(n - 1) + s = s + sc[k, i] * q[ax_u..., k] + end + Q[ax_u..., i - 1] = s + end + sc = sc[2:(end - 1), 2:(h - 1)] + M = transpose(sc) * sc + Q = reshape(Q, prod(size(u)[1:(end - 1)]), :) + P = (M \ Q')' + P = reshape(P, size(u)[1:(end - 1)]..., :) + c[ax_u..., 2:(end - 1)] = P + sc = zeros(eltype(t), h) + BSplineApprox( + u, t, d, h, p, k, c, sc, pVecType, knotVecType, extrapolate, assume_linear_t) +end """ CubicHermiteSpline(du, u, t; extrapolate = false, cache_parameters = false) @@ -796,7 +1012,8 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T} +struct CubicHermiteSpline{uType, tType, IType, duType, pType, T, N} <: + AbstractInterpolation{T, N} du::duType u::uType t::tType @@ -809,7 +1026,8 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte function CubicHermiteSpline( du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}( + N = get_output_dim(u) + new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u), N}( du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end @@ -875,8 +1093,8 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: - AbstractInterpolation{T} +struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T, N} <: + AbstractInterpolation{T, N} ddu::dduType du::duType u::uType @@ -890,8 +1108,9 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <: function QuinticHermiteSpline( ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) + N = get_output_dim(u) new{typeof(u), typeof(t), typeof(I), typeof(du), - typeof(ddu), typeof(p.c₁), eltype(u)}( + typeof(ddu), typeof(p.c₁), eltype(u), N}( ddu, du, u, t, I, p, extrapolate, Guesser(t), cache_parameters, linear_lookup) end end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index eb499353..f47f9e3c 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -33,11 +33,12 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues val end -function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) +function _interpolate(A::LinearInterpolation{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess) Δt = t - A.t[idx] slope = get_parameters(A, idx) - return A.u[:, idx] + slope * Δt + ax = axes(A.u)[1:(end - 1)] + return A.u[ax..., idx] + slope * Δt end # Quadratic Interpolation @@ -165,6 +166,18 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) I + C + D end +function _interpolate(A::CubicSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} + idx = get_idx(A, t, iguess) + Δt₁ = t - A.t[idx] + Δt₂ = A.t[idx + 1] - t + ax = axes(A.z)[1:(end - 1)] + I = (A.z[ax..., idx] * Δt₂^3 + A.z[ax..., idx + 1] * Δt₁^3) / (6A.h[idx + 1]) + c₁, c₂ = get_parameters(A, idx) + C = c₁ * Δt₁ + D = c₂ * Δt₂ + I + C + D +end + # BSpline Curve Interpolation function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, @@ -175,11 +188,30 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, idx = get_idx(A, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) n = length(A.t) - N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N - nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc + nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t) ucum = zero(eltype(A.u)) for i in nonzero_coefficient_idxs - ucum += N[i] * A.c[i] + ucum += sc[i] * A.c[i] + end + ucum +end + +function _interpolate(A::BSplineInterpolation{<:AbstractArray{T, N}}, + t::Number, + iguess) where {T <: Number, N} + ax_u = axes(A.u)[1:(end - 1)] + t < A.t[1] && return A.u[ax_u..., 1] + t > A.t[end] && return A.u[ax_u..., end] + # change t into param [0 1] + idx = get_idx(A, t, iguess) + t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) + n = length(A.t) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc + nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t) + ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...) + for i in nonzero_coefficient_idxs + ucum = ucum + (sc[i] * A.c[ax_u..., i]) end ucum end @@ -191,11 +223,28 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i # change t into param [0 1] idx = get_idx(A, t, iguess) t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) - N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N - nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc + nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t) ucum = zero(eltype(A.u)) for i in nonzero_coefficient_idxs - ucum += N[i] * A.c[i] + ucum += sc[i] * A.c[i] + end + ucum +end + +function _interpolate( + A::BSplineApprox{<:AbstractArray{T, N}}, t::Number, iguess) where {T <: Number, N} + ax_u = axes(A.u)[1:(end - 1)] + t < A.t[1] && return A.u[ax_u..., 1] + t > A.t[end] && return A.u[ax_u..., end] + # change t into param [0 1] + idx = get_idx(A, t, iguess) + t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx]) + sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc + nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t) + ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...) + for i in nonzero_coefficient_idxs + ucum = ucum + (sc[i] * A.c[ax_u..., i]) end ucum end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index cb205a1a..f085605b 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -59,7 +59,7 @@ function spline_coefficients!(N, d, k, u::AbstractVector) return nothing end -function quadratic_spline_params(t::AbstractVector, N::AbstractVector) +function quadratic_spline_params(t::AbstractVector, sc::AbstractVector) # Create knot vector # Don't use x[end-1] as knot to match number of degrees of freedom with data @@ -72,17 +72,17 @@ function quadratic_spline_params(t::AbstractVector, N::AbstractVector) # - A consists of basis function evaulations in t # - c are the 1D control points n = length(t) - dtype_N = typeof(t[1] / t[1]) + dtype_sc = typeof(t[1] / t[1]) - diag = Vector{dtype_N}(undef, n) - diag_hi = Vector{dtype_N}(undef, n - 1) - diag_lo = Vector{dtype_N}(undef, n - 1) + diag = Vector{dtype_sc}(undef, n) + diag_hi = Vector{dtype_sc}(undef, n - 1) + diag_lo = Vector{dtype_sc}(undef, n - 1) for (i, tᵢ) in enumerate(t) - spline_coefficients!(N, 2, k, tᵢ) - diag[i] = N[i] - (i > 1) && (diag_lo[i - 1] = N[i - 1]) - (i < n) && (diag_hi[i] = N[i + 1]) + spline_coefficients!(sc, 2, k, tᵢ) + diag[i] = sc[i] + (i > 1) && (diag_lo[i - 1] = sc[i - 1]) + (i < n) && (diag_hi[i] = sc[i + 1]) end A = Tridiagonal(diag_lo, diag, diag_hi) @@ -90,6 +90,19 @@ function quadratic_spline_params(t::AbstractVector, N::AbstractVector) return k, A end +# Get Output Dimension for parameterizing AbstractInterpolations +function get_output_dim(u::AbstractVector{<:Number}) + return (1,) +end + +function get_output_dim(u::AbstractVector) + return (length(first(u)),) +end + +function get_output_dim(u::AbstractArray) + return size(u)[1:(end - 1)] +end + # helper function for data manipulation function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) return u, t @@ -125,6 +138,21 @@ function munge_data(U::StridedMatrix, t::AbstractVector) return U, t end +function munge_data(U::AbstractArray{T, N}, t) where {T, N} + TU = Base.nonmissingtype(eltype(U)) + Tt = Base.nonmissingtype(eltype(t)) + @assert length(t) == size(U, ndims(U)) + ax = axes(U)[1:(end - 1)] + non_missing_indices = collect( + i for i in 1:length(t) + if !any(ismissing, U[ax..., i]) && !ismissing(t[i]) + ) + U = cat([TU.(U[ax..., i]) for i in non_missing_indices]...; dims = ndims(U)) + t = Tt.([t[i] for i in non_missing_indices]) + + return U, t +end + seems_linear(assume_linear_t::Bool, _) = assume_linear_t seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t) @@ -192,7 +220,7 @@ function get_parameters(A::QuadraticSpline, idx) if A.cache_parameters A.p.α[idx], A.p.β[idx] else - quadratic_spline_parameters(A.u, A.t, A.k, A.c, A.N, idx) + quadratic_spline_parameters(A.u, A.t, A.k, A.c, A.sc, idx) end end diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index fcf4b56d..2c787725 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -18,9 +18,11 @@ function safe_diff(b, a::T) where {T} b == a ? zero(T) : b - a end -function linear_interpolation_parameters(u::AbstractArray{T}, t, idx) where {T} - Δu = if u isa AbstractMatrix - [safe_diff(u[j, idx + 1], u[j, idx]) for j in 1:size(u)[1]] +function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where {T, N} + Δu = if N > 1 + ax = axes(u) + safe_diff.( + u[ax[1:(end - 1)]..., (idx + 1)], u[ax[1:(end - 1)]..., idx]) else safe_diff(u[idx + 1], u[idx]) end @@ -77,25 +79,25 @@ struct QuadraticSplineParameterCache{pType} β::pType end -function QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) +function QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) if cache_parameters parameters = quadratic_spline_parameters.( - Ref(u), Ref(t), Ref(k), Ref(c), Ref(N), 1:(length(t) - 1)) + Ref(u), Ref(t), Ref(k), Ref(c), Ref(sc), 1:(length(t) - 1)) α, β = collect.(eachrow(stack(collect.(parameters)))) QuadraticSplineParameterCache(α, β) else # Compute parameters once to infer types - α, β = quadratic_spline_parameters(u, t, k, c, N, 1) + α, β = quadratic_spline_parameters(u, t, k, c, sc, 1) QuadraticSplineParameterCache(typeof(α)[], typeof(β)[]) end end -function quadratic_spline_parameters(u, t, k, c, N, idx) +function quadratic_spline_parameters(u, t, k, c, sc, idx) tᵢ₊ = (t[idx] + t[idx + 1]) / 2 - nonzero_coefficient_idxs = spline_coefficients!(N, 2, k, tᵢ₊) + nonzero_coefficient_idxs = spline_coefficients!(sc, 2, k, tᵢ₊) uᵢ₊ = zero(first(u)) for j in nonzero_coefficient_idxs - uᵢ₊ += N[j] * c[j] + uᵢ₊ += sc[j] * c[j] end α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊ β = 4 * (uᵢ₊ - u[idx]) - (u[idx + 1] - u[idx]) @@ -121,12 +123,19 @@ function CubicSplineParameterCache(u, h, z, cache_parameters) end end -function cubic_spline_parameters(u, h, z, idx) +function cubic_spline_parameters(u::AbstractVector, h, z, idx) c₁ = (u[idx + 1] / h[idx + 1] - z[idx + 1] * h[idx + 1] / 6) c₂ = (u[idx] / h[idx + 1] - z[idx] * h[idx + 1] / 6) return c₁, c₂ end +function cubic_spline_parameters(u::AbstractArray, h, z, idx) + ax = axes(u)[1:(end - 1)] + c₁ = (u[ax..., idx + 1] / h[idx + 1] - z[ax..., idx + 1] * h[idx + 1] / 6) + c₂ = (u[ax..., idx] / h[idx + 1] - z[ax..., idx] * h[idx + 1] / 6) + return c₁, c₂ +end + struct CubicHermiteParameterCache{pType} c₁::pType c₂::pType diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 69d3700b..37595a4f 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -133,9 +133,11 @@ end @testset "Constant Interpolation" begin u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0] - t = collect(0.0:10.0) + t = collect(0.0:11.0) A = ConstantInterpolation(u, t) - @test all(derivative.(Ref(A), t) .== 0.0) + t2 = collect(0.0:10.0) + @test all(isnan, derivative.(Ref(A), t)) + @test all(derivative.(Ref(A), t2 .+ 0.1) .== 0.0) end @testset "Quadratic Spline" begin @@ -184,6 +186,44 @@ end :Uniform, :Uniform], name = "BSpline Approx (Uniform, Uniform)") + + f3d(t) = [sin(t) cos(t); + 0.0 cos(2t)] + + t3d = 0.1:0.1:1.0 |> collect + u3d = cat(f3d.(t3d)...; dims = 3) + test_derivatives(BSplineInterpolation; + args = [u3d, t3d, + 2, + :Uniform, + :Uniform], + name = "BSpline Interpolation (Uniform, Uniform): AbstractArray" + ) + + test_derivatives(BSplineInterpolation; + args = [u3d, t3d, + 2, + :ArcLen, + :Average], + name = "BSpline Interpolation (Arclen, Average): AbstractArray" + ) + + test_derivatives(BSplineApprox; + args = [u3d, t3d, + 3, + 4, + :Uniform, + :Uniform], + name = "BSpline Approx (Uniform, Uniform): AbstractArray") + + test_derivatives(BSplineApprox; + args = [u3d, t3d, + 3, + 4, + :ArcLen, + :Average], + name = "BSpline Approx (Arclen, Average): AbstractArray" + ) end @testset "Cubic Hermite Spline" begin diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 10b5ac0c..bbeb72f4 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -598,15 +598,35 @@ end A = CubicSpline(u, t) @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0) + + @testset "AbstractMatrix" begin + t = 0.1:0.1:1.0 + u = [sin.(t) cos.(t)]' |> collect + c = CubicSpline(u, t) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, c.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) + end + @testset "AbstractArray{T, 3}" begin + f3d(t) = [sin(t) cos(t); + 0.0 cos(2t)] + t = 0.1:0.1:1.0 + u3d = f3d.(t) + c = CubicSpline(u3d, t) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, c.(t_test)) + f_test = reduce(hcat, f3d.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + end end @testset "BSplines" begin # BSpline Interpolation and Approximation - t = [0, 62.25, 109.66, 162.66, 205.8, 252.3] - u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] - @testset "BSplineInterpolation" begin + t = [0, 62.25, 109.66, 162.66, 205.8, 252.3] + u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] test_interpolation_type(BSplineInterpolation) A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) @@ -642,10 +662,43 @@ end A = BSplineInterpolation(u, t, 2, :ArcLen, :Average) @test_throws DataInterpolations.ExtrapolationError A(-1.0) @test_throws DataInterpolations.ExtrapolationError A(300.0) + + @testset "AbstractMatrix" begin + t = 0.1:0.1:1.0 + u2d = [sin.(t) cos.(t)]' |> collect + A = BSplineInterpolation(u2d, t, 2, :Uniform, :Uniform) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) + + A = BSplineInterpolation(u2d, t, 2, :ArcLen, :Average) + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) + end + @testset "AbstractArray{T, 3}" begin + f3d(t) = [sin(t) cos(t); + 0.0 cos(2t)] + t = 0.1:0.1:1.0 + u3d = cat(f3d.(t)..., dims = 3) + A = BSplineInterpolation(u3d, t, 2, :Uniform, :Uniform) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + f_test = reduce(hcat, f3d.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + + A = BSplineInterpolation(u3d, t, 2, :ArcLen, :Average) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + end end @testset "BSplineApprox" begin test_interpolation_type(BSplineApprox) + t = [0, 62.25, 109.66, 162.66, 205.8, 252.3] + u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178] @@ -666,6 +719,37 @@ end A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) @test_throws DataInterpolations.ExtrapolationError A(-1.0) @test_throws DataInterpolations.ExtrapolationError A(300.0) + + @testset "AbstractMatrix" begin + t = 0.1:0.1:1.0 + u2d = [sin.(t) cos.(t)]' |> collect + A = BSplineApprox(u2d, t, 2, 5, :Uniform, :Uniform) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) + + A = BSplineApprox(u2d, t, 2, 5, :ArcLen, :Average) + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-2) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-2) + end + @testset "AbstractArray{T, 3}" begin + f3d(t) = [sin(t) cos(t); + 0.0 cos(2t)] + t = 0.1:0.1:1.0 + u3d = cat(f3d.(t)..., dims = 3) + A = BSplineApprox(u3d, t, 2, 6, :Uniform, :Uniform) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + f_test = reduce(hcat, f3d.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + + A = BSplineApprox(u3d, t, 2, 7, :ArcLen, :Average) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, A.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + end end end diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 32e8de10..5bc09005 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -4,13 +4,13 @@ using Zygote function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name::String) func = method(args..., u, t, args_after...; kwargs..., extrapolate = true) - (; u, t) = func trange = collect(range(minimum(t) - 5.0, maximum(t) + 5.0, step = 0.1)) trange_exclude = filter(x -> !in(x, t), trange) @testset "$name, derivatives w.r.t. input" begin for _t in trange_exclude adiff = DataInterpolations.derivative(func, _t) - zdiff = only(Zygote.gradient(func, _t)) + zdiff = u isa AbstractVector{<:Real} ? only(Zygote.gradient(func, _t)) : + only(Zygote.jacobian(func, _t)) isnothing(zdiff) && (zdiff = 0.0) @test adiff ≈ zdiff end @@ -19,15 +19,26 @@ function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name @testset "$name, derivatives w.r.t. u" begin function f(u) A = method(args..., u, t, args_after...; kwargs..., extrapolate = true) - out = zero(eltype(u)) + out = if u isa AbstractVector{<:Real} + zero(eltype(u)) + elseif u isa AbstractMatrix + zero(u[:, 1]) + else + zero(u[1]) + end + for _t in trange out += A(_t) end out end - zgrad = only(Zygote.gradient(f, u)) - fgrad = ForwardDiff.gradient(f, u) - @test zgrad ≈ fgrad + zgrad, fgrad = if u isa AbstractVector{<:Real} + Zygote.gradient(f, u), ForwardDiff.gradient(f, u) + elseif u isa AbstractMatrix + Zygote.jacobian(f, u), ForwardDiff.jacobian(f, u) + else + Zygote.jacobian(f, u), ForwardDiff.jacobian(f, hcat(u...)) + end end end end @@ -48,7 +59,15 @@ end @testset "Constant Interpolation" begin u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0] t = collect(0.0:10.0) - test_zygote(ConstantInterpolation, u, t; name = "Constant Interpolation") + test_zygote(ConstantInterpolation, u, t; name = "Constant Interpolation (vector)") + + t = [1.0, 4.0] + u = [1.0 2.0; 0.0 1.0; 1.0 2.0; 0.0 1.0] + test_zygote(ConstantInterpolation, u, t, name = "Constant Interpolation (matrix)") + + u = [[1.0, 2.0, 3.0, 4.0], [2.0, 3.0, 4.0, 5.0]] + test_zygote( + ConstantInterpolation, u, t, name = "Constant Interpolation (vector of vectors)") end @testset "Cubic Hermite Spline" begin