Add Regrid

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"

README.md

@@ -11,14 +11,15 @@ 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$) | ✅ |
 | Directional Derivative ($\nabla f \cdot v$) | ✅ |
 | 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.

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)
 ```

docs/src/index.md

@@ -11,14 +11,15 @@ 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$) | ✅ |
 | Directional Derivative ($\nabla f \cdot v$) | ✅ |
 | 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.

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

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))

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},

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"

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

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

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