From 4d78d5351cfc387baf6104901a9e50a447a9a2a2 Mon Sep 17 00:00:00 2001 From: Vikram Date: Tue, 14 Jun 2022 19:07:01 +0530 Subject: [PATCH 01/12] replace sklearn pls with julia version of modified pls --- src/GEKPLS.jl | 467 +++++++++++++++++++++++++++++++++++++++++++++++++ test/GEKPLS.jl | 96 ++++++++++ 2 files changed, 563 insertions(+) create mode 100644 src/GEKPLS.jl create mode 100644 test/GEKPLS.jl diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl new file mode 100644 index 000000000..fc62a63bb --- /dev/null +++ b/src/GEKPLS.jl @@ -0,0 +1,467 @@ +using LinearAlgebra +using Statistics + +mutable struct GEKPLS <: AbstractSurrogate + x::Matrix{Float64} #1 + y::Matrix{Float64} #2 + grads::Matrix{Float64} #3 + xl::Matrix{Float64} #xlimits #4 + delta:: Float64 #5 + extra_points::Int64 #6 + num_components::Int64 #7 + beta::Vector{Float64} #8 + gamma::Matrix{Float64} #9 + theta::Vector{Float64} #10 + reduced_likelihood_function_value:: Float64 #11 + X_offset::Matrix{Float64} #12 + X_scale::Matrix{Float64} #13 + X_after_std::Matrix{Float64} #14 - X after standardization + pls_mean::Matrix{Float64} #15 + y_mean::Float64 #16 + y_std::Float64 #17 +end + +#constructor for GEKPLS Struct +function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) + theta = [θ for i in 1:n_comp] + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) + X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, y_after_PLS) + D, ij = cross_distances(X_after_std) + pls_mean_reshaped = reshape(pls_mean, (size(X,2),n_comp)) + d = componentwise_distance_PLS(D, "squar_exp", n_comp, pls_mean_reshaped) + nt, nd = size(X_after_PLS) + beta, gamma, reduced_likelihood_function_value = _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) + return GEKPLS(X, y, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, reduced_likelihood_function_value, + X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) + println("struct created") +end + + +# predictor +function (g::GEKPLS)(X_test) + n_eval, n_features_x = size(X_test) + X_cont = (X_test .- g.X_offset) ./ g.X_scale + dx = differences(X_cont, g.X_after_std) + pred_d = componentwise_distance_PLS(dx, "squar_exp", g.num_components, g.pls_mean) + nt = size(g.X_after_std,1) + r = transpose(reshape(squar_exp(g.theta, pred_d),(nt,n_eval))) + f = ones(n_eval,1) + y_ = (f * g.beta) + (r * g.gamma) + y = g.y_mean .+ g.y_std * y_ + return y +end + +function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) + """ + Gradient-enhanced PLS-coefficients. + Parameters + ---------- + X: [n_obs,dim] - The input variables. + y: [n_obs,ny] - The output variable + n_comp: int - Number of principal components used. + gradients: - The gradient values. Matrix size (n_obs,dim) + delta_x: real - The step used in the First Order Taylor Approximation + xlimits: [dim, 2]- The upper and lower var bounds. + extra_points: int - The number of extra points per each training point. + Returns + ------- + Coeff_pls: [dim, n_comp] - The PLS-coefficients. + 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 + # 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 + nt, dim = size(X) + XX = zeros(0,dim) + yy = zeros(0,size(y)[2]) + #_pls = PLSRegression(n_comp) #todo: relic from sklearn versiom. REMOVE. + coeff_pls = zeros((dim, n_comp, nt)) + for i in 1:nt + if dim >= 3 #to do ... this is consistent with SMT but inefficient. Move outside for loop and evaluate dim >= 3 just once. + bb_vals = circshift(boxbehnken(dim, 1),1) + _X = zeros((size(bb_vals)[1], dim)) + _y = zeros((size(bb_vals)[1], 1)) + bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals + _X = X[i, :]' .+ bb_vals + bb_vals = bb_vals .* grads[i,:]' + _y = y[i, :] .+ sum(bb_vals, dims=2) + else + println("GEKPLS for less than 3 dimensions is coming soon") + end + #_pls.fit(_X, _y) #todo: relic from sklearn versiom. REMOVE. + #coeff_pls[:, :, i] = _pls.x_rotations_ #size of _pls.x_rotations_ is (dim, n_comp) #todo: relic from sklearn versiom. REMOVE. + coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) + if extra_points != 0 + start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) #todo: evaluate just once + max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] + for ii in max_coeff + XX = [XX; transpose(X[i, :])] + XX[end, ii] += delta_x * (xlimits[ii,2]-xlimits[ii,1]) + yy = [yy; y[i]] + yy[end] += grads[i,ii] * delta_x * (xlimits[ii,2]-xlimits[ii,1]) + end + end + end + if extra_points != 0 + X = [X; XX] + y = [y; yy] + end + pls_mean = mean(broadcast(abs,coeff_pls),dims=3) + return pls_mean, X, y +end + +######start of bbdesign###### + + + # + # Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' + # https://github.com/phrb/ExperimentalDesign.jl + + # MIT License + + # ExperimentalDesign.jl: Design of Experiments in Julia + # Copyright (C) 2019 Pedro Bruel + + # Permission is hereby granted, free of charge, to any person obtaining a copy of + # this software and associated documentation files (the "Software"), to deal in + # the Software without restriction, including without limitation the rights to + # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + # the Software, and to permit persons to whom the Software is furnished to do so, + # subject to the following conditions: + + # The above copyright notice and this permission notice (including the next + # paragraph) shall be included in all copies or substantial portions of the + # Software. + + # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS + # FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR + # 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) + end + + function boxbehnken(matrix_size::Int, center::Int) + @assert matrix_size>=3 + + A_fact = explicit_fullfactorial(Tuple([-1,1] for i = 1:2)) + + rows = floor(Int, (0.5*matrix_size*(matrix_size-1))*size(A_fact)[1]) + + A = zeros(rows, matrix_size) + + l = 0 + for i in 1:matrix_size-1 + for j in i+1:matrix_size + l = l +1 + A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] + A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] + end + end + + if center == matrix_size + if matrix_size <= 16 + points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] + center = points[matrix_size] + end + end + + A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) + end + + function explicit_fullfactorial(factors::Tuple) + explicit_fullfactorial(fullfactorial(factors)) + end + + function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) + hcat(vcat.(collect(iterator)...)...) + end + + function fullfactorial(factors::Tuple) + Base.Iterators.product(factors...) + end + + +######end of bb design###### + +function standardization(X,y) + """ + We substract the mean from each variable. Then, we divide the values of each + variable by its standard deviation. + + Parameters + ---------- + + X - The input variables. + y - The output variable. + + Returns + ------- + + X: [n_obs, dim] + The standardized input matrix. + + y: [n_obs, 1] + The standardized output vector. + + X_offset: The mean (or the min if scale_X_to_unit=True) of each input variable. + + y_mean: The mean of the output variable. + + X_scale: The standard deviation of each input variable. + + y_std: The standard deviation of the output variable. + + """ + #Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L21 + X_offset = mean(X, dims = 1) + X_scale = std(X, dims = 1) + X_scale = map(x -> (x==0.0) ? x=1 : x=x, X_scale) #to prevent division by 0 below + y_mean = mean(y) + y_std = std(y) + y_std = map(y -> (y==0) ? y=1 : y=y, y_std) #to prevent division by 0 below + X = (X.-X_offset) ./ X_scale + y = (y .- y_mean) ./ y_std + return X, y, X_offset, y_mean, X_scale, y_std + +end + +function cross_distances(X) + """ + Computes the nonzero componentwise cross-distances between the vectors + in X + + Parameters + ---------- + + X: [n_obs, dim] + + Returns + ------- + D: [n_obs * (n_obs - 1) / 2, dim] + - The cross-distances between the vectors in X. + + ij: [n_obs * (n_obs - 1) / 2, 2] + - The indices i and j of the vectors in X associated to the cross- + distances in D. + """ + # equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L86 + n_samples, n_features = size(X) + n_nonzero_cross_dist = ( n_samples * (n_samples - 1) ) ÷ 2 + ij = zeros((n_nonzero_cross_dist, 2)) + D = zeros((n_nonzero_cross_dist, n_features)) + ll_1 = 0 + + for k in 1:n_samples - 1 + ll_0 = ll_1 + 1 + ll_1 = ll_0 + n_samples - k - 1 + ij[ll_0:ll_1, 1] .= k + ij[ll_0:ll_1, 2] = k+1:1:n_samples + D[ll_0:ll_1, :] = -(X[(k + 1) : n_samples,:] .- X[k,:]') + + end + return D, Int.(ij) +end + +function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) + """ + Computes the nonzero componentwise cross-spatial-correlation-distance + between the vectors in X. + + Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L1257 + with some simplifications (removed theta and return_derivative as it's not required for GEKPLS) + + Parameters + ---------- + + D: [n_obs * (n_obs - 1) / 2, dim] + - The L1 cross-distances between the vectors in X. + + corr: str + - Name of the correlation function used. + squar_exp or abs_exp. + + n_comp: int + - Number of principal components used. + + coeff_pls: [dim, n_comp] + - The PLS-coefficients. + + Returns + ------- + + D_corr: [n_obs * (n_obs - 1) / 2, n_comp] + - The componentwise cross-spatial-correlation-distance between the + vectors in X. + + """ + #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L1257 + #todo + #figure out how to handle this computation in the case of very large matrices + #similar to what SMT has done + #at https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L1257 + D_corr = zeros((size(D)[1], n_comp)) + + if corr == "squar_exp" + D_corr = D.^2 * coeff_pls.^2 + else #abs_exp + D_corr = abs.(D) * abs.(coeff_pls) + end + return D_corr +end + +function squar_exp(theta, d) + """ + Squared exponential correlation model. + Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L604 + Parameters: + ----------- + theta : Hyperparameters of the correlation model + d: componentwise distances from componentwise_distance_PLS + + Returns: + -------- + r: array containing the values of the autocorrelation model + + """ + n_components = size(d)[2] + theta = reshape(theta, (1,n_components)) + return exp.(-sum(theta .* d, dims=2)) +end + +function differences(X, Y) + #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L392 + #code credit: Elias Carvalho - https://stackoverflow.com/questions/72392010/row-wise-operations-between-matrices-in-julia + Rx = repeat(X, inner=(size(Y, 1), 1)) + Ry = repeat(Y, size(X, 1)) + return Rx - Ry +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 + 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. + + Parameters + ---------- + theta: array containing the parameters at which the Gaussian Process model parameters should be determined. + kernel_type: name of the correlation function. + d: The componentwise cross-spatial-correlation-distance between the vectors in X. + nt: number of training points + ij: The indices i and j of the vectors in X associated to the cross-distances in D. + y_norma: Standardized y values + noise: noise hyperparameter - increasing noise reduces reduced_likelihood_function_value + + + Returns + ------- + reduced_likelihood_function_value: real + - The value of the reduced likelihood function associated with the given autocorrelation parameters theta. + beta: Generalized least-squares regression weights + gamma: Gaussian Process weights. + + """ + #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 + reduced_likelihood_function_value = -Inf + nugget = 1000000.0 * eps() #a jitter for numerical stability; reducing the multiple from 1000000.0 results in positive definite error for Cholesky decomposition; + if kernel_type =="squar_exp" #todo - add other kernel types like abs_exp, matern52 etc. + r = squar_exp(theta,d) + end + R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) + for k in 1:size(ij)[1] #todo - speed up using @inbounds / @simd where required + R[ij[k,1],ij[k,2]]=r[k] + R[ij[k,2],ij[k,1]]=r[k] + end + C = cholesky(R).L #todo - #values diverge at this point from SMT code; verify impact + F = ones(nt,1) #todo - examine if this should be a parameter for this function + Ft = C\F + Q, G = qr(Ft) + Q = Array(Q) + Yt = C\y_norma + #todo - in smt, they check if the matrix is ill-conditioned using SVD. Verify and include if necessary + beta = G \ [(transpose(Q) ⋅ Yt)] + rho = Yt .- (Ft .* beta) + gamma = transpose(C) \ rho + sigma2 = sum((rho).^2, dims=1)/nt + detR = prod(diag(C).^(2.0/nt)) + reduced_likelihood_function_value = - nt * log10(sum(sigma2)) - nt * log10(detR) + return beta, gamma, reduced_likelihood_function_value +end + + +### MODIFIED PLS BELOW ### + +function _center_scale(X,Y) + x_mean = mean(X, dims = 1) + X .-= x_mean + y_mean = mean(Y, dims = 1) + Y .-= y_mean + + x_std = std(X, dims = 1) + x_std[x_std .== 0] .= 1.0 + X ./= x_std + y_std = std(Y, dims = 1) + y_std[y_std .==0] .= 1.0 + Y ./= y_std + return X, Y, x_mean, y_mean, x_std, y_std +end + +function _svd_flip_1d(u,v) + # equivalent of https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L149 + biggest_abs_val_idx = findmax(abs.(vec(u)))[2] + sign_ = sign(u[biggest_abs_val_idx]) + u .*= sign_ + v .*= sign_ + +end + +function _get_first_singular_vectors_power_method(X,Y) + #todo - do a check for Y residual being constant + #as in https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L66 + x_weights_old = 100 + y_score = vec(Y) + x_weights = transpose(X)y_score / dot(y_score, y_score) + my_eps = eps() + x_weights ./= (sqrt(dot(x_weights,x_weights)) + my_eps) + x_score = X * x_weights + y_weights = transpose(Y)x_score / dot(x_score, x_score) + y_score = Y * y_weights / (dot(y_weights, y_weights) + my_eps) + x_weights_diff = x_weights .- x_weights_old + return x_weights, y_weights +end + +function _modified_pls(X,Y, n_components) + x_weights_ = zeros(size(X,2),n_components) + y_weights_ = zeros(1, n_components) + _x_scores = zeros(size(X,1), n_components) + _y_scores = zeros(size(Y,1), n_components) + x_loadings_ = zeros(size(X,2),n_components) + y_loadings_ = zeros(1, n_components) + + Xk, Yk, x_mean, y_mean, x_std, y_std = _center_scale(X,Y) + + for k in 1:n_components + x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) + _svd_flip_1d(x_weights, y_weights) + x_scores = Xk * x_weights + y_ss = dot(y_weights, y_weights) + y_scores = (Yk * y_weights) ./ y_ss + x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) + Xk = Xk - (x_scores * x_loadings) + y_loadings = transpose(x_scores)*Yk / dot(x_scores, x_scores) + Yk = Yk - x_scores*y_loadings + x_weights_[:,k] = x_weights + _x_scores[:,k] = x_scores + x_loadings_[:,k] = vec(x_loadings) + end + x_rotations_ = x_weights_ * pinv(transpose(x_loadings_)x_weights_) + return x_rotations_ +end + +### MODIFIED PLS ABOVE ### diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl new file mode 100644 index 000000000..5e0a6ab88 --- /dev/null +++ b/test/GEKPLS.jl @@ -0,0 +1,96 @@ +using Surrogates +using Zygote + + +function water_flow(x) + r_w = x[1] + r = x[2] + T_u = x[3] + H_u = x[4] + T_l = x[5] + H_l = x[6] + L = x[7] + K_w = x[8] + log_val = log(r/r_w) + return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) +end + + +function vector_of_tuples_to_matrix(v) + #convert training data generated by surrogate sampling into a matrix suitable for GEKPLS + num_rows = length(v) + num_cols = length(first(v)) + K = zeros(num_rows, num_cols) + for row in 1:num_rows + for col in 1:num_cols + K[row, col]=v[row][col] + end + end + return K +end + +function vector_of_tuples_to_matrix2(v) + #convert gradients into matrix form + num_rows = length(v) + num_cols = length(first(first(v))) + K = zeros(num_rows, num_cols) + for row in 1:num_rows + for col in 1:num_cols + K[row, col] = v[row][1][col] + end + end + return K +end + +n = 1000 +d = 8 +lb = [0.05,100,63070,990,63.1,700,1120,9855] +ub = [0.15,50000,115600,1110,116,820,1680,12045] +x = sample(n,lb,ub,SobolSample()) +X = vector_of_tuples_to_matrix(x) +grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x)) +y = reshape(water_flow.(x),(size(x,1),1)) +xlimits = hcat(lb, ub) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) +y_true = water_flow.(x_test) + +@testset "Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2; initial_theta=.01)" begin + n_comp = 2 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println() + println("rmse: ", rmse) #rmse: 0.028999713540341848 + #@test isapprox(rmse, 0.02, atol=0.01) #score tolerance with sklearn's PLS + @test isapprox(rmse, 0.02, atol=0.02) + +end + +@testset "Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2; initial_theta=.01)" begin + n_comp = 3 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #rmse: 0.024098165334772537 + @test isapprox(rmse, 0.02, atol=0.01) +end + +@testset "Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3; initial_theta=.01)" begin + n_comp = 3 + delta_x = 0.0001 + extra_points = 3 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #0.023712136706924902 + @test isapprox(rmse, 0.02, atol=0.01) +end From d6b1a544e3d3452b75197bbd827db0792242208b Mon Sep 17 00:00:00 2001 From: Vikram Date: Wed, 15 Jun 2022 20:46:04 +0530 Subject: [PATCH 02/12] format after rebase from master --- Project.toml | 1 + src/GEKPLS.jl | 1 + src/Surrogates.jl | 3 ++- test/GEKPLS.jl | 59 ++++++++++++++++++++++++++++++++++++++++++----- test/runtests.jl | 12 ++++++++++ 5 files changed, 69 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 79a481846..d87569e97 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index fc62a63bb..9be303012 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -460,6 +460,7 @@ function _modified_pls(X,Y, n_components) _x_scores[:,k] = x_scores x_loadings_[:,k] = vec(x_loadings) end + x_rotations_ = x_weights_ * pinv(transpose(x_loadings_)x_weights_) return x_rotations_ end diff --git a/src/Surrogates.jl b/src/Surrogates.jl index ef2137fdf..8d418ab3f 100644 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -12,12 +12,12 @@ include("Lobachevsky.jl") include("LinearSurrogate.jl") include("InverseDistanceSurrogate.jl") include("SecondOrderPolynomialSurrogate.jl") - include("Wendland.jl") include("MOE.jl") #rewrite gaussian mixture with own algorithm to fix deps issue include("VariableFidelity.jl") include("Earth.jl") include("GEK.jl") +include("GEKPLS.jl") current_surrogates = ["Kriging", "LinearSurrogate", "LobachevskySurrogate", "NeuralSurrogate", @@ -81,6 +81,7 @@ function PolyChaosStructure(; op) end export current_surrogates +export GEKPLS export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 5e0a6ab88..6f16134f3 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -64,9 +64,7 @@ y_true = water_flow.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println() - println("rmse: ", rmse) #rmse: 0.028999713540341848 - #@test isapprox(rmse, 0.02, atol=0.01) #score tolerance with sklearn's PLS + println("rmse: ", rmse) # ~.039 with custom Julia PLS regressor; ~.029 with sklearn's PLS regressor @test isapprox(rmse, 0.02, atol=0.02) end @@ -79,7 +77,7 @@ end g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #rmse: 0.024098165334772537 + println("rmse: ", rmse) # ~.027 with our custom PLS regressor; ~0.024 with sklearn's PLS regressor @test isapprox(rmse, 0.02, atol=0.01) end @@ -88,9 +86,58 @@ end delta_x = 0.0001 extra_points = 3 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #0.023712136706924902 + println("rmse: ", rmse) # ~.027 with our custom PLS regressor; ~0.024 with sklearn's PLS regressor @test isapprox(rmse, 0.02, atol=0.01) end + + +function welded_beam(x) + h = x[1] + l = x[2] + t = x[3] + a = 6000/(sqrt(2)*h*l) + b = (6000*(14+0.5*l)*sqrt(0.25*(l^2+(h+t)^2)))/(2*(0.707*h*l*(l^2/12 + 0.25*(h+t)^2))) + return (sqrt(a^2+b^2 + l*a*b))/(sqrt(0.25*(l^2+(h+t)^2))) +end + +n = 1000 +d = 3 +lb = [0.125,5.0,5.0] +ub = [1.,10.,10.] +x = sample(n,lb,ub,SobolSample()) +X = vector_of_tuples_to_matrix(x) +grads = vector_of_tuples_to_matrix2(gradient.(welded_beam, x)) +y = reshape(welded_beam.(x),(size(x,1),1)) +xlimits = hcat(lb, ub) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) +y_true = welded_beam.(x_test) + +@testset "Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2; initial_theta=.01)" begin + n_comp = 2 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #rmse: 39.48 + @test isapprox(rmse, 39.5, atol=0.5) +end + +#increasing extra points increases accuracy +@testset "Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4; initial_theta=.01)" begin + n_comp = 2 + delta_x = 0.0001 + extra_points = 4 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #rmse: 37.87 + @test isapprox(rmse, 37.5, atol=0.5) +end diff --git a/test/runtests.jl b/test/runtests.jl index 9a5353848..8fd87c66c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", end end +<<<<<<< HEAD @time @safetestset "Radials.jl" begin include("Radials.jl") end @time @safetestset "Kriging.jl" begin include("Kriging.jl") end @time @safetestset "Sampling" begin include("sampling.jl") end @@ -23,6 +24,17 @@ end @time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end @time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end @time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end +======= +@time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end +@time @safetestset "Radials.jl" begin include("Radials.jl") end +@time @safetestset "Kriging.jl" begin include("Kriging.jl") end +@time @safetestset "Sampling" begin include("sampling.jl") end +@time @safetestset "Optimization" begin include("optimization.jl") end +@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end +@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end +@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end +@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end +>>>>>>> 838e708 (add more tests for GEKPLS) # @time @safetestset "AD_Compatibility" begin include("AD_compatibility.jl") end @time @safetestset "Wendland" begin include("Wendland.jl") end @time @safetestset "VariableFidelity" begin include("VariableFidelity.jl") end From 6e8e3ade49d1d3621c74f5e9564bf665a5dd1ace Mon Sep 17 00:00:00 2001 From: Vikram Date: Sat, 18 Jun 2022 20:08:39 +0530 Subject: [PATCH 03/12] add 2D box behnken FOTA --- src/GEKPLS.jl | 42 ++++++++++++++++++----- test/GEKPLS.jl | 91 +++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 108 insertions(+), 25 deletions(-) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 9be303012..1c75e7c2a 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -77,21 +77,41 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) yy = zeros(0,size(y)[2]) #_pls = PLSRegression(n_comp) #todo: relic from sklearn versiom. REMOVE. coeff_pls = zeros((dim, n_comp, nt)) + for i in 1:nt - if dim >= 3 #to do ... this is consistent with SMT but inefficient. Move outside for loop and evaluate dim >= 3 just once. + if dim >= 3 bb_vals = circshift(boxbehnken(dim, 1),1) - _X = zeros((size(bb_vals)[1], dim)) - _y = zeros((size(bb_vals)[1], 1)) - bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals - _X = X[i, :]' .+ bb_vals - bb_vals = bb_vals .* grads[i,:]' - _y = y[i, :] .+ sum(bb_vals, dims=2) + # _X = zeros((size(bb_vals)[1], dim)) + # _y = zeros((size(bb_vals)[1], 1)) + # bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals + # _X = X[i, :]' .+ bb_vals + # bb_vals = bb_vals .* grads[i,:]' + # _y = y[i, :] .+ sum(bb_vals, dims=2) else - println("GEKPLS for less than 3 dimensions is coming soon") + bb_vals = [ + 0.0 0.0; #center + 1.0 0.0; #right + 0.0 1.0; #up + -1.0 0.0; #left + 0.0 -1.0; #down + 1.0 1.0; #right up + -1.0 1.0; #left up + -1.0 -1.0; #left down + 1.0 -1.0; #right down + ] end + _X = zeros((size(bb_vals)[1], dim)) + _y = zeros((size(bb_vals)[1], 1)) + bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals + _X = X[i, :]' .+ bb_vals + bb_vals = bb_vals .* grads[i,:]' + _y = y[i, :] .+ sum(bb_vals, dims=2) + #_pls.fit(_X, _y) #todo: relic from sklearn versiom. REMOVE. #coeff_pls[:, :, i] = _pls.x_rotations_ #size of _pls.x_rotations_ is (dim, n_comp) #todo: relic from sklearn versiom. REMOVE. + coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) + if extra_points != 0 start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) #todo: evaluate just once max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] @@ -107,6 +127,7 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) X = [X; XX] y = [y; yy] end + pls_mean = mean(broadcast(abs,coeff_pls),dims=3) return pls_mean, X, y end @@ -312,6 +333,7 @@ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) else #abs_exp D_corr = abs.(D) * abs.(coeff_pls) end + return D_corr end @@ -370,14 +392,16 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 reduced_likelihood_function_value = -Inf nugget = 1000000.0 * eps() #a jitter for numerical stability; reducing the multiple from 1000000.0 results in positive definite error for Cholesky decomposition; - if kernel_type =="squar_exp" #todo - add other kernel types like abs_exp, matern52 etc. + if kernel_type =="squar_exp" #todo - add other kernel type abs_exp etc. r = squar_exp(theta,d) end R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) + for k in 1:size(ij)[1] #todo - speed up using @inbounds / @simd where required R[ij[k,1],ij[k,2]]=r[k] R[ij[k,2],ij[k,1]]=r[k] end + C = cholesky(R).L #todo - #values diverge at this point from SMT code; verify impact F = ones(nt,1) #todo - examine if this should be a parameter for this function Ft = C\F diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 6f16134f3..8252cc185 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -2,20 +2,6 @@ using Surrogates using Zygote -function water_flow(x) - r_w = x[1] - r = x[2] - T_u = x[3] - H_u = x[4] - T_l = x[5] - H_l = x[6] - L = x[7] - K_w = x[8] - log_val = log(r/r_w) - return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) -end - - function vector_of_tuples_to_matrix(v) #convert training data generated by surrogate sampling into a matrix suitable for GEKPLS num_rows = length(v) @@ -42,6 +28,20 @@ function vector_of_tuples_to_matrix2(v) return K end +# water flow function tests +function water_flow(x) + r_w = x[1] + r = x[2] + T_u = x[3] + H_u = x[4] + T_l = x[5] + H_l = x[6] + L = x[7] + K_w = x[8] + log_val = log(r/r_w) + return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) +end + n = 1000 d = 8 lb = [0.05,100,63070,990,63.1,700,1120,9855] @@ -66,7 +66,6 @@ y_true = water_flow.(x_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) println("rmse: ", rmse) # ~.039 with custom Julia PLS regressor; ~.029 with sklearn's PLS regressor @test isapprox(rmse, 0.02, atol=0.02) - end @testset "Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2; initial_theta=.01)" begin @@ -93,7 +92,7 @@ end @test isapprox(rmse, 0.02, atol=0.01) end - +# welded beam tests function welded_beam(x) h = x[1] l = x[2] @@ -141,3 +140,63 @@ end println("rmse: ", rmse) #rmse: 37.87 @test isapprox(rmse, 37.5, atol=0.5) end + +#sphere function tests +function sphere_function(x) + return sum(x.^2) +end + +#3D +n = 100 +d = 3 +lb = [-5.0, -5.0, -5.0] +ub = [5.0, 5.0 ,5.0] +x = sample(n,lb,ub,SobolSample()) +X = vector_of_tuples_to_matrix(x) +grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) +y = reshape(sphere_function.(x),(size(x,1),1)) +xlimits = hcat(lb, ub) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) +y_true = sphere_function.(x_test) + +@testset "Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2; initial_theta = .01" begin + n_comp = 2 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #rmse: 0.0008 + @test isapprox(rmse, 0.001, atol=0.05) +end + +#2D +n = 20 # for vals > ~30 and < ~7 "PosDefException: matrix is not Hermitian; Cholesky factorization failed." +d = 2 +lb = [-10.0, -10.0] +ub = [10.0, 10.0] +x = sample(n,lb,ub,SobolSample()) +#x[1] = (-1.0, 2.0) #LINE INSERTED BY VIK FOR TEMPORARY TESTING OF 2D BB OUTPUTS; DELETE [ALSO, setting to (-1.0, 1.0) throws error] +X = vector_of_tuples_to_matrix(x) +grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) +y = reshape(sphere_function.(x),(size(x,1),1)) +xlimits = hcat(lb, ub) +n_test = 10 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) +y_true = sphere_function.(x_test) + +@testset "Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2; initial_theta = .01" begin + n_comp = 2 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) #rmse: 0.01 + @test isapprox(rmse, 0.1, atol=0.5) +end From bb8c534e77f1f70a11b77acf9bfcf38ddc21f37a Mon Sep 17 00:00:00 2001 From: Vikram Date: Sun, 19 Jun 2022 19:46:40 +0530 Subject: [PATCH 04/12] provide fix to prevent NaN values from creeping into matrices --- src/GEKPLS.jl | 35 +++++++++++++++++++---------------- test/GEKPLS.jl | 49 ++++++++++++++++++++++++++----------------------- 2 files changed, 45 insertions(+), 39 deletions(-) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 1c75e7c2a..c5ede2083 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -106,7 +106,7 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) _X = X[i, :]' .+ bb_vals bb_vals = bb_vals .* grads[i,:]' _y = y[i, :] .+ sum(bb_vals, dims=2) - + #_pls.fit(_X, _y) #todo: relic from sklearn versiom. REMOVE. #coeff_pls[:, :, i] = _pls.x_rotations_ #size of _pls.x_rotations_ is (dim, n_comp) #todo: relic from sklearn versiom. REMOVE. @@ -420,20 +420,24 @@ end ### MODIFIED PLS BELOW ### +# Note for developers: +# 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 + function _center_scale(X,Y) x_mean = mean(X, dims = 1) X .-= x_mean y_mean = mean(Y, dims = 1) Y .-= y_mean - x_std = std(X, dims = 1) x_std[x_std .== 0] .= 1.0 X ./= x_std y_std = std(Y, dims = 1) y_std[y_std .==0] .= 1.0 Y ./= y_std - return X, Y, x_mean, y_mean, x_std, y_std + return X,Y end function _svd_flip_1d(u,v) @@ -446,36 +450,35 @@ function _svd_flip_1d(u,v) end function _get_first_singular_vectors_power_method(X,Y) - #todo - do a check for Y residual being constant - #as in https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L66 - x_weights_old = 100 + my_eps = eps() y_score = vec(Y) x_weights = transpose(X)y_score / dot(y_score, y_score) - my_eps = eps() x_weights ./= (sqrt(dot(x_weights,x_weights)) + my_eps) x_score = X * x_weights y_weights = transpose(Y)x_score / dot(x_score, x_score) y_score = Y * y_weights / (dot(y_weights, y_weights) + my_eps) - x_weights_diff = x_weights .- x_weights_old + #Equivalent in intent to https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L66 + if any(isnan.(x_weights)) || any(isnan.(y_weights)) + return false, false + end return x_weights, y_weights end function _modified_pls(X,Y, n_components) x_weights_ = zeros(size(X,2),n_components) - y_weights_ = zeros(1, n_components) _x_scores = zeros(size(X,1), n_components) - _y_scores = zeros(size(Y,1), n_components) x_loadings_ = zeros(size(X,2),n_components) - y_loadings_ = zeros(1, n_components) + Xk, Yk = _center_scale(X,Y) - Xk, Yk, x_mean, y_mean, x_std, y_std = _center_scale(X,Y) - for k in 1:n_components - x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) + x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) + + if x_weights == false + break + end + _svd_flip_1d(x_weights, y_weights) x_scores = Xk * x_weights - y_ss = dot(y_weights, y_weights) - y_scores = (Yk * y_weights) ./ y_ss x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) Xk = Xk - (x_scores * x_loadings) y_loadings = transpose(x_scores)*Yk / dot(x_scores, x_scores) diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 8252cc185..91b34c209 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -56,7 +56,7 @@ x_test = sample(n_test,lb,ub,GoldenSample()) X_test = vector_of_tuples_to_matrix(x_test) y_true = water_flow.(x_test) -@testset "Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2; initial_theta=.01)" begin +@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 @@ -64,11 +64,10 @@ y_true = water_flow.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) # ~.039 with custom Julia PLS regressor; ~.029 with sklearn's PLS regressor - @test isapprox(rmse, 0.02, atol=0.02) + @test isapprox(rmse, 0.03, atol=0.02) #rmse: 0.039 end -@testset "Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2; initial_theta=.01)" begin +@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 @@ -76,11 +75,10 @@ end g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) # ~.027 with our custom PLS regressor; ~0.024 with sklearn's PLS regressor - @test isapprox(rmse, 0.02, atol=0.01) + @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 end -@testset "Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3; initial_theta=.01)" begin +@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin n_comp = 3 delta_x = 0.0001 extra_points = 3 @@ -88,8 +86,7 @@ end g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) # ~.027 with our custom PLS regressor; ~0.024 with sklearn's PLS regressor - @test isapprox(rmse, 0.02, atol=0.01) + @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 end # welded beam tests @@ -116,7 +113,18 @@ x_test = sample(n_test,lb,ub,GoldenSample()) X_test = vector_of_tuples_to_matrix(x_test) y_true = welded_beam.(x_test) -@testset "Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2; initial_theta=.01)" begin +@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin + n_comp = 3 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + y_pred = g(X_test) + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 39.0, atol=0.5) #rmse: 38.988 +end + +@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 @@ -124,12 +132,11 @@ y_true = welded_beam.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #rmse: 39.48 - @test isapprox(rmse, 39.5, atol=0.5) + @test isapprox(rmse, 39.5, atol=0.5) #rmse: 39.481 end #increasing extra points increases accuracy -@testset "Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4; initial_theta=.01)" begin +@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin n_comp = 2 delta_x = 0.0001 extra_points = 4 @@ -137,8 +144,7 @@ end g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #rmse: 37.87 - @test isapprox(rmse, 37.5, atol=0.5) + @test isapprox(rmse, 37.5, atol=0.5) #rmse: 37.87 end #sphere function tests @@ -161,7 +167,7 @@ x_test = sample(n_test,lb,ub,GoldenSample()) X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) -@testset "Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2; initial_theta = .01" begin +@testset "Test 7: Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 @@ -169,17 +175,15 @@ y_true = sphere_function.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #rmse: 0.0008 - @test isapprox(rmse, 0.001, atol=0.05) + @test isapprox(rmse, 0.001, atol=0.05) #rmse: 0.00083 end #2D -n = 20 # for vals > ~30 and < ~7 "PosDefException: matrix is not Hermitian; Cholesky factorization failed." +n = 50 d = 2 lb = [-10.0, -10.0] ub = [10.0, 10.0] x = sample(n,lb,ub,SobolSample()) -#x[1] = (-1.0, 2.0) #LINE INSERTED BY VIK FOR TEMPORARY TESTING OF 2D BB OUTPUTS; DELETE [ALSO, setting to (-1.0, 1.0) throws error] X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) y = reshape(sphere_function.(x),(size(x,1),1)) @@ -189,7 +193,7 @@ x_test = sample(n_test,lb,ub,GoldenSample()) X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) -@testset "Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2; initial_theta = .01" begin +@testset "Test 8: Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 @@ -197,6 +201,5 @@ y_true = sphere_function.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) #rmse: 0.01 - @test isapprox(rmse, 0.1, atol=0.5) + @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 end From 64f527b3ff5533c907edde1680cd5a40122e069e Mon Sep 17 00:00:00 2001 From: Vikram Date: Mon, 20 Jun 2022 19:35:36 +0530 Subject: [PATCH 05/12] fix formatting after git rebase --- docs/pages.jl | 77 ++++++++++++++++++++++-------------------- docs/src/gekpls.md | 84 ++++++++++++++++++++++++++++++++++++++++++++++ src/GEKPLS.jl | 12 +++---- test/runtests.jl | 11 ------ 4 files changed, 130 insertions(+), 54 deletions(-) create mode 100644 docs/src/gekpls.md diff --git a/docs/pages.jl b/docs/pages.jl index d4ce6dd70..5101a1e11 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,37 +1,40 @@ -pages = ["index.md" - "Tutorials" => [ - "Basics" => "tutorials.md", - "Radials" => "radials.md", - "Kriging" => "kriging.md", - "Gaussian Process" => "abstractgps.md", - "Lobachevsky" => "lobachevsky.md", - "Linear" => "LinearSurrogate.md", - "InverseDistance" => "InverseDistance.md", - "RandomForest" => "randomforest.md", - "SecondOrderPolynomial" => "secondorderpoly.md", - "NeuralSurrogate" => "neural.md", - "Wendland" => "wendland.md", - "Polynomial Chaos" => "polychaos.md", - "Variable Fidelity" => "variablefidelity.md", - "Gradient Enhanced Kriging" => "gek.md", - ] - "User guide" => [ - "Samples" => "samples.md", - "Surrogates" => "surrogate.md", - "Optimization" => "optimizations.md", - ] - "Benchmarks" => [ - "Sphere function" => "sphere_function.md", - "Lp norm" => "lp.md", - "Rosenbrock" => "rosenbrock.md", - "Tensor product" => "tensor_prod.md", - "Cantilever beam" => "cantilever.md", - "Water Flow function" => "water_flow.md", - "Welded beam function" => "welded_beam.md", - "Branin function" => "BraninFunction.md", - "Ackley function" => "ackley.md", - "Gramacy & Lee Function" => "gramacylee.md", - "Salustowicz Benchmark" => "Salustowicz.md", - "Multi objective optimization" => "multi_objective_opt.md", - ] - "Contributing" => "contributing.md"] +pages = [ + "index.md" + "Tutorials" => [ + "Basics" => "tutorials.md", + "Radials" => "radials.md", + "Kriging" => "kriging.md", + "Gaussian Process" =>"abstractgps.md", + "Lobachevsky" => "lobachevsky.md", + "Linear" => "LinearSurrogate.md", + "InverseDistance" => "InverseDistance.md", + "RandomForest" => "randomforest.md", + "SecondOrderPolynomial" => "secondorderpoly.md", + "NeuralSurrogate" => "neural.md", + "Wendland" => "wendland.md", + "Polynomial Chaos" => "polychaos.md", + "Variable Fidelity" => "variablefidelity.md", + "Gradient Enhanced Kriging" => "gek.md", + "GEKPLS" => "gekpls.md" + ] + "User guide" => [ + "Samples" => "samples.md", + "Surrogates" => "surrogate.md", + "Optimization" => "optimizations.md" + ] + "Benchmarks" => [ + "Sphere function" => "sphere_function.md", + "Lp norm" => "lp.md", + "Rosenbrock" => "rosenbrock.md", + "Tensor product" => "tensor_prod.md", + "Cantilever beam" => "cantilever.md", + "Water Flow function" => "water_flow.md", + "Welded beam function" => "welded_beam.md", + "Branin function" => "BraninFunction.md", + "Ackley function" => "ackley.md", + "Gramacy & Lee Function" => "gramacylee.md", + "Salustowicz Benchmark" => "Salustowicz.md", + "Multi objective optimization" => "multi_objective_opt.md" + ] + "Contributing" => "contributing.md" + ] diff --git a/docs/src/gekpls.md b/docs/src/gekpls.md new file mode 100644 index 000000000..7a5095d98 --- /dev/null +++ b/docs/src/gekpls.md @@ -0,0 +1,84 @@ +## GEKPLS Surrogate Tutorial + +Gradient Enhanced Kriging with Partial Least Squares Method (GEKPLS) is a surrogate modelling technique that brings down computation time and returns improved accuracy for high-dimensional problems. The Julia implementation of GEKPLS is adapted from the Python version by [SMT](https://github.com/SMTorg) which is based on this [paper](https://arxiv.org/pdf/1708.02663.pdf). + +The following are the inputs when building a GEKPLS surrogate: + +1. X - The matrix containing the training points +2. y - The vector containing the training outputs associated with each of the training points +3. grads - The gradients at each of the input X training points +4. n_comp - Number of components to retain for the partial least squares regression (PLS) +5. delta_x - The step size to use for the first order Taylor approximation +6. xlimits - The lower and upper bounds for the training points +7. extra_points - The number of additional points to use for the PLS +8. theta - The hyperparameter to use for the correlation model + +The following example illustrates how to use GEKPLS: + +```@example gekpls_water_flow + +using Surrogates +using Zygote + +function vector_of_tuples_to_matrix(v) + #helper function to convert training data generated by surrogate sampling into a matrix suitable for GEKPLS + num_rows = length(v) + num_cols = length(first(v)) + K = zeros(num_rows, num_cols) + for row in 1:num_rows + for col in 1:num_cols + K[row, col]=v[row][col] + end + end + return K +end + +function vector_of_tuples_to_matrix2(v) + #helper function to convert gradients into matrix form + num_rows = length(v) + num_cols = length(first(first(v))) + K = zeros(num_rows, num_cols) + for row in 1:num_rows + for col in 1:num_cols + K[row, col] = v[row][1][col] + end + end + return K +end + +function water_flow(x) + r_w = x[1] + r = x[2] + T_u = x[3] + H_u = x[4] + T_l = x[5] + H_l = x[6] + L = x[7] + K_w = x[8] + log_val = log(r/r_w) + return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) +end + +n = 1000 +d = 8 +lb = [0.05,100,63070,990,63.1,700,1120,9855] +ub = [0.15,50000,115600,1110,116,820,1680,12045] +x = sample(n,lb,ub,SobolSample()) +X = vector_of_tuples_to_matrix(x) +grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x)) +y = reshape(water_flow.(x),(size(x,1),1)) +xlimits = hcat(lb, ub) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) +y_true = water_flow.(x_test) +n_comp = 2 +delta_x = 0.0001 +extra_points = 2 +initial_theta = 0.01 +g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) +y_pred = g(X_test) +rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) #root mean squared error +println(rmse) +``` + diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index c5ede2083..76be055a9 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -107,13 +107,13 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) bb_vals = bb_vals .* grads[i,:]' _y = y[i, :] .+ sum(bb_vals, dims=2) - #_pls.fit(_X, _y) #todo: relic from sklearn versiom. REMOVE. - #coeff_pls[:, :, i] = _pls.x_rotations_ #size of _pls.x_rotations_ is (dim, n_comp) #todo: relic from sklearn versiom. REMOVE. + #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. + #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. - coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) + coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ if extra_points != 0 - start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) #todo: evaluate just once + start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] for ii in max_coeff XX = [XX; transpose(X[i, :])] @@ -420,7 +420,7 @@ end ### MODIFIED PLS BELOW ### -# Note for developers: + # 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 @@ -476,7 +476,7 @@ function _modified_pls(X,Y, n_components) if x_weights == false break end - + _svd_flip_1d(x_weights, y_weights) x_scores = Xk * x_weights x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) diff --git a/test/runtests.jl b/test/runtests.jl index 8fd87c66c..5f32d5058 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,16 +15,6 @@ for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", end end -<<<<<<< HEAD -@time @safetestset "Radials.jl" begin include("Radials.jl") end -@time @safetestset "Kriging.jl" begin include("Kriging.jl") end -@time @safetestset "Sampling" begin include("sampling.jl") end -@time @safetestset "Optimization" begin include("optimization.jl") end -@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end -@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end -@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end -@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end -======= @time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end @time @safetestset "Radials.jl" begin include("Radials.jl") end @time @safetestset "Kriging.jl" begin include("Kriging.jl") end @@ -34,7 +24,6 @@ end @time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end @time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end @time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end ->>>>>>> 838e708 (add more tests for GEKPLS) # @time @safetestset "AD_Compatibility" begin include("AD_compatibility.jl") end @time @safetestset "Wendland" begin include("Wendland.jl") end @time @safetestset "VariableFidelity" begin include("VariableFidelity.jl") end From f31d42c984791ee92c46aaa285da2ce1f422a399 Mon Sep 17 00:00:00 2001 From: Vikram Date: Tue, 21 Jun 2022 15:55:02 +0530 Subject: [PATCH 06/12] add add_point function, tests for add_point function and bounds checking --- src/GEKPLS.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ test/GEKPLS.jl | 48 +++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 87 insertions(+), 7 deletions(-) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 76be055a9..587d33c24 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -21,8 +21,28 @@ mutable struct GEKPLS <: AbstractSurrogate y_std::Float64 #17 end +function bounds_error(x, xl) + num_x_rows = size(x,1) + num_dim = size(xl, 1) + for i in 1:num_x_rows + for j in 1:num_dim + if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2]) + return true + end + end + end + return false +end + #constructor for GEKPLS Struct function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) + + #ensure that X values are within the upper and lower bounds + if bounds_error(X, xlimits) + println("X values outside bounds") + return + end + theta = [θ for i in 1:n_comp] pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, y_after_PLS) @@ -51,6 +71,32 @@ function (g::GEKPLS)(X_test) return y end + +function add_point!(g::GEKPLS, new_x, new_y, new_grads) + if new_x in g.x + println("Adding a sample that already exists. Cannot build GEKPLS") + return + end + + if bounds_error(new_x, g.xl) + println("x values outside bounds") + return + end + + + g.x = vcat(g.x, new_x) + g.y = vcat(g.y, new_y) + g.grads = vcat(g.grads, new_grads) + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) + g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, y_after_PLS) + D, ij = cross_distances(g.X_after_std) + #pls_mean_reshaped = reshape(pls_mean, (size(g.x, 2), g.num_components)) + g.pls_mean = reshape(pls_mean, (size(g.x, 2), g.num_components)) + d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) + nt, nd = size(X_after_PLS) + g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) +end + function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) """ Gradient-enhanced PLS-coefficients. diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 91b34c209..594f3487e 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -28,7 +28,7 @@ function vector_of_tuples_to_matrix2(v) return K end -# water flow function tests +# # water flow function tests function water_flow(x) r_w = x[1] r = x[2] @@ -61,7 +61,7 @@ y_true = water_flow.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) @test isapprox(rmse, 0.03, atol=0.02) #rmse: 0.039 @@ -89,7 +89,7 @@ end @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 end -# welded beam tests +## welded beam tests function welded_beam(x) h = x[1] l = x[2] @@ -135,7 +135,7 @@ end @test isapprox(rmse, 39.5, atol=0.5) #rmse: 39.481 end -#increasing extra points increases accuracy +## increasing extra points increases accuracy @testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin n_comp = 2 delta_x = 0.0001 @@ -147,12 +147,12 @@ end @test isapprox(rmse, 37.5, atol=0.5) #rmse: 37.87 end -#sphere function tests +## sphere function tests function sphere_function(x) return sum(x.^2) end -#3D +## 3D n = 100 d = 3 lb = [-5.0, -5.0, -5.0] @@ -178,7 +178,7 @@ y_true = sphere_function.(x_test) @test isapprox(rmse, 0.001, atol=0.05) #rmse: 0.00083 end -#2D +## 2D n = 50 d = 2 lb = [-10.0, -10.0] @@ -201,5 +201,39 @@ y_true = sphere_function.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + println("rmse: ", rmse) @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 end + +@testset "Test 9: Add Point Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin + #first we create a surrogate model with just 3 input points + initial_x_vec = [(1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0)] + initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec,1),1)) + initial_X = vector_of_tuples_to_matrix(initial_x_vec) + initial_grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, initial_x_vec)) + lb = [-5.0, -5.0, -5.0] + ub = [10.0, 10.0, 10.0] + xlimits = hcat(lb, ub) + n_comp = 2 + delta_x = 0.0001 + extra_points = 2 + initial_theta = 0.01 + g = GEKPLS(initial_X, initial_y, initial_grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + n_test = 100 + x_test = sample(n_test,lb,ub,GoldenSample()) + X_test = vector_of_tuples_to_matrix(x_test) + y_true = sphere_function.(x_test) + y_pred1 = g(X_test) + rmse1 = sqrt(sum(((y_pred1 - y_true).^2)/n_test)) #rmse1 = 31.91 + + #then we update the model with more points to see if performance improves + n = 100 + x = sample(n,lb,ub,SobolSample()) + X = vector_of_tuples_to_matrix(x) + grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) + y = reshape(sphere_function.(x),(size(x,1),1)) + add_point!(g, X, y, grads) + y_pred2 = g(X_test) + rmse2 = sqrt(sum(((y_pred2 - y_true).^2)/n_test)) #rmse2 = 0.0015 + @test (rmse2 < rmse1) +end From 2e2cd8815fc467b0a75204baa6550457be59af5f Mon Sep 17 00:00:00 2001 From: Vikram Date: Sat, 25 Jun 2022 19:56:16 +0530 Subject: [PATCH 07/12] change struct of GEKPLS --- src/GEKPLS.jl | 49 ++++++++++++++++++++----------------------------- test/GEKPLS.jl | 1 - 2 files changed, 20 insertions(+), 30 deletions(-) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 587d33c24..1c8345e73 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -1,24 +1,24 @@ using LinearAlgebra using Statistics -mutable struct GEKPLS <: AbstractSurrogate - x::Matrix{Float64} #1 - y::Matrix{Float64} #2 - grads::Matrix{Float64} #3 - xl::Matrix{Float64} #xlimits #4 - delta:: Float64 #5 - extra_points::Int64 #6 - num_components::Int64 #7 - beta::Vector{Float64} #8 - gamma::Matrix{Float64} #9 - theta::Vector{Float64} #10 - reduced_likelihood_function_value:: Float64 #11 - X_offset::Matrix{Float64} #12 - X_scale::Matrix{Float64} #13 - X_after_std::Matrix{Float64} #14 - X after standardization - pls_mean::Matrix{Float64} #15 - y_mean::Float64 #16 - y_std::Float64 #17 +mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate + x::Matrix{T} #1 + y::Matrix{T} #2 + grads::Matrix{T} #3 + xl::Matrix{T} #xlimits #4 + delta:: T #5 + extra_points::Int #6 + num_components::Int #7 + beta::Vector{T} #8 + gamma::Matrix{T} #9 + theta::Vector{T} #10 + reduced_likelihood_function_value:: T #11 + X_offset::Matrix{T} #12 + X_scale::Matrix{T} #13 + X_after_std::Matrix{T} #14 - X after standardization + pls_mean::Matrix{T} #15 + y_mean::T #16 + y_std::T #17 end function bounds_error(x, xl) @@ -90,7 +90,6 @@ function add_point!(g::GEKPLS, new_x, new_y, new_grads) pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, y_after_PLS) D, ij = cross_distances(g.X_after_std) - #pls_mean_reshaped = reshape(pls_mean, (size(g.x, 2), g.num_components)) g.pls_mean = reshape(pls_mean, (size(g.x, 2), g.num_components)) d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) nt, nd = size(X_after_PLS) @@ -121,18 +120,11 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) nt, dim = size(X) XX = zeros(0,dim) yy = zeros(0,size(y)[2]) - #_pls = PLSRegression(n_comp) #todo: relic from sklearn versiom. REMOVE. coeff_pls = zeros((dim, n_comp, nt)) for i in 1:nt if dim >= 3 bb_vals = circshift(boxbehnken(dim, 1),1) - # _X = zeros((size(bb_vals)[1], dim)) - # _y = zeros((size(bb_vals)[1], 1)) - # bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals - # _X = X[i, :]' .+ bb_vals - # bb_vals = bb_vals .* grads[i,:]' - # _y = y[i, :] .+ sum(bb_vals, dims=2) else bb_vals = [ 0.0 0.0; #center @@ -152,12 +144,11 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) _X = X[i, :]' .+ bb_vals bb_vals = bb_vals .* grads[i,:]' _y = y[i, :] .+ sum(bb_vals, dims=2) - + #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ - if extra_points != 0 start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] @@ -443,7 +434,7 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no end R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) - for k in 1:size(ij)[1] #todo - speed up using @inbounds / @simd where required + for k in 1:size(ij)[1] R[ij[k,1],ij[k,2]]=r[k] R[ij[k,2],ij[k,1]]=r[k] end diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 594f3487e..85a7134b9 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -201,7 +201,6 @@ y_true = sphere_function.(x_test) g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - println("rmse: ", rmse) @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 end From bb14ca445a111203d6369f6051feb05ce3a18646 Mon Sep 17 00:00:00 2001 From: Vikram Date: Sun, 26 Jun 2022 12:06:54 +0530 Subject: [PATCH 08/12] update files after rebase --- docs/pages.jl | 16 +-- src/GEKPLS.jl | 349 ++++++++++++++++++++++++---------------------- src/Surrogates.jl | 2 +- test/GEKPLS.jl | 164 ++++++++++++---------- test/runtests.jl | 36 +++-- 5 files changed, 310 insertions(+), 257 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index 5101a1e11..895a72e02 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -4,7 +4,7 @@ pages = [ "Basics" => "tutorials.md", "Radials" => "radials.md", "Kriging" => "kriging.md", - "Gaussian Process" =>"abstractgps.md", + "Gaussian Process" => "abstractgps.md", "Lobachevsky" => "lobachevsky.md", "Linear" => "LinearSurrogate.md", "InverseDistance" => "InverseDistance.md", @@ -15,13 +15,13 @@ pages = [ "Polynomial Chaos" => "polychaos.md", "Variable Fidelity" => "variablefidelity.md", "Gradient Enhanced Kriging" => "gek.md", - "GEKPLS" => "gekpls.md" - ] + "GEKPLS" => "gekpls.md", + ] "User guide" => [ "Samples" => "samples.md", "Surrogates" => "surrogate.md", - "Optimization" => "optimizations.md" - ] + "Optimization" => "optimizations.md", + ] "Benchmarks" => [ "Sphere function" => "sphere_function.md", "Lp norm" => "lp.md", @@ -34,7 +34,7 @@ pages = [ "Ackley function" => "ackley.md", "Gramacy & Lee Function" => "gramacylee.md", "Salustowicz Benchmark" => "Salustowicz.md", - "Multi objective optimization" => "multi_objective_opt.md" - ] - "Contributing" => "contributing.md" + "Multi objective optimization" => "multi_objective_opt.md", ] + "Contributing" => "contributing.md" +] diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 1c8345e73..facd34a9b 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -6,13 +6,13 @@ mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate y::Matrix{T} #2 grads::Matrix{T} #3 xl::Matrix{T} #xlimits #4 - delta:: T #5 + delta::T #5 extra_points::Int #6 num_components::Int #7 beta::Vector{T} #8 gamma::Matrix{T} #9 theta::Vector{T} #10 - reduced_likelihood_function_value:: T #11 + reduced_likelihood_function_value::T #11 X_offset::Matrix{T} #12 X_scale::Matrix{T} #13 X_after_std::Matrix{T} #14 - X after standardization @@ -22,11 +22,11 @@ mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate end function bounds_error(x, xl) - num_x_rows = size(x,1) + num_x_rows = size(x, 1) num_dim = size(xl, 1) - for i in 1:num_x_rows - for j in 1:num_dim - if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2]) + for i = 1:num_x_rows + for j = 1:num_dim + if (x[i, j] < xl[j, 1] || x[i, j] > xl[j, 2]) return true end end @@ -36,23 +36,43 @@ end #constructor for GEKPLS Struct function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) - + #ensure that X values are within the upper and lower bounds if bounds_error(X, xlimits) println("X values outside bounds") return end - theta = [θ for i in 1:n_comp] - pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) - X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, y_after_PLS) + theta = [θ for i = 1:n_comp] + pls_mean, X_after_PLS, y_after_PLS = + _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) + X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = + standardization(X_after_PLS, y_after_PLS) D, ij = cross_distances(X_after_std) - pls_mean_reshaped = reshape(pls_mean, (size(X,2),n_comp)) + pls_mean_reshaped = reshape(pls_mean, (size(X, 2), n_comp)) d = componentwise_distance_PLS(D, "squar_exp", n_comp, pls_mean_reshaped) nt, nd = size(X_after_PLS) - beta, gamma, reduced_likelihood_function_value = _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) - return GEKPLS(X, y, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, reduced_likelihood_function_value, - X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) + beta, gamma, reduced_likelihood_function_value = + _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) + return GEKPLS( + X, + y, + grads, + xlimits, + delta_x, + extra_points, + n_comp, + beta, + gamma, + theta, + reduced_likelihood_function_value, + X_offset, + X_scale, + X_after_std, + pls_mean_reshaped, + y_mean, + y_std, + ) println("struct created") end @@ -63,9 +83,9 @@ function (g::GEKPLS)(X_test) X_cont = (X_test .- g.X_offset) ./ g.X_scale dx = differences(X_cont, g.X_after_std) pred_d = componentwise_distance_PLS(dx, "squar_exp", g.num_components, g.pls_mean) - nt = size(g.X_after_std,1) - r = transpose(reshape(squar_exp(g.theta, pred_d),(nt,n_eval))) - f = ones(n_eval,1) + nt = size(g.X_after_std, 1) + r = transpose(reshape(squar_exp(g.theta, pred_d), (nt, n_eval))) + f = ones(n_eval, 1) y_ = (f * g.beta) + (r * g.gamma) y = g.y_mean .+ g.y_std * y_ return y @@ -80,20 +100,23 @@ function add_point!(g::GEKPLS, new_x, new_y, new_grads) if bounds_error(new_x, g.xl) println("x values outside bounds") - return + return end g.x = vcat(g.x, new_x) g.y = vcat(g.y, new_y) g.grads = vcat(g.grads, new_grads) - pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) - g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, y_after_PLS) + pls_mean, X_after_PLS, y_after_PLS = + _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) + g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = + standardization(X_after_PLS, y_after_PLS) D, ij = cross_distances(g.X_after_std) g.pls_mean = reshape(pls_mean, (size(g.x, 2), g.num_components)) d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) nt, nd = size(X_after_PLS) - g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) + g.beta, g.gamma, g.reduced_likelihood_function_value = + _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) end function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) @@ -118,45 +141,45 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) # 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 nt, dim = size(X) - XX = zeros(0,dim) - yy = zeros(0,size(y)[2]) + XX = zeros(0, dim) + yy = zeros(0, size(y)[2]) coeff_pls = zeros((dim, n_comp, nt)) - - for i in 1:nt - if dim >= 3 - bb_vals = circshift(boxbehnken(dim, 1),1) + + for i = 1:nt + if dim >= 3 + bb_vals = circshift(boxbehnken(dim, 1), 1) else bb_vals = [ - 0.0 0.0; #center - 1.0 0.0; #right - 0.0 1.0; #up - -1.0 0.0; #left - 0.0 -1.0; #down - 1.0 1.0; #right up - -1.0 1.0; #left up - -1.0 -1.0; #left down - 1.0 -1.0; #right down - ] + 0.0 0.0 #center + 1.0 0.0 #right + 0.0 1.0 #up + -1.0 0.0 #left + 0.0 -1.0 #down + 1.0 1.0 #right up + -1.0 1.0 #left up + -1.0 -1.0 #left down + 1.0 -1.0 #right down + ] end - _X = zeros((size(bb_vals)[1], dim)) - _y = zeros((size(bb_vals)[1], 1)) + _X = zeros((size(bb_vals)[1], dim)) + _y = zeros((size(bb_vals)[1], 1)) bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals - _X = X[i, :]' .+ bb_vals - bb_vals = bb_vals .* grads[i,:]' - _y = y[i, :] .+ sum(bb_vals, dims=2) - + _X = X[i, :]' .+ bb_vals + bb_vals = bb_vals .* grads[i, :]' + _y = y[i, :] .+ sum(bb_vals, dims = 2) + #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ if extra_points != 0 - start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) - max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] + start_index = max(1, length(coeff_pls[:, 1, i]) - extra_points + 1) + max_coeff = sortperm(broadcast(abs, coeff_pls[:, 1, i]))[start_index:end] for ii in max_coeff XX = [XX; transpose(X[i, :])] - XX[end, ii] += delta_x * (xlimits[ii,2]-xlimits[ii,1]) + XX[end, ii] += delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) yy = [yy; y[i]] - yy[end] += grads[i,ii] * delta_x * (xlimits[ii,2]-xlimits[ii,1]) + yy[end] += grads[i, ii] * delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) end end end @@ -165,89 +188,89 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) y = [y; yy] end - pls_mean = mean(broadcast(abs,coeff_pls),dims=3) + pls_mean = mean(broadcast(abs, coeff_pls), dims = 3) return pls_mean, X, y end ######start of bbdesign###### - # - # Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' - # https://github.com/phrb/ExperimentalDesign.jl +# +# Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' +# https://github.com/phrb/ExperimentalDesign.jl - # MIT License +# MIT License - # ExperimentalDesign.jl: Design of Experiments in Julia - # Copyright (C) 2019 Pedro Bruel +# ExperimentalDesign.jl: Design of Experiments in Julia +# Copyright (C) 2019 Pedro Bruel - # Permission is hereby granted, free of charge, to any person obtaining a copy of - # this software and associated documentation files (the "Software"), to deal in - # the Software without restriction, including without limitation the rights to - # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of - # the Software, and to permit persons to whom the Software is furnished to do so, - # subject to the following conditions: +# Permission is hereby granted, free of charge, to any person obtaining a copy of +# this software and associated documentation files (the "Software"), to deal in +# the Software without restriction, including without limitation the rights to +# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +# the Software, and to permit persons to whom the Software is furnished to do so, +# subject to the following conditions: - # The above copyright notice and this permission notice (including the next - # paragraph) shall be included in all copies or substantial portions of the - # Software. +# The above copyright notice and this permission notice (including the next +# paragraph) shall be included in all copies or substantial portions of the +# Software. - # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS - # FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR - # 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. - # +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +# 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) - end - - function boxbehnken(matrix_size::Int, center::Int) - @assert matrix_size>=3 - - A_fact = explicit_fullfactorial(Tuple([-1,1] for i = 1:2)) - - rows = floor(Int, (0.5*matrix_size*(matrix_size-1))*size(A_fact)[1]) - - A = zeros(rows, matrix_size) - - l = 0 - for i in 1:matrix_size-1 - for j in i+1:matrix_size - l = l +1 - A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] - A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] - end - end - - if center == matrix_size - if matrix_size <= 16 - points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] - center = points[matrix_size] - end +function boxbehnken(matrix_size::Int) + boxbehnken(matrix_size, matrix_size) +end + +function boxbehnken(matrix_size::Int, center::Int) + @assert matrix_size >= 3 + + A_fact = explicit_fullfactorial(Tuple([-1, 1] for i = 1:2)) + + rows = floor(Int, (0.5 * matrix_size * (matrix_size - 1)) * size(A_fact)[1]) + + A = zeros(rows, matrix_size) + + l = 0 + for i = 1:matrix_size-1 + for j = i+1:matrix_size + l = l + 1 + A[max(0, (l - 1) * size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] + A[max(0, (l - 1) * size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] end - - A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) - end - - function explicit_fullfactorial(factors::Tuple) - explicit_fullfactorial(fullfactorial(factors)) end - - function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) - hcat(vcat.(collect(iterator)...)...) - end - - function fullfactorial(factors::Tuple) - Base.Iterators.product(factors...) + + if center == matrix_size + if matrix_size <= 16 + points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] + center = points[matrix_size] + end end - + + A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) +end + +function explicit_fullfactorial(factors::Tuple) + explicit_fullfactorial(fullfactorial(factors)) +end + +function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) + hcat(vcat.(collect(iterator)...)...) +end + +function fullfactorial(factors::Tuple) + Base.Iterators.product(factors...) +end + ######end of bb design###### -function standardization(X,y) +function standardization(X, y) """ We substract the mean from each variable. Then, we divide the values of each variable by its standard deviation. @@ -279,11 +302,11 @@ function standardization(X,y) #Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L21 X_offset = mean(X, dims = 1) X_scale = std(X, dims = 1) - X_scale = map(x -> (x==0.0) ? x=1 : x=x, X_scale) #to prevent division by 0 below + X_scale = map(x -> (x == 0.0) ? x = 1 : x = x, X_scale) #to prevent division by 0 below y_mean = mean(y) y_std = std(y) - y_std = map(y -> (y==0) ? y=1 : y=y, y_std) #to prevent division by 0 below - X = (X.-X_offset) ./ X_scale + y_std = map(y -> (y == 0) ? y = 1 : y = y, y_std) #to prevent division by 0 below + X = (X .- X_offset) ./ X_scale y = (y .- y_mean) ./ y_std return X, y, X_offset, y_mean, X_scale, y_std @@ -310,18 +333,18 @@ function cross_distances(X) """ # equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L86 n_samples, n_features = size(X) - n_nonzero_cross_dist = ( n_samples * (n_samples - 1) ) ÷ 2 + n_nonzero_cross_dist = (n_samples * (n_samples - 1)) ÷ 2 ij = zeros((n_nonzero_cross_dist, 2)) D = zeros((n_nonzero_cross_dist, n_features)) ll_1 = 0 - - for k in 1:n_samples - 1 + + for k = 1:n_samples-1 ll_0 = ll_1 + 1 ll_1 = ll_0 + n_samples - k - 1 ij[ll_0:ll_1, 1] .= k ij[ll_0:ll_1, 2] = k+1:1:n_samples - D[ll_0:ll_1, :] = -(X[(k + 1) : n_samples,:] .- X[k,:]') - + D[ll_0:ll_1, :] = -(X[(k+1):n_samples, :] .- X[k, :]') + end return D, Int.(ij) end @@ -366,7 +389,7 @@ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) D_corr = zeros((size(D)[1], n_comp)) if corr == "squar_exp" - D_corr = D.^2 * coeff_pls.^2 + D_corr = D .^ 2 * coeff_pls .^ 2 else #abs_exp D_corr = abs.(D) * abs.(coeff_pls) end @@ -386,27 +409,27 @@ function squar_exp(theta, d) Returns: -------- r: array containing the values of the autocorrelation model - + """ n_components = size(d)[2] - theta = reshape(theta, (1,n_components)) - return exp.(-sum(theta .* d, dims=2)) + theta = reshape(theta, (1, n_components)) + return exp.(-sum(theta .* d, dims = 2)) end function differences(X, Y) #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L392 #code credit: Elias Carvalho - https://stackoverflow.com/questions/72392010/row-wise-operations-between-matrices-in-julia - Rx = repeat(X, inner=(size(Y, 1), 1)) + Rx = repeat(X, inner = (size(Y, 1), 1)) Ry = repeat(Y, size(X, 1)) return Rx - Ry end -function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise=0.0) +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 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. - + Parameters ---------- theta: array containing the parameters at which the Gaussian Process model parameters should be determined. @@ -429,29 +452,29 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 reduced_likelihood_function_value = -Inf nugget = 1000000.0 * eps() #a jitter for numerical stability; reducing the multiple from 1000000.0 results in positive definite error for Cholesky decomposition; - if kernel_type =="squar_exp" #todo - add other kernel type abs_exp etc. - r = squar_exp(theta,d) + if kernel_type == "squar_exp" #todo - add other kernel type abs_exp etc. + r = squar_exp(theta, d) end - R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) + R = (I + zeros(nt, nt)) .* (1.0 + nugget + noise) - for k in 1:size(ij)[1] - R[ij[k,1],ij[k,2]]=r[k] - R[ij[k,2],ij[k,1]]=r[k] + for k = 1:size(ij)[1] + R[ij[k, 1], ij[k, 2]] = r[k] + R[ij[k, 2], ij[k, 1]] = r[k] end C = cholesky(R).L #todo - #values diverge at this point from SMT code; verify impact - F = ones(nt,1) #todo - examine if this should be a parameter for this function - Ft = C\F + F = ones(nt, 1) #todo - examine if this should be a parameter for this function + Ft = C \ F Q, G = qr(Ft) Q = Array(Q) - Yt = C\y_norma + Yt = C \ y_norma #todo - in smt, they check if the matrix is ill-conditioned using SVD. Verify and include if necessary beta = G \ [(transpose(Q) ⋅ Yt)] rho = Yt .- (Ft .* beta) gamma = transpose(C) \ rho - sigma2 = sum((rho).^2, dims=1)/nt - detR = prod(diag(C).^(2.0/nt)) - reduced_likelihood_function_value = - nt * log10(sum(sigma2)) - nt * log10(detR) + sigma2 = sum((rho) .^ 2, dims = 1) / nt + detR = prod(diag(C) .^ (2.0 / nt)) + reduced_likelihood_function_value = -nt * log10(sum(sigma2)) - nt * log10(detR) return beta, gamma, reduced_likelihood_function_value end @@ -463,34 +486,34 @@ end # https://github.com/scikit-learn/scikit-learn/blob/80598905e/sklearn/cross_decomposition/_pls.py -function _center_scale(X,Y) +function _center_scale(X, Y) x_mean = mean(X, dims = 1) X .-= x_mean y_mean = mean(Y, dims = 1) Y .-= y_mean - x_std = std(X, dims = 1) - x_std[x_std .== 0] .= 1.0 + x_std = std(X, dims = 1) + x_std[x_std.==0] .= 1.0 X ./= x_std - y_std = std(Y, dims = 1) - y_std[y_std .==0] .= 1.0 + y_std = std(Y, dims = 1) + y_std[y_std.==0] .= 1.0 Y ./= y_std - return X,Y + return X, Y end -function _svd_flip_1d(u,v) +function _svd_flip_1d(u, v) # equivalent of https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L149 biggest_abs_val_idx = findmax(abs.(vec(u)))[2] sign_ = sign(u[biggest_abs_val_idx]) u .*= sign_ v .*= sign_ - + end -function _get_first_singular_vectors_power_method(X,Y) +function _get_first_singular_vectors_power_method(X, Y) my_eps = eps() y_score = vec(Y) x_weights = transpose(X)y_score / dot(y_score, y_score) - x_weights ./= (sqrt(dot(x_weights,x_weights)) + my_eps) + x_weights ./= (sqrt(dot(x_weights, x_weights)) + my_eps) x_score = X * x_weights y_weights = transpose(Y)x_score / dot(x_score, x_score) y_score = Y * y_weights / (dot(y_weights, y_weights) + my_eps) @@ -501,28 +524,28 @@ function _get_first_singular_vectors_power_method(X,Y) return x_weights, y_weights end -function _modified_pls(X,Y, n_components) - x_weights_ = zeros(size(X,2),n_components) - _x_scores = zeros(size(X,1), n_components) - x_loadings_ = zeros(size(X,2),n_components) - Xk, Yk = _center_scale(X,Y) +function _modified_pls(X, Y, n_components) + x_weights_ = zeros(size(X, 2), n_components) + _x_scores = zeros(size(X, 1), n_components) + x_loadings_ = zeros(size(X, 2), n_components) + Xk, Yk = _center_scale(X, Y) - for k in 1:n_components - x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) + for k = 1:n_components + x_weights, y_weights = _get_first_singular_vectors_power_method(Xk, Yk) if x_weights == false break - end + end - _svd_flip_1d(x_weights, y_weights) + _svd_flip_1d(x_weights, y_weights) x_scores = Xk * x_weights x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) Xk = Xk - (x_scores * x_loadings) - y_loadings = transpose(x_scores)*Yk / dot(x_scores, x_scores) - Yk = Yk - x_scores*y_loadings - x_weights_[:,k] = x_weights - _x_scores[:,k] = x_scores - x_loadings_[:,k] = vec(x_loadings) + y_loadings = transpose(x_scores) * Yk / dot(x_scores, x_scores) + Yk = Yk - x_scores * y_loadings + x_weights_[:, k] = x_weights + _x_scores[:, k] = x_scores + x_loadings_[:, k] = vec(x_loadings) end x_rotations_ = x_weights_ * pinv(transpose(x_loadings_)x_weights_) diff --git a/src/Surrogates.jl b/src/Surrogates.jl index 8d418ab3f..eaf50ed44 100644 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -81,7 +81,7 @@ function PolyChaosStructure(; op) end export current_surrogates -export GEKPLS +export GEKPLS export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 85a7134b9..d3e91f366 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -7,9 +7,9 @@ function vector_of_tuples_to_matrix(v) num_rows = length(v) num_cols = length(first(v)) K = zeros(num_rows, num_cols) - for row in 1:num_rows - for col in 1:num_cols - K[row, col]=v[row][col] + for row = 1:num_rows + for col = 1:num_cols + K[row, col] = v[row][col] end end return K @@ -20,8 +20,8 @@ function vector_of_tuples_to_matrix2(v) num_rows = length(v) num_cols = length(first(first(v))) K = zeros(num_rows, num_cols) - for row in 1:num_rows - for col in 1:num_cols + for row = 1:num_rows + for col = 1:num_cols K[row, col] = v[row][1][col] end end @@ -38,55 +38,56 @@ function water_flow(x) H_l = x[6] L = x[7] K_w = x[8] - log_val = log(r/r_w) - return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) + log_val = log(r / r_w) + return (2 * pi * T_u * (H_u - H_l)) / + (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l)) end n = 1000 d = 8 -lb = [0.05,100,63070,990,63.1,700,1120,9855] -ub = [0.15,50000,115600,1110,116,820,1680,12045] -x = sample(n,lb,ub,SobolSample()) +lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] +ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x)) -y = reshape(water_flow.(x),(size(x,1),1)) +y = reshape(water_flow.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = water_flow.(x_test) -@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin +@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.03, atol=0.02) #rmse: 0.039 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.03, atol = 0.02) #rmse: 0.039 end -@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin +@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 end -@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin +@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin n_comp = 3 delta_x = 0.0001 extra_points = 3 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 end ## welded beam tests @@ -94,77 +95,79 @@ function welded_beam(x) h = x[1] l = x[2] t = x[3] - a = 6000/(sqrt(2)*h*l) - b = (6000*(14+0.5*l)*sqrt(0.25*(l^2+(h+t)^2)))/(2*(0.707*h*l*(l^2/12 + 0.25*(h+t)^2))) - return (sqrt(a^2+b^2 + l*a*b))/(sqrt(0.25*(l^2+(h+t)^2))) + a = 6000 / (sqrt(2) * h * l) + b = + (6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))) / + (2 * (0.707 * h * l * (l^2 / 12 + 0.25 * (h + t)^2))) + return (sqrt(a^2 + b^2 + l * a * b)) / (sqrt(0.25 * (l^2 + (h + t)^2))) end n = 1000 d = 3 -lb = [0.125,5.0,5.0] -ub = [1.,10.,10.] -x = sample(n,lb,ub,SobolSample()) +lb = [0.125, 5.0, 5.0] +ub = [1.0, 10.0, 10.0] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(welded_beam, x)) -y = reshape(welded_beam.(x),(size(x,1),1)) +y = reshape(welded_beam.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = welded_beam.(x_test) -@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin +@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 39.0, atol=0.5) #rmse: 38.988 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 39.0, atol = 0.5) #rmse: 38.988 end -@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin +@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 39.5, atol=0.5) #rmse: 39.481 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 39.5, atol = 0.5) #rmse: 39.481 end ## increasing extra points increases accuracy -@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin +@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin n_comp = 2 delta_x = 0.0001 extra_points = 4 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 37.5, atol=0.5) #rmse: 37.87 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 37.5, atol = 0.5) #rmse: 37.87 end ## sphere function tests function sphere_function(x) - return sum(x.^2) + return sum(x .^ 2) end ## 3D n = 100 d = 3 lb = [-5.0, -5.0, -5.0] -ub = [5.0, 5.0 ,5.0] -x = sample(n,lb,ub,SobolSample()) +ub = [5.0, 5.0, 5.0] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x),(size(x,1),1)) +y = reshape(sphere_function.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 7: Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin @@ -172,25 +175,25 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.001, atol=0.05) #rmse: 0.00083 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.001, atol = 0.05) #rmse: 0.00083 end ## 2D -n = 50 +n = 50 d = 2 lb = [-10.0, -10.0] ub = [10.0, 10.0] -x = sample(n,lb,ub,SobolSample()) +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x),(size(x,1),1)) +y = reshape(sphere_function.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 10 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 10 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 8: Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2" begin @@ -198,16 +201,16 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.1, atol = 0.5) #rmse: 0.0022 end @testset "Test 9: Add Point Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin #first we create a surrogate model with just 3 input points initial_x_vec = [(1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0)] - initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec,1),1)) + initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec, 1), 1)) initial_X = vector_of_tuples_to_matrix(initial_x_vec) initial_grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, initial_x_vec)) lb = [-5.0, -5.0, -5.0] @@ -217,22 +220,31 @@ end delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(initial_X, initial_y, initial_grads, n_comp, delta_x, xlimits, extra_points, initial_theta) - n_test = 100 - x_test = sample(n_test,lb,ub,GoldenSample()) - X_test = vector_of_tuples_to_matrix(x_test) - y_true = sphere_function.(x_test) + g = GEKPLS( + initial_X, + initial_y, + initial_grads, + n_comp, + delta_x, + xlimits, + extra_points, + initial_theta, + ) + n_test = 100 + x_test = sample(n_test, lb, ub, GoldenSample()) + X_test = vector_of_tuples_to_matrix(x_test) + y_true = sphere_function.(x_test) y_pred1 = g(X_test) - rmse1 = sqrt(sum(((y_pred1 - y_true).^2)/n_test)) #rmse1 = 31.91 + rmse1 = sqrt(sum(((y_pred1 - y_true) .^ 2) / n_test)) #rmse1 = 31.91 #then we update the model with more points to see if performance improves n = 100 - x = sample(n,lb,ub,SobolSample()) + x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) - y = reshape(sphere_function.(x),(size(x,1),1)) + y = reshape(sphere_function.(x), (size(x, 1), 1)) add_point!(g, X, y, grads) y_pred2 = g(X_test) - rmse2 = sqrt(sum(((y_pred2 - y_true).^2)/n_test)) #rmse2 = 0.0015 + rmse2 = sqrt(sum(((y_pred2 - y_true) .^ 2) / n_test)) #rmse2 = 0.0015 @test (rmse2 < rmse1) end diff --git a/test/runtests.jl b/test/runtests.jl index 5f32d5058..7ff07219f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,15 +15,33 @@ for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", end end -@time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end -@time @safetestset "Radials.jl" begin include("Radials.jl") end -@time @safetestset "Kriging.jl" begin include("Kriging.jl") end -@time @safetestset "Sampling" begin include("sampling.jl") end -@time @safetestset "Optimization" begin include("optimization.jl") end -@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end -@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end -@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end -@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end +@time @safetestset "GEKPLS.jl" begin + include("GEKPLS.jl") +end +@time @safetestset "Radials.jl" begin + include("Radials.jl") +end +@time @safetestset "Kriging.jl" begin + include("Kriging.jl") +end +@time @safetestset "Sampling" begin + include("sampling.jl") +end +@time @safetestset "Optimization" begin + include("optimization.jl") +end +@time @safetestset "LinearSurrogate" begin + include("linearSurrogate.jl") +end +@time @safetestset "Lobachevsky" begin + include("lobachevsky.jl") +end +@time @safetestset "InverseDistanceSurrogate" begin + include("inverseDistanceSurrogate.jl") +end +@time @safetestset "SecondOrderPolynomialSurrogate" begin + include("secondOrderPolynomialSurrogate.jl") +end # @time @safetestset "AD_Compatibility" begin include("AD_compatibility.jl") end @time @safetestset "Wendland" begin include("Wendland.jl") end @time @safetestset "VariableFidelity" begin include("VariableFidelity.jl") end From 0b1984f74568c65a6aac78ceafc5396dffaa869f Mon Sep 17 00:00:00 2001 From: Vikram Date: Mon, 27 Jun 2022 20:00:16 +0530 Subject: [PATCH 09/12] remove unnecessary variables and comments --- src/GEKPLS.jl | 7 ++----- test/GEKPLS.jl | 4 ---- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index facd34a9b..48e1e801a 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -73,7 +73,7 @@ function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) y_mean, y_std, ) - println("struct created") + end @@ -140,6 +140,7 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) # 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 + nt, dim = size(X) XX = zeros(0, dim) yy = zeros(0, size(y)[2]) @@ -167,10 +168,6 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) _X = X[i, :]' .+ bb_vals bb_vals = bb_vals .* grads[i, :]' _y = y[i, :] .+ sum(bb_vals, dims = 2) - - #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. - #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. - coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ if extra_points != 0 start_index = max(1, length(coeff_pls[:, 1, i]) - extra_points + 1) diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index d3e91f366..fcda8b58b 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -44,7 +44,6 @@ function water_flow(x) end n = 1000 -d = 8 lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] x = sample(n, lb, ub, SobolSample()) @@ -103,7 +102,6 @@ function welded_beam(x) end n = 1000 -d = 3 lb = [0.125, 5.0, 5.0] ub = [1.0, 10.0, 10.0] x = sample(n, lb, ub, SobolSample()) @@ -157,7 +155,6 @@ end ## 3D n = 100 -d = 3 lb = [-5.0, -5.0, -5.0] ub = [5.0, 5.0, 5.0] x = sample(n, lb, ub, SobolSample()) @@ -183,7 +180,6 @@ end ## 2D n = 50 -d = 2 lb = [-10.0, -10.0] ub = [10.0, 10.0] x = sample(n, lb, ub, SobolSample()) From 45d864d0a5593b93848000b3620d4ce136e2b51c Mon Sep 17 00:00:00 2001 From: Vikram Date: Tue, 28 Jun 2022 18:22:38 +0530 Subject: [PATCH 10/12] fix changes introduced by rebase --- docs/pages.jl | 16 +-- src/GEKPLS.jl | 353 ++++++++++++++++++++++------------------------ src/Surrogates.jl | 16 +-- test/GEKPLS.jl | 167 +++++++++++----------- test/runtests.jl | 40 ++---- 5 files changed, 273 insertions(+), 319 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index 895a72e02..5101a1e11 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -4,7 +4,7 @@ pages = [ "Basics" => "tutorials.md", "Radials" => "radials.md", "Kriging" => "kriging.md", - "Gaussian Process" => "abstractgps.md", + "Gaussian Process" =>"abstractgps.md", "Lobachevsky" => "lobachevsky.md", "Linear" => "LinearSurrogate.md", "InverseDistance" => "InverseDistance.md", @@ -15,13 +15,13 @@ pages = [ "Polynomial Chaos" => "polychaos.md", "Variable Fidelity" => "variablefidelity.md", "Gradient Enhanced Kriging" => "gek.md", - "GEKPLS" => "gekpls.md", - ] + "GEKPLS" => "gekpls.md" + ] "User guide" => [ "Samples" => "samples.md", "Surrogates" => "surrogate.md", - "Optimization" => "optimizations.md", - ] + "Optimization" => "optimizations.md" + ] "Benchmarks" => [ "Sphere function" => "sphere_function.md", "Lp norm" => "lp.md", @@ -34,7 +34,7 @@ pages = [ "Ackley function" => "ackley.md", "Gramacy & Lee Function" => "gramacylee.md", "Salustowicz Benchmark" => "Salustowicz.md", - "Multi objective optimization" => "multi_objective_opt.md", - ] + "Multi objective optimization" => "multi_objective_opt.md" + ] "Contributing" => "contributing.md" -] + ] diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 48e1e801a..c228dea9b 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -6,13 +6,13 @@ mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate y::Matrix{T} #2 grads::Matrix{T} #3 xl::Matrix{T} #xlimits #4 - delta::T #5 + delta:: T #5 extra_points::Int #6 num_components::Int #7 beta::Vector{T} #8 gamma::Matrix{T} #9 theta::Vector{T} #10 - reduced_likelihood_function_value::T #11 + reduced_likelihood_function_value:: T #11 X_offset::Matrix{T} #12 X_scale::Matrix{T} #13 X_after_std::Matrix{T} #14 - X after standardization @@ -22,11 +22,11 @@ mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate end function bounds_error(x, xl) - num_x_rows = size(x, 1) + num_x_rows = size(x,1) num_dim = size(xl, 1) - for i = 1:num_x_rows - for j = 1:num_dim - if (x[i, j] < xl[j, 1] || x[i, j] > xl[j, 2]) + for i in 1:num_x_rows + for j in 1:num_dim + if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2]) return true end end @@ -36,44 +36,24 @@ end #constructor for GEKPLS Struct function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) - + #ensure that X values are within the upper and lower bounds if bounds_error(X, xlimits) println("X values outside bounds") return end - theta = [θ for i = 1:n_comp] - pls_mean, X_after_PLS, y_after_PLS = - _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) - X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = - standardization(X_after_PLS, y_after_PLS) + theta = [θ for i in 1:n_comp] + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) + X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, y_after_PLS) D, ij = cross_distances(X_after_std) - pls_mean_reshaped = reshape(pls_mean, (size(X, 2), n_comp)) + pls_mean_reshaped = reshape(pls_mean, (size(X,2),n_comp)) d = componentwise_distance_PLS(D, "squar_exp", n_comp, pls_mean_reshaped) nt, nd = size(X_after_PLS) - beta, gamma, reduced_likelihood_function_value = - _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) - return GEKPLS( - X, - y, - grads, - xlimits, - delta_x, - extra_points, - n_comp, - beta, - gamma, - theta, - reduced_likelihood_function_value, - X_offset, - X_scale, - X_after_std, - pls_mean_reshaped, - y_mean, - y_std, - ) - + beta, gamma, reduced_likelihood_function_value = _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) + return GEKPLS(X, y, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, reduced_likelihood_function_value, + X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) + println("struct created") end @@ -83,9 +63,9 @@ function (g::GEKPLS)(X_test) X_cont = (X_test .- g.X_offset) ./ g.X_scale dx = differences(X_cont, g.X_after_std) pred_d = componentwise_distance_PLS(dx, "squar_exp", g.num_components, g.pls_mean) - nt = size(g.X_after_std, 1) - r = transpose(reshape(squar_exp(g.theta, pred_d), (nt, n_eval))) - f = ones(n_eval, 1) + nt = size(g.X_after_std,1) + r = transpose(reshape(squar_exp(g.theta, pred_d),(nt,n_eval))) + f = ones(n_eval,1) y_ = (f * g.beta) + (r * g.gamma) y = g.y_mean .+ g.y_std * y_ return y @@ -100,23 +80,20 @@ function add_point!(g::GEKPLS, new_x, new_y, new_grads) if bounds_error(new_x, g.xl) println("x values outside bounds") - return + return end g.x = vcat(g.x, new_x) g.y = vcat(g.y, new_y) g.grads = vcat(g.grads, new_grads) - pls_mean, X_after_PLS, y_after_PLS = - _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) - g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = - standardization(X_after_PLS, y_after_PLS) + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) + g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, y_after_PLS) D, ij = cross_distances(g.X_after_std) g.pls_mean = reshape(pls_mean, (size(g.x, 2), g.num_components)) d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) nt, nd = size(X_after_PLS) - g.beta, g.gamma, g.reduced_likelihood_function_value = - _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) + g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) end function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) @@ -142,41 +119,45 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) # and https://github.com/SMTorg/smt/blob/f124c01ffa78c04b80221dded278a20123dac742/smt/surrogate_models/gekpls.py#L48 nt, dim = size(X) - XX = zeros(0, dim) - yy = zeros(0, size(y)[2]) + XX = zeros(0,dim) + yy = zeros(0,size(y)[2]) coeff_pls = zeros((dim, n_comp, nt)) - - for i = 1:nt - if dim >= 3 - bb_vals = circshift(boxbehnken(dim, 1), 1) + + for i in 1:nt + if dim >= 3 + bb_vals = circshift(boxbehnken(dim, 1),1) else bb_vals = [ - 0.0 0.0 #center - 1.0 0.0 #right - 0.0 1.0 #up - -1.0 0.0 #left - 0.0 -1.0 #down - 1.0 1.0 #right up - -1.0 1.0 #left up - -1.0 -1.0 #left down - 1.0 -1.0 #right down - ] + 0.0 0.0; #center + 1.0 0.0; #right + 0.0 1.0; #up + -1.0 0.0; #left + 0.0 -1.0; #down + 1.0 1.0; #right up + -1.0 1.0; #left up + -1.0 -1.0; #left down + 1.0 -1.0; #right down + ] end - _X = zeros((size(bb_vals)[1], dim)) - _y = zeros((size(bb_vals)[1], 1)) + _X = zeros((size(bb_vals)[1], dim)) + _y = zeros((size(bb_vals)[1], 1)) bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals - _X = X[i, :]' .+ bb_vals - bb_vals = bb_vals .* grads[i, :]' - _y = y[i, :] .+ sum(bb_vals, dims = 2) + _X = X[i, :]' .+ bb_vals + bb_vals = bb_vals .* grads[i,:]' + _y = y[i, :] .+ sum(bb_vals, dims=2) + + #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. + #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. + coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ if extra_points != 0 - start_index = max(1, length(coeff_pls[:, 1, i]) - extra_points + 1) - max_coeff = sortperm(broadcast(abs, coeff_pls[:, 1, i]))[start_index:end] + start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) + max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] for ii in max_coeff XX = [XX; transpose(X[i, :])] - XX[end, ii] += delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) + XX[end, ii] += delta_x * (xlimits[ii,2]-xlimits[ii,1]) yy = [yy; y[i]] - yy[end] += grads[i, ii] * delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) + yy[end] += grads[i,ii] * delta_x * (xlimits[ii,2]-xlimits[ii,1]) end end end @@ -185,89 +166,89 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) y = [y; yy] end - pls_mean = mean(broadcast(abs, coeff_pls), dims = 3) + pls_mean = mean(broadcast(abs,coeff_pls),dims=3) return pls_mean, X, y end ######start of bbdesign###### -# -# Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' -# https://github.com/phrb/ExperimentalDesign.jl + # + # Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' + # https://github.com/phrb/ExperimentalDesign.jl -# MIT License - -# ExperimentalDesign.jl: Design of Experiments in Julia -# Copyright (C) 2019 Pedro Bruel - -# Permission is hereby granted, free of charge, to any person obtaining a copy of -# this software and associated documentation files (the "Software"), to deal in -# the Software without restriction, including without limitation the rights to -# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of -# the Software, and to permit persons to whom the Software is furnished to do so, -# subject to the following conditions: - -# The above copyright notice and this permission notice (including the next -# paragraph) shall be included in all copies or substantial portions of the -# Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS -# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR -# 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) -end + # MIT License -function boxbehnken(matrix_size::Int, center::Int) - @assert matrix_size >= 3 + # ExperimentalDesign.jl: Design of Experiments in Julia + # Copyright (C) 2019 Pedro Bruel - A_fact = explicit_fullfactorial(Tuple([-1, 1] for i = 1:2)) + # Permission is hereby granted, free of charge, to any person obtaining a copy of + # this software and associated documentation files (the "Software"), to deal in + # the Software without restriction, including without limitation the rights to + # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + # the Software, and to permit persons to whom the Software is furnished to do so, + # subject to the following conditions: - rows = floor(Int, (0.5 * matrix_size * (matrix_size - 1)) * size(A_fact)[1]) + # The above copyright notice and this permission notice (including the next + # paragraph) shall be included in all copies or substantial portions of the + # Software. - A = zeros(rows, matrix_size) + # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS + # FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR + # 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. + # - l = 0 - for i = 1:matrix_size-1 - for j = i+1:matrix_size - l = l + 1 - A[max(0, (l - 1) * size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] - A[max(0, (l - 1) * size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] - end + function boxbehnken(matrix_size::Int) + boxbehnken(matrix_size, matrix_size) end - - if center == matrix_size - if matrix_size <= 16 - points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] - center = points[matrix_size] + + function boxbehnken(matrix_size::Int, center::Int) + @assert matrix_size>=3 + + A_fact = explicit_fullfactorial(Tuple([-1,1] for i = 1:2)) + + rows = floor(Int, (0.5*matrix_size*(matrix_size-1))*size(A_fact)[1]) + + A = zeros(rows, matrix_size) + + l = 0 + for i in 1:matrix_size-1 + for j in i+1:matrix_size + l = l +1 + A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] + A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] + end end + + if center == matrix_size + if matrix_size <= 16 + points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] + center = points[matrix_size] + end + end + + A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) end - - A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) -end - -function explicit_fullfactorial(factors::Tuple) - explicit_fullfactorial(fullfactorial(factors)) -end - -function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) - hcat(vcat.(collect(iterator)...)...) -end - -function fullfactorial(factors::Tuple) - Base.Iterators.product(factors...) -end - + + function explicit_fullfactorial(factors::Tuple) + explicit_fullfactorial(fullfactorial(factors)) + end + + function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) + hcat(vcat.(collect(iterator)...)...) + end + + function fullfactorial(factors::Tuple) + Base.Iterators.product(factors...) + end + ######end of bb design###### -function standardization(X, y) +function standardization(X,y) """ We substract the mean from each variable. Then, we divide the values of each variable by its standard deviation. @@ -299,11 +280,11 @@ function standardization(X, y) #Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L21 X_offset = mean(X, dims = 1) X_scale = std(X, dims = 1) - X_scale = map(x -> (x == 0.0) ? x = 1 : x = x, X_scale) #to prevent division by 0 below + X_scale = map(x -> (x==0.0) ? x=1 : x=x, X_scale) #to prevent division by 0 below y_mean = mean(y) y_std = std(y) - y_std = map(y -> (y == 0) ? y = 1 : y = y, y_std) #to prevent division by 0 below - X = (X .- X_offset) ./ X_scale + y_std = map(y -> (y==0) ? y=1 : y=y, y_std) #to prevent division by 0 below + X = (X.-X_offset) ./ X_scale y = (y .- y_mean) ./ y_std return X, y, X_offset, y_mean, X_scale, y_std @@ -330,18 +311,18 @@ function cross_distances(X) """ # equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L86 n_samples, n_features = size(X) - n_nonzero_cross_dist = (n_samples * (n_samples - 1)) ÷ 2 + n_nonzero_cross_dist = ( n_samples * (n_samples - 1) ) ÷ 2 ij = zeros((n_nonzero_cross_dist, 2)) D = zeros((n_nonzero_cross_dist, n_features)) ll_1 = 0 - - for k = 1:n_samples-1 + + for k in 1:n_samples - 1 ll_0 = ll_1 + 1 ll_1 = ll_0 + n_samples - k - 1 ij[ll_0:ll_1, 1] .= k ij[ll_0:ll_1, 2] = k+1:1:n_samples - D[ll_0:ll_1, :] = -(X[(k+1):n_samples, :] .- X[k, :]') - + D[ll_0:ll_1, :] = -(X[(k + 1) : n_samples,:] .- X[k,:]') + end return D, Int.(ij) end @@ -386,7 +367,7 @@ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) D_corr = zeros((size(D)[1], n_comp)) if corr == "squar_exp" - D_corr = D .^ 2 * coeff_pls .^ 2 + D_corr = D.^2 * coeff_pls.^2 else #abs_exp D_corr = abs.(D) * abs.(coeff_pls) end @@ -406,27 +387,27 @@ function squar_exp(theta, d) Returns: -------- r: array containing the values of the autocorrelation model - + """ n_components = size(d)[2] - theta = reshape(theta, (1, n_components)) - return exp.(-sum(theta .* d, dims = 2)) + theta = reshape(theta, (1,n_components)) + return exp.(-sum(theta .* d, dims=2)) end function differences(X, Y) #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L392 #code credit: Elias Carvalho - https://stackoverflow.com/questions/72392010/row-wise-operations-between-matrices-in-julia - Rx = repeat(X, inner = (size(Y, 1), 1)) + Rx = repeat(X, inner=(size(Y, 1), 1)) Ry = repeat(Y, size(X, 1)) return Rx - Ry end -function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise = 0.0) +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 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. - + Parameters ---------- theta: array containing the parameters at which the Gaussian Process model parameters should be determined. @@ -449,29 +430,29 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 reduced_likelihood_function_value = -Inf nugget = 1000000.0 * eps() #a jitter for numerical stability; reducing the multiple from 1000000.0 results in positive definite error for Cholesky decomposition; - if kernel_type == "squar_exp" #todo - add other kernel type abs_exp etc. - r = squar_exp(theta, d) + if kernel_type =="squar_exp" #todo - add other kernel type abs_exp etc. + r = squar_exp(theta,d) end - R = (I + zeros(nt, nt)) .* (1.0 + nugget + noise) + R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) - for k = 1:size(ij)[1] - R[ij[k, 1], ij[k, 2]] = r[k] - R[ij[k, 2], ij[k, 1]] = r[k] + for k in 1:size(ij)[1] + R[ij[k,1],ij[k,2]]=r[k] + R[ij[k,2],ij[k,1]]=r[k] end C = cholesky(R).L #todo - #values diverge at this point from SMT code; verify impact - F = ones(nt, 1) #todo - examine if this should be a parameter for this function - Ft = C \ F + F = ones(nt,1) #todo - examine if this should be a parameter for this function + Ft = C\F Q, G = qr(Ft) Q = Array(Q) - Yt = C \ y_norma + Yt = C\y_norma #todo - in smt, they check if the matrix is ill-conditioned using SVD. Verify and include if necessary beta = G \ [(transpose(Q) ⋅ Yt)] rho = Yt .- (Ft .* beta) gamma = transpose(C) \ rho - sigma2 = sum((rho) .^ 2, dims = 1) / nt - detR = prod(diag(C) .^ (2.0 / nt)) - reduced_likelihood_function_value = -nt * log10(sum(sigma2)) - nt * log10(detR) + sigma2 = sum((rho).^2, dims=1)/nt + detR = prod(diag(C).^(2.0/nt)) + reduced_likelihood_function_value = - nt * log10(sum(sigma2)) - nt * log10(detR) return beta, gamma, reduced_likelihood_function_value end @@ -483,34 +464,34 @@ end # https://github.com/scikit-learn/scikit-learn/blob/80598905e/sklearn/cross_decomposition/_pls.py -function _center_scale(X, Y) +function _center_scale(X,Y) x_mean = mean(X, dims = 1) X .-= x_mean y_mean = mean(Y, dims = 1) Y .-= y_mean - x_std = std(X, dims = 1) - x_std[x_std.==0] .= 1.0 + x_std = std(X, dims = 1) + x_std[x_std .== 0] .= 1.0 X ./= x_std - y_std = std(Y, dims = 1) - y_std[y_std.==0] .= 1.0 + y_std = std(Y, dims = 1) + y_std[y_std .==0] .= 1.0 Y ./= y_std - return X, Y + return X,Y end -function _svd_flip_1d(u, v) +function _svd_flip_1d(u,v) # equivalent of https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L149 biggest_abs_val_idx = findmax(abs.(vec(u)))[2] sign_ = sign(u[biggest_abs_val_idx]) u .*= sign_ v .*= sign_ - + end -function _get_first_singular_vectors_power_method(X, Y) +function _get_first_singular_vectors_power_method(X,Y) my_eps = eps() y_score = vec(Y) x_weights = transpose(X)y_score / dot(y_score, y_score) - x_weights ./= (sqrt(dot(x_weights, x_weights)) + my_eps) + x_weights ./= (sqrt(dot(x_weights,x_weights)) + my_eps) x_score = X * x_weights y_weights = transpose(Y)x_score / dot(x_score, x_score) y_score = Y * y_weights / (dot(y_weights, y_weights) + my_eps) @@ -521,28 +502,28 @@ function _get_first_singular_vectors_power_method(X, Y) return x_weights, y_weights end -function _modified_pls(X, Y, n_components) - x_weights_ = zeros(size(X, 2), n_components) - _x_scores = zeros(size(X, 1), n_components) - x_loadings_ = zeros(size(X, 2), n_components) - Xk, Yk = _center_scale(X, Y) +function _modified_pls(X,Y, n_components) + x_weights_ = zeros(size(X,2),n_components) + _x_scores = zeros(size(X,1), n_components) + x_loadings_ = zeros(size(X,2),n_components) + Xk, Yk = _center_scale(X,Y) - for k = 1:n_components - x_weights, y_weights = _get_first_singular_vectors_power_method(Xk, Yk) + for k in 1:n_components + x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) if x_weights == false break - end + end - _svd_flip_1d(x_weights, y_weights) + _svd_flip_1d(x_weights, y_weights) x_scores = Xk * x_weights x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) Xk = Xk - (x_scores * x_loadings) - y_loadings = transpose(x_scores) * Yk / dot(x_scores, x_scores) - Yk = Yk - x_scores * y_loadings - x_weights_[:, k] = x_weights - _x_scores[:, k] = x_scores - x_loadings_[:, k] = vec(x_loadings) + y_loadings = transpose(x_scores)*Yk / dot(x_scores, x_scores) + Yk = Yk - x_scores*y_loadings + x_weights_[:,k] = x_weights + _x_scores[:,k] = x_scores + x_loadings_[:,k] = vec(x_loadings) end x_rotations_ = x_weights_ * pinv(transpose(x_loadings_)x_weights_) diff --git a/src/Surrogates.jl b/src/Surrogates.jl index eaf50ed44..16335f97a 100644 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -31,11 +31,11 @@ function RadialBasisStructure(; radial_function, scale_factor, sparse) end #Kriging structure: -function KrigingStructure(; p, theta) +function KrigingStructure(;p,theta) return (name = "Kriging", p = p, theta = theta) end -function GEKStructure(; p, theta) +function GEKStructure(;p,theta) return (name = "GEK", p = p, theta = theta) end @@ -45,12 +45,12 @@ function LinearStructure() end #InverseDistance structure -function InverseDistanceStructure(; p) +function InverseDistanceStructure(;p) return (name = "InverseDistanceSurrogate", p = p) end #Lobachevsky structure -function LobachevskyStructure(; alpha, n, sparse) +function LobachevskyStructure(;alpha,n,sparse) return (name = "LobachevskySurrogate", alpha = alpha, n = n, sparse = sparse) end @@ -61,7 +61,7 @@ function NeuralStructure(; model, loss, opt, n_echos) end #Random forest structure -function RandomForestStructure(; num_round) +function RandomForestStructure(;num_round) return (name = "RandomForestSurrogate", num_round = num_round) end @@ -81,7 +81,7 @@ function PolyChaosStructure(; op) end export current_surrogates -export GEKPLS +export GEKPLS export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure @@ -89,7 +89,7 @@ export WendlandStructure export AbstractSurrogate, SamplingAlgorithm export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point # radial basis functions -export linearRadial, cubicRadial, multiquadricRadial, thinplateRadial +export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial # samplers export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, @@ -97,7 +97,7 @@ export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, export RandomSample, KroneckerSample, GoldenSample, SectionSample # Optimization algorithms -export SRBF, LCBS, EI, DYCORS, SOP, EGO, RTEA, SMB, surrogate_optimize +export SRBF,LCBS,EI,DYCORS,SOP,EGO,RTEA,SMB,surrogate_optimize export LobachevskySurrogate, lobachevsky_integral, lobachevsky_integrate_dimension export LinearSurrogate export SVMSurrogate diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index fcda8b58b..28c070989 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -7,9 +7,9 @@ function vector_of_tuples_to_matrix(v) num_rows = length(v) num_cols = length(first(v)) K = zeros(num_rows, num_cols) - for row = 1:num_rows - for col = 1:num_cols - K[row, col] = v[row][col] + for row in 1:num_rows + for col in 1:num_cols + K[row, col]=v[row][col] end end return K @@ -20,8 +20,8 @@ function vector_of_tuples_to_matrix2(v) num_rows = length(v) num_cols = length(first(first(v))) K = zeros(num_rows, num_cols) - for row = 1:num_rows - for col = 1:num_cols + for row in 1:num_rows + for col in 1:num_cols K[row, col] = v[row][1][col] end end @@ -38,55 +38,55 @@ function water_flow(x) H_l = x[6] L = x[7] K_w = x[8] - log_val = log(r / r_w) - return (2 * pi * T_u * (H_u - H_l)) / - (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l)) + log_val = log(r/r_w) + return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) end n = 1000 -lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] -ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] -x = sample(n, lb, ub, SobolSample()) +d = 8 +lb = [0.05,100,63070,990,63.1,700,1120,9855] +ub = [0.15,50000,115600,1110,116,820,1680,12045] +x = sample(n,lb,ub,SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x)) -y = reshape(water_flow.(x), (size(x, 1), 1)) +y = reshape(water_flow.(x),(size(x,1),1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test, lb, ub, GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = water_flow.(x_test) -@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin +@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 0.03, atol = 0.02) #rmse: 0.039 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 0.03, atol=0.02) #rmse: 0.039 end -@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin +@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 end -@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin +@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin n_comp = 3 delta_x = 0.0001 extra_points = 3 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 end ## welded beam tests @@ -94,77 +94,76 @@ function welded_beam(x) h = x[1] l = x[2] t = x[3] - a = 6000 / (sqrt(2) * h * l) - b = - (6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))) / - (2 * (0.707 * h * l * (l^2 / 12 + 0.25 * (h + t)^2))) - return (sqrt(a^2 + b^2 + l * a * b)) / (sqrt(0.25 * (l^2 + (h + t)^2))) + a = 6000/(sqrt(2)*h*l) + b = (6000*(14+0.5*l)*sqrt(0.25*(l^2+(h+t)^2)))/(2*(0.707*h*l*(l^2/12 + 0.25*(h+t)^2))) + return (sqrt(a^2+b^2 + l*a*b))/(sqrt(0.25*(l^2+(h+t)^2))) end n = 1000 -lb = [0.125, 5.0, 5.0] -ub = [1.0, 10.0, 10.0] -x = sample(n, lb, ub, SobolSample()) +d = 3 +lb = [0.125,5.0,5.0] +ub = [1.,10.,10.] +x = sample(n,lb,ub,SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(welded_beam, x)) -y = reshape(welded_beam.(x), (size(x, 1), 1)) +y = reshape(welded_beam.(x),(size(x,1),1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test, lb, ub, GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = welded_beam.(x_test) -@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin +@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 39.0, atol = 0.5) #rmse: 38.988 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 39.0, atol=0.5) #rmse: 38.988 end -@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin +@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 39.5, atol = 0.5) #rmse: 39.481 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 39.5, atol=0.5) #rmse: 39.481 end ## increasing extra points increases accuracy -@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin +@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin n_comp = 2 delta_x = 0.0001 extra_points = 4 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 37.5, atol = 0.5) #rmse: 37.87 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 37.5, atol=0.5) #rmse: 37.87 end ## sphere function tests function sphere_function(x) - return sum(x .^ 2) + return sum(x.^2) end ## 3D n = 100 lb = [-5.0, -5.0, -5.0] -ub = [5.0, 5.0, 5.0] -x = sample(n, lb, ub, SobolSample()) +ub = [5.0, 5.0 ,5.0] +x = sample(n,lb,ub,SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x), (size(x, 1), 1)) +y = reshape(sphere_function.(x),(size(x,1),1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test, lb, ub, GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 7: Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin @@ -172,24 +171,25 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 0.001, atol = 0.05) #rmse: 0.00083 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 0.001, atol=0.05) #rmse: 0.00083 end ## 2D -n = 50 +n = 50 +d = 2 lb = [-10.0, -10.0] ub = [10.0, 10.0] -x = sample(n, lb, ub, SobolSample()) +x = sample(n,lb,ub,SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x), (size(x, 1), 1)) +y = reshape(sphere_function.(x),(size(x,1),1)) xlimits = hcat(lb, ub) -n_test = 10 -x_test = sample(n_test, lb, ub, GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 10 +x_test = sample(n_test,lb,ub,GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 8: Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2" begin @@ -197,16 +197,16 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) - @test isapprox(rmse, 0.1, atol = 0.5) #rmse: 0.0022 + rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) + @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 end @testset "Test 9: Add Point Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin #first we create a surrogate model with just 3 input points initial_x_vec = [(1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0)] - initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec, 1), 1)) + initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec,1),1)) initial_X = vector_of_tuples_to_matrix(initial_x_vec) initial_grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, initial_x_vec)) lb = [-5.0, -5.0, -5.0] @@ -216,31 +216,22 @@ end delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS( - initial_X, - initial_y, - initial_grads, - n_comp, - delta_x, - xlimits, - extra_points, - initial_theta, - ) - n_test = 100 - x_test = sample(n_test, lb, ub, GoldenSample()) - X_test = vector_of_tuples_to_matrix(x_test) - y_true = sphere_function.(x_test) + g = GEKPLS(initial_X, initial_y, initial_grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + n_test = 100 + x_test = sample(n_test,lb,ub,GoldenSample()) + X_test = vector_of_tuples_to_matrix(x_test) + y_true = sphere_function.(x_test) y_pred1 = g(X_test) - rmse1 = sqrt(sum(((y_pred1 - y_true) .^ 2) / n_test)) #rmse1 = 31.91 + rmse1 = sqrt(sum(((y_pred1 - y_true).^2)/n_test)) #rmse1 = 31.91 #then we update the model with more points to see if performance improves n = 100 - x = sample(n, lb, ub, SobolSample()) + x = sample(n,lb,ub,SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) - y = reshape(sphere_function.(x), (size(x, 1), 1)) + y = reshape(sphere_function.(x),(size(x,1),1)) add_point!(g, X, y, grads) y_pred2 = g(X_test) - rmse2 = sqrt(sum(((y_pred2 - y_true) .^ 2) / n_test)) #rmse2 = 0.0015 + rmse2 = sqrt(sum(((y_pred2 - y_true).^2)/n_test)) #rmse2 = 0.0015 @test (rmse2 < rmse1) end diff --git a/test/runtests.jl b/test/runtests.jl index 7ff07219f..aa703cca4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,8 +4,8 @@ using SafeTestsets using Pkg function dev_subpkg(subpkg) - subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) - Pkg.develop(PackageSpec(path = subpkg_path)) + subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) + Pkg.develop(PackageSpec(path=subpkg_path)) end for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", "SurrogatesRandomForest", "SurrogatesSVM"] @@ -15,33 +15,15 @@ for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", end end -@time @safetestset "GEKPLS.jl" begin - include("GEKPLS.jl") -end -@time @safetestset "Radials.jl" begin - include("Radials.jl") -end -@time @safetestset "Kriging.jl" begin - include("Kriging.jl") -end -@time @safetestset "Sampling" begin - include("sampling.jl") -end -@time @safetestset "Optimization" begin - include("optimization.jl") -end -@time @safetestset "LinearSurrogate" begin - include("linearSurrogate.jl") -end -@time @safetestset "Lobachevsky" begin - include("lobachevsky.jl") -end -@time @safetestset "InverseDistanceSurrogate" begin - include("inverseDistanceSurrogate.jl") -end -@time @safetestset "SecondOrderPolynomialSurrogate" begin - include("secondOrderPolynomialSurrogate.jl") -end +@time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end +@time @safetestset "Radials.jl" begin include("Radials.jl") end +@time @safetestset "Kriging.jl" begin include("Kriging.jl") end +@time @safetestset "Sampling" begin include("sampling.jl") end +@time @safetestset "Optimization" begin include("optimization.jl") end +@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end +@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end +@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end +@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end # @time @safetestset "AD_Compatibility" begin include("AD_compatibility.jl") end @time @safetestset "Wendland" begin include("Wendland.jl") end @time @safetestset "VariableFidelity" begin include("VariableFidelity.jl") end From e39d42140c783df5f9e6e7a46e9ec73e58f437e2 Mon Sep 17 00:00:00 2001 From: Vikram Date: Tue, 28 Jun 2022 19:04:30 +0530 Subject: [PATCH 11/12] update code after rebase with master --- Project.toml | 1 + docs/pages.jl | 78 ++++++----- src/GEKPLS.jl | 336 +++++++++++++++++++++++----------------------- src/Surrogates.jl | 16 +-- test/GEKPLS.jl | 148 ++++++++++---------- test/runtests.jl | 22 +-- 6 files changed, 303 insertions(+), 298 deletions(-) diff --git a/Project.toml b/Project.toml index d87569e97..9893178ad 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/docs/pages.jl b/docs/pages.jl index 5101a1e11..b53c2d048 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,40 +1,38 @@ -pages = [ - "index.md" - "Tutorials" => [ - "Basics" => "tutorials.md", - "Radials" => "radials.md", - "Kriging" => "kriging.md", - "Gaussian Process" =>"abstractgps.md", - "Lobachevsky" => "lobachevsky.md", - "Linear" => "LinearSurrogate.md", - "InverseDistance" => "InverseDistance.md", - "RandomForest" => "randomforest.md", - "SecondOrderPolynomial" => "secondorderpoly.md", - "NeuralSurrogate" => "neural.md", - "Wendland" => "wendland.md", - "Polynomial Chaos" => "polychaos.md", - "Variable Fidelity" => "variablefidelity.md", - "Gradient Enhanced Kriging" => "gek.md", - "GEKPLS" => "gekpls.md" - ] - "User guide" => [ - "Samples" => "samples.md", - "Surrogates" => "surrogate.md", - "Optimization" => "optimizations.md" - ] - "Benchmarks" => [ - "Sphere function" => "sphere_function.md", - "Lp norm" => "lp.md", - "Rosenbrock" => "rosenbrock.md", - "Tensor product" => "tensor_prod.md", - "Cantilever beam" => "cantilever.md", - "Water Flow function" => "water_flow.md", - "Welded beam function" => "welded_beam.md", - "Branin function" => "BraninFunction.md", - "Ackley function" => "ackley.md", - "Gramacy & Lee Function" => "gramacylee.md", - "Salustowicz Benchmark" => "Salustowicz.md", - "Multi objective optimization" => "multi_objective_opt.md" - ] - "Contributing" => "contributing.md" - ] +pages = ["index.md" + "Tutorials" => [ + "Basics" => "tutorials.md", + "Radials" => "radials.md", + "Kriging" => "kriging.md", + "Gaussian Process" => "abstractgps.md", + "Lobachevsky" => "lobachevsky.md", + "Linear" => "LinearSurrogate.md", + "InverseDistance" => "InverseDistance.md", + "RandomForest" => "randomforest.md", + "SecondOrderPolynomial" => "secondorderpoly.md", + "NeuralSurrogate" => "neural.md", + "Wendland" => "wendland.md", + "Polynomial Chaos" => "polychaos.md", + "Variable Fidelity" => "variablefidelity.md", + "Gradient Enhanced Kriging" => "gek.md", + "GEKPLS" => "gekpls.md", + ] + "User guide" => [ + "Samples" => "samples.md", + "Surrogates" => "surrogate.md", + "Optimization" => "optimizations.md", + ] + "Benchmarks" => [ + "Sphere function" => "sphere_function.md", + "Lp norm" => "lp.md", + "Rosenbrock" => "rosenbrock.md", + "Tensor product" => "tensor_prod.md", + "Cantilever beam" => "cantilever.md", + "Water Flow function" => "water_flow.md", + "Welded beam function" => "welded_beam.md", + "Branin function" => "BraninFunction.md", + "Ackley function" => "ackley.md", + "Gramacy & Lee Function" => "gramacylee.md", + "Salustowicz Benchmark" => "Salustowicz.md", + "Multi objective optimization" => "multi_objective_opt.md", + ] + "Contributing" => "contributing.md"] diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index c228dea9b..255cbc588 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -1,18 +1,18 @@ using LinearAlgebra using Statistics -mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate +mutable struct GEKPLS{T <: AbstractFloat} <: AbstractSurrogate x::Matrix{T} #1 y::Matrix{T} #2 grads::Matrix{T} #3 xl::Matrix{T} #xlimits #4 - delta:: T #5 + delta::T #5 extra_points::Int #6 num_components::Int #7 beta::Vector{T} #8 gamma::Matrix{T} #9 theta::Vector{T} #10 - reduced_likelihood_function_value:: T #11 + reduced_likelihood_function_value::T #11 X_offset::Matrix{T} #12 X_scale::Matrix{T} #13 X_after_std::Matrix{T} #14 - X after standardization @@ -22,11 +22,11 @@ mutable struct GEKPLS{T<:AbstractFloat} <: AbstractSurrogate end function bounds_error(x, xl) - num_x_rows = size(x,1) + num_x_rows = size(x, 1) num_dim = size(xl, 1) for i in 1:num_x_rows for j in 1:num_dim - if (x[i, j] < xl[j,1] || x[i,j] > xl[j,2]) + if (x[i, j] < xl[j, 1] || x[i, j] > xl[j, 2]) return true end end @@ -36,7 +36,7 @@ end #constructor for GEKPLS Struct function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) - + #ensure that X values are within the upper and lower bounds if bounds_error(X, xlimits) println("X values outside bounds") @@ -44,34 +44,38 @@ function GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, θ) end theta = [θ for i in 1:n_comp] - pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) - X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, y_after_PLS) + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, + xlimits, extra_points) + X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, + y_after_PLS) D, ij = cross_distances(X_after_std) - pls_mean_reshaped = reshape(pls_mean, (size(X,2),n_comp)) + pls_mean_reshaped = reshape(pls_mean, (size(X, 2), n_comp)) d = componentwise_distance_PLS(D, "squar_exp", n_comp, pls_mean_reshaped) nt, nd = size(X_after_PLS) - beta, gamma, reduced_likelihood_function_value = _reduced_likelihood_function(theta, "squar_exp", d, nt, ij, y_after_std) - return GEKPLS(X, y, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, reduced_likelihood_function_value, - X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) + beta, gamma, reduced_likelihood_function_value = _reduced_likelihood_function(theta, + "squar_exp", + d, nt, ij, + y_after_std) + return GEKPLS(X, y, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, + reduced_likelihood_function_value, + X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) println("struct created") end - # predictor function (g::GEKPLS)(X_test) n_eval, n_features_x = size(X_test) X_cont = (X_test .- g.X_offset) ./ g.X_scale dx = differences(X_cont, g.X_after_std) pred_d = componentwise_distance_PLS(dx, "squar_exp", g.num_components, g.pls_mean) - nt = size(g.X_after_std,1) - r = transpose(reshape(squar_exp(g.theta, pred_d),(nt,n_eval))) - f = ones(n_eval,1) + nt = size(g.X_after_std, 1) + r = transpose(reshape(squar_exp(g.theta, pred_d), (nt, n_eval))) + f = ones(n_eval, 1) y_ = (f * g.beta) + (r * g.gamma) y = g.y_mean .+ g.y_std * y_ return y end - function add_point!(g::GEKPLS, new_x, new_y, new_grads) if new_x in g.x println("Adding a sample that already exists. Cannot build GEKPLS") @@ -80,20 +84,27 @@ function add_point!(g::GEKPLS, new_x, new_y, new_grads) if bounds_error(new_x, g.xl) println("x values outside bounds") - return + return end - g.x = vcat(g.x, new_x) g.y = vcat(g.y, new_y) g.grads = vcat(g.grads, new_grads) - pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, g.grads, g.delta, g.xl, g.extra_points) - g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, y_after_PLS) + pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(g.x, g.y, g.num_components, + g.grads, g.delta, g.xl, + g.extra_points) + g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, + y_after_PLS) D, ij = cross_distances(g.X_after_std) g.pls_mean = reshape(pls_mean, (size(g.x, 2), g.num_components)) d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) nt, nd = size(X_after_PLS) - g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, "squar_exp", d, nt, ij, y_after_std) + g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, + "squar_exp", + d, + nt, + ij, + y_after_std) end function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) @@ -119,45 +130,43 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) # and https://github.com/SMTorg/smt/blob/f124c01ffa78c04b80221dded278a20123dac742/smt/surrogate_models/gekpls.py#L48 nt, dim = size(X) - XX = zeros(0,dim) - yy = zeros(0,size(y)[2]) + XX = zeros(0, dim) + yy = zeros(0, size(y)[2]) coeff_pls = zeros((dim, n_comp, nt)) - + for i in 1:nt - if dim >= 3 - bb_vals = circshift(boxbehnken(dim, 1),1) + if dim >= 3 + bb_vals = circshift(boxbehnken(dim, 1), 1) else - bb_vals = [ - 0.0 0.0; #center - 1.0 0.0; #right - 0.0 1.0; #up - -1.0 0.0; #left - 0.0 -1.0; #down - 1.0 1.0; #right up - -1.0 1.0; #left up - -1.0 -1.0; #left down - 1.0 -1.0; #right down - ] + bb_vals = [0.0 0.0; #center + 1.0 0.0; #right + 0.0 1.0; #up + -1.0 0.0; #left + 0.0 -1.0; #down + 1.0 1.0; #right up + -1.0 1.0; #left up + -1.0 -1.0; #left down + 1.0 -1.0] end - _X = zeros((size(bb_vals)[1], dim)) - _y = zeros((size(bb_vals)[1], 1)) + _X = zeros((size(bb_vals)[1], dim)) + _y = zeros((size(bb_vals)[1], 1)) bb_vals = bb_vals .* (delta_x * (xlimits[:, 2] - xlimits[:, 1]))' #smt calls this sign. I've called it bb_vals - _X = X[i, :]' .+ bb_vals - bb_vals = bb_vals .* grads[i,:]' - _y = y[i, :] .+ sum(bb_vals, dims=2) - + _X = X[i, :]' .+ bb_vals + bb_vals = bb_vals .* grads[i, :]' + _y = y[i, :] .+ sum(bb_vals, dims = 2) + #_pls.fit(_X, _y) # relic from sklearn versiom; retained for future reference. #coeff_pls[:, :, i] = _pls.x_rotations_ #relic from sklearn versiom; retained for future reference. coeff_pls[:, :, i] = _modified_pls(_X, _y, n_comp) #_modified_pls returns the equivalent of SKLearn's _pls.x_rotations_ if extra_points != 0 - start_index = max(1, length(coeff_pls[:,1,i])-extra_points+1) - max_coeff = sortperm(broadcast(abs,coeff_pls[:,1,i]))[start_index:end] + start_index = max(1, length(coeff_pls[:, 1, i]) - extra_points + 1) + max_coeff = sortperm(broadcast(abs, coeff_pls[:, 1, i]))[start_index:end] for ii in max_coeff XX = [XX; transpose(X[i, :])] - XX[end, ii] += delta_x * (xlimits[ii,2]-xlimits[ii,1]) + XX[end, ii] += delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) yy = [yy; y[i]] - yy[end] += grads[i,ii] * delta_x * (xlimits[ii,2]-xlimits[ii,1]) + yy[end] += grads[i, ii] * delta_x * (xlimits[ii, 2] - xlimits[ii, 1]) end end end @@ -166,89 +175,89 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) y = [y; yy] end - pls_mean = mean(broadcast(abs,coeff_pls),dims=3) + pls_mean = mean(broadcast(abs, coeff_pls), dims = 3) return pls_mean, X, y end ######start of bbdesign###### +# +# Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' +# https://github.com/phrb/ExperimentalDesign.jl - # - # Adapted from 'ExperimentalDesign.jl: Design of Experiments in Julia' - # https://github.com/phrb/ExperimentalDesign.jl +# MIT License - # MIT License +# ExperimentalDesign.jl: Design of Experiments in Julia +# Copyright (C) 2019 Pedro Bruel - # ExperimentalDesign.jl: Design of Experiments in Julia - # Copyright (C) 2019 Pedro Bruel +# Permission is hereby granted, free of charge, to any person obtaining a copy of +# this software and associated documentation files (the "Software"), to deal in +# the Software without restriction, including without limitation the rights to +# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +# the Software, and to permit persons to whom the Software is furnished to do so, +# subject to the following conditions: - # Permission is hereby granted, free of charge, to any person obtaining a copy of - # this software and associated documentation files (the "Software"), to deal in - # the Software without restriction, including without limitation the rights to - # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of - # the Software, and to permit persons to whom the Software is furnished to do so, - # subject to the following conditions: +# The above copyright notice and this permission notice (including the next +# paragraph) shall be included in all copies or substantial portions of the +# Software. - # The above copyright notice and this permission notice (including the next - # paragraph) shall be included in all copies or substantial portions of the - # Software. +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +# 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. +# - # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS - # FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR - # 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) +end - function boxbehnken(matrix_size::Int) - boxbehnken(matrix_size, matrix_size) - end - - function boxbehnken(matrix_size::Int, center::Int) - @assert matrix_size>=3 - - A_fact = explicit_fullfactorial(Tuple([-1,1] for i = 1:2)) - - rows = floor(Int, (0.5*matrix_size*(matrix_size-1))*size(A_fact)[1]) - - A = zeros(rows, matrix_size) - - l = 0 - for i in 1:matrix_size-1 - for j in i+1:matrix_size - l = l +1 - A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], i] = A_fact[:, 1] - A[max(0, (l - 1)*size(A_fact)[1])+1:l*size(A_fact)[1], j] = A_fact[:, 2] - end - end - - if center == matrix_size - if matrix_size <= 16 - points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] - center = points[matrix_size] - end +function boxbehnken(matrix_size::Int, center::Int) + @assert matrix_size >= 3 + + A_fact = explicit_fullfactorial(Tuple([-1, 1] for i in 1:2)) + + rows = floor(Int, (0.5 * matrix_size * (matrix_size - 1)) * size(A_fact)[1]) + + A = zeros(rows, matrix_size) + + l = 0 + for i in 1:(matrix_size - 1) + for j in (i + 1):matrix_size + l = l + 1 + A[(max(0, (l - 1) * size(A_fact)[1]) + 1):(l * size(A_fact)[1]), i] = A_fact[:, + 1] + A[(max(0, (l - 1) * size(A_fact)[1]) + 1):(l * size(A_fact)[1]), j] = A_fact[:, + 2] end - - A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) - end - - function explicit_fullfactorial(factors::Tuple) - explicit_fullfactorial(fullfactorial(factors)) - end - - function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) - hcat(vcat.(collect(iterator)...)...) end - - function fullfactorial(factors::Tuple) - Base.Iterators.product(factors...) + + if center == matrix_size + if matrix_size <= 16 + points = [0, 0, 3, 3, 6, 6, 6, 8, 9, 10, 12, 12, 13, 14, 15, 16] + center = points[matrix_size] + end end - + + A = transpose(hcat(transpose(A), transpose(zeros(center, matrix_size)))) +end + +function explicit_fullfactorial(factors::Tuple) + explicit_fullfactorial(fullfactorial(factors)) +end + +function explicit_fullfactorial(iterator::Base.Iterators.ProductIterator) + hcat(vcat.(collect(iterator)...)...) +end + +function fullfactorial(factors::Tuple) + Base.Iterators.product(factors...) +end ######end of bb design###### -function standardization(X,y) +function standardization(X, y) """ We substract the mean from each variable. Then, we divide the values of each variable by its standard deviation. @@ -280,14 +289,13 @@ function standardization(X,y) #Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L21 X_offset = mean(X, dims = 1) X_scale = std(X, dims = 1) - X_scale = map(x -> (x==0.0) ? x=1 : x=x, X_scale) #to prevent division by 0 below + X_scale = map(x -> (x == 0.0) ? x = 1 : x = x, X_scale) #to prevent division by 0 below y_mean = mean(y) y_std = std(y) - y_std = map(y -> (y==0) ? y=1 : y=y, y_std) #to prevent division by 0 below - X = (X.-X_offset) ./ X_scale + y_std = map(y -> (y == 0) ? y = 1 : y = y, y_std) #to prevent division by 0 below + X = (X .- X_offset) ./ X_scale y = (y .- y_mean) ./ y_std return X, y, X_offset, y_mean, X_scale, y_std - end function cross_distances(X) @@ -311,18 +319,17 @@ function cross_distances(X) """ # equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L86 n_samples, n_features = size(X) - n_nonzero_cross_dist = ( n_samples * (n_samples - 1) ) ÷ 2 + n_nonzero_cross_dist = (n_samples * (n_samples - 1)) ÷ 2 ij = zeros((n_nonzero_cross_dist, 2)) D = zeros((n_nonzero_cross_dist, n_features)) ll_1 = 0 - - for k in 1:n_samples - 1 + + for k in 1:(n_samples - 1) ll_0 = ll_1 + 1 ll_1 = ll_0 + n_samples - k - 1 ij[ll_0:ll_1, 1] .= k - ij[ll_0:ll_1, 2] = k+1:1:n_samples - D[ll_0:ll_1, :] = -(X[(k + 1) : n_samples,:] .- X[k,:]') - + ij[ll_0:ll_1, 2] = (k + 1):1:n_samples + D[ll_0:ll_1, :] = -(X[(k + 1):n_samples, :] .- X[k, :]') end return D, Int.(ij) end @@ -367,7 +374,7 @@ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) D_corr = zeros((size(D)[1], n_comp)) if corr == "squar_exp" - D_corr = D.^2 * coeff_pls.^2 + D_corr = D .^ 2 * coeff_pls .^ 2 else #abs_exp D_corr = abs.(D) * abs.(coeff_pls) end @@ -387,27 +394,27 @@ function squar_exp(theta, d) Returns: -------- r: array containing the values of the autocorrelation model - + """ n_components = size(d)[2] - theta = reshape(theta, (1,n_components)) - return exp.(-sum(theta .* d, dims=2)) + theta = reshape(theta, (1, n_components)) + return exp.(-sum(theta .* d, dims = 2)) end function differences(X, Y) #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L392 #code credit: Elias Carvalho - https://stackoverflow.com/questions/72392010/row-wise-operations-between-matrices-in-julia - Rx = repeat(X, inner=(size(Y, 1), 1)) + Rx = repeat(X, inner = (size(Y, 1), 1)) Ry = repeat(Y, size(X, 1)) return Rx - Ry end -function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise=0.0) +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 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. - + Parameters ---------- theta: array containing the parameters at which the Gaussian Process model parameters should be determined. @@ -430,68 +437,65 @@ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, no #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 reduced_likelihood_function_value = -Inf nugget = 1000000.0 * eps() #a jitter for numerical stability; reducing the multiple from 1000000.0 results in positive definite error for Cholesky decomposition; - if kernel_type =="squar_exp" #todo - add other kernel type abs_exp etc. - r = squar_exp(theta,d) + if kernel_type == "squar_exp" #todo - add other kernel type abs_exp etc. + r = squar_exp(theta, d) end - R = (I + zeros(nt,nt)) .* (1.0 + nugget + noise) + R = (I + zeros(nt, nt)) .* (1.0 + nugget + noise) - for k in 1:size(ij)[1] - R[ij[k,1],ij[k,2]]=r[k] - R[ij[k,2],ij[k,1]]=r[k] + for k in 1:size(ij)[1] + R[ij[k, 1], ij[k, 2]] = r[k] + R[ij[k, 2], ij[k, 1]] = r[k] end C = cholesky(R).L #todo - #values diverge at this point from SMT code; verify impact - F = ones(nt,1) #todo - examine if this should be a parameter for this function - Ft = C\F + F = ones(nt, 1) #todo - examine if this should be a parameter for this function + Ft = C \ F Q, G = qr(Ft) Q = Array(Q) - Yt = C\y_norma + Yt = C \ y_norma #todo - in smt, they check if the matrix is ill-conditioned using SVD. Verify and include if necessary beta = G \ [(transpose(Q) ⋅ Yt)] rho = Yt .- (Ft .* beta) gamma = transpose(C) \ rho - sigma2 = sum((rho).^2, dims=1)/nt - detR = prod(diag(C).^(2.0/nt)) - reduced_likelihood_function_value = - nt * log10(sum(sigma2)) - nt * log10(detR) + sigma2 = sum((rho) .^ 2, dims = 1) / nt + detR = prod(diag(C) .^ (2.0 / nt)) + reduced_likelihood_function_value = -nt * log10(sum(sigma2)) - nt * log10(detR) return beta, gamma, reduced_likelihood_function_value end - ### MODIFIED PLS BELOW ### # 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 - -function _center_scale(X,Y) +function _center_scale(X, Y) x_mean = mean(X, dims = 1) X .-= x_mean y_mean = mean(Y, dims = 1) Y .-= y_mean - x_std = std(X, dims = 1) + x_std = std(X, dims = 1) x_std[x_std .== 0] .= 1.0 X ./= x_std - y_std = std(Y, dims = 1) - y_std[y_std .==0] .= 1.0 + y_std = std(Y, dims = 1) + y_std[y_std .== 0] .= 1.0 Y ./= y_std - return X,Y + return X, Y end -function _svd_flip_1d(u,v) +function _svd_flip_1d(u, v) # equivalent of https://github.com/scikit-learn/scikit-learn/blob/80598905e517759b4696c74ecc35c6e2eb508cff/sklearn/cross_decomposition/_pls.py#L149 biggest_abs_val_idx = findmax(abs.(vec(u)))[2] sign_ = sign(u[biggest_abs_val_idx]) u .*= sign_ v .*= sign_ - end -function _get_first_singular_vectors_power_method(X,Y) +function _get_first_singular_vectors_power_method(X, Y) my_eps = eps() y_score = vec(Y) x_weights = transpose(X)y_score / dot(y_score, y_score) - x_weights ./= (sqrt(dot(x_weights,x_weights)) + my_eps) + x_weights ./= (sqrt(dot(x_weights, x_weights)) + my_eps) x_score = X * x_weights y_weights = transpose(Y)x_score / dot(x_score, x_score) y_score = Y * y_weights / (dot(y_weights, y_weights) + my_eps) @@ -502,28 +506,28 @@ function _get_first_singular_vectors_power_method(X,Y) return x_weights, y_weights end -function _modified_pls(X,Y, n_components) - x_weights_ = zeros(size(X,2),n_components) - _x_scores = zeros(size(X,1), n_components) - x_loadings_ = zeros(size(X,2),n_components) - Xk, Yk = _center_scale(X,Y) +function _modified_pls(X, Y, n_components) + x_weights_ = zeros(size(X, 2), n_components) + _x_scores = zeros(size(X, 1), n_components) + x_loadings_ = zeros(size(X, 2), n_components) + Xk, Yk = _center_scale(X, Y) for k in 1:n_components - x_weights, y_weights = _get_first_singular_vectors_power_method(Xk,Yk) + x_weights, y_weights = _get_first_singular_vectors_power_method(Xk, Yk) if x_weights == false break - end + end - _svd_flip_1d(x_weights, y_weights) + _svd_flip_1d(x_weights, y_weights) x_scores = Xk * x_weights x_loadings = transpose(x_scores)Xk / dot(x_scores, x_scores) Xk = Xk - (x_scores * x_loadings) - y_loadings = transpose(x_scores)*Yk / dot(x_scores, x_scores) - Yk = Yk - x_scores*y_loadings - x_weights_[:,k] = x_weights - _x_scores[:,k] = x_scores - x_loadings_[:,k] = vec(x_loadings) + y_loadings = transpose(x_scores) * Yk / dot(x_scores, x_scores) + Yk = Yk - x_scores * y_loadings + x_weights_[:, k] = x_weights + _x_scores[:, k] = x_scores + x_loadings_[:, k] = vec(x_loadings) end x_rotations_ = x_weights_ * pinv(transpose(x_loadings_)x_weights_) diff --git a/src/Surrogates.jl b/src/Surrogates.jl index 16335f97a..eaf50ed44 100644 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -31,11 +31,11 @@ function RadialBasisStructure(; radial_function, scale_factor, sparse) end #Kriging structure: -function KrigingStructure(;p,theta) +function KrigingStructure(; p, theta) return (name = "Kriging", p = p, theta = theta) end -function GEKStructure(;p,theta) +function GEKStructure(; p, theta) return (name = "GEK", p = p, theta = theta) end @@ -45,12 +45,12 @@ function LinearStructure() end #InverseDistance structure -function InverseDistanceStructure(;p) +function InverseDistanceStructure(; p) return (name = "InverseDistanceSurrogate", p = p) end #Lobachevsky structure -function LobachevskyStructure(;alpha,n,sparse) +function LobachevskyStructure(; alpha, n, sparse) return (name = "LobachevskySurrogate", alpha = alpha, n = n, sparse = sparse) end @@ -61,7 +61,7 @@ function NeuralStructure(; model, loss, opt, n_echos) end #Random forest structure -function RandomForestStructure(;num_round) +function RandomForestStructure(; num_round) return (name = "RandomForestSurrogate", num_round = num_round) end @@ -81,7 +81,7 @@ function PolyChaosStructure(; op) end export current_surrogates -export GEKPLS +export GEKPLS export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure @@ -89,7 +89,7 @@ export WendlandStructure export AbstractSurrogate, SamplingAlgorithm export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point # radial basis functions -export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial +export linearRadial, cubicRadial, multiquadricRadial, thinplateRadial # samplers export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, @@ -97,7 +97,7 @@ export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, export RandomSample, KroneckerSample, GoldenSample, SectionSample # Optimization algorithms -export SRBF,LCBS,EI,DYCORS,SOP,EGO,RTEA,SMB,surrogate_optimize +export SRBF, LCBS, EI, DYCORS, SOP, EGO, RTEA, SMB, surrogate_optimize export LobachevskySurrogate, lobachevsky_integral, lobachevsky_integrate_dimension export LinearSurrogate export SVMSurrogate diff --git a/test/GEKPLS.jl b/test/GEKPLS.jl index 28c070989..26f33f39b 100644 --- a/test/GEKPLS.jl +++ b/test/GEKPLS.jl @@ -1,7 +1,6 @@ using Surrogates using Zygote - function vector_of_tuples_to_matrix(v) #convert training data generated by surrogate sampling into a matrix suitable for GEKPLS num_rows = length(v) @@ -9,7 +8,7 @@ function vector_of_tuples_to_matrix(v) K = zeros(num_rows, num_cols) for row in 1:num_rows for col in 1:num_cols - K[row, col]=v[row][col] + K[row, col] = v[row][col] end end return K @@ -38,55 +37,56 @@ function water_flow(x) H_l = x[6] L = x[7] K_w = x[8] - log_val = log(r/r_w) - return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) + log_val = log(r / r_w) + return (2 * pi * T_u * (H_u - H_l)) / + (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l)) end n = 1000 d = 8 -lb = [0.05,100,63070,990,63.1,700,1120,9855] -ub = [0.15,50000,115600,1110,116,820,1680,12045] -x = sample(n,lb,ub,SobolSample()) +lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] +ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x)) -y = reshape(water_flow.(x),(size(x,1),1)) +y = reshape(water_flow.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = water_flow.(x_test) -@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin +@testset "Test 1: Water Flow Function Test (dimensions = 8; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.03, atol=0.02) #rmse: 0.039 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.03, atol = 0.02) #rmse: 0.039 end -@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin +@testset "Test 2: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) #change hard-coded 2 param to variable y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 end -@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin +@testset "Test 3: Water Flow Function Test (dimensions = 8; n_comp = 3; extra_points = 3)" begin n_comp = 3 delta_x = 0.0001 extra_points = 3 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.02, atol=0.01) #rmse: 0.027 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.02, atol = 0.01) #rmse: 0.027 end ## welded beam tests @@ -94,76 +94,77 @@ function welded_beam(x) h = x[1] l = x[2] t = x[3] - a = 6000/(sqrt(2)*h*l) - b = (6000*(14+0.5*l)*sqrt(0.25*(l^2+(h+t)^2)))/(2*(0.707*h*l*(l^2/12 + 0.25*(h+t)^2))) - return (sqrt(a^2+b^2 + l*a*b))/(sqrt(0.25*(l^2+(h+t)^2))) + a = 6000 / (sqrt(2) * h * l) + b = (6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))) / + (2 * (0.707 * h * l * (l^2 / 12 + 0.25 * (h + t)^2))) + return (sqrt(a^2 + b^2 + l * a * b)) / (sqrt(0.25 * (l^2 + (h + t)^2))) end n = 1000 d = 3 -lb = [0.125,5.0,5.0] -ub = [1.,10.,10.] -x = sample(n,lb,ub,SobolSample()) +lb = [0.125, 5.0, 5.0] +ub = [1.0, 10.0, 10.0] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(welded_beam, x)) -y = reshape(welded_beam.(x),(size(x,1),1)) +y = reshape(welded_beam.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = welded_beam.(x_test) -@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin +@testset "Test 4: Welded Beam Function Test (dimensions = 3; n_comp = 3; extra_points = 2)" begin n_comp = 3 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 39.0, atol=0.5) #rmse: 38.988 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 39.0, atol = 0.5) #rmse: 38.988 end -@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin +@testset "Test 5: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin n_comp = 2 delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 39.5, atol=0.5) #rmse: 39.481 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 39.5, atol = 0.5) #rmse: 39.481 end ## increasing extra points increases accuracy -@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin +@testset "Test 6: Welded Beam Function Test (dimensions = 3; n_comp = 2; extra_points = 4)" begin n_comp = 2 delta_x = 0.0001 extra_points = 4 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 37.5, atol=0.5) #rmse: 37.87 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 37.5, atol = 0.5) #rmse: 37.87 end ## sphere function tests function sphere_function(x) - return sum(x.^2) + return sum(x .^ 2) end ## 3D n = 100 lb = [-5.0, -5.0, -5.0] -ub = [5.0, 5.0 ,5.0] -x = sample(n,lb,ub,SobolSample()) +ub = [5.0, 5.0, 5.0] +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x),(size(x,1),1)) +y = reshape(sphere_function.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 7: Sphere Function Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin @@ -171,25 +172,25 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.001, atol=0.05) #rmse: 0.00083 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.001, atol = 0.05) #rmse: 0.00083 end ## 2D -n = 50 +n = 50 d = 2 lb = [-10.0, -10.0] ub = [10.0, 10.0] -x = sample(n,lb,ub,SobolSample()) +x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) -y = reshape(sphere_function.(x),(size(x,1),1)) +y = reshape(sphere_function.(x), (size(x, 1), 1)) xlimits = hcat(lb, ub) -n_test = 10 -x_test = sample(n_test,lb,ub,GoldenSample()) -X_test = vector_of_tuples_to_matrix(x_test) +n_test = 10 +x_test = sample(n_test, lb, ub, GoldenSample()) +X_test = vector_of_tuples_to_matrix(x_test) y_true = sphere_function.(x_test) @testset "Test 8: Sphere Function Test (dimensions = 2; n_comp = 2; extra_points = 2" begin @@ -197,16 +198,16 @@ y_true = sphere_function.(x_test) delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) + g = GEKPLS(X, y, grads, n_comp, delta_x, xlimits, extra_points, initial_theta) y_pred = g(X_test) - rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) - @test isapprox(rmse, 0.1, atol=0.5) #rmse: 0.0022 + rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) + @test isapprox(rmse, 0.1, atol = 0.5) #rmse: 0.0022 end @testset "Test 9: Add Point Test (dimensions = 3; n_comp = 2; extra_points = 2)" begin #first we create a surrogate model with just 3 input points initial_x_vec = [(1.0, 2.0, 3.0), (4.0, 5.0, 6.0), (7.0, 8.0, 9.0)] - initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec,1),1)) + initial_y = reshape(sphere_function.(initial_x_vec), (size(initial_x_vec, 1), 1)) initial_X = vector_of_tuples_to_matrix(initial_x_vec) initial_grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, initial_x_vec)) lb = [-5.0, -5.0, -5.0] @@ -216,22 +217,23 @@ end delta_x = 0.0001 extra_points = 2 initial_theta = 0.01 - g = GEKPLS(initial_X, initial_y, initial_grads, n_comp, delta_x, xlimits, extra_points, initial_theta) - n_test = 100 - x_test = sample(n_test,lb,ub,GoldenSample()) - X_test = vector_of_tuples_to_matrix(x_test) - y_true = sphere_function.(x_test) + g = GEKPLS(initial_X, initial_y, initial_grads, n_comp, delta_x, xlimits, extra_points, + initial_theta) + n_test = 100 + x_test = sample(n_test, lb, ub, GoldenSample()) + X_test = vector_of_tuples_to_matrix(x_test) + y_true = sphere_function.(x_test) y_pred1 = g(X_test) - rmse1 = sqrt(sum(((y_pred1 - y_true).^2)/n_test)) #rmse1 = 31.91 + rmse1 = sqrt(sum(((y_pred1 - y_true) .^ 2) / n_test)) #rmse1 = 31.91 #then we update the model with more points to see if performance improves n = 100 - x = sample(n,lb,ub,SobolSample()) + x = sample(n, lb, ub, SobolSample()) X = vector_of_tuples_to_matrix(x) grads = vector_of_tuples_to_matrix2(gradient.(sphere_function, x)) - y = reshape(sphere_function.(x),(size(x,1),1)) + y = reshape(sphere_function.(x), (size(x, 1), 1)) add_point!(g, X, y, grads) y_pred2 = g(X_test) - rmse2 = sqrt(sum(((y_pred2 - y_true).^2)/n_test)) #rmse2 = 0.0015 + rmse2 = sqrt(sum(((y_pred2 - y_true) .^ 2) / n_test)) #rmse2 = 0.0015 @test (rmse2 < rmse1) end diff --git a/test/runtests.jl b/test/runtests.jl index aa703cca4..e435b486f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,8 +4,8 @@ using SafeTestsets using Pkg function dev_subpkg(subpkg) - subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) - Pkg.develop(PackageSpec(path=subpkg_path)) + subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) + Pkg.develop(PackageSpec(path = subpkg_path)) end for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", "SurrogatesRandomForest", "SurrogatesSVM"] @@ -15,15 +15,15 @@ for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos", end end -@time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end -@time @safetestset "Radials.jl" begin include("Radials.jl") end -@time @safetestset "Kriging.jl" begin include("Kriging.jl") end -@time @safetestset "Sampling" begin include("sampling.jl") end -@time @safetestset "Optimization" begin include("optimization.jl") end -@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end -@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end -@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end -@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end +@time @safetestset "GEKPLS.jl" begin include("GEKPLS.jl") end +@time @safetestset "Radials.jl" begin include("Radials.jl") end +@time @safetestset "Kriging.jl" begin include("Kriging.jl") end +@time @safetestset "Sampling" begin include("sampling.jl") end +@time @safetestset "Optimization" begin include("optimization.jl") end +@time @safetestset "LinearSurrogate" begin include("linearSurrogate.jl") end +@time @safetestset "Lobachevsky" begin include("lobachevsky.jl") end +@time @safetestset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end +@time @safetestset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end # @time @safetestset "AD_Compatibility" begin include("AD_compatibility.jl") end @time @safetestset "Wendland" begin include("Wendland.jl") end @time @safetestset "VariableFidelity" begin include("VariableFidelity.jl") end From cd72b3f13ff1ffc6087745c587dd9a1e6955856a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 29 Jun 2022 05:38:59 -0400 Subject: [PATCH 12/12] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9893178ad..d87569e97 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"