From 6d3d910fdb233033817b881ed5fdcd8c3a3b8072 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sun, 21 Apr 2024 13:49:54 -0400
Subject: [PATCH 1/9] add upwinding op

---
 src/RadialBasisFunctions.jl |  5 ++-
 src/operators/upwind.jl     | 69 +++++++++++++++++++++++++++++++++++++
 2 files changed, 73 insertions(+), 1 deletion(-)
 create mode 100644 src/operators/upwind.jl

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index 946e597..3de747e 100644
--- 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/upwind.jl")
+export Upwind, upwind
+
 include("operators/monomial.jl")
 
 include("operators/operator_combinations.jl")
diff --git a/src/operators/upwind.jl b/src/operators/upwind.jl
new file mode 100644
index 0000000..32bef0b
--- /dev/null
+++ b/src/operators/upwind.jl
@@ -0,0 +1,69 @@
+struct Upwind{L<:Function,T<:Int} <: AbstractOperator
+    ℒ::L
+    dim::T
+end
+
+# convienience constructors
+"""
+    function upwind(data, dim, basis; k=autoselect_k(data, basis))
+
+Builds a `RadialBasisOperator` where the operator is the partial derivative, `Partial`, of `order` with respect to `dim`.
+"""
+function upwind(
+    data::AbstractVector{D},
+    dim,
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+    tree = KDTree(data)
+    _, dists = knn(tree, data, k, true)
+    Δ = minimum(dists) do d
+        z = minimum(@view(d[2:end])) do i
+            abs(i - first(d))
+        end
+        return z
+    end
+
+    N = length(first(data))
+    dx = zeros(N)
+    dx[dim] = Δ
+
+    li = regrid(data, data .- Ref(dx))
+    ri = regrid(data, data .+ Ref(dx))
+    update_weights!(li)
+    update_weights!(ri)
+    one_typed = one(eltype(first(data)))
+    l = columnwise_div((sparse(one_typed * I, size(li.weights)...) .- li.weights), Δ)
+    r = columnwise_div((sparse(one_typed * I, size(ri.weights)...) .- ri.weights), Δ)
+    c = partial(data, 1, dim)
+
+    du = let l = l, r = r, c = c
+        (u, v, θ) -> begin
+            wl = max.(v, Ref(0)) .* (θ * (l * u) .+ (1 - θ) * c(u))
+            wr = min.(v, Ref(0)) .* (θ * (r * u) .+ (1 - θ) * c(u))
+            wl .+ wr
+        end
+    end
+
+    return du, l, r
+end
+
+function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)
+    I, J, V = findnz(A)
+    for idx in eachindex(V)
+        V[idx] /= B[I[idx]]
+    end
+    return sparse(I, J, V)
+end
+columnwise_div(A::SparseMatrixCSC, B::Number) = A ./ B
+
+separate(f, x) = (t = findall(f, x); (@view(x[t]), @view(x[setdiff(eachindex(x), t)])))
+
+function split_stencils(x, adjl, dim)
+    split = map(adjl) do a
+        center = first(a)
+        l, r = separate(i -> x[i][dim] < x[center][dim], @view(a[2:end]))
+        return vcat(center, l), vcat(center, r)
+    end
+    return getindex.(split, 1), getindex.(split, 2)
+end

From a067a41b52fc5d49f82bf8d88f4a98777b89dbfa Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Sun, 21 Apr 2024 14:38:28 -0400
Subject: [PATCH 2/9] fix directional deivative to be scalar output

---
 src/operators/directional.jl  | 46 +++++++++++++++++++++++++++--------
 test/operators/directional.jl | 19 +++++----------
 2 files changed, 42 insertions(+), 23 deletions(-)

diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index 837abac..4b1bc62 100644
--- 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
@@ -49,26 +49,52 @@ function directional(
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
 end
 
-function _update_weights!(
-    op::RadialBasisOperator{<:Directional}, weights::NTuple{N,AbstractMatrix}
-) where {N}
+function RadialBasisOperator(
+    ℒ::Directional,
+    data::AbstractVector{D},
+    basis::B=PHS(3; poly_deg=2);
+    k::T=autoselect_k(data, basis),
+) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+    adjl = find_neighbors(data, k)
+    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{D},
+    eval_points::AbstractVector{D},
+    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)
+    Na = length(adjl)
+    Nd = length(data)
+    weights = spzeros(eltype(D), Na, Nd)
+    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
+end
+
+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/test/operators/directional.jl b/test/operators/directional.jl
index bcae41d..69a425a 100644
--- a/test/operators/directional.jl
+++ b/test/operators/directional.jl
@@ -21,9 +21,8 @@ y = f.(x)
     v = SVector(2.0, 1.0)
     v /= norm(v)
     ∇v = directional(x, v, PHS3(2))
-    ∇vy = ∇v(y)
-    @test mean_percent_error(∇vy[1], df_dx.(x) .* v[1]) < 2
-    @test mean_percent_error(∇vy[2], df_dy.(x) .* v[2]) < 2
+    exact = map(x -> SVector(df_dx(x), df_dy(x)) ⋅ v, x)
+    @test mean_percent_error(∇v(y), exact) < 2
 end
 
 @testset "Direction Vector for Each Data Center" begin
