Merge pull request #27 from kylebeggs/feature/monomials

Slight refactor of monomial basis

.github/workflows/CI.yml

@@ -18,7 +18,6 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
           - '1.7'
           - '1.8'
           - '1.9'
@@ -28,7 +27,6 @@ jobs:
           - ubuntu-latest
         arch:
           - x64
-          - x86
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -1,7 +1,7 @@
 name = "RadialBasisFunctions"
 uuid = "79ee0514-adf7-4479-8807-6f72ea8967e8"
 authors = ["Kyle Beggs"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
@@ -22,7 +22,7 @@ NearestNeighbors = "0.4"
 PrecompileTools = "1"
 StaticArraysCore = "1.4"
 SymRCM = "0.2"
-julia = "1.6"
+julia = "1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/RadialBasisFunctions.jl

@@ -20,7 +20,7 @@ export MonomialBasis
 include("utils.jl")
 export find_neighbors, reorder_points!
 
-include("linalg/stencil.jl")
+include("solve.jl")
 
 include("operators/operators.jl")
 export RadialBasisOperator, update_weights!
@@ -40,7 +40,7 @@ export Directional, directional
 include("operators/virtual.jl")
 export ∂virtual
 
-include("operators/monomial.jl")
+include("operators/monomial/monomial.jl")
 
 include("operators/operator_combinations.jl")
 
@@ -63,7 +63,7 @@ using PrecompileTools
         basis_funcs = [
             IMQ(1),
             Gaussian(1),
-            PHS(1; poly_deg=-1),
+            PHS(1; poly_deg=0),
             PHS(3; poly_deg=0),
             PHS(5; poly_deg=1),
             PHS(7; poly_deg=2),

src/basis/basis.jl

@@ -1,29 +1,39 @@
 """
-    abstract type AbstractRadialBasis end
+    abstract type AbstractBasis end
 """
-abstract type AbstractRadialBasis end
+abstract type AbstractBasis end
+
+"""
+    abstract type AbstractRadialBasis <: AbstractBasis end
+"""
+abstract type AbstractRadialBasis <: AbstractBasis end
 
 struct ℒRadialBasisFunction{F<:Function}
     f::F
 end
 (ℒrbf::ℒRadialBasisFunction)(x, xᵢ) = ℒrbf.f(x, xᵢ)
 
-struct ℒMonomial{F<:Function}
+struct ℒMonomialBasis{Dim,Deg,F<:Function}
     f::F
+    function ℒMonomialBasis(dim::T, deg::T, f) where {T<:Int}
+        if deg < 0
+            throw(ArgumentError("Monomial basis must have non-negative degree"))
+        end
+        return new{dim,deg,typeof(f)}(f)
+    end
 end
-(ℒmon::ℒMonomial)(m, x) = ℒmon.f(m, x)
+function (ℒmon::ℒMonomialBasis{Dim,Deg})(x) where {Dim,Deg}
+    b = ones(_get_underlying_type(x), binomial(Dim + Deg, Dim))
+    ℒmon(b, x)
+    return b
+end
+(m::ℒMonomialBasis)(b, x) = m.f(b, x)
 
 include("polyharmonic_spline.jl")
 include("inverse_multiquadric.jl")
 include("gaussian.jl")
 include("monomial.jl")
 
-function ∂(basis::B, order::T, dim::T) where {T<:Int,B<:AbstractRadialBasis}
-    return ∂(basis, Val(order), dim)
-end
-∂(basis::B, dim::T) where {T<:Int,B} = ∂(basis, Val(1), dim)
-∂²(basis::B, dim::T) where {T<:Int,B} = ∂(basis, Val(2), dim)
-
 # pretty printing
 unicode_order(::Val{1}) = ""
 unicode_order(::Val{2}) = "²"

src/basis/gaussian.jl

@@ -18,7 +18,7 @@ end
 
 (rbf::Gaussian)(x, xᵢ) = exp(-(rbf.ε * euclidean(x, xᵢ))^2)
 
-function ∂(rbf::Gaussian, ::Val{1}, dim::Int)
+function ∂(rbf::Gaussian, dim::Int)
     function ∂ℒ(x, xᵢ)
         return -2 * rbf.ε^2 * (x[dim] - xᵢ[dim]) * exp(-rbf.ε^2 * sqeuclidean(x, xᵢ))
     end
@@ -30,7 +30,7 @@ function ∇(rbf::Gaussian)
     end
 end
 
-function ∂(rbf::Gaussian, ::Val{2}, dim::Int)
+function ∂²(rbf::Gaussian, dim::Int)
     function ∂²ℒ(x, xᵢ)
         ε2 = rbf.ε^2
         return (4 * ε2^2 * (x[dim] - xᵢ[dim])^2 - 2 * ε2) * exp(-ε2 * sqeuclidean(x, xᵢ))

src/basis/inverse_multiquadric.jl

@@ -15,7 +15,7 @@ end
 
 (rbf::IMQ)(x, xᵢ) = 1 / sqrt((euclidean(x, xᵢ) * rbf.ε)^2 + 1)
 
-function ∂(rbf::IMQ, ::Val{1}, dim::Int=1)
+function ∂(rbf::IMQ, dim::Int=1)
     function ∂ℒ(x, xᵢ)
         ε2 = rbf.ε .^ 2
         return (xᵢ[dim] - x[dim]) .* (ε2 / sqrt((ε2 * sqeuclidean(x, xᵢ) + 1)^3))
@@ -29,7 +29,7 @@ function ∇(rbf::IMQ)
     end
 end
 
-function ∂(rbf::IMQ, ::Val{2}, dim::Int=1)
+function ∂²(rbf::IMQ, dim::Int=1)
     function ∂²ℒ(x, xᵢ)
         ε2 = rbf.ε .^ 2
         ε4 = ε2^2

src/basis/monomial.jl

@@ -1,39 +1,100 @@
 """
-    struct MonomialBasis{T<:Int,B<:Function}
+    struct MonomialBasis{Dim,Deg} <: AbstractBasis
 
-Multivariate Monomial basis.
-n ∈ N: length of array, i.e., x ∈ Rⁿ
-deg ∈ N: degree
+`Dim` dimensional monomial basis of order `Deg`.
 """
-struct MonomialBasis{T<:Int,B}
-    n::T
-    deg::T
-    basis!::B
-    function MonomialBasis(n::T, deg::T) where {T<:Int}
+struct MonomialBasis{Dim,Deg,F<:Function} <: AbstractBasis
+    f::F
+    function MonomialBasis(dim::T, deg::T) where {T<:Int}
         if deg < 0
-            basis! = nothing
-        else
-            basis! = build_monomial_basis(n, deg)
+            throw(ArgumentError("Monomial basis must have non-negative degree"))
         end
-        return new{T,typeof(basis!)}(n, deg, basis!)
+        f = _get_monomial_basis(Val(dim), Val(deg))
+        return new{dim,deg,typeof(f)}(f)
     end
 end
 
-function (m::MonomialBasis)(x::AbstractVector{T}) where {T}
-    b = ones(T, binomial(m.n + m.deg, m.n))
-    m.basis!(b, x)
+function (m::MonomialBasis{Dim,Deg})(x) where {Dim,Deg}
+    b = ones(_get_underlying_type(x), binomial(Dim + Deg, Dim))
+    m.f(b, x)
     return b
 end
-(m::MonomialBasis)(b, x) = m.basis!(b, x)
+(m::MonomialBasis)(b, x) = m.f(b, x)
 
-function build_monomial_basis(n::T, deg::T) where {T<:Int}
-    if deg == 0
-        basis! = (b, x) -> b .= one(eltype(x))
-        return basis!
+for Dim in (:1, :2, :3)
+    @eval begin
+        function _get_monomial_basis(::Val{$Dim}, ::Val{0})
+            return function basis!(b, x)
+                b[1] = one(_get_underlying_type(x))
+                return nothing
+            end
+        end
+    end
+end
+
+function _get_monomial_basis(::Val{1}, ::Val{N}) where {N}
+    return function basis!(b, x)
+        b[1] = one(_get_underlying_type(x))
+        if N > 0
+            for n in 1:N
+                b[n + 1] = only(x)^n
+            end
+        end
+        return nothing
+    end
+end
+
+function _get_monomial_basis(::Val{2}, ::Val{1})
+    return function basis!(b, x)
+        b[1] = one(eltype(x))
+        b[2] = x[1]
+        b[3] = x[2]
+        return nothing
+    end
+end
+
+function _get_monomial_basis(::Val{2}, ::Val{2})
+    return function basis!(b, x)
+        b[1] = one(eltype(x))
+        b[2] = x[1]
+        b[3] = x[2]
+        b[4] = x[1] * x[2]
+        b[5] = x[1] * x[1]
+        b[6] = x[2] * x[2]
+        return nothing
+    end
+end
+
+function _get_monomial_basis(::Val{3}, ::Val{1})
+    return function basis!(b, x)
+        b[1] = one(eltype(x))
+        b[2] = x[1]
+        b[3] = x[2]
+        b[4] = x[3]
+        return nothing
     end
-    e = multiexponents(n + 1, deg)
-    e = map(i -> getindex.(e, i), 1:n)
-    ids = map(j -> map(i -> findall(x -> x >= i, e[j]), 1:deg), 1:n)
+end
+
+function _get_monomial_basis(::Val{3}, ::Val{2})
+    return function basis!(b, x)
+        b[1] = one(eltype(x))
+        b[2] = x[1]
+        b[3] = x[2]
+        b[4] = x[3]
+        b[5] = x[1] * x[2]
+        b[6] = x[1] * x[3]
+        b[7] = x[2] * x[3]
+        b[8] = x[1] * x[1]
+        b[9] = x[2] * x[2]
+        b[10] = x[3] * x[3]
+        return nothing
+    end
+end
+
+function _get_monomial_basis(::Val{Dim}, ::Val{Deg}) where {Dim,Deg}
+    e = multiexponents(Dim + 1, Deg)
+    e = map(i -> getindex.(e, i), 1:Dim)
+    ids = map(j -> map(i -> findall(x -> x >= i, e[j]), 1:Deg), 1:Dim)
     return build_monomial_basis(ids)
 end
 
@@ -49,7 +110,6 @@ function build_monomial_basis(ids::Vector{Vector{Vector{T}}}) where {T<:Int}
     return basis!
 end
 
-# pretty printing
-function Base.show(io::IO, pb::MonomialBasis)
-    return print(io, "MonomialBasis, deg=$(pb.deg) in $(pb.n) dimensions")
+function Base.show(io::IO, ::MonomialBasis{Dim,Deg}) where {Dim,Deg}
+    return print(io, "MonomialBasis of degree $(Deg) in $(Dim) dimensions")
 end

src/basis/polyharmonic_spline.jl

@@ -38,15 +38,15 @@ struct PHS1{T<:Int} <: AbstractPHS
 end
 
 (phs::PHS1)(x, xᵢ) = euclidean(x, xᵢ)
-function ∂(::PHS1, ::Val{1}, dim::Int)
+function ∂(::PHS1, dim::Int)
     ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + DIV0_OFFSET)
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS1)
     ∇ℒ(x, xᵢ) = (x .- xᵢ) / euclidean(x, xᵢ)
     return ℒRadialBasisFunction(∇ℒ)
 end
