Merge pull request #146 from SciML/sparse_loba

Sparse lobachesky

src/Lobachesky.jl

@@ -1,11 +1,13 @@
-mutable struct LobacheskySurrogate{X,Y,A,N,L,U,C} <: AbstractSurrogate
+using ExtendableSparse
+mutable struct LobacheskySurrogate{X,Y,A,N,L,U,C,S} <: AbstractSurrogate
     x::X
     y::Y
     alpha::A
     n::N
     lb::L
     ub::U
     coeff::C
+    sparse::S
 end
 
 function phi_nj1D(point,x,alpha,n)
@@ -24,9 +26,13 @@ function phi_nj1D(point,x,alpha,n)
     return val
 end
 
-function _calc_loba_coeff1D(x,y,alpha,n)
+function _calc_loba_coeff1D(x,y,alpha,n,sparse)
     dim = length(x)
-    D = zeros(eltype(x[1]), dim, dim)
+    if sparse
+        D = ExtendableSparseMatrix{eltype(x),Int}(dim,dim)
+    else
+        D = zeros(eltype(x[1]), dim, dim)
+    end
     for i = 1:dim
         for j = 1:dim
             D[i,j] = phi_nj1D(x[i],x[j],alpha,n)
@@ -38,15 +44,15 @@ end
 """
 Lobachesky interpolation, suggested parameters: 0 <= alpha <= 4, n must be even.
 """
-function LobacheskySurrogate(x,y,lb::Number,ub::Number; alpha::Number = 1.0,n::Int=4)
+function LobacheskySurrogate(x,y,lb::Number,ub::Number; alpha::Number = 1.0,n::Int=4, sparse = false)
     if alpha > 4 || alpha < 0
         error("Alpha must be between 0 and 4")
     end
     if n % 2 != 0
         error("Parameter n must be even")
     end
-    coeff = _calc_loba_coeff1D(x,y,alpha,n)
-    LobacheskySurrogate(x,y,alpha,n,lb,ub,coeff)
+    coeff = _calc_loba_coeff1D(x,y,alpha,n,sparse)
+    LobacheskySurrogate(x,y,alpha,n,lb,ub,coeff,sparse)
 end
 
 function (loba::LobacheskySurrogate)(val::Number)
@@ -57,9 +63,13 @@ function phi_njND(point,x,alpha,n)
     return prod(phi_nj1D(point[h],x[h],alpha[h],n) for h = 1:length(x))
 end
 
-function _calc_loba_coeffND(x,y,alpha,n)
+function _calc_loba_coeffND(x,y,alpha,n,sparse)
     dim = length(x)
-    D = zeros(eltype(x[1]), dim, dim)
+    if sparse
+        D = ExtendableSparseMatrix{eltype(x[1]),Int}(dim,dim)
+    else
+        D = zeros(eltype(x[1]), dim, dim)
+    end
     for i = 1:dim
         for j = 1:dim
             D[i,j] = phi_njND(x[i],x[j],alpha,n)
@@ -69,16 +79,16 @@ function _calc_loba_coeffND(x,y,alpha,n)
     return Sym\y
 end
 """
-LobacheskySurrogate(x,y,alpha,n::Int,lb,ub)
+LobacheskySurrogate(x,y,alpha,n::Int,lb,ub,sparse = false)
 
 Build the Lobachesky surrogate with parameters alpha and n.
 """
-function LobacheskySurrogate(x,y,lb,ub; alpha = collect(one.(x[1])),n::Int = 4)
+function LobacheskySurrogate(x,y,lb,ub; alpha = collect(one.(x[1])),n::Int = 4, sparse = false)
     if n % 2 != 0
         error("Parameter n must be even")
     end
-    coeff = _calc_loba_coeffND(x,y,alpha,n)
-    LobacheskySurrogate(x,y,alpha,n,lb,ub,coeff)
+    coeff = _calc_loba_coeffND(x,y,alpha,n,sparse)
+    LobacheskySurrogate(x,y,alpha,n,lb,ub,coeff,sparse)
 end
 
 function (loba::LobacheskySurrogate)(val)
@@ -90,12 +100,12 @@ function add_point!(loba::LobacheskySurrogate,x_new,y_new)
         #1D
         append!(loba.x,x_new)
         append!(loba.y,y_new)
-        loba.coeff = _calc_loba_coeff1D(loba.x,loba.y,loba.alpha,loba.n)
+        loba.coeff = _calc_loba_coeff1D(loba.x,loba.y,loba.alpha,loba.n,loba.sparse)
     else
         #ND
         loba.x = vcat(loba.x,x_new)
         loba.y = vcat(loba.y,y_new)
-        loba.coeff = _calc_loba_coeffND(loba.x,loba.y,loba.alpha,loba.n)
+        loba.coeff = _calc_loba_coeffND(loba.x,loba.y,loba.alpha,loba.n,loba.sparse)
     end
     nothing
 end
@@ -180,5 +190,5 @@ function lobachesky_integrate_dimension(loba::LobacheskySurrogate,lb,ub,dim::Int
     new_lb = deleteat!(lb,dim)
     new_ub = deleteat!(ub,dim)
     new_loba = deleteat!(loba.alpha,dim)
-    return LobacheskySurrogate(new_x,loba.y,loba.alpha,loba.n,new_lb,new_ub,new_coeff)
+    return LobacheskySurrogate(new_x,loba.y,loba.alpha,loba.n,new_lb,new_ub,new_coeff,loba.sparse)
 end

test/lobachesky.jl

@@ -54,3 +54,24 @@ y = obj.(x)
 n = 4
 my_loba_ND = LobacheskySurrogate(x,y,lb,ub)
 lobachesky_integrate_dimension(my_loba_ND,lb,ub,2)
+
+#Sparse
+#1D
+obj = x -> 3*x + log(x)
+a = 1.0
+b = 4.0
+x = sample(100,a,b,SobolSample())
+y = obj.(x)
+alpha = 2.0
+n = 6
+my_loba = LobacheskySurrogate(x,y,a,b, alpha = 2.0, n = 6,sparse=true)
+
+#ND
+obj = x -> x[1] + log(x[2])
+lb = [0.0,0.0]
+ub = [8.0,8.0]
+alpha = [2.4,2.4]
+n = 8
+x = sample(100,lb,ub,SobolSample())
+y = obj.(x)
+my_loba_ND = LobacheskySurrogate(x,y,lb,ub,alpha = [2.4,2.4],n=8,sparse=true)