@@ -32,11 +31,8 @@ end
         return v /= norm(v)
     end
     ∇v = directional(x, v, PHS3(2))
-    ∇vy = ∇v(y)
-    df_dx_v = map((df, v) -> df * v[1], df_dx.(x), v)
-    df_dy_v = map((df, v) -> df * v[2], df_dy.(x), v)
-    @test mean_percent_error(∇vy[1], df_dx_v) < 2
-    @test mean_percent_error(∇vy[2], df_dy_v) < 2
+    exact = map((x, vv) -> SVector(df_dx(x), df_dy(x)) ⋅ vv, x, v)
+    @test mean_percent_error(∇v(y), exact) < 2
 end
 
 @testset "Different Evaluation Points" begin
@@ -46,9 +42,6 @@ end
         return v /= norm(v)
     end
     ∇v = directional(x, x2, v, PHS3(2))
-    ∇vy = ∇v(y)
-    df_dx_v = map((df, v) -> df * v[1], df_dx.(x), v)
-    df_dy_v = map((df, v) -> df * v[2], df_dy.(x), v)
-    @test mean_percent_error(∇vy[1], df_dx_v) < 2
-    @test mean_percent_error(∇vy[2], df_dy_v) < 2
+    exact = map((x, vv) -> SVector(df_dx(x), df_dy(x)) ⋅ vv, x2, v)
+    @test mean_percent_error(∇v(y), exact) < 2
 end

From 24650deb6218b0889db5ca169a972337039c7e0f Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Mon, 22 Apr 2024 09:08:02 -0400
Subject: [PATCH 3/9] remove some type restrictions

---
 src/operators/directional.jl | 14 +++++++-------
 src/operators/gradient.jl    |  6 +++---
 src/operators/laplacian.jl   |  6 +++---
 src/operators/operators.jl   | 17 ++++++++---------
 src/operators/partial.jl     |  6 +++---
 src/operators/regridding.jl  |  6 +++---
 6 files changed, 27 insertions(+), 28 deletions(-)

diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index 4b1bc62..323f7f6 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -34,12 +34,12 @@ 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}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
@@ -64,15 +64,15 @@ end
 
 function RadialBasisOperator(
     ℒ::Directional,
-    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}
+) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     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
 
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index 692c15e..759784c 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -31,11 +31,11 @@ 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}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
             x -> ∂(x, 1, dim)
diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
index c493346..1ccc62e 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -16,11 +16,11 @@ function laplacian(
 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}
+) where {T<:Int,B<:AbstractRadialBasis}
     ℒ = Laplacian(∇²)
     return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
 end
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
index cc1e8a3..3af23fe 100644
--- a/src/operators/operators.jl
+++ b/src/operators/operators.jl
@@ -35,15 +35,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}
+) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     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
 
@@ -62,16 +62,15 @@ 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)
+) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     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
 
diff --git a/src/operators/partial.jl b/src/operators/partial.jl
index 5c26ead..df83402 100644
--- a/src/operators/partial.jl
+++ b/src/operators/partial.jl
@@ -35,13 +35,13 @@ end
 Builds a `RadialBasisOperator` where the operator is the partial derivative, `Partial`. The resulting operator will only evaluate at `eval_points`.
 """
 function partial(
-    data::AbstractVector{D},
-    eval_points::AbstractVector{D},
+    data::AbstractVector,
+    eval_points::AbstractVector,
     order::T,
     dim::T,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {T<:Int,B<:AbstractRadialBasis}
     f = let o = order, dim = dim
         x -> ∂(x, o, dim)
     end
diff --git a/src/operators/regridding.jl b/src/operators/regridding.jl
index 32582b3..ee9617d 100644
--- a/src/operators/regridding.jl
+++ b/src/operators/regridding.jl
@@ -15,11 +15,11 @@ end
 Builds a `RadialBasisOperator` where the operator is an regrid from one set of points to another, `data` -> `eval_points`.
 """
 function regrid(
-    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}
+) where {T<:Int,B<:AbstractRadialBasis}
     return RadialBasisOperator(Regrid(), data, eval_points, basis; k=k)
 end
 

From 05b75a1400add57aa14658d16d04594673fad4a2 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Mon, 22 Apr 2024 09:49:48 -0400
Subject: [PATCH 4/9] update stencil.jl with new looser type restrictions

---
 src/linalg/stencil.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl
index 699fbdf..eb17f9d 100644
--- a/src/linalg/stencil.jl
+++ b/src/linalg/stencil.jl
@@ -45,12 +45,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<:AbstractArray,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 +75,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<:AbstractArray,TE,B<:AbstractRadialBasis,K<:Int}
     # radial basis section
     @inbounds for i in eachindex(data)
         b[i] = ℒrbf(eval_point, data[i])

From e1c2f6b06f08fbd0afb088611a43f2e4b49cace3 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Thu, 25 Apr 2024 17:58:07 -0400
Subject: [PATCH 5/9] fix return on upwind()

