correctly handle data centers

Project.toml

@@ -1,7 +1,7 @@
 name = "RadialBasisFunctions"
 uuid = "79ee0514-adf7-4479-8807-6f72ea8967e8"
 authors = ["Kyle Beggs"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/interpolation.jl

@@ -17,15 +17,14 @@ end
 
 Construct a radial basis interpolator.
 """
-function RadialBasisInterp(x, c, y, basis::B=PHS()) where {B<:AbstractRadialBasis}
+function RadialBasisInterp(x, y, basis::B=PHS()) where {B<:AbstractRadialBasis}
  dim = length(first(x))
  k = length(x) # number of data in influence/support domain
  npoly = binomial(dim + basis.poly_deg, basis.poly_deg)
  n = k + npoly
  mon = MonomialBasis(dim, basis.poly_deg)
  A = Symmetric(zeros(eltype(first(x)), n, n))
- _build_collocation_matrix!(A, x, c, basis, mon, k)
- return A
+ _build_collocation_matrix!(A, x, basis, mon, k)
  b = zeros(eltype(first(x)), n)
  b[1:k] .= y
  w = A \ b

src/linalg/stencil.jl

@@ -1,5 +1,9 @@
 function _build_weightmx(
- ℒ, data::AbstractVector{D}, centers, adjl::Vector{Vector{T}}, basis::B
+ ℒ,
+ data::AbstractVector{D},
+ centers::AbstractVector{D},
+ adjl::Vector{Vector{T}},
+ basis::B,
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
  TD = eltype(first(data))
  dim = length(first(data)) # dimension of data
@@ -42,7 +46,11 @@ function _build_weightmx(
 end
 
 function _build_weight_vec(
- ℒ, data::AbstractVector{D}, adjl::Vector{Vector{T}}, basis::B
+ ℒ,
+ data::AbstractVector{D},
+ centers::AbstractVector{D},
+ adjl::Vector{Vector{T}},
+ basis::B,
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
  TD = eltype(first(data))
  dim = length(first(data)) # dimension of data
@@ -70,7 +78,7 @@ function _build_weight_vec(
  Threads.@threads for i in eachindex(adjl)
  d[Threads.threadid()] = data[adjl[i]]
  V[i] = @views _build_stencil!(
- A, b, Threads.threadid(), ℒrbf, ℒmon, d, basis, mon, k
+ A, b, Threads.threadid(), ℒrbf, ℒmon, d, centers[i], basis, mon, k
  )
  end
 

test/linalg/stencil.jl

@@ -42,11 +42,12 @@ end
 
 @testset "Right-hand side" begin
  b = zeros(n)
- _build_rhs!(b, Lrb, Lmb, x, rb, k)
+ center = SVector(0.0, 0.0)
+ _build_rhs!(b, Lrb, Lmb, x, center, rb, k)
  @testset "RBFs" begin
- @test b[1] ≈ Lrb(x[1], x[1])
- @test b[2] ≈ Lrb(x[1], x[2])
- @test b[3] ≈ Lrb(x[1], x[3])
+ @test b[1] ≈ Lrb(center, x[1])
+ @test b[2] ≈ Lrb(center, x[2])
+ @test b[3] ≈ Lrb(center, x[3])
  end
  @testset "Monomials" begin
  bb = zeros(3)

test/operators/laplacian.jl

@@ -26,8 +26,8 @@ y = f.(x)
  @test mean_percent_error(∇²(y), ∇²f.(x)) < 2
 end
 
-#@testset "Different Data Centers" begin
-# x2 = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), (N + 1):(N + 1000))
-# ∇² = laplacian(x, x2, PHS(3; poly_deg=4))
-# @test mean_percent_error(∇²(y), ∇²f.(x2)) < 2
-#end
+@testset "Different Data Centers" begin
+ x2 = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), (N + 1):(N + 1000))
+ ∇² = laplacian(x, x2, PHS(3; poly_deg=4))
+ @test mean_percent_error(∇²(y), ∇²f.(x2)) < 2
+end