Added sparse radials

Project.toml

@@ -5,12 +5,14 @@ version = "1.1.2"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Stheno = "8188c328-b5d6-583d-959b-9690869a5511"
 XGBoost = "009559a3-9522-5dbb-924b-0b6ed2b22bb9"
 

src/Radials.jl

@@ -1,6 +1,7 @@
 using LinearAlgebra
+using ExtendableSparse
 
-mutable struct RadialBasis{F,Q,X,Y,L,U,C,S} <: AbstractSurrogate
+mutable struct RadialBasis{F,Q,X,Y,L,U,C,S,D} <: AbstractSurrogate
     phi::F
     dim_poly::Q
     x::X
@@ -9,6 +10,7 @@ mutable struct RadialBasis{F,Q,X,Y,L,U,C,S} <: AbstractSurrogate
     ub::U
     coeff::C
     scale_factor::S
+    sparse::D
 end
 
 mutable struct RadialFunction{Q,P}
@@ -26,60 +28,66 @@ thinplateRadial = RadialFunction(2,z->norm(z)^2*log(norm(z)))
 
 Constructor for Radial basis surrogate.
 """
-function RadialBasis(x, y, lb::Number, ub::Number; rad::RadialFunction=linearRadial, scale_factor::Real=1.0)
+function RadialBasis(x, y, lb::Number, ub::Number; rad::RadialFunction=linearRadial, scale_factor::Real=1.0, sparse = false)
     q = rad.q
     phi = rad.phi
-    coeff = _calc_coeffs(x, y, lb, ub, phi, q,scale_factor)
-    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor)
+    coeff = _calc_coeffs(x, y, lb, ub, phi, q,scale_factor, sparse)
+    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor,sparse)
 end
 
 """
 RadialBasis(x,y,lb,ub,rad::RadialFunction, scale_factor::Float = 1.0)
 
 Constructor for RadialBasis surrogate
 """
-function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0)
+function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0, sparse = false)
     q = rad.q
     phi = rad.phi
-    coeff = _calc_coeffs(x, y, lb, ub, phi, q, scale_factor)
-    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor)
+    coeff = _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse)
+    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor, sparse)
 end
 
-function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor)
-    D = _construct_rbf_interp_matrix(x, first(x), lb, ub, phi, q, scale_factor)
+function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse)
+    D = _construct_rbf_interp_matrix(x, first(x), lb, ub, phi, q, scale_factor, sparse)
     Y = _construct_rbf_y_matrix(y, first(y), length(y) + q)
     coeff = D \ Y
     return coeff
 end
 
-function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_factor)
+function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_factor, sparse)
     n = length(x)
     m = n + q
-    D = zeros(eltype(x_el), m, m)
-
+    if sparse
+        D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m)
+    else
+        D = zeros(eltype(x_el), m, m)
+    end
     @inbounds for i = 1:n
         for j = 1:n
             D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor )
         end
         if i <= n
             for k = 1:q
-                D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub)
+                    D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub)
             end
         end
     end
     D_sym = Symmetric(D, :U)
     return D_sym
 end
 
-function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor)
+function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor,sparse)
     n = length(x)
     nd = length(x_el)
     central_point = _center_bounds(x_el, lb, ub)
     sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:nd)
     mean_half_diameter = sum_half_diameter / nd
     m = n+q
-    D = zeros(eltype(x_el), m, m)
-
+    if sparse
+        D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m)
+    else
+        D = zeros(eltype(x_el), m, m)
+    end
     @inbounds for i = 1:n
         for j = 1:n
             D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor)
@@ -171,6 +179,6 @@ function add_point!(rad::RadialBasis,new_x,new_y)
         append!(rad.x,new_x)
         append!(rad.y,new_y)
     end
-    rad.coeff = _calc_coeffs(rad.x,rad.y,rad.lb,rad.ub,rad.phi,rad.dim_poly, rad.scale_factor)
+    rad.coeff = _calc_coeffs(rad.x,rad.y,rad.lb,rad.ub,rad.phi,rad.dim_poly, rad.scale_factor, rad.sparse)
     nothing
 end

test/Radials.jl

@@ -75,7 +75,7 @@ my_rad_ND = RadialBasis(x,y,lb,ub,rad = cubicRadial)
 my_rad_ND = RadialBasis(x,y,lb,ub, rad = multiquadricRadial)
 prediction = my_rad_ND((1.0,1.0,1.0))
 
-# #100
+
 f = x -> x[1]*x[2]
 lb = [1.0, 2.0]
 ub = [10.0, 8.5]
@@ -87,7 +87,6 @@ my_radial_basis = RadialBasis(x, y, lb, ub, rad = linearRadial)
 my_radial_basis = RadialBasis(x, y, lb, ub, rad =linearRadial)
 @test my_radial_basis((1.0, 2.0)) ≈ 2
 
-# #100
 f = x -> x[1]*x[2]
 lb = [1.0, 2.0]
 ub = [10.0, 8.5]
@@ -97,7 +96,7 @@ y = f.(x)
 my_radial_basis = RadialBasis(x, y, lb,ub, rad = linearRadial)
 @test my_radial_basis((1.0, 2.0)) ≈ 2
 
-# Multi-output #98
+# Multi-output
 f  = x -> [x^2, x]
 lb = 1.0
 ub = 10.0
@@ -122,3 +121,22 @@ x_new = (2.0, 2.0)
 y_new = f(x_new)
 add_point!(my_radial_basis, x_new, y_new)
 @test my_radial_basis(x_new) ≈ y_new
+
+
+
+#sparse
+
+#1D
+lb = 0.0
+ub = 4.0
+x = [1.0,2.0,3.0]
+y = [4.0,5.0,6.0]
+my_rad = RadialBasis(x, y, lb, ub, rad = linearRadial, sparse = true)
+
+#ND
+x = [(1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0)]
+y = [4.0, 5.0, 6.0]
+lb = [0.0,3.0,6.0]
+ub = [4.0,7.0,10.0]
+#bounds = [[0.0, 3.0, 6.0], [4.0, 7.0, 10.0]]
+my_rad = RadialBasis(x, y, lb, ub, sparse = true)