Merge pull request #26 from kylebeggs/feature/virtual

Feature/virtual
kylebeggs · Jul 17, 2024 · 5ca0c01 · 5ca0c01 · kylebeggs · Jul 17, 2024
2 parents 56a7811 + 4a96b39
commit 5ca0c01
Show file tree

Hide file tree

Showing 20 changed files with 280 additions and 104 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.2.1"
+version = "0.2.2"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -1,12 +1,12 @@
-# Exported Functions
+## Exported Functions
 
 ```@autodocs
 Modules = [RadialBasisFunctions]
 Private = false
 Order   = [:function, :type]
 ```
 
-# Private
+## Private
 
 ```@autodocs
 Modules = [RadialBasisFunctions]

diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
@@ -50,8 +50,8 @@ abs.(y_true .- y_new)
 This package also provides an API for operators. There is support for several built-in operators along with support for user-defined operators. Currently, we have implementations for
 
 - partial derivative (1st and 2nd order)
-- gradient
 - laplacian
+- gradient
 
 but we plan to add more in the future. Please make and issue or pull request for additional operators.
 

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
@@ -23,7 +23,7 @@ export find_neighbors, reorder_points!
 include("linalg/stencil.jl")
 
 include("operators/operators.jl")
-export RadialBasisOperator
+export RadialBasisOperator, update_weights!
 
 include("operators/partial.jl")
 export Partial, partial
@@ -37,6 +37,9 @@ export Gradient, gradient
 include("operators/directional.jl")
 export Directional, directional
 
+include("operators/virtual.jl")
+export ∂virtual
+
 include("operators/monomial.jl")
 
 include("operators/operator_combinations.jl")

diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl
@@ -1,19 +1,28 @@
 function _build_weights(ℒ, op; nchunks=Threads.nthreads())
+    data = op.data
+    basis = op.basis
+    dim = length(first(data)) # dimension of data
+
+    # build monomial basis and operator
+    mon = MonomialBasis(dim, basis.poly_deg)
+    ℒmon = ℒ(mon)
+    ℒrbf = ℒ(basis)
+
+    return _build_weights(op, ℒrbf, ℒmon, mon; nchunks=nchunks)
+end
+
+function _build_weights(op, ℒrbf, ℒmon, mon; nchunks=Threads.nthreads())
     data = op.data
     eval_points = op.eval_points
     adjl = op.adjl
     basis = op.basis
+
     TD = eltype(first(data))
     dim = length(first(data)) # dimension of data
     nmon = binomial(dim + basis.poly_deg, basis.poly_deg)
     k = length(first(adjl))  # number of data in influence/support domain
     sizes = (k, nmon)
 
-    # build monomial basis and operator
-    mon = MonomialBasis(dim, basis.poly_deg)
-    ℒmon = ℒ(mon)
-    ℒrbf = ℒ(basis)
-
     # allocate arrays to build sparse matrix
     Na = length(adjl)
     I = zeros(Int, k * Na)
@@ -45,12 +54,12 @@ function _build_stencil!(
     b::Vector,
     ℒrbf,
     ℒmon,
-    data::AbstractVector{D},
-    eval_point::D,
+    data::AbstractVector{TD},
+    eval_point::TE,
     basis::B,
     mon::MonomialBasis,
     k::Int,
-) where {D<:AbstractArray,B<:AbstractRadialBasis}
+) 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]
@@ -75,8 +84,8 @@ function _build_collocation_matrix!(
 end
 
 function _build_rhs!(
-    b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{D}, eval_point::D, basis::B, k::K
-) where {D<:AbstractArray,B<:AbstractRadialBasis,K<:Int}
+    b::AbstractVector, ℒ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)
         b[i] = ℒrbf(eval_point, data[i])

diff --git a/src/operators/directional.jl b/src/operators/directional.jl
@@ -1,9 +1,9 @@
 """
-    Directional <: VectorValuedOperator
+    Directional <: ScalarValuedOperator
 
 Operator for the directional derivative, or the inner product of the gradient and a direction vector.
 """
-struct Directional{L<:NTuple,T} <: VectorValuedOperator
+struct Directional{L<:NTuple,T} <: ScalarValuedOperator
     ℒ::L
     v::T
 end
@@ -18,14 +18,15 @@ function directional(
     v::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
         end
     end
     ℒ = Directional(f, v)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 """
@@ -34,41 +35,68 @@ end
 Builds a `RadialBasisOperator` where the operator is the directional derivative, `Directional`.
 """
 function directional(
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector,
+    eval_points::AbstractVector,
     v::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+    adjl=find_neighbors(data, eval_points, k),
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
         end
     end
     ℒ = Directional(f, v)
-    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
+    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
+end
+
+function RadialBasisOperator(
+    ℒ::Directional,
+    data::AbstractVector{D},
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
+) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+    Na = length(adjl)
+    Nd = length(data)
+    weights = spzeros(eltype(D), Na, Nd)
+    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}
+    Na = length(adjl)
+    Nd = length(data)
+    weights = spzeros(eltype(TD), Na, Nd)
+    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end
 
-function _update_weights!(
-    op::RadialBasisOperator{<:Directional}, weights::NTuple{N,AbstractMatrix}
-) where {N}
+function update_weights!(op::RadialBasisOperator{<:Directional})
     v = op.ℒ.v
+    N = length(first(op.data))
     @assert length(v) == N || length(v) == size(op)[1] "wrong size for v"
     if length(v) == N
-        for (i, ℒ) in enumerate(op.ℒ.ℒ)
-            weights[i] .= _build_weights(ℒ, op) * v[i]
+        op.weights .= mapreduce(+, enumerate(op.ℒ.ℒ)) do (i, ℒ)
+            _build_weights(ℒ, op) * v[i]
         end
     else
         vv = ntuple(i -> getindex.(v, i), N)
-        for (i, ℒ) in enumerate(op.ℒ.ℒ)
-            weights[i] .= Diagonal(vv[i]) * _build_weights(ℒ, op)
+        op.weights .= mapreduce(+, enumerate(op.ℒ.ℒ)) do (i, ℒ)
+            Diagonal(vv[i]) * _build_weights(ℒ, op)
         end
     end
     validate_cache(op)
     return nothing
 end
 
-Base.size(op::RadialBasisOperator{<:Directional}) = size(first(op.weights))
+Base.size(op::RadialBasisOperator{<:Directional}) = size(op.weights)
 
 # pretty printing
 print_op(op::Directional) = "Directional Gradient (∇f⋅v)"
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
@@ -14,15 +14,18 @@ end
 Builds a `RadialBasisOperator` where the operator is the gradient, `Gradient`.
 """
 function gradient(
-    data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis)
+    data::AbstractVector{D},
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
         end
     end
     ℒ = Gradient(f)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 """
@@ -31,18 +34,19 @@ end
 Builds a `RadialBasisOperator` where the operator is the gradient, `Gradient`. The resulting operator will only evaluate at `eval_points`.
 """
 function gradient(
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector,
+    eval_points::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+    adjl=find_neighbors(data, eval_points, k),
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
         end
     end
     ℒ = Gradient(f)
-    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
+    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end
 
 Base.size(op::RadialBasisOperator{<:Gradient}) = size(first(op.weights))

diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
@@ -9,20 +9,24 @@ end
 
 # convienience constructors
 function laplacian(
-    data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis)
+    data::AbstractVector{D},
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
     ℒ = Laplacian(∇²)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 function laplacian(
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector,
+    eval_points::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+    adjl=find_neighbors(data, eval_points, k),
+) where {T<:Int,B<:AbstractRadialBasis}
     ℒ = Laplacian(∇²)
-    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
+    return RadialBasisOperator(ℒ, data, eval_points, basis; k=k, adjl=adjl)
 end
 
 # pretty printing

diff --git a/src/operators/operator_combinations.jl b/src/operators/operator_combinations.jl
@@ -36,7 +36,7 @@ for op in (:+, :-, :*, :/)
         k2 = length(first((op2.adjl)))
         k = k1 > k2 ? k1 : k2
         ℒ = Base.$op(op1.ℒ, op2.ℒ)
-        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k)
+        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
     end
 end
 

diff --git a/src/operators/operators.jl b/src/operators/operators.jl
@@ -24,9 +24,12 @@ end
 
 # convienience constructors
 function RadialBasisOperator(
-    ℒ, data::AbstractVector{D}, basis::B=PHS(3; poly_deg=2); k::T=autoselect_k(data, basis)
+    ℒ,
+    data::AbstractVector{D},
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    adjl = find_neighbors(data, k)
     Na = length(adjl)
     Nd = length(data)
     weights = spzeros(eltype(D), Na, Nd)
@@ -35,15 +38,15 @@ end
 
 function RadialBasisOperator(
     ℒ,
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector{TD},
+    eval_points::AbstractVector{TE},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    adjl = find_neighbors(data, eval_points, k)
+    adjl=find_neighbors(data, eval_points, k),
+) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
     Na = length(adjl)
     Nd = length(data)
-    weights = spzeros(eltype(D), Na, Nd)
+    weights = spzeros(eltype(TD), Na, Nd)
     return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end
 
@@ -52,26 +55,25 @@ function RadialBasisOperator(
     data::AbstractVector{D},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, k),
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
     TD = eltype(D)
-    adjl = find_neighbors(data, k)
     N = length(adjl)
     weights = ntuple(_ -> _allocate_weights(TD, N, N, k), length(ℒ.ℒ))
     return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
 end
 
 function RadialBasisOperator(
     ℒ::VectorValuedOperator,
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector{TD},
+    eval_points::AbstractVector{TE},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    TD = eltype(D)
-    adjl = find_neighbors(data, eval_points, k)
+    adjl=find_neighbors(data, eval_points, k),
+) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
     Na = length(adjl)
     Nd = length(data)
-    weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k), length(ℒ.ℒ))
+    weights = ntuple(_ -> _allocate_weights(eltype(TD), Na, Nd, k), length(ℒ.ℒ))
     return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end