Add generalized RBFs and refactor the RBF code (#7)

* Add generalized RBFs.

* Drop 0.6 support.

* Relax type constraints and refactor interpolation/evaluation.

* Relax a few more type constraints.

* Iteration fix.

* Add documentation and tests.

* Tweak tests to not raise SingularException on OSX.

* Another test tweak to pass on Openblas.

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 0.6
   - 1.0
   - nightly
 matrix:

REQUIRE

@@ -1,3 +1,4 @@
-julia 0.6
+julia 0.7
 Distances
-NearestNeighbors
+NearestNeighbors
+Combinatorics

docs/src/api.md

@@ -21,6 +21,8 @@ Gaussian
 InverseQuadratic
 Polyharmonic
 ThinPlate
+GeneralizedMultiquadratic
+GeneralizedPolyharmonic
 ```
 
 ### Inverse Distance Weighting (Shepard)

docs/src/methods.md

@@ -15,6 +15,15 @@ where ``||\mathbf{x} - \mathbf{x_i}||`` is the distance between ``\mathbf{x}`` a
 To use radial basis function interpolation, pass one of the available basis functions as 
 `method` to `interpolate`.
 
+If a `GeneralizedRadialBasisFunction` is used, an additional polynomial term is added in
+order for the resulting matrix to be positive definite:
+```math
+u(\mathbf{x}) = \displaystyle \sum_{i = 1}^{N}{ w_i \phi(||\mathbf{x} - \mathbf{x_i}||)} + \mathbf{P}(\mathbf{x})\mathbf{λ}
+```
+where ``\mathbf{P}(\mathbf{x})`` is the matrix defining a complete homogeneous symmetric
+polynomial of degree `degree`, and ``\mathbf{λ}`` is a vector containing the polynomial
+coefficients.
+
 ### Available basis functions
 
   * [`Multiquadratic`](@ref)
@@ -54,10 +63,40 @@ To use radial basis function interpolation, pass one of the available basis func
     ```
 
   * [`ThinPlate`](@ref) spline
-    
+
     A thin plate spline is the special case ``k = 2`` of the polyharmonic splines.
     `ThinPlate()` is a shorthand for `Polyharmonic(2)`.
 
+  * [`GeneralizedMultiquadratic`](@ref)
+
+    ```math
+    ϕ(r) = \left(1 + (ɛr)^2\right)^\beta
+    ```
+    The generalzized multiquadratic results in a positive definite system for polynomials of
+    `degree` ``m \geq \lceil\beta\rceil``.
+
+  * [`GeneralizedPolyharmonic`](@ref) spline
+
+    ```math
+    ϕ(r) = 
+    \begin{cases}
+        \begin{align*}
+            &r^k                    &   k = 1, 3, 5, ... \\
+            &r^k \mathrm{ln}(r)     &   k = 2, 4, 6, ...
+        \end{align*}
+    \end{cases}
+    ```
+    The generalized polyharmonic spline results in a positive definite system for
+    polynomials of `degree`
+    ```math
+    \begin{cases}
+        \begin{align*}
+            m &\geq \left\lceil\frac{\beta}{2}\right\rceil          &   k = 1, 3, 5, ... \\
+            m &= \beta + 1                                          &   k = 2, 4, 6, ...
+        \end{align*}
+    \end{cases}
+    ```
+
 ## Inverse Distance Weighting
 Also called Shepard interpolation, the basic version computes the interpolated value at
 some point ``\mathbf{x}`` by

src/ScatteredInterpolation.jl

@@ -1,7 +1,9 @@
 module ScatteredInterpolation
 
 using Distances,
-    NearestNeighbors
+    NearestNeighbors,
+    Combinatorics,
+    LinearAlgebra
 
 export interpolate,
     evaluate

src/rbf.jl

@@ -1,11 +1,18 @@
-abstract type RadialBasisFunction <: InterpolationMethod end
+abstract type AbstractRadialBasisFunction <: InterpolationMethod end
+abstract type RadialBasisFunction <: AbstractRadialBasisFunction end
+abstract type GeneralizedRadialBasisFunction <: AbstractRadialBasisFunction end
+
+Base.iterate(x::AbstractRadialBasisFunction) = (x, nothing)
+Base.iterate(x::AbstractRadialBasisFunction, ::Any) = nothing
 
 export  Gaussian,
         Multiquadratic,
         InverseQuadratic,
         InverseMultiquadratic,
         Polyharmonic,
-        ThinPlate
+        ThinPlate,
+        GeneralizedMultiquadratic,
+        GeneralizedPolyharmonic
 
 # Define types for the different kinds of radial basis functions
 for rbf in (:Gaussian,
@@ -23,6 +30,18 @@ for rbf in (:Gaussian,
     end
 end
 
+# Generalized RBF:s
+struct GeneralizedMultiquadratic{T<:Real, S<:Real, U<:Integer} <: GeneralizedRadialBasisFunction
+    ε::T
+    β::S
+    degree::U
+end
+
+struct GeneralizedPolyharmonic{S<:Integer, U<:Integer} <: GeneralizedRadialBasisFunction
+    k::S
+    degree::U
+end
+
 @doc "
     Gaussian(ɛ = 1)
 
@@ -104,81 +123,162 @@ This is a shorthand for `Polyharmonic(2)`.
 " ThinPlate
 ThinPlate() = Polyharmonic(2)
 
+@doc "
+    GeneralizedMultiquadratic(ɛ, β, degree)
+
+Define a generalized Multiquadratic Radial Basis Function
+
+```math
+ϕ(r) = (1 + (ɛ*r)^2)^β
+```
+Results in a positive definite system for a 'degree' of ⌈β⌉ or higher.
+
+" GeneralizedMultiquadratic
+(rbf::GeneralizedMultiquadratic)(r) = (1 + (rbf.ε*r)^2)^rbf.β
+
+
+@doc "
+    GeneralizedPolyharmonic(k, degree)
+
+Define a generalized Polyharmonic Radial Basis Function
+
+```math
+ϕ(r) = r^k, k = 1, 3, 5, ...
+\\\\
+ϕ(r) = r^k ln(r), k = 2, 4, 6, ...
+```
+Results in a positive definite system for a 'degree' of ⌈k/2⌉ or higher for k = 1, 3, 5, ...
+and of exactly k + 1 for k = 2, 4, 6, ...
+" GeneralizedPolyharmonic
+function (rbf::GeneralizedPolyharmonic)(r)
+    @assert rbf.k > 0
+
+    # Distinguish odd and even cases
+    expr = if rbf.k % 2 == 0
+        (r > 0 ? r^rbf.k*log(r) : 0.0)
+    else
+        (r^rbf.k)
+    end
+
+    expr
+end
+
+abstract type RadialBasisInterpolant <: ScatteredInterpolant end
+
+struct RBFInterpolant{F, T, A, N, M} <: RadialBasisInterpolant where A <: AbstractArray{<:Real, 2}
+
+    w::Array{T,N}
+    points::A
+    rbf::F
+    metric::M
+end
 
-struct RBFInterpolant{F, T, N, M} <: ScatteredInterpolant
+struct GeneralizedRBFInterpolant{F, T, A, N, M} <: RadialBasisInterpolant where A <: AbstractArray{<:Real, 2}
 
     w::Array{T,N}
-    points::Array{T,2}
+    λ::Array{T,N}
+    points::A
     rbf::F
     metric::M
 end
 
-function interpolate(rbf::RadialBasisFunction,
+function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialBasisFunction,
                      points::AbstractArray{<:Real,2},
                      samples::AbstractArray{<:Number,N};
                      metric = Euclidean(), returnRBFmatrix::Bool = false) where {N}
 
     # Compute pairwise distances and apply the Radial Basis Function
     A = pairwise(metric, points)
-    @. A = rbf(A)
+    evaluateRBF!(A, rbf)
 
     # Solve for the weights
-    w = A\samples
+    itp = solveForWeights(A, points, samples, rbf, metric)
 
     # Create and return an interpolation object
     if returnRBFmatrix    # Return matrix A
-        return RBFInterpolant(w, points, rbf, metric), A
+        return itp, A
     else
-        return RBFInterpolant(w, points, rbf, metric)
+        return itp
     end
 
 end
 
-function interpolate(rbfs::Vector{T} where T <: ScatteredInterpolation.RadialBasisFunction,
-                     points::AbstractArray{<:Real,2},
-                     samples::AbstractArray{<:Number,N};
-                     metric = Euclidean(), returnRBFmatrix::Bool = false) where {N}
-
-    # Compute pairwise distances and apply the Radial Basis Function
-    A = pairwise(metric, points)
-
-    n = size(points,2)
+@inline evaluateRBF!(A, rbf) = A .= rbf.(A)
+@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction})
     for (j, rbf) in enumerate(rbfs)
-        for i = 1:n
-            A[i,j] = rbf(A[i,j])
-        end
+        A[:,j] .= rbf.(@view A[:,j])
     end
+end
 
-    # Solve for the weights
+@inline function solveForWeights(A, points, samples,
+                                    rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction,
+                                    metric)
     w = A\samples
-
-    # Create and return an interpolation object
-    if returnRBFmatrix    # Return matrix A
-        return RBFInterpolant(w, points, rbfs, metric), A
-    else
-        return RBFInterpolant(w, points, rbfs, metric)
-    end
-
+    RBFInterpolant(w, points, rbf, metric)
+end
+@inline function solveForWeights(A, points, samples,
+                                    rbf::Union{T, AbstractVector{T}} where T <: Union{GeneralizedRadialBasisFunction, RadialBasisFunction},
+                                    metric)
+    # Use the maximum degree among the generalized RBF:s
+    P = getPolynomial(rbf, points)
+
+    # Solve for the weights and polynomial coefficients
+    # We end up with a blocked system, so we don't have to form the full matrix
+    Af = factorize(A)
+    B = -P'*(Af\P)
+    E = B\(P'*(Af\samples))
+
+    w = A\(I*samples + P*E)
+    λ = -E
+
+    GeneralizedRBFInterpolant(w, λ, points, rbf, metric)
 end
 
-function evaluate(itp::RBFInterpolant, points::AbstractArray{<:Real, 2})
+function evaluate(itp::RadialBasisInterpolant, points::AbstractArray{<:Real, 2})
 
-    # Compute distance matrix
+    # Compute distance matrix and evaluate the RBF
     A = pairwise(itp.metric, points, itp.points)
-    @. A = itp.rbf(A)
+    evaluateRBF!(A, itp.rbf)
+
+    # Compute polynomial matrix for generalized RBF:s
+    P = getPolynomial(itp.rbf, points)
 
     # Compute the interpolated values
-    return A*itp.w
+    return computeInterpolatedValues(A, P, itp)
 end
 
-function evaluate(itp::RBFInterpolant{S,T,U,V}, points::AbstractArray{<:Real, 2}) where {S <: Vector, T, U, V}
+@inline getPolynomial(rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction, points) = nothing
+@inline function getPolynomial(rbf::Union{T, AbstractVector{T}} where T <: Union{GeneralizedRadialBasisFunction, RadialBasisFunction}, points)
 
-    # Compute distance matrix
-    A = pairwise(itp.metric, points, itp.points)
-    for (j, rbf) in enumerate(itp.rbf)
-        @. A[:,j] = rbf(A[:,j])
+    # Use the maximum degree among the generalized RBF:s
+    degree = maximum(x isa GeneralizedRadialBasisFunction ? x.degree : 0 for x in rbf)
+    P = generateMultivariatePolynomial(points, degree)
+end
+
+@inline computeInterpolatedValues(A, P, itp::RBFInterpolant) = A*itp.w
+@inline computeInterpolatedValues(A, P, itp::GeneralizedRBFInterpolant) = A*itp.w + P*itp.λ
+
+# Helper function to generate matrices defining complete homogenic symmetric polynomials
+function generateMultivariatePolynomial(points::AbstractArray{<:Real, 2}, degree::Integer)
+
+    # How big should the matrix be?
+    nDimensions, nPoints = size(points)
+    nTerms = binomial(degree + nDimensions, degree)
+
+    P = ones(nPoints, nTerms)
+
+    # Start with the lowest orders and work upwards
+    position = 2
+    while position <= nTerms
+        for order = 1:degree
+            for combination in with_replacement_combinations(1:nDimensions, order)
+                for var in combination
+                    P[:, position] .*= points[var, :]
+                end
+                position += 1
+            end
+        end
     end
 
-    # Compute the interpolated values
-    return A*itp.w
-end
+    P
+end

test/rbf.jl

@@ -1,9 +1,16 @@
 
 # Define some points and data in 2D
-points = permutedims([0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0], (2,1))
-data = [0.0; 0.5; 0.5; 1.0]
-
-radialBasisFunctions = (Gaussian(2), Multiquadratic(2), InverseQuadratic(2), InverseMultiquadratic(2), Polyharmonic(2), ThinPlate())
+points = permutedims([0.0 0.0; 0.0 1.0; 0.5 0.5; 1.0 0.0; 1.0 1.0], (2,1))
+data = [0.0; 0.5; 0.5; 0.5; 1.0]
+
+radialBasisFunctions = (Gaussian(2),
+                        Multiquadratic(2),
+                        InverseQuadratic(2),
+                        InverseMultiquadratic(2),
+                        Polyharmonic(2),
+                        ThinPlate(),
+                        GeneralizedMultiquadratic(1, 1/2, 2),
+                        GeneralizedPolyharmonic(3, 2))
 
 @testset "RBF" begin
 
@@ -23,6 +30,10 @@ radialBasisFunctions = (Gaussian(2), Multiquadratic(2), InverseQuadratic(2), Inv
             1/sqrt(1 + (2x)^2)
         elseif isa(r, Polyharmonic)
             x > 0.0 ? x^2*log(x) : 0.0
+        elseif isa(r, GeneralizedMultiquadratic)
+            (1 + (1x)^2)^(1/2)
+        elseif isa(r, GeneralizedPolyharmonic)
+            x^3
         end
 
         @test r.(data) ≈ f.(data)
@@ -60,4 +71,11 @@ radialBasisFunctions = (Gaussian(2), Multiquadratic(2), InverseQuadratic(2), Inv
         @test_throws TypeError interpolate(r, points, data; returnRBFmatrix = "true")
     end
 
+    @testset "Generalized RBF:s" begin
+        points = [1. 2; 3 4]
+        @test ScatteredInterpolation.generateMultivariatePolynomial(points, 2) ≈ 
+                                                                            [1 1 3 1 3 9;
+                                                                             1 2 4 4 8 16]
+    end
+
 end