From 56d0965cb89fbaf30c2bbb07cfb1c216032a4952 Mon Sep 17 00:00:00 2001 From: MRIDUL JAIN <105979087+Spinachboul@users.noreply.github.com> Date: Sun, 25 Feb 2024 14:29:28 +0530 Subject: [PATCH] Optimize GEK.jl --- src/GEK.jl | 188 +++++++++++++---------------------------------------- 1 file changed, 46 insertions(+), 142 deletions(-) diff --git a/src/GEK.jl b/src/GEK.jl index 790b2ba0..3ff0af53 100644 --- a/src/GEK.jl +++ b/src/GEK.jl @@ -1,3 +1,5 @@ +using LinearAlgebra + mutable struct GEK{X, Y, L, U, P, T, M, B, S, R} <: AbstractSurrogate x::X y::Y @@ -11,149 +13,45 @@ mutable struct GEK{X, Y, L, U, P, T, M, B, S, R} <: AbstractSurrogate inverse_of_R::R end -function _calc_gek_coeffs(x, y, p::Number, theta::Number) - nd1 = length(y) #2n - n = length(x) - R = zeros(eltype(x[1]), nd1, nd1) - - #top left - @inbounds for i in 1:n - for j in 1:n - R[i, j] = exp(-theta * abs(x[i] - x[j])^p) - end - end - #top right - @inbounds for i in 1:n - for j in (n + 1):nd1 - R[i, j] = 2 * theta * (x[i] - x[j - n]) * exp(-theta * abs(x[i] - x[j - n])^p) - end - end - #bottom left - @inbounds for i in (n + 1):nd1 - for j in 1:n - R[i, j] = -2 * theta * (x[i - n] - x[j]) * exp(-theta * abs(x[i - n] - x[j])^p) - end - end - #bottom right - @inbounds for i in (n + 1):nd1 - for j in (n + 1):nd1 - R[i, j] = -4 * theta * (x[i - n] - x[j - n])^2 * - exp(-theta * abs(x[i - n] - x[j - n])^p) - end - end - one = ones(eltype(x[1]), nd1, 1) - for i in (n + 1):nd1 - one[i] = zero(eltype(x[1])) - end - one_t = one' - inverse_of_R = inv(R) - mu = (one_t * inverse_of_R * y) / (one_t * inverse_of_R * one) - b = inverse_of_R * (y - one * mu) - sigma = ((y - one * mu)' * inverse_of_R * (y - one * mu)) / n - mu[1], b, sigma[1], inverse_of_R -end - -function std_error_at_point(k::GEK, val::Number) - # Check to make sure dimensions of input matches expected dimension of surrogate - _check_dimension(k, val) - - phi(z) = exp(-(abs(z))^k.p) - nd1 = length(k.y) - n = length(k.x) - r = zeros(eltype(k.x[1]), nd1, 1) - @inbounds for i in 1:n - r[i] = phi(val - k.x[i]) - end - one = ones(eltype(k.x[1]), nd1, 1) - @inbounds for i in (n + 1):nd1 - one[i] = zero(eltype(k.x[1])) - end - one_t = one' - a = r' * k.inverse_of_R * r - a = a[1] - b = one_t * k.inverse_of_R * one - b = b[1] - mean_squared_error = k.sigma * (1 - a + (1 - a)^2 / (b)) - return sqrt(abs(mean_squared_error)) -end - -function (k::GEK)(val::Number) - # Check to make sure dimensions of input matches expected dimension of surrogate - _check_dimension(k, val) - - phi = z -> exp(-(abs(z))^k.p) - n = length(k.x) - prediction = zero(eltype(k.x[1])) - @inbounds for i in 1:n - prediction = prediction + k.b[i] * phi(val - k.x[i]) - end - prediction = k.mu + prediction - return prediction -end - -function GEK(x, y, lb::Number, ub::Number; p = 1.0, theta = 1.0) - if length(x) != length(unique(x)) - println("There exists a repetition in the samples, cannot build Kriging.") - return - end - mu, b, sigma, inverse_of_R = _calc_gek_coeffs(x, y, p, theta) - return GEK(x, y, lb, ub, p, theta, mu, b, sigma, inverse_of_R) -end - function _calc_gek_coeffs(x, y, p, theta) nd = length(y) n = length(x) d = length(x[1]) - R = zeros(eltype(x[1]), nd, nd) + R = similar(zeros(eltype(x[1]), nd, nd)) - #top left - @inbounds for i in 1:n - for j in 1:n - R[i, j] = prod([exp(-theta[l] * (x[i][l] - x[j][l])) for l in 1:d]) - end + @inbounds for i in 1:n, j in 1:n + R[i, j] = prod(exp.(-theta .* (x[i] .- x[j]))) end - #top right - @inbounds for i in 1:n - jr = 1 - for j in (n + 1):d:nd - for l in 1:d - R[i, j + l - 1] = +2 * theta[l] * (x[i][l] - x[jr][l]) * R[i, jr] - end - jr = jr + 1 + jr = 1 + @inbounds for i in 1:n, j in (n + 1):d:nd + for l in 1:d + R[i, j + l - 1] = +2 * theta[l] * (x[i][l] - x[jr][l]) * R[i, jr] end + jr += 1 end - #bottom left - @inbounds for j in 1:n - ir = 1 - for i in (n + 1):d:nd - for l in 1:d - R[i + l - 1, j] = -2 * theta[l] * (x[ir][l] - x[j][l]) * R[ir, j] - end - ir = ir + 1 + ir = 1 + @inbounds for j in 1:n, i in (n + 1):d:nd + for l in 1:d + R[i + l - 1, j] = -2 * theta[l] * (x[ir][l] - x[j][l]) * R[ir, j] end + ir += 1 end - #bottom right ir = 1 - @inbounds for i in (n + 1):d:nd - for j in (n + 1):d:nd - jr = 1 - for l in 1:d - for k in 1:d - R[i + l - 1, j + k - 1] = -4 * theta[l] * theta[k] * - (x[ir][l] - x[jr][l]) * - (x[ir][k] - x[jr][k]) * R[ir, jr] - end - end - jr = jr + 1 + @inbounds for i in (n + 1):d:nd, j in (n + 1):d:nd + jr = 1 + for l in 1:d, k in 1:d + R[i + l - 1, j + k - 1] = -4 * theta[l] * theta[k] * + (x[ir][l] - x[jr][l]) * + (x[ir][k] - x[jr][k]) * R[ir, jr] end - ir = ir + +1 + jr += 1 end one = ones(eltype(x[1]), nd, 1) - for i in (n + 1):nd + @inbounds for i in (n + 1):nd one[i] = zero(eltype(x[1])) end one_t = one' @@ -161,22 +59,19 @@ function _calc_gek_coeffs(x, y, p, theta) mu = (one_t * inverse_of_R * y) / (one_t * inverse_of_R * one) b = inverse_of_R * (y - one * mu) sigma = ((y - one * mu)' * inverse_of_R * (y - one * mu)) / n - mu[1], b, sigma[1], inverse_of_R + return mu[1], b, sigma[1], inverse_of_R end function std_error_at_point(k::GEK, val) - - # Check to make sure dimensions of input matches expected dimension of surrogate _check_dimension(k, val) - nd1 = length(k.y) n = length(k.x) d = length(k.x[1]) - r = zeros(eltype(k.x[1]), nd1, 1) + r = similar(zeros(eltype(k.x[1]), nd1, 1)) @inbounds for i in 1:n sum = zero(eltype(k.x[1])) for l in 1:d - sum = sum + k.theta[l] * norm(val[l] - k.x[i][l])^(k.p[l]) + sum += k.theta[l] * norm(val[l] - k.x[i][l])^k.p[l] end r[i] = exp(-sum) end @@ -194,15 +89,19 @@ function std_error_at_point(k::GEK, val) end function (k::GEK)(val) - # Check to make sure dimensions of input matches expected dimension of surrogate _check_dimension(k, val) - d = length(val) n = length(k.x) - return k.mu + - sum(k.b[i] * - exp(-sum(k.theta[j] * norm(val[j] - k.x[i][j])^k.p[j] for j in 1:d)) - for i in 1:n) + phi = z -> exp(-abs(z)^k.p) + prediction = zero(eltype(k.x[1])) + @inbounds for i in 1:n + sum = zero(eltype(k.x[1])) + for j in 1:d + sum += k.theta[j] * norm(val[j] - k.x[i][j])^k.p[j] + end + prediction += k.b[i] * phi(sum) + end + return k.mu + prediction end function GEK(x, y, lb, ub; p = collect(one.(x[1])), theta = collect(one.(x[1]))) @@ -211,7 +110,7 @@ function GEK(x, y, lb, ub; p = collect(one.(x[1])), theta = collect(one.(x[1]))) return end mu, b, sigma, inverse_of_R = _calc_gek_coeffs(x, y, p, theta) - GEK(x, y, lb, ub, p, theta, mu, b, sigma, inverse_of_R) + return GEK(x, y, lb, ub, p, theta, mu, b, sigma, inverse_of_R) end function add_point!(k::GEK, new_x, new_y) @@ -219,16 +118,21 @@ function add_point!(k::GEK, new_x, new_y) println("Adding a sample that already exists, cannot build Kriging.") return end - if (length(new_x) == 1 && length(new_x[1]) == 1) || - (length(new_x) > 1 && length(new_x[1]) == 1 && length(k.theta) > 1) - n = length(k.x) + n = length(k.x) + if length(new_x) == 1 || length(new_x[1]) == 1 k.x = insert!(k.x, n + 1, new_x) k.y = insert!(k.y, n + 1, new_y) else - n = length(k.x) k.x = insert!(k.x, n + 1, new_x) k.y = insert!(k.y, n + 1, new_y) end k.mu, k.b, k.sigma, k.inverse_of_R = _calc_gek_coeffs(k.x, k.y, k.p, k.theta) nothing end + +function _check_dimension(k::GEK, val) + d = length(val) + if length(k.x[1]) != d + throw(DimensionMismatch("Input dimension does not match surrogate dimension")) + end +end