Remove type instabilities and allow general floating point types (#121)

Author: JoshuaLampert

CHANGELOG.md

@@ -5,6 +5,20 @@ KernelInterpolation.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.3 from v0.2.x
+
+#### Added
+
+- General floating point support ([#121]).
+
+#### Changed
+
+- The functions `random_hypersphere` and `random_hypersphere_boundary` not require a `Tuple` for
+  the argument `center`. Before, e.g., a `Vector` was allowed ([#121]).
+- The element type of `NodeSet`s will now always be converted to a floating point type, i.e., also when
+  integer values are passed. This is more consistent for an interpolation framework makes many things easier.
+  A similar approach is also used in the Meshes.jl/CoordRefSystems.jl ecosystem ([#121]).
+
 ## Changes in the v0.2 lifecycle
 
 #### Added

Project.toml

@@ -28,19 +28,19 @@ Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 KernelInterpolationMeshesExt = "Meshes"
 
 [compat]
-DiffEqCallbacks = "3, 4"
+DiffEqCallbacks = "4"
 ForwardDiff = "0.10.36"
 LinearAlgebra = "1"
 Meshes = "0.52.1, 0.53"
 Printf = "1"
 Random = "1"
 ReadVTK = "0.2"
 RecipesBase = "1.3.4"
-Reexport = "1.2"
-SciMLBase = "2.56"
+Reexport = "1.2.2"
+SciMLBase = "2.78"
 SimpleUnPack = "1.1"
 SpecialFunctions = "2"
-StaticArrays = "1.9"
+StaticArrays = "1.9.7"
 TimerOutputs = "0.5.23"
 TrixiBase = "0.1.3"
 TypedPolynomials = "0.4.1"

examples/interpolation/interpolation_2d_sphere.jl

@@ -5,7 +5,7 @@ using Plots
 f(x) = x[1] * x[2]
 
 r = 2.0
-center = [-1.0, 2.0]
+center = (-1.0, 2.0)
 n = 40
 nodeset = random_hypersphere(n, r, center)
 nodeset_boundary = random_hypersphere_boundary(20, r, center)

src/KernelInterpolation.jl

@@ -13,7 +13,7 @@ module KernelInterpolation
 
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 using ForwardDiff: ForwardDiff
-using LinearAlgebra: Symmetric, I, norm, tr, muladd, dot, diagind
+using LinearAlgebra: Symmetric, I, norm, tr, dot, diagind
 using Printf: @sprintf
 using Random: Random
 using ReadVTK: VTKFile, get_points, get_point_data, get_data

src/basis.jl

@@ -59,14 +59,14 @@ The basis functions are given by
 
 where `K` is the kernel and `x_j` are the nodes in `centers`.
 """
-struct StandardBasis{Kernel} <: AbstractBasis
-    centers::NodeSet
+struct StandardBasis{Dim, RealT, Kernel} <: AbstractBasis
+    centers::NodeSet{Dim, RealT}
     kernel::Kernel
     function StandardBasis(centers::NodeSet, kernel::Kernel) where {Kernel}
         if dim(kernel) != dim(centers)
             throw(DimensionMismatch("The dimension of the kernel and the centers must be the same"))
         end
-        new{typeof(kernel)}(centers, kernel)
+        new{dim(centers), eltype(centers), typeof(kernel)}(centers, kernel)
     end
 end
 
@@ -85,36 +85,35 @@ already includes polynomial augmentation of degree `m` defaulting to `order(kern
 which means that the [`kernel_matrix`](@ref) of this basis is the identity matrix making it suitable for interpolation. Since the
 basis already includes polynomials no additional polynomial augmentation is needed for interpolation with this basis.
 """
-struct LagrangeBasis{Kernel, I <: AbstractInterpolation, Monomials, PolyVars} <:
+struct LagrangeBasis{Dim, RealT, Kernel, I <: AbstractInterpolation, Monomials, PolyVars} <:
        AbstractBasis
-    centers::NodeSet
+    centers::NodeSet{Dim, RealT}
     kernel::Kernel
     basis_functions::Vector{I}
     ps::Monomials
     xx::PolyVars
-    function LagrangeBasis(centers::NodeSet, kernel::Kernel;
-                           m = order(kernel)) where {Kernel}
+    function LagrangeBasis(centers::NodeSet, kernel::AbstractKernel;
+                           m = order(kernel))
         if dim(kernel) != dim(centers)
             throw(DimensionMismatch("The dimension of the kernel and the centers must be the same"))
         end
+        RealT = eltype(centers)
         K = length(centers)
-        values = zeros(K)
-        values[1] = 1.0
+        values = zeros(RealT, K)
+        values[1] = one(RealT)
         b = interpolate(centers, values, kernel; m = m)
         basis_functions = Vector{typeof(b)}(undef, K)
         basis_functions[1] = b
         for i in 2:K
-            values[i - 1] = 0.0
-            values[i] = 1.0
+            values[i - 1] = zero(RealT)
+            values[i] = one(RealT)
             basis_functions[i] = interpolate(centers, values, kernel; m = m)
         end
         # All basis functions have same polynomials
         ps = first(basis_functions).ps
         xx = first(basis_functions).xx
-        new{typeof(kernel), eltype(basis_functions), typeof(ps), typeof(xx)}(centers,
-                                                                             kernel,
-                                                                             basis_functions,
-                                                                             ps, xx)
+        new{dim(centers), eltype(centers), typeof(kernel), eltype(basis_functions),
+            typeof(ps), typeof(xx)}(centers, kernel, basis_functions, ps, xx)
     end
 end
 

src/discretization.jl

@@ -176,7 +176,8 @@ function rhs!(dc, c, semi, t)
                           boundary_condition.(Ref(t), nodeset_boundary)]
         end
         # dc = -pde_boundary_matrix * c + rhs_vector
-        @trixi_timeit timer() "muladd" dc[:]=muladd(pde_boundary_matrix, -c, rhs_vector)
+        @trixi_timeit timer() "muladd" dc[:]=Base.muladd(pde_boundary_matrix, -c,
+                                                         rhs_vector)
     end
     return nothing
 end

src/interpolation.jl

@@ -159,7 +159,7 @@ Otherwise, `nodeset` is set to `centers(basis)` or `centers`.
 A regularization can be applied to the kernel matrix using the `regularization` argument, cf. [`regularize!`](@ref).
 """
 function interpolate(basis::AbstractBasis, values::Vector{RealT},
-                     nodeset::NodeSet{Dim} = centers(basis);
+                     nodeset::NodeSet{Dim, RealT} = centers(basis);
                      m = order(basis),
                      regularization = NoRegularization()) where {Dim, RealT}
     @assert dim(basis) == Dim
@@ -174,27 +174,29 @@ function interpolate(basis::AbstractBasis, values::Vector{RealT},
     else
         system_matrix = least_squares_matrix(basis, nodeset, ps, regularization)
     end
-    b = [values; zeros(q)]
+    b = [values; zeros(RealT, q)]
     c = system_matrix \ b
-    return Interpolation(basis, nodeset, c, system_matrix, ps, xx)
+    return Interpolation{typeof(basis), dim(basis), eltype(nodeset), typeof(system_matrix),
+                         typeof(ps), typeof(xx)}(basis, nodeset, c, system_matrix, ps, xx)
 end
 function interpolate(centers::NodeSet{Dim, RealT}, nodeset::NodeSet{Dim, RealT},
-                     values::AbstractVector{RealT}, kernel = GaussKernel{Dim}();
+                     values::AbstractVector{RealT}, kernel = GaussKernel{Dim, RealT}();
                      kwargs...) where {Dim, RealT}
     interpolate(StandardBasis(centers, kernel), values, nodeset; kwargs...)
 end
 
 function interpolate(centers::NodeSet{Dim, RealT},
-                     values::AbstractVector{RealT}, kernel = GaussKernel{Dim}();
+                     values::AbstractVector{RealT},
+                     kernel = GaussKernel{Dim}(; shape_parameter = RealT(1.0));
                      kwargs...) where {Dim, RealT}
     interpolate(StandardBasis(centers, kernel), values; kwargs...)
 end
 
 # Evaluate interpolant
 function (itp::Interpolation)(x)
-    s = 0
     bas = basis(itp)
     c = kernel_coefficients(itp)
+    s = zero(eltype(x))
     for j in eachindex(c)
         s += c[j] * bas[j](x)
     end

src/kernel_matrices.jl

@@ -76,7 +76,7 @@ function interpolation_matrix(basis::AbstractBasis, ps,
     regularize!(k_matrix, regularization)
     p_matrix = polynomial_matrix(centers(basis), ps)
     system_matrix = [k_matrix p_matrix
-                     p_matrix' zeros(q, q)]
+                     p_matrix' zeros(eltype(k_matrix), q, q)]
     return Symmetric(system_matrix)
 end
 
@@ -112,7 +112,7 @@ function least_squares_matrix(basis::AbstractBasis, nodeset::NodeSet, ps,
     p_matrix1 = polynomial_matrix(nodeset, ps)
     p_matrix2 = polynomial_matrix(centers(basis), ps)
     system_matrix = [k_matrix p_matrix1
-                     p_matrix2' zeros(q, q)]
+                     p_matrix2' zeros(eltype(k_matrix), q, q)]
     return system_matrix
 end
 

src/kernels/radialsymmetric_kernel.jl

@@ -278,21 +278,21 @@ function Base.show(io::IO, kernel::WendlandKernel{Dim}) where {Dim}
           kernel.shape_parameter, ", d = ", kernel.d, ")")
 end
 
-function phi(kernel::WendlandKernel, r::Real)
+function phi(kernel::WendlandKernel, r::RealT) where {RealT <: Real}
     a_r = kernel.shape_parameter * r
     if a_r >= 1
-        return 0.0
+        return RealT(0.0)
     end
-    l = floor(kernel.d / 2) + kernel.k + 1
+    l = floor(Int, kernel.d / 2) + kernel.k + 1
     if kernel.k == 0
         return (1 - a_r)^l
     elseif kernel.k == 1
         return (1 - a_r)^(l + 1) * ((l + 1) * a_r + 1)
     elseif kernel.k == 2
-        return 1 / 3 * (1 - a_r)^(l + 2) *
+        return 1 // 3 * (1 - a_r)^(l + 2) *
                ((l^2 + 4 * l + 3) * a_r^2 + (3 * l + 6) * a_r + 3)
     elseif kernel.k == 3
-        return 1 / 15 * (1 - a_r)^(l + 3) *
+        return 1 // 15 * (1 - a_r)^(l + 3) *
                ((l^3 + 9 * l^2 + 23 * l + 15) * a_r^3 + (6 * l^2 + 36 * l + 45) * a_r^2 +
                 (15 * l + 45) * a_r + 15)
     end
@@ -346,10 +346,10 @@ function Base.show(io::IO, kernel::WuKernel{Dim}) where {Dim}
           ", shape_parameter = ", kernel.shape_parameter, ")")
 end
 
-function phi(kernel::WuKernel, r::Real)
+function phi(kernel::WuKernel, r::RealT) where {RealT <: Real}
     a_r = kernel.shape_parameter * r
     if a_r >= 1
-        return 0.0
+        return RealT(0.0)
     end
     if kernel.l == 0
         # k = 0
@@ -358,29 +358,29 @@ function phi(kernel::WuKernel, r::Real)
         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)
+            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)
+            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)
+            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 *
+            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 *
+            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 *
+            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)
+            return 1 // 16 * (1 - a_r)^4 * (5 * a_r^3 + 20 * a_r^2 + 29 * a_r + 16)
         end
     end
 end
@@ -544,8 +544,8 @@ function Base.show(io::IO, kernel::Matern32Kernel{Dim}) where {Dim}
                  kernel.shape_parameter, ")")
 end
 
-function phi(kernel::Matern32Kernel, r::Real)
-    y = sqrt(3) * kernel.shape_parameter * r
+function phi(kernel::Matern32Kernel, r::RealT) where {RealT <: Real}
+    y = RealT(sqrt(3)) * kernel.shape_parameter * r
     return (1 + y) * exp(-y)
 end
 order(::Matern32Kernel) = 0
@@ -580,9 +580,9 @@ function Base.show(io::IO, kernel::Matern52Kernel{Dim}) where {Dim}
     print(io, "Matern52Kernel{", Dim, "}(shape_parameter = ", kernel.shape_parameter, ")")
 end
 
-function phi(kernel::Matern52Kernel, r::Real)
-    y = sqrt(5) * kernel.shape_parameter * r
-    return 1 / 3 * (3 + 3 * y + y^2) * exp(-y)
+function phi(kernel::Matern52Kernel, r::RealT) where {RealT <: Real}
+    y = RealT(sqrt(5)) * kernel.shape_parameter * r
+    return 1 // 3 * (3 + 3 * y + y^2) * exp(-y)
 end
 order(::Matern52Kernel) = 0
 
@@ -616,8 +616,8 @@ function Base.show(io::IO, kernel::Matern72Kernel{Dim}) where {Dim}
     print(io, "Matern72Kernel{", Dim, "}(shape_parameter = ", kernel.shape_parameter, ")")
 end
 
-function phi(kernel::Matern72Kernel, r::Real)
-    y = sqrt(7) * kernel.shape_parameter * r
+function phi(kernel::Matern72Kernel, r::RealT) where {RealT <: Real}
+    y = RealT(sqrt(7)) * kernel.shape_parameter * r
     return (1 + y + 6 * y^2 / 15 + y^3 / 15) * exp(-y)
 end
 order(::Matern72Kernel) = 0

src/kernels/special_kernel.jl

@@ -57,7 +57,7 @@ end
 
 function (kernel::ProductKernel)(x, y)
     @assert length(x) == length(y)
-    res = 1.0
+    res = eltype(x)(1.0)
     for k in kernel.kernels
         res *= k(x, y)
     end