SOP ND

src/Optimization.jl

@@ -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
@@ -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)
@@ -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

test/optimization.jl

@@ -4,7 +4,6 @@ using Flux
 using Flux: @epochs
 
 #######SRBF############
-
 ##### 1D #####
 lb = 0.0
 ub = 15.0
@@ -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)
@@ -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)