From e555198c26d1cb7d3606c8d20895990fac1d34b9 Mon Sep 17 00:00:00 2001 From: ArnoStrouwen Date: Sat, 23 Sep 2023 05:55:00 +0200 Subject: [PATCH] format --- .JuliaFormatter.toml | 2 + LICENSE.md | 37 +- README.md | 54 +- docs/make.jl | 12 +- docs/src/index.md | 52 +- docs/src/interface.md | 8 +- docs/src/methods.md | 70 +- ext/DataInterpolationsChainRulesCoreExt.jl | 18 +- ext/DataInterpolationsOptimExt.jl | 48 +- ...ataInterpolationsRegularizationToolsExt.jl | 100 ++- ext/DataInterpolationsSymbolicsExt.jl | 6 +- src/DataInterpolations.jl | 60 +- src/derivatives.jl | 283 ++++---- src/integrals.jl | 113 +-- src/interpolation_caches.jl | 647 ++++++++++-------- src/interpolation_methods.jl | 263 +++---- src/interpolation_utils.jl | 250 +++---- src/online.jl | 44 +- src/plot_rec.jl | 489 +++++++------ test/derivative_tests.jl | 90 +-- test/integral_tests.jl | 59 +- test/interface.jl | 12 +- test/interpolation_tests.jl | 282 ++++---- test/online_tests.jl | 4 +- test/regularization.jl | 171 +++-- test/runtests.jl | 24 +- 26 files changed, 1692 insertions(+), 1506 deletions(-) create mode 100644 .JuliaFormatter.toml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 00000000..9c793591 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,2 @@ +style = "sciml" +format_markdown = true \ No newline at end of file diff --git a/LICENSE.md b/LICENSE.md index dcdb711f..395353a6 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,40 +1,33 @@ The DataInterpolations.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2018: University of Maryland, Center for Translational Medicine. -> -> +> > Permission is hereby granted, free of charge, to any person obtaining a copy -> +> > of this software and associated documentation files (the "Software"), to deal -> +> > in the Software without restriction, including without limitation the rights -> +> > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -> +> > copies of the Software, and to permit persons to whom the Software is -> +> > furnished to do so, subject to the following conditions: -> -> -> +> > The above copyright notice and this permission notice shall be included in all -> +> > copies or substantial portions of the Software. -> -> -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -> +> > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -> +> > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -> +> > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -> +> > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -> +> > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -> +> > SOFTWARE. -> -> diff --git a/README.md b/README.md index 27b2b8b7..d9a64551 100644 --- a/README.md +++ b/README.md @@ -6,10 +6,9 @@ [![codecov](https://codecov.io/gh/SciML/DataInterpolations.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SciML/DataInterpolations.jl) [![CI](https://github.com/SciML/DataInterpolations.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/SciML/DataInterpolations.jl/actions/workflows/CI.yml) -[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![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) - 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 some methods are mixtures of regressions with interpolations (i.e. do not hit the data @@ -23,7 +22,7 @@ All interpolation objects act as functions. Thus for example, using an interpola ```julia u = rand(5) t = 0:4 -interp = LinearInterpolation(u,t) +interp = LinearInterpolation(u, t) interp(3.5) # Gives the linear interpolation value at t=3.5 ``` @@ -52,45 +51,40 @@ interp[4] # Gives the 4th value of u In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` of timepoints corresponding to `(u,t)` pairs. -- `ConstantInterpolation(u,t)` - A piecewise constant interpolation. - -- `LinearInterpolation(u,t)` - A linear interpolation. - -- `QuadraticInterpolation(u,t)` - A quadratic interpolation. - -- `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`. - -- `QuadraticSpline(u,t)` - A quadratic spline interpolation. - -- `CubicSpline(u,t)` - A cubic spline interpolation. - -- `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient. - -- `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are: - - `d` - degree of B-spline - - `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method. - - `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector. - -- `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h "methods.md", "Interface" => "interface.md"], -) + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/DataInterpolations/stable/"), + pages = ["index.md", "Methods" => "methods.md", "Interface" => "interface.md"]) deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true) diff --git a/docs/src/index.md b/docs/src/index.md index c5e04a24..7e31d620 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -6,10 +6,9 @@ some methods are mixtures of regressions with interpolations (i.e. do not hit th points exactly, smoothing out the lines). This library can be used to fill in intermediate data points in applications like timeseries data. - ## Installation -To install LinearSolve.jl, use the Julia package manager: +To install DataInterpolations.jl, use the Julia package manager: ```julia using Pkg @@ -21,46 +20,41 @@ Pkg.add("DataInterpolations") In all cases, `u` an `AbstractVector` of values and `t` is an `AbstractVector` of timepoints corresponding to `(u,t)` pairs. -- `ConstantInterpolation(u,t)` - A piecewise constant interpolation. - -- `LinearInterpolation(u,t)` - A linear interpolation. - -- `QuadraticInterpolation(u,t)` - A quadratic interpolation. - -- `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`. - -- `QuadraticSpline(u,t)` - A quadratic spline interpolation. - -- `CubicSpline(u,t)` - A cubic spline interpolation. + - `ConstantInterpolation(u,t)` - A piecewise constant interpolation. -- `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient. - -- `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are: - - `d` - degree of B-spline - - `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method. - - `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector. - -- `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `herrfun(t,u,p), p0, alg) + mfit = optimize(p -> errfun(t, u, p), p0, alg) else if lb === nothing || ub === nothing error("lower or upper bound should not be nothing") end - od = OnceDifferentiable(p->errfun(t,u,p), p0, autodiff=:finite) + od = OnceDifferentiable(p -> errfun(t, u, p), p0, autodiff = :finite) mfit = optimize(od, lb, ub, p0, Fminbox(alg)) end pmin = Optim.minimizer(mfit) @@ -44,11 +64,15 @@ function Curvefit(u, t, model, p0, alg, box=false, lb=nothing, ub=nothing) end # Curvefit -function _interpolate(A::CurvefitCache{<:AbstractVector{<:Number}}, t::Union{AbstractVector{<:Number},Number}) +function _interpolate(A::CurvefitCache{<:AbstractVector{<:Number}}, + t::Union{AbstractVector{<:Number}, Number}) A.m(t, A.pmin) end -_interpolate(A::CurvefitCache{<:AbstractVector{<:Number}}, t::Union{AbstractVector{<:Number},Number}, i) = - _interpolate(A, t), i - +function _interpolate(A::CurvefitCache{<:AbstractVector{<:Number}}, + t::Union{AbstractVector{<:Number}, Number}, + i) + _interpolate(A, t), i +end + end # module diff --git a/ext/DataInterpolationsRegularizationToolsExt.jl b/ext/DataInterpolationsRegularizationToolsExt.jl index 4fa42d74..15fd0aeb 100644 --- a/ext/DataInterpolationsRegularizationToolsExt.jl +++ b/ext/DataInterpolationsRegularizationToolsExt.jl @@ -31,7 +31,6 @@ end # - validate data and t̂ # x unit tests - const LA = LinearAlgebra """ @@ -67,14 +66,14 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, - wls::AbstractVector, wr::AbstractVector, d::Int=2; - λ::Real=1.0, alg::Symbol=:gcv_svd) + wls::AbstractVector, wr::AbstractVector, d::Int = 2; + λ::Real = 1.0, alg::Symbol = :gcv_svd) u, t = munge_data(u, t) M = _mapping_matrix(t̂, t) Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, wls, wr, d, λ, alg, Aitp) end """ Direct smoothing, no `t̂` or weights @@ -82,16 +81,17 @@ Direct smoothing, no `t̂` or weights A = RegularizationSmooth(u, t, d; λ=[1.0], alg=[:gcv_svd]) ``` """ -function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int=2; λ::Real=1.0, - alg::Symbol=:gcv_svd) +function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2; + λ::Real = 1.0, + alg::Symbol = :gcv_svd) u, t = munge_data(u, t) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) Wls½ = Array{Float64}(LA.I, N, N) - Wr½ = Array{Float64}(LA.I, N-d, N-d) + Wr½ = Array{Float64}(LA.I, N - d, N - d) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, LA.diag(Wls½), LA.diag(Wr½), d, λ, alg, Aitp) end """ `t̂` provided, no weights @@ -100,14 +100,14 @@ A = RegularizationSmooth(u, t, t̂, d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, - d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd) + d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd) u, t = munge_data(u, t) N, N̂ = length(t), length(t̂) M = _mapping_matrix(t̂, t) Wls½ = Array{Float64}(LA.I, N, N) - Wr½ = Array{Float64}(LA.I, N̂-d, N̂-d) + Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, LA.diag(Wls½), LA.diag(Wr½), d, λ, alg, Aitp) end """ `t̂` and `wls` provided @@ -116,15 +116,15 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, - wls::AbstractVector, d::Int=2; λ::Real=1.0, - alg::Symbol=:gcv_svd) + wls::AbstractVector, d::Int = 2; λ::Real = 1.0, + alg::Symbol = :gcv_svd) u, t = munge_data(u, t) N, N̂ = length(t), length(t̂) M = _mapping_matrix(t̂, t) Wls½ = LA.diagm(sqrt.(wls)) - Wr½ = Array{Float64}(LA.I, N̂-d, N̂-d) + Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, wls, LA.diag(Wr½), d, λ, alg, Aitp) end """ `wls` provided, no `t̂` @@ -133,16 +133,16 @@ A = RegularizationSmooth(u, t, nothing, wls,d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, - wls::AbstractVector, d::Int=2; λ::Real=1.0, - alg::Symbol=:gcv_svd) + wls::AbstractVector, d::Int = 2; λ::Real = 1.0, + alg::Symbol = :gcv_svd) u, t = munge_data(u, t) t̂ = t N = length(t) M = Array{Float64}(LA.I, N, N) Wls½ = LA.diagm(sqrt.(wls)) - Wr½ = Array{Float64}(LA.I, N-d, N-d) + Wr½ = Array{Float64}(LA.I, N - d, N - d) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, wls, LA.diag(Wr½), d, λ, alg, Aitp) end """ `wls` and `wr` provided, no `t̂` @@ -151,8 +151,8 @@ A = RegularizationSmooth(u, t, nothing, wls, wr, d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, - wls::AbstractVector, wr::AbstractVector, d::Int=2; - λ::Real=1.0, alg::Symbol=:gcv_svd) + wls::AbstractVector, wr::AbstractVector, d::Int = 2; + λ::Real = 1.0, alg::Symbol = :gcv_svd) u, t = munge_data(u, t) t̂ = t N = length(t) @@ -160,7 +160,7 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂,wls, LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, wls, LA.diag(Wr½), d, λ, alg, Aitp) end """ Keyword provided for `wls`, no `t̂` @@ -169,7 +169,7 @@ A = RegularizationSmooth(u, t, nothing, :midpoint, d; λ=[1.0], alg=[:gcv_svd]) ``` """ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing, - wls::Symbol, d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd) + wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd) u, t = munge_data(u, t) t̂ = t N = length(t) @@ -178,20 +178,19 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing Wls½ = LA.diagm(sqrt.(wls)) Wr½ = LA.diagm(sqrt.(wr)) û, λ, Aitp = _reg_smooth_solve(u, t̂, d, M, Wls½, Wr½, λ, alg) - RegularizationSmooth{true}(u,û,t,t̂, LA.diag(Wls½), LA.diag(Wr½), d,λ,alg,Aitp) + RegularizationSmooth{true}(u, û, t, t̂, LA.diag(Wls½), LA.diag(Wr½), d, λ, alg, Aitp) end # """ t̂ provided and keyword for wls _TBD_ """ # function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector, # wls::Symbol, d::Int=2; λ::Real=1.0, alg::Symbol=:gcv_svd) - """ Solve for the smoothed depedent variables and create spline interpolator """ function _reg_smooth_solve(u::AbstractVector, t̂::AbstractVector, d::Int, M::AbstractMatrix, - Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol) + Wls½::AbstractMatrix, Wr½::AbstractMatrix, λ::Real, alg::Symbol) λ = float(λ) # `float` expected by RT - D = _derivative_matrix(t̂,d) - Ψ = RT.setupRegularizationProblem(Wls½*M, Wr½*D) - Wls½u = Wls½*u + D = _derivative_matrix(t̂, d) + Ψ = RT.setupRegularizationProblem(Wls½ * M, Wr½ * D) + Wls½u = Wls½ * u if alg == :fixed b̄ = RT.to_standard_form(Ψ, Wls½u) # via b̄ @@ -200,18 +199,17 @@ function _reg_smooth_solve(u::AbstractVector, t̂::AbstractVector, d::Int, M::Ab else # the provided λ (a scalar) is used as an initial guess; using bounds for Brent() # method is TBD, JJS 12/21/21 - result = RT.solve(Ψ, Wls½u; alg=alg, method=RT.NelderMead(), λ₀=λ) + result = RT.solve(Ψ, Wls½u; alg = alg, method = RT.NelderMead(), λ₀ = λ) û = result.x λ = result.λ end - Aitp = CubicSpline(û,t̂) + Aitp = CubicSpline(û, t̂) # It seems logical to use B-Spline of order d+1, but I am unsure if theory supports the # extra computational cost, JJS 12/25/21 #Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average) return û, λ, Aitp end - """ Order d derivative matrix for the provided t vector """ @@ -220,28 +218,28 @@ function _derivative_matrix(t::AbstractVector, d::Int) if d == 0 return Array{Float64}(LA.I, (N, N)) end - dt = t[d+1:end] - t[1:end-d] + dt = t[(d + 1):end] - t[1:(end - d)] V = LA.diagm(1 ./ dt) - Ddm1 = diff(_derivative_matrix(t, d-1), dims=1) - D = d * V*Ddm1 + Ddm1 = diff(_derivative_matrix(t, d - 1), dims = 1) + D = d * V * Ddm1 return D end """Linear interpolation mapping matrix, which maps `û` to `u`.""" -function _mapping_matrix(t̂::AbstractVector,t::AbstractVector) +function _mapping_matrix(t̂::AbstractVector, t::AbstractVector) N = length(t) N̂ = length(t̂) # map the scattered points to the appropriate index of the smoothed points - idx = searchsortedlast.(Ref(t̂),t) + idx = searchsortedlast.(Ref(t̂), t) # allow for "extrapolation"; i.e. for t̂ extremum that are interior to t - idx[idx.==0] .+= 1 - idx[idx.==N̂] .+= -1 + idx[idx .== 0] .+= 1 + idx[idx .== N̂] .+= -1 # create the linear interpolation matrix - m2 = @. (t - t̂[idx])/(t̂[idx+1] - t̂[idx]) - M = zeros(eltype(t), (N,N̂)) + m2 = @. (t - t̂[idx]) / (t̂[idx + 1] - t̂[idx]) + M = zeros(eltype(t), (N, N̂)) for i in 1:N - M[i,idx[i]] = 1 - m2[i] - M[i,idx[i]+1] = m2[i] + M[i, idx[i]] = 1 - m2[i] + M[i, idx[i] + 1] = m2[i] end return M end @@ -253,15 +251,15 @@ function _weighting_by_kw(t::AbstractVector, d::Int, wls::Symbol) if wls == :midpoint bmp = zeros(N) bmp[1] = -t[1] + t[2] - for i in 2:N-1 - bmp[i] = -t[i-1] + t[i+1] + for i in 2:(N - 1) + bmp[i] = -t[i - 1] + t[i + 1] end - bmp[N] = -t[N-1] + t[N] + bmp[N] = -t[N - 1] + t[N] # divide by 2 doesn't matter in the minimize step, but keeping for correctness if # used as a template elsewhere - bmp = bmp/2 - start = floor(Int, d/2) + 1 - final = iseven(d) ? N - (start - 1) : N - start + bmp = bmp / 2 + start = floor(Int, d / 2) + 1 + final = iseven(d) ? N - (start - 1) : N - start b̃mp = bmp[start:final] return bmp, b̃mp else @@ -269,8 +267,8 @@ function _weighting_by_kw(t::AbstractVector, d::Int, wls::Symbol) end end -_interpolate(A::RegularizationSmooth{<:AbstractVector{<:Number}}, t::Number) = +function _interpolate(A::RegularizationSmooth{<:AbstractVector{<:Number}}, t::Number) _interpolate(A.Aitp, t) +end end # module - diff --git a/ext/DataInterpolationsSymbolicsExt.jl b/ext/DataInterpolationsSymbolicsExt.jl index d45d96c1..2341745b 100644 --- a/ext/DataInterpolationsSymbolicsExt.jl +++ b/ext/DataInterpolationsSymbolicsExt.jl @@ -16,10 +16,12 @@ end SymbolicUtils.promote_symtype(t::AbstractInterpolation, _...) = Real Base.nameof(interp::AbstractInterpolation) = :Interpolation -derivative(interp::AbstractInterpolation, t::Num) = SymbolicUtils.term(derivative, interp, unwrap(t)) +function derivative(interp::AbstractInterpolation, t::Num) + SymbolicUtils.term(derivative, interp, unwrap(t)) +end SymbolicUtils.promote_symtype(::typeof(derivative), _...) = Real -function Symbolics.derivative(interp::AbstractInterpolation, args::NTuple{1,Any}, ::Val{1}) +function Symbolics.derivative(interp::AbstractInterpolation, args::NTuple{1, Any}, ::Val{1}) derivative(interp, Num(args[1])) end diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index c0c6130c..83f6c4d4 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -2,16 +2,18 @@ module DataInterpolations ### Interface Functionality -abstract type AbstractInterpolation{FT,T} <: AbstractVector{T} end +abstract type AbstractInterpolation{FT, T} <: AbstractVector{T} end Base.size(A::AbstractInterpolation) = size(A.u) Base.size(A::AbstractInterpolation{true}) = length(A.u) .+ size(A.t) -Base.getindex(A::AbstractInterpolation,i) = A.u[i] -Base.getindex(A::AbstractInterpolation{true},i) = - i<=length(A.u) ? A.u[i] : A.t[i-length(A.u)] -Base.setindex!(A::AbstractInterpolation,x,i) = A.u[i] = x -Base.setindex!(A::AbstractInterpolation{true},x,i) = - i <= length(A.u) ? (A.u[i] = x) : (A.t[i-length(A.u)] = x) +Base.getindex(A::AbstractInterpolation, i) = A.u[i] +function Base.getindex(A::AbstractInterpolation{true}, i) + i <= length(A.u) ? A.u[i] : A.t[i - length(A.u)] +end +Base.setindex!(A::AbstractInterpolation, x, i) = A.u[i] = x +function Base.setindex!(A::AbstractInterpolation{true}, x, i) + i <= length(A.u) ? (A.u[i] = x) : (A.t[i - length(A.u)] = x) +end using LinearAlgebra, RecursiveArrayTools, RecipesBase @@ -27,11 +29,11 @@ include("online.jl") (interp::AbstractInterpolation)(t::Number, i::Integer) = _interpolate(interp, t, i) (interp::AbstractInterpolation)(t::AbstractVector) = interp(similar(t, eltype(interp)), t) function (interp::AbstractInterpolation)(u::AbstractVector, t::AbstractVector) - iguess = firstindex(interp.t) - @inbounds for i in eachindex(u, t) - u[i], iguess = interp(t[i], iguess) - end - u + iguess = firstindex(interp.t) + @inbounds for i in eachindex(u, t) + u[i], iguess = interp(t[i], iguess) + end + u end export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, @@ -40,7 +42,7 @@ export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, # added for RegularizationSmooth, JJS 11/27/21 ### Regularization data smoothing and interpolation -struct RegularizationSmooth{uType, tType, FT, T, T2} <: AbstractInterpolation{FT,T} +struct RegularizationSmooth{uType, tType, FT, T, T2} <: AbstractInterpolation{FT, T} u::uType û::uType t::tType @@ -50,9 +52,19 @@ struct RegularizationSmooth{uType, tType, FT, T, T2} <: AbstractInterpolation{FT d::Int # derivative degree used to calculate the roughness λ::T2 # regularization parameter alg::Symbol # how to determine λ: `:fixed`, `:gcv_svd`, `:gcv_tr`, `L_curve` - Aitp::AbstractInterpolation{FT,T} - RegularizationSmooth{FT}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp) where FT = - new{typeof(u), typeof(t), FT, eltype(u), typeof(λ)}(u,û,t,t̂,wls,wr,d,λ,alg,Aitp) + Aitp::AbstractInterpolation{FT, T} + function RegularizationSmooth{FT}(u, û, t, t̂, wls, wr, d, λ, alg, Aitp) where {FT} + new{typeof(u), typeof(t), FT, eltype(u), typeof(λ)}(u, + û, + t, + t̂, + wls, + wr, + d, + λ, + alg, + Aitp) + end end export RegularizationSmooth @@ -63,10 +75,18 @@ end @static if !isdefined(Base, :get_extension) function __init__() - Requires.@require ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" begin include("../ext/DataInterpolationsChainRulesCoreExt.jl") end - Requires.@require Optim = "429524aa-4258-5aef-a3af-852621145aeb" begin include("../ext/DataInterpolationsOptimExt.jl") end - Requires.@require RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182" begin include("../ext/DataInterpolationsRegularizationToolsExt.jl") end - Requires.@require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" begin include("../ext/DataInterpolationsSymbolicsExt.jl") end + Requires.@require ChainRulesCore="d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" begin + include("../ext/DataInterpolationsChainRulesCoreExt.jl") + end + Requires.@require Optim="429524aa-4258-5aef-a3af-852621145aeb" begin + include("../ext/DataInterpolationsOptimExt.jl") + end + Requires.@require RegularizationTools="29dad682-9a27-4bc3-9c72-016788665182" begin + include("../ext/DataInterpolationsRegularizationToolsExt.jl") + end + Requires.@require Symbolics="0c5d862f-8b57-4792-8d23-62f2024744c7" begin + include("../ext/DataInterpolationsSymbolicsExt.jl") + end end end diff --git a/src/derivatives.jl b/src/derivatives.jl index a856fbf4..e428149f 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,183 +1,190 @@ -derivative(A, t) = derivative(A, t, firstindex(A.t)-1)[1] +derivative(A, t) = derivative(A, t, firstindex(A.t) - 1)[1] function derivative(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess) - idx = searchsortedfirstcorrelated(A.t, t, iguess) - if A.t[idx] >= t - idx -= 1 - end - idx == 0 ? idx += 1 : nothing - θ = 1 / (A.t[idx+1] - A.t[idx]) - (A.u[idx+1] - A.u[idx]) / (A.t[idx+1] - A.t[idx]), idx + idx = searchsortedfirstcorrelated(A.t, t, iguess) + if A.t[idx] >= t + idx -= 1 + end + idx == 0 ? idx += 1 : nothing + θ = 1 / (A.t[idx + 1] - A.t[idx]) + (A.u[idx + 1] - A.u[idx]) / (A.t[idx + 1] - A.t[idx]), idx end function derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) - idx = searchsortedfirstcorrelated(A.t, t, iguess) - if A.t[idx] >= t - idx -= 1 - end - idx == 0 ? idx += 1 : nothing - θ = 1 / (A.t[idx+1] - A.t[idx]) - (@views @. (A.u[:, idx+1] - A.u[:, idx]) / (A.t[idx+1] - A.t[idx])), idx + idx = searchsortedfirstcorrelated(A.t, t, iguess) + if A.t[idx] >= t + idx -= 1 + end + idx == 0 ? idx += 1 : nothing + θ = 1 / (A.t[idx + 1] - A.t[idx]) + (@views @. (A.u[:, idx + 1] - A.u[:, idx]) / (A.t[idx + 1] - A.t[idx])), idx end function derivative(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess) - i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂, i₀ + i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) + dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) + dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) + dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) + A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂, i₀ end function derivative(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess) - idx = searchsortedfirstcorrelated(A.t, t, iguess) - if A.t[idx] >= t - idx -= 1 - end - idx == 0 ? idx += 1 : nothing - if idx == length(A.t) - 1 - i₀ = idx - 1; i₁ = idx; i₂ = i₁ + 1; - else - i₀ = idx; i₁ = i₀ + 1; i₂ = i₁ + 1; - end - dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) - dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) - dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) - (@views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂), idx + idx = searchsortedfirstcorrelated(A.t, t, iguess) + if A.t[idx] >= t + idx -= 1 + end + idx == 0 ? idx += 1 : nothing + if idx == length(A.t) - 1 + i₀ = idx - 1 + i₁ = idx + i₂ = i₁ + 1 + else + i₀ = idx + i₁ = i₀ + 1 + i₂ = i₁ + 1 + end + dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) + dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) + dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) + (@views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂), idx end function derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) - idxs = findRequiredIdxs(A, t) - if A.t[idxs[1]] == t - return zero(A.u[idxs[1]]) - end - G = zero(A.u[1]); F = zero(A.t[1]) - DG = zero(A.u[1]); DF = zero(A.t[1]) - tmp = G - for i = 1:length(idxs) - if isnan(A.bcache[idxs[i]]) - mult = one(A.t[1]) - for j = 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - for j = (i+1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - A.bcache[idxs[i]] = mult - else - mult = A.bcache[idxs[i]] + idxs = findRequiredIdxs(A, t) + if A.t[idxs[1]] == t + return zero(A.u[idxs[1]]) end - wi = inv(mult) - tti = t - A.t[idxs[i]] - tmp = wi / (t - A.t[idxs[i]]) - g = tmp * A.u[idxs[i]] - G += g - DG -= g / (t - A.t[idxs[i]]) - F += tmp - DF -= tmp / (t - A.t[idxs[i]]) - end - (DG * F - G * DF) / (F ^ 2) + G = zero(A.u[1]) + F = zero(A.t[1]) + DG = zero(A.u[1]) + DF = zero(A.t[1]) + tmp = G + for i in 1:length(idxs) + if isnan(A.bcache[idxs[i]]) + mult = one(A.t[1]) + for j in 1:(i - 1) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + for j in (i + 1):length(idxs) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + A.bcache[idxs[i]] = mult + else + mult = A.bcache[idxs[i]] + end + wi = inv(mult) + tti = t - A.t[idxs[i]] + tmp = wi / (t - A.t[idxs[i]]) + g = tmp * A.u[idxs[i]] + G += g + DG -= g / (t - A.t[idxs[i]]) + F += tmp + DF -= tmp / (t - A.t[idxs[i]]) + end + (DG * F - G * DF) / (F^2) end function derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) - idxs = findRequiredIdxs(A, t) - if A.t[idxs[1]] == t - return zero(A.u[:, idxs[1]]) - end - G = zero(A.u[:, 1]); F = zero(A.t[1]) - DG = zero(A.u[:, 1]); DF = zero(A.t[1]) - tmp = G - for i = 1:length(idxs) - if isnan(A.bcache[idxs[i]]) - mult = one(A.t[1]) - for j = 1:(i - 1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - for j = (i+1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - A.bcache[idxs[i]] = mult - else - mult = A.bcache[idxs[i]] + idxs = findRequiredIdxs(A, t) + if A.t[idxs[1]] == t + return zero(A.u[:, idxs[1]]) + end + G = zero(A.u[:, 1]) + F = zero(A.t[1]) + DG = zero(A.u[:, 1]) + DF = zero(A.t[1]) + tmp = G + for i in 1:length(idxs) + if isnan(A.bcache[idxs[i]]) + mult = one(A.t[1]) + for j in 1:(i - 1) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + for j in (i + 1):length(idxs) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + A.bcache[idxs[i]] = mult + else + mult = A.bcache[idxs[i]] + end + wi = inv(mult) + tti = t - A.t[idxs[i]] + tmp = wi / (t - A.t[idxs[i]]) + g = tmp * A.u[:, idxs[i]] + @. G += g + @. DG -= g / (t - A.t[idxs[i]]) + F += tmp + DF -= tmp / (t - A.t[idxs[i]]) end - wi = inv(mult) - tti = t - A.t[idxs[i]] - tmp = wi / (t - A.t[idxs[i]]) - g = tmp * A.u[:, idxs[i]] - @. G += g - @. DG -= g / (t - A.t[idxs[i]]) - F += tmp - DF -= tmp / (t - A.t[idxs[i]]) - end - @. (DG * F - G * DF) / (F ^ 2) + @. (DG * F - G * DF) / (F^2) end -derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = - derivative(A, t), i -derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) = - derivative(A, t), i +derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = derivative(A, t), i +derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) = derivative(A, t), i function derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) - t < A.t[1] && return zero(A.u[1]), 1 - t > A.t[end] && return zero(A.u[end]), lastindex(t) - i = searchsortedlastcorrelated(A.t, t, iguess) - j = min(i, length(A.c)) # for smooth derivative at A.t[end] - wj = t - A.t[i] - (@evalpoly wj A.b[i] 2A.c[j] 3A.d[j]), i + t < A.t[1] && return zero(A.u[1]), 1 + t > A.t[end] && return zero(A.u[end]), lastindex(t) + i = searchsortedlastcorrelated(A.t, t, iguess) + j = min(i, length(A.c)) # for smooth derivative at A.t[end] + wj = t - A.t[i] + (@evalpoly wj A.b[i] 2A.c[j] 3A.d[j]), i end function derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number) - return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) + return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) end function derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number) - return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1] + return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1] end # QuadraticSpline Interpolation function derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - i = searchsortedfirstcorrelated(A.t, t, iguess) - i == 1 ? i += 1 : nothing - σ = 1//2 * (A.z[i] - A.z[i - 1]) / (A.t[i] - A.t[i - 1]) - A.z[i-1] + 2σ * (t - A.t[i-1]), i + i = searchsortedfirstcorrelated(A.t, t, iguess) + i == 1 ? i += 1 : nothing + σ = 1 // 2 * (A.z[i] - A.z[i - 1]) / (A.t[i] - A.t[i - 1]) + A.z[i - 1] + 2σ * (t - A.t[i - 1]), i end # CubicSpline Interpolation function derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess) - i = searchsortedfirstcorrelated(A.t, t, iguess) - isnothing(i) ? i = length(A.t) - 1 : i -= 1 - i == 0 ? i += 1 : nothing - dI = -3A.z[i] * (A.t[i + 1] - t)^2 / (6A.h[i + 1]) + 3A.z[i + 1] * (t - A.t[i])^2 / (6A.h[i + 1]) - dC = A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6 - dD = -(A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6) - dI + dC + dD, i + i = searchsortedfirstcorrelated(A.t, t, iguess) + isnothing(i) ? i = length(A.t) - 1 : i -= 1 + i == 0 ? i += 1 : nothing + dI = -3A.z[i] * (A.t[i + 1] - t)^2 / (6A.h[i + 1]) + + 3A.z[i + 1] * (t - A.t[i])^2 / (6A.h[i + 1]) + dC = A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6 + dD = -(A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6) + dI + dC + dD, i end function derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) - # change t into param [0 1] - idx = searchsortedlastcorrelated(A.t,t, iguess) - idx == length(A.t) ? idx -= 1 : nothing - 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 = DataInterpolations.spline_coefficients(n, A.d-1, A.k, t_) - ducum = zero(eltype(A.u)) - for i = 1:(n - 1) - ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) - end - ducum * A.d * scale, idx + # change t into param [0 1] + idx = searchsortedlastcorrelated(A.t, t, iguess) + idx == length(A.t) ? idx -= 1 : nothing + 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 = DataInterpolations.spline_coefficients(n, A.d - 1, A.k, t_) + ducum = zero(eltype(A.u)) + 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]) + end + ducum * A.d * scale, idx end # BSpline Curve Approx function derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess) - # change t into param [0 1] - idx = searchsortedlastcorrelated(A.t,t, iguess) - idx == 0 ? idx += 1 : nothing - 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 = spline_coefficients(A.h, A.d-1, A.k, t_) - ducum = zero(eltype(A.u)) - for i = 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]) - end - ducum * A.d * scale, idx + # change t into param [0 1] + idx = searchsortedlastcorrelated(A.t, t, iguess) + idx == 0 ? idx += 1 : nothing + 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 = spline_coefficients(A.h, A.d - 1, A.k, t_) + ducum = zero(eltype(A.u)) + 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]) + end + ducum * A.d * scale, idx end diff --git a/src/integrals.jl b/src/integrals.jl index d15dcc28..1c42a241 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -1,85 +1,98 @@ function integral(A::AbstractInterpolation, t::Number) bw, fw = samples(A) - idx = max(1+bw, min(searchsortedlast(A.t, t), length(A.t) - fw)) + idx = max(1 + bw, min(searchsortedlast(A.t, t), length(A.t) - fw)) _integral(A, idx, t) end function integral(A::AbstractInterpolation, t1::Number, t2::Number) bw, fw = samples(A) # the index less than or equal to t1 - idx1 = max(1+bw, min(searchsortedlast(A.t, t1), length(A.t) - fw)) + idx1 = max(1 + bw, min(searchsortedlast(A.t, t1), length(A.t) - fw)) # the index less than t2 - idx2 = max(2+bw, min(searchsortedlast(A.t, t2), length(A.t) - fw)) + idx2 = max(2 + bw, min(searchsortedlast(A.t, t2), length(A.t) - fw)) if A.t[idx2] == t2 - idx2-=1 + idx2 -= 1 end total = zero(eltype(A)) for idx in idx1:idx2 lt1 = idx == idx1 ? t1 : A.t[idx] - lt2 = idx == idx2 ? t2 : A.t[idx+1] - total += _integral(A, idx, lt2)-_integral(A, idx, lt1) + lt2 = idx == idx2 ? t2 : A.t[idx + 1] + total += _integral(A, idx, lt2) - _integral(A, idx, lt1) end total end samples(A::LinearInterpolation{<:AbstractVector}) = (0, 1) -function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx+1] - u1 = A.u[idx] - u2 = A.u[idx+1] - t^2*(u1 - u2)/(2*t1 - 2*t2) + t*(t1*u2 - t2*u1)/(t1 - t2) +function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, + idx::Number, + t::Number) + t1 = A.t[idx] + t2 = A.t[idx + 1] + u1 = A.u[idx] + u2 = A.u[idx + 1] + t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2) end samples(A::ConstantInterpolation{<:AbstractVector}) = (0, 1) function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number) - if A.dir === :left - # :left means that value to the left is used for interpolation - return A.u[idx]*t - else - # :right means that value to the right is used for interpolation - return A.u[idx+1]*t - end + if A.dir === :left + # :left means that value to the left is used for interpolation + return A.u[idx] * t + else + # :right means that value to the right is used for interpolation + return A.u[idx + 1] * t + end end samples(A::QuadraticInterpolation{<:AbstractVector}) = (0, 1) -function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) - A.mode == :Backward && idx > 1 && (idx -= 1) - idx = min(length(A.t) - 2, idx) - t1 = A.t[idx] - t2 = A.t[idx+1] - t3 = A.t[idx+2] - u1 = A.u[idx] - u2 = A.u[idx+1] - u3 = A.u[idx+2] - (t^3*(-t1*u2 + t1*u3 + t2*u1 - t2*u3 - t3*u1 + t3*u2)/ - (3*t1^2*t2 - 3*t1^2*t3 - 3*t1*t2^2 + 3*t1*t3^2 + 3*t2^2*t3 - 3*t2*t3^2) + - t^2*(t1^2*u2 - t1^2*u3 - t2^2*u1 + t2^2*u3 + t3^2*u1 - t3^2*u2)/ - (2*t1^2*t2 - 2*t1^2*t3 - 2*t1*t2^2 + 2*t1*t3^2 + 2*t2^2*t3 - 2*t2*t3^2) + - t*(t1^2*t2*u3 - t1^2*t3*u2 - t1*t2^2*u3 + t1*t3^2*u2 + t2^2*t3*u1 - t2*t3^2*u1)/ - (t1^2*t2 - t1^2*t3 - t1*t2^2 + t1*t3^2 + t2^2*t3 - t2*t3^2)) +function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, + idx::Number, + t::Number) + A.mode == :Backward && idx > 1 && (idx -= 1) + idx = min(length(A.t) - 2, idx) + t1 = A.t[idx] + t2 = A.t[idx + 1] + t3 = A.t[idx + 2] + u1 = A.u[idx] + u2 = A.u[idx + 1] + u3 = A.u[idx + 2] + (t^3 * (-t1 * u2 + t1 * u3 + t2 * u1 - t2 * u3 - t3 * u1 + t3 * u2) / + (3 * t1^2 * t2 - 3 * t1^2 * t3 - 3 * t1 * t2^2 + 3 * t1 * t3^2 + 3 * t2^2 * t3 - + 3 * t2 * t3^2) + + t^2 * (t1^2 * u2 - t1^2 * u3 - t2^2 * u1 + t2^2 * u3 + t3^2 * u1 - t3^2 * u2) / + (2 * t1^2 * t2 - 2 * t1^2 * t3 - 2 * t1 * t2^2 + 2 * t1 * t3^2 + 2 * t2^2 * t3 - + 2 * t2 * t3^2) + + t * + (t1^2 * t2 * u3 - t1^2 * t3 * u2 - t1 * t2^2 * u3 + t1 * t3^2 * u2 + t2^2 * t3 * u1 - + t2 * t3^2 * u1) / + (t1^2 * t2 - t1^2 * t3 - t1 * t2^2 + t1 * t3^2 + t2^2 * t3 - t2 * t3^2)) end samples(A::QuadraticSpline{<:AbstractVector{<:Number}}) = (0, 1) function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx+1] - u1 = A.u[idx] - z1 = A.z[idx] - z2 = A.z[idx+1] - t^3*(z1 - z2)/(6*t1 - 6*t2) + t^2*(t1*z2 - t2*z1)/(2*t1 - 2*t2) + t*(-t1^2*z1 - t1^2*z2 + 2*t1*t2*z1 + 2*t1*u1 - 2*t2*u1)/(2*t1 - 2*t2) + t1 = A.t[idx] + t2 = A.t[idx + 1] + u1 = A.u[idx] + z1 = A.z[idx] + z2 = A.z[idx + 1] + t^3 * (z1 - z2) / (6 * t1 - 6 * t2) + t^2 * (t1 * z2 - t2 * z1) / (2 * t1 - 2 * t2) + + t * (-t1^2 * z1 - t1^2 * z2 + 2 * t1 * t2 * z1 + 2 * t1 * u1 - 2 * t2 * u1) / + (2 * t1 - 2 * t2) end samples(A::CubicSpline{<:AbstractVector{<:Number}}) = (0, 1) function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx+1] - u1 = A.u[idx] - u2 = A.u[idx+1] - z1 = A.z[idx] - z2 = A.z[idx+1] - h2 = A.h[idx+1] - (t^4*(-z1 + z2)/(24*h2) + t^3*(-t1*z2 + t2*z1)/(6*h2) + - t^2*(h2^2*z1 - h2^2*z2 + 3*t1^2*z2 - 3*t2^2*z1 - 6*u1 + 6*u2)/(12*h2) + - t*(h2^2*t1*z2 - h2^2*t2*z1 - t1^3*z2 - 6*t1*u2 + t2^3*z1 + 6*t2*u1)/(6*h2)) + t1 = A.t[idx] + t2 = A.t[idx + 1] + u1 = A.u[idx] + u2 = A.u[idx + 1] + z1 = A.z[idx] + z2 = A.z[idx + 1] + h2 = A.h[idx + 1] + (t^4 * (-z1 + z2) / (24 * h2) + t^3 * (-t1 * z2 + t2 * z1) / (6 * h2) + + t^2 * (h2^2 * z1 - h2^2 * z2 + 3 * t1^2 * z2 - 3 * t2^2 * z1 - 6 * u1 + 6 * u2) / + (12 * h2) + + t * + (h2^2 * t1 * z2 - h2^2 * t2 * z1 - t1^3 * z2 - 6 * t1 * u2 + t2^3 * z1 + 6 * t2 * u1) / + (6 * h2)) end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index fe21ef4e..ea39ab7f 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -1,364 +1,413 @@ ### Linear Interpolation -struct LinearInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - LinearInterpolation{FT}(u,t) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t) +struct LinearInterpolation{uType, tType, FT, T} <: AbstractInterpolation{FT, T} + u::uType + t::tType + function LinearInterpolation{FT}(u, t) where {FT} + new{typeof(u), typeof(t), FT, eltype(u)}(u, t) + end end -function LinearInterpolation(u,t) - u, t = munge_data(u, t) - LinearInterpolation{true}(u,t) +function LinearInterpolation(u, t) + u, t = munge_data(u, t) + LinearInterpolation{true}(u, t) end ### Quadratic Interpolation -struct QuadraticInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - mode::Symbol - function QuadraticInterpolation{FT}(u,t,mode) where FT - mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") - new{typeof(u),typeof(t),FT,eltype(u)}(u,t,mode) - end +struct QuadraticInterpolation{uType, tType, FT, T} <: AbstractInterpolation{FT, T} + u::uType + t::tType + mode::Symbol + function QuadraticInterpolation{FT}(u, t, mode) where {FT} + mode ∈ (:Forward, :Backward) || + error("mode should be :Forward or :Backward for QuadraticInterpolation") + new{typeof(u), typeof(t), FT, eltype(u)}(u, t, mode) + end end -function QuadraticInterpolation(u,t,mode) - u, t = munge_data(u, t) - QuadraticInterpolation{true}(u,t,mode) +function QuadraticInterpolation(u, t, mode) + u, t = munge_data(u, t) + QuadraticInterpolation{true}(u, t, mode) end -QuadraticInterpolation(u,t) = QuadraticInterpolation(u,t,:Forward) +QuadraticInterpolation(u, t) = QuadraticInterpolation(u, t, :Forward) ### Lagrange Interpolation -struct LagrangeInterpolation{uType,tType,FT,T,bcacheType} <: AbstractInterpolation{FT,T} - u::uType - t::tType - n::Int - bcache::bcacheType - function LagrangeInterpolation{FT}(u,t,n) where FT - bcache = zeros(eltype(u[1]),n+1) - fill!(bcache, NaN) - new{typeof(u),typeof(t),FT,eltype(u),typeof(bcache)}(u,t,n,bcache) - end +struct LagrangeInterpolation{uType, tType, FT, T, bcacheType} <: + AbstractInterpolation{FT, T} + u::uType + t::tType + n::Int + bcache::bcacheType + function LagrangeInterpolation{FT}(u, t, n) where {FT} + bcache = zeros(eltype(u[1]), n + 1) + fill!(bcache, NaN) + new{typeof(u), typeof(t), FT, eltype(u), typeof(bcache)}(u, t, n, bcache) + end end -function LagrangeInterpolation(u,t,n=nothing) - u, t = munge_data(u, t) - if isnothing(n) - n = length(t) - 1 # degree - end - if n != length(t) - 1 - error("Currently only n=length(t) - 1 is supported") - end - LagrangeInterpolation{true}(u,t,n) +function LagrangeInterpolation(u, t, n = nothing) + u, t = munge_data(u, t) + if isnothing(n) + n = length(t) - 1 # degree + end + if n != length(t) - 1 + error("Currently only n=length(t) - 1 is supported") + end + LagrangeInterpolation{true}(u, t, n) end ### Akima Interpolation -struct AkimaInterpolation{uType,tType,bType,cType,dType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - b::bType - c::cType - d::dType - AkimaInterpolation{FT}(u,t,b,c,d) where FT = new{typeof(u),typeof(t),typeof(b),typeof(c), - typeof(d),FT,eltype(u)}(u,t,b,c,d) +struct AkimaInterpolation{uType, tType, bType, cType, dType, FT, T} <: + AbstractInterpolation{FT, T} + u::uType + t::tType + b::bType + c::cType + d::dType + function AkimaInterpolation{FT}(u, t, b, c, d) where {FT} + new{typeof(u), typeof(t), typeof(b), typeof(c), + typeof(d), FT, eltype(u)}(u, + t, + b, + c, + d) + end end function AkimaInterpolation(u, t) - u, t = munge_data(u, t) - n = length(t) - dt = diff(t) - m = Array{eltype(u)}(undef, n + 3) - m[3:end-2] = diff(u) ./ dt - m[2] = 2m[3] - m[4] - m[1] = 2m[2] - m[3] - m[end-1] = 2m[end-2] - m[end-3] - m[end] = 2m[end-1] - m[end-2] - - b = 0.5 .* (m[4:end] .+ m[1:end-3]) - dm = abs.(diff(m)) - f1 = dm[3:n+2] - f2 = dm[1:n] - f12 = f1 + f2 - ind = findall(f12 .> 1e-9 * maximum(f12)) - b[ind] = (f1[ind] .* m[ind .+ 1] .+ - f2[ind] .* m[ind .+ 2]) ./ f12[ind] - c = (3.0 .* m[3:end-2] .- 2.0 .* b[1:end-1] .- b[2:end]) ./ dt - d = (b[1:end-1] .+ b[2:end] .- 2.0 .* m[3:end-2]) ./ dt.^2 - - AkimaInterpolation{true}(u,t,b,c,d) + u, t = munge_data(u, t) + n = length(t) + dt = diff(t) + m = Array{eltype(u)}(undef, n + 3) + m[3:(end - 2)] = diff(u) ./ dt + m[2] = 2m[3] - m[4] + m[1] = 2m[2] - m[3] + m[end - 1] = 2m[end - 2] - m[end - 3] + m[end] = 2m[end - 1] - m[end - 2] + + b = 0.5 .* (m[4:end] .+ m[1:(end - 3)]) + dm = abs.(diff(m)) + f1 = dm[3:(n + 2)] + f2 = dm[1:n] + f12 = f1 + f2 + ind = findall(f12 .> 1e-9 * maximum(f12)) + b[ind] = (f1[ind] .* m[ind .+ 1] .+ + f2[ind] .* m[ind .+ 2]) ./ f12[ind] + c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt + d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 + + AkimaInterpolation{true}(u, t, b, c, d) end - ### ConstantInterpolation Interpolation -struct ConstantInterpolation{uType,tType,dirType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - dir::Symbol # indicates if value to the $dir should be used for the interpolation - ConstantInterpolation{FT}(u,t,dir) where FT = new{typeof(u),typeof(t),typeof(dir),FT,eltype(u)}(u,t,dir) +struct ConstantInterpolation{uType, tType, dirType, FT, T} <: AbstractInterpolation{FT, T} + u::uType + t::tType + dir::Symbol # indicates if value to the $dir should be used for the interpolation + function ConstantInterpolation{FT}(u, t, dir) where {FT} + new{typeof(u), typeof(t), typeof(dir), FT, eltype(u)}(u, t, dir) + end end -function ConstantInterpolation(u,t;dir=:left) - u, t = munge_data(u, t) - ConstantInterpolation{true}(u,t,dir) +function ConstantInterpolation(u, t; dir = :left) + u, t = munge_data(u, t) + ConstantInterpolation{true}(u, t, dir) end Base.@deprecate_binding ZeroSpline ConstantInterpolation ### QuadraticSpline Interpolation -struct QuadraticSpline{uType,tType,tAType,dType,zType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - tA::tAType - d::dType - z::zType - QuadraticSpline{FT}(u,t,tA,d,z) where FT = new{typeof(u),typeof(t),typeof(tA), - typeof(d),typeof(z),FT,eltype(u)}(u,t,tA,d,z) +struct QuadraticSpline{uType, tType, tAType, dType, zType, FT, T} <: + AbstractInterpolation{FT, T} + u::uType + t::tType + tA::tAType + d::dType + z::zType + function QuadraticSpline{FT}(u, t, tA, d, z) where {FT} + new{typeof(u), typeof(t), typeof(tA), + typeof(d), typeof(z), FT, eltype(u)}(u, + t, + tA, + d, + z) + end end -function QuadraticSpline(u::uType,t) where {uType<:AbstractVector{<:Number}} - u, t = munge_data(u, 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) +function QuadraticSpline(u::uType, t) where {uType <: AbstractVector{<:Number}} + u, t = munge_data(u, 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])) + # 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 - QuadraticSpline{true}(u,t,tA,d,z) + d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) + z = tA \ d + QuadraticSpline{true}(u, t, tA, d, z) end -function QuadraticSpline(u::uType,t) where {uType<:AbstractVector} - u, t = munge_data(u, 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_))] - QuadraticSpline{true}(u,t,tA,d,z) +function QuadraticSpline(u::uType, t) where {uType <: AbstractVector} + u, t = munge_data(u, 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_))] + QuadraticSpline{true}(u, t, tA, d, z) end # Cubic Spline Interpolation -struct CubicSpline{uType,tType,hType,zType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - h::hType - z::zType - CubicSpline{FT}(u,t,h,z) where FT = new{typeof(u),typeof(t),typeof(h),typeof(z),FT,eltype(u)}(u,t,h,z) +struct CubicSpline{uType, tType, hType, zType, FT, T} <: AbstractInterpolation{FT, T} + u::uType + t::tType + h::hType + z::zType + function CubicSpline{FT}(u, t, h, z) where {FT} + new{typeof(u), typeof(t), typeof(h), typeof(z), FT, eltype(u)}(u, t, h, z) + end end -function CubicSpline(u::uType,t) where {uType<:AbstractVector{<:Number}} - 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 = h[2:n+1] - d_tmp = 2 .* (h[1:n+1] .+ h[2:n+2]) - du = h[2:n+1] - tA = Tridiagonal(dl,d_tmp,du) - - # zero for element type of d, which we don't know yet - typed_zero = zero(6(u[begin+2] - u[begin+1]) / h[begin+2] - 6(u[begin+1] - u[begin]) / h[begin+1]) - - d = map(i -> i == 1 || i == n + 1 ? typed_zero : 6(u[i+1] - u[i]) / h[i+1] - 6(u[i] - u[i-1]) / h[i], 1:n+1) - z = tA\d - CubicSpline{true}(u,t,h[1:n+1],z) +function CubicSpline(u::uType, t) where {uType <: AbstractVector{<:Number}} + 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 = h[2:(n + 1)] + d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) + du = h[2:(n + 1)] + tA = Tridiagonal(dl, d_tmp, du) + + # zero for element type of d, which we don't know yet + typed_zero = zero(6(u[begin + 2] - u[begin + 1]) / h[begin + 2] - + 6(u[begin + 1] - u[begin]) / h[begin + 1]) + + d = map(i -> i == 1 || i == n + 1 ? typed_zero : + 6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i], + 1:(n + 1)) + z = tA \ d + CubicSpline{true}(u, t, h[1:(n + 1)], z) end -function CubicSpline(u::uType,t) where {uType<:AbstractVector} - 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 = h[2:n+1] - d_tmp = 2 .* (h[1:n+1] .+ h[2:n+2]) - du = h[2:n+1] - tA = Tridiagonal(dl,d_tmp,du) - d_ = map(i -> i == 1 || i == n + 1 ? zeros(eltype(t),size(u[1])) : 6(u[i+1] - u[i]) / h[i+1] - 6(u[i] - u[i-1]) / h[i], 1:n+1) - d = transpose(reshape(reduce(hcat, d_), :, n+1)) - z_ = reshape(transpose(tA\d), size(u[1])...,:) - z = [z_s for z_s in eachslice(z_, dims=ndims(z_))] - CubicSpline{true}(u,t,h[1:n+1],z) +function CubicSpline(u::uType, t) where {uType <: AbstractVector} + 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 = h[2:(n + 1)] + d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) + du = h[2:(n + 1)] + tA = Tridiagonal(dl, d_tmp, du) + d_ = map(i -> i == 1 || i == n + 1 ? zeros(eltype(t), size(u[1])) : + 6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i], + 1:(n + 1)) + d = transpose(reshape(reduce(hcat, d_), :, n + 1)) + z_ = reshape(transpose(tA \ d), size(u[1])..., :) + z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] + CubicSpline{true}(u, t, h[1:(n + 1)], z) end ### BSpline Curve Interpolation -struct BSplineInterpolation{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - d::Int # degree - p::pType # params vector - k::kType # knot vector - c::cType # control points - pVecType::Symbol - knotVecType::Symbol - BSplineInterpolation{FT}(u,t,d,p,k,c,pVecType,knotVecType) where FT = new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,p,k,c,pVecType,knotVecType) +struct BSplineInterpolation{uType, tType, pType, kType, cType, FT, T} <: + AbstractInterpolation{FT, T} + u::uType + t::tType + d::Int # degree + p::pType # params vector + k::kType # knot vector + c::cType # control points + pVecType::Symbol + knotVecType::Symbol + function BSplineInterpolation{FT}(u, t, d, p, k, c, pVecType, knotVecType) where {FT} + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), FT, eltype(u)}(u, + t, + d, + p, + k, + c, + pVecType, + knotVecType) + end end -function BSplineInterpolation(u,t,d,pVecType,knotVecType) - u, t = munge_data(u, t) - n = length(t) - 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)) - - for i = 2:n - s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2) - l[i-1] = s - end - if pVecType == :Uniform - for i = 2:(n-1) - p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1) +function BSplineInterpolation(u, t, d, pVecType, knotVecType) + u, t = munge_data(u, t) + n = length(t) + 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)) + + for i in 2:n + s += √((t[i] - t[i - 1])^2 + (u[i] - u[i - 1])^2) + l[i - 1] = s end - elseif pVecType == :ArcLen - for i = 2:(n-1) - p[i] = p[1] + l[i-1]/s * (p[end]-p[1]) + 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 - 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 = 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 = (d+2):n - k[i] = k[1] + (i-d-1)//(n-d) * (k[end]-k[1]) + + 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 - elseif knotVecType == :Average - # average spaced knot vector - idx = 1 - if d+2 <= n - k[d+2] = 1//d * ps[d] + + ps = zeros(eltype(t), n - 2) + s = zero(eltype(t)) + for i in 2:(n - 1) + s += p[i] + ps[i - 1] = s end - for i = (d+3):n - k[i] = 1//d * (ps[idx+d] - ps[idx]) - idx += 1 + + 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 - end - # control points - N = spline_coefficients(n,d,k,p) - c = vec(N\u[:,:]) - BSplineInterpolation{true}(u,t,d,p,k,c,pVecType,knotVecType) + # control points + N = spline_coefficients(n, d, k, p) + c = vec(N \ u[:, :]) + BSplineInterpolation{true}(u, t, d, p, k, c, pVecType, knotVecType) end ### BSpline Curve Approx -struct BSplineApprox{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T} - u::uType - t::tType - d::Int # degree - h::Int # number of control points (n => h >= d >= 1) - p::pType # params vector - k::kType # knot vector - c::cType # control points - pVecType::Symbol - knotVecType::Symbol - BSplineApprox{FT}(u,t,d,h,p,k,c,pVecType,knotVecType) where FT = new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,h,p,k,c,pVecType,knotVecType) +struct BSplineApprox{uType, tType, pType, kType, cType, FT, T} <: + AbstractInterpolation{FT, T} + u::uType + t::tType + d::Int # degree + h::Int # number of control points (n => h >= d >= 1) + p::pType # params vector + k::kType # knot vector + c::cType # control points + pVecType::Symbol + knotVecType::Symbol + function BSplineApprox{FT}(u, t, d, h, p, k, c, pVecType, knotVecType) where {FT} + new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), FT, eltype(u)}(u, + t, + d, + h, + p, + k, + c, + pVecType, + knotVecType) + end end -function BSplineApprox(u,t,d,h,pVecType,knotVecType) - u, t = munge_data(u, t) - n = length(t) - 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)) - - for i = 2:n - s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2) - l[i-1] = s - end - if pVecType == :Uniform - for i = 2:(n-1) - p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1) +function BSplineApprox(u, t, d, h, pVecType, knotVecType) + u, t = munge_data(u, t) + n = length(t) + 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)) + + for i in 2:n + s += √((t[i] - t[i - 1])^2 + (u[i] - u[i - 1])^2) + l[i - 1] = s end - elseif pVecType == :ArcLen - for i = 2:(n-1) - p[i] = p[1] + l[i-1]/s * (p[end]-p[1]) + 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 - 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 = 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 = (d+2):h - k[i] = k[1] + (i-d-1)//(h-d) * (k[end]-k[1]) + + 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 - 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] + # control points + c = zeros(eltype(u), h) + c[1] = u[1] + c[end] = u[end] + q = zeros(eltype(u), n) + N = zeros(eltype(t), n, h) + for i in 1:n + N[i, :] .= spline_coefficients(h, d, k, p[i]) end - for i = (d+3):h - k[i] = 1//d * (ps[idx+d] - ps[idx]) - idx += 1 + for k in 2:(n - 1) + q[k] = u[k] - N[k, 1] * u[1] - N[k, h] * u[end] end - end - # control points - c = zeros(eltype(u), h) - c[1] = u[1] - c[end] = u[end] - q = zeros(eltype(u), n) - N = zeros(eltype(t), n, h) - for i = 1:n - N[i, :] .= spline_coefficients(h,d,k,p[i]) - end - for k = 2:n-1 - q[k] = u[k] - N[k, 1]*u[1] - N[k, h]*u[end] - end - Q = Matrix{eltype(u)}(undef, h-2, 1) - for i = 2:h-1 - s = 0.0 - for k = 2:n-1 - s += N[k, i]*q[k] + 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] + end + Q[i - 1] = s end - Q[i-1] = s - end - N = N[2:end-1,2:h-1] - M = transpose(N) * N - P = M\Q - c[2:end-1] .= vec(P) - BSplineApprox{true}(u,t,d,h,p,k,c,pVecType,knotVecType) + N = N[2:(end - 1), 2:(h - 1)] + M = transpose(N) * N + P = M \ Q + c[2:(end - 1)] .= vec(P) + BSplineApprox{true}(u, t, d, h, p, k, c, pVecType, knotVecType) end diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index ba543cfb..093aaa6d 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1,19 +1,19 @@ -_interpolate(interp, t) = _interpolate(interp, t, firstindex(interp.t)-1)[1] +_interpolate(interp, t) = _interpolate(interp, t, firstindex(interp.t) - 1)[1] # Linear Interpolation function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess) if isnan(t) # For correct derivative with NaN - idx = firstindex(A)-1 + idx = firstindex(A) - 1 t1 = t2 = one(eltype(A.t)) u1 = u2 = one(eltype(A.u)) else idx = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) - t1, t2 = A.t[idx], A.t[idx+1] - u1, u2 = A.u[idx], A.u[idx+1] + t1, t2 = A.t[idx], A.t[idx + 1] + u1, u2 = A.u[idx], A.u[idx + 1] end - θ = (t - t1)/(t2 - t1) - val = (1 - θ)*u1 + θ*u2 + θ = (t - t1) / (t2 - t1) + val = (1 - θ) * u1 + θ * u2 # Note: The following is limited to when val is NaN as to not change the derivative of exact points. t == t1 && any(isnan, val) && return oftype(val, u1), idx # Return exact value if no interpolation needed (eg when NaN at t2) t == t2 && any(isnan, val) && return oftype(val, u2), idx # ... (eg when NaN at t1) @@ -21,170 +21,179 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues end function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess) - idx = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) - θ = (t - A.t[idx])/(A.t[idx + 1] - A.t[idx]) - return (1 - θ)*A.u[:,idx] + θ*A.u[:,idx + 1], idx + idx = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) + θ = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) + return (1 - θ) * A.u[:, idx] + θ * A.u[:, idx + 1], idx end # Quadratic Interpolation -_quad_interp_indices(A, t) = _quad_interp_indices(A, t, firstindex(A.t)-1) +_quad_interp_indices(A, t) = _quad_interp_indices(A, t, firstindex(A.t) - 1) function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess) - inner_idx = searchsortedlastcorrelated(A.t, t, iguess) - A.mode == :Backward && (inner_idx -= 1) - idx = max(1, min(inner_idx, length(A.t) - 2)) - idx, idx + 1, idx + 2 + inner_idx = searchsortedlastcorrelated(A.t, t, iguess) + A.mode == :Backward && (inner_idx -= 1) + idx = max(1, min(inner_idx, length(A.t) - 2)) + idx, idx + 1, idx + 2 end function _interpolate(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess) - i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - l₀ = ((t - A.t[i₁])*(t - A.t[i₂]))/((A.t[i₀] - A.t[i₁])*(A.t[i₀] - A.t[i₂])) - l₁ = ((t - A.t[i₀])*(t - A.t[i₂]))/((A.t[i₁] - A.t[i₀])*(A.t[i₁] - A.t[i₂])) - l₂ = ((t - A.t[i₀])*(t - A.t[i₁]))/((A.t[i₂] - A.t[i₀])*(A.t[i₂] - A.t[i₁])) - return A.u[i₀]*l₀ + A.u[i₁]*l₁ + A.u[i₂]*l₂, i₀ + i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) + l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) + l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) + l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) + return A.u[i₀] * l₀ + A.u[i₁] * l₁ + A.u[i₂] * l₂, i₀ end function _interpolate(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess) - i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) - l₀ = ((t - A.t[i₁])*(t - A.t[i₂]))/((A.t[i₀] - A.t[i₁])*(A.t[i₀] - A.t[i₂])) - l₁ = ((t - A.t[i₀])*(t - A.t[i₂]))/((A.t[i₁] - A.t[i₀])*(A.t[i₁] - A.t[i₂])) - l₂ = ((t - A.t[i₀])*(t - A.t[i₁]))/((A.t[i₂] - A.t[i₀])*(A.t[i₂] - A.t[i₁])) - return A.u[:,i₀]*l₀ + A.u[:,i₁]*l₁ + A.u[:,i₂]*l₂, i₀ + i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess) + l₀ = ((t - A.t[i₁]) * (t - A.t[i₂])) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂])) + l₁ = ((t - A.t[i₀]) * (t - A.t[i₂])) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂])) + l₂ = ((t - A.t[i₀]) * (t - A.t[i₁])) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁])) + return A.u[:, i₀] * l₀ + A.u[:, i₁] * l₁ + A.u[:, i₂] * l₂, i₀ end # Lagrange Interpolation function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number) - idxs = findRequiredIdxs(A,t) - if A.t[idxs[1]] == t - return A.u[idxs[1]] - end - N = zero(A.u[1]); D = zero(A.t[1]); tmp = N - for i = 1:length(idxs) - if isnan(A.bcache[idxs[i]]) - mult = one(A.t[1]) - for j = 1:(i-1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - for j = (i+1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - A.bcache[idxs[i]] = mult - else - mult = A.bcache[idxs[i]] + idxs = findRequiredIdxs(A, t) + if A.t[idxs[1]] == t + return A.u[idxs[1]] + end + N = zero(A.u[1]) + D = zero(A.t[1]) + tmp = N + for i in 1:length(idxs) + if isnan(A.bcache[idxs[i]]) + mult = one(A.t[1]) + for j in 1:(i - 1) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + for j in (i + 1):length(idxs) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + A.bcache[idxs[i]] = mult + else + mult = A.bcache[idxs[i]] + end + tmp = inv((t - A.t[idxs[i]]) * mult) + D += tmp + N += (tmp * A.u[idxs[i]]) end - tmp = inv((t - A.t[idxs[i]]) * mult) - D += tmp - N += (tmp * A.u[idxs[i]]) - end - N/D + N / D end function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) - idxs = findRequiredIdxs(A,t) - if A.t[idxs[1]] == t - return A.u[:,idxs[1]] - end - N = zero(A.u[:,1]); D = zero(A.t[1]); tmp = D - for i = 1:length(idxs) - if isnan(A.bcache[idxs[i]]) - mult = one(A.t[1]) - for j = 1:(i-1) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - for j = (i+1):length(idxs) - mult *= (A.t[idxs[i]] - A.t[idxs[j]]) - end - A.bcache[idxs[i]] = mult - else - mult = A.bcache[idxs[i]] + idxs = findRequiredIdxs(A, t) + if A.t[idxs[1]] == t + return A.u[:, idxs[1]] end - tmp = inv((t - A.t[idxs[i]]) * mult) - D += tmp - @. N += (tmp * A.u[:,idxs[i]]) - end - N/D + N = zero(A.u[:, 1]) + D = zero(A.t[1]) + tmp = D + for i in 1:length(idxs) + if isnan(A.bcache[idxs[i]]) + mult = one(A.t[1]) + for j in 1:(i - 1) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + for j in (i + 1):length(idxs) + mult *= (A.t[idxs[i]] - A.t[idxs[j]]) + end + A.bcache[idxs[i]] = mult + else + mult = A.bcache[idxs[i]] + end + tmp = inv((t - A.t[idxs[i]]) * mult) + D += tmp + @. N += (tmp * A.u[:, idxs[i]]) + end + N / D end -_interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) = - _interpolate(A, t), i -_interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) = - _interpolate(A, t), i +function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, i) + _interpolate(A, t), i +end +function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, i) + _interpolate(A, t), i +end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) - i = searchsortedlastcorrelated(A.t, t, iguess) - i == 0 && return A.u[1], i - i == length(A.t) && return A.u[end], i - wj = t - A.t[i] - (@evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]), i + i = searchsortedlastcorrelated(A.t, t, iguess) + i == 0 && return A.u[1], i + i == length(A.t) && return A.u[end], i + wj = t - A.t[i] + (@evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]), i end # ConstantInterpolation Interpolation function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess) - if A.dir === :left - # :left means that value to the left is used for interpolation - i = max(1, searchsortedlastcorrelated(A.t, t, iguess)) - return A.u[i], i - else - # :right means that value to the right is used for interpolation - i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess)) - return A.u[i], i - end + if A.dir === :left + # :left means that value to the left is used for interpolation + i = max(1, searchsortedlastcorrelated(A.t, t, iguess)) + return A.u[i], i + else + # :right means that value to the right is used for interpolation + i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess)) + return A.u[i], i + end end function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) - if A.dir === :left - # :left means that value to the left is used for interpolation - i = max(1, searchsortedlastcorrelated(A.t, t, iguess)) - return A.u[:, i], i - else - # :right means that value to the right is used for interpolation - i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess)) - return A.u[:, i], i - end + if A.dir === :left + # :left means that value to the left is used for interpolation + i = max(1, searchsortedlastcorrelated(A.t, t, iguess)) + return A.u[:, i], i + else + # :right means that value to the right is used for interpolation + i = min(length(A.t), searchsortedfirstcorrelated(A.t, t, iguess)) + return A.u[:, i], i + end end # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - i = min(max(2, searchsortedfirstcorrelated(A.t, t, iguess)), length(A.t)) - Cᵢ = A.u[i-1] - σ = 1//2 * (A.z[i] - A.z[i-1])/(A.t[i] - A.t[i-1]) - return A.z[i-1] * (t - A.t[i-1]) + σ * (t - A.t[i-1])^2 + Cᵢ, i + i = min(max(2, searchsortedfirstcorrelated(A.t, t, iguess)), length(A.t)) + Cᵢ = A.u[i - 1] + σ = 1 // 2 * (A.z[i] - A.z[i - 1]) / (A.t[i] - A.t[i - 1]) + return A.z[i - 1] * (t - A.t[i - 1]) + σ * (t - A.t[i - 1])^2 + Cᵢ, i end # CubicSpline Interpolation function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) - i = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) - I = A.z[i] * (A.t[i+1] - t)^3 / (6A.h[i+1]) + A.z[i+1] * (t - A.t[i])^3 / (6A.h[i+1]) - C = (A.u[i+1]/A.h[i+1] - A.z[i+1]*A.h[i+1]/6)*(t - A.t[i]) - D = (A.u[i]/A.h[i+1] - A.z[i]*A.h[i+1]/6)*(A.t[i+1] - t) - I + C + D, i + i = max(1, min(searchsortedlastcorrelated(A.t, t, iguess), length(A.t) - 1)) + I = A.z[i] * (A.t[i + 1] - t)^3 / (6A.h[i + 1]) + + A.z[i + 1] * (t - A.t[i])^3 / (6A.h[i + 1]) + C = (A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6) * (t - A.t[i]) + D = (A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6) * (A.t[i + 1] - t) + I + C + D, i end # BSpline Curve Interpolation -function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) - # change t into param [0 1] - idx = searchsortedlastcorrelated(A.t,t, iguess) - idx == length(A.t) ? idx -= 1 : nothing - 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 = spline_coefficients(n,A.d,A.k,t) - ucum = zero(eltype(A.u)) - for i = 1:n - ucum += N[i] * A.c[i] - end - ucum, idx +function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, + t::Number, + iguess) + # change t into param [0 1] + idx = searchsortedlastcorrelated(A.t, t, iguess) + idx == length(A.t) ? idx -= 1 : nothing + 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 = spline_coefficients(n, A.d, A.k, t) + ucum = zero(eltype(A.u)) + for i in 1:n + ucum += N[i] * A.c[i] + end + ucum, idx end # BSpline Curve Approx function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess) - # change t into param [0 1] - idx = searchsortedlastcorrelated(A.t,t, iguess) - idx == length(A.t) ? idx -= 1 : nothing - 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 = spline_coefficients(A.h,A.d,A.k,t) - ucum = zero(eltype(A.u)) - for i = 1:A.h - ucum += N[i] * A.c[i] - end - ucum, idx + # change t into param [0 1] + idx = searchsortedlastcorrelated(A.t, t, iguess) + idx == length(A.t) ? idx -= 1 : nothing + 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 = spline_coefficients(A.h, A.d, A.k, t) + ucum = zero(eltype(A.u)) + for i in 1:(A.h) + ucum += N[i] * A.c[i] + end + ucum, idx end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index c3f6e825..17bae864 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -1,125 +1,125 @@ -function findRequiredIdxs(A::LagrangeInterpolation, t) - idxs = sortperm(A.t,by=x->abs(t-x)) - idxs[1:(A.n+1)] -end - -function spline_coefficients(n, d, k, u::Number) - N = zeros(n) - if u == k[1] - N[1] = one(u) - elseif u == k[end] - N[end] = one(u) - else - i = findfirst(x->x>u,k) - 1 - N[i] = one(u) - for deg = 1:d - N[i-deg] = (k[i+1]-u)/(k[i+1]-k[i-deg+1]) * N[i-deg+1] - for j = (i-deg+1):(i-1) - N[j] = (u-k[j])/(k[j+deg]-k[j]) * N[j] + (k[j+deg+1]-u)/(k[j+deg+1]-k[j+1]) * N[j+1] - end - N[i] = (u-k[i])/(k[i+deg]-k[i]) * N[i] - end - end - N -end - -function spline_coefficients(n, d, k, u::AbstractVector) - N = zeros(eltype(u),n,n) - for i = 1:n - N[i,:] .= spline_coefficients(n,d,k,u[i]) - end - N -end - -# helper function for data manipulation -function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) - return u, t -end - -function munge_data(u::AbstractVector, t::AbstractVector) - Tu = Base.nonmissingtype(eltype(u)) - Tt = Base.nonmissingtype(eltype(t)) - @assert length(t) == length(u) - non_missing_indices = collect(i for i in 1:length(t) if !ismissing(u[i]) && !ismissing(t[i])) - newu = Tu.([u[i] for i in non_missing_indices]) - newt = Tt.([t[i] for i in non_missing_indices]) - - return newu, newt -end - -function munge_data(U::StridedMatrix, t::AbstractVector) - TU = Base.nonmissingtype(eltype(U)) - Tt = Base.nonmissingtype(eltype(t)) - @assert length(t) == size(U,2) - non_missing_indices = collect(i for i in 1:length(t) if !any(ismissing,U[:,i]) && !ismissing(t[i])) - newUs = [TU.(U[:,i]) for i in non_missing_indices] - newt= Tt.([t[i] for i in non_missing_indices]) - - return hcat(newUs...), newt -end - - -""" - bracketstrictlymontonic(v, x, guess; lt=, by=, rev=false) - -Starting from an initial `guess` index, find indices `(lo, hi)` such that `v[lo] ≤ x ≤ -v[hi]` according to the specified order, assuming that `x` is actually within the range of -values found in `v`. If `x` is outside that range, either `lo` will be `firstindex(v)` or -`hi` will be `lastindex(v)`. - -Note that the results will not typically satify `lo ≤ guess ≤ hi`. If `x` is precisely -equal to a value that is not unique in the input `v`, there is no guarantee that `(lo, hi)` -will encompass *all* indices corresponding to that value. - -This algorithm is essentially an expanding binary search, which can be used as a precursor -to `searchsorted` and related functions, which can take `lo` and `hi` as arguments. The -purpose of using this function first would be to accelerate convergence in those functions -by using correlated `guess`es for repeated calls. The best `guess` for the next call of -this function would be the index returned by the previous call to `searchsorted`. - -See [`sort!`](@ref) for an explanation of the keyword arguments `by`, `lt` and `rev`. -""" -function bracketstrictlymontonic( - v::AbstractVector, - x, - guess::T, - o::Base.Order.Ordering -)::NTuple{2,keytype(v)} where {T<:Integer} - bottom = firstindex(v) - top = lastindex(v) - if guess < bottom || guess > top - return bottom, top - # # NOTE: for cache efficiency in repeated calls, we avoid accessing the first and last elements of `v` - # # on each call to this function. This should only result in significant slow downs for calls with - # # out-of-bounds values of `x` *and* bad `guess`es. - # elseif lt(o, x, v[bottom]) - # return bottom, bottom - # elseif lt(o, v[top], x) - # return top, top - else - u = T(1) - lo, hi = guess, min(guess + u, top) - @inbounds if Base.Order.lt(o, x, v[lo]) - while lo > bottom && Base.Order.lt(o, x, v[lo]) - lo, hi = max(bottom, lo - u), lo - u += u - end - else - while hi < top && !Base.Order.lt(o, x, v[hi]) - lo, hi = hi, min(top, hi + u) - u += u - end - end - end - return lo, hi -end - -function searchsortedfirstcorrelated(v::AbstractVector, x, guess) - lo, hi = bracketstrictlymontonic(v, x, guess, Base.Order.Forward) - searchsortedfirst(v, x, lo, hi, Base.Order.Forward) -end - -function searchsortedlastcorrelated(v::AbstractVector, x, guess) - lo, hi = bracketstrictlymontonic(v, x, guess, Base.Order.Forward) - searchsortedlast(v, x, lo, hi, Base.Order.Forward) -end +function findRequiredIdxs(A::LagrangeInterpolation, t) + idxs = sortperm(A.t, by = x -> abs(t - x)) + idxs[1:(A.n + 1)] +end + +function spline_coefficients(n, d, k, u::Number) + N = zeros(n) + if u == k[1] + N[1] = one(u) + elseif u == k[end] + N[end] = one(u) + else + i = findfirst(x -> x > u, k) - 1 + N[i] = one(u) + for deg in 1:d + N[i - deg] = (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) * N[i - deg + 1] + for j in (i - deg + 1):(i - 1) + N[j] = (u - k[j]) / (k[j + deg] - k[j]) * N[j] + + (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * N[j + 1] + end + N[i] = (u - k[i]) / (k[i + deg] - k[i]) * N[i] + end + end + N +end + +function spline_coefficients(n, d, k, u::AbstractVector) + N = zeros(eltype(u), n, n) + for i in 1:n + N[i, :] .= spline_coefficients(n, d, k, u[i]) + end + N +end + +# helper function for data manipulation +function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real}) + return u, t +end + +function munge_data(u::AbstractVector, t::AbstractVector) + Tu = Base.nonmissingtype(eltype(u)) + Tt = Base.nonmissingtype(eltype(t)) + @assert length(t) == length(u) + non_missing_indices = collect(i for i in 1:length(t) + if !ismissing(u[i]) && !ismissing(t[i])) + newu = Tu.([u[i] for i in non_missing_indices]) + newt = Tt.([t[i] for i in non_missing_indices]) + + return newu, newt +end + +function munge_data(U::StridedMatrix, t::AbstractVector) + TU = Base.nonmissingtype(eltype(U)) + Tt = Base.nonmissingtype(eltype(t)) + @assert length(t) == size(U, 2) + non_missing_indices = collect(i for i in 1:length(t) + if !any(ismissing, U[:, i]) && !ismissing(t[i])) + newUs = [TU.(U[:, i]) for i in non_missing_indices] + newt = Tt.([t[i] for i in non_missing_indices]) + + return hcat(newUs...), newt +end + +""" + bracketstrictlymontonic(v, x, guess; lt=, by=, rev=false) + +Starting from an initial `guess` index, find indices `(lo, hi)` such that `v[lo] ≤ x ≤ +v[hi]` according to the specified order, assuming that `x` is actually within the range of +values found in `v`. If `x` is outside that range, either `lo` will be `firstindex(v)` or +`hi` will be `lastindex(v)`. + +Note that the results will not typically satify `lo ≤ guess ≤ hi`. If `x` is precisely +equal to a value that is not unique in the input `v`, there is no guarantee that `(lo, hi)` +will encompass *all* indices corresponding to that value. + +This algorithm is essentially an expanding binary search, which can be used as a precursor +to `searchsorted` and related functions, which can take `lo` and `hi` as arguments. The +purpose of using this function first would be to accelerate convergence in those functions +by using correlated `guess`es for repeated calls. The best `guess` for the next call of +this function would be the index returned by the previous call to `searchsorted`. + +See [`sort!`](@ref) for an explanation of the keyword arguments `by`, `lt` and `rev`. +""" +function bracketstrictlymontonic(v::AbstractVector, + x, + guess::T, + o::Base.Order.Ordering)::NTuple{2, keytype(v)} where {T <: Integer} + bottom = firstindex(v) + top = lastindex(v) + if guess < bottom || guess > top + return bottom, top + # # NOTE: for cache efficiency in repeated calls, we avoid accessing the first and last elements of `v` + # # on each call to this function. This should only result in significant slow downs for calls with + # # out-of-bounds values of `x` *and* bad `guess`es. + # elseif lt(o, x, v[bottom]) + # return bottom, bottom + # elseif lt(o, v[top], x) + # return top, top + else + u = T(1) + lo, hi = guess, min(guess + u, top) + @inbounds if Base.Order.lt(o, x, v[lo]) + while lo > bottom && Base.Order.lt(o, x, v[lo]) + lo, hi = max(bottom, lo - u), lo + u += u + end + else + while hi < top && !Base.Order.lt(o, x, v[hi]) + lo, hi = hi, min(top, hi + u) + u += u + end + end + end + return lo, hi +end + +function searchsortedfirstcorrelated(v::AbstractVector, x, guess) + lo, hi = bracketstrictlymontonic(v, x, guess, Base.Order.Forward) + searchsortedfirst(v, x, lo, hi, Base.Order.Forward) +end + +function searchsortedlastcorrelated(v::AbstractVector, x, guess) + lo, hi = bracketstrictlymontonic(v, x, guess, Base.Order.Forward) + searchsortedlast(v, x, lo, hi, Base.Order.Forward) +end diff --git a/src/online.jl b/src/online.jl index cdf09c0d..4ad06f3a 100644 --- a/src/online.jl +++ b/src/online.jl @@ -1,40 +1,40 @@ import Base: append!, push! function push!(di::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di + push!(di.u, u) + push!(di.t, t) + di end function push!(di::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di + push!(di.u, u) + push!(di.t, t) + di end function push!(di::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T} - push!(di.u, u) - push!(di.t, t) - di + push!(di.u, u) + push!(di.t, t) + di end function append!(di::LinearInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di + u, t = munge_data(u, t) + append!(di.u, u) + append!(di.t, t) + di end function append!(di::QuadraticInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di + u, t = munge_data(u, t) + append!(di.u, u) + append!(di.t, t) + di end function append!(di::ConstantInterpolation{U, T}, u::U, t::T) where {U, T} - u, t = munge_data(u, t) - append!(di.u, u) - append!(di.t, t) - di -end \ No newline at end of file + u, t = munge_data(u, t) + append!(di.u, u) + append!(di.t, t) + di +end diff --git a/src/plot_rec.jl b/src/plot_rec.jl index 28f5f8b8..d801bfa5 100644 --- a/src/plot_rec.jl +++ b/src/plot_rec.jl @@ -1,251 +1,238 @@ -################################################################################ -# Type recipes # -################################################################################ - -function to_plottable(A::AbstractInterpolation; plotdensity = 10_000, denseplot = true) - t = sort(A.t) - start = t[1]; stop = t[end] - if denseplot - plott = collect(range(start,stop=stop,length=plotdensity)) - else - plott = t - end - output = A.(plott) - plott, output -end - -@recipe function f(A::AbstractInterpolation; plotdensity = 10_000, denseplot = true) - to_plottable(A; plotdensity = plotdensity, denseplot=denseplot) -end - -################################################################################ -# Series recipes # -################################################################################ - -############################################################ -# Interpolations # -############################################################ - -######################################## -# Linear Interpolation # -######################################## - -@recipe function f(::Type{Val{:linear_interp}}, x, y, z; plotdensity = 10_000, denseplot = true) - - seriestype := :path - - label --> "Linear fit" - - nx, ny = to_plottable(LinearInterpolation(y, x); plotdensity = plotdensity, denseplot = denseplot) - - x := nx - y := ny -end - -######################################## -# Quadratic Interpolation # -######################################## - -@recipe function f(::Type{Val{:quadratic_interp}}, x, y, z; plotdensity = 10_000, denseplot = true) - - seriestype := :path - - label --> "Quadratic fit" - - nx, ny = to_plottable( - QuadraticInterpolation( - T.(y), - T.(x) - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - x := nx - y := ny -end - -######################################## -# Quadratic Spline # -######################################## - -@recipe function f(::Type{Val{:quadratic_spline}}, x, y, z; plotdensity = 10_000, denseplot = true) - - seriestype := :path - - label --> "Quadratic Spline" - - T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - QuadraticSpline( - T.(y), - T.(x) - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - - x := nx - y := ny -end - -######################################## -# Lagrange Interpolation # -######################################## - -@recipe function f(::Type{Val{:lagrange_interp}}, - x, y, z; - n = length(x) - 1, - plotdensity = 10_000, - denseplot = true - ) - - seriestype := :path - - label --> "Lagrange Fit" - - T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - LagrangeInterpolation( - T.(y), - T.(x), - n - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - - x := nx - y := ny -end - -######################################## -# Cubic Spline # -######################################## - -@recipe function f(::Type{Val{:cubic_spline}}, x, y, z; plotdensity = 10_000, denseplot = true) - - seriestype := :path - - label --> "Cubic Spline" - - T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - CubicSpline( - T.(y), - T.(x) - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - x := nx - y := ny -end - - -@recipe function f(::Type{Val{:bspline_interp}}, - x, y, z; - d = 5, - pVec=:ArcLen, - knotVec = :Average, - plotdensity = length(x) * 6, - denseplot = true - ) - - seriestype := :path - - label --> "B-Spline" - - @show x y eltype(x) - - # T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - BSplineInterpolation( - T.(y), - T.(x), - d, - pVec, - knotVec - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - x := nx - y := ny -end - -######################################## -# B-spline (approximation) # -######################################## - -@recipe function f(::Type{Val{:bspline_approx}}, - x, y, z; - d = 5, - h = length(x) - 1, - pVec = :ArcLen, - knotVec = :Average, - plotdensity = length(x) * 6, - denseplot = true - ) - - seriestype := :path - - label --> "B-Spline" - - T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - BSplineApprox( - T.(y), - T.(x), - d, - h, - pVec, - knotVec - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - x := nx - y := ny -end - -######################################## -# Akima intepolation # -######################################## - -@recipe function f(::Type{Val{:akima}}, x, y, z; plotdensity = length(x) * 6, denseplot = true) - - seriestype := :path - - label --> "Akima" - - T = promote_type(eltype(y), eltype(x)) - - nx, ny = to_plottable( - AkimaInterpolation( - T.(y), - T.(x) - ); - plotdensity = plotdensity, - denseplot = denseplot - ) - x := nx - y := ny -end - - -################################################################################ -# Shorthands # -################################################################################ - -""" - akima(u, t) - akima!(u, t) - -Plot the Akima interpolation of the given data. -""" -@shorthands akima +################################################################################ +# Type recipes # +################################################################################ + +function to_plottable(A::AbstractInterpolation; plotdensity = 10_000, denseplot = true) + t = sort(A.t) + start = t[1] + stop = t[end] + if denseplot + plott = collect(range(start, stop = stop, length = plotdensity)) + else + plott = t + end + output = A.(plott) + plott, output +end + +@recipe function f(A::AbstractInterpolation; plotdensity = 10_000, denseplot = true) + to_plottable(A; plotdensity = plotdensity, denseplot = denseplot) +end + +################################################################################ +# Series recipes # +################################################################################ + +############################################################ +# Interpolations # +############################################################ + +######################################## +# Linear Interpolation # +######################################## + +@recipe function f(::Type{Val{:linear_interp}}, + x, + y, + z; + plotdensity = 10_000, + denseplot = true) + seriestype := :path + + label --> "Linear fit" + + nx, ny = to_plottable(LinearInterpolation(y, x); + plotdensity = plotdensity, + denseplot = denseplot) + + x := nx + y := ny +end + +######################################## +# Quadratic Interpolation # +######################################## + +@recipe function f(::Type{Val{:quadratic_interp}}, + x, + y, + z; + plotdensity = 10_000, + denseplot = true) + seriestype := :path + + label --> "Quadratic fit" + + nx, ny = to_plottable(QuadraticInterpolation(T.(y), + T.(x)); + plotdensity = plotdensity, + denseplot = denseplot) + x := nx + y := ny +end + +######################################## +# Quadratic Spline # +######################################## + +@recipe function f(::Type{Val{:quadratic_spline}}, + x, + y, + z; + plotdensity = 10_000, + denseplot = true) + seriestype := :path + + label --> "Quadratic Spline" + + T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(QuadraticSpline(T.(y), + T.(x)); + plotdensity = plotdensity, + denseplot = denseplot) + + x := nx + y := ny +end + +######################################## +# Lagrange Interpolation # +######################################## + +@recipe function f(::Type{Val{:lagrange_interp}}, + x, y, z; + n = length(x) - 1, + plotdensity = 10_000, + denseplot = true) + seriestype := :path + + label --> "Lagrange Fit" + + T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(LagrangeInterpolation(T.(y), + T.(x), + n); + plotdensity = plotdensity, + denseplot = denseplot) + + x := nx + y := ny +end + +######################################## +# Cubic Spline # +######################################## + +@recipe function f(::Type{Val{:cubic_spline}}, + x, + y, + z; + plotdensity = 10_000, + denseplot = true) + seriestype := :path + + label --> "Cubic Spline" + + T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(CubicSpline(T.(y), + T.(x)); + plotdensity = plotdensity, + denseplot = denseplot) + x := nx + y := ny +end + +@recipe function f(::Type{Val{:bspline_interp}}, + x, y, z; + d = 5, + pVec = :ArcLen, + knotVec = :Average, + plotdensity = length(x) * 6, + denseplot = true) + seriestype := :path + + label --> "B-Spline" + + @show x y eltype(x) + + # T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(BSplineInterpolation(T.(y), + T.(x), + d, + pVec, + knotVec); + plotdensity = plotdensity, + denseplot = denseplot) + x := nx + y := ny +end + +######################################## +# B-spline (approximation) # +######################################## + +@recipe function f(::Type{Val{:bspline_approx}}, + x, y, z; + d = 5, + h = length(x) - 1, + pVec = :ArcLen, + knotVec = :Average, + plotdensity = length(x) * 6, + denseplot = true) + seriestype := :path + + label --> "B-Spline" + + T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(BSplineApprox(T.(y), + T.(x), + d, + h, + pVec, + knotVec); + plotdensity = plotdensity, + denseplot = denseplot) + x := nx + y := ny +end + +######################################## +# Akima intepolation # +######################################## + +@recipe function f(::Type{Val{:akima}}, + x, + y, + z; + plotdensity = length(x) * 6, + denseplot = true) + seriestype := :path + + label --> "Akima" + + T = promote_type(eltype(y), eltype(x)) + + nx, ny = to_plottable(AkimaInterpolation(T.(y), + T.(x)); + plotdensity = plotdensity, + denseplot = denseplot) + x := nx + y := ny +end + +################################################################################ +# Shorthands # +################################################################################ + +""" + akima(u, t) + akima!(u, t) + +Plot the Akima interpolation of the given data. +""" +@shorthands akima diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 61afa197..3b3a1518 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -3,86 +3,86 @@ using FiniteDifferences using DataInterpolations: derivative function test_derivatives(func, tspan, name::String) - trange = range(minimum(tspan), maximum(tspan), length=32)[2:end-1] - @testset "$name" begin - for t in trange - # Linearly spaced points might lead to evaluations outside - # trange - cdiff = central_fdm(5, 1; geom=true)(_t -> func(_t), t) - adiff = derivative(func, t) - @test isapprox(cdiff, adiff, atol=1e-8) + trange = range(minimum(tspan), maximum(tspan), length = 32)[2:(end - 1)] + @testset "$name" begin + for t in trange + # Linearly spaced points might lead to evaluations outside + # trange + cdiff = central_fdm(5, 1; geom = true)(_t -> func(_t), t) + adiff = derivative(func, t) + @test isapprox(cdiff, adiff, atol = 1e-8) + end end - end end # Linear Interpolation u = 2.0collect(1:10) t = 1.0collect(1:10) -A = LinearInterpolation(u,t) +A = LinearInterpolation(u, t) test_derivatives(A, t, "Linear Interpolation (Vector)") u = vcat(2.0collect(1:10)', 3.0collect(1:10)') -A = LinearInterpolation(u,t) +A = LinearInterpolation(u, t) test_derivatives(A, t, "Linear Interpolation (Matrix)") # Quadratic Interpolation u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] -A = QuadraticInterpolation(u,t) +A = QuadraticInterpolation(u, t) test_derivatives(A, t, "Quadratic Interpolation (Vector)") -Ab = QuadraticInterpolation(u,t,:Backward) +Ab = QuadraticInterpolation(u, t, :Backward) test_derivatives(Ab, t, "Quadratic Interpolation (Vector), backward") u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0] -A = QuadraticInterpolation(u,t) +A = QuadraticInterpolation(u, t) test_derivatives(A, t, "Quadratic Interpolation (Matrix)") @testset "Backward Quadratic Interpolation" begin - u = [0.5, 0.0, 0.5, 0.0] - t = [1.0, 2.0, 3.0, 4.0] - A_f = QuadraticInterpolation(u,t) - A_b = QuadraticInterpolation(u,t,:Backward) - @test derivative(A_f, 1.5) ≈ -0.5 - @test derivative(A_b, 1.5) ≈ -0.5 - @test derivative(A_f, 2.25) ≈ 0.75 - @test derivative(A_b, 2.25) ≈ 0.25 - @test derivative(A_f, 2.75) ≈ 0.25 - @test derivative(A_b, 2.75) ≈ 0.75 - @test derivative(A_f, 3.5) ≈ -0.5 - @test derivative(A_b, 3.5) ≈ -0.5 + u = [0.5, 0.0, 0.5, 0.0] + t = [1.0, 2.0, 3.0, 4.0] + A_f = QuadraticInterpolation(u, t) + A_b = QuadraticInterpolation(u, t, :Backward) + @test derivative(A_f, 1.5) ≈ -0.5 + @test derivative(A_b, 1.5) ≈ -0.5 + @test derivative(A_f, 2.25) ≈ 0.75 + @test derivative(A_b, 2.25) ≈ 0.25 + @test derivative(A_f, 2.75) ≈ 0.25 + @test derivative(A_b, 2.75) ≈ 0.75 + @test derivative(A_f, 3.5) ≈ -0.5 + @test derivative(A_b, 3.5) ≈ -0.5 end # Lagrange Interpolation u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] -A = LagrangeInterpolation(u,t) +A = LagrangeInterpolation(u, t) test_derivatives(A, t, "Lagrange Interpolation (Vector)") # Lagrange Interpolation u = [1.0 4.0 9.0; 1.0 2.0 3.0] t = [1.0, 2.0, 3.0] -A = LagrangeInterpolation(u,t) +A = LagrangeInterpolation(u, t) test_derivatives(A, t, "Lagrange Interpolation (Matrix)") # Lagrange Interpolation u = [[1.0, 4.0, 9.0], [3.0, 7.0, 4.0], [5.0, 4.0, 1.0]] t = [1.0, 2.0, 3.0] -A = LagrangeInterpolation(u,t) +A = LagrangeInterpolation(u, t) test_derivatives(A, t, "Lagrange Interpolation (Vector of Vectors)") # Lagrange Interpolation u = [[3.0 1.0 4.0; 1.0 5.0 9.0], [2.0 6.0 5.0; 3.0 5.0 8.0], [9.0 7.0 9.0; 3.0 2.0 3.0]] t = [1.0, 2.0, 3.0] -A = LagrangeInterpolation(u,t) +A = LagrangeInterpolation(u, t) test_derivatives(A, t, "Lagrange Interpolation (Vector of Matrices)") @@ -102,7 +102,7 @@ end u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] -A = QuadraticSpline(u,t) +A = QuadraticSpline(u, t) test_derivatives(A, t, "Quadratic Interpolation (Vector)") @@ -110,7 +110,7 @@ test_derivatives(A, t, "Quadratic Interpolation (Vector)") u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]] t = [-1.0, 0.0, 1.0] -A = QuadraticSpline(u,t) +A = QuadraticSpline(u, t) test_derivatives(A, t, "Quadratic Interpolation (Vector of Vectors)") @@ -118,7 +118,7 @@ test_derivatives(A, t, "Quadratic Interpolation (Vector of Vectors)") u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]] t = [-1.0, 0.0, 1.0] -A = QuadraticSpline(u,t) +A = QuadraticSpline(u, t) test_derivatives(A, t, "Quadratic Interpolation (Vector of Matrices)") @@ -126,7 +126,7 @@ test_derivatives(A, t, "Quadratic Interpolation (Vector of Matrices)") u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] -A = CubicSpline(u,t) +A = CubicSpline(u, t) test_derivatives(A, t, "Cubic Spline Interpolation (Vector)") @@ -134,7 +134,7 @@ test_derivatives(A, t, "Cubic Spline Interpolation (Vector)") u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]] t = [-1.0, 0.0, 1.0] -A = CubicSpline(u,t) +A = CubicSpline(u, t) test_derivatives(A, t, "Cubic Spline Interpolation (Vector of Vectors)") @@ -142,23 +142,23 @@ test_derivatives(A, t, "Cubic Spline Interpolation (Vector of Vectors)") u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]] t = [-1.0, 0.0, 1.0] -A = CubicSpline(u,t) +A = CubicSpline(u, t) test_derivatives(A, t, "Cubic Spline Interpolation (Vector of Matrices)") # 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] +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 = BSplineInterpolation(u,t,2,:Uniform,:Uniform) +A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) test_derivatives(A, t, "BSpline Interpolation (Uniform, Uniform)") -A = BSplineInterpolation(u,t,2,:ArcLen,:Average) +A = BSplineInterpolation(u, t, 2, :ArcLen, :Average) test_derivatives(A, t, "BSpline Interpolation (Arclen, Average)") -A = BSplineApprox(u,t,3,4,:Uniform,:Uniform) +A = BSplineApprox(u, t, 3, 4, :Uniform, :Uniform) test_derivatives(A, t, "BSpline Approx (Uniform, Uniform)") @@ -170,8 +170,8 @@ A = QuadraticSpline(u, t) @variables τ, ω(τ) D = Symbolics.Differential(τ) expr = A(ω) -@test isequal(Symbolics.derivative(expr, τ), D(ω)*DataInterpolations.derivative(A, ω)) +@test isequal(Symbolics.derivative(expr, τ), D(ω) * DataInterpolations.derivative(A, ω)) -derivexpr = expand_derivatives(substitute(D(A(ω)), Dict(ω=>0.5τ))) -symfunc = Symbolics.build_function(derivexpr, τ; expression=Val{false}) -@test symfunc(0.5) == 0.5*3 +derivexpr = expand_derivatives(substitute(D(A(ω)), Dict(ω => 0.5τ))) +symfunc = Symbolics.build_function(derivexpr, τ; expression = Val{false}) +@test symfunc(0.5) == 0.5 * 3 diff --git a/test/integral_tests.jl b/test/integral_tests.jl index 3927494a..d781cb6d 100644 --- a/test/integral_tests.jl +++ b/test/integral_tests.jl @@ -2,58 +2,55 @@ using DataInterpolations, Test using QuadGK using DataInterpolations: integral - function test_integral(func, tspan, name::String) - t1 = minimum(tspan) - t2 = maximum(tspan) - @testset "$name" begin - qint, err = quadgk(func, t1, t2) - aint = integral(func, t1, t2) - @test isapprox(qint, aint, atol=1e-8) - end + t1 = minimum(tspan) + t2 = maximum(tspan) + @testset "$name" begin + qint, err = quadgk(func, t1, t2) + aint = integral(func, t1, t2) + @test isapprox(qint, aint, atol = 1e-8) + end end # Linear Interpolation u = 2.0collect(1:10) t = 1.0collect(1:10) -A = LinearInterpolation(u,t) +A = LinearInterpolation(u, t) test_integral(A, t, "Linear Interpolation (Vector)") # Quadratic Interpolation u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] -A = QuadraticInterpolation(u,t) +A = QuadraticInterpolation(u, t) test_integral(A, t, "Quadratic Interpolation (Vector)") # Quadratic forward/backward Interpolation @testset "QuadraticInterpolation - forward/backward modes" begin - u = [3.0, 0.0, 3.0, 0.0] - t = [1.0, 2.0, 3.0, 4.0] - A_f = QuadraticInterpolation(u,t) - A_b = QuadraticInterpolation(u,t,:Backward) - @test integral(A_f, 1.0, 2.0) ≈ 1.0 - @test integral(A_f, 2.0, 3.0) ≈ 2.0 - @test integral(A_f, 3.0, 4.0) ≈ 2.0 - @test integral(A_f, 1.0, 3.0) ≈ 3.0 - @test integral(A_f, 2.0, 4.0) ≈ 4.0 - @test integral(A_f, 1.0, 4.0) ≈ 5.0 - @test integral(A_b, 1.0, 2.0) ≈ 1.0 - @test integral(A_b, 2.0, 3.0) ≈ 1.0 - @test integral(A_b, 3.0, 4.0) ≈ 2.0 - @test integral(A_b, 1.0, 3.0) ≈ 2.0 - @test integral(A_b, 2.0, 4.0) ≈ 3.0 - @test integral(A_b, 1.0, 4.0) ≈ 4.0 - + u = [3.0, 0.0, 3.0, 0.0] + t = [1.0, 2.0, 3.0, 4.0] + A_f = QuadraticInterpolation(u, t) + A_b = QuadraticInterpolation(u, t, :Backward) + @test integral(A_f, 1.0, 2.0) ≈ 1.0 + @test integral(A_f, 2.0, 3.0) ≈ 2.0 + @test integral(A_f, 3.0, 4.0) ≈ 2.0 + @test integral(A_f, 1.0, 3.0) ≈ 3.0 + @test integral(A_f, 2.0, 4.0) ≈ 4.0 + @test integral(A_f, 1.0, 4.0) ≈ 5.0 + @test integral(A_b, 1.0, 2.0) ≈ 1.0 + @test integral(A_b, 2.0, 3.0) ≈ 1.0 + @test integral(A_b, 3.0, 4.0) ≈ 2.0 + @test integral(A_b, 1.0, 3.0) ≈ 2.0 + @test integral(A_b, 2.0, 4.0) ≈ 3.0 + @test integral(A_b, 1.0, 4.0) ≈ 4.0 end - # QuadraticSpline Interpolation u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] -A = QuadraticSpline(u,t) +A = QuadraticSpline(u, t) test_integral(A, t, "Quadratic Spline (Vector)") @@ -61,6 +58,6 @@ test_integral(A, t, "Quadratic Spline (Vector)") u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] -A = CubicSpline(u,t) +A = CubicSpline(u, t) -test_integral(A, t, "Cubic Spline Interpolation (Vector)") \ No newline at end of file +test_integral(A, t, "Cubic Spline Interpolation (Vector)") diff --git a/test/interface.jl b/test/interface.jl index ef503448..3299104c 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -1,26 +1,26 @@ u = 2.0collect(1:10) t = 1.0collect(1:10) -A = LinearInterpolation(u,t) +A = LinearInterpolation(u, t) @test length(A) == 20 for i in 1:10 - @test u[i] == A[i] + @test u[i] == A[i] end for i in 11:20 - @test t[i-10] == A[i] + @test t[i - 10] == A[i] end -A = LinearInterpolation{false}(u,t) +A = LinearInterpolation{false}(u, t) @test length(A) == 10 for i in 1:10 - @test u[i] == A[i] + @test u[i] == A[i] end using Symbolics u = 2.0collect(1:10) t = 1.0collect(1:10) -A = LinearInterpolation(u,t) +A = LinearInterpolation(u, t) @variables t x(t) substitute(A(t), Dict(t => x)) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 2dbc256e..6ce3fefb 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -6,41 +6,41 @@ import ForwardDiff @testset "Linear Interpolation" begin u = 2.0collect(1:10) t = 1.0collect(1:10) - A = LinearInterpolation(u,t) + A = LinearInterpolation(u, t) for (_t, _u) in zip(t, u) @test A(_t) == _u end - @test A(0) == 0.0 - @test A(5.5) == 11.0 - @test A(11) == 22 + @test A(0) == 0.0 + @test A(5.5) == 11.0 + @test A(11) == 22 u = vcat(2.0collect(1:10)', 3.0collect(1:10)') - A = LinearInterpolation(u,t) + A = LinearInterpolation(u, t) for (_t, _u) in zip(t, eachcol(u)) @test A(_t) == _u end - @test A(0) == [0.0, 0.0] - @test A(5.5) == [11.0, 16.5] - @test A(11) == [22, 33] + @test A(0) == [0.0, 0.0] + @test A(5.5) == [11.0, 16.5] + @test A(11) == [22, 33] x = 1:10 y = 2:4 - u_= x' .* y - u = [u_[:,i] for i = 1:size(u_,2)] - A = LinearInterpolation(u,t) - @test A(0) == [0.0, 0.0, 0.0] - @test A(5.5) == [11.0, 16.5, 22.0] - @test A(11) == [22.0, 33.0, 44.0] - - u = [u_[:,i:i+1] for i = 1:2:10] + u_ = x' .* y + u = [u_[:, i] for i in 1:size(u_, 2)] + A = LinearInterpolation(u, t) + @test A(0) == [0.0, 0.0, 0.0] + @test A(5.5) == [11.0, 16.5, 22.0] + @test A(11) == [22.0, 33.0, 44.0] + + u = [u_[:, i:(i + 1)] for i in 1:2:10] t = 1.0collect(2:2:10) - A = LinearInterpolation(u,t) + A = LinearInterpolation(u, t) - @test A(0) == [-2.0 0.0; -3.0 0.0; -4.0 0.0] - @test A(3) == [4.0 6.0; 6.0 9.0; 8.0 12.0] - @test A(5) == [8.0 10.0; 12.0 15.0; 16.0 20.0] + @test A(0) == [-2.0 0.0; -3.0 0.0; -4.0 0.0] + @test A(3) == [4.0 6.0; 6.0 9.0; 8.0 12.0] + @test A(5) == [8.0 10.0; 12.0 15.0; 16.0 20.0] # with NaNs (#113) u = [NaN, 1.0, 2.0, 3.0] @@ -91,19 +91,19 @@ import ForwardDiff u = 1:5 t = 1:5 A2 = LinearInterpolation(u, t) - u = [1//i for i in 1:5] + u = [1 // i for i in 1:5] t = (1:5) A3 = LinearInterpolation(u, t) - u = [1//i for i in 1:5] - t = [1//(6-i) for i in 1:5] + u = [1 // i for i in 1:5] + t = [1 // (6 - i) for i in 1:5] A4 = LinearInterpolation(u, t) F32 = Float32(1) F64 = Float64(1) I32 = Int32(1) I64 = Int64(1) - R32 = Int32(1)//Int32(1) - R64 = 1//1 + R32 = Int32(1) // Int32(1) + R64 = 1 // 1 for A in Any[A1, A2, A3, A4] @test @inferred(A(F32)) === A(F32) @test @inferred(A(F64)) === A(F64) @@ -114,20 +114,20 @@ import ForwardDiff end # Nan time value: - t = 0.0:3; # Floats - u = [0, -2, -1, -2]; - A = LinearInterpolation(u, t); + t = 0.0:3 # Floats + u = [0, -2, -1, -2] + A = LinearInterpolation(u, t) dA = t -> ForwardDiff.derivative(A, t) @test isnan(dA(NaN)) - t = 0:3; # Integers - u = [0, -2, -1, -2]; - A = LinearInterpolation(u, t); + t = 0:3 # Integers + u = [0, -2, -1, -2] + A = LinearInterpolation(u, t) dA = t -> ForwardDiff.derivative(A, t) @test isnan(dA(NaN)) # Test derivative at point gives derivative to the right (except last is to left): - ts = t[begin:end-1] + ts = t[begin:(end - 1)] @test dA.(ts) == dA.(ts .+ 0.5) # Test last derivitive is to the left: @test dA(last(t)) == dA(last(t) - 0.5) @@ -135,16 +135,16 @@ import ForwardDiff # Test array-valued interpolation u = collect.(2.0collect(1:10)) t = 1.0collect(1:10) - A = LinearInterpolation(u,t) - @test A(0) == fill(0.0) - @test A(5.5) == fill(11.0) - @test A(11) == fill(22) + A = LinearInterpolation(u, t) + @test A(0) == fill(0.0) + @test A(5.5) == fill(11.0) + @test A(11) == fill(22) end @testset "Quadratic Interpolation" begin u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u,t) + A = QuadraticInterpolation(u, t) for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -158,7 +158,7 @@ end # backward-looking interpolation u = [1.0, 4.0, 9.0, 16.0] t = [1.0, 2.0, 3.0, 4.0] - A = QuadraticInterpolation(u,t,:Backward) + A = QuadraticInterpolation(u, t, :Backward) for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -172,8 +172,8 @@ end # Test both forward and backward-looking quadratic interpolation u = [1.0, 4.5, 6.0, 2.0] t = [1.0, 2.0, 3.0, 4.0] - A_f = QuadraticInterpolation(u,t,:Forward) - A_b = QuadraticInterpolation(u,t,:Backward) + A_f = QuadraticInterpolation(u, t, :Forward) + A_b = QuadraticInterpolation(u, t, :Backward) for (_t, _u) in zip(t, u) @test A_f(_t) == _u @@ -181,40 +181,40 @@ end end l₀, l₁, l₂ = 0.375, 0.75, -0.125 # In the first subinterval they're the same (no other option) - @test A_f(1.5) ≈ l₀*u[1] + l₁*u[2] + l₂*u[3] - @test A_b(1.5) ≈ l₀*u[1] + l₁*u[2] + l₂*u[3] + @test A_f(1.5) ≈ l₀ * u[1] + l₁ * u[2] + l₂ * u[3] + @test A_b(1.5) ≈ l₀ * u[1] + l₁ * u[2] + l₂ * u[3] # In the second subinterval they should be different - @test A_f(2.5) ≈ l₀*u[2] + l₁*u[3] + l₂*u[4] - @test A_b(2.5) ≈ l₂*u[1] + l₁*u[2] + l₀*u[3] + @test A_f(2.5) ≈ l₀ * u[2] + l₁ * u[3] + l₂ * u[4] + @test A_b(2.5) ≈ l₂ * u[1] + l₁ * u[2] + l₀ * u[3] # In the last subinterval they should be the same again - @test A_f(3.5) ≈ l₂*u[2] + l₁*u[3] + l₀*u[4] - @test A_b(3.5) ≈ l₂*u[2] + l₁*u[3] + l₀*u[4] + @test A_f(3.5) ≈ l₂ * u[2] + l₁ * u[3] + l₀ * u[4] + @test A_b(3.5) ≈ l₂ * u[2] + l₁ * u[3] + l₀ * u[4] # Matrix interpolation test u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0] - A = QuadraticInterpolation(u,t) + A = QuadraticInterpolation(u, t) for (_t, _u) in zip(t, eachcol(u)) @test A(_t) == _u end - @test A(0.0) == [0.0 , 0.0 ] - @test A(1.5) == [2.25 , 2.25] - @test A(2.5) == [6.25 , 6.25] + @test A(0.0) == [0.0, 0.0] + @test A(1.5) == [2.25, 2.25] + @test A(2.5) == [6.25, 6.25] @test A(3.5) == [12.25, 12.25] - @test A(5.0) == [25.0 , 25.0 ] + @test A(5.0) == [25.0, 25.0] - u_= [1.0, 4.0, 9.0, 16.0]' .* ones(5) - u = [u_[:,i] for i = 1:size(u_,2)] - A = QuadraticInterpolation(u,t) - @test A(0) == zeros(5) + u_ = [1.0, 4.0, 9.0, 16.0]' .* ones(5) + u = [u_[:, i] for i in 1:size(u_, 2)] + A = QuadraticInterpolation(u, t) + @test A(0) == zeros(5) @test A(1.5) == 2.25 * ones(5) @test A(2.5) == 6.25 * ones(5) @test A(3.5) == 12.25 * ones(5) @test A(5.0) == 25.0 * ones(5) - u = [repeat(u[i], 1, 3) for i=1:4] - A = QuadraticInterpolation(u,t) - @test A(0) == zeros(5, 3) + u = [repeat(u[i], 1, 3) for i in 1:4] + A = QuadraticInterpolation(u, t) + @test A(0) == zeros(5, 3) @test A(1.5) == 2.25 * ones(5, 3) @test A(2.5) == 6.25 * ones(5, 3) @test A(3.5) == 12.25 * ones(5, 3) @@ -224,45 +224,45 @@ end @testset "Lagrange Interpolation" begin u = [1.0, 4.0, 9.0] t = [1.0, 2.0, 3.0] - A = LagrangeInterpolation(u,t) + A = LagrangeInterpolation(u, t) @test A(2.0) == 4.0 @test A(1.5) == 2.25 u = [1.0, 8.0, 27.0, 64.0] t = [1.0, 2.0, 3.0, 4.0] - A = LagrangeInterpolation(u,t) + A = LagrangeInterpolation(u, t) @test A(2.0) == 8.0 @test A(1.5) ≈ 3.375 @test A(3.5) ≈ 42.875 u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0] - A = LagrangeInterpolation(u,t) + A = LagrangeInterpolation(u, t) - @test A(2.0) == [4.0,4.0] - @test A(1.5) ≈ [2.25,2.25] - @test A(3.5) ≈ [12.25,12.25] + @test A(2.0) == [4.0, 4.0] + @test A(1.5) ≈ [2.25, 2.25] + @test A(3.5) ≈ [12.25, 12.25] - u_= [1.0, 4.0, 9.0]' .* ones(4) - u = [u_[:,i] for i = 1:size(u_,2)] + u_ = [1.0, 4.0, 9.0]' .* ones(4) + u = [u_[:, i] for i in 1:size(u_, 2)] t = [1.0, 2.0, 3.0] - A = LagrangeInterpolation(u,t) + A = LagrangeInterpolation(u, t) @test A(2.0) == 4.0 * ones(4) @test A(1.5) == 2.25 * ones(4) - u_= [1.0, 8.0, 27.0, 64.0]' .* ones(4) - u = [u_[:,i] for i = 1:size(u_,2)] + u_ = [1.0, 8.0, 27.0, 64.0]' .* ones(4) + u = [u_[:, i] for i in 1:size(u_, 2)] t = [1.0, 2.0, 3.0, 4.0] - A = LagrangeInterpolation(u,t) + A = LagrangeInterpolation(u, t) @test A(2.0) == 8.0 * ones(4) @test A(1.5) ≈ 3.375 * ones(4) @test A(3.5) ≈ 42.875 * ones(4) - u = [repeat(u[i], 1, 3) for i=1:4] - A = LagrangeInterpolation(u,t) + u = [repeat(u[i], 1, 3) for i in 1:4] + A = LagrangeInterpolation(u, t) @test A(2.0) == 8.0 * ones(4, 3) @test A(1.5) ≈ 3.375 * ones(4, 3) @@ -290,13 +290,10 @@ end end @testset "ConstantInterpolation" begin - t = [1.0, 2.0, 3.0, 4.0] - @testset "Vector case" for u in - [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]] - - A = ConstantInterpolation(u, t, dir=:right) + @testset "Vector case" for u in [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]] + A = ConstantInterpolation(u, t, dir = :right) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -319,37 +316,37 @@ end @test A(4.5) == u[1] end - @testset "Matrix case" for u in - [[1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], ["B" "C" "A" "B"; "B" "C" "A" "B"]] - - A = ConstantInterpolation(u, t, dir=:right) - @test A(0.5) == u[:,1] - @test A(1.0) == u[:,1] - @test A(1.5) == u[:,2] - @test A(2.0) == u[:,2] - @test A(2.5) == u[:,3] - @test A(3.0) == u[:,3] - @test A(3.5) == u[:,1] - @test A(4.0) == u[:,1] - @test A(4.5) == u[:,1] + @testset "Matrix case" for u in [ + [1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], + ["B" "C" "A" "B"; "B" "C" "A" "B"], + ] + A = ConstantInterpolation(u, t, dir = :right) + @test A(0.5) == u[:, 1] + @test A(1.0) == u[:, 1] + @test A(1.5) == u[:, 2] + @test A(2.0) == u[:, 2] + @test A(2.5) == u[:, 3] + @test A(3.0) == u[:, 3] + @test A(3.5) == u[:, 1] + @test A(4.0) == u[:, 1] + @test A(4.5) == u[:, 1] A = ConstantInterpolation(u, t) # dir=:left is default - @test A(0.5) == u[:,1] - @test A(1.0) == u[:,1] - @test A(1.5) == u[:,1] - @test A(2.0) == u[:,2] - @test A(2.5) == u[:,2] - @test A(3.0) == u[:,3] - @test A(3.5) == u[:,3] - @test A(4.0) == u[:,1] - @test A(4.5) == u[:,1] + @test A(0.5) == u[:, 1] + @test A(1.0) == u[:, 1] + @test A(1.5) == u[:, 1] + @test A(2.0) == u[:, 2] + @test A(2.5) == u[:, 2] + @test A(3.0) == u[:, 3] + @test A(3.5) == u[:, 3] + @test A(4.0) == u[:, 1] + @test A(4.5) == u[:, 1] end - @testset "Vector of Vectors case" for u in - [[[1.0, 2.0], [0.0, 1.0], [1.0, 2.0], [0.0, 1.0]], + @testset "Vector of Vectors case" for u in [ + [[1.0, 2.0], [0.0, 1.0], [1.0, 2.0], [0.0, 1.0]], [["B", "C"], ["A", "B"], ["B", "C"], ["A", "B"]]] - - A = ConstantInterpolation(u, t, dir=:right) + A = ConstantInterpolation(u, t, dir = :right) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -372,11 +369,10 @@ end @test A(4.5) == u[4] end - @testset "Vector of Matrices case" for u in - [[[1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0], [1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0]], + @testset "Vector of Matrices case" for u in [ + [[1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0], [1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0]], [["B" "C"; "B" "C"], ["A" "B"; "A" "B"], ["B" "C"; "B" "C"], ["A" "B"; "A" "B"]]] - - A = ConstantInterpolation(u, t, dir=:right) + A = ConstantInterpolation(u, t, dir = :right) @test A(0.5) == u[1] @test A(1.0) == u[1] @test A(1.5) == u[2] @@ -404,42 +400,41 @@ end u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = QuadraticSpline(u,t) + A = QuadraticSpline(u, t) # Solution P₁ = x -> (x + 1)^2 # for x ∈ [-1, 0] - P₂ = x -> 2*x + 1 # for x ∈ [ 0, 1] + P₂ = x -> 2 * x + 1 # for x ∈ [ 0, 1] 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) - u_= [0.0, 1.0, 3.0]' .* ones(4) - u = [u_[:,i] for i = 1:size(u_,2)] - A = QuadraticSpline(u,t) + u_ = [0.0, 1.0, 3.0]' .* ones(4) + u = [u_[:, i] for i in 1:size(u_, 2)] + A = QuadraticSpline(u, t) @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=1:3] - A = QuadraticSpline(u,t) + u = [repeat(u[i], 1, 3) for i in 1:3] + A = QuadraticSpline(u, t) @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) end - @testset "CubicSpline Interpolation" begin u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] - A = CubicSpline(u,t) + A = CubicSpline(u, t) # Solution P₁ = x -> 1 + 1.5x + x^2 + 0.5x^3 # for x ∈ [-1.0, 0.0] @@ -455,9 +450,9 @@ end @test A(x) ≈ P₂(x) end - u_= [0.0, 1.0, 3.0]' .* ones(4) - u = [u_[:,i] for i = 1:size(u_,2)] - A = CubicSpline(u,t) + u_ = [0.0, 1.0, 3.0]' .* ones(4) + u = [u_[:, i] for i in 1:size(u_, 2)] + A = CubicSpline(u, t) for x in (-1.5, -0.5, -0.7) @test A(x) ≈ P₁(x) * ones(4) end @@ -465,8 +460,8 @@ end @test A(x) ≈ P₂(x) * ones(4) end - u = [repeat(u[i], 1, 3) for i=1:3] - A = CubicSpline(u,t) + u = [repeat(u[i], 1, 3) for i in 1:3] + A = CubicSpline(u, t) for x in (-1.5, -0.5, -0.7) @test A(x) ≈ P₁(x) * ones(4, 3) end @@ -475,41 +470,45 @@ end end end - # 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] +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 = BSplineInterpolation(u,t,2,:Uniform,:Uniform) +A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform) @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845] @test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251] @test [A(t[1]), A(t[end])] == [u[1], u[end]] -A = BSplineInterpolation(u,t,2,:ArcLen,:Average) +A = BSplineInterpolation(u, t, 2, :ArcLen, :Average) @test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609] @test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463] @test [A(t[1]), A(t[end])] == [u[1], u[end]] -A = BSplineApprox(u,t,2,4,:Uniform,:Uniform) +A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform) @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178] @test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575] @test [A(t[1]), A(t[end])] ≈ [u[1], u[end]] - # Curvefit Interpolation rng = StableRNG(12345) -model(x, p) = @. p[1]/(1 + exp(x - p[2])) -t = range(-10, stop=10, length=40) -u = model(t, [1.0, 2.0]) + 0.01*randn(rng, length(t)) +model(x, p) = @. p[1] / (1 + exp(x - p[2])) +t = range(-10, stop = 10, length = 40) +u = model(t, [1.0, 2.0]) + 0.01 * randn(rng, length(t)) p0 = [0.5, 0.5] A = Curvefit(u, t, model, p0, LBFGS()) ts = [-7.0, -2.0, 0.0, 2.5, 5.0] -vs = [1.0013468217936277, 0.9836755196317837, 0.8833959853995836, 0.3810348276782708, 0.048062978598861855] +vs = [ + 1.0013468217936277, + 0.9836755196317837, + 0.8833959853995836, + 0.3810348276782708, + 0.048062978598861855, +] us = A.(ts) @test vs ≈ us @@ -517,7 +516,7 @@ us = A.(ts) # missing values handling tests u = [1.0, 4.0, 9.0, 16.0, 25.0, missing, missing] t = [1.0, 2.0, 3.0, 4.0, missing, 6.0, missing] -A = QuadraticInterpolation(u,t) +A = QuadraticInterpolation(u, t) @test A(2.0) == 4.0 @test A(1.5) == 2.25 @@ -525,14 +524,13 @@ A = QuadraticInterpolation(u,t) @test A(2.5) == 6.25 u = copy(hcat(u, u)') -A = QuadraticInterpolation(u,t) +A = QuadraticInterpolation(u, t) @test A(2.0) == [4.0, 4.0] @test A(1.5) == [2.25, 2.25] -@test A(3.5) == [12.25,12.25] +@test A(3.5) == [12.25, 12.25] @test A(2.5) == [6.25, 6.25] - # ForwardDiff compatibility with respect to cofficients function square(INTERPOLATION_TYPE, c) # elaborate way to write f(x) = x² diff --git a/test/online_tests.jl b/test/online_tests.jl index 8b7cf949..a7da6d86 100644 --- a/test/online_tests.jl +++ b/test/online_tests.jl @@ -1,6 +1,6 @@ using DataInterpolations, Test -t = [1,2,3] +t = [1, 2, 3] u = [0, 1, 0] for di in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation] @@ -11,4 +11,4 @@ for di in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation] li = di(copy(u), copy(t)) push!(li, 1, 4) @test li == di(vcat(u, 1), vcat(t, 4)) -end \ No newline at end of file +end diff --git a/test/regularization.jl b/test/regularization.jl index 55357393..a2b540e5 100644 --- a/test/regularization.jl +++ b/test/regularization.jl @@ -4,17 +4,17 @@ using RegularizationTools # create scattered data npts = 50 xmin = 0.0 -xspan = 3/2*π -x = collect(range(xmin, xmin+xspan, length=npts)) +xspan = 3 / 2 * π +x = collect(range(xmin, xmin + xspan, length = npts)) rng = StableRNG(655) -x = x + xspan/npts*(rand(rng,npts) .- 0.5) +x = x + xspan / npts * (rand(rng, npts) .- 0.5) # select a subset randomly idx = unique(rand(rng, collect(eachindex(x)), 20)) t = x[unique(idx)] npts = length(t) ut = sin.(t) -stdev = 1e-1*maximum(ut) -u = ut + stdev*randn(rng, npts) +stdev = 1e-1 * maximum(ut) +u = ut + stdev * randn(rng, npts) # data must be ordered if t̂ is not provided idx = sortperm(t) tₒ = t[idx] @@ -24,63 +24,148 @@ tolerance = 1e-3 @testset "Direct smoothing" begin # fixed with default λ = 1.0 - A = RegularizationSmooth(uₒ,tₒ; alg=:fixed) + A = RegularizationSmooth(uₒ, tₒ; alg = :fixed) ans = [0.6456173647252937 0.663974701324226 0.7631218523665086 0.778654700697601 0.7489958320589535 0.7319087707475104 0.6807082599508811 0.6372557895089508 0.5832859790765743 0.5021274805916013 0.3065928203396211 0.1353332321156384 -0.3260000640060584 -0.6557906092739154 -0.9204882447932498]' - @test isapprox(A.û, ans, rtol=tolerance) + @test isapprox(A.û, ans, rtol = tolerance) # non-default d and λ - A = RegularizationSmooth(uₒ,tₒ, 4; λ=1e-2, alg=:fixed) + A = RegularizationSmooth(uₒ, tₒ, 4; λ = 1e-2, alg = :fixed) ans = [0.19865190868740357 0.2885349151737291 0.6756699442978945 0.9165887141895426 0.9936113717653254 1.0042825002191034 0.9768118192829827 0.9184595331808411 0.8214983284892922 0.6538356458824783 0.28295521578898 0.018060767871253963 -0.5301723647977373 -0.8349855890541111 -1.1085048455468356]' - @test isapprox(A.û, ans, rtol=tolerance) + @test isapprox(A.û, ans, rtol = tolerance) # GCV (default) to determine λ - A = RegularizationSmooth(uₒ,tₒ) - @test isapprox(A.λ, 0.12788440382063268, rtol=tolerance) + A = RegularizationSmooth(uₒ, tₒ) + @test isapprox(A.λ, 0.12788440382063268, rtol = tolerance) ans = [0.21974931848164914 0.2973284508009968 0.6908546278415386 0.9300465474303226 0.9741453042418977 0.9767572556868123 0.9432951659303452 0.8889834700087442 0.804842790047182 0.6603217445567791 0.30341652659101737 0.05924456463634589 -0.5239939779242144 -0.8421768233191822 -1.107517099580091]' - @test isapprox(A.û, ans, rtol=tolerance) + @test isapprox(A.û, ans, rtol = tolerance) # L-curve to determine λ - A = RegularizationSmooth(uₒ,tₒ; alg=:L_curve) - @test isapprox(A.λ, 0.9536286111306728, rtol=tolerance) - ans = [0.6261657429321232, 0.6470204841904836, 0.7599270022828396, 0.7835725197598925, 0.7567105872094757, 0.7404815750685363, 0.6906841961987067, 0.647628105931872, 0.5937273796308717, 0.512087780658067, 0.3136272387739983, 0.1392761732695201, -0.3312498167413961, -0.6673268474631847, -0.9370342562716745] - @test isapprox(A.û, ans, rtol=tolerance) + A = RegularizationSmooth(uₒ, tₒ; alg = :L_curve) + @test isapprox(A.λ, 0.9536286111306728, rtol = tolerance) + ans = [ + 0.6261657429321232, + 0.6470204841904836, + 0.7599270022828396, + 0.7835725197598925, + 0.7567105872094757, + 0.7404815750685363, + 0.6906841961987067, + 0.647628105931872, + 0.5937273796308717, + 0.512087780658067, + 0.3136272387739983, + 0.1392761732695201, + -0.3312498167413961, + -0.6673268474631847, + -0.9370342562716745, + ] + @test isapprox(A.û, ans, rtol = tolerance) end @testset "Smoothing with weights" begin # midpoint rule integration - A = RegularizationSmooth(uₒ,tₒ, nothing, :midpoint) - @test isapprox(A.λ, 0.10787235405005478, rtol=tolerance) - ans = [0.3068904607028622, 0.3637388879266782, 0.6654462500501238, 0.9056440536733456, 0.9738150157541853, 0.9821315604309402, 0.9502526946446999, 0.8953643918063283, 0.8024431779821514, 0.6415812230114304, 0.2834706832220367, 0.05281575111822609, -0.5333542714497277, -0.8406745098604134, -1.0983391396173634] - @test isapprox(A.û, ans, rtol=tolerance) + A = RegularizationSmooth(uₒ, tₒ, nothing, :midpoint) + @test isapprox(A.λ, 0.10787235405005478, rtol = tolerance) + ans = [ + 0.3068904607028622, + 0.3637388879266782, + 0.6654462500501238, + 0.9056440536733456, + 0.9738150157541853, + 0.9821315604309402, + 0.9502526946446999, + 0.8953643918063283, + 0.8024431779821514, + 0.6415812230114304, + 0.2834706832220367, + 0.05281575111822609, + -0.5333542714497277, + -0.8406745098604134, + -1.0983391396173634, + ] + @test isapprox(A.û, ans, rtol = tolerance) # arbitrary weights for wls (and fixed λ, GCV not working well for some of these) - A = RegularizationSmooth(uₒ,tₒ, nothing, collect(1:npts); λ=1e-1, alg=:fixed) - ans = [0.24640196218427968, 0.3212059975226125, 0.6557626475144205, 0.9222911426465459, 0.9913331910731215, 1.0072241662103494, 0.9757899817730779, 0.935880516370941, 0.8381074902073471, 0.6475589703422522, 0.2094170714475404, 0.09102085384961625, -0.5640882848240228, -0.810519277110118, -1.1159124134900906] - @test isapprox(A.û, ans, rtol=tolerance) + A = RegularizationSmooth(uₒ, tₒ, nothing, collect(1:npts); λ = 1e-1, alg = :fixed) + ans = [ + 0.24640196218427968, + 0.3212059975226125, + 0.6557626475144205, + 0.9222911426465459, + 0.9913331910731215, + 1.0072241662103494, + 0.9757899817730779, + 0.935880516370941, + 0.8381074902073471, + 0.6475589703422522, + 0.2094170714475404, + 0.09102085384961625, + -0.5640882848240228, + -0.810519277110118, + -1.1159124134900906, + ] + @test isapprox(A.û, ans, rtol = tolerance) # arbitrary weights for wls and wr - nhalf = Int(floor(npts/2)) - wls = vcat(ones(nhalf),10*ones(npts-nhalf)) - wr = collect(1:npts-2) - A = RegularizationSmooth(uₒ,tₒ, nothing, wls, wr; λ=1e-1, alg=:fixed) - ans = [0.21878709713242372, 0.3118480645325099, 0.7669822464946172, 1.0232343854914931, 1.0526513115274412, 1.0469579284244412, 0.9962426294084775, 0.9254407155702626, 0.8204764044515936, 0.6514510142804217, 0.27796896299068763, 0.04756024028728636, -0.5301034620974782, -0.8408107101140526, -1.1058428573417736] - @test isapprox(A.û, ans, rtol=tolerance) + nhalf = Int(floor(npts / 2)) + wls = vcat(ones(nhalf), 10 * ones(npts - nhalf)) + wr = collect(1:(npts - 2)) + A = RegularizationSmooth(uₒ, tₒ, nothing, wls, wr; λ = 1e-1, alg = :fixed) + ans = [ + 0.21878709713242372, + 0.3118480645325099, + 0.7669822464946172, + 1.0232343854914931, + 1.0526513115274412, + 1.0469579284244412, + 0.9962426294084775, + 0.9254407155702626, + 0.8204764044515936, + 0.6514510142804217, + 0.27796896299068763, + 0.04756024028728636, + -0.5301034620974782, + -0.8408107101140526, + -1.1058428573417736, + ] + @test isapprox(A.û, ans, rtol = tolerance) end @testset "Smoothing with t̂ provided" begin N̂ = 20 - t̂ = collect(range(xmin, xmin+xspan, length=N̂)) + t̂ = collect(range(xmin, xmin + xspan, length = N̂)) # with t̂, no weights - A = RegularizationSmooth(u,t, t̂) - @test isapprox(A.λ, 0.138273889585313, rtol=tolerance) + A = RegularizationSmooth(u, t, t̂) + @test isapprox(A.λ, 0.138273889585313, rtol = tolerance) ans = [0.21626377852882872 0.39235926952322575 0.5573848799950002 0.7072474496656729 0.8361906119247042 0.9313473799797176 0.9809844353757837 0.9750833208625507 0.9096038940899813 0.7816929736202427 0.6052694276527628 0.4015903497629387 0.1913719025253403 -0.01979786871512895 -0.23400354001942947 -0.44481229967011127 -0.6457913359497256 -0.8405146928672158 -1.0367229293434395 -1.2334090099343238]' - @test isapprox(A.û, ans, rtol=tolerance) + @test isapprox(A.û, ans, rtol = tolerance) # t̂ and wls - A = RegularizationSmooth(u,t, t̂, collect(1:npts)) - @test isapprox(A.λ, 0.26746430253489195, rtol=tolerance) + A = RegularizationSmooth(u, t, t̂, collect(1:npts)) + @test isapprox(A.λ, 0.26746430253489195, rtol = tolerance) ans = [0.3118247878815087 0.44275860852897864 0.5705834985506882 0.6979119448253899 0.8234189540704866 0.9289458273102476 0.9970803470992273 1.0071205506077525 0.9443157518324818 0.7954860908242515 0.5847385548859145 0.34813493129868633 0.1237494751337505 -0.0823517516424196 -0.28265170846635246 -0.4760833187699964 -0.6615795059853024 -0.844779821396189 -1.0341162283806349 -1.225270266213379]' - @test isapprox(A.û, ans, rtol=tolerance) + @test isapprox(A.û, ans, rtol = tolerance) # t̂, wls, and wr - nhalf = Int(floor(npts/2)) - wls = vcat(ones(nhalf),10*ones(npts-nhalf)) - wr = collect(1:N̂-2) - A = RegularizationSmooth(u,t, t̂, wls, wr) - @test isapprox(A.λ, 0.04555080890920959, rtol=tolerance) - ans = [0.2799800686914433, 0.4627548444527547, 0.5611922868318674, 0.6647761469309206, 0.7910803348948329, 0.9096001134420562, 1.0067644979677808, 1.0541868144785513, 0.9889720466386331, 0.8088479651943575, 0.5677592185997403, 0.31309698432269184, 0.08587106716465115, -0.11476265128730469, -0.30749376694236485, -0.4942769809562725, -0.676806367664006, -0.8587832527770329, -1.0443430843364814, -1.2309001260104093] - @test isapprox(A.û, ans, rtol=tolerance) + nhalf = Int(floor(npts / 2)) + wls = vcat(ones(nhalf), 10 * ones(npts - nhalf)) + wr = collect(1:(N̂ - 2)) + A = RegularizationSmooth(u, t, t̂, wls, wr) + @test isapprox(A.λ, 0.04555080890920959, rtol = tolerance) + ans = [ + 0.2799800686914433, + 0.4627548444527547, + 0.5611922868318674, + 0.6647761469309206, + 0.7910803348948329, + 0.9096001134420562, + 1.0067644979677808, + 1.0541868144785513, + 0.9889720466386331, + 0.8088479651943575, + 0.5677592185997403, + 0.31309698432269184, + 0.08587106716465115, + -0.11476265128730469, + -0.30749376694236485, + -0.4942769809562725, + -0.676806367664006, + -0.8587832527770329, + -1.0443430843364814, + -1.2309001260104093, + ] + @test isapprox(A.û, ans, rtol = tolerance) end diff --git a/test/runtests.jl b/test/runtests.jl index a63c1ae8..0c971d40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,20 @@ using DataInterpolations, Test -@testset "Interface" begin include("interface.jl") end -@testset "Interpolation Tests" begin include("interpolation_tests.jl") end -@testset "Derivative Tests" begin include("derivative_tests.jl") end -@testset "Integral Tests" begin include("integral_tests.jl") end -@testset "Online Tests" begin include("online_tests.jl") end -@testset "Regularization Smoothing" begin include("regularization.jl") end +@testset "Interface" begin + include("interface.jl") +end +@testset "Interpolation Tests" begin + include("interpolation_tests.jl") +end +@testset "Derivative Tests" begin + include("derivative_tests.jl") +end +@testset "Integral Tests" begin + include("integral_tests.jl") +end +@testset "Online Tests" begin + include("online_tests.jl") +end +@testset "Regularization Smoothing" begin + include("regularization.jl") +end