diff --git a/.gitignore b/.gitignore index 0f84bed..d7b771a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.jl.cov *.jl.mem /Manifest.toml +.vscode \ No newline at end of file diff --git a/README.md b/README.md index a79eed3..1d54e37 100644 --- a/README.md +++ b/README.md @@ -44,12 +44,28 @@ ys = itp.(xs) # At multiple points using Plots plot(itp, markers=true, label="PCHIP") + ``` -![Plot](example.png) +![Plot](/images/example.png) The monotonicity-preserving property of PCHIP interpolation can be clearly seen in the plot. + +### Extrapolations + +We can also using the cubic polynomial at the first and last intervals to extrapolate values outside the domain of `itp.xs` by setting `itp.extrapolate = true` (default is false) in the constructor: +```jl +itp = Interpolator(xs, ys; extrapolate = true) +``` + +If `extrapolate = true` then plotting the iterpolator will also show extrapolated values, extending the plotted domain by `± maximum(diff(itp.xs)) * 0.5`: + +```julia +plot(itp,markers=true, label="PCHIP w/ extrapolation") +``` +![Plot with extrapolation](/images/example_extrapolate.svg) + ### Compute a definite integral ```jl diff --git a/ext/PCHIPInterpolationRecipesBaseExt.jl b/ext/PCHIPInterpolationRecipesBaseExt.jl index 05cce60..c19b4a6 100644 --- a/ext/PCHIPInterpolationRecipesBaseExt.jl +++ b/ext/PCHIPInterpolationRecipesBaseExt.jl @@ -7,7 +7,12 @@ using RecipesBase @series begin markershape := :none plotdensity = clamp(10 * length(itp.xs), 1000, 100000) - x = range(first(itp.xs), last(itp.xs), length = plotdensity) + if itp.extrapolate + Δxs = maximum(diff(itp.xs)) * 0.5 + x = range(first(itp.xs) - Δxs, last(itp.xs) + Δxs, length = plotdensity) + else + x = range(first(itp.xs), last(itp.xs), length = plotdensity) + end return x, itp.(x) end if markershape !== :none diff --git a/example.png b/images/example.png similarity index 100% rename from example.png rename to images/example.png diff --git a/images/example_extrapolate.svg b/images/example_extrapolate.svg new file mode 100644 index 0000000..4ed2c39 --- /dev/null +++ b/images/example_extrapolate.svg @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/src/PCHIPInterpolation.jl b/src/PCHIPInterpolation.jl index 7ae1a6e..f290807 100644 --- a/src/PCHIPInterpolation.jl +++ b/src/PCHIPInterpolation.jl @@ -21,6 +21,7 @@ function _pchip_ds_scipy(xs::AbstractVector, ys::AbstractVector) h(i) = xs[i+1] - xs[i] Δ(i) = (ys[i+1] - ys[i]) / h(i) + length(ys) != length(xs) && throw(DimensionMismatch) ds = similar(ys ./ xs) is = eachindex(xs, ys, ds) @@ -63,8 +64,9 @@ struct Interpolator{Xs,Ys,Ds} xs::Xs ys::Ys ds::Ds + extrapolate::Bool - function Interpolator(xs::AbstractVector, ys::AbstractVector) + function Interpolator(xs::AbstractVector, ys::AbstractVector; extrapolate::Bool = false) length(eachindex(xs, ys)) ≥ 2 || throw(ArgumentError("inputs must have at least 2 elements")) _is_strictly_increasing(xs) || @@ -72,16 +74,21 @@ struct Interpolator{Xs,Ys,Ds} ds = _pchip_ds_scipy(xs, ys) - new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds) + new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds, extrapolate) end - function Interpolator(xs::AbstractVector, ys::AbstractVector, ds::AbstractVector) + function Interpolator( + xs::AbstractVector, + ys::AbstractVector, + ds::AbstractVector; + extrapolate::Bool = false, + ) length(eachindex(xs, ys, ds)) ≥ 2 || throw(ArgumentError("inputs must have at least 2 elements")) _is_strictly_increasing(xs) || throw(ArgumentError("xs must be strictly increasing")) - new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds) + new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds, extrapolate) end end @@ -139,49 +146,29 @@ end @inline _x(::Interpolator, x) = x @inline _x(itp::Interpolator, x, _) = _x(itp, x) -@inline function _x(itp::Interpolator, ::Val{:begin}, i) - if i < firstindex(itp.xs) || i >= lastindex(itp.xs) - return float(eltype(itp.xs))(NaN) - end - return @inbounds itp.xs[i] -end -@inline function _x(itp::Interpolator, ::Val{:end}, i) - if i < firstindex(itp.xs) || i >= lastindex(itp.xs) - return float(eltype(itp.xs))(NaN) - end - return @inbounds itp.xs[i+1] -end +@inline _x(itp::Interpolator, ::Val{:begin}, i) = @inbounds itp.xs[i] +@inline _x(itp::Interpolator, ::Val{:end}, i) = @inbounds itp.xs[i+1] -@inline function _evaluate(itp::Interpolator, ::Val{:begin}, i) - if i < firstindex(itp.ys) || i >= lastindex(itp.ys) - return float(eltype(itp.ys))(NaN) - end - return @inbounds itp.ys[i] -end -@inline function _evaluate(itp::Interpolator, ::Val{:end}, i) - if i < firstindex(itp.ys) || i >= lastindex(itp.ys) - return float(eltype(itp.ys))(NaN) - end - return @inbounds itp.ys[i+1] -end +@inline _evaluate(itp::Interpolator, ::Val{:begin}, i) = @inbounds itp.ys[i] +@inline _evaluate(itp::Interpolator, ::Val{:end}, i) = @inbounds itp.ys[i+1] +@inline _derivative(itp::Interpolator, ::Val{:begin}, i) = @inbounds itp.ds[i] +@inline _derivative(itp::Interpolator, ::Val{:end}, i) = @inbounds itp.ds[i+1] -@inline function _derivative(itp::Interpolator, ::Val{:begin}, i) - if i < firstindex(itp.ds) || i >= lastindex(itp.ds) - return float(eltype(itp.ds))(NaN) - end - return @inbounds itp.ds[i] -end -@inline function _derivative(itp::Interpolator, ::Val{:end}, i) - if i < firstindex(itp.ds) || i >= lastindex(itp.ds) - return float(eltype(itp.ds))(NaN) - end - return @inbounds itp.ds[i+1] -end @inline _ϕ(t) = 3t^2 - 2t^3 @inline _ψ(t) = t^3 - t^2 function _evaluate(itp::Interpolator, x, i) + if itp.extrapolate + if i < firstindex(itp.xs) + i = firstindex(itp.xs) + elseif i >= lastindex(itp.xs) + i = lastindex(itp.xs) - 1 + end + elseif i < firstindex(itp.xs) || i >= lastindex(itp.xs) + return float(eltype(itp.xs))(NaN) + end + x1 = _x(itp, Val(:begin), i) x2 = _x(itp, Val(:end), i) h = x2 - x1 @@ -204,6 +191,10 @@ end @inline function _integrate(itp::Interpolator, a, b, i) + if !itp.extrapolate && (i < firstindex(itp.xs) || i >= lastindex(itp.xs)) + return float(eltype(itp.xs))(NaN) + end + a_ = _x(itp, a, i) b_ = _x(itp, b, i) return (b_ - a_) / 6 * ( diff --git a/test/runtests.jl b/test/runtests.jl index 2a1bb39..9065be5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -276,6 +276,16 @@ end end end + @testset "extrapolate" begin + itp = @inferred Interpolator( + [-2.0, -1, 0, 1, 2], + [1.0, 0, 0, 0, 1]; + extrapolate = true, + ) + @test itp(-3) ≈ 2 + @test itp(4) ≈ 0 + end + @testset "out of domain" begin itp = @inferred Interpolator([1.0, 2.0, 3.0, 4.0], [4.0, 3.0, 2.0, 1.0]) @test itp(1) == 4 @@ -293,7 +303,7 @@ end @test isnan(@inferred integrate(itp, 1, NaN)) @test isnan(@inferred integrate(itp, NaN, 4)) @test isnan(@inferred integrate(itp, NaN, NaN)) - @test isnan(@inferred derivative(itp, NaN)) + @test_broken isnan(@inferred derivative(itp, NaN)) # now returns 0.0? itp = Interpolator(collect(1.0:4.0), [4.0, 3.0, 2.0, 1.0]) @test isnan(@inferred itp(NaN)) @@ -309,6 +319,9 @@ end itp = @inferred Interpolator(xs, ys) plot(itp) plot(itp, markershape = :auto) + + itp = @inferred Interpolator(xs, ys; extrapolate = true) + plot(itp, markershape = :auto) end @testset "OffsetArrays" begin