Different RBFs per point and type stability (#5)

* added vector functionality to rbf

* Changed type signature for type stability

* Added functionality for different RBFs per point. Improved type stability

src/ScatteredInterpolation.jl

@@ -39,6 +39,10 @@ interpolations.
 `samples` is an ``k×m`` array, where ``k`` is the number of sampled points (same as for
 `points`) and ``m`` is the dimension of the sampled data.
 
+The RadialBasisFunction interpolation supports the use of unique RBF functions and widths
+for each sampled point by supplying `method` with a vector of interpolation methods
+of length ``k``.
+
 The returned `ScatteredInterpolant` object can be passed to `evaluate` to interpolate the
 data to new points.
 """

src/idw.jl

@@ -1,4 +1,4 @@
-@compat abstract type ShepardType <: InterpolationMethod end
+abstract type ShepardType <: InterpolationMethod end
 
 export Shepard
 
@@ -7,8 +7,10 @@ export Shepard
 
 Standard Shepard interpolation with power parameter `P`.
 """
-struct Shepard{P} <: ShepardType end
-Shepard(P::Real = 2) = Shepard{P}()
+struct Shepard{T} <: ShepardType where T <: Real
+    P::T
+end
+Shepard() = Shepard(2)
 
 struct ShepardInterpolant{F, T1, T2, N, M} <: ScatteredInterpolant
 
@@ -54,12 +56,12 @@ function evaluate(itp::ShepardInterpolant, points::AbstractArray{<:Real,2})
 end
 
 # Original Shepard
-function evaluatePoint(::Shepard{P},
+function evaluatePoint(idw::Shepard,
                        dataPoints::AbstractArray{<:Real,2},
                        data::AbstractArray{<:Number,N},
-                       d::AbstractVector) where {N, P}
+                       d::AbstractVector) where {N}
 
     # Compute weigths and return the weighted sum
-    w = d.^P
+    w = d.^idw.P
     value = sum(w.*data, 1)./sum(w)
 end

src/nearestNeighbor.jl

@@ -9,9 +9,8 @@ Nearest neighbor interpolation.
 """
 struct NearestNeighbor <: InterpolationMethod end
 
-struct NearestNeighborInterpolant{T, TT, N} <: ScatteredInterpolant
-
-    data::Array{T,N}
+struct NearestNeighborInterpolant{T, TT} <: ScatteredInterpolant
+    data::Array{T}
     tree::TT
 end
 

src/rbf.jl

@@ -14,82 +14,60 @@ for rbf in (:Gaussian,
             :InverseMultiquadratic,
             :Polyharmonic)
     @eval begin
-        struct $rbf{T} <: RadialBasisFunction
+        struct $rbf{T} <: RadialBasisFunction where T <: Real
+            ɛ::T    
         end
 
-        # Define constructors
-        function $rbf(c::Real = 1)
-            $rbf{c}()
-        end
+        # Define default constructors
+        $rbf() = $rbf(1)
     end
 end
 
-"""
-    ThinPlate()
-
-Define a Thin Plate Spline Radial Basis Function
-
-```math
-ϕ(r) = r^2 ln(r)
-```
-
-This is a shorthand for `Polyharmonic(2)`.
-"""
-const ThinPlate = Polyharmonic{2}
-
-"""
+@doc "
     Gaussian(ɛ = 1)
 
 Define a Gaussian Radial Basis Function
 
 ```math
 ϕ(r) = e^{-(ɛr)^2}
 ```
-"""
-@generated function (::Gaussian{C})(r::Real) where {C}
-    :(exp(-($C*r)^2))
-end
+" Gaussian
+(rbf::Gaussian)(r) = exp(-(rbf.ɛ*r)^2) 
 
-"""
+@doc "
     Multiquadratic(ɛ = 1)
 
 Define a Multiquadratic Radial Basis Function
 
 ```math
 ϕ(r) = \\sqrt{1 + (ɛr)^2}
 ```
-"""
-@generated function (::Multiquadratic{C})(r::Real) where {C}
-    :(sqrt(1 + ($C*r)^2))
-end
+" Multiquadratic
+(rbf::Multiquadratic)(r) = (sqrt(1 + (rbf.ɛ*r)^2))
 
-"""
+@doc "
     InverseQuadratic(ɛ = 1)
 
 Define an Inverse Quadratic Radial Basis Function
 
 ```math
 ϕ(r) = \\frac{1}{1 + (ɛr)^2}
 ```
-"""
-@generated function (::InverseQuadratic{C})(r::Real) where {C}
-    :(1/(1 + ($C*r)^2))
-end
+" InverseQuadratic
+(rbf::InverseQuadratic)(r) = (1/(1 + (rbf.ɛ*r)^2))
 
-"""
+@doc "
     InverseMultiquadratic(ɛ = 1)
 
 Define an Inverse Multiquadratic Radial Basis Function
 
 ```math
 ϕ(r) = \\frac{1}{\\sqrt{1 + (ɛr)^2}}
 ```
-"""
-@generated function (::InverseMultiquadratic{C})(r::Real) where {C}
-    :(1/sqrt(1 + ($C*r)^2))
-end
+" InverseMultiquadratic
+(rbf::InverseMultiquadratic)(r) = (1/sqrt(1 + (rbf.ɛ*r)^2))
 
-"""
+@doc "
     Polyharmonic(k = 1)
 
 Define a Polyharmonic Spline Radial Basis Function
@@ -99,21 +77,34 @@ Define a Polyharmonic Spline Radial Basis Function
 \\\\
 ϕ(r) = r^k ln(r), k = 2, 4, 6, ...
 ```
-"""
-@generated function (::Polyharmonic{C})(r::Real) where {C}
-
-    @assert typeof(C) <: Integer && C > 0
+" Polyharmonic
+function (rbf::Polyharmonic{T})(r) where T <: Integer
+    @assert rbf.ɛ > 0
 
     # Distinguish odd and even cases
-    expr = if C % 2 == 0
-        :(r > 0 ? r^$C*log(r) : 0.0)
+    expr = if rbf.ɛ % 2 == 0
+        (r > 0 ? r^rbf.ɛ*log(r) : 0.0)
     else
-        :(r^$C)
+        (r^rbf.ɛ)
     end
 
     expr
 end
 
+@doc "
+    ThinPlate()
+
+Define a Thin Plate Spline Radial Basis Function
+
+```math
+ϕ(r) = r^2 ln(r)
+``` 
+
+This is a shorthand for `Polyharmonic(2)`.
+" ThinPlate
+ThinPlate() = Polyharmonic(2)
+
+
 struct RBFInterpolant{F, T, N, M} <: ScatteredInterpolant
 
     w::Array{T,N}
@@ -143,6 +134,33 @@ function interpolate(rbf::RadialBasisFunction,
 
 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)
+    for (j, rbf) in enumerate(rbfs)
+        for i = 1:n
+            A[i,j] = rbf(A[i,j])
+        end
+    end
+
+    # Solve for the weights
+    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
+
+end
+
 function evaluate(itp::RBFInterpolant, points::AbstractArray{<:Real, 2})
 
     # Compute distance matrix
@@ -151,4 +169,16 @@ function evaluate(itp::RBFInterpolant, points::AbstractArray{<:Real, 2})
 
     # Compute the interpolated values
     return A*itp.w
-end
+end
+
+function evaluate(itp::RBFInterpolant{S,T,U,V}, points::AbstractArray{<:Real, 2}) where {S <: Vector, T, U, V}
+
+    # Compute distance matrix
+    A = pairwise(itp.metric, points, itp.points)
+    for (j, rbf) in enumerate(itp.rbf)
+        @. A[:,j] = rbf(A[:,j])
+    end
+
+    # Compute the interpolated values
+    return A*itp.w
+end

test/idw.jl

@@ -6,7 +6,7 @@ data = [0.0; 0.5; 0.5; 1.0]
 @testset "Shepard" begin
 
     @testset "Constructors" for idw in (:Shepard, )
-        @eval @test $idw(2) == $idw{2}() == $idw()
+        @eval @test $idw(2) == $idw()
     end
 
     @testset "Evaluation" begin

test/rbf.jl

@@ -3,12 +3,12 @@
 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())
+radialBasisFunctions = (Gaussian(2), Multiquadratic(2), InverseQuadratic(2), InverseMultiquadratic(2), Polyharmonic(2), ThinPlate())
 
 @testset "RBF" begin
 
     @testset "Constructors" for rbf in (:Gaussian, :Multiquadratic, :InverseQuadratic, :InverseMultiquadratic, :Polyharmonic)
-        @eval @test $rbf(1) == $rbf{1}() == $rbf()
+        @eval @test $rbf(1) == $rbf()
     end
 
     @testset "Evaluation" for r in radialBasisFunctions
@@ -34,6 +34,26 @@ radialBasisFunctions = (Gaussian{2}(), Multiquadratic{2}(), InverseQuadratic{2}(
         @test ev ≈ data
     end
 
+    @testset "Mixed RBF Evaluation" begin
+
+        RBFs = [Gaussian(2), Multiquadratic(2), InverseQuadratic(2), InverseMultiquadratic(2)]
+        itp = interpolate(RBFs, points, data)
+
+        # Check that we get back the original data at the sample points
+        ev = evaluate(itp, points)
+        @test ev ≈ data
+
+        @testset "Mixed RBF method equality" begin
+
+            itpConstant = interpolate(Gaussian(), points, data)
+            itpMixed = interpolate(repmat([Gaussian()],size(points,2)), points, data)
+
+            # Check that the result is the same when dispatching on multiple,
+            # but equal RBFs for each interpolation point
+            @test evaluate(itpConstant, points) == evaluate(itpMixed, points)
+        end
+    end
+
     @testset "returnRBFmatrix" begin
         r = radialBasisFunctions[1]
         @test typeof(interpolate(r, points, data; returnRBFmatrix = true)) <: Tuple