add SumKernel (#34)

examples/interpolation_2d_sum_kernel.jl

@@ -0,0 +1,27 @@
+using KernelInterpolation
+using Plots
+
+# function to interpolate
+f(x) = sinpi(x[1]) + sinpi(x[2])
+
+x_min = -1.0
+x_max = 1.0
+n = 100
+nodeset = random_hypercube(n, x_min, x_max; dim = 2)
+values = f.(nodeset)
+
+proj1(x) = [x[1]]
+kernel1 = TransformationKernel{2}(Matern12Kernel{1}(), proj1)
+proj2(x) = [x[2]]
+kernel2 = TransformationKernel{2}(RadialCharacteristicKernel{1}(shape_parameter = 2.0),
+                                  proj2)
+sum_kernel = SumKernel{2}([kernel1, kernel2])
+itp = interpolate(nodeset, values, sum_kernel)
+
+many_nodes = homogeneous_hypercube(20, x_min, x_max; dim = 2)
+
+plot(layout = (1, 2))
+plot!(sum_kernel, subplot = 1)
+
+plot!(many_nodes, itp, subplot = 2)
+plot!(many_nodes, f, subplot = 2)

src/KernelInterpolation.jl

@@ -17,7 +17,7 @@ export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
        PolyharmonicSplineKernel, ThinPlateSplineKernel, WendlandKernel,
        RadialCharacteristicKernel, MaternKernel, Matern12Kernel, Matern32Kernel,
        Matern52Kernel, Matern72Kernel, RieszKernel,
-       TransformationKernel, ProductKernel
+       TransformationKernel, ProductKernel, SumKernel
 export phi, Phi, order
 export NodeSet, separation_distance, dim, values_along_dim, random_hypercube,
        random_hypercube_boundary, homogeneous_hypercube, homogeneous_hypercube_boundary,

src/interpolation.jl

@@ -54,7 +54,7 @@ dim(itp::Interpolation{Kernel, Dim, RealT, A}) where {Kernel, Dim, RealT, A} = D
 """
     coefficients(itp::Interpolation)
 
-Obtain all the coefficients of the linear combination for the interpolant, i.e. both
+Obtain all the coefficients of the linear combination for the interpolant, i.e., both
 the coefficients for the kernel part and for the polynomial part.
 
 See also [`kernel_coefficients`](@ref) and [`polynomial_coefficients`](@ref).
@@ -102,15 +102,15 @@ polyvars(itp::Interpolation) = itp.xx
 """
     order(itp)
 
-Return the order ``m`` of the polynomial used for the interpolation, i.e.
+Return the order ``m`` of the polynomial used for the interpolation, i.e.,
 the polynomial degree plus 1. If ``m = 0``, no polynomial is added.
 """
 order(itp::Interpolation) = maximum(degree.(itp.ps), init = -1) + 1
 
 @doc raw"""
     system_matrix(itp::Interpolation)
 
-Return the system matrix, i.e. the matrix
+Return the system matrix, i.e., the matrix
 ```math
 \begin{pmatrix}
 A & P \\
@@ -129,7 +129,7 @@ system_matrix(itp::Interpolation) = itp.symmetric_system_matrix
     interpolate(nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel))
 
 Interpolate the `values` evaluated at the nodes in the `nodeset` to a function using the kernel `kernel`
-and polynomials up to a degree `polynomial_degree`, i.e. determine the coefficients `c_j` and `d_k` in the expansion
+and polynomials up to a degree `polynomial_degree`, i.e., determine the coefficients `c_j` and `d_k` in the expansion
 ```math
     s(x) = \sum_{j = 1}^n c_jK(x, x_j) + \sum_{k = 1}^q d_kp_k(x),
 ```

src/kernels/radialsymmetric_kernel.jl

@@ -188,7 +188,7 @@ Thin plate spline kernel function with
 ```math
     \phi(r) = r^2\log(r),
 ```
-i.e. [`PolyharmonicSplineKernel`](@ref) with ``k = 2``.
+i.e., [`PolyharmonicSplineKernel`](@ref) with ``k = 2``.
 The thin plate spline is conditionally positive definite of order ``m = 2``.
 See Wendland (2004), p. 112.
 
@@ -357,7 +357,7 @@ order(kernel::MaternKernel) = 0
 @doc raw"""
     Matern12Kernel{Dim}(; shape_parameter = 1.0)
 
-Matern kernel with ``\nu = 1/2``, i.e.
+Matern kernel with ``\nu = 1/2``, i.e.,
 ```math
     \phi(r) = \exp(-\varepsilon r),
 ```
@@ -388,7 +388,7 @@ order(kernel::Matern12Kernel) = 0
 @doc raw"""
     Matern32Kernel{Dim}(; shape_parameter = 1.0)
 
-Matern kernel with ``\nu = 3/2``, i.e.
+Matern kernel with ``\nu = 3/2``, i.e.,
 ```math
     \phi(r) =  (1 + \sqrt{3}\varepsilon r)\exp(-\sqrt{3}\varepsilon r),
 ```
@@ -420,7 +420,7 @@ order(kernel::Matern32Kernel) = 0
 @doc raw"""
     Matern52Kernel{Dim}(; shape_parameter = 1.0)
 
-Matern kernel with ``\nu = 5/2``, i.e.
+Matern kernel with ``\nu = 5/2``, i.e.,
 ```math
     \phi(r) =  (1 + \sqrt{5}\varepsilon r + 5\cdot(\varepsilon r)^2/3)\exp(-\sqrt{5}\varepsilon r),
 ```
@@ -451,7 +451,7 @@ order(kernel::Matern52Kernel) = 0
 @doc raw"""
     Matern72Kernel{Dim}(; shape_parameter = 1.0)
 
-Matern kernel with ``\nu = 7/2``, i.e.
+Matern kernel with ``\nu = 7/2``, i.e.,
 ```math
     \phi(r) =  (1 + \sqrt{7}\varepsilon r + 12\cdot(\varepsilon r)^2/5 + 7\cdot(\varepsilon r)^3/15)\exp(-\sqrt{7}\varepsilon r),
 ```

src/kernels/special_kernel.jl

@@ -75,3 +75,47 @@ end
 
 # TODO: Is that correct in general?
 order(kernel::ProductKernel) = maximum(order.(kernel.kernels))
+
+@doc raw"""
+    SumKernel(kernels)
+
+Given a vector of `kernels`, construct a new kernel that sums the
+results of the component kernels, i.e., the new kernel ``K`` is given by
+```math
+    K(x, y) = \sum_{i = 1}^n K_i(x, y),
+```
+where ``K_i`` are the component kernels and ``n`` the number of kernels.
+Note that all component kernels need to have the same [`dim`](@ref).
+"""
+struct SumKernel{Dim} <: AbstractKernel{Dim}
+    kernels::Vector{AbstractKernel}
+
+    function SumKernel{Dim}(kernels) where {Dim}
+        @assert all(dim.(kernels) .== Dim)
+        new(kernels)
+    end
+end
+
+function (kernel::SumKernel)(x, y)
+    @assert length(x) == length(y)
+    res = 0.0
+    for k in kernel.kernels
+        res += k(x, y)
+    end
+    return res
+end
+
+function Base.show(io::IO, kernel::SumKernel{Dim}) where {Dim}
+    print(io, "SumKernel{", Dim, "}(kernels = [")
+    for (i, k) in enumerate(kernel.kernels)
+        if i < length(kernel.kernels)
+            print(io, k, ", ")
+        else
+            print(io, k)
+        end
+    end
+    print("])")
+end
+
+# TODO: Is that correct in general?
+order(kernel::SumKernel) = minimum(order.(kernel.kernels))

test/test_examples.jl

@@ -56,6 +56,11 @@ EXAMPLES_DIR = examples_dir()
                               l2=1.6308820716679489, linf=0.4234550933525697)
     end
 
+    @ki_testset "interpolation_2d_sum_kernel.jl" begin
+        @test_include_example(joinpath(EXAMPLES_DIR, "interpolation_2d_sum_kernel.jl"),
+                              l2=0.684013657976209, linf=0.16498541680450884)
+    end
+
     @ki_testset "interpolation_5d.jl" begin
         @test_include_example(joinpath(EXAMPLES_DIR, "interpolation_5d.jl"),
                               l2=0.4308925377778874, linf=0.06402624845465965)

test/test_unit.jl

@@ -137,6 +137,12 @@ using Plots
         @test_nowarn display(kernel11)
         @test order(kernel11) == 1
         @test isapprox(kernel11(x, y), kernel1(x, y) * kernel2(x, y))
+
+        kernel12 = @test_nowarn SumKernel{2}([kernel1, kernel2])
+        @test_nowarn println(kernel12)
+        @test_nowarn display(kernel12)
+        @test order(kernel12) == 0
+        @test isapprox(kernel12(x, y), kernel1(x, y) + kernel2(x, y))
     end
 
     @testset "NodeSet" begin