Add smoothing RBF:s using ridge regression (#8)

* Added option to smooth RBF interpolants using ridge regression

* Bugfix, smoothing now applied after evaluateRBF!

* Refactored smoothing, added test for calling method with true

* Documentation update

src/ScatteredInterpolation.jl

@@ -27,7 +27,7 @@ function evaluate(itp::ScatteredInterpolant, points::AbstractArray{<:Real, 1})
 end
 
 """
-    interpolate(method, points, samples; metric = Euclidean(), returnRBFmatrix = false)
+    interpolate(method, points, samples; metric = Euclidean(), returnRBFmatrix = false; smooth = false)
 
 Create an interpolation of the data in `samples` sampled at the locations defined in
 `points` based on the interpolation method `method`. `metric` is any of the metrics defined
@@ -45,6 +45,11 @@ The RadialBasisFunction interpolation supports the use of unique RBF functions a
 for each sampled point by supplying `method` with a vector of interpolation methods
 of length ``k``.
 
+The RadialBasisFunction interpolation also supports smoothing of the data points using
+ridge regression. All points can be smoothed equally supplying a scalar value,
+alternatively each point can be smoothed independently by supplying a vector of smoothing
+values. Note that it is no longer interpolating when using smoothing. 
+
 The returned `ScatteredInterpolant` object can be passed to `evaluate` to interpolate the
 data to new points.
 """

src/rbf.jl

@@ -185,11 +185,17 @@ end
 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}
+                     metric = Euclidean(), returnRBFmatrix::Bool = false,
+                     smooth::Union{S, AbstractVector{S}} = false) where {N} where {S<:Number}
 
-    # Compute pairwise distances and apply the Radial Basis Function
+    #hinder smooth from being set to true and interpreted as the value 1 
+    @assert smooth != true "set the smoothing value as a number or vector of numbers"
+
+    # Compute pairwise distances, apply the Radial Basis Function
+    # and optional smoothing (ridge regression)
     A = pairwise(metric, points)
-    evaluateRBF!(A, rbf)
+
+    evaluateRBF!(A, rbf, smooth)
 
     # Solve for the weights
     itp = solveForWeights(A, points, samples, rbf, metric)
@@ -203,12 +209,50 @@ function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialB
 
 end
 
-@inline evaluateRBF!(A, rbf) = A .= rbf.(A)
+@inline function evaluateRBF!(A, rbf)
+    A .= rbf.(A)
+end
 @inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction})
     for (j, rbf) in enumerate(rbfs)
         A[:,j] .= rbf.(@view A[:,j])
     end
 end
+@inline function evaluateRBF!(A, rbf, smooth::Bool)
+    A .= rbf.(A)
+end
+@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction}, smooth::Bool)
+    for (j, rbf) in enumerate(rbfs)
+        A[:,j] .= rbf.(@view A[:,j])
+    end
+end
+@inline function evaluateRBF!(A, rbf, smooth::Vector{T}) where T <: Number
+    A .= rbf.(A)
+    for i = 1:size(A,1)
+        A[i,i] += smooth[i]
+    end
+end
+@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction}, smooth::Vector{T}) where T <: Number
+    for (j, rbf) in enumerate(rbfs)
+        A[:,j] .= rbf.(@view A[:,j])
+    end
+    for i = 1:size(A,1)
+        A[i,i] += smooth[i]
+    end
+end
+@inline function evaluateRBF!(A, rbf, smooth::T) where T <: Number
+    A .= rbf.(A)
+    for i = 1:size(A,1)
+        A[i,i] += smooth
+    end
+end
+@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction}, smooth::T) where T <: Number
+    for (j, rbf) in enumerate(rbfs)
+        A[:,j] .= rbf.(@view A[:,j])
+    end
+    for i = 1:size(A,1)
+        A[i,i] += smooth
+    end
+end
 
 @inline function solveForWeights(A, points, samples,
                                     rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction,

test/rbf.jl

@@ -71,6 +71,22 @@ radialBasisFunctions = (Gaussian(2),
         @test_throws TypeError interpolate(r, points, data; returnRBFmatrix = "true")
     end
 
+    @testset "Smooth RBF" begin
+        r = radialBasisFunctions[1]
+        @test interpolate(r, points, data; smooth = false).w ≈ interpolate(r, points, data).w
+
+
+        itp = interpolate(r, points, data; smooth = 1000.0)
+        @test itp.w ≈ interpolate(r, points, data; smooth = [1000.0 for i = 1:size(points,2)]).w
+
+        # Can not call the method with true causing true to be silently interpreted as 1 
+        @test_throws AssertionError interpolate(r, points, data; smooth = true)
+
+        # The method is no longer interpolating using smoothing
+        ev = evaluate(itp, points)
+        @test !(ev ≈ data)
+    end
+
     @testset "Generalized RBF:s" begin
         points = [1. 2; 3 4]
         @test ScatteredInterpolation.generateMultivariatePolynomial(points, 2) ≈