add some tests for basis functions

src/basis/basis.jl

@@ -8,30 +8,6 @@ abstract type AbstractBasis end
 """
 abstract type AbstractRadialBasis <: AbstractBasis end
 
-struct ℒRadialBasisFunction{F<:Function}
-    f::F
-end
-(ℒrbf::ℒRadialBasisFunction)(x, xᵢ) = ℒrbf.f(x, xᵢ)
-
-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
-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)
-
-degree(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Deg
-dim(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Dim
-
 include("polyharmonic_spline.jl")
 include("inverse_multiquadric.jl")
 include("gaussian.jl")

src/basis/monomial.jl

@@ -21,6 +21,9 @@ function (m::MonomialBasis{Dim,Deg})(x) where {Dim,Deg}
 end
 (m::MonomialBasis)(b, x) = m.f(b, x)
 
+degree(::MonomialBasis{Dim,Deg}) where {Dim,Deg} = Deg
+dim(::MonomialBasis{Dim,Deg}) where {Dim,Deg} = Dim
+
 for Dim in (:1, :2, :3)
     @eval begin
         function _get_monomial_basis(::Val{$Dim}, ::Val{0})

src/basis/polyharmonic_spline.jl

@@ -123,13 +123,13 @@ function ∇(::PHS5)
     return ℒRadialBasisFunction(∇ℒ)
 end
 function ∂²(::PHS5, dim::Int)
-    return function ∂²ℒ(x, xᵢ)
+    function ∂²ℒ(x, xᵢ)
         return 5 * euclidean(x, xᵢ) * (3 * (x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ))
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS5)
-    return function ∇²ℒ(x, xᵢ)
+    function ∇²ℒ(x, xᵢ)
         return sum(5 * euclidean(x, xᵢ) * (3 * (x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)))
     end
     return ℒRadialBasisFunction(∇²ℒ)

src/operators/monomial/monomial.jl

@@ -1,3 +1,28 @@
+struct ℒRadialBasisFunction{F<:Function}
+    f::F
+end
+(ℒrbf::ℒRadialBasisFunction)(x, xᵢ) = ℒrbf.f(x, xᵢ)
+
+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
+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)
+
+degree(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Deg
+dim(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Dim
+
+# include partial definitions for monomials
 include("partial.jl")
 ∂(mb::MonomialBasis, differentiation_dim::Int) = ∂(mb, Val(differentiation_dim))
 ∂²(mb::MonomialBasis, differentiation_dim::Int) = ∂²(mb, Val(differentiation_dim))
@@ -77,3 +102,7 @@
 function monomial_recursive_list(m::MonomialBasis, me::Vector{<:Monomial})
     return [monomial_recursive_list(m, me[i]) for i in eachindex(me)]
 end
+
+function Base.show(io::IO, ::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg}
+    return print(io, "ℒMonomialBasis of degree $(Deg) in $(Dim) dimensions")
+end

src/operators/virtual.jl

@@ -26,12 +26,12 @@ function ∂virtual(
     return if backward
         li = regrid(data, eval_points .- Ref(dx), basis; k=k)
         update_weights!(li)
-        w = columnwise_div(self.weights .- li.weights, Δ)
+        w = (self.weights .- li.weights) ./ Δ
         x -> w * x
     else # forward difference
         ri = regrid(data, eval_points .+ Ref(dx), basis; k=k)
         update_weights!(ri)
-        w = columnwise_div((ri.weights .- self.weights), Δ)
+        w = (ri.weights .- self.weights) ./ Δ
         x -> w * x
     end
 end

src/utils.jl

@@ -10,24 +10,6 @@ function find_neighbors(data::AbstractVector, eval_points::AbstractVector, k::In
     return adjl
 end
 
-function get_num_params(f)
-    return length.(getfield.(getfield.(methods(f), :sig), :parameters)) .- 1
-end
-
-function check_num_params(ℒ, basis::B) where {B<:AbstractRadialBasis}
-    num_params_ϕ = get_num_params(ℒ(basis.ϕ))
-    !any(num_params_ϕ .== 2) && throw(ArgumentError("ϕ must have 2 input arguments."))
-    #num_params_poly = get_num_params(ℒ(basis.poly))
-    #!any(num_params_poly .== 1) && throw(ArgumentError("poly must have 1 input argument."))
-    return nothing
-end
-
-function check_if_deg_odd(deg::Int)
-    if isodd(deg)
-        @warn "Monomial degree is recommended to be even, unless it is -1 which indicates no Monomials. (Flyer, 2016 - https://doi.org/10.1016/j.jcp.2016.05.026)"
-    end
-end
-
 """
     autoselect_k(data::Vector, basis<:AbstractRadialBasis)
 
@@ -58,36 +40,5 @@ function check_poly_deg(poly_deg)
     return nothing
 end
 
-function scale_cloud(data)
-    furthest_point = maximum(p -> euclidean(first(data), p), data)
-    return data ./ furthest_point
-end
-
-_allocate_weights(m, n, k) = _allocate_weights(Float64, m, n, k)
-function _allocate_weights(T, m, n, k)
-    return spzeros(T, m, n)
-end
-
-function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)
-    I, J, V = findnz(A)
-    for idx in eachindex(V)
-        V[idx] /= B[I[idx]]
-    end
-    return sparse(I, J, V)
-end
-columnwise_div(A::SparseMatrixCSC, B::Number) = A ./ B
-
-function _find_smallest_dist(data, k)
-    tree = KDTree(data)
-    _, dists = knn(tree, data, k, true)
-    Δ = minimum(dists) do d
-        z = minimum(@view(d[2:end])) do i
-            abs(i - first(d))
-        end
-        return z
-    end
-    return Δ
-end
-
 _get_underlying_type(x::AbstractVector) = eltype(x)
 _get_underlying_type(x::Number) = typeof(x)

