diff --git a/config_model.yml b/config_model.yml
index 8c5dbe22f4525d569ca72e8d3f8018ea30104695..8ef6219103fc88c9f6b975d62e9861af034801f3 100644
--- a/config_model.yml
+++ b/config_model.yml
@@ -33,13 +33,13 @@ siting_params:
threads: 36
MIRSA:
# Simulated Annealing with Local Search
- set: True
+ set: False
neighborhood: 1
no_iterations: 2000
no_epochs: 500
initial_temp: 200.
no_runs: 50
- algorithm: 'SALS' # 'GLS' (Greedy Local Search), 'RSSA', 'SALS', 'STGHSA'
+ algorithm: 'SALS' # 'GLS' (Greedy Local Search), 'RSSA', 'SALS', 'STGHSAS', 'STGHSAP'
GRED_DET:
set: False
no_runs: 10
@@ -57,3 +57,14 @@ siting_params:
no_epochs: 1000
no_runs: 200
algorithm: 'RS'
+ SGHLS:
+ # Stochastic greedy & simulated annealing
+ set: True
+ neighborhood: 1
+ no_iterations: 2000
+ no_epochs: 500
+ initial_temp: 200.
+ no_runs_LS: 50
+ algorithm: 'STGHSAP' # 'STGHSAP'
+ p: 27 # 10% of 10150 locations
+ no_runs_SGH: 200
diff --git a/src/jl/MCP_heuristics.jl b/src/jl/MCP_heuristics.jl
index eedd0fbcd675b9ee162f4d8951ccc69592333b72..cfc35ad4b5e89c73608be90923d390f95c5cf384 100644
--- a/src/jl/MCP_heuristics.jl
+++ b/src/jl/MCP_heuristics.jl
@@ -439,4 +439,50 @@ function stochastic_threshold_greedy_algorithm(D::Array{Float64,2}, c::Float64,
x_incumbent[ind_incumbent] .= 1.
return x_incumbent, obj_incumbent
+end
+
+function stochastic_threshold_greedy_algorithm_trajectory(D::Array{Float64,2}, c::Float64, n::Float64, p::Float64)
+ W, L = size(D)
+ s = convert(Int64, round((L/n)*p))
+ random_ind_set = Vector{Int64}(undef, s)
+ ind_compl_incumbent = [i for i in 1:L]
+ ind_incumbent = zeros(Int64, 0)
+ Dx_incumbent = zeros(Float64, W)
+ y_incumbent = zeros(Float64, W)
+ obj_incumbent = 0
+ obj_incumbent_vec = zeros(Float64, convert(Int64,n))
+ Dx_tmp = Vector{Float64}(undef, W)
+ y_tmp = Vector{Float64}(undef, W)
+ locations_added, threshold = 0, 0
+ @inbounds while locations_added < n
+ if locations_added < c
+ threshold = locations_added + 1
+ obj_candidate = 0
+ else
+ obj_candidate = obj_incumbent
+ end
+ ind_candidate_list = Vector{Int64}(undef, 0)
+ random_ind_set .= sample(ind_compl_incumbent, s, replace=false)
+ @inbounds for ind in random_ind_set
+ Dx_tmp .= Dx_incumbent .+ view(D, :, ind)
+ y_tmp .= Dx_tmp .>= threshold
+ obj_tmp = sum(y_tmp)
+ if obj_tmp > obj_candidate
+ ind_candidate_list = [ind]
+ obj_candidate = obj_tmp
+ elseif obj_tmp == obj_candidate
+ ind_candidate_list = union(ind, ind_candidate_list)
+ end
+ end
+ ind_candidate = sample(ind_candidate_list)
+ ind_incumbent = union(ind_incumbent, ind_candidate)
+ ind_compl_incumbent = setdiff(ind_compl_incumbent, ind_candidate)
+ Dx_incumbent .= Dx_incumbent .+ view(D, :, ind_candidate)
+ y_incumbent .= Dx_incumbent .>= threshold
+ obj_incumbent = sum(y_incumbent)
+ obj_incumbent_vec[locations_added+1] = obj_incumbent
+ locations_added += 1
+ end
+ return ind_incumbent, obj_incumbent
+
end
\ No newline at end of file
diff --git a/src/jl/SitingHeuristics.jl b/src/jl/SitingHeuristics.jl
index d36e740984d9a2fe11b3089fcf45108023d78753..d76cc1b4a9b83066d567574e02049306e5960899 100644
--- a/src/jl/SitingHeuristics.jl
+++ b/src/jl/SitingHeuristics.jl
@@ -62,6 +62,69 @@ function main_MIRSA(index_dict, deployment_dict, D, c, N, I, E, T_init, R, run)
end
+function main_SGHLS(index_dict, deployment_dict, D, c, N, I, E, T_init, p, R_SGH, R_LS, run)
+
+ index_dict = Dict([(convert(Int64, k), convert(Int64, index_dict[k])) for k in keys(index_dict)])
+ deployment_dict = Dict([(convert(Int64, k), convert(Int64, deployment_dict[k])) for k in keys(deployment_dict)])
+ n = convert(Float64, deployment_dict[1])
+ D = convert.(Float64, D)
+
+ c = convert(Float64, c)
+ N = convert(Int64, N)
+ I = convert(Int64, I)
+ E = convert(Int64, E)
+ T_init = convert(Float64, T_init)
+ p = convert(Float64, p)
+ R_SGH = convert(Int64, R_SGH)
+ R_LS = convert(Int64, R_LS)
+ run = string(run)
+
+ n_partitions = [deployment_dict[i] for i in 1:maximum(values(index_dict))]
+
+ W, L = size(D)
+
+ if run == "STGHSAS"
+
+ x_init, LB_init = Array{Float64, 2}(undef, R_SGH, L), Array{Float64, 1}(undef, R_SGH)
+ x_sol, LB_sol, obj_sol = Array{Float64, 2}(undef, R_LS, L), Array{Float64, 1}(undef, R_LS), Array{Float64, 2}(undef, R_LS, I)
+
+ for r = 1:R_SGH
+ x_init[r, :], LB_init[r] = stochastic_threshold_greedy_algorithm(D, c, n, p)
+ end
+
+ LB_init_best = argmax(LB_init)
+ x_init_best = x_init[LB_init_best, :]
+
+ for r = 1:R_LS
+ x_sol[r, :], LB_sol[r], obj_sol[r, :] = simulated_annealing_local_search_partition(D, c, n_partitions, N, I, E, x_init_best, T_init, index_dict)
+ end
+
+ elseif run == "STGHSAP"
+
+ x_sol, LB_sol, obj_sol = Array{Float64, 2}(undef, R_LS, L), Array{Float64, 1}(undef, R_LS), Array{Float64, 2}(undef, R_LS, I)
+
+ for r = 1:R_LS
+
+ x_init, LB_init = Array{Float64, 2}(undef, 1, L), Array{Float64, 1}(undef, 1)
+ x_init[1, :], LB_init[1] = stochastic_threshold_greedy_algorithm(D, c, n, p)
+
+ x_init_best = x_init[1, :]
+
+ x_sol[r, :], LB_sol[r], obj_sol[r, :] = simulated_annealing_local_search_partition(D, c, n_partitions, N, I, E, x_init_best, T_init, index_dict)
+
+ end
+
+ else
+ println("No such run available.")
+ throw(ArgumentError)
+ end
+
+ return x_sol, LB_sol, obj_sol
+
+end
+
+
+
function main_GRED(deployment_dict, D, c, R, p, run)
deployment_dict = Dict([(convert(Int64, k), convert(Int64, deployment_dict[k])) for k in keys(deployment_dict)])
diff --git a/src/main.py b/src/main.py
index e223b0780aac4fda719ca4140114e79dc3397c97..1704c6c9b37a6adb104367a6b4258ae4e230d7eb 100644
--- a/src/main.py
+++ b/src/main.py
@@ -22,6 +22,7 @@ def parse_args():
parser.add_argument('--run_GRED_DET', type=bool, default=False)
parser.add_argument('--run_GRED_STO', type=bool, default=False)
parser.add_argument('--run_RAND', type=bool, default=False)
+ parser.add_argument('--run_SGHLS', type=bool, default=False)
parsed_args = vars(parser.parse_args())
@@ -88,10 +89,12 @@ if __name__ == '__main__':
siting_parameters['solution_method']['GRED_DET']['set'] = args['run_GRED_DET']
siting_parameters['solution_method']['GRED_STO']['set'] = args['run_GRED_STO']
siting_parameters['solution_method']['RAND']['set'] = args['run_RAND']
+ siting_parameters['solution_method']['SGHLS']['set'] = args['run_SGHLS']
c = args['c']
- if not single_true([args['run_BB'], args['run_MIRSA'], args['run_GRED_DET'], args['run_GRED_STO'], args['run_RAND']]):
+ if not single_true([args['run_BB'], args['run_MIRSA'], args['run_GRED_DET'], args['run_GRED_STO'],
+ args['run_RAND'], args['run_SGHLS']]):
raise ValueError(' More than one run selected in the argparser.')
if siting_parameters['solution_method']['BB']['set']:
@@ -167,8 +170,7 @@ if __name__ == '__main__':
j = julia.Julia(compiled_modules=False)
from julia import Main
Main.include("jl/SitingHeuristics.jl")
-
-
+
jl_selected, jl_objective = Main.main_RAND(jl_dict['deployment_dict'], criticality_data,
c, params['no_iterations'], params['no_runs'],
params['algorithm'])
@@ -224,5 +226,33 @@ if __name__ == '__main__':
pickle.dump(jl_selected, open(join(output_folder, 'solution_matrix.p'), 'wb'))
pickle.dump(jl_objective, open(join(output_folder, 'objective_vector.p'), 'wb'))
+ elif siting_parameters['solution_method']['SGHLS']['set']:
+
+ params = siting_parameters['solution_method']['SGHLS']
+ custom_log(f" SHGLS_{params['algorithm']} chosen to solve the IP. Opening a Julia instance.")
+
+ jl_dict = generate_jl_input(deployment_dict, site_coordinates)
+
+ import julia
+ j = julia.Julia(compiled_modules=False)
+ from julia import Main
+ Main.include("jl/SitingHeuristics.jl")
+
+ start = time.time()
+ jl_selected, jl_objective, jl_traj = Main.main_SGHLS(jl_dict['index_dict'], jl_dict['deployment_dict'],
+ criticality_data, c, params['neighborhood'],
+ params['no_iterations'], params['no_epochs'],
+ params['initial_temp'], params['p'],
+ params['no_runs_SGH'], params['no_runs_LS'],
+ params['algorithm'])
+ end = time.time()
+ print(f"Average CPU time for c={c}: {round((end-start)/params['no_runs_LS'], 1)} s")
+
+ output_folder = init_folder(model_parameters, c, suffix=f"_SGHLS_{params['algorithm']}_p{params['p']}")
+
+ pickle.dump(jl_selected, open(join(output_folder, 'solution_matrix.p'), 'wb'))
+ pickle.dump(jl_objective, open(join(output_folder, 'objective_vector.p'), 'wb'))
+ pickle.dump(jl_traj, open(join(output_folder, 'trajectory_matrix.p'), 'wb'))
+
else:
raise ValueError(' This solution method is not available. ')