diff --git a/Project.toml b/Project.toml index 696a63d7..fa8bbe62 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RadialBasisFunctions" uuid = "79ee0514-adf7-4479-8807-6f72ea8967e8" authors = ["Kyle Beggs"] -version = "0.1.5" +version = "0.2.0" [deps] ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" diff --git a/README.md b/README.md index 42f0e9a6..0903a68e 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,7 @@ This package intends to provide tools for all things regarding Radial Basis Func | Feature | Status | | ------- | ------ | | Interpolation | ✅ | +| Regridding | ✅ | | Partial derivative ($\partial f$) | ✅ | | Laplacian ($\nabla^2 f$, $\Delta f$) | ✅ | | Gradient ($\nabla f$) | ✅ | @@ -18,7 +19,7 @@ This package intends to provide tools for all things regarding Radial Basis Func | Custom / user supplied ($\mathcal{L} f$) | ✅ | | divergence ($\textrm{div} \mathbf{F}$ or $\nabla \cdot \mathbf{F}$) | ❌ | | curl ($\nabla \times \mathbf{F}$) | ❌ | -| Reduced Order Models | ❌ | +| Reduced Order Models (i.e. POD) | ❌ | Currently, we support the following types of RBFs (all have polynomial augmentation by default, but is optional) @@ -47,4 +48,4 @@ Simply install the latest stable release using Julia's package manager: That said, we currently only support the second option here (`Vector{AbstractVector}`), but plan to support matrix inputs in the future. -2. `RadialBasisInterp` uses all points, but there are plans to support local collocation / subdomains like the operators use. +2. `Interpolator` uses all points, but there are plans to support local collocation / subdomains like the operators use. diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index 7a2b0145..31d8f157 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -20,7 +20,7 @@ y = f.(x) and now we can build the interpolator ```@example overview -interp = RadialBasisInterp(x, y) +interp = Interpolator(x, y) ``` and evaluate it at a new point @@ -37,10 +37,10 @@ and compare the error abs.(y_true .- y_new) ``` -Wow! The error is numerically zero! Well... we set ourselves up for success here. `RadialBasisInterp` (along with `RadialBasisOperator`) has an optional argument to provide the type of radial basis including the degree of polynomial augmentation. The default basis is a cubic polyharmonic spline with 2nd degree polynomial augmentation (which the constructor is `PHS(3, poly_deg=2)`) and given the underlying function we are interpolating is a 2nd order polynomial itself, we are able to represent it exactly (up to machine precision). Let's see what happens when we only use 1st order polynomial augmentation +Wow! The error is numerically zero! Well... we set ourselves up for success here. `Interpolator` (along with `RadialBasisOperator`) has an optional argument to provide the type of radial basis including the degree of polynomial augmentation. The default basis is a cubic polyharmonic spline with 2nd degree polynomial augmentation (which the constructor is `PHS(3, poly_deg=2)`) and given the underlying function we are interpolating is a 2nd order polynomial itself, we are able to represent it exactly (up to machine precision). Let's see what happens when we only use 1st order polynomial augmentation ```@example overview -interp = RadialBasisInterp(x, y, PHS(3, poly_deg=1)) +interp = Interpolator(x, y, PHS(3, poly_deg=1)) y_new = interp(x_new) abs.(y_true .- y_new) ``` diff --git a/docs/src/index.md b/docs/src/index.md index da116098..683c8cd6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,6 +11,7 @@ This package intends to provide tools for all things regarding Radial Basis Func | Feature | Status | | ------- | ------ | | Interpolation | ✅ | +| Regridding | ✅ | | Partial derivative ($\partial f$) | ✅ | | Laplacian ($\nabla^2 f$, $\Delta f$) | ✅ | | Gradient ($\nabla f$) | ✅ | @@ -18,7 +19,7 @@ This package intends to provide tools for all things regarding Radial Basis Func | Custom / user supplied ($\mathcal{L} f$) | ✅ | | divergence ($\textrm{div} \mathbf{F}$ or $\nabla \cdot \mathbf{F}$) | ❌ | | curl ($\nabla \times \mathbf{F}$) | ❌ | -| Reduced Order Models | ❌ | +| Reduced Order Models (i.e. POD) | ❌ | Currently, we support the following types of RBFs (all have polynomial augmentation by default, but is optional) @@ -39,13 +40,6 @@ Simply install the latest stable release using Julia's package manager: ] add RadialBasisFunctions ``` -## Planned Features - -* Adaptive operators and interpolation. Adding / removing / modifying points and automatically updating the weights without a complete recalculation. -* Add more built-in operator combinations that will allow you to lazily construct operators such as - * divergence ($\textrm{div} \mathbf{F}$ or $\nabla \cdot \mathbf{F}$) - * curl ($\nabla \times \mathbf{F}$) - ## Current Limitations 1. A critical dependency of this package is [NearestNeighbors.jl](https://github.com/KristofferC/NearestNeighbors.jl) which requires that the dimension of each data point is inferrable. To quote from NearestNeighbors.jl: @@ -55,4 +49,4 @@ Simply install the latest stable release using Julia's package manager: That said, we currently only support the second option here (`Vector{AbstractVector}`), but plan to support matrix inputs in the future. -2. `RadialBasisInterp` uses all points, but there are plans to support local collocation / subdomains like the operators use. +2. `Interpolator` uses all points, but there are plans to support local collocation / subdomains like the operators use. diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl index 41f59b54..946e597a 100644 --- a/src/RadialBasisFunctions.jl +++ b/src/RadialBasisFunctions.jl @@ -42,7 +42,10 @@ include("operators/monomial.jl") include("operators/operator_combinations.jl") include("interpolation.jl") -export RadialBasisInterp +export Interpolator + +include("operators/regridding.jl") +export Regrid, regrid const Δ = ∇² # some people like this notation for the Laplacian const DIV0_OFFSET = 1e-8 @@ -77,7 +80,7 @@ using PrecompileTools ∇²z = ∇²(z) # interpolation - interp = RadialBasisInterp(x, z, b) + interp = Interpolator(x, z, b) zz = interp([SVector{2}(rand(2)), SVector{2}(rand(2))]) end end diff --git a/src/interpolation.jl b/src/interpolation.jl index c535d6a8..c1c64703 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -1,9 +1,9 @@ """ - struct RadialBasisInterp + struct Interpolator Construct a radial basis interpolation. """ -struct RadialBasisInterp{X,Y,R,M,RB,MB} +struct Interpolator{X,Y,R,M,RB,MB} x::X y::Y rbf_weights::R @@ -13,25 +13,25 @@ struct RadialBasisInterp{X,Y,R,M,RB,MB} end """ - function RadialBasisInterp(x, y, basis::B=PHS()) + function Interpolator(x, y, basis::B=PHS()) Construct a radial basis interpolator. """ -function RadialBasisInterp(x, y, basis::B=PHS()) where {B<:AbstractRadialBasis} +function Interpolator(x, y, basis::B=PHS()) where {B<:AbstractRadialBasis} dim = length(first(x)) k = length(x) # number of data in influence/support domain npoly = binomial(dim + basis.poly_deg, basis.poly_deg) n = k + npoly mon = MonomialBasis(dim, basis.poly_deg) - A = Symmetric(zeros(eltype(first(x)), n, n)) + data_type = promote_type(eltype(first(x)), eltype(y)) + A = Symmetric(zeros(data_type, n, n)) _build_collocation_matrix!(A, x, basis, mon, k) - b = zeros(eltype(first(x)), n) - b[1:k] .= y + b = data_type[i < k ? y[i] : 0 for i in 1:n] w = A \ b - return RadialBasisInterp(x, y, w[1:k], w[(k + 1):end], basis, mon) + return Interpolator(x, y, w[1:k], w[(k + 1):end], basis, mon) end -function (rbfi::RadialBasisInterp)(x::T) where {T} +function (rbfi::Interpolator)(x::T) where {T} rbf = zero(eltype(T)) for i in eachindex(rbfi.rbf_weights) rbf += rbfi.rbf_weights[i] * rbfi.rbf_basis(x, rbfi.x[i]) @@ -47,11 +47,11 @@ function (rbfi::RadialBasisInterp)(x::T) where {T} return rbf + poly end -(rbfi::RadialBasisInterp)(x::Vector{<:AbstractVector}) = [rbfi(val) for val in x] +(rbfi::Interpolator)(x::Vector{<:AbstractVector}) = [rbfi(val) for val in x] # pretty printing -function Base.show(io::IO, op::RadialBasisInterp) - println(io, "RadialBasisInterp") +function Base.show(io::IO, op::Interpolator) + println(io, "Interpolator") println(io, "├─Input type: ", typeof(first(op.x))) println(io, "├─Output type: ", typeof(first(op.y))) println(io, "├─Number of points: ", length(op.x)) diff --git a/src/operators/operators.jl b/src/operators/operators.jl index 2eb7dd0a..e8b522bb 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -125,7 +125,6 @@ function LinearAlgebra.:⋅( return sum(op(x)) end -# TODO function LinearAlgebra.mul!( y::AbstractVector{<:Real}, op::RadialBasisOperator{<:VectorValuedOperator}, diff --git a/src/operators/regridding.jl b/src/operators/regridding.jl new file mode 100644 index 00000000..32582b35 --- /dev/null +++ b/src/operators/regridding.jl @@ -0,0 +1,27 @@ +""" + Regrid + +Builds an operator for interpolating from one set of points to another. +""" +struct Regrid + ℒ + Regrid() = new(identity) +end + +# convienience constructors +""" + function regrid(data, eval_points, order, dim, basis; k=autoselect_k(data, basis)) + +Builds a `RadialBasisOperator` where the operator is an regrid from one set of points to another, `data` -> `eval_points`. +""" +function regrid( + data::AbstractVector{D}, + eval_points::AbstractVector{D}, + basis::B=PHS(3; poly_deg=2); + k::T=autoselect_k(data, basis), +) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis} + return RadialBasisOperator(Regrid(), data, eval_points, basis; k=k) +end + +# pretty printing +print_op(op::Regrid) = "regrid" diff --git a/test/operators/interpolation.jl b/test/operators/interpolation.jl index db576f78..cbe42679 100644 --- a/test/operators/interpolation.jl +++ b/test/operators/interpolation.jl @@ -19,8 +19,8 @@ N = 1000 x = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), 1:N) y = franke.(x) -interp = RadialBasisInterp(x, y, PHS(3; poly_deg=2)) -@test interp isa RadialBasisInterp +interp = Interpolator(x, y, PHS(3; poly_deg=2)) +@test interp isa Interpolator xnew = SVector(0.5, 0.5) @test abs(interp(xnew) - franke(xnew)) < 1e-5 diff --git a/test/operators/regrid.jl b/test/operators/regrid.jl new file mode 100644 index 00000000..ceb5982c --- /dev/null +++ b/test/operators/regrid.jl @@ -0,0 +1,16 @@ +using RadialBasisFunctions +using StaticArrays +using Statistics +using Random + +rsme(test, correct) = sqrt(sum((test - correct) .^ 2) / sum(correct .^ 2)) +mean_percent_error(test, correct) = mean(abs.((test .- correct) ./ correct)) * 100 + +f(x) = 1 + sin(4 * x[1]) + cos(3 * x[1]) + sin(2 * x[2]) +N = 10_000 +x = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), 1:N) +y = f.(x) + +x2 = map(x -> SVector{2}(rand(2)), 1:100) +r = regrid(x, x2, PHS(3; poly_deg=2)) +@test mean_percent_error(r(y), f.(x2)) < 0.1 diff --git a/test/runtests.jl b/test/runtests.jl index 61ce900a..7b391d2c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,6 +36,10 @@ end include("operators/interpolation.jl") end +@safetestset "Regridding" begin + include("operators/regrid.jl") +end + @safetestset "Stencil" begin include("linalg/stencil.jl") end