Feature/clean up

src/RadialBasisFunctions.jl

@@ -16,6 +16,7 @@ export AbstractPHS, PHS, PHS1, PHS3, PHS5, PHS7
 export IMQ
 export Gaussian
 export MonomialBasis
+export degree, dim
 
 include("utils.jl")
 export find_neighbors, reorder_points!
@@ -42,7 +43,7 @@ export ∂virtual
 
 include("operators/monomial/monomial.jl")
 
-include("operators/operator_combinations.jl")
+include("operators/operator_algebra.jl")
 
 include("interpolation.jl")
 export Interpolator
@@ -52,7 +53,7 @@ export Regrid, regrid
 
 # Some consts and aliases
 const Δ = ∇² # some people like this notation for the Laplacian
-const AVOID_NAN = 1e-16
+const AVOID_INF = 1e-16
 
 using PrecompileTools
 @setup_workload begin

src/basis/basis.jl

@@ -29,6 +29,9 @@ function (ℒmon::ℒ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
+
 include("polyharmonic_spline.jl")
 include("inverse_multiquadric.jl")
 include("gaussian.jl")

src/basis/polyharmonic_spline.jl

@@ -39,7 +39,7 @@ end
 
 (phs::PHS1)(x, xᵢ) = euclidean(x, xᵢ)
 function ∂(::PHS1, dim::Int)
-    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + AVOID_NAN)
+    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + AVOID_INF)
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS1)
@@ -49,14 +49,14 @@ end
 function ∂²(::PHS1, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) /
-               (euclidean(x, xᵢ)^3 + AVOID_NAN)
+               (euclidean(x, xᵢ)^3 + AVOID_INF)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS1)
     function ∇²ℒ(x, xᵢ)
         return sum(
-            (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_NAN)
+            (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_INF)
         )
     end
     return ℒRadialBasisFunction(∇²ℒ)
@@ -87,14 +87,14 @@ end
 function ∂²(::PHS3, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) /
-               (euclidean(x, xᵢ) + AVOID_NAN)
+               (euclidean(x, xᵢ) + AVOID_INF)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS3)
     function ∇²ℒ(x, xᵢ)
         return sum(
-            3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + AVOID_NAN)
+            3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + AVOID_INF)
         )
     end
     return ℒRadialBasisFunction(∇²ℒ)

src/operators/directional.jl

@@ -14,12 +14,12 @@ end
 Builds a `RadialBasisOperator` where the operator is the directional derivative, `Directional`.
 """
 function directional(
-    data::AbstractVector{D},
+    data::AbstractVector,
     v::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Directional(f, v)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
@@ -45,14 +45,12 @@ end
 
 function RadialBasisOperator(
     ℒ::Directional,
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(D), Na, Nd)
+) where {T<:Int,B<:AbstractRadialBasis}
+    weights = _build_weights(ℒ, data, data, adjl, basis)
     return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
 end
 
@@ -64,12 +62,26 @@ function RadialBasisOperator(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(TD), Na, Nd)
+    weights = _build_weights(ℒ, data, eval_points, adjl, basis)
     return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end
 
+function _build_weights(ℒ::Directional, data, eval_points, adjl, basis)
+    v = ℒ.v
+    N = length(first(data))
+    @assert length(v) == N || length(v) == length(data) "wrong size for v"
+    if length(v) == N
+        return mapreduce(+, enumerate(ℒ.ℒ)) do (i, ℒ)
+            _build_weights(ℒ, data, eval_points, adjl, basis) * v[i]
+        end
+    else
+        vv = ntuple(i -> getindex.(v, i), N)
+        return mapreduce(+, enumerate(ℒ.ℒ)) do (i, ℒ)
+            Diagonal(vv[i]) * _build_weights(ℒ, data, eval_points, adjl, basis)
+        end
+    end
+end
+
 function update_weights!(op::RadialBasisOperator{<:Directional})
     v = op.ℒ.v
     N = length(first(op.data))
@@ -88,7 +100,5 @@ function update_weights!(op::RadialBasisOperator{<:Directional})
     return nothing
 end
 
-Base.size(op::RadialBasisOperator{<:Directional}) = size(op.weights)
-
 # pretty printing
 print_op(op::Directional) = "Directional Gradient (∇f⋅v)"

src/operators/gradient.jl

@@ -14,11 +14,11 @@ end
 Builds a `RadialBasisOperator` where the operator is the gradient, `Gradient`.
 """
 function gradient(
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Gradient(f)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)

src/operators/laplacian.jl

@@ -9,11 +9,11 @@ end
 
 # convienience constructors
 function laplacian(
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {T<:Int,B<:AbstractRadialBasis}
     ℒ = Laplacian(∇²)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end

src/operators/monomial/monomial.jl

@@ -15,10 +15,9 @@ end
 function ∇²(m::MonomialBasis{Dim,Deg}) where {Dim,Deg}
     ∂² = ntuple(dim -> ∂(m, 2, dim), Dim)
     function basis!(b, x)
-        cache = ones(size(b))
-        b .= 0
+        cache = ones(eltype(x), size(b))
+        b .= zero(eltype(x))
         for ∂²! in ∂²
-            # use mapreduce here instead?
             ∂²!(cache, x)
             b .+= cache
         end
@@ -34,8 +33,8 @@ end
 
 function build_monomial_basis(ids::Vector{Vector{Vector{T}}}, c::Vector{T}) where {T<:Int}
     function basis!(db::AbstractVector{B}, x::AbstractVector) where {B}
-        db .= 1
-        # TODO flatten loop - why does it allocate here
+        db .= one(eltype(x))
+        # TODO optimize - allocations
         @views @inbounds for i in eachindex(ids), j in eachindex(ids[i])
             db[ids[i][j]] *= x[i]
         end

src/operators/operator_algebra.jl

@@ -0,0 +1,45 @@
+for op in (:+, :-)
+    @eval function Base.$op(a::ℒRadialBasisFunction, b::ℒRadialBasisFunction)
+        additive_ℒRBF(x, xᵢ) = Base.$op(a(x, xᵢ), b(x, xᵢ))
+        return ℒRadialBasisFunction(additive_ℒRBF)
+    end
+end
+
+for op in (:+, :-)
+    @eval function Base.$op(
+        a::ℒMonomialBasis{Dim,Deg}, b::ℒMonomialBasis{Dim,Deg}
+    ) where {Dim,Deg}
+        function additive_ℒMon(m, x)
+            m .= Base.$op.(a(x), b(x))
+            return nothing
+        end
+        return ℒMonomialBasis(Dim, Deg, additive_ℒMon)
+    end
+end
+
+for op in (:+, :-)
+    @eval function Base.$op(op1::RadialBasisOperator, op2::RadialBasisOperator)
+        _check_compatible(op1, op2)
+        k = _update_stencil(op1, op2)
+        ℒ(x) = Base.$op(op1.ℒ(x), op2.ℒ(x))
+        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
+    end
+end
+
+function _check_compatible(op1::RadialBasisOperator, op2::RadialBasisOperator)
+    if !all(op1.data .≈ op2.data)
+        throw(
+            ArgumentError("Can not add operators that were not built with the same data.")
+        )
+    end
+    if !all(op1.adjl .≈ op2.adjl)
+        throw(ArgumentError("Can not add operators that do not have the same stencils."))
+    end
+end
+
+function _update_stencil(op1::RadialBasisOperator, op2::RadialBasisOperator)
+    k1 = length(first((op1.adjl)))
+    k2 = length(first((op2.adjl)))
+    k = k1 > k2 ? k1 : k2
+    return k
+end