---
 src/operators/upwind.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/operators/upwind.jl b/src/operators/upwind.jl
index 32bef0b..b95f884 100644
--- a/src/operators/upwind.jl
+++ b/src/operators/upwind.jl
@@ -28,24 +28,24 @@ function upwind(
     dx = zeros(N)
     dx[dim] = Δ
 
-    li = regrid(data, data .- Ref(dx))
-    ri = regrid(data, data .+ Ref(dx))
+    li = regrid(data, data .- Ref(dx); k=k)
+    ri = regrid(data, data .+ Ref(dx); k=k)
     update_weights!(li)
     update_weights!(ri)
     one_typed = one(eltype(first(data)))
     l = columnwise_div((sparse(one_typed * I, size(li.weights)...) .- li.weights), Δ)
     r = columnwise_div((sparse(one_typed * I, size(ri.weights)...) .- ri.weights), Δ)
-    c = partial(data, 1, dim)
+    c = partial(data, 1, dim; k=k)
 
     du = let l = l, r = r, c = c
-        (u, v, θ) -> begin
-            wl = max.(v, Ref(0)) .* (θ * (l * u) .+ (1 - θ) * c(u))
-            wr = min.(v, Ref(0)) .* (θ * (r * u) .+ (1 - θ) * c(u))
+        (ϕ, v, θ) -> begin
+            wl = max.(v, Ref(0)) .* (θ * (l * ϕ) .+ (1 - θ) * c(ϕ))
+            wr = min.(v, Ref(0)) .* (θ * (r * ϕ) .+ (1 - θ) * c(ϕ))
             wl .+ wr
         end
     end
 
-    return du, l, r
+    return du
 end
 
 function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)

From 7849e2fb52d8a06e2b9b05bd8b2bc6081dd21f11 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 30 Apr 2024 10:47:20 -0400
Subject: [PATCH 6/9] fix upwinding, turning wave now works

---
 src/linalg/stencil.jl        |  4 ++--
 src/operators/directional.jl |  2 +-
 src/operators/operators.jl   |  4 ++--
 src/operators/upwind.jl      | 34 ++++++++++++++++++++--------------
 4 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl
