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