From 4188da5c0f071e8c48a09e0e2effadc4ebfa4671 Mon Sep 17 00:00:00 2001
From: Emil Ljungskog <emil.ljungskog@gmail.com>
Date: Wed, 14 Aug 2019 08:40:31 +0200
Subject: [PATCH] Restructure evaluateRBF! and exploit RBF matrix symmetry.

---
 src/rbf.jl  | 47 +++++++++++++++++------------------------------
 test/rbf.jl |  2 +-
 2 files changed, 18 insertions(+), 31 deletions(-)

diff --git a/src/rbf.jl b/src/rbf.jl
index ff3f674..30a1dce 100755
--- a/src/rbf.jl
+++ b/src/rbf.jl
@@ -195,7 +195,7 @@ function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialB
     # and optional smoothing (ridge regression)
     A = pairwise(metric, points;dims=2)
     
-    evaluateRBF!(A, rbf, smooth)
+    A = evaluateRBF!(A, rbf, smooth)
 
     # Solve for the weights
     itp = solveForWeights(A, points, samples, rbf, metric)
@@ -211,47 +211,34 @@ end
 
 @inline function evaluateRBF!(A, rbf)
     A .= rbf.(A)
+    return Symmetric(A)
 end
 @inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction})
     for (j, rbf) in enumerate(rbfs)
         A[:,j] .= rbf.(@view A[:,j])
     end
+    return A
 end
-@inline function evaluateRBF!(A, rbf, smooth::Bool)
-    A .= rbf.(A)
+@inline function evaluateRBF!(A, rbf, smooth)
+    A = evaluateRBF!(A, rbf)
+    A = addSmoothing!(A, smooth)
+    return A
 end
-@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction}, smooth::Bool)
-    for (j, rbf) in enumerate(rbfs)
-        A[:,j] .= rbf.(@view A[:,j])
-    end
+
+@inline function addSmoothing!(A, smooth::Bool)
+    return A
 end
-@inline function evaluateRBF!(A, rbf, smooth::Vector{T}) where T <: Number
-    A .= rbf.(A)
-    for i = 1:size(A,1)
+@inline function addSmoothing!(A, smooth::Vector{T}) where T <: Number
+    for i = 1:size(A, 1)
         A[i,i] += smooth[i]
     end
+    return A
 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)
+@inline function addSmoothing!(A, smooth::T) where T <: Number
+    for i = 1:size(A, 1)
         A[i,i] += smooth
     end
+    return A
 end
 
 @inline function solveForWeights(A, points, samples,
@@ -282,7 +269,7 @@ function evaluate(itp::RadialBasisInterpolant, points::AbstractArray{<:Real, 2})
 
     # Compute distance matrix and evaluate the RBF
     A = pairwise(itp.metric, points, itp.points;dims=2)
-    evaluateRBF!(A, itp.rbf)
+    A = evaluateRBF!(A, itp.rbf)
 
     # Compute polynomial matrix for generalized RBF:s
     P = getPolynomial(itp.rbf, points)
diff --git a/test/rbf.jl b/test/rbf.jl
index 00513bf..0814748 100755
--- a/test/rbf.jl
+++ b/test/rbf.jl
@@ -67,7 +67,7 @@ radialBasisFunctions = (Gaussian(2),
 
             # 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)
+            @test evaluate(itpConstant, points) ≈ evaluate(itpMixed, points)
         end
     end