From 270568291aafedd91e63ceec254c47b5f993fdd7 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Fri, 7 Mar 2025 00:10:05 -0500
Subject: [PATCH 01/13] remove unnecesary fields of operators and make the
 operators functors

---
 src/basis/polyharmonic_spline.jl   | 32 +++++++++++++++---------------
 src/operators/directional.jl       |  8 ++++----
 src/operators/laplacian.jl         | 11 ++++------
 src/operators/monomial/monomial.jl |  5 -----
 src/operators/operator_algebra.jl  | 12 +++++------
 src/operators/operators.jl         |  6 ++----
 src/operators/partial.jl           | 17 +++++++---------
 test/Project.toml                  |  1 +
 8 files changed, 40 insertions(+), 52 deletions(-)

diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index fa8af001..b1015dd2 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -40,18 +40,18 @@ end
 (phs::PHS1)(x, xᵢ) = euclidean(x, xᵢ)
 function ∂(::PHS1, dim::Int)
     ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + AVOID_INF)
-    return ℒRadialBasisFunction(∂ℒ)
+    return ∂ℒ
 end
 function ∇(::PHS1)
     ∇ℒ(x, xᵢ) = (x .- xᵢ) / euclidean(x, xᵢ)
-    return ℒRadialBasisFunction(∇ℒ)
+    return ∇ℒ
 end
 function ∂²(::PHS1, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) /
                (euclidean(x, xᵢ)^3 + AVOID_INF)
     end
-    return ℒRadialBasisFunction(∂²ℒ)
+    return ∂²ℒ
 end
 function ∇²(::PHS1)
     function ∇²ℒ(x, xᵢ)
@@ -59,7 +59,7 @@ function ∇²(::PHS1)
             (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_INF)
         )
     end
-    return ℒRadialBasisFunction(∇²ℒ)
+    return ∇²ℒ
 end
 
 """
@@ -78,18 +78,18 @@ end
 (phs::PHS3)(x, xᵢ) = euclidean(x, xᵢ)^3
 function ∂(::PHS3, dim::Int)
     ∂ℒ(x, xᵢ) = 3 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)
-    return ℒRadialBasisFunction(∂ℒ)
+    return ∂ℒ
 end
 function ∇(::PHS3)
     ∇ℒ(x, xᵢ) = 3 * (x .- xᵢ) * euclidean(x, xᵢ)
-    return ℒRadialBasisFunction(∇ℒ)
+    return ∇ℒ
 end
 function ∂²(::PHS3, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) /
                (euclidean(x, xᵢ) + AVOID_INF)
     end
-    return ℒRadialBasisFunction(∂²ℒ)
+    return ∂²ℒ
 end
 function ∇²(::PHS3)
     function ∇²ℒ(x, xᵢ)
@@ -97,7 +97,7 @@ function ∇²(::PHS3)
             3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + AVOID_INF)
         )
     end
-    return ℒRadialBasisFunction(∇²ℒ)
+    return ∇²ℒ
 end
 
 """
@@ -116,23 +116,23 @@ end
 (phs::PHS5)(x, xᵢ) = euclidean(x, xᵢ)^5
 function ∂(::PHS5, dim::Int)
     ∂ℒ(x, xᵢ) = 5 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)^3
-    return ℒRadialBasisFunction(∂ℒ)
+    return ∂ℒ
 end
 function ∇(::PHS5)
     ∇ℒ(x, xᵢ) = 5 * (x .- xᵢ) * euclidean(x, xᵢ)^3