index eb17f9d..5ddf7b6 100644
--- a/src/linalg/stencil.jl
+++ b/src/linalg/stencil.jl
@@ -50,7 +50,7 @@ function _build_stencil!(
     basis::B,
     mon::MonomialBasis,
     k::Int,
-) where {TD<:AbstractArray,TE,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]
@@ -76,7 +76,7 @@ end
 
 function _build_rhs!(
     b::AbstractVector, ℒrbf, ℒmon, data::AbstractVector{TD}, eval_point::TE, basis::B, k::K
-) where {TD<:AbstractArray,TE,B<:AbstractRadialBasis,K<:Int}
+) 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
index 323f7f6..296616c 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -68,7 +68,7 @@ function RadialBasisOperator(
     eval_points::AbstractVector{TE},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
index 3af23fe..9d1b3e6 100644
--- a/src/operators/operators.jl
+++ b/src/operators/operators.jl
@@ -39,7 +39,7 @@ function RadialBasisOperator(
     eval_points::AbstractVector{TE},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
@@ -66,7 +66,7 @@ function RadialBasisOperator(
     eval_points::AbstractVector{TE},
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
-) where {TD<:AbstractArray,TE<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
     adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
diff --git a/src/operators/upwind.jl b/src/operators/upwind.jl
index b95f884..2286650 100644
--- a/src/operators/upwind.jl
+++ b/src/operators/upwind.jl
@@ -12,17 +12,11 @@ Builds a `RadialBasisOperator` where the operator is the partial derivative, `Pa
 function upwind(
     data::AbstractVector{D},
     dim,
+    Δ=nothing,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
 ) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    tree = KDTree(data)
-    _, dists = knn(tree, data, k, true)
-    Δ = minimum(dists) do d
-        z = minimum(@view(d[2:end])) do i
-            abs(i - first(d))
-        end
-        return z
-    end
+    Δ === nothing && (Δ = _find_smallest_dist(data, k))
 
     N = length(first(data))
     dx = zeros(N)
@@ -33,14 +27,14 @@ function upwind(
     update_weights!(li)
     update_weights!(ri)
     one_typed = one(eltype(first(data)))
-    l = columnwise_div((sparse(one_typed * I, size(li.weights)...) .- li.weights), Δ)
-    r = columnwise_div((sparse(one_typed * I, size(ri.weights)...) .- ri.weights), Δ)
-    c = partial(data, 1, dim; k=k)
+    backward = columnwise_div((sparse(one_typed * I, size(li.weights)...) .- li.weights), Δ)
+    forward = columnwise_div((-sparse(one_typed * I, size(ri.weights)...) .+ ri.weights), Δ)
+    center = partial(data, 1, dim; k=k)
 
-    du = let l = l, r = r, c = c
+    du = let backward = backward, forward = forward, center = center
         (ϕ, v, θ) -> begin
-            wl = max.(v, Ref(0)) .* (θ * (l * ϕ) .+ (1 - θ) * c(ϕ))
-            wr = min.(v, Ref(0)) .* (θ * (r * ϕ) .+ (1 - θ) * c(ϕ))
+            wl = max.(v, Ref(0)) .* (θ * (backward * ϕ) .+ (1 - θ) * center(ϕ))
+            wr = min.(v, Ref(0)) .* (θ * (forward * ϕ) .+ (1 - θ) * center(ϕ))
             wl .+ wr
         end
     end
@@ -57,6 +51,18 @@ function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)
 end
 columnwise_div(A::SparseMatrixCSC, B::Number) = A ./ B
 
+function _find_smallest_dist(data, k)
+    tree = KDTree(data)
+    _, dists = knn(tree, data, k, true)
+    Δ = minimum(dists) do d
+        z = minimum(@view(d[2:end])) do i
+            abs(i - first(d))
+        end
+        return z
+    end
+    return Δ
+end
+
 separate(f, x) = (t = findall(f, x); (@view(x[t]), @view(x[setdiff(eachindex(x), t)])))
 
 function split_stencils(x, adjl, dim)

From 87072e3f6903e0e1be0f94aec0233d8e325f5c61 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 16 Jul 2024 21:35:40 -0400
Subject: [PATCH 7/9] add tests for virtual

---
 src/RadialBasisFunctions.jl            |  4 +-
 src/linalg/stencil.jl                  | 19 +++++--
 src/operators/directional.jl           | 10 ++--
 src/operators/gradient.jl              | 10 ++--
 src/operators/laplacian.jl             | 10 ++--
 src/operators/operator_combinations.jl |  2 +-
 src/operators/operators.jl             | 13 +++--
 src/operators/partial.jl               |  6 ++-
 src/operators/regridding.jl            |  3 +-
 src/operators/upwind.jl                | 75 --------------------------
 src/operators/virtual.jl               | 56 +++++++++++++++++++
 src/utils.jl                           | 27 ++++++++--
 test/operators/virtual.jl              | 49 +++++++++++++++++
 test/runtests.jl                       |  4 ++
 14 files changed, 184 insertions(+), 104 deletions(-)
 delete mode 100644 src/operators/upwind.jl
 create mode 100644 src/operators/virtual.jl
 create mode 100644 test/operators/virtual.jl

diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index 3de747e..5ee2c88 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -37,8 +37,8 @@ export Gradient, gradient
 include("operators/directional.jl")
 export Directional, directional
 
-include("operators/upwind.jl")
-export Upwind, upwind
+include("operators/virtual.jl")
+export ∂virtual
 
 include("operators/monomial.jl")
 
diff --git a/src/linalg/stencil.jl b/src/linalg/stencil.jl
index 5ddf7b6..4f4fd42 100644
--- 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)
diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index 296616c..0b71639 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -18,6 +18,7 @@ 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
@@ -25,7 +26,7 @@ function directional(
         end
     end
     ℒ = Directional(f, v)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 """
@@ -39,6 +40,7 @@ function directional(
     v::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, eval_points, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
@@ -46,7 +48,7 @@ function directional(
         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(
@@ -54,8 +56,8 @@ 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}
-    adjl = find_neighbors(data, k)
     Na = length(adjl)
     Nd = length(data)
     weights = spzeros(eltype(D), Na, Nd)
@@ -68,8 +70,8 @@ function RadialBasisOperator(
     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}
-    adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
     weights = spzeros(eltype(TD), Na, Nd)
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index 759784c..0c01746 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -14,7 +14,10 @@ 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
@@ -22,7 +25,7 @@ function gradient(
         end
     end
     ℒ = Gradient(f)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 """
@@ -35,6 +38,7 @@ function gradient(
     eval_points::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, eval_points, k),
 ) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(length(first(data))) do dim
         return let dim = dim
@@ -42,7 +46,7 @@ function gradient(
         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
index 1ccc62e..8fedacf 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -9,10 +9,13 @@ 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(
@@ -20,9 +23,10 @@ function laplacian(
     eval_points::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     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)
+    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
index 2b7cbf2..6f03f86 100644
--- 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
index 9d1b3e6..031100b 100644
--- 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)
@@ -39,8 +42,8 @@ function RadialBasisOperator(
     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}
-    adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
     weights = spzeros(eltype(TD), Na, Nd)
@@ -52,9 +55,9 @@ 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)
@@ -66,8 +69,8 @@ function RadialBasisOperator(
     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}
-    adjl = find_neighbors(data, eval_points, k)
     Na = length(adjl)
     Nd = length(data)
     weights = ntuple(_ -> _allocate_weights(eltype(TD), Na, Nd, k), length(ℒ.ℒ))
diff --git a/src/operators/partial.jl b/src/operators/partial.jl
index df83402..5dac278 100644
--- a/src/operators/partial.jl
+++ b/src/operators/partial.jl
@@ -21,12 +21,13 @@ function partial(
     dim::T,
     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}
     f = let o = order, dim = dim
         x -> ∂(x, o, dim)
     end
     ℒ = Partial(f, order, dim)
-    return RadialBasisOperator(ℒ, data, basis; k=k)
+    return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
 
 """
@@ -41,12 +42,13 @@ function partial(
     dim::T,
     basis::B=PHS(3; poly_deg=2);
     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)
-    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/regridding.jl b/src/operators/regridding.jl
index ee9617d..7d138fc 100644
--- a/src/operators/regridding.jl
+++ b/src/operators/regridding.jl
@@ -19,8 +19,9 @@ function regrid(
     eval_points::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
+    adjl=find_neighbors(data, eval_points, k),
 ) where {T<:Int,B<:AbstractRadialBasis}
-    return RadialBasisOperator(Regrid(), data, eval_points, basis; k=k)
+    return RadialBasisOperator(Regrid(), data, eval_points, basis; k=k, adjl=adjl)
 end
 
 # pretty printing
diff --git a/src/operators/upwind.jl b/src/operators/upwind.jl
deleted file mode 100644
index 2286650..0000000
--- a/src/operators/upwind.jl
+++ /dev/null
@@ -1,75 +0,0 @@
-struct Upwind{L<:Function,T<:Int} <: AbstractOperator
-    ℒ::L
-    dim::T
-end
-
-# convienience constructors
-"""
-    function upwind(data, dim, basis; k=autoselect_k(data, basis))
-
-Builds a `RadialBasisOperator` where the operator is the partial derivative, `Partial`, of `order` with respect to `dim`.
-"""
-function upwind(
-    data::AbstractVector{D},
-    dim,
-    Δ=nothing,
-    basis::B=PHS(3; poly_deg=2);
-    k::T=autoselect_k(data, basis),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    Δ === nothing && (Δ = _find_smallest_dist(data, k))
-
-    N = length(first(data))
-    dx = zeros(N)
-    dx[dim] = Δ
-
-    li = regrid(data, data .- Ref(dx); k=k)
-    ri = regrid(data, data .+ Ref(dx); k=k)
-    update_weights!(li)
-    update_weights!(ri)
-    one_typed = one(eltype(first(data)))
-    backward = columnwise_div((sparse(one_typed * I, size(li.weights)...) .- li.weights), Δ)
-    forward = columnwise_div((-sparse(one_typed * I, size(ri.weights)...) .+ ri.weights), Δ)
-    center = partial(data, 1, dim; k=k)
-
-    du = let backward = backward, forward = forward, center = center
-        (ϕ, v, θ) -> begin
-            wl = max.(v, Ref(0)) .* (θ * (backward * ϕ) .+ (1 - θ) * center(ϕ))
-            wr = min.(v, Ref(0)) .* (θ * (forward * ϕ) .+ (1 - θ) * center(ϕ))
-            wl .+ wr
-        end
-    end
-
-    return du
-end
-
-function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)
-    I, J, V = findnz(A)
-    for idx in eachindex(V)
-        V[idx] /= B[I[idx]]
-    end
-    return sparse(I, J, V)
-end
-columnwise_div(A::SparseMatrixCSC, B::Number) = A ./ B
-
-function _find_smallest_dist(data, k)
-    tree = KDTree(data)
-    _, dists = knn(tree, data, k, true)
-    Δ = minimum(dists) do d
-        z = minimum(@view(d[2:end])) do i
-            abs(i - first(d))
-        end
-        return z
-    end
-    return Δ
-end
-
-separate(f, x) = (t = findall(f, x); (@view(x[t]), @view(x[setdiff(eachindex(x), t)])))
-
-function split_stencils(x, adjl, dim)
-    split = map(adjl) do a
-        center = first(a)
-        l, r = separate(i -> x[i][dim] < x[center][dim], @view(a[2:end]))
-        return vcat(center, l), vcat(center, r)
-    end
-    return getindex.(split, 1), getindex.(split, 2)
-end
diff --git a/src/operators/virtual.jl b/src/operators/virtual.jl
new file mode 100644
index 0000000..16c8afd
--- /dev/null
+++ b/src/operators/virtual.jl
@@ -0,0 +1,56 @@
+# convienience constructors
+"""
+    function ∂virtual(data, eval_points, dim, Δ, basis; k=autoselect_k(data, basis))
+
+Builds a virtual `RadialBasisOperator` whichi will be evaluated at `eval_points` where the
+operator is the partial derivative with respect to `dim`. Virtual operators interpolate the
+data to structured points at a distance `Δ` for which standard finite difference formulas
+can be applied.
+"""
+function ∂virtual(
+    data::AbstractVector,
+    eval_points::AbstractVector,
+    dim,
+    Δ,
+    basis::B=PHS(3; poly_deg=2);
+    backward=false,
+    k::T=autoselect_k(data, basis),
+) where {T<:Int,B<:AbstractRadialBasis}
+    N = length(first(data))
+    dx = zeros(N)
+    dx[dim] = Δ
+
+    self = regrid(data, eval_points, basis; k=k)
+    update_weights!(self)
+
+    return if backward
+        li = regrid(data, eval_points .- Ref(dx), basis; k=k)
+        update_weights!(li)
+        w = columnwise_div(self.weights .- li.weights, Δ)
+        x -> w * x
+    else # forward difference
+        ri = regrid(data, eval_points .+ Ref(dx), basis; k=k)
+        update_weights!(ri)
+        w = columnwise_div((ri.weights .- self.weights), Δ)
+        x -> w * x
+    end
+end
+
+"""
+    function ∂virtual(data, dim, Δ, basis; k=autoselect_k(data, basis))
+
+Builds a virtual `RadialBasisOperator` whichi will be evaluated at the input points (`data`)
+where the operator is the partial derivative with respect to `dim`. Virtual operators
+interpolate the data to structured points at a distance `Δ` for which standard finite
+difference formulas can be applied.
+"""
+function ∂virtual(
+    data::AbstractVector,
+    dim,
+    Δ,
+    basis::B=PHS(3; poly_deg=2);
+    backward=true,
+    k::T=autoselect_k(data, basis),
+) where {T<:Int,B<:AbstractRadialBasis}
+    return ∂virtual(data, data, dim, Δ, basis; backward=backward, k=k)
+end
diff --git a/src/utils.jl b/src/utils.jl
index dd61b17..c47dbca 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -40,14 +40,14 @@ function autoselect_k(data::Vector, basis::B) where {B<:AbstractRadialBasis}
 end
 
 function reorder_points!(
-    x::AbstractVector{D}, adjl::Vector{Vector{T}}, k::T
-) where {D,T<:Int}
+    x::AbstractVector, adjl::AbstractVector{AbstractVector{T}}, k::T
+) where {T<:Int}
     i = symrcm(adjl, ones(T, length(x)) .* k)
     permute!(x, i)
     return nothing
 end
 
-function reorder_points!(x::AbstractVector{D}, k::T) where {D,T<:Int}
+function reorder_points!(x::AbstractVector, k::T) where {T<:Int}
     return reorder_points!(x, find_neighbors(x, k), k)
 end
 
@@ -67,3 +67,24 @@ _allocate_weights(m, n, k) = _allocate_weights(Float64, m, n, k)
 function _allocate_weights(T, m, n, k)
     return spzeros(T, m, n)
 end
+
+function columnwise_div(A::SparseMatrixCSC, B::AbstractVector)
+    I, J, V = findnz(A)
+    for idx in eachindex(V)
+        V[idx] /= B[I[idx]]
+    end
+    return sparse(I, J, V)
+end
+columnwise_div(A::SparseMatrixCSC, B::Number) = A ./ B
+
+function _find_smallest_dist(data, k)
+    tree = KDTree(data)
+    _, dists = knn(tree, data, k, true)
+    Δ = minimum(dists) do d
+        z = minimum(@view(d[2:end])) do i
+            abs(i - first(d))
+        end
+        return z
+    end
+    return Δ
+end
diff --git a/test/operators/virtual.jl b/test/operators/virtual.jl
new file mode 100644
index 0000000..5d0b405
--- /dev/null
+++ b/test/operators/virtual.jl
@@ -0,0 +1,49 @@
+using RadialBasisFunctions
+using StaticArrays
+using Statistics
+using Random
+
+rsme(test, correct) = sqrt(sum((test - correct) .^ 2) / sum(correct .^ 2))
+mean_percent_error(test, correct) = mean(abs.((test .- correct) ./ correct)) * 100
+
+f(x) = 1 + sin(4 * x[1]) + cos(3 * x[1]) + sin(2 * x[2])
+df_dx(x) = 4 * cos(4 * x[1]) - 3 * sin(3 * x[1])
+df_dy(x) = 2 * cos(2 * x[2])
+d2f_dxx(x) = -16 * sin(4 * x[1]) - 9 * cos(3 * x[1])
+d2f_dyy(x) = -4 * sin(2 * x[2])
+
+N = 10_000
+Δ = 0.001
+x = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), 1:N)
+y = f.(x)
+
+@testset "First Derivative Partials" begin
+    @testset "Polyharmonic Splines" begin
+        ∂x = ∂virtual(x, 1, Δ, PHS(3; poly_deg=2))
+        ∂y = ∂virtual(x, 2, Δ, PHS(3; poly_deg=2))
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+    end
+
+    @testset "Inverse Multiquadrics" begin
+        ∂x = ∂virtual(x, 1, Δ, IMQ(1; poly_deg=2))
+        ∂y = ∂virtual(x, 2, Δ, IMQ(1; poly_deg=2))
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+    end
+
+    @testset "Gaussian" begin
+        ∂x = ∂virtual(x, 1, Δ, Gaussian(1; poly_deg=2))
+        ∂y = ∂virtual(x, 2, Δ, Gaussian(1; poly_deg=2))
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+    end
+end
+
+@testset "Different Evaluation Points" begin
+    x2 = map(x -> SVector{2}(rand(2)), 1:100)
+    ∂x = ∂virtual(x, x2, 1, Δ, PHS(3; poly_deg=2))
+    ∂y = ∂virtual(x, x2, 2, Δ, PHS(3; poly_deg=2))
+    @test mean_percent_error(∂x(y), df_dx.(x2)) < 2
+    @test mean_percent_error(∂y(y), df_dy.(x2)) < 2
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7b391d2..078de4f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -40,6 +40,10 @@ end
     include("operators/regrid.jl")
 end
 
+@safetestset "Virtual" begin
+    include("operators/virtual.jl")
+end
+
 @safetestset "Stencil" begin
     include("linalg/stencil.jl")
 end

From 7f9ad31b1628bf9a7269b03fb698418079d899fd Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 16 Jul 2024 21:41:37 -0400
Subject: [PATCH 8/9] bump version

---
 Project.toml                | 2 +-
 docs/src/api.md             | 4 ++--
 docs/src/getting_started.md | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index f20891e..ade5c5f 100644
--- 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
index 4cfad03..671e6c2 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -1,4 +1,4 @@
-# Exported Functions
+## Exported Functions
 
 ```@autodocs
 Modules = [RadialBasisFunctions]
@@ -6,7 +6,7 @@ Private = false
 Order   = [:function, :type]
 ```
 
-# Private
+## Private
 
 ```@autodocs
 Modules = [RadialBasisFunctions]
diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
index 31d8f15..161174b 100644
--- 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.
 

From 4a96b39f0b6d12f685da240007c9b7f25dc28364 Mon Sep 17 00:00:00 2001
From: Kyle Beggs <beggskw@gmail.com>
Date: Tue, 16 Jul 2024 21:50:56 -0400
Subject: [PATCH 9/9] loosen test bounds

---
 test/operators/directional.jl |  6 +++---
 test/operators/gradient.jl    |  8 ++++----
 test/operators/laplacian.jl   |  8 ++++----
 test/operators/partial.jl     | 28 ++++++++++++++--------------
 test/operators/virtual.jl     | 16 ++++++++--------
 5 files changed, 33 insertions(+), 33 deletions(-)

diff --git a/test/operators/directional.jl b/test/operators/directional.jl
index 69a425a..0824a09 100644
--- a/test/operators/directional.jl
+++ b/test/operators/directional.jl
@@ -22,7 +22,7 @@ y = f.(x)
     v /= norm(v)
     ∇v = directional(x, v, PHS3(2))
     exact = map(x -> SVector(df_dx(x), df_dy(x)) ⋅ v, x)
-    @test mean_percent_error(∇v(y), exact) < 2
+    @test mean_percent_error(∇v(y), exact) < 5
 end
 
 @testset "Direction Vector for Each Data Center" begin
@@ -32,7 +32,7 @@ end
     end
     ∇v = directional(x, v, PHS3(2))
     exact = map((x, vv) -> SVector(df_dx(x), df_dy(x)) ⋅ vv, x, v)
-    @test mean_percent_error(∇v(y), exact) < 2
+    @test mean_percent_error(∇v(y), exact) < 5
 end
 
 @testset "Different Evaluation Points" begin
@@ -43,5 +43,5 @@ end
     end
     ∇v = directional(x, x2, v, PHS3(2))
     exact = map((x, vv) -> SVector(df_dx(x), df_dy(x)) ⋅ vv, x2, v)
-    @test mean_percent_error(∇v(y), exact) < 2
+    @test mean_percent_error(∇v(y), exact) < 5
 end
diff --git a/test/operators/gradient.jl b/test/operators/gradient.jl
index e609dea..eceeeef 100644
--- a/test/operators/gradient.jl
+++ b/test/operators/gradient.jl
@@ -19,14 +19,14 @@ y = f.(x)
 @testset "First Derivative gradients" begin
     ∇ = gradient(x, PHS(3; poly_deg=2))
     ∇y = ∇(y)
-    @test mean_percent_error(∇y[1], df_dx.(x)) < 2
-    @test mean_percent_error(∇y[2], df_dy.(x)) < 2
+    @test mean_percent_error(∇y[1], df_dx.(x)) < 5
+    @test mean_percent_error(∇y[2], df_dy.(x)) < 5
 end
 
 @testset "Different Evaluation Points" begin
     x2 = map(x -> SVector{2}(rand(2)), 1:100)
     ∇ = gradient(x, x2, PHS(3; poly_deg=2))
     ∇y = ∇(y)
-    @test mean_percent_error(∇y[1], df_dx.(x2)) < 2
-    @test mean_percent_error(∇y[2], df_dy.(x2)) < 2
+    @test mean_percent_error(∇y[1], df_dx.(x2)) < 5
+    @test mean_percent_error(∇y[2], df_dy.(x2)) < 5
 end
diff --git a/test/operators/laplacian.jl b/test/operators/laplacian.jl
index 7affb57..57b2e29 100644
--- a/test/operators/laplacian.jl
+++ b/test/operators/laplacian.jl
@@ -17,17 +17,17 @@ y = f.(x)
 
 @testset "Laplacian" begin
     ∇² = laplacian(x, PHS(3; poly_deg=4))
-    @test mean_percent_error(∇²(y), ∇²f.(x)) < 2
+    @test mean_percent_error(∇²(y), ∇²f.(x)) < 5
 
     ∇² = laplacian(x, IMQ(1; poly_deg=4))
-    @test mean_percent_error(∇²(y), ∇²f.(x)) < 2
+    @test mean_percent_error(∇²(y), ∇²f.(x)) < 5
 
     ∇² = laplacian(x, Gaussian(1; poly_deg=4))
-    @test mean_percent_error(∇²(y), ∇²f.(x)) < 2
+    @test mean_percent_error(∇²(y), ∇²f.(x)) < 5
 end
 
 @testset "Different Evaluation Points" 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
+    @test mean_percent_error(∇²(y), ∇²f.(x2)) < 5
 end
diff --git a/test/operators/partial.jl b/test/operators/partial.jl
index aaac527..86bb900 100644
--- a/test/operators/partial.jl
+++ b/test/operators/partial.jl
@@ -20,22 +20,22 @@ y = f.(x)
     @testset "Polyharmonic Splines" begin
         ∂x = partial(x, 1, 1, PHS(3; poly_deg=2))
         ∂y = partial(x, 1, 2, PHS(3; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 
     @testset "Inverse Multiquadrics" begin
         ∂x = partial(x, 1, 1, IMQ(1; poly_deg=2))
         ∂y = partial(x, 1, 2, IMQ(1; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 
     @testset "Gaussian" begin
         ∂x = partial(x, 1, 1, Gaussian(1; poly_deg=2))
         ∂y = partial(x, 1, 2, Gaussian(1; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 end
 
@@ -43,22 +43,22 @@ end
     @testset "Polyharmonic Splines" begin
         ∂2x = partial(x, 2, 1, PHS(3; poly_deg=4))
         ∂2y = partial(x, 2, 2, PHS(3; poly_deg=4))
-        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 2
-        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 2
+        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 5
+        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 5
     end
 
     @testset "Inverse Multiquadrics" begin
         ∂2x = partial(x, 2, 1, IMQ(1; poly_deg=4))
         ∂2y = partial(x, 2, 2, IMQ(1; poly_deg=4))
-        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 2
-        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 2
+        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 5
+        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 5
     end
 
     @testset "Gaussian" begin
         ∂2x = partial(x, 2, 1, Gaussian(1; poly_deg=4))
         ∂2y = partial(x, 2, 2, Gaussian(1; poly_deg=4))
-        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 2
-        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 2
+        @test mean_percent_error(∂2x(y), d2f_dxx.(x)) < 5
+        @test mean_percent_error(∂2y(y), d2f_dyy.(x)) < 5
     end
 end
 
@@ -66,6 +66,6 @@ end
     x2 = map(x -> SVector{2}(rand(2)), 1:100)
     ∂x = partial(x, x2, 1, 1, PHS(3; poly_deg=2))
     ∂y = partial(x, x2, 1, 2, PHS(3; poly_deg=2))
-    @test mean_percent_error(∂x(y), df_dx.(x2)) < 2
-    @test mean_percent_error(∂y(y), df_dy.(x2)) < 2
+    @test mean_percent_error(∂x(y), df_dx.(x2)) < 5
+    @test mean_percent_error(∂y(y), df_dy.(x2)) < 5
 end
diff --git a/test/operators/virtual.jl b/test/operators/virtual.jl
index 5d0b405..079accc 100644
--- a/test/operators/virtual.jl
+++ b/test/operators/virtual.jl
@@ -21,22 +21,22 @@ y = f.(x)
     @testset "Polyharmonic Splines" begin
         ∂x = ∂virtual(x, 1, Δ, PHS(3; poly_deg=2))
         ∂y = ∂virtual(x, 2, Δ, PHS(3; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 
     @testset "Inverse Multiquadrics" begin
         ∂x = ∂virtual(x, 1, Δ, IMQ(1; poly_deg=2))
         ∂y = ∂virtual(x, 2, Δ, IMQ(1; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 
     @testset "Gaussian" begin
         ∂x = ∂virtual(x, 1, Δ, Gaussian(1; poly_deg=2))
         ∂y = ∂virtual(x, 2, Δ, Gaussian(1; poly_deg=2))
-        @test mean_percent_error(∂x(y), df_dx.(x)) < 2
-        @test mean_percent_error(∂y(y), df_dy.(x)) < 2
+        @test mean_percent_error(∂x(y), df_dx.(x)) < 5
+        @test mean_percent_error(∂y(y), df_dy.(x)) < 5
     end
 end
 
@@ -44,6 +44,6 @@ end
     x2 = map(x -> SVector{2}(rand(2)), 1:100)
     ∂x = ∂virtual(x, x2, 1, Δ, PHS(3; poly_deg=2))
     ∂y = ∂virtual(x, x2, 2, Δ, PHS(3; poly_deg=2))
-    @test mean_percent_error(∂x(y), df_dx.(x2)) < 2
-    @test mean_percent_error(∂y(y), df_dy.(x2)) < 2
+    @test mean_percent_error(∂x(y), df_dx.(x2)) < 5
+    @test mean_percent_error(∂y(y), df_dy.(x2)) < 5
 end