Use FastInterpolations.jl as the pdf interpolation backend

Project.toml

@@ -7,13 +7,15 @@ 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"
 
 [compat]
 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"

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.

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

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

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

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)

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