diff --git a/Project.toml b/Project.toml index fa75d048..eb908588 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ authors = ["Simon Byrne and various contributors"] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FFTA = "b86e33f2-c0db-4aa1-a6e0-ab43e668529e" +FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -14,6 +15,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Distributions = "0.23, 0.24, 0.25" DocStringExtensions = "0.8, 0.9" FFTA = "0.3" +FastInterpolations = "0.4.14" Interpolations = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16" StatsBase = "0.33, 0.34" julia = "1.10" diff --git a/README.md b/README.md index c671a9c0..a0b71177 100644 --- a/README.md +++ b/README.md @@ -78,26 +78,39 @@ The `BivariateKDE` object `B` contains gridded coordinates (`B.x` and `B.y`) and estimate (`B.density`). ### Interpolation - The KDE objects are stored as gridded density values, with attached coordinates. These are typically sufficient for plotting (see above), but -intermediate values can be interpolated using the -[Interpolations.jl](https://github.com/tlycken/Interpolations.jl) package via the `pdf` method -(extended from Distributions.jl). +intermediate values can be evaluated with the `pdf` method (extended from +Distributions.jl), which by default interpolates through the +[FastInterpolations.jl](https://github.com/ProjectTorreyPines/FastInterpolations.jl) backend. -``` +```julia pdf(k::UnivariateKDE, x) pdf(k::BivariateKDE, x, y) ``` where `x` and `y` are real numbers or arrays. -If you are making multiple calls to `pdf`, it will be more efficient to -construct an intermediate `InterpKDE` to store the interpolation structure: +The default is a quadratic, C¹-continuous interpolation; you can select a different +scheme — and its boundary condition — with the `method` keyword: +```julia +import FastInterpolations as FI +pdf(k, x; method = FI.LinearInterp()) +pdf(k, x; method = FI.CubicInterp()) +pdf(k, x; method = FI.CubicInterp(bc = FI.ZeroCurvBC())) # with a boundary condition of your choice ``` -ik = InterpKDE(k) + +See the [FastInterpolations.jl docs](https://projecttorreypines.github.io/FastInterpolations.jl/stable/boundary-conditions/overview/) for the available boundary conditions. + +For repeated calls — or to use the Interpolations.jl backend instead — it is +more efficient to construct an interpolation object once and reuse it: + +```julia +ik = InterpKDE(k) # Interpolations.jl backend +ik = FastInterpKDE(k) # FastInterpolations.jl backend pdf(ik, x) ``` -`InterpKDE` will pass any extra arguments to `interpolate`. +- `InterpKDE` ([Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl)) passes any extra arguments to `interpolate`. +- `FastInterpKDE` ([FastInterpolations.jl](https://github.com/ProjectTorreyPines/FastInterpolations.jl)) takes a `method` keyword to select the interpolation scheme. \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 71b83a3e..36ce5ff0 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -81,9 +81,9 @@ estimate (`B.density`). The KDE objects are stored as gridded density values, with attached coordinates. These are typically sufficient for plotting (see above), but -intermediate values can be interpolated using the -[Interpolations.jl](https://github.com/tlycken/Interpolations.jl) package via the `pdf` method -(extended from Distributions.jl). +intermediate values can be evaluated with the `pdf` method (extended from +Distributions.jl), which by default interpolates through the +[FastInterpolations.jl](https://github.com/ProjectTorreyPines/FastInterpolations.jl) backend. ```julia pdf(k::UnivariateKDE, x) @@ -92,16 +92,29 @@ pdf(k::BivariateKDE, x, y) where `x` and `y` are real numbers or arrays. -If you are making multiple calls to `pdf`, it will be more efficient to -construct an intermediate `InterpKDE` to store the interpolation structure: +The default is a quadratic, C¹-continuous interpolation; you can select a different +scheme — and its boundary condition — with the `method` keyword: ```julia -ik = InterpKDE(k) -pdf(ik, x) +import FastInterpolations as FI +pdf(k, x; method = FI.LinearInterp()) +pdf(k, x; method = FI.CubicInterp()) +pdf(k, x; method = FI.CubicInterp(bc = FI.ZeroCurvBC())) # with a boundary condition of your choice ``` -`InterpKDE` will pass any extra arguments to `interpolate`. +See the [FastInterpolations.jl docs](https://projecttorreypines.github.io/FastInterpolations.jl/stable/boundary-conditions/overview/) for the available boundary conditions. + +For repeated calls — or to use the Interpolations.jl backend instead — it is +more efficient to construct an interpolation object once and reuse it: + +```julia +ik = InterpKDE(k) # Interpolations.jl backend +ik = FastInterpKDE(k) # FastInterpolations.jl backend (faster) +pdf(ik, x) +``` +- `InterpKDE` ([Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl)) passes any extra arguments to `interpolate`. +- `FastInterpKDE` ([FastInterpolations.jl](https://github.com/ProjectTorreyPines/FastInterpolations.jl)) takes a `method` keyword to select the interpolation scheme. ## API Reference ```@autodocs diff --git a/src/KernelDensity.jl b/src/KernelDensity.jl index 42d2f0b8..a1804faa 100644 --- a/src/KernelDensity.jl +++ b/src/KernelDensity.jl @@ -4,11 +4,12 @@ using DocStringExtensions: TYPEDEF, FIELDS using StatsBase using Distributions using Interpolations +using FastInterpolations import Distributions: twoπ, pdf import FFTA: rfft, irfft -export kde, kde_lscv, UnivariateKDE, BivariateKDE, InterpKDE, pdf +export kde, kde_lscv, UnivariateKDE, BivariateKDE, InterpKDE, FastInterpKDE, pdf abstract type AbstractKDE end @@ -18,5 +19,6 @@ Base.Broadcast.broadcastable(x::AbstractKDE) = Ref(x) include("univariate.jl") include("bivariate.jl") include("interp.jl") +include("fast_interp.jl") end # module diff --git a/src/fast_interp.jl b/src/fast_interp.jl new file mode 100644 index 00000000..2613a743 --- /dev/null +++ b/src/fast_interp.jl @@ -0,0 +1,39 @@ +const FI_DefaultInterpMethod = QuadraticInterp(bc = MinCurvFit()) + +mutable struct FastInterpKDE{K,I} <: AbstractKDE + kde::K + itp::I + FastInterpKDE{K,I}(kde::K, itp::I) where {K,I} = new{K,I}(kde, itp) +end + +function FastInterpKDE(kde::UnivariateKDE; method::AbstractInterpMethod = FI_DefaultInterpMethod) + itp = interp(kde.x, kde.density; method, extrap = FillExtrap(0.0)) + FastInterpKDE{typeof(kde), typeof(itp)}(kde, itp) +end + +function FastInterpKDE(kde::BivariateKDE; method::AbstractInterpMethod = FI_DefaultInterpMethod) + itp = interp((kde.x, kde.y), kde.density; method, extrap = FillExtrap(0.0)) + FastInterpKDE{typeof(kde), typeof(itp)}(kde, itp) +end + +pdf(fik::FastInterpKDE, x::Real...) = fik.itp(x...) +pdf(fik::FastInterpKDE, xs::AbstractVector) = fik.itp(xs) +pdf(fik::FastInterpKDE, xs::AbstractVector, ys::AbstractVector) = [fik.itp(x, y) for x in xs, y in ys] + + +# One-shot API: a direct query builds no persistent interpolant. +# 1D handles a scalar or a whole vector in a single one-shot call. +function pdf(k::UnivariateKDE, x; method::AbstractInterpMethod = FI_DefaultInterpMethod) + return interp(k.x, k.density, x; method = method, extrap = FillExtrap(0.0)) +end + +# 2D single point: one-shot (a local build, no full-grid interpolant) — the fast path. +function pdf(k::BivariateKDE, x::Real, y::Real; method::AbstractInterpMethod = FI_DefaultInterpMethod) + return interp((k.x, k.y), k.density, (x, y); method = method, extrap = FillExtrap(0.0)) +end + +# 2D grid: build the interpolant once, then evaluate the xs × ys grid (the ND one-shot +# has no efficient outer-product grid form). +function pdf(k::BivariateKDE, xs::AbstractVector, ys::AbstractVector; method::AbstractInterpMethod = FI_DefaultInterpMethod) + return pdf(FastInterpKDE(k; method = method), xs, ys) +end \ No newline at end of file diff --git a/src/interp.jl b/src/interp.jl index 05e648e6..8cfa6419 100644 --- a/src/interp.jl +++ b/src/interp.jl @@ -28,5 +28,3 @@ pdf(ik::InterpKDE,x::Real...) = ik.itp(x...) pdf(ik::InterpKDE,xs::AbstractVector) = [ik.itp(x) for x in xs] pdf(ik::InterpKDE,xs::AbstractVector,ys::AbstractVector) = [ik.itp(x,y) for x in xs, y in ys] -pdf(k::UnivariateKDE,x) = pdf(InterpKDE(k),x) -pdf(k::BivariateKDE,x,y) = pdf(InterpKDE(k),x,y) diff --git a/test/interp.jl b/test/interp.jl index ed93d914..eadc8219 100644 --- a/test/interp.jl +++ b/test/interp.jl @@ -1,5 +1,8 @@ using Test using KernelDensity +using Distributions +import FastInterpolations: QuadraticInterp, CubicInterp +import Interpolations X = randn(100) Y = randn(100) @@ -23,3 +26,116 @@ k = kde(([0.0],[0.0]), bandwidth=(1.0, 1.0)) @test pdf(k, +10.0, 0.0) ≤ pdf(MvNormal(2, 1.0), [+10.0, 0.0]) @test pdf(k, 0.0, -10.0) ≤ pdf(MvNormal(2, 1.0), [0.0, -10.0]) @test pdf(k, 0.0, +10.0) ≤ pdf(MvNormal(2, 1.0), [0.0, +10.0]) + +@testset "FastInterpKDE" begin + X = randn(100); Y = randn(100) + + # nodal exactness (1D) + k = kde(X) + fik = FastInterpKDE(k) + @test pdf(fik, k.x) ≈ k.density + + # nodal exactness (2D) + k2 = kde((X, Y)) + fik2 = FastInterpKDE(k2) + @test pdf(fik2, k2.x, k2.y) ≈ k2.density + + # tail safety + extrapolation → exactly 0 (1D) + kc = kde([0.0]; bandwidth = 1.0) + fkc = FastInterpKDE(kc) + @test pdf(fkc, kc.x) ≈ kc.density + @test pdf(fkc, -10.0) == 0.0 + @test pdf(fkc, +10.0) == 0.0 + @test pdf(fkc, -10.0) ≤ pdf(Normal(), -10.0) + @test pdf(fkc, +10.0) ≤ pdf(Normal(), +10.0) + + # tail safety + extrapolation → exactly 0 (2D) + kd = kde(([0.0], [0.0]); bandwidth = (1.0, 1.0)) + fkd = FastInterpKDE(kd) + @test pdf(fkd, kd.x, kd.y) ≈ kd.density + @test pdf(fkd, -10.0, 0.0) == 0.0 + @test pdf(fkd, 0.0, +10.0) == 0.0 + + # method swap stays nodal-exact + @test pdf(FastInterpKDE(k; method = CubicInterp()), k.x) ≈ k.density +end + +@testset "one-shot pdf" begin + X = randn(100); Y = randn(100) + k = kde(X); k2 = kde((X, Y)) + + # direct queries reproduce the grid nodes + @test pdf(k, k.x) ≈ k.density # 1D vector (one-shot) + @test pdf(k, 0.37) isa Real # 1D scalar (one-shot) + @test pdf(k2, k2.x, k2.y) ≈ k2.density # 2D grid (build-once) + @test pdf(k2, 0.1, 0.2) isa Real # 2D scalar (one-shot fast path) + + # extrapolation outside the grid → exactly 0 + @test pdf(k, first(k.x) - 1.0) == 0.0 + @test pdf(k2, first(k2.x) - 1.0, 0.0) == 0.0 + + # method keyword threads through + @test pdf(k, k.x; method = QuadraticInterp()) ≈ k.density +end + +@testset "type stability and parity with master (InterpKDE)" begin + k = kde(randn(80)); k2 = kde((randn(80), randn(80))) + + # 1D: scalar → Float64; vector/range (incl. length-1) → Vector{Float64} + @test (@inferred pdf(k, 0.5)) isa Float64 + @test (@inferred pdf(k, [0.5])) isa Vector{Float64} + @test (@inferred pdf(k, k.x)) isa Vector{Float64} + + # 2D: scalar → Float64; vectors/ranges → Matrix{Float64} + @test (@inferred pdf(k2, 0.1, 0.2)) isa Float64 + @test (@inferred pdf(k2, [0.1], [0.2])) isa Matrix{Float64} + @test (@inferred pdf(k2, k2.x, k2.y)) isa Matrix{Float64} + + # the same types come out of the InterpKDE path, so the migration is type-transparent + @test pdf(k, 0.5) isa Float64 + @test pdf(InterpKDE(k), 0.5) isa Float64 + @test pdf(k2, k2.x, k2.y) isa Matrix{Float64} + + # explicit FastInterpKDE wrapper is type-stable too + @test (@inferred pdf(FastInterpKDE(k), 0.5)) isa Float64 + @test (@inferred pdf(FastInterpKDE(k2), 0.1, 0.2)) isa Float64 +end + +# The bare pdf(k, …) convenience now routes to the FastInterpolations path, so these +# tests pin the original Interpolations.jl backend explicitly via InterpKDE to guard +# backward compatibility (it is unchanged in src/interp.jl). +@testset "InterpKDE (Interpolations.jl) backward compatibility" begin + X = randn(100); Y = randn(100) + + # 1D: nodal exactness + tail safety (the original suite's assertions) + k = kde(X) + ik = InterpKDE(k) + @test pdf(ik, k.x) ≈ k.density + + kc = kde([0.0]; bandwidth = 1.0) + ikc = InterpKDE(kc) + @test pdf(ikc, kc.x) ≈ kc.density + @test pdf(ikc, -10.0) ≤ pdf(Normal(), -10.0) + @test pdf(ikc, +10.0) ≤ pdf(Normal(), +10.0) + + # 2D: nodal exactness + tail safety + k2 = kde((X, Y)) + ik2 = InterpKDE(k2) + @test pdf(ik2, k2.x, k2.y) ≈ k2.density + + kd = kde(([0.0], [0.0]); bandwidth = (1.0, 1.0)) + ikd = InterpKDE(kd) + @test pdf(ikd, kd.x, kd.y) ≈ kd.density + @test pdf(ikd, -10.0, 0.0) ≤ pdf(MvNormal(2, 1.0), [-10.0, 0.0]) + @test pdf(ikd, 0.0, +10.0) ≤ pdf(MvNormal(2, 1.0), [0.0, +10.0]) + + # documented API: extra args are forwarded to Interpolations.interpolate + bspline = Interpolations.BSpline(Interpolations.Quadratic(Interpolations.Line(Interpolations.OnGrid()))) + @test pdf(InterpKDE(k, bspline), k.x) ≈ k.density + + # return types unchanged: scalar / vector / matrix + @test pdf(ik, 0.5) isa Float64 + @test pdf(ik, k.x) isa Vector{Float64} + @test pdf(ik2, 0.1, 0.2) isa Float64 + @test pdf(ik2, k2.x, k2.y) isa Matrix{Float64} +end