test/basis/basis.jl

@@ -0,0 +1,17 @@
+using RadialBasisFunctions
+import RadialBasisFunctions as RBF
+
+@testset "Pretty Printing" begin
+    superscripts = ("", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹")
+    @test all(RBF.unicode_order.(Val.(1:9)) .== superscripts)
+end
+
+@testset "Basis - General Utils" begin
+    m = MonomialBasis(2, 3)
+    @test dim(m) == 2
+    @test degree(m) == 3
+
+    ∂m = RBF.∂(m, 1)
+    @test dim(∂m) == 2
+    @test degree(∂m) == 3
+end

test/basis/monomial.jl

@@ -2,52 +2,57 @@ using RadialBasisFunctions
 import RadialBasisFunctions as RBF
 using StaticArrays
 
-@testset "dim=1, deg=0, vector input" begin
-    x = SVector(2.0)
-
-    m = MonomialBasis(1, 0)
-    @test m isa MonomialBasis
-    @test typeof(m) <: MonomialBasis{1,0}
-
-    # standard evaluation
-    @test all(isapprox.(m(x), [1]))
-
-    # derivatives
-    @test all(isapprox.(RBF.∂(m, 1)(x), [0]))
-    @test all(isapprox.(RBF.∂²(m, 1)(x), [0]))
-    @test all(isapprox.(RBF.∇²(m)(x), [0]))
+@testset "dim=1, deg=0" begin
+    inputs = (SVector(2.0), 2.0)
+    foreach(inputs) do x
+        m = MonomialBasis(1, 0)
+        @test m isa MonomialBasis
+        @test typeof(m) <: MonomialBasis{1,0}
+
+        # standard evaluation
+        @test all(isapprox.(m(x), [1]))
+
+        # derivatives
+        @test all(isapprox.(RBF.∂(m, 1)(x), [0]))
+        @test all(isapprox.(RBF.∂²(m, 1)(x), [0]))
+        @test all(isapprox.(RBF.∇²(m)(x), [0]))
+    end
 end
 
 @testset "dim=1, deg=1" begin
-    x = SVector(2.0)
+    inputs = (SVector(2.0), 2.0)
 
-    m = MonomialBasis(1, 1)
-    @test m isa MonomialBasis
-    @test typeof(m) <: MonomialBasis{1,1}
+    foreach(inputs) do x
+        m = MonomialBasis(1, 1)
+        @test m isa MonomialBasis
+        @test typeof(m) <: MonomialBasis{1,1}
 
-    # standard evaluation
-    @test all(isapprox.(m(x), [1, 2]))
+        # standard evaluation
+        @test all(isapprox.(m(x), [1, 2]))
 
-    # derivatives
-    @test all(isapprox.(RBF.∂(m, 1)(x), [0, 1]))
-    @test all(isapprox.(RBF.∂²(m, 1)(x), [0, 0]))
-    @test all(isapprox.(RBF.∇²(m)(x), [0, 0]))
+        # derivatives
+        @test all(isapprox.(RBF.∂(m, 1)(x), [0, 1]))
+        @test all(isapprox.(RBF.∂²(m, 1)(x), [0, 0]))
+        @test all(isapprox.(RBF.∇²(m)(x), [0, 0]))
+    end
 end
 
 @testset "dim=1, deg=2" begin
-    x = SVector(2.0)
+    inputs = (SVector(2.0), 2.0)
 
-    m = MonomialBasis(1, 2)
-    @test m isa MonomialBasis
-    @test typeof(m) <: MonomialBasis{1,2}
+    foreach(inputs) do x
+        m = MonomialBasis(1, 2)
+        @test m isa MonomialBasis
+        @test typeof(m) <: MonomialBasis{1,2}
 
-    # standard evaluation
-    @test all(isapprox.(m(x), [1, 2, 4]))
+        # standard evaluation
+        @test all(isapprox.(m(x), [1, 2, 4]))
 
-    # derivatives
-    @test all(isapprox.(RBF.∂(m, 1)(x), [0, 1, 4]))
-    @test all(isapprox.(RBF.∂²(m, 1)(x), [0, 0, 2]))
-    @test all(isapprox.(RBF.∇²(m)(x), [0, 0, 2]))
+        # derivatives
+        @test all(isapprox.(RBF.∂(m, 1)(x), [0, 1, 4]))
+        @test all(isapprox.(RBF.∂²(m, 1)(x), [0, 0, 2]))
+        @test all(isapprox.(RBF.∇²(m)(x), [0, 0, 2]))
+    end
 end
 
 @testset "dim=2, deg=0" begin

test/runtests.jl

@@ -1,5 +1,9 @@
 using SafeTestsets
 
+@safetestset "Basis - General Utils" begin
+    include("basis/basis.jl")
+end
+
 @safetestset "Polyharmonic Splines" begin
     include("basis/polyharmonic_spline.jl")
 end