From 47f6c859c441b45381e73630b2468c090791561a Mon Sep 17 00:00:00 2001 From: Aaron Date: Sun, 24 Aug 2025 16:06:08 +1000 Subject: [PATCH] Started `DimensionalDataInterpolationsExt.jl` for `AbstractDimVector`s (#420) --- Project.toml | 6 +- docs/src/.vitepress/config.mts | 1 + docs/src/interpolations.md | 453 +++++++++++ ext/DimensionalDataDataInterpolationsExt.jl | 180 +++++ src/interpolations.jl | 3 + test/interpolation.jl | 822 ++++++++++++++++++++ test/runtests.jl | 2 +- 7 files changed, 1465 insertions(+), 2 deletions(-) create mode 100644 docs/src/interpolations.md create mode 100644 ext/DimensionalDataDataInterpolationsExt.jl create mode 100644 src/interpolations.jl create mode 100644 test/interpolation.jl diff --git a/Project.toml b/Project.toml index 595268942..900cf6f9e 100644 --- a/Project.toml +++ b/Project.toml @@ -28,6 +28,7 @@ AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67" CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" @@ -41,6 +42,7 @@ DimensionalDataAlgebraOfGraphicsExt = "AlgebraOfGraphics" DimensionalDataCategoricalArraysExt = "CategoricalArrays" DimensionalDataChainRulesCoreExt = "ChainRulesCore" DimensionalDataDiskArraysExt = "DiskArrays" +DimensionalDataDataInterpolationsExt = "DataInterpolations" DimensionalDataMakie = "Makie" DimensionalDataNearestNeighborsExt = "NearestNeighbors" DimensionalDataPythonCall = "PythonCall" @@ -63,6 +65,7 @@ ConstructionBase = "1" CoordinateTransformations = "0.6" DataAPI = "1.16" DataFrames = "1" +DataInterpolations = "8.5.0" Dates = "1" DiskArrays = "0.3, 0.4" Distributions = "0.25" @@ -111,6 +114,7 @@ ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" @@ -133,4 +137,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["AlgebraOfGraphics", "Aqua", "ArrayInterface", "BenchmarkTools", "CairoMakie", "CategoricalArrays", "ChainRulesCore", "ColorTypes", "Combinatorics", "CoordinateTransformations", "DataFrames", "DiskArrays", "Distributions", "Documenter", "FFTW", "GPUArrays", "ImageFiltering", "ImageTransformations", "JLArrays", "NearestNeighbors", "OffsetArrays", "OpenSSL", "Plots", "PythonCall", "Random", "SafeTestsets", "SparseArrays", "StatsBase", "StatsPlots", "Test", "Unitful"] +test = ["AlgebraOfGraphics", "Aqua", "ArrayInterface", "BenchmarkTools", "CairoMakie", "CategoricalArrays", "ChainRulesCore", "ColorTypes", "Combinatorics", "CoordinateTransformations", "DataFrames", "DataInterpolations", "DiskArrays", "Distributions", "Documenter", "FFTW", "GPUArrays", "ImageFiltering", "ImageTransformations", "JLArrays", "NearestNeighbors", "OffsetArrays", "OpenSSL", "Plots", "PythonCall", "Random", "SafeTestsets", "SparseArrays", "StatsBase", "StatsPlots", "Test", "Unitful"] diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index f0503bf04..c10c387bf 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -30,6 +30,7 @@ const navTemp = { { text: 'CUDA and GPUs', link: '/cuda' }, { text: 'DiskArrays', link: '/diskarrays' }, { text: 'Xarray', link: '/xarray' }, + { text: 'DataInterpolations', link: '/interpolation' }, { text: 'Extending DimensionalData', link: '/extending_dd' }, ], }, diff --git a/docs/src/interpolations.md b/docs/src/interpolations.md new file mode 100644 index 000000000..2a290974c --- /dev/null +++ b/docs/src/interpolations.md @@ -0,0 +1,453 @@ +```@meta +Description = "Interpolate and extrapolate `DimArray` data - interoperability with `DataInterpolations.jl` and `DataInterpolationsND.jl`" +``` + +# Interpolation of `DimensionalData.jl` Objects + +The following functionalities are available +for DimensionalData version 0.29.25 +with DataInterpolations version +but are experimental and under development. +Breaking changes may come in future patches. + +## 1D Interpolation with `DataInterpolations.jl` + +Convenience interpolation methods of `DimensionalData` objects +has been implemented as the package extension +`DimensionalDataDataInterpolationsExt`, +which is loaded simply by loading the two dependent packages: + +```@ansi Interpolation1D +using DimensionalData +using DataInterpolations +``` + +All interpolation methods provided by `DataInterpolations` +are supported. Explicitly: + +* `LinearInterpolation` +* `QuadraticInterpolation` +* `LagrangeInterpolation` +* `AkimaInterpolation` +* `ConstantInterpolation` +* `SmoothedConstantInterpolation` +* `QuadraticSpline` +* `CubicSpline` +* `BSplineInterpolation` +* `CubicHermiteSpline` +* `PCHIPInterpolation` +* `QuinticHermiteSpline` + +All extrapolation types provided by `DataInterpolations` +are also supported. Explicitly: + +* `ExtrapolationType.None` +* `ExtrapolationType.Constant` +* `ExtrapolationType.Linear` +* `ExtrapolationType.Extension` +* `ExtrapolationType.Periodic` +* `ExtrapolationType.Reflective` + +The extension methods are defined such that +the two positional arguments of `u` and `t` are replaced with +the single positional argument of an `AbstractDimVector`. +The `parent` vector is defined as `u` +and the `Dimension` vector is defined as `t`. + +For the following examples, the following `DimArray` is used: + +```@ansi Interpolation1D +A = DimArray( + [14.7, 11.51, 10.41, 14.95, 12.24, 11.22], + [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] |> X +) +dA = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0011] +ddA = [0.0, -0.00033, 0.0051, -0.0067, 0.0029, 0.0] + +nothing # hide +``` + +### Linear Interpolation Examples + +::: tabs + + +== Interpolation + +```@ansi Interpolation1D +itp = LinearInterpolation(A) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(0 : 10.0 : 100) +``` + + +== Extrapolation + +```@ansi Interpolation1D +itp = LinearInterpolation( + A; + extrapolation = ExtrapolationType.Periodic +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + + +== Left Extrapolation + +```@ansi Interpolation1D +itp = LinearInterpolation( + A; + extrapolation_left = ExtrapolationType.Extension +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 250) +``` + + +== Right Extrapolation + +```@ansi Interpolation1D +itp = LinearInterpolation( + A; + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(50 : 100.0 : 250) +``` + + +== Mixed Extrapolation + +```@ansi Interpolation1D +itp = LinearInterpolation( + A; + extrapolation_left = ExtrapolationType.Constant, + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + +### Quadratic Interpolation Examples + +::: tabs + + +== Interpolation + +```@ansi Interpolation1D +itp = QuadraticInterpolation(A) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(0 : 10.0 : 100) +``` + + +== Backward + +```@ansi Interpolation1D +itp = QuadraticInterpolation(A, :Backward) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(0 : 10.0 : 100) +``` + + +== Extrapolation + +```@ansi Interpolation1D +itp = QuadraticInterpolation( + A; + extrapolation = ExtrapolationType.Periodic +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + + +== Left Extrapolation + +```@ansi Interpolation1D +itp = QuadraticInterpolation( + A, :Backward; + extrapolation_left = ExtrapolationType.Extension +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 250) +``` + + +== Right Extrapolation + +```@ansi Interpolation1D +itp = QuadraticInterpolation( + A; + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(50 : 100.0 : 250) +``` + + +== Mixed Extrapolation + +```@ansi Interpolation1D +itp = QuadraticInterpolation( + A, :Backward; + extrapolation_left = ExtrapolationType.Constant, + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + +### Quintic Hermite Spline Examples + +::: tabs + + +== Interpolation + +```@ansi Interpolation1D +itp = QuinticHermiteSpline(ddA, dA, A) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(0 : 10.0 : 100) +``` + + +== Extrapolation + +```@ansi Interpolation1D +itp = QuinticHermiteSpline( + ddA, dA, A; + extrapolation = ExtrapolationType.Periodic +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + + +== Left Extrapolation + +```@ansi Interpolation1D +itp = QuinticHermiteSpline( + ddA, dA, A; + extrapolation_left = ExtrapolationType.Extension +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 250) +``` + + +== Right Extrapolation + +```@ansi Interpolation1D +itp = QuinticHermiteSpline( + ddA, dA, A; + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(50 : 100.0 : 250) +``` + + +== Mixed Extrapolation + +```@ansi Interpolation1D +itp = QuinticHermiteSpline( + ddA, dA, A; + extrapolation_left = ExtrapolationType.Constant, + extrapolation_right = ExtrapolationType.Reflective +) +``` + +```@ansi Interpolation1D +itp(12.34) +``` + +```@ansi Interpolation1D +itp(-12.34) +``` + +```@ansi Interpolation1D +itp(1234.56) +``` + +```@ansi Interpolation1D +itp(-250 : 100.0 : 450) +``` + +## ND Interpolation with `DataInterpolationsND.jl` + +Not yet implemented. + +## Roadmap and Discussion Prompts + +Feel free to follow, contribute to, and veto features in +the discussion on this functionality at +[Issue 420](https://github.com/rafaqz/DimensionalData.jl/issues/420), +especially regarding the following planned features. + +For `itp::AbstractInterpolation` +(including `itp = PCHIPInterpolation(...)`): + +* Enforce `<:Number` on array and dimension element types via stricter method signatures. +* `itp(::Dimension)` and `itp(::DimArray)` return a `DimArray` instead of a raw `Array`. +* `itp` construction requiring/accepting additional numerical information (such as derivatives at data points), define signatures for specification of the derivatives in a `DimStack` and `DimTree`. +* Expansion of functionality of 1D interpolation to receive `DimArray`s and specify the `Dimension` for the interpolation. +* Check downstream functionalities of `DimensionalData` extensions such as `Rasters.jl`. +* Package extension for `DataInterpolationsND.jl`. +* Package extension for `Interpolations.jl`. diff --git a/ext/DimensionalDataDataInterpolationsExt.jl b/ext/DimensionalDataDataInterpolationsExt.jl new file mode 100644 index 000000000..11bade642 --- /dev/null +++ b/ext/DimensionalDataDataInterpolationsExt.jl @@ -0,0 +1,180 @@ +module DimensionalDataDataInterpolationsExt + +using DimensionalData +using DataInterpolations + +function (Itp::Type{<:DataInterpolations.AbstractInterpolation})( + data::AbstractDimVector, + args...; + kw... +) + return Itp( + data |> parent, + dims(data) |> only |> parent |> parent, + args...; + kw... + ) +end + +doctargets = [ + :LinearInterpolation, + :QuadraticInterpolation, + :LagrangeInterpolation, + :AkimaInterpolation, + :ConstantInterpolation, + :QuadraticSpline, + :CubicSpline, + :BSplineInterpolation, + :PCHIPInterpolation +] + +for doctarget in doctargets + @eval begin + """ + ``` + $($doctarget)( + ut::DimensionalData.AbstractDimVector, + ...; ... + ) + ``` + + The two positional arguments `u, t` + in the original method definitions + are assigned the parent array of `ut` + and the `Dimension` of `ut` respectively. + + Remaining positional and keyword arguments are as per the original + `DataInterpolations.$($doctarget)` method definitions. + + This interoperability between + DimensionalData.jl and DataInterpolations.jl + is experimental and under development. + """ + $(doctarget) + end +end + +function DataInterpolations.QuadraticInterpolation( + data::AbstractDimVector, + mode::Symbol, + args...; + kw... +) + return QuadraticInterpolation( + data |> parent, + dims(data) |> only |> parent |> parent, + mode, + args...; + kw... + ) +end + +""" +``` +PCHIPInterpolation( + ut::DimensionalData.AbstractDimVector, + ...; ... +) +``` + +The two positional arguments `u, t` +in the original method definitions +are assigned the parent array of `ut` +and the `Dimension` of `ut` respectively. + +Remaining positional and keyword arguments are as per the original +`DataInterpolationsQuadraticInterpolation` method definitions. + +This interoperability between +DimensionalData.jl and DataInterpolations.jl +is experimental and under development. +""" +function DataInterpolations.PCHIPInterpolation( + data::AbstractDimVector, + args...; + kw... +) + return PCHIPInterpolation( + data |> parent, + dims(data) |> only |> parent |> parent, + args...; + kw... + ) +end + +""" +``` +CubicHermiteSpline( + du::AbstractVector, + ut::DimensionalData.AbstractDimVector, + ...; ... +) +``` + +The two positional arguments `u, t` +in the original method definitions +are assigned the parent array of `ut` +and the `Dimension` of `ut` respectively. + +Remaining positional and keyword arguments are as per the original +`DataInterpolationsQuadraticInterpolation` method definitions. + +This interoperability between +DimensionalData.jl and DataInterpolations.jl +is experimental and under development. +""" +function DataInterpolations.CubicHermiteSpline( + du::AbstractVector, + data::AbstractDimVector, + args...; + kw... +) + return CubicHermiteSpline( + du, + data |> parent, + dims(data) |> only |> parent |> parent, + args...; + kw... + ) +end + +""" +``` +QuinticHermiteSpline( + ddu::AbstractVector, + du::AbstractVector, + ut::DimensionalData.AbstractDimVector, + ...; ... +) +``` + +The two positional arguments `u, t` +in the original method definitions +are assigned the parent array of `ut` +and the `Dimension` of `ut` respectively. + +Remaining positional and keyword arguments are as per the original +`DataInterpolationsQuadraticInterpolation` method definitions. + +This interoperability between +DimensionalData.jl and DataInterpolations.jl +is experimental and under development. +""" +function DataInterpolations.QuinticHermiteSpline( + ddu::AbstractVector, + du::AbstractVector, + data::AbstractDimVector, + args...; + kw... +) + return QuinticHermiteSpline( + ddu, + du, + data |> parent, + dims(data) |> only |> parent |> parent, + args...; + kw... + ) +end + +end diff --git a/src/interpolations.jl b/src/interpolations.jl new file mode 100644 index 000000000..e6dfb195a --- /dev/null +++ b/src/interpolations.jl @@ -0,0 +1,3 @@ +function define_interpolation(Itp, data::AbstractDimVector) + +end diff --git a/test/interpolation.jl b/test/interpolation.jl new file mode 100644 index 000000000..fec75743d --- /dev/null +++ b/test/interpolation.jl @@ -0,0 +1,822 @@ +using DimensionalData +using DataInterpolations +using Test + +x_vec = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3] + +x_vec_dense = range(extrema(x_vec)..., 123) + +x_vec_wider = [ + minimum(x_vec) - 100 : 50 : maximum(x_vec) + 100 + extrema(x_vec)... +] |> unique! |> sort! + +x_vec_left = filter(≤(x_vec |> maximum), x_vec_wider) + +x_vec_right = filter(≥(x_vec |> minimum), x_vec_wider) + +A_vec = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22] + +dA_vec = [-0.047, -0.058, 0.054, 0.012, -0.068, 0.0011] + +ddA_vec = [0.0, -0.00033, 0.0051, -0.0067, 0.0029, 0.0] + +A_dimvec = DimArray(A_vec, X(x_vec)) + +reproducers = [ + LinearInterpolation + QuadraticInterpolation + LagrangeInterpolation + AkimaInterpolation + ConstantInterpolation + QuadraticSpline + CubicSpline + BSplineInterpolation + CubicHermiteSpline + PCHIPInterpolation + QuinticHermiteSpline +] + +interpolators = [ + reproducers... + SmoothedConstantInterpolation +] + +extrapolators = [ + ExtrapolationType.Constant + ExtrapolationType.Linear + ExtrapolationType.Extension + ExtrapolationType.Periodic + ExtrapolationType.Reflective +] + +@testset "Reproduces Data" begin + @testset "`$Itp`" for Itp in reproducers + itp = if Itp == BSplineInterpolation + BSplineInterpolation(A_dimvec, 3, :ArcLen, :Average) + elseif Itp == CubicHermiteSpline + CubicHermiteSpline(dA_vec, A_dimvec) + elseif Itp == QuinticHermiteSpline + QuinticHermiteSpline(ddA_vec, dA_vec, A_dimvec) + else + Itp(A_dimvec) + end + + if Itp in [CubicSpline, QuinticHermiteSpline, BSplineInterpolation] + @testset "$x" for x in x_vec + @test isapprox( + A_dimvec[X = At(x)], + itp(x); + atol = 1e-14 + ) + end + + @test isapprox( + A_dimvec |> parent, + itp(x_vec) + ) + else + @testset "$x" for x in x_vec + @test isequal( + A_dimvec[X = At(x)], + itp(x) + ) + end + + @test isequal( + A_dimvec |> parent, + itp(x_vec) + ) + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp = QuadraticInterpolation(A_dimvec, :Backward) + @testset "$x" for x in x_vec + @test isequal( + A_dimvec[X = At(x)], + itp(x) + ) + @test isequal( + A_vec, + itp(x_vec) + ) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp = ConstantInterpolation(A_dimvec; dir = :Right) + @testset "$x" for x in x_vec + @test isequal( + A_dimvec[X = At(x)], + itp(x) + ) + @test isequal( + A_vec, + itp(x_vec) + ) + end + end +end + +@testset "Preserves Interpolation" begin + @testset "$Itp" for Itp in interpolators + (; itp_raw, itp_dim) = if Itp == BSplineInterpolation + itp_raw = BSplineInterpolation(A_vec, x_vec, 3, :ArcLen, :Average) + itp_dim = BSplineInterpolation(A_dimvec, 3, :ArcLen, :Average) + (; itp_raw, itp_dim) + elseif Itp == CubicHermiteSpline + itp_raw = CubicHermiteSpline(dA_vec, A_vec, x_vec) + itp_dim = CubicHermiteSpline(dA_vec, A_dimvec) + (; itp_raw, itp_dim) + elseif Itp == QuinticHermiteSpline + itp_raw = QuinticHermiteSpline(ddA_vec, dA_vec, A_vec, x_vec) + itp_dim = QuinticHermiteSpline(ddA_vec, dA_vec, A_dimvec) + (; itp_raw, itp_dim) + else + itp_raw = Itp(A_vec, x_vec) + itp_dim = Itp(A_dimvec) + (; itp_raw, itp_dim) + end + + @testset "$x" for x in x_vec_dense + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec), + itp_dim(x_vec) + ) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_wider) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_left) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_right) + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp_raw = QuadraticInterpolation(A_vec, x_vec, :Backward) + itp_dim = QuadraticInterpolation(A_dimvec, :Backward) + @testset "$x" for x in x_vec + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec), + itp_dim(x_vec) + ) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_wider) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_left) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_right) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp_raw = ConstantInterpolation(A_vec, x_vec; dir = :Right) + itp_dim = ConstantInterpolation(A_dimvec; dir = :Right) + @testset "$x" for x in x_vec + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec), + itp_dim(x_vec) + ) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_wider) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_left) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_right) + end + end +end + +@testset "Preserves Extrapolation" begin + @testset "$Xtp" for Xtp in extrapolators + @testset "$Itp" for Itp in interpolators + + (; itp_raw, itp_dim) = if Itp == BSplineInterpolation + itp_raw = BSplineInterpolation( + A_vec, x_vec, 3, :ArcLen, :Average; + extrapolation = Xtp + ) + itp_dim = BSplineInterpolation( + A_dimvec, 3, :ArcLen, :Average; + extrapolation = Xtp + ) + (; itp_raw, itp_dim) + + elseif Itp == CubicHermiteSpline + itp_raw = CubicHermiteSpline( + dA_vec, A_vec, x_vec; + extrapolation = Xtp + ) + itp_dim = CubicHermiteSpline( + dA_vec, A_dimvec; + extrapolation = Xtp + ) + (; itp_raw, itp_dim) + + elseif Itp == QuinticHermiteSpline + itp_raw = QuinticHermiteSpline( + ddA_vec, dA_vec, A_vec, x_vec; + extrapolation = Xtp + ) + itp_dim = QuinticHermiteSpline( + ddA_vec, dA_vec, A_dimvec; + extrapolation = Xtp + ) + (; itp_raw, itp_dim) + + else + itp_raw = Itp( + A_vec, x_vec; + extrapolation = Xtp + ) + itp_dim = Itp( + A_dimvec; + extrapolation = Xtp + ) + (; itp_raw, itp_dim) + end + + if (Itp, Xtp) == ( + SmoothedConstantInterpolation, + ExtrapolationType.Linear + ) + @testset "$x" for x in x_vec_wider + if x > maximum(x_vec) + @test_broken isequal( + itp_raw(x), + itp_dim(x) + ) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test_broken isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + @test isequal( + itp_raw(x_vec_left), + itp_dim(x_vec_left) + ) + end + end + + else + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp_raw = QuadraticInterpolation( + A_vec, x_vec, :Backward; + extrapolation = Xtp + ) + itp_dim = QuadraticInterpolation( + A_dimvec, :Backward; + extrapolation = Xtp + ) + + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Inputs" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp_raw = ConstantInterpolation( + A_vec, x_vec; dir = :Right, + extrapolation = Xtp + ) + itp_dim = ConstantInterpolation( + A_dimvec; dir = :Right, + extrapolation = Xtp + ) + + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Inputs" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + end +end + +@testset "Preserves Left Extrapolation" begin + @testset "$Ext" for Ext in extrapolators + @testset "$Itp" for Itp in interpolators + + (; itp_raw, itp_dim) = if Itp == BSplineInterpolation + itp_raw = BSplineInterpolation( + A_vec, x_vec, 3, :ArcLen, :Average; + extrapolation_left = Ext + ) + itp_dim = BSplineInterpolation( + A_dimvec, 3, :ArcLen, :Average; + extrapolation_left = Ext + ) + (; itp_raw, itp_dim) + + elseif Itp == CubicHermiteSpline + itp_raw = CubicHermiteSpline( + dA_vec, A_vec, x_vec; + extrapolation_left = Ext + ) + itp_dim = CubicHermiteSpline( + dA_vec, A_dimvec; + extrapolation_left = Ext + ) + (; itp_raw, itp_dim) + + elseif Itp == QuinticHermiteSpline + itp_raw = QuinticHermiteSpline( + ddA_vec, dA_vec, A_vec, x_vec; + extrapolation_left = Ext + ) + itp_dim = QuinticHermiteSpline( + ddA_vec, dA_vec, A_dimvec; + extrapolation_left = Ext + ) + (; itp_raw, itp_dim) + + else + itp_raw = Itp( + A_vec, x_vec; + extrapolation_left = Ext + ) + itp_dim = Itp( + A_dimvec; + extrapolation_left = Ext + ) + (; itp_raw, itp_dim) + end + + @testset "$x" for x in x_vec_wider + if x > maximum(x_vec) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_left), + itp_dim(x_vec_left) + ) + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp_raw = QuadraticInterpolation( + A_vec, x_vec, :Backward; + extrapolation_left = Ext + ) + itp_dim = QuadraticInterpolation( + A_dimvec, :Backward; + extrapolation_left = Ext + ) + + @testset "$x" for x in x_vec_wider + if x > maximum(x_vec) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_left), + itp_dim(x_vec_left) + ) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp_raw = ConstantInterpolation( + A_vec, x_vec; dir = :Right, + extrapolation_left = Ext + ) + itp_dim = ConstantInterpolation( + A_dimvec; dir = :Right, + extrapolation_left = Ext + ) + + @testset "$x" for x in x_vec_wider + if x > maximum(x_vec) + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.RightExtrapolationError itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_left), + itp_dim(x_vec_left) + ) + end + end + end +end + +@testset "Preserves Right Extrapolation" begin + @testset "$Xtp" for Xtp in extrapolators + @testset "$Itp" for Itp in interpolators + + (; itp_raw, itp_dim) = if Itp == BSplineInterpolation + itp_raw = BSplineInterpolation( + A_vec, x_vec, 3, :ArcLen, :Average; + extrapolation_right = Xtp + ) + itp_dim = BSplineInterpolation( + A_dimvec, 3, :ArcLen, :Average; + extrapolation_right = Xtp + ) + (; itp_raw, itp_dim) + + elseif Itp == CubicHermiteSpline + itp_raw = CubicHermiteSpline( + dA_vec, A_vec, x_vec; + extrapolation_right = Xtp + ) + itp_dim = CubicHermiteSpline( + dA_vec, A_dimvec; + extrapolation_right = Xtp + ) + (; itp_raw, itp_dim) + + elseif Itp == QuinticHermiteSpline + itp_raw = QuinticHermiteSpline( + ddA_vec, dA_vec, A_vec, x_vec; + extrapolation_right = Xtp + ) + itp_dim = QuinticHermiteSpline( + ddA_vec, dA_vec, A_dimvec; + extrapolation_right = Xtp + ) + (; itp_raw, itp_dim) + + else + itp_raw = Itp( + A_vec, x_vec; + extrapolation_right = Xtp + ) + itp_dim = Itp( + A_dimvec; + extrapolation_right = Xtp + ) + (; itp_raw, itp_dim) + end + + if (Itp, Xtp) == ( + SmoothedConstantInterpolation, + ExtrapolationType.Linear + ) + @testset "$x" for x in x_vec_wider + if x < minimum(x_vec) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x) + elseif x > maximum(x_vec) + @test_broken isequal( + itp_raw(x), + itp_dim(x) + ) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + end + + @testset "Vector Input" begin + @test_broken isequal( + itp_raw(x_vec_right), + itp_dim(x_vec_right) + ) + end + + else + @testset "$x" for x in x_vec_wider + if x < minimum(x_vec) + @test_throws DataInterpolations.LeftExtrapolationError itp_raw(x) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.LeftExtrapolationError itp_raw(x_vec_wider) + @test isequal( + itp_raw(x_vec_right), + itp_dim(x_vec_right) + ) + end + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp_raw = QuadraticInterpolation( + A_vec, x_vec, :Backward; + extrapolation_right = Xtp + ) + itp_dim = QuadraticInterpolation( + A_dimvec, :Backward; + extrapolation_right = Xtp + ) + + @testset "$x" for x in x_vec_wider + if x < minimum(x_vec) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_right), + itp_dim(x_vec_right) + ) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp_raw = ConstantInterpolation( + A_vec, x_vec; dir = :Right, + extrapolation_right = Xtp + ) + itp_dim = ConstantInterpolation( + A_dimvec; dir = :Right, + extrapolation_right = Xtp + ) + + @testset "$x" for x in x_vec_wider + if x < minimum(x_vec) + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_throws DataInterpolations.LeftExtrapolationError itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_right), + itp_dim(x_vec_right) + ) + end + end + end +end + +@testset "Preserves Mixed Extrapolation" begin + @testset "$XtpLeft" for XtpLeft in extrapolators + @testset "$XtpRight" for XtpRight in extrapolators + @testset "$Itp" for Itp in interpolators + + (; itp_raw, itp_dim) = if Itp == BSplineInterpolation + itp_raw = BSplineInterpolation( + A_vec, x_vec, 3, :ArcLen, :Average; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = BSplineInterpolation( + A_dimvec, 3, :ArcLen, :Average; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + (; itp_raw, itp_dim) + + elseif Itp == CubicHermiteSpline + itp_raw = CubicHermiteSpline( + dA_vec, A_vec, x_vec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = CubicHermiteSpline( + dA_vec, A_dimvec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + (; itp_raw, itp_dim) + + elseif Itp == QuinticHermiteSpline + itp_raw = QuinticHermiteSpline( + ddA_vec, dA_vec, A_vec, x_vec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = QuinticHermiteSpline( + ddA_vec, dA_vec, A_dimvec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + (; itp_raw, itp_dim) + + else + itp_raw = Itp( + A_vec, x_vec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = Itp( + A_dimvec; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + (; itp_raw, itp_dim) + end + + @testset "$x" for x in x_vec_wider + if (Itp, XtpRight) == ( + SmoothedConstantInterpolation, + ExtrapolationType.Linear + ) && x > maximum(x_vec) + @test_broken isequal( + itp_raw(x), + itp_dim(x) + ) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + if (Itp, XtpRight) == ( + SmoothedConstantInterpolation, + ExtrapolationType.Linear + ) + @testset "$x" for x in x_vec_wider + if x > maximum(x_vec) + @test_broken isequal( + itp_raw(x), + itp_dim(x) + ) + else + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + end + + @testset "Vector Input" begin + @test_broken itp_dim(x_vec_wider) + + @test isequal( + itp_raw(x_vec_left), + itp_dim(x_vec_left) + ) + end + + else + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + end + + @testset "Backward `QuadraticInterpolation`" begin + itp_raw = QuadraticInterpolation( + A_vec, x_vec, :Backward; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = QuadraticInterpolation( + A_dimvec, :Backward; + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + + @testset "Right `ConstantInterpolation`" begin + itp_raw = ConstantInterpolation( + A_vec, x_vec; dir = :Right, + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + itp_dim = ConstantInterpolation( + A_dimvec; dir = :Right, + extrapolation_left = XtpLeft, + extrapolation_right = XtpRight + ) + + @testset "$x" for x in x_vec_wider + @test isequal( + itp_raw(x), + itp_dim(x) + ) + end + + @testset "Vector Input" begin + @test isequal( + itp_raw(x_vec_wider), + itp_dim(x_vec_wider) + ) + end + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index bcada46ad..ddf8907f7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -43,7 +43,7 @@ end @time @safetestset "categorical" begin include("categorical.jl") end @time @safetestset "xarray" begin include("xarray.jl") end @time @safetestset "chainrules" begin include("chainrules.jl") end - +@time @safetestset "interpolation" begin include("interpolation.jl") end if Sys.islinux() # Unfortunately this can hang on other platforms.