-function ∂(::PHS1, ::Val{2}, dim::Int)
+function ∂²(::PHS1, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) /
                (euclidean(x, xᵢ)^3 + DIV0_OFFSET)
@@ -76,15 +76,15 @@ struct PHS3{T<:Int} <: AbstractPHS
 end
 
 (phs::PHS3)(x, xᵢ) = euclidean(x, xᵢ)^3
-function ∂(::PHS3, ::Val{1}, dim::Int)
+function ∂(::PHS3, dim::Int)
     ∂ℒ(x, xᵢ) = 3 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS3)
     ∇ℒ(x, xᵢ) = 3 * (x .- xᵢ) * euclidean(x, xᵢ)
     return ℒRadialBasisFunction(∇ℒ)
 end
-function ∂(::PHS3, ::Val{2}, dim::Int)
+function ∂²(::PHS3, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) /
                (euclidean(x, xᵢ) + DIV0_OFFSET)
@@ -114,15 +114,15 @@ struct PHS5{T<:Int} <: AbstractPHS
 end
 
 (phs::PHS5)(x, xᵢ) = euclidean(x, xᵢ)^5
-function ∂(::PHS5, ::Val{1}, dim::Int)
+function ∂(::PHS5, dim::Int)
     ∂ℒ(x, xᵢ) = 5 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)^3
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS5)
     ∇ℒ(x, xᵢ) = 5 * (x .- xᵢ) * euclidean(x, xᵢ)^3
     return ℒRadialBasisFunction(∇ℒ)
 end
