From 2eb4353eaa4dd317944e1687861d0efd2fecd955 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 20 May 2023 18:13:11 -0500 Subject: [PATCH 01/16] Added 6 parallelization strategies for Kriging --- src/Kriging.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/Kriging.jl b/src/Kriging.jl index e0b3a3671..6dfedf38a 100644 --- a/src/Kriging.jl +++ b/src/Kriging.jl @@ -261,3 +261,39 @@ function add_point!(k::Kriging, new_x, new_y) k.mu, k.b, k.sigma, k.inverse_of_R = _calc_kriging_coeffs(k.x, k.y, k.p, k.theta) nothing end + +# Minimum Constant Liar +function CLmin!(tmp_k::Kriging, k::Kriging, new_x) + new_y = minimum(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Maximum Constant Liar +function CLmax!(tmp_k::Kriging, k::Kriging, new_x) + new_y = maximum(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Mean Constant Liar +function CLmean!(tmp_k::Kriging, k::Kriging, new_x) + new_y = mean(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer +function KB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer Upper Bound +function KBUB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) + 3 * std_error_at_point(k, new_x) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer Lower Bound +function KBLB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) - 3 * std_error_at_point(k, new_x) + add_point!(tmp_k, new_x, new_y) +end From b31e29cb0cc751cd7d27331d922d5e3e82065a84 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 20 May 2023 18:13:35 -0500 Subject: [PATCH 02/16] Added Ask function for Kriging --- src/Optimization.jl | 64 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/src/Optimization.jl b/src/Optimization.jl index a2fb005d2..77a12a2f0 100644 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -569,6 +569,70 @@ function surrogate_optimize(obj::Function, ::EI, lb::Number, ub::Number, krig, println("Completed maximum number of iterations") end +function Ask(::EI, lb::Number, ub::Number, krig, + sample_type::SamplingAlgorithm, n_parallel, strategy!; + num_new_samples = 100) + + dtol = 1e-3 * norm(ub - lb) + eps = 0.01 + + tmp_krig = deepcopy(krig) # Temporary copy of the kriging model to store virtual points + + new_x_max = zeros(eltype(tmp_krig.x[1]), n_parallel) # New x point + new_EI_max = zeros(eltype(tmp_krig.x[1]), n_parallel) # EI at new x point + + for i in 1:n_parallel + # Sample lots of points from the design space -- we will evaluate the EI function at these points + new_sample = sample(num_new_samples, lb, ub, sample_type) + + # Find the best point so far + f_min = minimum(tmp_krig.y) + + # Allocate some arrays + evaluations = zeros(eltype(tmp_krig.x[1]), num_new_samples) # Holds EI function evaluations + point_found = false # Whether we have found a new point to test + while point_found == false + # For each point in the sample set, evaluate the Expected Improvement function + for j in 1:length(new_sample) + std = std_error_at_point(tmp_krig, new_sample[j]) + u = tmp_krig(new_sample[j]) + if abs(std) > 1e-6 + z = (f_min - u - eps) / std + else + z = 0 + end + # Evaluate EI at point new_sample[j] + evaluations[j] = (f_min - u - eps) * cdf(Normal(), z) + + std * pdf(Normal(), z) + end + # find the sample which maximizes the EI function + index_max = argmax(evaluations) + x_new = new_sample[index_max] # x point which maximized EI + y_new = maximum(evaluations) # EI at the new point + diff_x = abs.(tmp_krig.x .- x_new) + bit_x = diff_x .> dtol + #new_min_x has to have some distance from tmp_krig.x + if false in bit_x + #The new_point is not actually that new, discard it! + deleteat!(evaluations, index_max) + deleteat!(new_sample, index_max) + if length(new_sample) == 0 + println("Out of sampling points") + index = argmin(tmp_krig.y) + return (tmp_krig.x[index], tmp_krig.y[index]) + end + else + point_found = true + new_x_max[i] = x_new + new_EI_max[i] = y_new + strategy!(tmp_krig, krig, x_new) + end + end + end + + return (new_x_max, new_EI_max) +end + """ This is an implementation of Expected Improvement (EI), arguably the most popular acquisition function in Bayesian optimization. From cd1973c13f8ab2b77a518505c0bbd08e5dfa82d4 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 20 May 2023 18:18:31 -0500 Subject: [PATCH 03/16] Exported relevant functions --- src/Surrogates.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Surrogates.jl b/src/Surrogates.jl index eaf50ed44..20f94577b 100644 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -88,6 +88,10 @@ export LobachevskyStructure, NeuralStructure, RandomForestStructure, export WendlandStructure export AbstractSurrogate, SamplingAlgorithm export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point +# Parallelization Strategies +export KB!, KBLB!, KBUB!, CLmax!, CLmean!, CLmin! +export Ask + # radial basis functions export linearRadial, cubicRadial, multiquadricRadial, thinplateRadial From 8528a7332284c247238fd59475b46f67057f71b8 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 10 Jun 2023 20:32:50 -0500 Subject: [PATCH 04/16] Removed unnecessary arguments --- src/Optimization.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Optimization.jl b/src/Optimization.jl index 77a12a2f0..37c29599a 100644 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -569,16 +569,18 @@ function surrogate_optimize(obj::Function, ::EI, lb::Number, ub::Number, krig, println("Completed maximum number of iterations") end -function Ask(::EI, lb::Number, ub::Number, krig, - sample_type::SamplingAlgorithm, n_parallel, strategy!; +function Ask(::EI, krig, sample_type::SamplingAlgorithm, n_parallel::Number, strategy!; num_new_samples = 100) + lb = krig.lb + ub = krig.ub + dtol = 1e-3 * norm(ub - lb) eps = 0.01 tmp_krig = deepcopy(krig) # Temporary copy of the kriging model to store virtual points - new_x_max = zeros(eltype(tmp_krig.x[1]), n_parallel) # New x point + new_x_max = Vector{typeof(tmp_krig.x[1])}(undef, n_parallel) # New x point new_EI_max = zeros(eltype(tmp_krig.x[1]), n_parallel) # EI at new x point for i in 1:n_parallel @@ -593,7 +595,7 @@ function Ask(::EI, lb::Number, ub::Number, krig, point_found = false # Whether we have found a new point to test while point_found == false # For each point in the sample set, evaluate the Expected Improvement function - for j in 1:length(new_sample) + for j in eachindex(new_sample) std = std_error_at_point(tmp_krig, new_sample[j]) u = tmp_krig(new_sample[j]) if abs(std) > 1e-6 @@ -609,8 +611,8 @@ function Ask(::EI, lb::Number, ub::Number, krig, index_max = argmax(evaluations) x_new = new_sample[index_max] # x point which maximized EI y_new = maximum(evaluations) # EI at the new point - diff_x = abs.(tmp_krig.x .- x_new) - bit_x = diff_x .> dtol + diff_x = [abs.(prev_point .- x_new) for prev_point in tmp_krig.x] + bit_x = [diff_x_point .> dtol for diff_x_point in diff_x] #new_min_x has to have some distance from tmp_krig.x if false in bit_x #The new_point is not actually that new, discard it! From de1f4a371f13328cc02622e3c525cc19189229ec Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 10 Jun 2023 20:33:05 -0500 Subject: [PATCH 05/16] Added tests for Ask --- test/Ask.jl | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 test/Ask.jl diff --git a/test/Ask.jl b/test/Ask.jl new file mode 100644 index 000000000..169b42eab --- /dev/null +++ b/test/Ask.jl @@ -0,0 +1,41 @@ +using Surrogates +using Test + + +#1D +lb = 0.0 +ub = 10.0 +f = x -> log(x) * exp(x) +x = sample(5, lb, ub, SobolSample()) +y = f.(x) + + +# Test lengths of new_x and EI (1D) + +my_k = Kriging(x, y, lb, ub) +new_x, eis = Ask(EI(), my_k, SobolSample(), 3, CLmean!) + +@test length(new_x) == 3 +@test length(eis) == 3 + + +# Test lengths of new_x and EI (ND) + +lb = [0.0, 0.0, 1.0] +ub = [5.0, 7.5, 10.0] +x = sample(5, lb, ub, SobolSample()) +f = x -> x[1] + x[2] * x[3] +y = f.(x) + +my_k = Kriging(x, y, lb, ub) + +new_x, eis = Ask(EI(), my_k, SobolSample(), 5, CLmean!) + +@test length(new_x) == 5 +@test length(eis) == 5 + +@test length(new_x[1]) == 3 + +# # Check hyperparameter validation for Ask +@test_throws ArgumentError new_x, eis = Ask(EI(), my_k, SobolSample(), -1, CLmean!) + From d114a20855cfb9a02a5e4f34034eb637f704e03d Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Fri, 11 Aug 2023 19:10:40 -0500 Subject: [PATCH 06/16] Moved parallel strategies to their own file --- src/Kriging.jl | 36 ------------------------------------ src/VirtualStrategy.jl | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 36 deletions(-) create mode 100644 src/VirtualStrategy.jl diff --git a/src/Kriging.jl b/src/Kriging.jl index 6dfedf38a..e0b3a3671 100644 --- a/src/Kriging.jl +++ b/src/Kriging.jl @@ -261,39 +261,3 @@ function add_point!(k::Kriging, new_x, new_y) k.mu, k.b, k.sigma, k.inverse_of_R = _calc_kriging_coeffs(k.x, k.y, k.p, k.theta) nothing end - -# Minimum Constant Liar -function CLmin!(tmp_k::Kriging, k::Kriging, new_x) - new_y = minimum(k.y) - add_point!(tmp_k, new_x, new_y) -end - -# Maximum Constant Liar -function CLmax!(tmp_k::Kriging, k::Kriging, new_x) - new_y = maximum(k.y) - add_point!(tmp_k, new_x, new_y) -end - -# Mean Constant Liar -function CLmean!(tmp_k::Kriging, k::Kriging, new_x) - new_y = mean(k.y) - add_point!(tmp_k, new_x, new_y) -end - -# Kriging Believer -function KB!(tmp_k::Kriging, k::Kriging, new_x) - new_y = k(new_x) - add_point!(tmp_k, new_x, new_y) -end - -# Kriging Believer Upper Bound -function KBUB!(tmp_k::Kriging, k::Kriging, new_x) - new_y = k(new_x) + 3 * std_error_at_point(k, new_x) - add_point!(tmp_k, new_x, new_y) -end - -# Kriging Believer Lower Bound -function KBLB!(tmp_k::Kriging, k::Kriging, new_x) - new_y = k(new_x) - 3 * std_error_at_point(k, new_x) - add_point!(tmp_k, new_x, new_y) -end diff --git a/src/VirtualStrategy.jl b/src/VirtualStrategy.jl new file mode 100644 index 000000000..e2558cb45 --- /dev/null +++ b/src/VirtualStrategy.jl @@ -0,0 +1,35 @@ +# Minimum Constant Liar +function CLmin!(tmp_k::Kriging, k::Kriging, new_x) + new_y = minimum(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Maximum Constant Liar +function CLmax!(tmp_k::Kriging, k::Kriging, new_x) + new_y = maximum(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Mean Constant Liar +function CLmean!(tmp_k::Kriging, k::Kriging, new_x) + new_y = mean(k.y) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer +function KB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer Upper Bound +function KBUB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) + 3 * std_error_at_point(k, new_x) + add_point!(tmp_k, new_x, new_y) +end + +# Kriging Believer Lower Bound +function KBLB!(tmp_k::Kriging, k::Kriging, new_x) + new_y = k(new_x) - 3 * std_error_at_point(k, new_x) + add_point!(tmp_k, new_x, new_y) +end \ No newline at end of file From 3ad117f314d2e07983c0b7c1e8d9106e0d72aef8 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Fri, 11 Aug 2023 20:17:20 -0500 Subject: [PATCH 07/16] Fixed argument types for virtual point functions --- src/VirtualStrategy.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/VirtualStrategy.jl b/src/VirtualStrategy.jl index e2558cb45..3c0a14a2e 100644 --- a/src/VirtualStrategy.jl +++ b/src/VirtualStrategy.jl @@ -1,19 +1,19 @@ # Minimum Constant Liar -function CLmin!(tmp_k::Kriging, k::Kriging, new_x) - new_y = minimum(k.y) - add_point!(tmp_k, new_x, new_y) +function CLmin!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) + new_y = minimum(surr.y) + add_point!(tmp_surr, new_x, new_y) end # Maximum Constant Liar -function CLmax!(tmp_k::Kriging, k::Kriging, new_x) - new_y = maximum(k.y) - add_point!(tmp_k, new_x, new_y) +function CLmax!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) + new_y = maximum(surr.y) + add_point!(tmp_surr, new_x, new_y) end # Mean Constant Liar -function CLmean!(tmp_k::Kriging, k::Kriging, new_x) - new_y = mean(k.y) - add_point!(tmp_k, new_x, new_y) +function CLmean!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) + new_y = mean(surr.y) + add_point!(tmp_surr, new_x, new_y) end # Kriging Believer From a7c3308aa5433ebbb717cd6b2dc555de26540cb4 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Thu, 14 Sep 2023 19:28:36 -0500 Subject: [PATCH 08/16] Added Ask-Tell for SRBF merit function --- src/Optimization.jl | 214 +++++++++++++++++++++++++++++++++++++++++++- src/Surrogates.jl | 1 + test/Ask.jl | 24 ++++- 3 files changed, 235 insertions(+), 4 deletions(-) mode change 100644 => 100755 src/Optimization.jl mode change 100644 => 100755 src/Surrogates.jl mode change 100644 => 100755 test/Ask.jl diff --git a/src/Optimization.jl b/src/Optimization.jl old mode 100644 new mode 100755 index 37c29599a..6bd966da3 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -375,6 +375,218 @@ function surrogate_optimize(obj::Function, ::SRBF, lb::Number, ub::Number, end end +# SRBF ND +function Ask(::SRBF, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgorithm, n_parallel, strategy!; + num_new_samples = 500) + + scale = 0.2 + success = 0 + failure = 0 + w_range = [0.3, 0.5, 0.7, 0.95] + w_cycle = Iterators.cycle(w_range) + + w, state = iterate(w_cycle) + + #Vector containing size in each direction + box_size = lb - ub + success = 0 + failures = 0 + dtol = 1e-3 * norm(ub - lb) + d = length(surr.x) + incumbent_x = surr.x[argmin(surr.y)] + + new_lb = incumbent_x .- 3 * scale * norm(incumbent_x .- lb) + new_ub = incumbent_x .+ 3 * scale * norm(incumbent_x .- ub) + + @inbounds for i in 1:length(new_lb) + if new_lb[i] < lb[i] + new_lb = collect(new_lb) + new_lb[i] = lb[i] + end + if new_ub[i] > ub[i] + new_ub = collect(new_ub) + new_ub[i] = ub[i] + end + end + + new_sample = sample(num_new_samples, new_lb, new_ub, sample_type) + s = zeros(eltype(surr.x[1]), num_new_samples) + for j in 1:num_new_samples + s[j] = surr(new_sample[j]) + end + s_max = maximum(s) + s_min = minimum(s) + + d_min = norm(box_size .+ 1) + d_max = 0.0 + for r in 1:length(surr.x) + for c in 1:num_new_samples + distance_rc = norm(surr.x[r] .- new_sample[c]) + if distance_rc > d_max + d_max = distance_rc + end + if distance_rc < d_min + d_min = distance_rc + end + end + end + + tmp_surr = deepcopy(surr) + + + new_addition = 0 + diff_x = zeros(eltype(surr.x[1]), d) + + evaluation_of_merit_function = zeros(float(eltype(surr.x[1])), num_new_samples) + proposed_points_x = Vector{typeof(surr.x[1])}(undef, n_parallel) + merit_of_proposed_points = zeros(Float64, n_parallel) + + while new_addition < n_parallel + #find minimum + + @inbounds for r in 1:num_new_samples + evaluation_of_merit_function[r] = merit_function(new_sample[r], w, tmp_surr, + s_max, s_min, d_max, d_min, + box_size) + end + + min_index = argmin(evaluation_of_merit_function) + new_min_x = new_sample[min_index] + min_x_merit = evaluation_of_merit_function[min_index] + + for l in 1:d + diff_x[l] = norm(surr.x[l] .- new_min_x) + end + bit_x = diff_x .> dtol + #new_min_x has to have some distance from krig.x + if false in bit_x + #The new_point is not actually that new, discard it! + + deleteat!(evaluation_of_merit_function, min_index) + deleteat!(new_sample, min_index) + + if length(new_sample) == 0 + println("Out of sampling points") + index = argmin(surr.y) + return (surr.x[index], surr.y[index]) + end + else + new_addition += 1 + proposed_points_x[new_addition] = new_min_x + merit_of_proposed_points[new_addition] = min_x_merit + + # Update temporary surrogate using provided strategy + strategy!(tmp_surr, surr, new_min_x) + end + + #4) Update w + w, state = iterate(w_cycle, state) + end + + return (proposed_points_x, merit_of_proposed_points) +end + +# Ask SRBF 1D +function Ask(::SRBF, lb::Number, ub::Number, surr::AbstractSurrogate, + sample_type::SamplingAlgorithm, n_parallel, strategy!; + num_new_samples = 500) + scale = 0.2 + success = 0 + w_range = [0.3, 0.5, 0.7, 0.95] + w_cycle = Iterators.cycle(w_range) + + w, state = iterate(w_cycle) + + box_size = lb - ub + success = 0 + failures = 0 + dtol = 1e-3 * norm(ub - lb) + num_of_iterations = 0 + + #1) Sample near incumbent (the 2 fraction is arbitrary here) + incumbent_x = surr.x[argmin(surr.y)] + + new_lb = incumbent_x - scale * norm(incumbent_x - lb) + new_ub = incumbent_x + scale * norm(incumbent_x - ub) + if new_lb < lb + new_lb = lb + end + if new_ub > ub + new_ub = ub + end + + new_sample = sample(num_new_samples, new_lb, new_ub, sample_type) + + #2) Create merit function + s = zeros(eltype(surr.x[1]), num_new_samples) + for j in 1:num_new_samples + s[j] = surr(new_sample[j]) + end + s_max = maximum(s) + s_min = minimum(s) + + d_min = box_size + 1 + d_max = 0.0 + for r in 1:length(surr.x) + for c in 1:num_new_samples + distance_rc = norm(surr.x[r] - new_sample[c]) + if distance_rc > d_max + d_max = distance_rc + end + if distance_rc < d_min + d_min = distance_rc + end + end + end + + new_addition = 0 + proposed_points_x = zeros(eltype(new_sample[1]), n_parallel) + merit_of_proposed_points = zeros(eltype(new_sample[1]), n_parallel) + + # Temporary surrogate for virtual points + tmp_surr = deepcopy(surr) + + # Loop until we have n_parallel new points + while new_addition < n_parallel + + #3) Evaluate merit function at the sampled points in parallel + evaluation_of_merit_function = merit_function.(new_sample, w, tmp_surr, s_max, + s_min, d_max, d_min, box_size) + + #find minimum + min_index = argmin(evaluation_of_merit_function) + new_min_x = new_sample[min_index] + min_x_merit = evaluation_of_merit_function[min_index] + + diff_x = abs.(tmp_surr.x .- new_min_x) + bit_x = diff_x .> dtol + #new_min_x has to have some distance from krig.x + if false in bit_x + #The new_point is not actually that new, discard it! + deleteat!(evaluation_of_merit_function, min_index) + deleteat!(new_sample, min_index) + if length(new_sample) == 0 + println("Out of sampling points") + return (proposed_points_x[1:new_addition], + merit_of_proposed_points[1:new_addition]) + end + else + new_addition += 1 + proposed_points_x[new_addition] = new_min_x + merit_of_proposed_points[new_addition] = min_x_merit + + # Update temporary surrogate using provided strategy + strategy!(tmp_surr, surr, new_min_x) + end + + #4) Update w + w, state = iterate(w_cycle, state) + end + + return (proposed_points_x, merit_of_proposed_points) +end + + """ This is an implementation of Lower Confidence Bound (LCB), a popular acquisition function in Bayesian optimization. @@ -569,7 +781,7 @@ function surrogate_optimize(obj::Function, ::EI, lb::Number, ub::Number, krig, println("Completed maximum number of iterations") end -function Ask(::EI, krig, sample_type::SamplingAlgorithm, n_parallel::Number, strategy!; +function Ask(::EI, lb, ub, krig, sample_type::SamplingAlgorithm, n_parallel::Number, strategy!; num_new_samples = 100) lb = krig.lb diff --git a/src/Surrogates.jl b/src/Surrogates.jl old mode 100644 new mode 100755 index 20f94577b..6c6a7706e --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -18,6 +18,7 @@ include("VariableFidelity.jl") include("Earth.jl") include("GEK.jl") include("GEKPLS.jl") +include("VirtualStrategy.jl") current_surrogates = ["Kriging", "LinearSurrogate", "LobachevskySurrogate", "NeuralSurrogate", diff --git a/test/Ask.jl b/test/Ask.jl old mode 100644 new mode 100755 index 169b42eab..ec7abebf1 --- a/test/Ask.jl +++ b/test/Ask.jl @@ -1,5 +1,6 @@ using Surrogates using Test +using Revise #1D @@ -13,11 +14,18 @@ y = f.(x) # Test lengths of new_x and EI (1D) my_k = Kriging(x, y, lb, ub) -new_x, eis = Ask(EI(), my_k, SobolSample(), 3, CLmean!) +new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) @test length(new_x) == 3 @test length(eis) == 3 +# Test lengths of new_x and SRBF (1D) + +my_surr = RadialBasis(x, y, lb, ub) +new_x, eis = Ask(SRBF(), lb, ub, my_surr, SobolSample(), 3, CLmean!) +@test length(new_x) == 3 +@test length(eis) == 3 + # Test lengths of new_x and EI (ND) @@ -29,7 +37,17 @@ y = f.(x) my_k = Kriging(x, y, lb, ub) -new_x, eis = Ask(EI(), my_k, SobolSample(), 5, CLmean!) +new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 5, CLmean!) + +@test length(new_x) == 5 +@test length(eis) == 5 + +@test length(new_x[1]) == 3 + +# Test lengths of new_x and SRBF (ND) + +my_surr = RadialBasis(x, y, lb, ub) +new_x, eis = Ask(SRBF(), lb, ub, my_surr, SobolSample(), 5, CLmean!) @test length(new_x) == 5 @test length(eis) == 5 @@ -37,5 +55,5 @@ new_x, eis = Ask(EI(), my_k, SobolSample(), 5, CLmean!) @test length(new_x[1]) == 3 # # Check hyperparameter validation for Ask -@test_throws ArgumentError new_x, eis = Ask(EI(), my_k, SobolSample(), -1, CLmean!) +@test_throws ArgumentError new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), -1, CLmean!) From 95f8b5a77fd2ff11fa818f94caef5aa3e9d8d0b7 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Sat, 16 Sep 2023 22:11:35 -0500 Subject: [PATCH 09/16] Wrote documentation --- docs/src/parallel.md | 56 ++++++++++++++++++++++++++++++++++++++++++++ src/Optimization.jl | 3 ++- 2 files changed, 58 insertions(+), 1 deletion(-) create mode 100755 docs/src/parallel.md diff --git a/docs/src/parallel.md b/docs/src/parallel.md new file mode 100755 index 000000000..99d0004df --- /dev/null +++ b/docs/src/parallel.md @@ -0,0 +1,56 @@ +# Parallel Optimization + +There are some situations where it can be beneficial to run multiple optimizations in parallel. For example, if your objective function is very expensive to evaluate, you may want to run multiple evaluations in parallel. + +## Ask-Tell Interface + +To enable parallel optimization, we make use of an Ask-Tell interface. The user will construct the initial surrogate model the same way as for non-parallel surrogate models, but instead of using `surrogate_optimize`, the user will use `Ask`. This will return the coordinates of points that the optimizer has determined are most useful to evaluate next. How the user evaluates these points is up to them. The Ask-Tell interface requires more manual control than `surrogate_optimize`, but it allows for more flexibility. After the point has been evaluated, the user will *tell* the surrogate model the new points with the `add_point!` function. + +## Virtual Points + +To ensure that points of interest returned by `Ask` are sufficiently far from each other, the function makes use of *virtual points*. They are used as follows: +1. `Ask` is told to return `n` points. +2. The point with the highest merit function value is selected. +3. This point is now treated as a virtual point and is assigned a temporary value that changes the landscape of the merit function. How the the temporary value is chosen depends on the strategy used. (see below) +4. The point with the new highest merit is selected. +5. The process is repeated until `n` points have been selected. + +The following strategies are available for virtual point selection for all optimization algorithms: + +- "Minimum Constant Liar (CLmin)": + - The virtual point is assigned using the lowest known value of the merit function across all evaluated points. +- "Mean Constant Liar (CLmean)": + - The virtual point is assigned using the mean of the merit function across all evaluated points. +- "Maximum Constant Liar (CLmax)": + - The virtual point is assigned using the great known value of the merit function across all evaluated points. + +For Kriging surrogates, specifically, the above and follow strategies are available: + +- "Kriging Believer (KB)": + - The virtual point is assigned using the mean of the Kriging surrogate at the virtual point. +- "Kriging Believer Upper Bound (KBUB)": + - The virtual point is assigned using 3$\sigma$ above the mean of the Kriging surrogate at the virtual point. +- "Kriging Believer Lower Bound (KBLB)": + - The virtual point is assigned using 3$\sigma$ below the mean of the Kriging surrogate at the virtual point. + + +In general, CLmin and KBLB tend to favor exploitation while CLmax and KBUB tend to favor exploration. CLmean and KB tend to be a compromise between the two. + +## Examples + +```jl +using Surrogates + +lb = 0.0 +ub = 10.0 +f = x -> log(x) * exp(x) +x = sample(5, lb, ub, SobolSample()) +y = f.(x) + +my_k = Kriging(x, y, lb, ub) + +for _ in 1:10 + new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) + add_point!(my_k, new_x, f.(new_x)) +end +``` diff --git a/src/Optimization.jl b/src/Optimization.jl index 6bd966da3..0cc60a15e 100755 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -375,7 +375,7 @@ function surrogate_optimize(obj::Function, ::SRBF, lb::Number, ub::Number, end end -# SRBF ND +# Ask SRBF ND function Ask(::SRBF, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgorithm, n_parallel, strategy!; num_new_samples = 500) @@ -781,6 +781,7 @@ function surrogate_optimize(obj::Function, ::EI, lb::Number, ub::Number, krig, println("Completed maximum number of iterations") end +# Ask EI 1D & ND function Ask(::EI, lb, ub, krig, sample_type::SamplingAlgorithm, n_parallel::Number, strategy!; num_new_samples = 100) From 48d4601d330b44cf318f26c5318d397e25dc5ca6 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 17:49:55 -0500 Subject: [PATCH 10/16] Renamed test file for parallelization to parallel.jl --- test/{Ask.jl => parallel.jl} | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) rename test/{Ask.jl => parallel.jl} (57%) diff --git a/test/Ask.jl b/test/parallel.jl similarity index 57% rename from test/Ask.jl rename to test/parallel.jl index ec7abebf1..7cc133e52 100755 --- a/test/Ask.jl +++ b/test/parallel.jl @@ -13,8 +13,10 @@ y = f.(x) # Test lengths of new_x and EI (1D) +# TODO + my_k = Kriging(x, y, lb, ub) -new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) +new_x, eis = potential_optimal_points(EI(), MeanConstantLiar(), lb, ub, my_k, SobolSample(), 3) @test length(new_x) == 3 @test length(eis) == 3 @@ -22,7 +24,7 @@ new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) # Test lengths of new_x and SRBF (1D) my_surr = RadialBasis(x, y, lb, ub) -new_x, eis = Ask(SRBF(), lb, ub, my_surr, SobolSample(), 3, CLmean!) +new_x, eis = potential_optimal_points(SRBF(), MeanConstantLiar(), lb, ub, my_surr, SobolSample(), 3) @test length(new_x) == 3 @test length(eis) == 3 @@ -37,7 +39,7 @@ y = f.(x) my_k = Kriging(x, y, lb, ub) -new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 5, CLmean!) +new_x, eis = potential_optimal_points(EI(), MeanConstantLiar(), lb, ub, my_k, SobolSample(), 5) @test length(new_x) == 5 @test length(eis) == 5 @@ -47,13 +49,13 @@ new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 5, CLmean!) # Test lengths of new_x and SRBF (ND) my_surr = RadialBasis(x, y, lb, ub) -new_x, eis = Ask(SRBF(), lb, ub, my_surr, SobolSample(), 5, CLmean!) +new_x, eis = potential_optimal_points(SRBF(), MeanConstantLiar(), lb, ub, my_surr, SobolSample(), 5) @test length(new_x) == 5 @test length(eis) == 5 @test length(new_x[1]) == 3 -# # Check hyperparameter validation for Ask -@test_throws ArgumentError new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), -1, CLmean!) +# # Check hyperparameter validation for potential_optimal_points +@test_throws ArgumentError new_x, eis = potential_optimal_points(EI(), MeanConstantLiar(), lb, ub, my_k, SobolSample(), -1) From 6c4a39361f180ae67477fb4319e58b306b27696b Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 17:50:54 -0500 Subject: [PATCH 11/16] Renamed functions to be more explicit and changed signatures to dispatch on types --- src/Optimization.jl | 22 +++++++++++++++------- src/Surrogates.jl | 5 +++-- src/VirtualStrategy.jl | 12 ++++++------ 3 files changed, 24 insertions(+), 15 deletions(-) diff --git a/src/Optimization.jl b/src/Optimization.jl index 0cc60a15e..65e8fc77f 100755 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -3,6 +3,14 @@ using Distributions using Zygote abstract type SurrogateOptimizationAlgorithm end +abstract type ParallelStrategy end + +struct KrigingBeliever <: ParallelStrategy end +struct KrigingBelieverUpperBound <: ParallelStrategy end +struct KrigingBelieverLowerBound <: ParallelStrategy end +struct MinimumConstantLiar <: ParallelStrategy end +struct MaximumConstantLiar <: ParallelStrategy end +struct MeanConstantLiar <: ParallelStrategy end #single objective optimization struct SRBF <: SurrogateOptimizationAlgorithm end @@ -376,7 +384,7 @@ function surrogate_optimize(obj::Function, ::SRBF, lb::Number, ub::Number, end # Ask SRBF ND -function Ask(::SRBF, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgorithm, n_parallel, strategy!; +function potential_optimal_points(::SRBF, strategy, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgorithm, n_parallel; num_new_samples = 500) scale = 0.2 @@ -476,7 +484,7 @@ function Ask(::SRBF, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgor merit_of_proposed_points[new_addition] = min_x_merit # Update temporary surrogate using provided strategy - strategy!(tmp_surr, surr, new_min_x) + calculate_liars(strategy, tmp_surr, surr, new_min_x) end #4) Update w @@ -487,8 +495,8 @@ function Ask(::SRBF, lb, ub, surr::AbstractSurrogate, sample_type::SamplingAlgor end # Ask SRBF 1D -function Ask(::SRBF, lb::Number, ub::Number, surr::AbstractSurrogate, - sample_type::SamplingAlgorithm, n_parallel, strategy!; +function potential_optimal_points(::SRBF, strategy, lb::Number, ub::Number, surr::AbstractSurrogate, + sample_type::SamplingAlgorithm, n_parallel; num_new_samples = 500) scale = 0.2 success = 0 @@ -576,7 +584,7 @@ function Ask(::SRBF, lb::Number, ub::Number, surr::AbstractSurrogate, merit_of_proposed_points[new_addition] = min_x_merit # Update temporary surrogate using provided strategy - strategy!(tmp_surr, surr, new_min_x) + calculate_liars(strategy, tmp_surr, surr, new_min_x) end #4) Update w @@ -782,7 +790,7 @@ function surrogate_optimize(obj::Function, ::EI, lb::Number, ub::Number, krig, end # Ask EI 1D & ND -function Ask(::EI, lb, ub, krig, sample_type::SamplingAlgorithm, n_parallel::Number, strategy!; +function potential_optimal_points(::EI, strategy, lb, ub, krig, sample_type::SamplingAlgorithm, n_parallel::Number; num_new_samples = 100) lb = krig.lb @@ -840,7 +848,7 @@ function Ask(::EI, lb, ub, krig, sample_type::SamplingAlgorithm, n_parallel::Num point_found = true new_x_max[i] = x_new new_EI_max[i] = y_new - strategy!(tmp_krig, krig, x_new) + calculate_liars(strategy, tmp_krig, krig, x_new) end end end diff --git a/src/Surrogates.jl b/src/Surrogates.jl index 6c6a7706e..45df925b3 100755 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -90,8 +90,9 @@ export WendlandStructure export AbstractSurrogate, SamplingAlgorithm export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point # Parallelization Strategies -export KB!, KBLB!, KBUB!, CLmax!, CLmean!, CLmin! -export Ask +export potential_optimal_points +export MinimumConstantLiar, MaximumConstantLiar, MeanConstantLiar, KrigingBeliever, + KrigingBelieverUpperBound, KrigingBelieverLowerBound # radial basis functions export linearRadial, cubicRadial, multiquadricRadial, thinplateRadial diff --git a/src/VirtualStrategy.jl b/src/VirtualStrategy.jl index 3c0a14a2e..cc63a0c19 100644 --- a/src/VirtualStrategy.jl +++ b/src/VirtualStrategy.jl @@ -1,35 +1,35 @@ # Minimum Constant Liar -function CLmin!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) +function calculate_liars(::MinimumConstantLiar, tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) new_y = minimum(surr.y) add_point!(tmp_surr, new_x, new_y) end # Maximum Constant Liar -function CLmax!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) +function calculate_liars(::MaximumConstantLiar, tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) new_y = maximum(surr.y) add_point!(tmp_surr, new_x, new_y) end # Mean Constant Liar -function CLmean!(tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) +function calculate_liars(::MeanConstantLiar, tmp_surr::AbstractSurrogate, surr::AbstractSurrogate, new_x) new_y = mean(surr.y) add_point!(tmp_surr, new_x, new_y) end # Kriging Believer -function KB!(tmp_k::Kriging, k::Kriging, new_x) +function calculate_liars(::KrigingBeliever, tmp_k::Kriging, k::Kriging, new_x) new_y = k(new_x) add_point!(tmp_k, new_x, new_y) end # Kriging Believer Upper Bound -function KBUB!(tmp_k::Kriging, k::Kriging, new_x) +function calculate_liars(::KrigingBelieverUpperBound, tmp_k::Kriging, k::Kriging, new_x) new_y = k(new_x) + 3 * std_error_at_point(k, new_x) add_point!(tmp_k, new_x, new_y) end # Kriging Believer Lower Bound -function KBLB!(tmp_k::Kriging, k::Kriging, new_x) +function calculate_liars(::KrigingBelieverLowerBound, tmp_k::Kriging, k::Kriging, new_x) new_y = k(new_x) - 3 * std_error_at_point(k, new_x) add_point!(tmp_k, new_x, new_y) end \ No newline at end of file From 1b27538bce657410e4e41e87275deef1eeac58a2 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 20:31:49 -0500 Subject: [PATCH 12/16] Fixed indexing mistake --- src/Optimization.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/Optimization.jl b/src/Optimization.jl index 65e8fc77f..1130dd473 100755 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -388,8 +388,6 @@ function potential_optimal_points(::SRBF, strategy, lb, ub, surr::AbstractSurrog num_new_samples = 500) scale = 0.2 - success = 0 - failure = 0 w_range = [0.3, 0.5, 0.7, 0.95] w_cycle = Iterators.cycle(w_range) @@ -397,8 +395,6 @@ function potential_optimal_points(::SRBF, strategy, lb, ub, surr::AbstractSurrog #Vector containing size in each direction box_size = lb - ub - success = 0 - failures = 0 dtol = 1e-3 * norm(ub - lb) d = length(surr.x) incumbent_x = surr.x[argmin(surr.y)] @@ -452,7 +448,7 @@ function potential_optimal_points(::SRBF, strategy, lb, ub, surr::AbstractSurrog while new_addition < n_parallel #find minimum - @inbounds for r in 1:num_new_samples + @inbounds for r in eachindex(evaluation_of_merit_function) evaluation_of_merit_function[r] = merit_function(new_sample[r], w, tmp_surr, s_max, s_min, d_max, d_min, box_size) From 837002d544d4472d9d1d4aa198a772120efab455 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 21:34:17 -0500 Subject: [PATCH 13/16] Added documentation to pages.jl --- docs/pages.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/pages.jl b/docs/pages.jl index e870f3b69..ce76c3fdb 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -16,6 +16,7 @@ pages = ["index.md" "Gradient Enhanced Kriging" => "gek.md", "GEKPLS" => "gekpls.md", "MOE" => "moe.md", + "Parallel Optimization" => "parallel.md" ] "User guide" => [ "Samples" => "samples.md", From 7c2398428f5e725d1035352cf8c7670f29fd6a11 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 21:48:46 -0500 Subject: [PATCH 14/16] Updated documentation --- docs/src/parallel.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/parallel.md b/docs/src/parallel.md index 99d0004df..994d3ea32 100755 --- a/docs/src/parallel.md +++ b/docs/src/parallel.md @@ -4,12 +4,12 @@ There are some situations where it can be beneficial to run multiple optimizatio ## Ask-Tell Interface -To enable parallel optimization, we make use of an Ask-Tell interface. The user will construct the initial surrogate model the same way as for non-parallel surrogate models, but instead of using `surrogate_optimize`, the user will use `Ask`. This will return the coordinates of points that the optimizer has determined are most useful to evaluate next. How the user evaluates these points is up to them. The Ask-Tell interface requires more manual control than `surrogate_optimize`, but it allows for more flexibility. After the point has been evaluated, the user will *tell* the surrogate model the new points with the `add_point!` function. +To enable parallel optimization, we make use of an Ask-Tell interface. The user will construct the initial surrogate model the same way as for non-parallel surrogate models, but instead of using `surrogate_optimize`, the user will use `potential_optimal_points`. This will return the coordinates of points that the optimizer has determined are most useful to evaluate next. How the user evaluates these points is up to them. The Ask-Tell interface requires more manual control than `surrogate_optimize`, but it allows for more flexibility. After the point has been evaluated, the user will *tell* the surrogate model the new points with the `add_point!` function. ## Virtual Points -To ensure that points of interest returned by `Ask` are sufficiently far from each other, the function makes use of *virtual points*. They are used as follows: -1. `Ask` is told to return `n` points. +To ensure that points of interest returned by `potential_optimal_points` are sufficiently far from each other, the function makes use of *virtual points*. They are used as follows: +1. `potential_optimal_points` is told to return `n` points. 2. The point with the highest merit function value is selected. 3. This point is now treated as a virtual point and is assigned a temporary value that changes the landscape of the merit function. How the the temporary value is chosen depends on the strategy used. (see below) 4. The point with the new highest merit is selected. @@ -50,7 +50,7 @@ y = f.(x) my_k = Kriging(x, y, lb, ub) for _ in 1:10 - new_x, eis = Ask(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) + new_x, eis = potential_optimal_points(EI(), lb, ub, my_k, SobolSample(), 3, CLmean!) add_point!(my_k, new_x, f.(new_x)) end ``` From bcc5451102e2e5ed5db7952212ba735e312ef060 Mon Sep 17 00:00:00 2001 From: Dreycen Foiles Date: Tue, 19 Sep 2023 22:44:11 -0500 Subject: [PATCH 15/16] Fix bug in parallel EI --- src/Optimization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Optimization.jl b/src/Optimization.jl index 1130dd473..81e8e3523 100755 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -828,7 +828,7 @@ function potential_optimal_points(::EI, strategy, lb, ub, krig, sample_type::Sam index_max = argmax(evaluations) x_new = new_sample[index_max] # x point which maximized EI y_new = maximum(evaluations) # EI at the new point - diff_x = [abs.(prev_point .- x_new) for prev_point in tmp_krig.x] + diff_x = [norm(prev_point .- x_new) for prev_point in tmp_krig.x] bit_x = [diff_x_point .> dtol for diff_x_point in diff_x] #new_min_x has to have some distance from tmp_krig.x if false in bit_x From ea4eef23c7ece2a210d15050d21805c4aa32f55b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 09:26:29 -0400 Subject: [PATCH 16/16] Update docs/src/parallel.md --- docs/src/parallel.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/parallel.md b/docs/src/parallel.md index 994d3ea32..2388e1eec 100755 --- a/docs/src/parallel.md +++ b/docs/src/parallel.md @@ -38,7 +38,7 @@ In general, CLmin and KBLB tend to favor exploitation while CLmax and KBUB tend ## Examples -```jl +```@example using Surrogates lb = 0.0