Skip to content

Commit

Permalink
Merge a44ca7c into 3bd8db5
Browse files Browse the repository at this point in the history
  • Loading branch information
@ludoro
ludoro committed Aug 20, 2020
2 parents 3bd8db5 + a44ca7c commit 9567392
Show file tree
Hide file tree
Showing 2 changed files with 263 additions and 5 deletions.
255 changes: 252 additions & 3 deletions src/Optimization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -970,7 +970,7 @@ SOP Surrogate optimization method, following closely the following papers:
- Multiobjective Optimization Using Evolutionary Algorithms by Kalyan Deb
#Suggested number of new_samples = min(500*d,5000)
"""
function surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSOP::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100)
function surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSOP::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=min(500*1,5000))
d = length(lb)
N_fail = 3
N_tenure = 5
Expand Down Expand Up @@ -1111,15 +1111,15 @@ function surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSO
N_failures[i] = N_failures_global[i]
end

f_1 = new_points[i,1]
f_1 = obj(new_points[i,1])
f_2 = obj2_1D(f_1,surrSOP.x)

l = length(Fronts[1])
Pareto_set = zeros(eltype(surrSOP.x[1]),l,2)

for j = 1:l
val = obj2_1D(Fronts[1][j],surrSOP.x)
Pareto_set[j,1] = Fronts[1][j]
Pareto_set[j,1] = obj(Fronts[1][j])
Pareto_set[j,2] = val
end
if (Hypervolume_Pareto_improving(f_1,f_2,Pareto_set)<tau)
Expand All @@ -1142,3 +1142,252 @@ function surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSO
index = argmin(surrSOP.y)
return (surrSOP.x[index],surrSOP.y[index])
end


function obj2_ND(value,points)
min = +Inf
my_p = filter(x->norm(x.-value)>10^-6,points)
for i = 1:length(my_p)
new_val = norm(my_p[i].-value)
if new_val < min
min = new_val
end
end
return min
end

function I_tier_ranking_ND(P,surrSOPD::AbstractSurrogate)
#obj1 = objective_function
#obj2 = obj2_1D
Fronts = Dict{Int,Array{eltype(surrSOPD.x),1}}()
i = 1
while true
F = Array{eltype(surrSOPD.x),1}()
j = 1
for p in P
n_p = 0
k = 1
for q in P
#I use equality with floats because p and q are in surrSOP.x
#for sure at this stage
p_index = j
q_index = k
val1_p = surrSOPD.y[p_index]
val2_p = obj2_ND(p,P)
val1_q = surrSOPD.y[q_index]
val2_q = obj2_ND(q,P)
p_dominates_q = (val1_p < val1_q || abs(val1_p-val1_q) <= 10^-5) &&
(val2_p < val2_q || abs(val2_p-val2_q) <= 10^-5) &&
((val1_p < val1_q) || (val2_p < val2_q))

q_dominates_p = (val1_p < val1_q || abs(val1_p-val1_q) < 10^-5) &&
(val2_p < val2_q || abs(val2_p-val2_q) < 10^-5) &&
((val1_p < val1_q) || (val2_p < val2_q))
if q_dominates_p
n_p += 1
end
k = k + 1
end
if n_p == 0
# no individual dominates p
push!(F,p)
end
j = j + 1
end
if length(F) > 0
Fronts[i] = F
P = setdiff(P,F)
i = i + 1
else
return Fronts
end
end
return F
end

function II_tier_ranking_ND(D::Dict,srgD::AbstractSurrogate)
for i = 1:length(D)
pos = []
yn = []
for j = 1:length(D[i])
push!(pos,findall(e->e==D[i][j],srgD.x))
push!(yn,srgD.y[pos[j]])
end
D[i] = D[i][sortperm(D[i])]
end
return D
end

function surrogate_optimize(obj::Function,sopd::SOP,lb,ub,surrSOPD::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=min(500*length(lb),5000))
d = length(lb)
N_fail = 3
N_tenure = 5
tau = 10^-5
num_P = sopd.p
centers_global = surrSOPD.x
r_centers_global = 0.2*norm(ub.-lb)*ones(length(surrSOPD.x))
N_failures_global = zeros(length(surrSOPD.x))
tabu = []
N_tenures_tabu = []
for k = 1:maxiters
N_tenures_tabu .+= 1
#deleting points that have been in tabu for too long
del = N_tenures_tabu .> N_tenure

if length(del) > 0
for i = 1:length(del)
if del[i]
del[i] = i
end
end
deleteat!(N_tenures_tabu,del)
deleteat!(tabu,del)
end


##### P CENTERS ######
C = Array{eltype(surrSOPD.x),1}()

#S(x) set of points already evaluated
#Rank points in S with:
#1) Non dominated sorting
Fronts_I = I_tier_ranking_ND(centers_global,surrSOPD)
#2) Second tier ranking
Fronts = II_tier_ranking_ND(Fronts_I,surrSOPD)
ranked_list = Array{eltype(surrSOPD.x),1}()
for i = 1:length(Fronts)
for j = 1:length(Fronts[i])
push!(ranked_list,Fronts[i][j])
end
end

centers_full = 0
i = 1
while i <= length(ranked_list) && centers_full == 0
flag = 0
for j = 1:length(ranked_list)
for m = 1:length(tabu)
if norm(ranked_list[j].-tabu[m]) < tau
flag = 1
end
end
for l = 1:length(centers_global)
if norm(ranked_list[j].-centers_global[l]) < tau
flag = 1
end
end
end
if flag == 1
skip
else
push!(C,ranked_list[i])
if length(C) == num_P
centers_full = 1
end
end
i = i + 1
end

#I examined all the points in the ranked list but num_selected < num_p
#I just iterate again using only radius rule
if length(C) < num_P
i = 1
while i <= length(ranked_list) && centers_full == 0
flag = 0
for j = 1:length(ranked_list)
for m = 1:length(centers_global)
if norm(centers_global[j] .- ranked_list[m]) < tau
flag = 1
end
end
end
if flag == 1
skip
else
push!(C,ranked_list[i])
if length(C) == num_P
centers_full = 1
end
end
i = i + 1
end
end

#If I still have num_selected < num_P, I double down on some centers iteratively
if length(C) < num_P
i = 1
while i <= length(ranked_list)
push!(C,ranked_list[i])
if length(C) == num_P
centers_full = 1
end
i = i + 1
end
end

#Here I have selected C = [(1.0,2.0),(3.0,4.0),.....] containing the centers
r_centers = 0.2*norm(ub.-lb)*ones(num_P)
N_failures = zeros(num_P)
#2.3 Candidate search
new_points_x = Array{eltype(surrSOPD.x),1}()
new_points_y = zeros(eltype(surrSOPD.y[1]),num_P)
for i = 1:num_P
N_candidates = zeros(eltype(surrSOPD.x[1]),num_new_samples,d)
#Using phi(n) just like DYCORS, merit function = surrogate
#Like in DYCORS, I_perturb = 1 always
evaluations = zeros(eltype(surrSOPD.y[1]),num_new_samples)
for j = 1:num_new_samples
for k = 1:d
a = lb[k] - C[i][k]
b = ub[k] - C[i][k]
N_candidates[j,k] = C[i][k] + rand(truncated(Normal(0,r_centers[i]), a,b))
end
evaluations[j] = surrSOPD(Tuple(N_candidates[j,:]))
end
x_best = Tuple(N_candidates[argmin(evaluations),:])
y_best = minimum(evaluations)
push!(new_points_x,x_best)
new_points_y[i] = y_best
end

#new_points[i] is splitted in new_points_x and new_points_y now contains:
#[x_1,y_1; x_2,y_2,...,x_{num_new_samples},y_{num_new_samples}]


#2.4 Adaptive learning and tabu archive
for i=1:num_P
if new_points_x[i] in centers_global
r_centers[i] = r_centers_global[i]
N_failures[i] = N_failures_global[i]
end

f_1 = obj(Tuple(new_points_x[i]))
f_2 = obj2_ND(f_1,surrSOPD.x)

l = length(Fronts[1])
Pareto_set = zeros(eltype(surrSOPD.x[1]),l,2)
for j = 1:l
val = obj2_ND(Fronts[1][j],surrSOPD.x)
Pareto_set[j,1] = obj(Tuple(Fronts[1][j]))
Pareto_set[j,2] = val
end
if (Hypervolume_Pareto_improving(f_1,f_2,Pareto_set)<tau)#check this
#failure
r_centers[i] = r_centers[i]/2
N_failures[i] += 1
if N_failures[i] > N_fail
push!(tabu,C[i])
push!(N_tenures_tabu,0)
end
else
#P_i is success
#Adaptive_learning
add_point!(surrSOPD,new_points_x[i],new_points_y[i])
push!(r_centers_global,r_centers[i])
push!(N_failures_global,N_failures[i])
end
end
end
index = argmin(surrSOPD.y)
return (surrSOPD.x[index],surrSOPD.y[index])
end
13 changes: 11 additions & 2 deletions test/optimization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ using Flux
using Flux: @epochs

#######SRBF############

##### 1D #####
lb = 0.0
ub = 15.0
Expand Down Expand Up @@ -220,7 +219,6 @@ p = [1.5,1.5]
theta = [2.0,2.0]
lb = [1.0,1.0]
ub = [6.0,6.0]
bounds = [lb,ub]


my_k_DYCORSN = Kriging(x,y,lb,ub)
Expand All @@ -243,3 +241,14 @@ ub = 6.0
num_centers = 2
my_k_SOP1 = Kriging(x,y,lb,ub,p=1.9)
surrogate_optimize(objective_function,SOP(num_centers),lb,ub,my_k_SOP1,SobolSample(),maxiters=60)
#ND
objective_function_ND = z -> 2*norm(z)+1
x = [(2.3,2.2),(1.4,1.5)]
y = objective_function_ND.(x)
p = [1.5,1.5]
theta = [2.0,2.0]
lb = [1.0,1.0]
ub = [6.0,6.0]
my_k_SOPND = Kriging(x,y,lb,ub)
num_centers = 2
surrogate_optimize(objective_function_ND,SOP(num_centers),lb,ub,my_k_SOPND,SobolSample(),maxiters=20)

0 comments on commit 9567392

Please sign in to comment.