Skip to content

Commit

Permalink
Add compactly supported Wu kernels (#64)
Browse files Browse the repository at this point in the history 
* add compactly supported Wu kernels

* format

* format

* increase coverage

* fix some docs

* fix formatting of docstring
  • Loading branch information
@JoshuaLampert
JoshuaLampert authored Jul 14, 2024
1 parent 2ef5169 commit 146763a
Show file tree
Hide file tree
Showing 5 changed files with 158 additions and 13 deletions.
3 changes: 2 additions & 1 deletion docs/src/interpolation.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,7 +244,8 @@ kernels already defined, which can be used in an analogous way. For an overview
| [`PolyharmonicSplineKernel`](@ref) | ``\phi_k(r) = \begin{cases} r^k, &\text{ if } k \text{ odd}\\ r^k\log{r}, &\text{ if } k \text{ even} \end{cases}, k\in\mathbb{N}`` | ``\left\lceil{\frac{k}{2}}\right\rceil`` for odd ``k`` and ``\frac{k}{2} + 1`` for even ``k`` | ``C^{k - 1}`` for odd ``k`` and ``C^k`` for even ``k``
| [`ThinPlateSplineKernel`](@ref) | ``\phi(r) = r^2\log{r}`` | 2 | ``C^2``
| [`WendlandKernel`](@ref) | ``\phi_{d,k}(r) = \begin{cases}p_{d,k}(r), &\text{ if } 0\le r\le 1\\0, &\text{ else}\end{cases}, d, k\in\mathbb{N}`` for some polynomial ``p_{d,k}``| ``0`` | ``C^{2k}``
| [`RadialCharacteristicKernel`](@ref) | ``\phi(r) = (1 - r^2)^\beta_+, \beta\ge(d + 1)/2`` | ``0`` | ``C^0``
| [`WuKernel`](@ref) | ``\phi_{l,k}(r) = \begin{cases}p_{l,k}(r), &\text{ if } 0\le r\le 1\\0, &\text{ else}\end{cases}, l, k\in\mathbb{N}`` for some polynomial ``p_{l,k}``| ``0`` | ``C^{2(l - k)}``
| [`RadialCharacteristicKernel`](@ref) | ``\phi(r) = (1 - r)^\beta_+, \beta\ge(d + 1)/2`` | ``0`` | ``C^0``
| [`MaternKernel`](@ref) | ``\phi_{\nu}(r) = \frac{2^{1 - \nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}r\right)^\nu K_{\nu}\left(\sqrt{2\nu}r\right), \nu > 0`` | ``0`` | ``C^{2(\nu - 1)}``
| [`RieszKernel`](@ref) | ``\phi(r) = -r^\beta, 0 < \beta < 2`` | ``1`` | ``C^\infty``

Expand Down
2 changes: 1 addition & 1 deletion src/KernelInterpolation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ include("util.jl")

export get_name
export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
PolyharmonicSplineKernel, ThinPlateSplineKernel, WendlandKernel,
PolyharmonicSplineKernel, ThinPlateSplineKernel, WendlandKernel, WuKernel,
RadialCharacteristicKernel, MaternKernel, Matern12Kernel, Matern32Kernel,
Matern52Kernel, Matern72Kernel, RieszKernel,
TransformationKernel, ProductKernel, SumKernel
Expand Down
2 changes: 1 addition & 1 deletion src/kernel_matrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ end
Compute the operator matrix ``L`` discretizing ``\mathcal{L}`` for a given kernel. The operator matrix is defined as
```math
L = A_\mathcal{L} * A^{-1},
L = A_\mathcal{L} A^{-1},
```
where ``A_\mathcal{L}`` is the matrix of the differential operator (defined by the `equations`), and ``A`` the kernel matrix.
Expand Down
110 changes: 100 additions & 10 deletions src/kernels/radialsymmetric_kernel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ function Phi(kernel::RadialSymmetricKernel{Dim}, x) where {Dim}
end

function (kernel::RadialSymmetricKernel)(x, y)
@assert length(x) == length(y)
@assert length(x)==length(y) "x and y must have the same length"
return Phi(kernel, x .- y)
end

Expand Down Expand Up @@ -241,16 +241,20 @@ Wendland kernel with
0, \text{ if } \varepsilon r > 1
\end{cases},
```
where ``\varepsilon`` is the shape parameter and ``p`` is a polynomial. The
Wendland kernel is positive definite and compactly supported.
See Wendland (2004), p. 129.
where ``\varepsilon`` is the shape parameter and ``p`` is a polynomial with
minimal degree. The Wendland kernel is positive definite for `d\le Dim` and
compactly supported. See Wendland (2004), p. 129 or Fasshauer (2007), pp. 87--88.
See also [`RadialSymmetricKernel`](@ref).
- Holger Wendland (2004)
Scattered Data Approximation
Cambridge University Press
[DOI: 10.1017/CBO9780511617539](https://doi.org/10.1017/CBO9780511617539)
- Gregory Fasshauer (2007)
Meshfree Approximation Methods with MATLAB
World Scientific
[DOI: 10.1142/6437](https://doi.org/10.1142/6437)
"""
struct WendlandKernel{Dim, RealT} <: RadialSymmetricKernel{Dim}
k::Int
Expand All @@ -259,8 +263,8 @@ struct WendlandKernel{Dim, RealT} <: RadialSymmetricKernel{Dim}
end

function WendlandKernel{Dim}(k::Int; shape_parameter = 1.0, d::Int = Dim) where {Dim}
@assert d <= Dim
@assert k in 0:3
@assert d<=Dim "d has to be smaller or equal to Dim"
@assert k in 0:3 "kernel only implemented for k in 0:3"
WendlandKernel{Dim, typeof(shape_parameter)}(k, shape_parameter, d)
end

Expand Down Expand Up @@ -295,12 +299,99 @@ function phi(kernel::WendlandKernel, r::Real)
end
order(::WendlandKernel) = 0

@doc raw"""
WuKernel{Dim}(l, k; shape_parameter = 1.0)
Wu kernel with
```math
\phi_{l,k}(r) = \begin{cases}
p_{l,k}(\varepsilon r), \text{ if } 0\le \varepsilon r\le 1\\
0, \text{ if } \varepsilon r > 1
\end{cases},
```
where ``\varepsilon`` is the shape parameter, ``k\le l``, and ``p`` is a polynomial
of degree ``4l - 2k + 1``. The Wu kernel is positive definite for ``Dim\le 2k + 1``
and compactly supported. See Fasshauer (2007), pp. 88--90 and Wu (1995).
See also [`RadialSymmetricKernel`](@ref).
- Gregory Fasshauer (2007)
Meshfree Approximation Methods with MATLAB
World Scientific
[DOI: 10.1142/6437](https://doi.org/10.1142/6437)
- Zongmin Wu (1995)
Compactly supported positive definite radial functions
Advances in Computational Mathematics
[DOI: 10.1007/BF03177517](https://doi.org/10.1007/BF03177517)
"""
struct WuKernel{Dim, RealT} <: RadialSymmetricKernel{Dim}
l::Int
k::Int
shape_parameter::RealT
end

function WuKernel{Dim}(l::Int, k::Int; shape_parameter = 1.0) where {Dim}
@assert l>=k "l has to be bigger or equal to k"
@assert l in 0:3 "kernel only implemented for l in 0:3"
WuKernel{Dim, typeof(shape_parameter)}(l, k, shape_parameter)
end

function get_name(kernel::WuKernel)
string(nameof(typeof(kernel))) * string(kernel.l) * "," * string(kernel.k) * "{" *
string(dim(kernel)) * "}"
end

function Base.show(io::IO, kernel::WuKernel{Dim}) where {Dim}
print(io, "WuKernel{", Dim, "}(l = ", kernel.l, ", k = ", kernel.k,
", shape_parameter = ", kernel.shape_parameter, ")")
end

function phi(kernel::WuKernel, r::Real)
a_r = kernel.shape_parameter * r
if a_r >= 1
return 0.0
end
if kernel.l == 0
# k = 0
return 1 - a_r
elseif kernel.l == 1
if kernel.k == 0
return (1 - a_r)^3 * (a_r^2 + 3 * a_r + 1)
elseif kernel.k == 1
return 1 / 2 * (1 - a_r)^2 * (a_r + 2)
end
elseif kernel.l == 2
if kernel.k == 0
return (1 - a_r)^5 * (a_r^4 + 5 * a_r^3 + 9 * a_r^2 + 5 * a_r + 1)
elseif kernel.k == 1
return 1 / 4 * (1 - a_r)^4 * (3 * a_r^3 + 12 * a_r^2 + 16 * a_r + 4)
elseif kernel.k == 2
return 1 / 8 * (1 - a_r)^3 * (3 * a_r^2 + 9 * a_r + 8)
end
elseif kernel.l == 3
if kernel.k == 0
return 1 / 5 * (1 - a_r)^7 *
(5 * a_r^6 + 35 * a_r^5 + 101 * a_r^4 + 147 * a_r^3 + 101 * a_r^2 +
35 * a_r + 5)
elseif kernel.k == 1
return 1 / 6 * (1 - a_r)^6 *
(5 * a_r^5 + 30 * a_r^4 + 72 * a_r^3 + 82 * a_r^2 + 36 * a_r + 6)
elseif kernel.k == 2
return 1 / 8 * (1 - a_r)^5 *
(5 * a_r^4 + 25 * a_r^3 + 48 * a_r^2 + 40 * a_r + 8)
elseif kernel.k == 3
return 1 / 16 * (1 - a_r)^4 * (5 * a_r^3 + 20 * a_r^2 + 29 * a_r + 16)
end
end
end
order(::WuKernel) = 0

@doc raw"""
RadialCharacteristicKernel{Dim}(beta = 2.0; shape_parameter = 1.0)
Radial characteristic function kernel function with
Radial characteristic function (or also called truncated power or Askey) kernel function with
```math
\phi(r) = (1 - (\varepsilon r)^2)^\beta_+,
\phi(r) = (1 - \varepsilon r)^\beta_+,
```
where ``\varepsilon`` is the shape parameter. The radial characteristic function is
positive definite if ``\beta\ge (d + 1)/2``. It is compactly supported.
Expand Down Expand Up @@ -549,8 +640,7 @@ struct RieszKernel{Dim, RealT} <: RadialSymmetricKernel{Dim}
end

function RieszKernel{Dim}(beta; shape_parameter = 1.0) where {Dim}
@assert beta > 0
@assert beta < 2
@assert 0<beta<2 "beta has to be in (0, 2)"
RieszKernel{Dim, typeof(shape_parameter)}(beta, shape_parameter)
end

Expand Down
54 changes: 54 additions & 0 deletions test/test_unit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ include("test_util.jl")
kernel1 = @test_nowarn GaussKernel{2}(shape_parameter = 2.0)
@test_nowarn println(kernel1)
@test_nowarn display(kernel1)
@test_nowarn get_name(kernel1)
@test dim(kernel1) == 2
@test order(kernel1) == 0
@test isapprox(phi(kernel1, 0.0), 1.0)
Expand All @@ -32,6 +33,7 @@ include("test_util.jl")
kernel2 = @test_nowarn MultiquadricKernel{2}()
@test_nowarn println(kernel2)
@test_nowarn display(kernel2)
@test_nowarn get_name(kernel2)
@test order(kernel2) == 1
@test isapprox(phi(kernel2, 0.0), 1.0)
@test isapprox(phi(kernel2, 0.5), 1.118033988749895)
Expand All @@ -40,6 +42,7 @@ include("test_util.jl")
kernel3 = @test_nowarn InverseMultiquadricKernel{2}()
@test_nowarn println(kernel3)
@test_nowarn display(kernel3)
@test_nowarn get_name(kernel3)
@test order(kernel3) == 0
@test isapprox(phi(kernel3, 0.0), 1.0)
@test isapprox(phi(kernel3, 0.5), 0.8944271909999159)
Expand All @@ -48,13 +51,15 @@ include("test_util.jl")
kernel4 = @test_nowarn PolyharmonicSplineKernel{2}(3)
@test_nowarn println(kernel4)
@test_nowarn display(kernel4)
@test_nowarn get_name(kernel4)
@test order(kernel4) == 2
@test isapprox(phi(kernel4, 0.5), 0.125)
@test isapprox(kernel4(x, y), 0.02778220597956396)

kernel5 = @test_nowarn ThinPlateSplineKernel{2}()
@test_nowarn println(kernel5)
@test_nowarn display(kernel5)
@test_nowarn get_name(kernel5)
@test order(kernel5) == 2
@test isapprox(phi(kernel5, 0.5), -0.17328679513998632)
@test isapprox(kernel5(x, y), -0.10956712895893082)
Expand All @@ -72,15 +77,58 @@ include("test_util.jl")
kernel6 = @test_nowarn WendlandKernel{2}(k)
@test_nowarn println(kernel6)
@test_nowarn display(kernel6)
@test_nowarn get_name(kernel6)
@test order(kernel6) == 0
@test isapprox(phi(kernel6, 0.0), 1.0)
@test isapprox(phi(kernel6, 1.1), 0.0)
@test isapprox(phi(kernel6, 0.5), expected_values[k + 1])
@test isapprox(kernel6(x, y), expected_differences[k + 1])
end

expected_values = [
0.5,
0.34375,
0.3125,
0.201171875,
0.240234375,
0.20703125,
0.1150146484375,
0.14461263020833331,
0.169677734375,
0.14111328125
]
expected_differences = [
0.6971304755630894,
0.6777128254016545,
0.5595868163344161,
0.5741858708746038,
0.5922403769292889,
0.46589172653038635,
0.47839230403264094,
0.5106135476269533,
0.5197783607481103,
0.39497468985502254
]
i = 1
for l in 0:3
for k in 0:l
kernel6_1 = @test_nowarn WuKernel{2}(l, k)
@test_nowarn println(kernel6_1)
@test_nowarn display(kernel6_1)
@test_nowarn get_name(kernel6_1)
@test order(kernel6_1) == 0
@test isapprox(phi(kernel6_1, 0.0), 1.0)
@test isapprox(phi(kernel6_1, 1.1), 0.0)
@test isapprox(phi(kernel6_1, 0.5), expected_values[i])
@test isapprox(kernel6_1(x, y), expected_differences[i])
i += 1
end
end

kernel7 = @test_nowarn RadialCharacteristicKernel{2}()
@test_nowarn println(kernel7)
@test_nowarn display(kernel7)
@test_nowarn get_name(kernel7)
@test order(kernel7) == 0
@test isapprox(phi(kernel7, 0.0), 1.0)
@test isapprox(phi(kernel7, 0.5), 0.25)
Expand All @@ -105,6 +153,7 @@ include("test_util.jl")
kernel8 = @test_nowarn MaternKernel{2}(nus[i])
@test_nowarn println(kernel8)
@test_nowarn display(kernel8)
@test_nowarn get_name(kernel8)
@test order(kernel8) == 0
@test isapprox(phi(kernel8, 0.0), 1.0)
@test isapprox(phi(kernel8, 0.5), expected_values[i])
Expand All @@ -113,6 +162,7 @@ include("test_util.jl")
kernel8_1 = @test_nowarn kernels[i]{2}()
@test_nowarn println(kernel8_1)
@test_nowarn display(kernel8_1)
@test_nowarn get_name(kernel8_1)
@test order(kernel8_1) == 0

@test isapprox(phi(kernel8, 0.5), phi(kernel8_1, 0.5))
Expand All @@ -122,6 +172,7 @@ include("test_util.jl")
kernel9 = @test_nowarn RieszKernel{2}(1.1)
@test_nowarn println(kernel9)
@test_nowarn display(kernel9)
@test_nowarn get_name(kernel9)
@test order(kernel9) == 1
@test isapprox(phi(kernel9, 0.0), 0.0)
@test isapprox(phi(kernel9, 0.5), -0.4665164957684037)
Expand All @@ -131,6 +182,7 @@ include("test_util.jl")
kernel10 = @test_nowarn TransformationKernel{3}(kernel1, trafo)
@test_nowarn println(kernel10)
@test_nowarn display(kernel10)
@test_nowarn get_name(kernel10)
@test order(kernel10) == 0
x3 = [-1.0, 2.0, pi / 8]
y3 = [2.3, 4.2, -12.3]
Expand All @@ -139,13 +191,15 @@ include("test_util.jl")
kernel11 = @test_nowarn ProductKernel{2}([kernel1, kernel2])
@test_nowarn println(kernel11)
@test_nowarn display(kernel11)
@test_nowarn get_name(kernel11)
@test order(kernel11) == 1
@test isapprox(kernel11(x, y), kernel1(x, y) * kernel2(x, y))
@test isapprox((kernel1 * kernel2)(x, y), kernel1(x, y) * kernel2(x, y))

kernel12 = @test_nowarn SumKernel{2}([kernel1, kernel2])
@test_nowarn println(kernel12)
@test_nowarn display(kernel12)
@test_nowarn get_name(kernel12)
@test order(kernel12) == 0
@test isapprox(kernel12(x, y), kernel1(x, y) + kernel2(x, y))
@test isapprox((kernel1 + kernel2)(x, y), kernel1(x, y) + kernel2(x, y))
Expand Down

0 comments on commit 146763a

Please sign in to comment.