Merge pull request #23 from kylebeggs/feature/directional

Feature/directional

Project.toml

@@ -1,7 +1,7 @@
 name = "RadialBasisFunctions"
 uuid = "79ee0514-adf7-4479-8807-6f72ea8967e8"
 authors = ["Kyle Beggs"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

src/operators/directional.jl

@@ -28,24 +28,32 @@ function directional(
  return RadialBasisOperator(ℒ, data, basis; k=k)
 end
 
-function RadialBasisOperator(
- ℒ::Directional,
+"""
+ function directional(data, eval_points, v, basis; k=autoselect_k(data, basis))
+
+Builds a `RadialBasisOperator` where the operator is the directional derivative, `Directional`.
+"""
+function directional(
  data::AbstractVector{D},
+ eval_points::AbstractVector{D},
+ v::AbstractVector,
  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, k)
- N = length(adjl)
- weights = ntuple(_ -> _allocate_weights(TD, N, N, k), length(ℒ.ℒ))
- return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
+) 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, eval_points, basis; k=k)
 end
 
 function _update_weights!(
  op::RadialBasisOperator{<:Directional}, weights::NTuple{N,AbstractMatrix}
 ) where {N}
  v = op.ℒ.v
- @assert length(v) == N || length(v) == size(op)[2] "wrong size for v"
+ @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]

src/operators/gradient.jl

@@ -45,34 +45,6 @@ function gradient(
  return RadialBasisOperator(ℒ, data, eval_points, basis; k=k)
 end
 
-function RadialBasisOperator(
- ℒ::Gradient,
- data::AbstractVector{D},
- 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, k)
- N = length(adjl)
- weights = ntuple(_ -> _allocate_weights(TD, N, N, k), length(ℒ.ℒ))
- return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
-end
-
-function RadialBasisOperator(
- ℒ::Gradient,
- 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}
- TD = eltype(D)
- adjl = find_neighbors(data, eval_points, k)
- Na = length(adjl)
- Nd = length(data)
- weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k), length(ℒ.ℒ))
- return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
-end
-
 Base.size(op::RadialBasisOperator{<:Gradient}) = size(first(op.weights))
 
 # pretty printing

src/operators/operators.jl

@@ -47,6 +47,34 @@ function RadialBasisOperator(
  return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end
 
+function RadialBasisOperator(
+ ℒ::VectorValuedOperator,
+ data::AbstractVector{D},
+ 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, 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},
+ 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)
+ Na = length(adjl)
+ Nd = length(data)
+ weights = ntuple(_ -> _allocate_weights(TD, Na, Nd, k), length(ℒ.ℒ))
+ return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
+end
+
 # extend Base methods
 Base.length(op::RadialBasisOperator) = length(op.adjl)
 Base.size(op::RadialBasisOperator) = size(op.weights)

test/operators/directional.jl

@@ -38,3 +38,17 @@ end
  @test mean_percent_error(∇vy[1], df_dx_v) < 2
  @test mean_percent_error(∇vy[2], df_dy_v) < 2
 end
+
+@testset "Different Evaluation Points" begin
+ x2 = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), 1:N)
+ v = map(1:length(x2)) do i
+ v = SVector{2}(rand(2))
+ 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
+end