Refactor of operator typing

.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.11"
+            - 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

.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')

benchmark/Project.toml

@@ -0,0 +1,11 @@
+[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"
+
+[compat]
+BenchmarkTools = "1.5"
+RadialBasisFunctions = "0.2.5"
+julia = "1.9"

benchmark/benchmarks.jl

@@ -0,0 +1,80 @@
+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
+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))
+
+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
+
+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

src/RadialBasisFunctions.jl

@@ -21,12 +21,15 @@ 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/custom.jl")
+export Custom, custom
+
 include("operators/partial.jl")
 export Partial, partial
 

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})

src/basis/polyharmonic_spline.jl

@@ -40,26 +40,26 @@ 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ᵢ)
         return sum(
             (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_INF)
         )
     end
-    return ℒRadialBasisFunction(∇²ℒ)
+    return ∇²ℒ
 end
 
 """
@@ -78,26 +78,26 @@ 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ᵢ)
         return sum(
             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,30 @@ 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
+
+# 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}

src/operators/custom.jl

@@ -0,0 +1,12 @@
+"""
+    Custom <: ScalarValuedOperator
+
+Builds an operator for a first order partial derivative.
+"""
+struct Custom{F<:Function} <: AbstractOperator
+    ℒ::F
+end
+(op::Custom)(basis) = op.ℒ(basis)
+
+# pretty printing
+print_op(op::Custom) = "Custom Operator"