-    return ℒRadialBasisFunction(∇ℒ)
+    return ∇ℒ
 end
 function ∂²(::PHS5, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 5 * euclidean(x, xᵢ) * (3 * (x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ))
     end
-    return ℒRadialBasisFunction(∂²ℒ)
+    return ∂²ℒ
 end
 function ∇²(::PHS5)
     function ∇²ℒ(x, xᵢ)
         return sum(5 * euclidean(x, xᵢ) * (3 * (x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)))
     end
-    return ℒRadialBasisFunction(∇²ℒ)
+    return ∇²ℒ
 end
 
 """
@@ -151,23 +151,23 @@ end
 (phs::PHS7)(x, xᵢ) = euclidean(x, xᵢ)^7
 function ∂(::PHS7, dim::Int)
     ∂ℒ(x, xᵢ) = 7 * (x[dim] - xᵢ[dim]) * euclidean(x, xᵢ)^5
-    return ℒRadialBasisFunction(∂ℒ)
+    return ∂ℒ
 end
 function ∇(::PHS7)
     ∇ℒ(x, xᵢ) = 7 * (x .- xᵢ) * euclidean(x, xᵢ)^5
-    return ℒRadialBasisFunction(∇ℒ)
+    return ∇ℒ
 end
 function ∂²(::PHS7, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 7 * euclidean(x, xᵢ)^3 * (5 * (x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ))
     end
-    return ℒRadialBasisFunction(∂²ℒ)
+    return ∂²ℒ
 end
 function ∇²(::PHS7)
     function ∇²ℒ(x, xᵢ)
         return sum(7 * euclidean(x, xᵢ)^3 * (5 * (x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)))
     end
-    return ℒRadialBasisFunction(∇²ℒ)
+    return ∇²ℒ
 end
 
 function Base.show(io::IO, rbf::R) where {R<:AbstractPHS}
diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index b36654e8..ac89db5e 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -71,12 +71,12 @@ function _build_weights(ℒ::Directional, data, eval_points, adjl, basis)
     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, ℒ)
+        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, ℒ)
+        return mapreduce(+, enumerate(ℒ)) do (i, ℒ)
             Diagonal(vv[i]) * _build_weights(ℒ, data, eval_points, adjl, basis)
         end
     end
@@ -86,12 +86,12 @@ function update_weights!(op::RadialBasisOperator{<:Directional})
     v = op.ℒ.v
     N = length(first(op.data))
     if length(v) == N
-        op.weights .= mapreduce(+, enumerate(op.ℒ.ℒ)) do (i, ℒ)
+        op.weights .= mapreduce(+, enumerate(op.ℒ)) do (i, ℒ)
             _build_weights(ℒ, op) * v[i]
         end
     else
         vv = ntuple(i -> getindex.(v, i), N)
-        op.weights .= mapreduce(+, enumerate(op.ℒ.ℒ)) do (i, ℒ)
+        op.weights .= mapreduce(+, enumerate(op.ℒ)) do (i, ℒ)
             Diagonal(vv[i]) * _build_weights(ℒ, op)
         end
     end
diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
index 137677d9..711a6d9d 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -3,9 +3,8 @@
 
 Operator for the sum of the second derivatives w.r.t. each independent variable.
 """
-struct Laplacian{L<:Function} <: ScalarValuedOperator
-    ℒ::L
-end
+struct Laplacian <: ScalarValuedOperator end
+(::Laplacian)(x) = ∇²(x)
 
 # convienience constructors
 function laplacian(
@@ -14,8 +13,7 @@ function laplacian(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
 ) where {T<:Int,B<:AbstractRadialBasis}
-    ℒ = Laplacian(∇²)
-    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
+    return RadialBasisOperator(Laplacian(), data, basis; k=k, adjl=adjl)
 end
 
 function laplacian(
@@ -25,8 +23,7 @@ function laplacian(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {T<:Int,B<:AbstractRadialBasis}
-    ℒ = Laplacian(∇²)
-    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
+    return RadialBasisOperator(Laplacian(), data, eval_points, basis; k=k, adjl=adjl)
 end
 
 # pretty printing
diff --git a/src/operators/monomial/monomial.jl b/src/operators/monomial/monomial.jl
index 616daf98..4d763e6d 100644
--- a/src/operators/monomial/monomial.jl
+++ b/src/operators/monomial/monomial.jl
@@ -1,8 +1,3 @@
-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}
diff --git a/src/operators/operator_algebra.jl b/src/operators/operator_algebra.jl
index 98d8fef7..cec96a74 100644
--- a/src/operators/operator_algebra.jl
+++ b/src/operators/operator_algebra.jl
@@ -1,9 +1,9 @@
-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::ℒ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(
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
index 13cefc3b..7dbcb931 100644
--- a/src/operators/operators.jl
+++ b/src/operators/operators.jl
@@ -2,8 +2,6 @@ abstract type AbstractOperator end
 abstract type ScalarValuedOperator <: AbstractOperator end
 abstract type VectorValuedOperator <: AbstractOperator end
 
-(op::AbstractOperator)(x) = op.ℒ(x)
-
 """
     struct RadialBasisOperator
 
@@ -131,13 +129,13 @@ end
 
 # update weights
 function update_weights!(op::RadialBasisOperator)
-    op.weights .= _build_weights(op.ℒ.ℒ, op)
+    op.weights .= _build_weights(op.ℒ, op)
     validate_cache!(op)
     return nothing
 end
 
 function update_weights!(op::RadialBasisOperator{<:VectorValuedOperator})
-    for (i, ℒ) in enumerate(op.ℒ.ℒ)
+    for (i, ℒ) in enumerate(op.ℒ)
         op.weights[i] .= _build_weights(ℒ, op)
     end
     validate_cache!(op)
diff --git a/src/operators/partial.jl b/src/operators/partial.jl
index 491d1bce..c8409cb5 100644
--- a/src/operators/partial.jl
+++ b/src/operators/partial.jl
@@ -3,12 +3,15 @@
 
 Builds an operator for a first order partial derivative.
 """
-struct Partial{L<:Function,T<:Int} <: ScalarValuedOperator
-    ℒ::L
+struct Partial{T<:Int} <: ScalarValuedOperator
     order::T
     dim::T
 end
 
+function (op::Partial)(x)
+    return ∂(x, op.order, op.dim)
+end
+
 # convienience constructors
 """
     function partial(data, order, dim, basis; k=autoselect_k(data, basis))
@@ -23,10 +26,7 @@ function partial(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
 ) where {T<:Int,B<:AbstractRadialBasis}
-    f = let o = order, dim = dim
-        x -> ∂(x, o, dim)
-    end
-    ℒ = Partial(f, order, dim)
+    ℒ = Partial(order, dim)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
@@ -44,10 +44,7 @@ function partial(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {T<:Int,B<:AbstractRadialBasis}
-    f = let o = order, dim = dim
-        x -> ∂(x, o, dim)
-    end
-    ℒ = Partial(f, order, dim)
+    ℒ = Partial(order, dim)
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end
 
diff --git a/test/Project.toml b/test/Project.toml
index 4b975f5d..db9f57bb 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,3 +6,4 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

From 089f9decef27e44124d7df5ac8c8370c98f0cf3a Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sat, 8 Mar 2025 08:59:18 -0500
Subject: [PATCH 02/13] refactor operators to remove one layer of type wrapping

---
 src/RadialBasisFunctions.jl |  8 +++----
 src/operators/gradient.jl   | 14 ++++++-----
 src/operators/laplacian.jl  |  2 +-
 src/operators/operators.jl  | 34 ++++----------------------
 src/operators/partial.jl    |  5 +---
 src/solve.jl                | 48 +++++++++++++++++++++++++++++++------
 test/operators/gradient.jl  |  2 +-
 test/operators/partial.jl   |  2 +-
 8 files changed, 62 insertions(+), 53 deletions(-)

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index 84c2fa03..562d895d 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -21,12 +21,12 @@ export degree, dim
 include("utils.jl")
 export find_neighbors, reorder_points!
 
-include("solve.jl")
-
 include("operators/operators.jl")
 export RadialBasisOperator, ScalarValuedOperator, VectorValuedOperator
 export update_weights!, is_cache_valid
 
+include("solve.jl")
+
 include("operators/partial.jl")
 export Partial, partial
 
@@ -78,8 +78,8 @@ using PrecompileTools
             ∂x = ∂(z)
 
             # gradient
-            ∇ = gradient(x, b)
-            ∇z = ∇(z)
+            #∇ = gradient(x, b)
+            #∇z = ∇(z)
 
             # laplacian
             ∇² = laplacian(x, b)
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index 9a0d9a39..b3e7bd93 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -3,8 +3,10 @@
 
 Builds an operator for the gradient of a function.
 """
-struct Gradient{L<:NTuple} <: VectorValuedOperator
-    ℒ::L
+struct Gradient{Dim} <: VectorValuedOperator{Dim} end
+
+function (op::Gradient{Dim})(basis) where {Dim}
+    return ntuple(dim -> ∂(basis, dim), Dim)
 end
 
 # convienience constructors
@@ -19,8 +21,8 @@ function gradient(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
-    ℒ = Gradient(f)
+    Dim = length(first(data))
+    ℒ = Gradient{Dim}()
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
@@ -36,8 +38,8 @@ function gradient(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
-    ℒ = Gradient(f)
+    Dim = length(first(data))
+    ℒ = Gradient{Dim}()
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end
 
diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
index 711a6d9d..110b5967 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -4,7 +4,7 @@
 Operator for the sum of the second derivatives w.r.t. each independent variable.
 """
 struct Laplacian <: ScalarValuedOperator end
-(::Laplacian)(x) = ∇²(x)
+(::Laplacian)(basis) = ∇²(basis)
 
 # convienience constructors
 function laplacian(
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
index 7dbcb931..2a5c3898 100644
--- a/src/operators/operators.jl
+++ b/src/operators/operators.jl
@@ -1,6 +1,6 @@
 abstract type AbstractOperator end
 abstract type ScalarValuedOperator <: AbstractOperator end
-abstract type VectorValuedOperator <: AbstractOperator end
+abstract type VectorValuedOperator{Dim} <: AbstractOperator end
 
 """
     struct RadialBasisOperator
@@ -54,31 +54,6 @@ function RadialBasisOperator(
     return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis, true)
 end
 
-function RadialBasisOperator(
-    ℒ::VectorValuedOperator,
-    data::AbstractVector,
-    basis::B=PHS(3; poly_deg=2);
-    k::T=autoselect_k(data, basis),
-    adjl=find_neighbors(data, k),
-) where {T<:Int,B<:AbstractRadialBasis}
-    weights = ntuple(i -> _build_weights(ℒ.ℒ[i], data, data, adjl, basis), length(ℒ.ℒ))
-    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis, true)
-end
-
-function RadialBasisOperator(
-    ℒ::VectorValuedOperator,
-    data::AbstractVector{TD},
-    eval_points::AbstractVector{TE},
-    basis::B=PHS(3; poly_deg=2);
-    k::T=autoselect_k(data, basis),
-    adjl=find_neighbors(data, eval_points, k),
-) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    weights = ntuple(length(ℒ.ℒ)) do i
-        _build_weights(ℒ.ℒ[i], data, eval_points, adjl, basis)
-    end
-    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis, true)
-end
-
 dim(op::RadialBasisOperator) = length(first(op.data))
 
 # caching
@@ -134,9 +109,10 @@ function update_weights!(op::RadialBasisOperator)
     return nothing
 end
 
-function update_weights!(op::RadialBasisOperator{<:VectorValuedOperator})
-    for (i, ℒ) in enumerate(op.ℒ)
-        op.weights[i] .= _build_weights(ℒ, op)
+function update_weights!(op::RadialBasisOperator{<:VectorValuedOperator{Dim}}) where {Dim}
+    new_weights = _build_weights(op.ℒ, op)
+    for i in 1:Dim
+        op.weights[i] .= new_weights[i]
     end
     validate_cache!(op)
     return nothing
diff --git a/src/operators/partial.jl b/src/operators/partial.jl
index c8409cb5..72b3c12d 100644
--- a/src/operators/partial.jl
+++ b/src/operators/partial.jl
@@ -7,10 +7,7 @@ struct Partial{T<:Int} <: ScalarValuedOperator
     order::T
     dim::T
 end
-
-function (op::Partial)(x)
-    return ∂(x, op.order, op.dim)
-end
+(op::Partial)(basis) = ∂(basis, op.order, op.dim)
 
 # convienience constructors
 """
diff --git a/src/solve.jl b/src/solve.jl
index 39402976..3dbdc20c 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -29,12 +29,12 @@ function _build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon)
     Na = length(adjl)
     I = zeros(Int, k * Na)
     J = reduce(vcat, adjl)
-    V = zeros(TD, k * Na)
+    V = zeros(TD, k * Na, _num_ops(ℒrbf))
 
     # create work arrays
     n = sum(sizes)
     A = Symmetric[Symmetric(zeros(TD, n, n), :U) for _ in 1:nchunks]
-    b = Vector[zeros(TD, n) for _ in 1:nchunks]
+    b = [_prepare_b(ℒrbf, TD, n) for _ in 1:nchunks]
     d = Vector{Vector{eltype(data)}}(undef, nchunks)
 
     # build stencil for each data point and store in global weight matrix
@@ -42,18 +42,29 @@ function _build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon)
         for i in xrange
             I[((i - 1) * k + 1):(i * k)] .= i
             d[ichunk] = data[adjl[i]]
-            V[((i - 1) * k + 1):(i * k)] = @views _build_stencil!(
+            V[((i - 1) * k + 1):(i * k), :] = @views _build_stencil!(
                 A[ichunk], b[ichunk], ℒrbf, ℒmon, d[ichunk], eval_points[i], basis, mon, k
             )
         end
     end
 
-    return sparse(I, J, V, length(adjl), length(data))
+    nrows = length(adjl)
+    ncols = length(data)
+    if size(V, 2) == 1
+        return sparse(I, J, V[:, 1], nrows, ncols)
+    else
+        return ntuple(i -> sparse(I, J, V[:, i], nrows, ncols), size(V, 2))
+    end
 end
 
+_num_ops(_) = 1
+_num_ops(ℒ::Tuple) = length(ℒ)
+_prepare_b(_, T, n) = zeros(T, n)
+_prepare_b(ℒ::Tuple, T, n) = zeros(T, n, length(ℒ))
+
 function _build_stencil!(
     A::Symmetric,
-    b::Vector,
+    b,
     ℒrbf,
     ℒmon,
     data::AbstractVector{TD},
@@ -64,7 +75,7 @@ function _build_stencil!(
 ) where {TD,TE,B<:AbstractRadialBasis}
     _build_collocation_matrix!(A, data, basis, mon, k)
     _build_rhs!(b, ℒrbf, ℒmon, data, eval_point, basis, k)
-    return (A \ b)[1:k]
+    return (A \ b)[1:k, :]
 end
 
 function _build_collocation_matrix!(
@@ -89,7 +100,7 @@ function _build_collocation_matrix!(
 end
 
 function _build_rhs!(
-    b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{TD}, eval_point::TE, basis::B, k::K
+    b, ℒrbf, ℒmon, data::AbstractVector{TD}, eval_point::TE, basis::B, k::K
 ) where {TD,TE,B<:AbstractRadialBasis,K<:Int}
     # radial basis section
     @inbounds for i in eachindex(data)
@@ -105,3 +116,26 @@ function _build_rhs!(
 
     return nothing
 end
+
+function _build_rhs!(
+    b, ℒrbf::Tuple, ℒmon::Tuple, data::AbstractVector{TD}, eval_point::TE, basis::B, k::K
+) where {TD,TE,B<:AbstractRadialBasis,K<:Int}
+    @assert size(b, 2) == length(ℒrbf) == length(ℒmon) "b, ℒrbf, ℒmon must have the same length"
+    # radial basis section
+    for (j, ℒ) in enumerate(ℒrbf)
+        @inbounds for i in eachindex(data)
+            b[i, j] = ℒ(eval_point, data[i])
+        end
+    end
+
+    # monomial augmentation
+    if basis.poly_deg > -1
+        N = size(b, 1)
+        for (j, ℒ) in enumerate(ℒmon)
+            bmono = view(b, (k + 1):N, j)
+            ℒ(bmono, eval_point)
+        end
+    end
+
+    return nothing
+end
diff --git a/test/operators/gradient.jl b/test/operators/gradient.jl
index e3e5d436..8e21b071 100644
--- a/test/operators/gradient.jl
+++ b/test/operators/gradient.jl
@@ -32,6 +32,6 @@ end
 end
 
 @testset "Printing" begin
-    ∇ = Gradient((1, 2))
+    ∇ = Gradient{2}()
     @test RadialBasisFunctions.print_op(∇) == "Gradient (∇f)"
 end
diff --git a/test/operators/partial.jl b/test/operators/partial.jl
index ad3ae6b4..eedaedde 100644
--- a/test/operators/partial.jl
+++ b/test/operators/partial.jl
@@ -75,6 +75,6 @@ end
 end
 
 @testset "Printing" begin
-    ∂ = Partial(identity, 1, 2)
+    ∂ = Partial(1, 2)
     @test RadialBasisFunctions.print_op(∂) == "∂ⁿf/∂xᵢ (n = 1, i = 2)"
 end

From fb70e6a421f182382a4cc5090b1efabbe991b133 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sat, 8 Mar 2025 09:11:32 -0500
Subject: [PATCH 03/13] uncomment precompile for Gradient

---
 src/RadialBasisFunctions.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index 562d895d..d71ba6d7 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -78,8 +78,8 @@ using PrecompileTools
             ∂x = ∂(z)
 
             # gradient
-            #∇ = gradient(x, b)
-            #∇z = ∇(z)
+            ∇ = gradient(x, b)
+            ∇z = ∇(z)
 
             # laplacian
             ∇² = laplacian(x, b)

From 7df48da41c8082aff1c10fe31c3bc03ee6308a2e Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Mon, 10 Mar 2025 21:49:28 -0400
Subject: [PATCH 04/13] fix Directional, add benchmark and formatting to
 workflows

---
 .github/workflows/AirspeedVelocity.yml | 78 +++++++++++++++++++++++++
 .github/workflows/FormatPR.yml         | 11 ++++
 benchmark/Project.toml                 |  8 +++
 benchmark/benchmarks.jl                |  4 ++
 src/operators/directional.jl           | 79 +++++++++-----------------
 test/operators/directional.jl          |  2 +-
 6 files changed, 128 insertions(+), 54 deletions(-)
 create mode 100644 .github/workflows/AirspeedVelocity.yml
 create mode 100644 .github/workflows/FormatPR.yml
 create mode 100644 benchmark/Project.toml
 create mode 100644 benchmark/benchmarks.jl

diff --git a/.github/workflows/AirspeedVelocity.yml b/.github/workflows/AirspeedVelocity.yml
new file mode 100644
index 00000000..81c3af42
--- /dev/null
+++ b/.github/workflows/AirspeedVelocity.yml
@@ -0,0 +1,78 @@
+name: Benchmark a Pull Request
+
+on:
+  pull_request:
+    branches:
+      - main
+
+permissions:
+  pull-requests: write
+
+jobs:
+    benchmark:
+        runs-on: ubuntu-latest
+
+        steps:
+            - uses: actions/checkout@v4
+            - uses: julia-actions/setup-julia@v2
+              with:
+                version: "1.10"
+            - uses: julia-actions/cache@v2
+            - name: Extract Package Name from Project.toml
+              id: extract-package-name
+              run: |
+                PACKAGE_NAME=$(grep "^name" Project.toml | sed 's/^name = "\(.*\)"$/\1/')
+                echo "package_name=$PACKAGE_NAME" >> $GITHUB_OUTPUT
+            - name: Build AirspeedVelocity
+              env:
+                JULIA_NUM_THREADS: 2
+              run: |
+                # Lightweight build step, as sometimes the runner runs out of memory:
+                julia -e 'ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0; import Pkg; Pkg.add(;url="https://github.com/MilesCranmer/AirspeedVelocity.jl.git")'
+                julia -e 'ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0; import Pkg; Pkg.build("AirspeedVelocity")'
+            - name: Add ~/.julia/bin to PATH
+              run: |
+                echo "$HOME/.julia/bin" >> $GITHUB_PATH
+            - name: Run benchmarks
+              run: |
+                echo $PATH
+                ls -l ~/.julia/bin
+                mkdir results
+                benchpkg ${{ steps.extract-package-name.outputs.package_name }} --rev="${{github.event.repository.default_branch}},${{github.event.pull_request.head.sha}}" --url=${{ github.event.repository.clone_url }} --bench-on="${{github.event.pull_request.head.sha}}" --output-dir=results/ --tune
+            - name: Create plots from benchmarks
+              run: |
+                mkdir -p plots
+                benchpkgplot ${{ steps.extract-package-name.outputs.package_name }} --rev="${{github.event.repository.default_branch}},${{github.event.pull_request.head.sha}}" --npart=10 --format=png --input-dir=results/ --output-dir=plots/
+            - name: Upload plot as artifact
+              uses: actions/upload-artifact@v4
+              with:
+                name: plots
+                path: plots
+            - name: Create markdown table from benchmarks
+              run: |
+                benchpkgtable ${{ steps.extract-package-name.outputs.package_name }} --rev="${{github.event.repository.default_branch}},${{github.event.pull_request.head.sha}}" --input-dir=results/ --ratio > table.md
+                echo '### Benchmark Results' > body.md
+                echo '' >> body.md
+                echo '' >> body.md
+                cat table.md >> body.md
+                echo '' >> body.md
+                echo '' >> body.md
+                echo '### Benchmark Plots' >> body.md
+                echo 'A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.' >> body.md
+                echo 'Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).' >> body.md
+
+            - name: Find Comment
+              uses: peter-evans/find-comment@v3
+              id: fcbenchmark
+              with:
+                issue-number: ${{ github.event.pull_request.number }}
+                comment-author: 'github-actions[bot]'
+                body-includes: Benchmark Results
+
+            - name: Comment on PR
+              uses: peter-evans/create-or-update-comment@v4
+              with:
+                comment-id: ${{ steps.fcbenchmark.outputs.comment-id }}
+                issue-number: ${{ github.event.pull_request.number }}
+                body-path: body.md
+                edit-mode: replace
\ No newline at end of file
diff --git a/.github/workflows/FormatPR.yml b/.github/workflows/FormatPR.yml
new file mode 100644
index 00000000..8ce8a2d4
--- /dev/null
+++ b/.github/workflows/FormatPR.yml
@@ -0,0 +1,11 @@
+name: Format suggestions
+on:
+  pull_request
+
+jobs:
+  code-style:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: julia-actions/julia-format@v3
+        with:
+          version: '1' # Set `version` to '1.0.54' if you need to use JuliaFormatter.jl v1.0.54 (default: '1')
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 00000000..9e70f370
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+RadialBasisFunctions = "79ee0514-adf7-4479-8807-6f72ea8967e8"
+
+[compat]
+BenchmarkTools = "1.5"
+RadialBasisFunctions = "0.2.5"
+julia = "1.9"
\ No newline at end of file
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
new file mode 100644
index 00000000..aa278e4f
--- /dev/null
+++ b/benchmark/benchmarks.jl
@@ -0,0 +1,4 @@
+using BenchmarkTools
+using RadialBasisFunctions
+
+const SUITE = BenchmarkGroup()
diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index ac89db5e..ec36dbb5 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -3,10 +3,14 @@
 
 Operator for the directional derivative, or the inner product of the gradient and a direction vector.
 """
-struct Directional{L<:NTuple,T} <: ScalarValuedOperator
-    ℒ::L
+struct Directional{Dim,T} <: ScalarValuedOperator
     v::T
 end
+Directional{Dim}(v) where {Dim} = Directional{Dim,typeof(v)}(v)
+
+function (op::Directional{Dim})(basis) where {Dim}
+    return ntuple(dim -> ∂(basis, dim), Dim)
+end
 
 """
     function directional(data, v, basis; k=autoselect_k(data, basis))
@@ -20,8 +24,8 @@ function directional(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
-    ℒ = Directional(f, v)
+    Dim = length(first(data))
+    ℒ = Directional{Dim}(v)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
@@ -38,65 +42,34 @@ function directional(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
-    f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
-    ℒ = Directional(f, v)
+    Dim = length(first(data))
+    ℒ = Directional{Dim}(v)
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end
 
-function RadialBasisOperator(
-    ℒ::Directional,
-    data::AbstractVector,
-    basis::B=PHS(3; poly_deg=2);
-    k::T=autoselect_k(data, basis),
-    adjl=find_neighbors(data, k),
-) where {T<:Int,B<:AbstractRadialBasis}
-    weights = _build_weights(ℒ, data, data, adjl, basis)
-    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
-end
-
-function RadialBasisOperator(
-    ℒ::Directional,
-    data::AbstractVector{TD},
-    eval_points::AbstractVector{TE},
-    basis::B=PHS(3; poly_deg=2);
-    k::T=autoselect_k(data, basis),
-    adjl=find_neighbors(data, eval_points, k),
-) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    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)
+function _build_weights(ℒ::Directional{Dim}, data, eval_points, adjl, basis) where {Dim}
     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
+    if !(length(v) == Dim || length(v) == length(data))
+        throw(
+            DomainError(
+                "The direction vector for Directional() should match either the dimension of the input or the number of input points. The direction vector length is $(length(v)) while there are $(length(data)) points with a dimension of $Dim",
+            ),
+        )
     end
-end
+    weights = _build_weights(Gradient{Dim}(), data, eval_points, adjl, basis)
 
-function update_weights!(op::RadialBasisOperator{<:Directional})
-    v = op.ℒ.v
-    N = length(first(op.data))
-    if length(v) == N
-        op.weights .= mapreduce(+, enumerate(op.ℒ)) do (i, ℒ)
-            _build_weights(ℒ, op) * v[i]
+    if length(v) == Dim
+        return mapreduce(+, zip(weights, v)) do zipped
+            w, vᵢ = zipped
+            w * vᵢ
         end
     else
-        vv = ntuple(i -> getindex.(v, i), N)
-        op.weights .= mapreduce(+, enumerate(op.ℒ)) do (i, ℒ)
-            Diagonal(vv[i]) * _build_weights(ℒ, op)
+        vv = ntuple(i -> getindex.(v, i), Dim)
+        return mapreduce(+, zip(weights, vv)) do zipped
+            w, vᵢ = zipped
+            Diagonal(vᵢ) * w
         end
     end
-    validate_cache!(op)
-    return nothing
 end
 
 # pretty printing
diff --git a/test/operators/directional.jl b/test/operators/directional.jl
index 7db218a2..abe81623 100644
--- a/test/operators/directional.jl
+++ b/test/operators/directional.jl
@@ -47,6 +47,6 @@ end
 end
 
 @testset "Printing" begin
-    ∇v = Directional((1, 2), (3, 4))
+    ∇v = Directional{2}(SVector(1.0, 1.0))
     @test RadialBasisFunctions.print_op(∇v) == "Directional Derivative (∇f⋅v)"
 end

From 0dd5c7fb144198ec5010ecccbdf2ffa093b37d2f Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Mon, 10 Mar 2025 22:11:32 -0400
Subject: [PATCH 05/13] initial changes to operator algebra

---
 src/RadialBasisFunctions.jl       |  3 +++
 src/operators/custom.jl           |  9 +++++++++
 src/operators/operator_algebra.jl | 15 ++++++++-------
 test/operators/laplacian.jl       |  2 +-
 4 files changed, 21 insertions(+), 8 deletions(-)
 create mode 100644 src/operators/custom.jl

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index d71ba6d7..e741f691 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -27,6 +27,9 @@ export update_weights!, is_cache_valid
 
 include("solve.jl")
 
+include("operators/custom.jl")
+export Custom
+
 include("operators/partial.jl")
 export Partial, partial
 
diff --git a/src/operators/custom.jl b/src/operators/custom.jl
new file mode 100644
index 00000000..f07a4c1e
--- /dev/null
+++ b/src/operators/custom.jl
@@ -0,0 +1,9 @@
+"""
+    Custom <: ScalarValuedOperator
+
+Builds an operator for a first order partial derivative.
+"""
+struct Custom{F<:Function} <: AbstractOperator
+    ℒ::F
+end
+(op::Custom)(basis) = op.ℒ(basis)
diff --git a/src/operators/operator_algebra.jl b/src/operators/operator_algebra.jl
index cec96a74..721a1e99 100644
--- a/src/operators/operator_algebra.jl
+++ b/src/operators/operator_algebra.jl
@@ -1,9 +1,9 @@
-#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(op1::AbstractOperator, op2::AbstractOperator)
+        additive_ℒ(basis) = Base.$op(op1(basis), op2(basis))
+        return Custom(additive_ℒ)
+    end
+end
 
 for op in (:+, :-)
     @eval function Base.$op(
@@ -21,7 +21,8 @@ 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))
+        ℒ = Base.$op(op1.ℒ, op2.ℒ)
+        @show ℒ
         return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
     end
 end
diff --git a/test/operators/laplacian.jl b/test/operators/laplacian.jl
index 606dbbbe..8c733e28 100644
--- a/test/operators/laplacian.jl
+++ b/test/operators/laplacian.jl
@@ -33,6 +33,6 @@ end
 end
 
 @testset "Printing" begin
-    ∇ = Laplacian(identity)
+    ∇ = Laplacian()
     @test RadialBasisFunctions.print_op(∇) == "Laplacian (∇²f)"
 end

From 3e2f4be2eaf364896270240dc0623827e08980d8 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 18 Mar 2025 23:34:57 -0400
Subject: [PATCH 06/13] fix operator algebra and add benchmarks

---
 benchmark/Project.toml             |  4 +-
 benchmark/benchmarks.jl            | 66 ++++++++++++++++++++++++++++++
 src/RadialBasisFunctions.jl        |  2 +-
 src/basis/monomial.jl              |  4 +-
 src/basis/polyharmonic_spline.jl   |  6 +++
 src/operators/custom.jl            |  3 ++
 src/operators/operator_algebra.jl  | 29 ++++++++-----
 test/basis/monomial.jl             |  5 +++
 test/basis/polyharmonic_spline.jl  | 20 +++++++++
 test/operators/custom.jl           | 11 +++++
 test/operators/operator_algebra.jl |  3 ++
 test/runtests.jl                   |  4 ++
 12 files changed, 143 insertions(+), 14 deletions(-)
 create mode 100644 test/operators/custom.jl

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 9e70f370..2cc4e33d 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,8 +1,10 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+HaltonSequences = "13907d55-377f-55d6-a9d6-25ac19e11b95"
 RadialBasisFunctions = "79ee0514-adf7-4479-8807-6f72ea8967e8"
+StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
 BenchmarkTools = "1.5"
 RadialBasisFunctions = "0.2.5"
-julia = "1.9"
\ No newline at end of file
+julia = "1.9"
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index aa278e4f..39fd1d49 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -1,4 +1,70 @@
 using BenchmarkTools
 using RadialBasisFunctions
+using StaticArraysCore
+using HaltonSequences
+
+f(x) = 1 + sin(4 * x[1]) + cos(3 * x[1]) + sin(2 * x[2])
+N = 100_000
+x = SVector{3}.(HaltonPoint(3)[1:N])
+y = f.(x)
 
 const SUITE = BenchmarkGroup()
+
+basis = PHS(3; poly_deg=2)
+
+∂x = partial(x, 1, 1, basis)
+SUITE["Partial"] = let s = BenchmarkGroup()
+    s["build weights"] = @benchmarkable update_weights!($∂x)
+    s["eval"] = @benchmarkable ∂x($y)
+end
+
+grad = gradient(x, basis)
+SUITE["Gradient"] = let s = BenchmarkGroup()
+    s["build weights"] = @benchmarkable update_weights!($grad)
+    s["eval"] = @benchmarkable grad($y)
+end
+
+v = SVector(2.0, 1.0, 0.5)
+v /= norm(v)
+∇v = directional(x, v, basis)
+SUITE["Directional"] = let s = BenchmarkGroup()
+    s["build weights"] = @benchmarkable update_weights!($∇v)
+    s["eval"] = @benchmarkable ∇v($y)
+end
+
+v = map(1:length(x)) do i
+    v = SVector{3}(rand(3))
+    return v /= norm(v)
+end
+∇v = directional(x, v, basis)
+SUITE["Directional (per point)"] = let s = BenchmarkGroup()
+    s["build weights"] = @benchmarkable update_weights!($∇v)
+    s["eval"] = @benchmarkable ∇v($y)
+end
+
+x1 = SVector{3}(rand(3))
+x2 = SVector{3}(rand(3))
+for basis in (IMQ, Gaussian, PHS1, PHS3, PHS5, PHS7)
+    for poly_deg in 0:2
+        b = basis(; poly_deg=poly_deg)
+        SUITE["RBF"]["$basis"]["$poly_deg"]["∂"] = @benchmarkable rbf($x1, $x2) setup = (
+            rbf = RadialBasisFunctions.∂($b, 1)
+        )
+        SUITE["RBF"]["$basis"]["$poly_deg"]["∂²"] = @benchmarkable rbf($x1, $x2) setup = (
+            rbf = RadialBasisFunctions.∂²($b, 1)
+        )
+        SUITE["RBF"]["$basis"]["$poly_deg"]["∇"] = @benchmarkable rbf($x1, $x2) setup = (
+            rbf = RadialBasisFunctions.∇($b)
+        )
+        SUITE["RBF"]["$basis"]["$poly_deg"]["∇²"] = @benchmarkable rbf($x1, $x2) setup = (
+            rbf = RadialBasisFunctions.∇²($b)
+        )
+    end
+end
+
+for dim in 1:3, deg in 0:2
+    b = zeros(binomial(dim + deg, dim))
+    SUITE["MonomialBasis"]["dim=$dim"]["deg=$deg"] = @benchmarkable mon($b, $x1) setup = (
+        mon = MonomialBasis($dim, $deg)
+    )
+end
diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index e741f691..02316e38 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -28,7 +28,7 @@ export update_weights!, is_cache_valid
 include("solve.jl")
 
 include("operators/custom.jl")
-export Custom
+export Custom, custom
 
 include("operators/partial.jl")
 export Partial, partial
diff --git a/src/basis/monomial.jl b/src/basis/monomial.jl
index fe033a62..f6802ebf 100644
--- a/src/basis/monomial.jl
+++ b/src/basis/monomial.jl
@@ -36,7 +36,7 @@ for Dim in (:1, :2, :3)
 end
 
 function _get_monomial_basis(::Val{1}, ::Val{N}) where {N}
-    return function basis!(b, x)
+    function basis!(b, x)
         b[1] = one(_get_underlying_type(x))
         if N > 0
             for n in 1:N
@@ -45,6 +45,8 @@ function _get_monomial_basis(::Val{1}, ::Val{N}) where {N}
         end
         return nothing
     end
+    basis!(b, x::AbstractVector) = basis!(b, x[1])
+    return basis!
 end
 
 function _get_monomial_basis(::Val{2}, ::Val{1})
diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index b1015dd2..0befe1aa 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -24,6 +24,12 @@ function PHS(n::T=3; poly_deg::T=2) where {T<:Int}
     return PHS7(poly_deg)
 end
 
+for phs in (:PHS1, :PHS3, :PHS5, :PHS7)
+    @eval function $phs(; poly_deg::Int=2)
+        return $phs(poly_deg)
+    end
+end
+
 """
     struct PHS1{T<:Int} <: AbstractPHS
 
diff --git a/src/operators/custom.jl b/src/operators/custom.jl
index f07a4c1e..b09c2aa4 100644
--- a/src/operators/custom.jl
+++ b/src/operators/custom.jl
@@ -7,3 +7,6 @@ struct Custom{F<:Function} <: AbstractOperator
     ℒ::F
 end
 (op::Custom)(basis) = op.ℒ(basis)
+
+# pretty printing
+print_op(op::Custom) = "Custom Operator"
diff --git a/src/operators/operator_algebra.jl b/src/operators/operator_algebra.jl
index 721a1e99..00c58bee 100644
--- a/src/operators/operator_algebra.jl
+++ b/src/operators/operator_algebra.jl
@@ -1,6 +1,23 @@
+for op in (:+, :-)
+    @eval function Base.$op(op1::RadialBasisOperator, op2::RadialBasisOperator)
+        _check_compatible(op1, op2)
+        k = _update_stencil(op1, op2)
+        ℒ = Base.$op(op1.ℒ, op2.ℒ)
+        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
+    end
+end
+
 for op in (:+, :-)
     @eval function Base.$op(op1::AbstractOperator, op2::AbstractOperator)
-        additive_ℒ(basis) = Base.$op(op1(basis), op2(basis))
+        function additive_ℒ(basis)
+            return additive_ℒrbf(x1, x2) = Base.$op(op1(basis)(x1, x2), op2(basis)(x1, x2))
+        end
+        function additive_ℒ(basis::MonomialBasis)
+            function additive_ℒMon(b, x)
+                b .= Base.$op(op1(basis)(x), op2(basis)(x))
+                return nothing
+            end
+        end
         return Custom(additive_ℒ)
     end
 end
@@ -17,16 +34,6 @@ for op in (:+, :-)
     end
 end
 
-for op in (:+, :-)
-    @eval function Base.$op(op1::RadialBasisOperator, op2::RadialBasisOperator)
-        _check_compatible(op1, op2)
-        k = _update_stencil(op1, op2)
-        ℒ = Base.$op(op1.ℒ, op2.ℒ)
-        @show ℒ
-        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(
diff --git a/test/basis/monomial.jl b/test/basis/monomial.jl
index fab7092f..a22b1cce 100644
--- a/test/basis/monomial.jl
+++ b/test/basis/monomial.jl
@@ -36,6 +36,11 @@ end
         # standard evaluation
         @test all(isapprox.(m(x), [1, 2]))
 
+        # in-place evaluation
+        b = ones(_get_underlying_type(x), 2)
+        m(b, x)
+        @test all(isapprox.(b, [1, 2]))
+
         # derivatives
         @test all(isapprox.(RBF.∂(m, 1)(x), [0, 1]))
         @test all(isapprox.(RBF.∂²(m, 1)(x), [0, 0]))
diff --git a/test/basis/polyharmonic_spline.jl b/test/basis/polyharmonic_spline.jl
index a621d0e3..ded95a4f 100644
--- a/test/basis/polyharmonic_spline.jl
+++ b/test/basis/polyharmonic_spline.jl
@@ -22,6 +22,11 @@ end
     x₁ = SVector(1.0, 2)
     x₂ = SVector(2.0, 4)
     phs = PHS(1; poly_deg=-1)
+    @test phs isa PHS1
+    @test phs.poly_deg == -1
+    phs = PHS1(; poly_deg=-1)
+    @test phs isa PHS1
+    @test phs.poly_deg == -1
 
     @testset "Distances" begin
         @test phs(x₁, x₂) ≈ sqrt((x₁[1] - x₂[1])^2 + (x₁[2] - x₂[2])^2)^1
@@ -43,6 +48,11 @@ end
     x₁ = SVector(1.0, 2)
     x₂ = SVector(2.0, 4)
     phs = PHS(3; poly_deg=-1)
+    @test phs isa PHS3
+    @test phs.poly_deg == -1
+    phs = PHS3(; poly_deg=-1)
+    @test phs isa PHS3
+    @test phs.poly_deg == -1
 
     @testset "Distances" begin
         @test phs(x₁, x₂) ≈ sqrt((x₁[1] - x₂[1])^2 + (x₁[2] - x₂[2])^2)^3
@@ -64,6 +74,11 @@ end
     x₁ = SVector(1.0, 2)
     x₂ = SVector(2.0, 4)
     phs = PHS(5; poly_deg=-1)
+    @test phs isa PHS5
+    @test phs.poly_deg == -1
+    phs = PHS5(; poly_deg=-1)
+    @test phs isa PHS5
+    @test phs.poly_deg == -1
 
     @testset "Distances" begin
         @test phs(x₁, x₂) ≈ sqrt((x₁[1] - x₂[1])^2 + (x₁[2] - x₂[2])^2)^5
@@ -85,6 +100,11 @@ end
     x₁ = SVector(1.0, 2)
     x₂ = SVector(2.0, 4)
     phs = PHS(7; poly_deg=-1)
+    @test phs isa PHS7
+    @test phs.poly_deg == -1
+    phs = PHS7(; poly_deg=-1)
+    @test phs isa PHS7
+    @test phs.poly_deg == -1
     @testset "Distances" begin
         @test phs(x₁, x₂) ≈ sqrt((x₁[1] - x₂[1])^2 + (x₁[2] - x₂[2])^2)^7
     end
diff --git a/test/operators/custom.jl b/test/operators/custom.jl
new file mode 100644
index 00000000..89f4c6f8
--- /dev/null
+++ b/test/operators/custom.jl
@@ -0,0 +1,11 @@
+using RadialBasisFunctions
+
+@testset "Constructors" begin
+    c = Custom(identity)
+    @test c(1) == 1
+end
+
+@testset "Printing" begin
+    c = Custom(identity)
+    @test RadialBasisFunctions.print_op(c) == "Custom Operator"
+end
diff --git a/test/operators/operator_algebra.jl b/test/operators/operator_algebra.jl
index 19227705..82a09862 100644
--- a/test/operators/operator_algebra.jl
+++ b/test/operators/operator_algebra.jl
@@ -18,6 +18,9 @@ dx = partial(x, 1, 1)
 dy = partial(x, 1, 2)
 
 dxdy = dx + dy
+
+##
+
 @test mean_percent_error(dxdy(y), df_dx.(x) .+ df_dy.(x)) < 1e-6
 
 dxdy = dx - dy
diff --git a/test/runtests.jl b/test/runtests.jl
index 6018a719..915ba979 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -40,6 +40,10 @@ end
     include("operators/laplacian.jl")
 end
 
+@safetestset "Custom" begin
+    include("operators/custom.jl")
+end
+
 @safetestset "Interpolation" begin
     include("operators/interpolation.jl")
 end

From 901a5a110e8ccb8eefd38a3ae05975bc3e232ab7 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 18 Mar 2025 23:56:05 -0400
Subject: [PATCH 07/13] fix PHS keyword constructors

---
 src/basis/polyharmonic_spline.jl | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index 0befe1aa..572f53e6 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -24,12 +24,6 @@ function PHS(n::T=3; poly_deg::T=2) where {T<:Int}
     return PHS7(poly_deg)
 end
 
-for phs in (:PHS1, :PHS3, :PHS5, :PHS7)
-    @eval function $phs(; poly_deg::Int=2)
-        return $phs(poly_deg)
-    end
-end
-
 """
     struct PHS1{T<:Int} <: AbstractPHS
 
@@ -176,6 +170,13 @@ function ∇²(::PHS7)
     return ∇²ℒ
 end
 
+# convient constructors using keyword arguments
+for phs in (:PHS1, :PHS3, :PHS5, :PHS7)
+    @eval function $phs(; poly_deg::Int=2)
+        return $phs(poly_deg)
+    end
+end
+
 function Base.show(io::IO, rbf::R) where {R<:AbstractPHS}
     print(io, print_basis(rbf))
     print(io, "\n└─Polynomial augmentation: degree $(rbf.poly_deg)")

From ca5f78adc8da504e3b9864acc3716054875adbcd Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 18 Mar 2025 23:57:15 -0400
Subject: [PATCH 08/13] fix monomial test

---
 test/basis/monomial.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/basis/monomial.jl b/test/basis/monomial.jl
index a22b1cce..98c2998d 100644
--- a/test/basis/monomial.jl
+++ b/test/basis/monomial.jl
@@ -37,7 +37,7 @@ end
         @test all(isapprox.(m(x), [1, 2]))
 
         # in-place evaluation
-        b = ones(_get_underlying_type(x), 2)
+        b = ones(2)
         m(b, x)
         @test all(isapprox.(b, [1, 2]))
 

From 9b69b6152155ecab9f8756a220528b9c05a7a407 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Wed, 19 Mar 2025 00:32:49 -0400
Subject: [PATCH 09/13] fix benchmarks and add some tests

---
 benchmark/Project.toml             |  1 +
 benchmark/benchmarks.jl            |  1 +
 src/operators/directional.jl       |  8 ++------
 src/operators/gradient.jl          |  2 +-
 src/operators/operator_algebra.jl  | 16 ++--------------
 test/operators/directional.jl      | 15 +++++++++++++++
 test/operators/operator_algebra.jl | 16 +++++++++++++---
 7 files changed, 35 insertions(+), 24 deletions(-)

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 2cc4e33d..9cf457e1 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 HaltonSequences = "13907d55-377f-55d6-a9d6-25ac19e11b95"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RadialBasisFunctions = "79ee0514-adf7-4479-8807-6f72ea8967e8"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 39fd1d49..1e55631c 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -2,6 +2,7 @@ using BenchmarkTools
 using RadialBasisFunctions
 using StaticArraysCore
 using HaltonSequences
+using LinearAlgebra
 
 f(x) = 1 + sin(4 * x[1]) + cos(3 * x[1]) + sin(2 * x[2])
 N = 100_000
diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index ec36dbb5..d572eec8 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -1,5 +1,5 @@
 """
-    Directional <: ScalarValuedOperator
+    Directional{Dim,T} <: ScalarValuedOperator
 
 Operator for the directional derivative, or the inner product of the gradient and a direction vector.
 """
@@ -8,10 +8,6 @@ struct Directional{Dim,T} <: ScalarValuedOperator
 end
 Directional{Dim}(v) where {Dim} = Directional{Dim,typeof(v)}(v)
 
-function (op::Directional{Dim})(basis) where {Dim}
-    return ntuple(dim -> ∂(basis, dim), Dim)
-end
-
 """
     function directional(data, v, basis; k=autoselect_k(data, basis))
 
@@ -52,7 +48,7 @@ function _build_weights(ℒ::Directional{Dim}, data, eval_points, adjl, basis) w
     if !(length(v) == Dim || length(v) == length(data))
         throw(
             DomainError(
-                "The direction vector for Directional() should match either the dimension of the input or the number of input points. The direction vector length is $(length(v)) while there are $(length(data)) points with a dimension of $Dim",
+                "The geometrical vector for Directional() should match either the dimension of the input or the number of input points. The geometrical vector length is $(length(v)) while there are $(length(data)) points with a dimension of $Dim",
             ),
         )
     end
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index b3e7bd93..7115c406 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -1,5 +1,5 @@
 """
-    Gradient <: VectorValuedOperator
+    Gradient{Dim} <: VectorValuedOperator
 
 Builds an operator for the gradient of a function.
 """
diff --git a/src/operators/operator_algebra.jl b/src/operators/operator_algebra.jl
index 00c58bee..c0585732 100644
--- a/src/operators/operator_algebra.jl
+++ b/src/operators/operator_algebra.jl
@@ -22,25 +22,13 @@ for op in (:+, :-)
     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
-
 function _check_compatible(op1::RadialBasisOperator, op2::RadialBasisOperator)
-    if !all(op1.data .≈ op2.data)
+    if (length(op1.data) != length(op2.data)) || !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)
+    if op1.adjl != op2.adjl
         throw(ArgumentError("Can not add operators that do not have the same stencils."))
     end
 end
diff --git a/test/operators/directional.jl b/test/operators/directional.jl
index abe81623..5ca36fc7 100644
--- a/test/operators/directional.jl
+++ b/test/operators/directional.jl
@@ -46,6 +46,21 @@ end
     @test mean_percent_error(∇v(y), exact) < 10
 end
 
+@testset "Incompatible Geometrical  Vector" begin
+    # wrong geometrical vector dimension
+    v = SVector(2.0, 1.0, 0.5)
+    v /= norm(v)
+    @test_throws DomainError directional(x, v)
+
+    # number of geometrical vectors don't match number of points when using a different
+    # geometrical vector for each point
+    v = map(1:(length(x) - 1)) do i
+        v = SVector{2}(rand(2))
+        return v /= norm(v)
+    end
+    @test_throws DomainError directional(x, v)
+end
+
 @testset "Printing" begin
     ∇v = Directional{2}(SVector(1.0, 1.0))
     @test RadialBasisFunctions.print_op(∇v) == "Directional Derivative (∇f⋅v)"
diff --git a/test/operators/operator_algebra.jl b/test/operators/operator_algebra.jl
index 82a09862..321932a7 100644
--- a/test/operators/operator_algebra.jl
+++ b/test/operators/operator_algebra.jl
@@ -18,10 +18,20 @@ dx = partial(x, 1, 1)
 dy = partial(x, 1, 2)
 
 dxdy = dx + dy
-
-##
-
 @test mean_percent_error(dxdy(y), df_dx.(x) .+ df_dy.(x)) < 1e-6
 
 dxdy = dx - dy
 @test mean_percent_error(dxdy(y), df_dx.(x) .- df_dy.(x)) < 1e-6
+
+# test compatibility with other operators
+dy = partial(x[1:500], 1, 2)
+@test_throws ArgumentError dx + dy
+
+dy = partial(x, 1, 2; adjl=dx.adjl[1:500])
+@test_throws ArgumentError dx + dy
+
+adjl = copy(dx.adjl)
+adjl[1] = dx.adjl[2]
+adjl[2] = dx.adjl[1]
+dy = partial(x, 1, 2; adjl=adjl)
+@test_throws ArgumentError dx + dy

From d2eed3b4cca05f2cdd6c955e28e146fa52458b6b Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Wed, 19 Mar 2025 00:49:22 -0400
Subject: [PATCH 10/13] update benchmark CI

---
 .github/workflows/AirspeedVelocity.yml | 2 +-
 src/basis/polyharmonic_spline.jl       | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/AirspeedVelocity.yml b/.github/workflows/AirspeedVelocity.yml
index 81c3af42..dd5953f6 100644
--- a/.github/workflows/AirspeedVelocity.yml
+++ b/.github/workflows/AirspeedVelocity.yml
@@ -16,7 +16,7 @@ jobs:
             - uses: actions/checkout@v4
             - uses: julia-actions/setup-julia@v2
               with:
-                version: "1.10"
+                version: "1.11"
             - uses: julia-actions/cache@v2
             - name: Extract Package Name from Project.toml
               id: extract-package-name
diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index 572f53e6..4abb744e 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -31,7 +31,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r``
 """
 struct PHS1{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS1(poly_deg::T) where {T<:Int}
+    function PHS1(poly_deg::T=2) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -69,7 +69,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^3``
 """
 struct PHS3{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS3(poly_deg::T) where {T<:Int}
+    function PHS3(poly_deg::T=2) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -107,7 +107,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^5``
 """
 struct PHS5{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS5(poly_deg::T) where {T<:Int}
+    function PHS5(poly_deg::T=2) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -142,7 +142,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^7``
 """
 struct PHS7{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS7(poly_deg::T) where {T<:Int}
+    function PHS7(poly_deg::T=2) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end

From 0b3c5f967cf04c1ae671529001a350b9c95e7dc6 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Wed, 19 Mar 2025 01:08:00 -0400
Subject: [PATCH 11/13] rm Zygote as test dep

---
 src/basis/polyharmonic_spline.jl | 8 ++++----
 test/Project.toml                | 1 -
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index 4abb744e..572f53e6 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -31,7 +31,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r``
 """
 struct PHS1{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS1(poly_deg::T=2) where {T<:Int}
+    function PHS1(poly_deg::T) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -69,7 +69,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^3``
 """
 struct PHS3{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS3(poly_deg::T=2) where {T<:Int}
+    function PHS3(poly_deg::T) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -107,7 +107,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^5``
 """
 struct PHS5{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS5(poly_deg::T=2) where {T<:Int}
+    function PHS5(poly_deg::T) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
@@ -142,7 +142,7 @@ Polyharmonic spline radial basis function:``ϕ(r) = r^7``
 """
 struct PHS7{T<:Int} <: AbstractPHS
     poly_deg::T
-    function PHS7(poly_deg::T=2) where {T<:Int}
+    function PHS7(poly_deg::T) where {T<:Int}
         check_poly_deg(poly_deg)
         return new{T}(poly_deg)
     end
diff --git a/test/Project.toml b/test/Project.toml
index db9f57bb..4b975f5d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,4 +6,3 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

From 4e8f48ddb7bff935a07a2049bd9979c2d46b63b4 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Wed, 19 Mar 2025 01:16:34 -0400
Subject: [PATCH 12/13] remove utils.jl tests because I haven;t written any yet

---
 test/runtests.jl | 4 ----
 test/utils.jl    | 0
 2 files changed, 4 deletions(-)
 delete mode 100644 test/utils.jl

diff --git a/test/runtests.jl b/test/runtests.jl
index 915ba979..8b3502e9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -63,7 +63,3 @@ end
 @safetestset "Stencil" begin
     include("solve.jl")
 end
-
-@safetestset "Utils" begin
-    include("utils.jl")
-end
diff --git a/test/utils.jl b/test/utils.jl
deleted file mode 100644
index e69de29b..00000000

From e0c76b928cf67c7cc2e4717d776fee4c96be7e1e Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Wed, 19 Mar 2025 21:04:52 -0400
Subject: [PATCH 13/13] refactor benchmarks to work with old API

---
 benchmark/benchmarks.jl | 39 ++++++++++++++++++++++++---------------
 1 file changed, 24 insertions(+), 15 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 1e55631c..48d4dc19 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -45,21 +45,30 @@ end
 
 x1 = SVector{3}(rand(3))
 x2 = SVector{3}(rand(3))
-for basis in (IMQ, Gaussian, PHS1, PHS3, PHS5, PHS7)
-    for poly_deg in 0:2
-        b = basis(; poly_deg=poly_deg)
-        SUITE["RBF"]["$basis"]["$poly_deg"]["∂"] = @benchmarkable rbf($x1, $x2) setup = (
-            rbf = RadialBasisFunctions.∂($b, 1)
-        )
-        SUITE["RBF"]["$basis"]["$poly_deg"]["∂²"] = @benchmarkable rbf($x1, $x2) setup = (
-            rbf = RadialBasisFunctions.∂²($b, 1)
-        )
-        SUITE["RBF"]["$basis"]["$poly_deg"]["∇"] = @benchmarkable rbf($x1, $x2) setup = (
-            rbf = RadialBasisFunctions.∇($b)
-        )
-        SUITE["RBF"]["$basis"]["$poly_deg"]["∇²"] = @benchmarkable rbf($x1, $x2) setup = (
-            rbf = RadialBasisFunctions.∇²($b)
-        )
+
+function benchmark_basis(SUITE, basis, poly_deg, x1, x2)
+    SUITE["RBF"]["$basis"]["$poly_deg"]["∂"] = @benchmarkable rbf($x1, $x2) setup = (
+        rbf = RadialBasisFunctions.∂($basis, 1)
+    )
+    SUITE["RBF"]["$basis"]["$poly_deg"]["∂²"] = @benchmarkable rbf($x1, $x2) setup = (
+        rbf = RadialBasisFunctions.∂²($basis, 1)
+    )
+    SUITE["RBF"]["$basis"]["$poly_deg"]["∇"] = @benchmarkable rbf($x1, $x2) setup = (
+        rbf = RadialBasisFunctions.∇($basis)
+    )
+    return SUITE["RBF"]["$basis"]["$poly_deg"]["∇²"] = @benchmarkable rbf($x1, $x2) setup = (
+        rbf = RadialBasisFunctions.∇²($basis)
+    )
+end
+
+for poly_deg in 0:2
+    for basis in (IMQ, Gaussian)
+        rbf = basis(; poly_deg=poly_deg)
+        benchmark_basis(SUITE, rbf, poly_deg, x1, x2)
+    end
+    for basis in (PHS1, PHS3, PHS5, PHS7)
+        rbf = basis(poly_deg)
+        benchmark_basis(SUITE, rbf, poly_deg, x1, x2)
     end
 end