-function ∂(::PHS5, ::Val{2}, dim::Int)
+function ∂²(::PHS5, dim::Int)
     return function ∂²ℒ(x, xᵢ)
         return 5 * euclidean(x, xᵢ) * (3 * (x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ))
     end
@@ -149,15 +149,15 @@ struct PHS7{T<:Int} <: AbstractPHS
 end
 
 (phs::PHS7)(x, xᵢ) = euclidean(x, xᵢ)^7
-function ∂(::PHS7, ::Val{1}, dim::Int)
+function ∂(::PHS7, dim::Int)
     ∂ℒ(x, xᵢ) = 7 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)^5
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS7)
     ∇ℒ(x, xᵢ) = 7 * (x .- xᵢ) * euclidean(x, xᵢ)^5
     return ℒRadialBasisFunction(∇ℒ)
 end
-function ∂(::PHS7, ::Val{2}, dim::Int)
+function ∂²(::PHS7, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 7 * euclidean(x, xᵢ)^3 * (5 * (x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ))
     end

src/interpolation.jl

@@ -40,8 +40,8 @@ function (rbfi::Interpolator)(x::T) where {T}
     poly = zero(eltype(T))
     if !isempty(rbfi.monomial_weights)
         val_poly = rbfi.monomial_basis(x)
-        for i in eachindex(rbfi.monomial_weights)
-            poly += rbfi.monomial_weights[i] * val_poly[i]
+        for (i, val) in enumerate(val_poly)
+            poly += rbfi.monomial_weights[i] * val
         end
     end
     return rbf + poly
@@ -59,6 +59,7 @@ function Base.show(io::IO, op::Interpolator)
         io,
         "└─Basis: ",
         print_basis(op.rbf_basis),
-        " with degree $(op.monomial_basis.deg) Monomial",
+        " with degree $(_get_deg(op.monomial_basis)) Monomial",
     )
 end
+_get_deg(::MonomialBasis{Dim,Deg}) where {Dim,Deg} = Deg

src/operators/directional.jl

@@ -20,11 +20,7 @@ function directional(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(length(first(data))) do dim
-        return let dim = dim
-            x -> ∂(x, 1, dim)
-        end
-    end
+    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Directional(f, v)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
@@ -42,11 +38,7 @@ function directional(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(length(first(data))) do dim
-        return let dim = dim
-            x -> ∂(x, 1, dim)
-        end
-    end
+    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Directional(f, v)
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end