Merge pull request #14 from kylebeggs/data-centers

Data centers
JuliaMeshless · Oct 6, 2023 · 42f33ed · 42f33ed · kylebeggs · Oct 6, 2023
2 parents dfa0983 + 61ffa96
commit 42f33ed
Show file tree

Hide file tree

Showing 17 changed files with 152 additions and 64 deletions.
diff --git a/Project.toml b/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"

diff --git a/docs/make.jl b/docs/make.jl
@@ -20,6 +20,7 @@ makedocs(;
         "Home" => "index.md",
         "Getting Started" => "getting_started.md",
         "Theory" => "theory.md",
+        "API" => "api.md",
     ],
 )
 

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -0,0 +1,4 @@
+```@autodocs
+Modules = [RadialBasisFunctions]
+Order   = [:function, :type]
+```
diff --git a/docs/src/theory.md b/docs/src/theory.md
@@ -96,7 +96,7 @@ In the RBF-FD, a stencil is built to approximate derivatives using the same neig
 
 ```math
 \mathcal{L}u(\mathbf{x}_{c})=\sum_{i=1}^{N}w_{i}u(\mathbf{x}_{i})
-```m
+```
 
 After some algebraic manipulation and substitution using equations, one arrives at a linear system for the weights ``\mathbf{w}`` as
 

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
@@ -24,6 +24,7 @@ include("operators/operator_combinations.jl")
 include("interpolation.jl")
 
 const Δ = ∇² # some people like this notation for the Laplacian
+const DIV0_OFFSET = 1e-8
 
 # operators
 export RadialBasisOperator

diff --git a/src/basis/monomial.jl b/src/basis/monomial.jl
@@ -5,40 +5,48 @@ Multivariate Monomial basis.
 n ∈ N: length of array, i.e., x ∈ Rⁿ
 deg ∈ N: degree
 """
-struct MonomialBasis{T<:Int,B<:Function}
+struct MonomialBasis{T<:Int,B}
     n::T
     deg::T
-    basis::B
+    basis!::B
     function MonomialBasis(n::T, deg::T) where {T<:Int}
-        basis = build_monomial_basis(n, deg)
-        return new{T,typeof(basis)}(n, deg, basis)
+        if deg < 0
+            basis! = nothing
+        else
+            basis! = build_monomial_basis(n, deg)
+        end
+        return new{T,typeof(basis!)}(n, deg, basis!)
     end
 end
 
 function (m::MonomialBasis)(x::AbstractVector{T}) where {T}
     b = ones(T, binomial(m.n + m.deg, m.n))
-    m.basis(b, x)
+    m.basis!(b, x)
     return b
 end
-(m::MonomialBasis)(b, x) = m.basis(b, x)
+(m::MonomialBasis)(b, x) = m.basis!(b, x)
 
 function build_monomial_basis(n::T, deg::T) where {T<:Int}
+    if deg == 0
+        basis! = (b, x) -> b .= one(eltype(x))
+        return basis!
+    end
     e = multiexponents(n + 1, deg)
     e = map(i -> getindex.(e, i), 1:n)
     ids = map(j -> map(i -> findall(x -> x >= i, e[j]), 1:deg), 1:n)
     return build_monomial_basis(ids)
 end
 
 function build_monomial_basis(ids::Vector{Vector{Vector{T}}}) where {T<:Int}
-    function basis(b::AbstractVector{B}, x::AbstractVector) where {B}
+    function basis!(b::AbstractVector{B}, x::AbstractVector) where {B}
         b .= one(B)
         # TODO flatten loop - why does it allocate here
         @views @inbounds for i in eachindex(ids), k in eachindex(ids[i])
             b[ids[i][k]] *= x[i]
         end
         return nothing
     end
-    return basis
+    return basis!
 end
 
 # pretty printing

diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
@@ -39,7 +39,7 @@ end
 
 (phs::PHS1)(x, xᵢ) = euclidean(x, xᵢ)
 function ∂(::PHS1, ::Val{1}, dim::Int)
-    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / euclidean(x, xᵢ)
+    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + DIV0_OFFSET)
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS1)
@@ -48,13 +48,16 @@ function ∇(::PHS1)
 end
 function ∂(::PHS1, ::Val{2}, dim::Int)
     function ∂²ℒ(x, xᵢ)
-        return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + 1e-8)
+        return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) /
+               (euclidean(x, xᵢ)^3 + DIV0_OFFSET)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS1)
     function ∇²ℒ(x, xᵢ)
-        return sum((-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + 1e-8))
+        return sum(
+            (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + DIV0_OFFSET)
+        )
     end
     return ℒRadialBasisFunction(∇²ℒ)
 end
@@ -83,13 +86,16 @@ function ∇(::PHS3)
 end
 function ∂(::PHS3, ::Val{2}, dim::Int)
     function ∂²ℒ(x, xᵢ)
-        return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) / (euclidean(x, xᵢ) + 1e-8)
+        return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) /
+               (euclidean(x, xᵢ) + DIV0_OFFSET)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS3)
     function ∇²ℒ(x, xᵢ)
-        return sum(3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + 1e-8))
+        return sum(
+            3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + DIV0_OFFSET)
+        )
     end
     return ℒRadialBasisFunction(∇²ℒ)
 end

diff --git a/src/interpolation.jl b/src/interpolation.jl
@@ -37,10 +37,12 @@ function (rbfi::RadialBasisInterp)(x::T) where {T}
         rbf += rbfi.rbf_weights[i] * rbfi.rbf_basis(x, rbfi.x[i])
     end
 
-    val_poly = rbfi.monomial_basis(x)
     poly = zero(eltype(T))
-    for i in eachindex(rbfi.monomial_weights)
-        poly += rbfi.monomial_weights[i] * val_poly[i]
+    if !isempty(rbfi.monomial_weights)
+        val_poly = rbfi.monomial_basis(x)
+        for i in eachindex(rbfi.monomial_weights)
+            poly += rbfi.monomial_weights[i] * val_poly[i]
+        end
     end
     return rbf + poly
 end

diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl
@@ -1,5 +1,9 @@
 function _build_weightmx(
-    ℒ, 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
@@ -18,6 +22,7 @@ function _build_weightmx(
     V = zeros(TD, k * length(adjl))
 
     n_threads = Threads.nthreads()
+    n_threads = 1
 
     # create work arrays
     n = sum(sizes)
@@ -27,20 +32,25 @@ function _build_weightmx(
     d = Vector{Vector{eltype(data)}}(undef, n_threads)
 
     # build stencil for each data point and store in global weight matrix
-    Threads.@threads for i in eachindex(adjl)
+    #Threads.@threads for i in eachindex(adjl)
+    for i in eachindex(adjl)
         range[Threads.threadid()] = ((i - 1) * k + 1):(i * k)
         @turbo I[range[Threads.threadid()]] .= i
         d[Threads.threadid()] = data[adjl[i]]
         V[range[Threads.threadid()]] = @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
 
     return sparse(I, J, V, length(adjl), length(data))
 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
@@ -68,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
 
@@ -82,12 +92,13 @@ function _build_stencil!(
     ℒrbf,
     ℒmon,
     data::AbstractVector{D},
+    center,
     basis::B,
     mon::MonomialBasis,
     k::Int,
 ) where {D<:AbstractArray,B<:AbstractRadialBasis}
     _build_collocation_matrix!(A[id], data[id], basis, mon, k)
-    _build_rhs!(b[id], ℒrbf, ℒmon, data[id], basis, k)
+    _build_rhs!(b[id], ℒrbf, ℒmon, data[id], center, basis, k)
     return (A[id] \ b[id])[1:k]
 end
 
@@ -110,18 +121,18 @@ function _build_collocation_matrix!(
 end
 
 function _build_rhs!(
-    b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{D}, basis::B, k::K
+    b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{D}, center, basis::B, k::K
 ) where {D<:AbstractArray,B<:AbstractRadialBasis,K<:Int}
     # radial basis section
     @inbounds for i in eachindex(data)
-        b[i] = ℒrbf(first(data), data[i])
+        b[i] = ℒrbf(center, data[i])
     end
 
     # monomial augmentation
     if basis.poly_deg > -1
         N = length(b)
         bmono = view(b, (k + 1):N)
-        ℒmon(bmono, first(data))
+        ℒmon(bmono, center)
     end
 
     return nothing

diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
@@ -50,7 +50,7 @@ function RadialBasisOperator(
     adjl = find_neighbors(data, k)
     N = length(adjl)
     weights = ntuple(_ -> _allocate_weights(TD, N, N, k; sparse=sparse), length(ℒ.ℒ))
-    return RadialBasisOperator(ℒ, weights, data, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
 end
 
 function RadialBasisOperator(
@@ -66,7 +66,7 @@ function RadialBasisOperator(
     Na = length(adjl)
     Nd = length(data)
     weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k; sparse=sparse), length(ℒ.ℒ))
-    return RadialBasisOperator(ℒ, weights, data, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, centers, adjl, basis)
 end
 
 Base.size(op::RadialBasisOperator{<:Gradient}) = size(first(op.weights))

diff --git a/src/operators/operators.jl b/src/operators/operators.jl
@@ -7,17 +7,18 @@ abstract type VectorValuedOperator <: AbstractOperator end
 
 Operator of data using a radial basis with potential monomial augmentation.
 """
-struct RadialBasisOperator{L,W,D,A,B<:AbstractRadialBasis}
+struct RadialBasisOperator{L,W,D,C,A,B<:AbstractRadialBasis}
     ℒ::L
     weights::W
     data::D
+    centers::C
     adjl::A
     basis::B
     valid_cache::Base.RefValue{Bool}
     function RadialBasisOperator(
-        ℒ::L, weights::W, data::D, adjl::A, basis::B
-    ) where {L,W,D,A,B<:AbstractRadialBasis}
-        return new{L,W,D,A,B}(ℒ, weights, data, adjl, basis, Ref(false))
+        ℒ::L, weights::W, data::D, centers::C, adjl::A, basis::B
+    ) where {L,W,D,C,A,B<:AbstractRadialBasis}
+        return new{L,W,D,C,A,B}(ℒ, weights, data, centers, adjl, basis, Ref(false))
     end
 end
 
@@ -34,7 +35,7 @@ function RadialBasisOperator(
     Na = length(adjl)
     Nd = length(data)
     weights = _allocate_weights(TD, Na, Nd, k; sparse=sparse)
-    return RadialBasisOperator(ℒ, weights, data, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
 end
 
 function RadialBasisOperator(
@@ -50,7 +51,7 @@ function RadialBasisOperator(
     Na = length(adjl)
     Nd = length(data)
     weights = _allocate_weights(TD, Na, Nd, k; sparse=sparse)
-    return RadialBasisOperator(ℒ, weights, data, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, centers, adjl, basis)
 end
 
 # extend Base methods
@@ -151,15 +152,15 @@ function _build_weights(ℒ, op::RadialBasisOperator)
     if op.weights isa Vector{<:Vector}
         return _build_weight_vec(ℒ, op.data, op.adjl, op.basis)
     else
-        return _build_weightmx(ℒ, op.data, op.adjl, op.basis)
+        return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis)
     end
 end
 
 function _build_weights(ℒ, op::RadialBasisOperator{<:VectorValuedOperator})
     if first(op.weights) isa Vector{<:Vector}
         return _build_weight_vec(ℒ, op.data, op.adjl, op.basis)
     else
-        return _build_weightmx(ℒ, op.data, op.adjl, op.basis)
+        return _build_weightmx(ℒ, op.data, op.centers, op.adjl, op.basis)
     end
 end
 

diff --git a/src/utils.jl b/src/utils.jl
@@ -58,6 +58,11 @@ function check_poly_deg(poly_deg)
     return nothing
 end
 
+function scale_cloud(data)
+    furthest_point = maximum(p -> euclidean(first(data), p), data)
+    return data ./ furthest_point
+end
+
 _allocate_weights(m, n, k; sparse=true) = _allocate_weights(Float64, m, n, k; sparse=sparse)
 function _allocate_weights(T, m, n, k; sparse=true)
     return sparse ? spzeros(T, m, n) : [zeros(T, k) for _ in 1:m]

diff --git a/test/linalg/stencil.jl b/test/linalg/stencil.jl
@@ -14,12 +14,6 @@ Lmb = L(mb)
 k = length(x)
 n = k + 3
 
-A = Symmetric(zeros(n, n))
-b = zeros(n)
-
-_build_collocation_matrix!(A, x, rb, mb, k)
-_build_rhs!(b, Lrb, Lmb, x, rb, k)
-
 function unordered_approx(A::AbstractVector, B::AbstractVector)
     length(A) != length(B) && return false
     b = deepcopy(B)
@@ -32,6 +26,8 @@ function unordered_approx(A::AbstractVector, B::AbstractVector)
 end
 
 @testset "Coefficient Matrix" begin
+    A = Symmetric(zeros(n, n))
+    _build_collocation_matrix!(A, x, rb, mb, k)
     @testset "RBFs" begin
         @test A[1, 2] ≈ (sqrt(sum((x[1] .- x[2]) .^ 2)))^3
         @test A[1, 3] ≈ (sqrt(sum((x[1] .- x[3]) .^ 2)))^3
@@ -45,10 +41,13 @@ end
 end
 
 @testset "Right-hand side" begin
+    b = zeros(n)
+    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)
-Original file line number
+Diff line change
@@ Expand Up @@
     ```math
     \mathcal{L}u(\mathbf{x}_{c})=\sum_{i=1}^{N}w_{i}u(\mathbf{x}_{i})
-    ```m
+    ```
     After some algebraic manipulation and substitution using equations, one arrives at a linear system for the weights ``\mathbf{w}`` as
@@ Expand Down @@