diff --git a/.github/workflows/AirspeedVelocity.yml b/.github/workflows/AirspeedVelocity.yml new file mode 100644 index 00000000..dd5953f6 --- /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.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 \ 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..9cf457e1 --- /dev/null +++ b/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" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 00000000..48d4dc19 --- /dev/null +++ b/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 diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl index 84c2fa03..02316e38 100644 --- a/src/RadialBasisFunctions.jl +++ b/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 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 fa8af001..572f53e6 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,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} diff --git a/src/operators/custom.jl b/src/operators/custom.jl new file mode 100644 index 00000000..b09c2aa4 --- /dev/null +++ b/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" diff --git a/src/operators/directional.jl b/src/operators/directional.jl index b36654e8..d572eec8 100644 --- a/src/operators/directional.jl +++ b/src/operators/directional.jl @@ -1,12 +1,12 @@ """ - Directional <: ScalarValuedOperator + Directional{Dim,T} <: ScalarValuedOperator 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 directional(data, v, basis; k=autoselect_k(data, basis)) @@ -20,8 +20,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 +38,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 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 -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/src/operators/gradient.jl b/src/operators/gradient.jl index 9a0d9a39..7115c406 100644 --- a/src/operators/gradient.jl +++ b/src/operators/gradient.jl @@ -1,10 +1,12 @@ """ - Gradient <: VectorValuedOperator + Gradient{Dim} <: VectorValuedOperator 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 137677d9..110b5967 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)(basis) = ∇²(basis) # 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..c0585732 100644 --- a/src/operators/operator_algebra.jl +++ b/src/operators/operator_algebra.jl @@ -1,38 +1,34 @@ 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) + @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( - 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 + @eval function Base.$op(op1::AbstractOperator, op2::AbstractOperator) + function additive_ℒ(basis) + return additive_ℒrbf(x1, x2) = Base.$op(op1(basis)(x1, x2), op2(basis)(x1, x2)) 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) + 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 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/src/operators/operators.jl b/src/operators/operators.jl index 13cefc3b..2a5c3898 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -1,8 +1,6 @@ abstract type AbstractOperator end abstract type ScalarValuedOperator <: AbstractOperator end -abstract type VectorValuedOperator <: AbstractOperator end - -(op::AbstractOperator)(x) = op.ℒ(x) +abstract type VectorValuedOperator{Dim} <: AbstractOperator end """ struct RadialBasisOperator @@ -56,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 @@ -131,14 +104,15 @@ 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.ℒ.ℒ) - 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 491d1bce..72b3c12d 100644 --- a/src/operators/partial.jl +++ b/src/operators/partial.jl @@ -3,11 +3,11 @@ 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 +(op::Partial)(basis) = ∂(basis, op.order, op.dim) # convienience constructors """ @@ -23,10 +23,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 +41,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/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/basis/monomial.jl b/test/basis/monomial.jl index fab7092f..98c2998d 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(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/directional.jl b/test/operators/directional.jl index 7db218a2..5ca36fc7 100644 --- a/test/operators/directional.jl +++ b/test/operators/directional.jl @@ -46,7 +46,22 @@ 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((1, 2), (3, 4)) + ∇v = Directional{2}(SVector(1.0, 1.0)) @test RadialBasisFunctions.print_op(∇v) == "Directional Derivative (∇f⋅v)" 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/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 diff --git a/test/operators/operator_algebra.jl b/test/operators/operator_algebra.jl index 19227705..321932a7 100644 --- a/test/operators/operator_algebra.jl +++ b/test/operators/operator_algebra.jl @@ -22,3 +22,16 @@ dxdy = dx + dy 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 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 diff --git a/test/runtests.jl b/test/runtests.jl index 6018a719..8b3502e9 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 @@ -59,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