Merge pull request #376 from archermarx/radial_input_checking

Input dimension checking when evaluating surrogates

lib/SurrogatesAbstractGPs/src/SurrogatesAbstractGPs.jl

@@ -1,6 +1,6 @@
 module SurrogatesAbstractGPs
 
-import Surrogates: add_point!, AbstractSurrogate, std_error_at_point
+import Surrogates: add_point!, AbstractSurrogate, std_error_at_point, _check_dimension
 export AbstractGPSurrogate, var_at_point, logpdf_surrogate
 
 using AbstractGPs
@@ -18,14 +18,17 @@ function AbstractGPSurrogate(x, y; gp = GP(Matern52Kernel()), Σy = 0.1)
     AbstractGPSurrogate(x, y, gp, posterior(gp(x, Σy), y), Σy)
 end
 
-# predictor 
+# predictor
 function (g::AbstractGPSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(g, val)
+
     return only(mean(g.gp_posterior([val])))
 end
 
 # for add point
-# copies of x and y need to be made because we get 
-#"Error: cannot resize array with shared data " if we push! directly to x and y  
+# copies of x and y need to be made because we get
+#"Error: cannot resize array with shared data " if we push! directly to x and y
 function add_point!(g::AbstractGPSurrogate, new_x, new_y)
     if new_x in g.x
         println("Adding a sample that already exists, cannot build AbstracgGPSurrogate.")

lib/SurrogatesFlux/src/SurrogatesFlux.jl

@@ -1,6 +1,6 @@
 module SurrogatesFlux
 
-import Surrogates: add_point!, AbstractSurrogate
+import Surrogates: add_point!, AbstractSurrogate, _check_dimension
 export NeuralSurrogate
 
 using Flux
@@ -37,6 +37,8 @@ function NeuralSurrogate(x, y, lb, ub; model = Chain(Dense(length(x[1]), 1), fir
 end
 
 function (my_neural::NeuralSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(my_neural, val)
     v = [val...]
     out = my_neural.model(v)
     if length(out) == 1

lib/SurrogatesPolyChaos/src/SurrogatesPolyChaos.jl

@@ -1,6 +1,6 @@
 module SurrogatesPolyChaos
 
-import Surrogates: AbstractSurrogate, add_point!
+import Surrogates: AbstractSurrogate, add_point!, _check_dimension
 export PolynomialChaosSurrogate
 
 using PolyChaos
@@ -37,6 +37,9 @@ function PolynomialChaosSurrogate(x, y, lb::Number, ub::Number;
 end
 
 function (pc::PolynomialChaosSurrogate)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(pc, val)
+
     return sum([pc.coeff[i] * PolyChaos.evaluate(val, pc.ortopolys)[i]
                 for i in 1:(pc.num_of_multi_indexes)])
 end
@@ -71,6 +74,9 @@ function PolynomialChaosSurrogate(x, y, lb, ub;
 end
 
 function (pcND::PolynomialChaosSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(pcND, val)
+
     sum = zero(eltype(val[1]))
     for i in 1:(pcND.num_of_multi_indexes)
         sum = sum +

lib/SurrogatesRandomForest/src/SurrogatesRandomForest.jl

@@ -1,6 +1,6 @@
 module SurrogatesRandomForest
 
-import Surrogates: add_point!, AbstractSurrogate
+import Surrogates: add_point!, AbstractSurrogate, _check_dimension
 export RandomForestSurrogate
 
 using XGBoost
@@ -19,6 +19,8 @@ function RandomForestSurrogate(x, y, lb::Number, ub::Number; num_round::Int = 1)
 end
 
 function (rndfor::RandomForestSurrogate)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(rndfor, val)
     return XGBoost.predict(rndfor.bst, reshape([val], 1, 1))[1]
 end
 
@@ -38,6 +40,8 @@ function RandomForestSurrogate(x, y, lb, ub; num_round::Int = 1)
 end
 
 function (rndfor::RandomForestSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(rndfor, val)
     return XGBoost.predict(rndfor.bst, reshape(collect(val), 1, length(val)))[1]
 end
 

src/Earth.jl

@@ -161,6 +161,8 @@ function EarthSurrogate(x, y, lb::Number, ub::Number; penalty::Number = 2.0,
 end
 
 function (earth::EarthSurrogate)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(earth, val)
     return sum([earth.coeff[i] * earth.basis[i](val) for i in 1:length(earth.coeff)]) +
            earth.intercept
 end
@@ -328,6 +330,8 @@ function EarthSurrogate(x, y, lb, ub; penalty::Number = 2.0, n_min_terms::Int =
 end
 
 function (earth::EarthSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(earth, val)
     return sum([earth.coeff[i] * prod([earth.basis[i][j](val[j]) for j in 1:length(val)])
                 for i in 1:length(earth.coeff)]) + earth.intercept
 end

src/GEK.jl

@@ -54,6 +54,9 @@ function _calc_gek_coeffs(x, y, p::Number, theta::Number)
 end
 
 function std_error_at_point(k::GEK, val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     phi(z) = exp(-(abs(z))^k.p)
     nd1 = length(k.y)
     n = length(k.x)
@@ -75,6 +78,9 @@ function std_error_at_point(k::GEK, val::Number)
 end
 
 function (k::GEK)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     phi = z -> exp(-(abs(z))^k.p)
     n = length(k.x)
     prediction = zero(eltype(k.x[1]))
@@ -159,6 +165,10 @@ function _calc_gek_coeffs(x, y, p, theta)
 end
 
 function std_error_at_point(k::GEK, val)
+
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     nd1 = length(k.y)
     n = length(k.x)
     d = length(k.x[1])
@@ -184,6 +194,9 @@ function std_error_at_point(k::GEK, val)
 end
 
 function (k::GEK)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     d = length(val)
     n = length(k.x)
     return k.mu +

src/GEKPLS.jl

@@ -61,7 +61,7 @@ function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, theta)
     println("struct created")
 end
 
-# predictor 
+# predictor
 function (g::GEKPLS)(X_test)
     n_eval, n_features_x = size(X_test)
     X_cont = (X_test .- g.X_offset) ./ g.X_scale
@@ -124,7 +124,7 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points)
         X: Concatenation of XX: [extra_points*nt, dim] - Extra points added (when extra_points > 0) and X
         y: Concatenation of yy[extra_points*nt, 1]- Extra points added (when extra_points > 0) and y
         """
-    # this function is equivalent to a combination of 
+    # this function is equivalent to a combination of
     # https://github.com/SMTorg/smt/blob/f124c01ffa78c04b80221dded278a20123dac742/smt/utils/kriging_utils.py#L1036
     # and https://github.com/SMTorg/smt/blob/f124c01ffa78c04b80221dded278a20123dac742/smt/surrogate_models/gekpls.py#L48
 
@@ -180,7 +180,7 @@ end
 
 ######start of bbdesign######
 
-# 
+#
 # Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia'
 # https://github.com/phrb/ExperimentalDesign.jl
 
@@ -206,7 +206,7 @@ end
 # COPYRIGHT HOLDERS BE  LIABLE FOR ANY CLAIM, DAMAGES OR  OTHER LIABILITY, WHETHER
 # IN  AN ACTION  OF  CONTRACT, TORT  OR  OTHERWISE,  ARISING FROM,  OUT  OF OR  IN
 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-# 
+#
 
 function boxbehnken(matrix_size::Int)
     boxbehnken(matrix_size, matrix_size)
@@ -259,7 +259,7 @@ end
 function standardization(X, y)
     """
     We substract the mean from each variable. Then, we divide the values of each
-    variable by its standard deviation. 
+    variable by its standard deviation.
 
     Parameters
     ----------
@@ -410,7 +410,7 @@ end
 
 function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise = 0.0)
     """
-    This function is a loose translation of SMT code from 
+    This function is a loose translation of SMT code from
     https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247
     It  determines the BLUP parameters and evaluates the reduced likelihood function for the given theta.
 
@@ -428,7 +428,7 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no
     Returns
     -------
     reduced_likelihood_function_value: real
-        - The value of the reduced likelihood function associated with the given autocorrelation parameters theta. 
+        - The value of the reduced likelihood function associated with the given autocorrelation parameters theta.
     beta:  Generalized least-squares regression weights
     gamma: Gaussian Process weights.
 
@@ -464,7 +464,7 @@ end
 
 ### MODIFIED PLS BELOW ###
 
-# The code below is a simplified version of 
+# The code below is a simplified version of
 # SKLearn's PLS
 # https://github.com/scikit-learn/scikit-learn/blob/80598905e/sklearn/cross_decomposition/_pls.py
 

src/InverseDistanceSurrogate.jl

@@ -20,6 +20,9 @@ function InverseDistanceSurrogate(x, y, lb, ub; p::Number = 1.0)
 end
 
 function (inverSurr::InverseDistanceSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(inverSurr, val)
+
     if val in inverSurr.x
         return inverSurr.y[findfirst(x -> x == val, inverSurr.x)]
     else

src/Kriging.jl

@@ -21,8 +21,13 @@ end
 Gives the current estimate for array 'val' with respect to the Kriging object k.
 """
 function (k::Kriging)(val)
-    d = length(val)
+
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     n = length(k.x)
+    d = length(val)
+
     return k.mu +
            sum(k.b[i] *
                exp(-sum(k.theta[j] * norm(val[j] - k.x[i][j])^k.p[j] for j in 1:d))
@@ -33,6 +38,10 @@ end
     Returns sqrt of expected mean_squared_error error at the point.
 """
 function std_error_at_point(k::Kriging, val)
+
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
+
     n = length(k.x)
     d = length(k.x[1])
     r = zeros(eltype(k.x[1]), n, 1)
@@ -58,6 +67,8 @@ end
 Gives the current estimate for 'val' with respect to the Kriging object k.
 """
 function (k::Kriging)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
     n = length(k.x)
     return k.mu + sum(k.b[i] * exp(-sum(k.theta * abs(val - k.x[i])^k.p)) for i in 1:n)
 end
@@ -66,6 +77,8 @@ end
     Returns sqrt of expected mean_squared_error error at the point.
 """
 function std_error_at_point(k::Kriging, val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(k, val)
     n = length(k.x)
     r = [exp(-k.theta * abs(val - k.x[i])^k.p) for i in 1:n]
     one = ones(eltype(k.x), n)
@@ -138,7 +151,7 @@ function _calc_kriging_coeffs(x, y, p::Number, theta::Number)
 
     mu = (one_t * inverse_of_R * y) / (one_t * inverse_of_R * one)
     b = inverse_of_R * (y - one * mu)
-    sigma = ((y - one * mu)' * inverse_of_R * (y - one * mu)) / n
+    sigma = ((y - one * mu)' * b) / n
     mu[1], b, sigma[1], inverse_of_R
 end
 
@@ -148,13 +161,13 @@ end
 Constructor for Kriging surrogate.
 
 - (x,y): sampled points
-- p: array of values 0<=p<=2 modeling the
+- p: array of values 0<=p<2 modeling the
      smoothness of the function being approximated in the i-th variable.
      low p -> rough, high p -> smooth
 - theta: array of values > 0 modeling how much the function is
           changing in the i-th variable.
 """
-function Kriging(x, y, lb, ub; p = 2 .* collect(one.(x[1])),
+function Kriging(x, y, lb, ub; p = 2.0 .* collect(one.(x[1])),
                  theta = [0.5 / max(1e-6 * norm(ub .- lb), std(x_i[i] for x_i in x))^p[i]
                           for i in 1:length(x[1])])
     if length(x) != length(unique(x))
@@ -219,7 +232,7 @@ function _calc_kriging_coeffs(x, y, p, theta)
 
     b = inverse_of_R * y_minus_1μ
 
-    sigma = (y_minus_1μ' * inverse_of_R * y_minus_1μ) / n
+    sigma = (y_minus_1μ' * b) / n
 
     mu[1], b, sigma[1], inverse_of_R
 end

src/LinearSurrogate.jl

@@ -51,6 +51,8 @@ function add_point!(my_linear::LinearSurrogate, new_x, new_y)
 end
 
 function (lin::LinearSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(lin, val)
     return lin.coeff' * [val...]
 end
 

src/Lobachevsky.jl

@@ -57,6 +57,9 @@ function LobachevskySurrogate(x, y, lb::Number, ub::Number; alpha::Number = 1.0,
 end
 
 function (loba::LobachevskySurrogate)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(loba, val)
+
     return sum(loba.coeff[j] * phi_nj1D(val, loba.x[j], loba.alpha, loba.n)
                for j in 1:length(loba.x))
 end
@@ -95,6 +98,8 @@ function LobachevskySurrogate(x, y, lb, ub; alpha = collect(one.(x[1])), n::Int
 end
 
 function (loba::LobachevskySurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(loba, val)
     return sum(loba.coeff[j] * phi_njND(val, loba.x[j], loba.alpha, loba.n)
                for j in 1:length(loba.x))
 end

src/PolynomialChaos.jl

@@ -32,6 +32,8 @@ function PolynomialChaosSurrogate(x, y, lb::Number, ub::Number;
 end
 
 function (pc::PolynomialChaosSurrogate)(val::Number)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(pc, val)
     return sum([pc.coeff[i] * PolyChaos.evaluate(val, pc.ortopolys)[i]
                 for i in 1:(pc.num_of_multi_indexes)])
 end
@@ -66,6 +68,8 @@ function PolynomialChaosSurrogate(x, y, lb, ub;
 end
 
 function (pcND::PolynomialChaosSurrogate)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(pcND, val)
     sum = zero(eltype(val[1]))
     for i in 1:(pcND.num_of_multi_indexes)
         sum = sum +

src/Radials.jl

@@ -63,7 +63,7 @@ function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse)
     if (typeof(y) == Vector{Float64}) #single output case
         coeff = copy(transpose(D \ y))
     else
-        coeff = copy(transpose(D \ Y[1:size(D)[1], :])) #if y is multi output; 
+        coeff = copy(transpose(D \ Y[1:size(D)[1], :])) #if y is multi output;
     end
     return coeff
 end
@@ -155,6 +155,9 @@ end
 Calculates current estimate of value 'val' with respect to the RadialBasis object.
 """
 function (rad::RadialBasis)(val)
+    # Check to make sure dimensions of input matches expected dimension of surrogate
+    _check_dimension(rad, val)
+
     approx = _approx_rbf(val, rad)
     return _match_container(approx, first(rad.y))
 end
@@ -171,9 +174,11 @@ function _approx_rbf(val::Number, rad::R) where {R}
     end
     return approx
 end
+
 function _approx_rbf(val, rad::R) where {R}
     n = length(rad.x)
     d = length(rad.x[1])
+
     q = rad.dim_poly
     num_poly_terms = binomial(q + d, q)
     lb = rad.lb