diff --git a/src/basis/basis.jl b/src/basis/basis.jl index 47858d67..a0894003 100644 --- a/src/basis/basis.jl +++ b/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") diff --git a/src/basis/monomial.jl b/src/basis/monomial.jl index 9b55d760..fe033a62 100644 --- a/src/basis/monomial.jl +++ b/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}) diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl index d1b2b609..07aa59f6 100644 --- a/src/basis/polyharmonic_spline.jl +++ b/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(∇²ℒ) diff --git a/src/operators/monomial/monomial.jl b/src/operators/monomial/monomial.jl index 10685899..616daf98 100644 --- a/src/operators/monomial/monomial.jl +++ b/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 @@ end 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 diff --git a/src/operators/virtual.jl b/src/operators/virtual.jl index 16c8afde..b704c2f4 100644 --- a/src/operators/virtual.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 3a8886f0..60e0bdfd 100644 --- a/src/utils.jl +++ b/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) diff --git a/test/basis/basis.jl b/test/basis/basis.jl new file mode 100644 index 00000000..9f93aba5 --- /dev/null +++ b/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 diff --git a/test/basis/monomial.jl b/test/basis/monomial.jl index 568ab253..b6a356b4 100644 --- a/test/basis/monomial.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 759eb83c..95e446d3 100644 --- a/test/runtests.jl +++ b/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