From 82a990966c289eff75d664f355b6778cd7661fd7 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 10:42:03 +0100 Subject: [PATCH 1/8] Update QuadraticSpline --- src/derivatives.jl | 8 ++-- src/integrals.jl | 7 ++-- src/interpolation_caches.jl | 77 +++++++++++++++++------------------- src/interpolation_methods.jl | 8 ++-- src/interpolation_utils.jl | 35 +++++++++++++++- src/parameter_caches.jl | 28 ++++++++----- test/derivative_tests.jl | 4 +- test/interpolation_tests.jl | 19 +++++---- test/parameter_tests.jl | 3 +- test/zygote_tests.jl | 2 +- 10 files changed, 116 insertions(+), 75 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index ff5b9850..c0835060 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -126,9 +126,11 @@ end # QuadraticSpline Interpolation function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first) - σ = get_parameters(A, idx - 1) - A.z[idx - 1] + 2σ * (t - A.t[idx - 1]) + idx = get_idx(A, t, iguess) + α, β = get_parameters(A, idx) + Δt = t - A.t[idx] + Δt_full = A.t[idx + 1] - A.t[idx] + 2α * Δt / Δt_full^2 + β / Δt_full end # CubicSpline Interpolation diff --git a/src/integrals.jl b/src/integrals.jl index 3ae893f7..5e2e5907 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -71,10 +71,11 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, end function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - Cᵢ = A.u[idx] + α, β = get_parameters(A, idx) + uᵢ = A.u[idx] Δt = t - A.t[idx] - σ = get_parameters(A, idx) - return A.z[idx] * Δt^2 / 2 + σ * Δt^3 / 3 + Cᵢ * Δt + Δt_full = (A.t[idx + 1] - A.t[idx]) + Δt * (α * Δt^2 / (3Δt_full^2) + β * Δt / (2Δt_full) + uᵢ) end function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index d7dce336..309a17a6 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -299,30 +299,30 @@ 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, tAType, dType, zType, T} <: +struct QuadraticSpline{uType, tType, IType, pType, kType, cType, NType, T} <: AbstractInterpolation{T} u::uType t::tType I::IType p::QuadraticSplineParameterCache{pType} - tA::tAType - d::dType - z::zType + k::kType # knot vector + c::cType # 1D B-spline control points + N::NType # Spline coefficients (preallocated memory) extrapolate::Bool iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticSpline( - u, t, I, p, tA, d, z, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) - new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA), - typeof(d), typeof(z), eltype(u)}(u, + new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k), + typeof(c), typeof(N), eltype(u)}(u, t, I, p, - tA, - d, - z, + k, + c, + N, extrapolate, Guesser(t), cache_parameters, @@ -336,24 +336,18 @@ function QuadraticSpline( cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector{<:Number}} u, t = munge_data(u, t) - linear_lookup = seems_linear(assume_linear_t, t) - s = length(t) - dl = ones(eltype(t), s - 1) - d_tmp = ones(eltype(t), s) - du = zeros(eltype(t), s - 1) - tA = Tridiagonal(dl, d_tmp, du) - - # zero for element type of d, which we don't know yet - typed_zero = zero(2 // 1 * (u[begin + 1] - u[begin]) / (t[begin + 1] - t[begin])) - d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) - z = tA \ d + n = length(t) + dtype_N = typeof(t[1] / t[1]) + N = zeros(dtype_N, n) + k, A = quadratic_spline_params(t, N) + c = A \ u - p = QuadraticSplineParameterCache(z, t, cache_parameters) + p = QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + QuadraticSpline(u, t, I, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) end function QuadraticSpline( @@ -361,25 +355,28 @@ function QuadraticSpline( assume_linear_t = 1e-2) where {uType <: AbstractVector} u, t = munge_data(u, t) - linear_lookup = seems_linear(assume_linear_t, t) - s = length(t) - dl = ones(eltype(t), s - 1) - d_tmp = ones(eltype(t), s) - du = zeros(eltype(t), s - 1) - tA = Tridiagonal(dl, d_tmp, du) - d_ = map( - i -> i == 1 ? zeros(eltype(t), size(u[1])) : - 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), - 1:s) - d = transpose(reshape(reduce(hcat, d_), :, s)) - z_ = reshape(transpose(tA \ d), size(u[1])..., :) - z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - p = QuadraticSplineParameterCache(z, t, cache_parameters) + n = length(t) + dtype_N = typeof(t[1] / t[1]) + N = zeros(dtype_N, n) + k, A = quadratic_spline_params(t, N) + + eltype_c_prototype = one(dtype_N) * first(u) + c = [similar(eltype_c_prototype) for _ in 1:n] + + # Assuming u contains arrays of equal shape + for j in eachindex(eltype_c_prototype) + c_dim = A \ [u_[j] for u_ in u] + for (i, c_dim_) in enumerate(c_dim) + c[i][j] = c_dim_ + end + end + + p = QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + QuadraticSpline(u, t, I, p, k, c, N, extrapolate, cache_parameters, assume_linear_t) end """ diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 697f7375..eb499353 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -147,10 +147,10 @@ end # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) - Cᵢ = A.u[idx] - Δt = t - A.t[idx] - σ = get_parameters(A, idx) - return A.z[idx] * Δt + σ * Δt^2 + Cᵢ + α, β = get_parameters(A, idx) + uᵢ = A.u[idx] + Δt_scaled = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) + Δt_scaled * (α * Δt_scaled + β) + uᵢ end # CubicSpline Interpolation diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 0df5e484..cb205a1a 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -59,6 +59,37 @@ function spline_coefficients!(N, d, k, u::AbstractVector) return nothing end +function quadratic_spline_params(t::AbstractVector, N::AbstractVector) + + # Create knot vector + # Don't use x[end-1] as knot to match number of degrees of freedom with data + k = zeros(length(t) + 3) + k[1:3] .= t[1] + k[(end - 2):end] .= t[end] + k[4:(end - 3)] .= t[2:(end - 2)] + + # Create and solve linear system Ac = u, where: + # - A consists of basis function evaulations in t + # - c are the 1D control points + n = length(t) + dtype_N = 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) + + 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]) + end + + A = Tridiagonal(diag_lo, diag, diag_hi) + + return k, A +end + # helper function for data manipulation function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) return u, t @@ -159,9 +190,9 @@ end function get_parameters(A::QuadraticSpline, idx) if A.cache_parameters - A.p.σ[idx] + A.p.α[idx], A.p.β[idx] else - quadratic_spline_parameters(A.z, A.t, idx) + quadratic_spline_parameters(A.u, A.t, A.k, A.c, A.N, idx) end end diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 0701b3a2..fcf4b56d 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -73,23 +73,33 @@ function quadratic_interpolation_parameters(u, t, idx) end struct QuadraticSplineParameterCache{pType} - σ::pType + α::pType + β::pType end -function QuadraticSplineParameterCache(z, t, cache_parameters) +function QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) if cache_parameters - σ = quadratic_spline_parameters.(Ref(z), Ref(t), 1:(length(t) - 1)) - QuadraticSplineParameterCache(σ) + parameters = quadratic_spline_parameters.( + Ref(u), Ref(t), Ref(k), Ref(c), Ref(N), 1:(length(t) - 1)) + α, β = collect.(eachrow(stack(collect.(parameters)))) + QuadraticSplineParameterCache(α, β) else # Compute parameters once to infer types - σ = quadratic_spline_parameters(z, t, 1) - QuadraticSplineParameterCache(typeof(σ)[]) + α, β = quadratic_spline_parameters(u, t, k, c, N, 1) + QuadraticSplineParameterCache(typeof(α)[], typeof(β)[]) end end -function quadratic_spline_parameters(z, t, idx) - σ = 1 // 2 * (z[idx + 1] - z[idx]) / (t[idx + 1] - t[idx]) - return σ +function quadratic_spline_parameters(u, t, k, c, N, idx) + tᵢ₊ = (t[idx] + t[idx + 1]) / 2 + nonzero_coefficient_idxs = spline_coefficients!(N, 2, k, tᵢ₊) + uᵢ₊ = zero(first(u)) + for j in nonzero_coefficient_idxs + uᵢ₊ += N[j] * c[j] + end + α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊ + β = 4 * (uᵢ₊ - u[idx]) - (u[idx + 1] - u[idx]) + return α, β end struct CubicSplineParameterCache{pType} diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index a4eb349b..69d3700b 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -259,8 +259,8 @@ end derivexpr2 = expand_derivatives(substitute(D2(A(ω)), Dict(ω => 0.5τ))) symfunc1 = Symbolics.build_function(derivexpr1, τ; expression = Val{false}) symfunc2 = Symbolics.build_function(derivexpr2, τ; expression = Val{false}) - @test symfunc1(0.5) == 0.5 * 3 - @test symfunc2(0.5) == 0.5 * 6 + @test symfunc1(0.5) == 1.5 + @test symfunc2(0.5) == -3.0 u = [0.0, 1.5, 0.0] t = [0.0, 0.5, 1.0] diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 3ebc1741..10b5ac0c 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -510,16 +510,15 @@ end A = QuadraticSpline(u, t; extrapolate = true) # Solution - P₁ = x -> (x + 1)^2 # for x ∈ [-1, 0] - P₂ = x -> 2 * x + 1 # for x ∈ [ 0, 1] + P₁ = x -> 0.5 * (x + 1) * (x + 2) for (_t, _u) in zip(t, u) @test A(_t) == _u end @test A(-2.0) == P₁(-2.0) @test A(-0.5) == P₁(-0.5) - @test A(0.7) == P₂(0.7) - @test A(2.0) == P₂(2.0) + @test A(0.7) == P₁(0.7) + @test A(2.0) == P₁(2.0) test_cached_index(A) u_ = [0.0, 1.0, 3.0]' .* ones(4) @@ -527,22 +526,22 @@ end A = QuadraticSpline(u, t; extrapolate = true) @test A(-2.0) == P₁(-2.0) * ones(4) @test A(-0.5) == P₁(-0.5) * ones(4) - @test A(0.7) == P₂(0.7) * ones(4) - @test A(2.0) == P₂(2.0) * ones(4) + @test A(0.7) == P₁(0.7) * ones(4) + @test A(2.0) == P₁(2.0) * ones(4) u = [repeat(u[i], 1, 3) for i in 1:3] A = QuadraticSpline(u, t; extrapolate = true) @test A(-2.0) == P₁(-2.0) * ones(4, 3) @test A(-0.5) == P₁(-0.5) * ones(4, 3) - @test A(0.7) == P₂(0.7) * ones(4, 3) - @test A(2.0) == P₂(2.0) * ones(4, 3) + @test A(0.7) == P₁(0.7) * ones(4, 3) + @test A(2.0) == P₁(2.0) * ones(4, 3) # Test extrapolation u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] A = QuadraticSpline(u, t; extrapolate = true) - @test A(-2.0) == 1.0 - @test A(2.0) == 5.0 + @test A(-2.0) == 0.0 + @test A(2.0) == 6.0 A = QuadraticSpline(u, t) @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0) diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 2e84b98d..af734768 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -20,7 +20,8 @@ end u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = QuadraticSpline(u, t; cache_parameters = true) - @test A.p.σ ≈ [4.0, -10.0, 13.0, -14.0] + @test A.p.α ≈ [-9.5, 3.5, -0.5, -0.5] + @test A.p.β ≈ [13.5, -5.5, 1.5, 0.5] end @testset "Cubic Spline" begin diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 1a7fc447..d706c1af 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -15,7 +15,7 @@ function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name @test adiff ≈ zdiff end end - if method ∉ [LagrangeInterpolation, BSplineInterpolation, BSplineApprox] + if method ∉ [LagrangeInterpolation, BSplineInterpolation, BSplineApprox, QuadraticSpline] @testset "$name, derivatives w.r.t. u" begin function f(u) A = method(args..., u, t, args_after...; kwargs..., extrapolate = true) From 8498a29a500e6392d8a617b92ed17fcaff1c7824 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 10:47:22 +0100 Subject: [PATCH 2/8] Remove QuadraticSpline Zygote test --- test/zygote_tests.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index d706c1af..32e8de10 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -67,12 +67,6 @@ end QuinticHermiteSpline, u, t, args = [ddu, du], name = "Quintic Hermite Spline") end -@testset "Quadratic Spline" begin - u = [1.0, 4.0, 9.0, 16.0] - t = [1.0, 2.0, 3.0, 4.0] - test_zygote(QuadraticSpline, u, t, name = "Quadratic Spline") -end - @testset "Lagrange Interpolation" begin u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] From 2cca6312629fb9241df1dfbe0762b637a52bfe6d Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 11:06:34 +0100 Subject: [PATCH 3/8] 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 From e58c72eda8c57c49a676d9aabd43be2369bfd23b Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 10:42:03 +0100 Subject: [PATCH 4/8] Update QuadraticSpline --- src/derivatives.jl | 8 ++-- src/integrals.jl | 7 ++-- src/interpolation_caches.jl | 77 +++++++++++++++++------------------- src/interpolation_methods.jl | 8 ++-- src/interpolation_utils.jl | 35 +++++++++++++++- src/parameter_caches.jl | 28 ++++++++----- test/derivative_tests.jl | 4 +- test/interpolation_tests.jl | 19 +++++---- test/parameter_tests.jl | 3 +- test/zygote_tests.jl | 2 +- 10 files changed, 116 insertions(+), 75 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 77224706..f7a2b0e1 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -126,9 +126,11 @@ end # QuadraticSpline Interpolation function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first) - σ = get_parameters(A, idx - 1) - A.z[idx - 1] + 2σ * (t - A.t[idx - 1]) + idx = get_idx(A, t, iguess) + α, β = get_parameters(A, idx) + Δt = t - A.t[idx] + Δt_full = A.t[idx + 1] - A.t[idx] + 2α * Δt / Δt_full^2 + β / Δt_full end # CubicSpline Interpolation diff --git a/src/integrals.jl b/src/integrals.jl index 3ae893f7..5e2e5907 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -71,10 +71,11 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, end function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - Cᵢ = A.u[idx] + α, β = get_parameters(A, idx) + uᵢ = A.u[idx] Δt = t - A.t[idx] - σ = get_parameters(A, idx) - return A.z[idx] * Δt^2 / 2 + σ * Δt^3 / 3 + Cᵢ * Δt + Δt_full = (A.t[idx + 1] - A.t[idx]) + Δt * (α * Δt^2 / (3Δt_full^2) + β * Δt / (2Δt_full) + uᵢ) end function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 37da05f2..9e375f33 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -305,31 +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, tAType, dType, zType, T, N} <: +struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType p::QuadraticSplineParameterCache{pType} - tA::tAType - d::dType - z::zType + k::kType # knot vector + c::cType # 1D B-spline control points + sc::scType # Spline coefficients (preallocated memory) extrapolate::Bool iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function QuadraticSpline( - u, t, I, p, tA, d, z, 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(tA), - typeof(d), typeof(z), eltype(u), N}(u, + new{typeof(u), typeof(t), typeof(I), typeof(p.α), typeof(k), + typeof(c), typeof(sc), eltype(u), N}(u, t, I, p, - tA, - d, - z, + k, + c, + sc, extrapolate, Guesser(t), cache_parameters, @@ -343,24 +343,18 @@ function QuadraticSpline( cache_parameters = false, assume_linear_t = 1e-2) where {uType <: AbstractVector{<:Number}} u, t = munge_data(u, t) - linear_lookup = seems_linear(assume_linear_t, t) - s = length(t) - dl = ones(eltype(t), s - 1) - d_tmp = ones(eltype(t), s) - du = zeros(eltype(t), s - 1) - tA = Tridiagonal(dl, d_tmp, du) - - # zero for element type of d, which we don't know yet - typed_zero = zero(2 // 1 * (u[begin + 1] - u[begin]) / (t[begin + 1] - t[begin])) - d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) - z = tA \ d + n = length(t) + dtype_sc = typeof(t[1] / t[1]) + sc = zeros(dtype_sc, n) + k, A = quadratic_spline_params(t, sc) + c = A \ u - p = QuadraticSplineParameterCache(z, t, cache_parameters) + p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + QuadraticSpline(u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) end function QuadraticSpline( @@ -368,25 +362,28 @@ function QuadraticSpline( assume_linear_t = 1e-2) where {uType <: AbstractVector} u, t = munge_data(u, t) - linear_lookup = seems_linear(assume_linear_t, t) - s = length(t) - dl = ones(eltype(t), s - 1) - d_tmp = ones(eltype(t), s) - du = zeros(eltype(t), s - 1) - tA = Tridiagonal(dl, d_tmp, du) - d_ = map( - i -> i == 1 ? zeros(eltype(t), size(u[1])) : - 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), - 1:s) - d = transpose(reshape(reduce(hcat, d_), :, s)) - z_ = reshape(transpose(tA \ d), size(u[1])..., :) - z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - p = QuadraticSplineParameterCache(z, t, cache_parameters) + n = length(t) + dtype_sc = typeof(t[1] / t[1]) + sc = zeros(dtype_sc, n) + k, A = quadratic_spline_params(t, sc) + + eltype_c_prototype = one(dtype_sc) * first(u) + c = [similar(eltype_c_prototype) for _ in 1:n] + + # Assuming u contains arrays of equal shape + for j in eachindex(eltype_c_prototype) + c_dim = A \ [u_[j] for u_ in u] + for (i, c_dim_) in enumerate(c_dim) + c[i][j] = c_dim_ + end + end + + p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters) A = QuadraticSpline( - u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) I = cumulative_integral(A, cache_parameters) - QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + QuadraticSpline(u, t, I, p, k, c, sc, extrapolate, cache_parameters, assume_linear_t) end """ diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index c429ff1f..f47f9e3c 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -148,10 +148,10 @@ end # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) - Cᵢ = A.u[idx] - Δt = t - A.t[idx] - σ = get_parameters(A, idx) - return A.z[idx] * Δt + σ * Δt^2 + Cᵢ + α, β = get_parameters(A, idx) + uᵢ = A.u[idx] + Δt_scaled = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) + Δt_scaled * (α * Δt_scaled + β) + uᵢ end # CubicSpline Interpolation diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 57302257..748cf662 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -72,6 +72,37 @@ function get_output_dim(u::AbstractArray) return size(u)[1:(end - 1)] end +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 + k = zeros(length(t) + 3) + k[1:3] .= t[1] + k[(end - 2):end] .= t[end] + k[4:(end - 3)] .= t[2:(end - 2)] + + # Create and solve linear system Ac = u, where: + # - A consists of basis function evaulations in t + # - c are the 1D control points + n = length(t) + dtype_sc = typeof(t[1] / t[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!(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) + + return k, A +end + # helper function for data manipulation function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) return u, t @@ -187,9 +218,9 @@ end function get_parameters(A::QuadraticSpline, idx) if A.cache_parameters - A.p.σ[idx] + A.p.α[idx], A.p.β[idx] else - quadratic_spline_parameters(A.z, A.t, 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 7c2d6f0b..301df9f5 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -75,23 +75,33 @@ function quadratic_interpolation_parameters(u, t, idx) end struct QuadraticSplineParameterCache{pType} - σ::pType + α::pType + β::pType end -function QuadraticSplineParameterCache(z, t, cache_parameters) +function QuadraticSplineParameterCache(u, t, k, c, N, cache_parameters) if cache_parameters - σ = quadratic_spline_parameters.(Ref(z), Ref(t), 1:(length(t) - 1)) - QuadraticSplineParameterCache(σ) + parameters = quadratic_spline_parameters.( + Ref(u), Ref(t), Ref(k), Ref(c), Ref(N), 1:(length(t) - 1)) + α, β = collect.(eachrow(stack(collect.(parameters)))) + QuadraticSplineParameterCache(α, β) else # Compute parameters once to infer types - σ = quadratic_spline_parameters(z, t, 1) - QuadraticSplineParameterCache(typeof(σ)[]) + α, β = quadratic_spline_parameters(u, t, k, c, N, 1) + QuadraticSplineParameterCache(typeof(α)[], typeof(β)[]) end end -function quadratic_spline_parameters(z, t, idx) - σ = 1 // 2 * (z[idx + 1] - z[idx]) / (t[idx + 1] - t[idx]) - return σ +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] + end + α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊ + β = 4 * (uᵢ₊ - u[idx]) - (u[idx + 1] - u[idx]) + return α, β end struct CubicSplineParameterCache{pType} diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 0acb36cb..37595a4f 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -299,8 +299,8 @@ end derivexpr2 = expand_derivatives(substitute(D2(A(ω)), Dict(ω => 0.5τ))) symfunc1 = Symbolics.build_function(derivexpr1, τ; expression = Val{false}) symfunc2 = Symbolics.build_function(derivexpr2, τ; expression = Val{false}) - @test symfunc1(0.5) == 0.5 * 3 - @test symfunc2(0.5) == 0.5 * 6 + @test symfunc1(0.5) == 1.5 + @test symfunc2(0.5) == -3.0 u = [0.0, 1.5, 0.0] t = [0.0, 0.5, 1.0] diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 9af43d37..bbeb72f4 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -510,16 +510,15 @@ end A = QuadraticSpline(u, t; extrapolate = true) # Solution - P₁ = x -> (x + 1)^2 # for x ∈ [-1, 0] - P₂ = x -> 2 * x + 1 # for x ∈ [ 0, 1] + P₁ = x -> 0.5 * (x + 1) * (x + 2) for (_t, _u) in zip(t, u) @test A(_t) == _u end @test A(-2.0) == P₁(-2.0) @test A(-0.5) == P₁(-0.5) - @test A(0.7) == P₂(0.7) - @test A(2.0) == P₂(2.0) + @test A(0.7) == P₁(0.7) + @test A(2.0) == P₁(2.0) test_cached_index(A) u_ = [0.0, 1.0, 3.0]' .* ones(4) @@ -527,22 +526,22 @@ end A = QuadraticSpline(u, t; extrapolate = true) @test A(-2.0) == P₁(-2.0) * ones(4) @test A(-0.5) == P₁(-0.5) * ones(4) - @test A(0.7) == P₂(0.7) * ones(4) - @test A(2.0) == P₂(2.0) * ones(4) + @test A(0.7) == P₁(0.7) * ones(4) + @test A(2.0) == P₁(2.0) * ones(4) u = [repeat(u[i], 1, 3) for i in 1:3] A = QuadraticSpline(u, t; extrapolate = true) @test A(-2.0) == P₁(-2.0) * ones(4, 3) @test A(-0.5) == P₁(-0.5) * ones(4, 3) - @test A(0.7) == P₂(0.7) * ones(4, 3) - @test A(2.0) == P₂(2.0) * ones(4, 3) + @test A(0.7) == P₁(0.7) * ones(4, 3) + @test A(2.0) == P₁(2.0) * ones(4, 3) # Test extrapolation u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] A = QuadraticSpline(u, t; extrapolate = true) - @test A(-2.0) == 1.0 - @test A(2.0) == 5.0 + @test A(-2.0) == 0.0 + @test A(2.0) == 6.0 A = QuadraticSpline(u, t) @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0) diff --git a/test/parameter_tests.jl b/test/parameter_tests.jl index 2e84b98d..af734768 100644 --- a/test/parameter_tests.jl +++ b/test/parameter_tests.jl @@ -20,7 +20,8 @@ end u = [1.0, 5.0, 3.0, 4.0, 4.0] t = collect(1:5) A = QuadraticSpline(u, t; cache_parameters = true) - @test A.p.σ ≈ [4.0, -10.0, 13.0, -14.0] + @test A.p.α ≈ [-9.5, 3.5, -0.5, -0.5] + @test A.p.β ≈ [13.5, -5.5, 1.5, 0.5] end @testset "Cubic Spline" begin diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 09b7be17..2e1fc400 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -15,7 +15,7 @@ function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name @test adiff ≈ zdiff end end - if method ∉ [LagrangeInterpolation, BSplineInterpolation, BSplineApprox] + if method ∉ [LagrangeInterpolation, BSplineInterpolation, BSplineApprox, QuadraticSpline] @testset "$name, derivatives w.r.t. u" begin function f(u) A = method(args..., u, t, args_after...; kwargs..., extrapolate = true) From 26467e1df6966a6f3626388c182d181290cd7a64 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 10:47:22 +0100 Subject: [PATCH 5/8] Remove QuadraticSpline Zygote test --- test/zygote_tests.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 2e1fc400..5bc09005 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -86,12 +86,6 @@ end QuinticHermiteSpline, u, t, args = [ddu, du], name = "Quintic Hermite Spline") end -@testset "Quadratic Spline" begin - u = [1.0, 4.0, 9.0, 16.0] - t = [1.0, 2.0, 3.0, 4.0] - test_zygote(QuadraticSpline, u, t, name = "Quadratic Spline") -end - @testset "Lagrange Interpolation" begin u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] From 0658d1601937732c56c5af2ca9540bb45b5d6e18 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 11:32:38 +0100 Subject: [PATCH 6/8] nit --- src/interpolation_caches.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 9e375f33..d0bfe26e 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -312,7 +312,7 @@ struct QuadraticSpline{uType, tType, IType, pType, kType, cType, scType, T, N} < I::IType p::QuadraticSplineParameterCache{pType} k::kType # knot vector - c::cType # 1D B-spline control points + c::cType # B-spline control points sc::scType # Spline coefficients (preallocated memory) extrapolate::Bool iguesser::Guesser{tType} From c2627e626a23f38930b465778fa8136d6ad5b97e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 11:36:04 +0100 Subject: [PATCH 7/8] Formatting --- test/zygote_tests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl index 5bc09005..7af735af 100644 --- a/test/zygote_tests.jl +++ b/test/zygote_tests.jl @@ -15,7 +15,8 @@ function test_zygote(method, u, t; args = [], args_after = [], kwargs = [], name @test adiff ≈ zdiff end end - if method ∉ [LagrangeInterpolation, BSplineInterpolation, BSplineApprox, QuadraticSpline] + if method ∉ + [LagrangeInterpolation, BSplineInterpolation, BSplineApprox, QuadraticSpline] @testset "$name, derivatives w.r.t. u" begin function f(u) A = method(args..., u, t, args_after...; kwargs..., extrapolate = true) From e5bd8dc73d3f45ca0f12ff617de7d095d32f4799 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 12 Nov 2024 11:41:34 +0100 Subject: [PATCH 8/8] Last nits --- src/integrals.jl | 2 +- src/interpolation_utils.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index 5e2e5907..de1b4633 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -74,7 +74,7 @@ function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, α, β = get_parameters(A, idx) uᵢ = A.u[idx] Δt = t - A.t[idx] - Δt_full = (A.t[idx + 1] - A.t[idx]) + Δt_full = A.t[idx + 1] - A.t[idx] Δt * (α * Δt^2 / (3Δt_full^2) + β * Δt / (2Δt_full) + uᵢ) end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 748cf662..2d1c6432 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -81,9 +81,9 @@ function quadratic_spline_params(t::AbstractVector, sc::AbstractVector) k[(end - 2):end] .= t[end] k[4:(end - 3)] .= t[2:(end - 2)] - # Create and solve linear system Ac = u, where: + # Create linear system Ac = u, where: # - A consists of basis function evaulations in t - # - c are the 1D control points + # - c are 1D control points n = length(t) dtype_sc = typeof(t[1] / t[1])