diff --git a/config_model.yml b/config_model.yml
index 9323bfcbb1994f857f436e6fe6cd1190a090b8c3..97953475d3d9b7395064cededd06a0f7cf643bb6 100644
--- a/config_model.yml
+++ b/config_model.yml
@@ -37,8 +37,8 @@ siting_params:
       no_iterations: 2000
       no_epochs: 500
       initial_temp: 100.
-      no_runs: 20
-      algorithm: 'MIR' # 'SGH', 'RS'
+      no_runs: 100
+      algorithm: 'SGH' # 'SGH', 'RS'
     GRED_DET:
       set: False
       no_runs: 10
@@ -62,7 +62,7 @@ siting_params:
       no_years: 10
       no_years_train: 5
       no_years_test: 5
-      no_experiments: 10
+      no_experiments: 50
       no_runs_per_experiment: 10
       criterion: 'max'
       algorithm: 'custom' # 'k_fold'
diff --git a/src/jl/SitingHeuristics.jl b/src/jl/SitingHeuristics.jl
index af2a2e3e8b977ea2e269b8818042226408312574..adc1318a398cdd6911737ca7838c60c1b41d1d92 100644
--- a/src/jl/SitingHeuristics.jl
+++ b/src/jl/SitingHeuristics.jl
@@ -21,7 +21,7 @@ function myunpickle(filename)
 
 end
 
-function main_MIRSA(index_dict, deployment_dict, D, c, N, I, E, T_init, R, run, p, data_path)
+function main_MIRSA(index_dict, deployment_dict, D, c, N, I, E, T_init, R, run, p)
 
   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)])
@@ -34,8 +34,7 @@ function main_MIRSA(index_dict, deployment_dict, D, c, N, I, E, T_init, R, run,
   T_init = convert(Float64, T_init)
   R = convert(Int64, R)
   run = string(run)
-  p = string(p)
-  data_path = string(data_path)
+  p = convert(Float64, p)
   legacy_index = Vector{Int64}(undef, 0)
 
   W, L = size(D)
@@ -74,16 +73,17 @@ function main_MIRSA(index_dict, deployment_dict, D, c, N, I, E, T_init, R, run,
 
   elseif run == "SGH"
 
-    x_sol, LB_sol, obj_sol = Array{Float64, 2}(undef, R, L), Array{Float64, 1}(undef, R), Array{Float64, 2}(undef, R, I)
-
-    # Read LB_init from file
-    LB_init = myunpickle(joinpath(data_path, "objective_vector.p"))
-    # Read x_init from file
-    x_init = myunpickle(joinpath(data_path, "solution_matrix.p"))
+    R_init = 10
+    x_init, LB_init = Array{Float64, 2}(undef, R_init, L), Array{Float64, 1}(undef, R_init)
+    for r = 1:R_init
+      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, :]
+    println("Initial solution computed. Starting LS.")
 
+    x_sol, LB_sol, obj_sol = Array{Float64, 2}(undef, R, L), Array{Float64, 1}(undef, R), Array{Float64, 2}(undef, R, I)
     for r = 1:R
       x_sol[r, :], LB_sol[r], obj_sol[r, :] = simulated_annealing_local_search(D, c, n, N, I, E, x_init_best, T_init, legacy_index)
     end
diff --git a/src/main.py b/src/main.py
index 397855515dddc72e7dfc80a734691afe1b052551..43061f4d2359630693153ca0f535bffb2480d7d2 100644
--- a/src/main.py
+++ b/src/main.py
@@ -26,18 +26,11 @@ def parse_args():
     parser.add_argument('--run_RAND', type=bool, default=False)
     parser.add_argument('--run_CROSS', type=bool, default=False)
     parser.add_argument('--LS_init_algorithm', type=str, default=None)
-    parser.add_argument('--init_sol_folder', type=str, default=None)
 
     parsed_args = vars(parser.parse_args())
 
     return parsed_args
 
-
-def single_true(iterable):
-    i = iter(iterable)
-    return any(i) and not any(i)
-
-
 if __name__ == '__main__':
 
     args = parse_args()
@@ -99,10 +92,6 @@ if __name__ == '__main__':
 
     c = args['c']
 
-    if not single_true([args['run_BB'], args['run_LS'], args['run_GRED_DET'], args['run_GRED_STO'],
-                        args['run_RAND'], args['run_CROSS']]):
-        raise ValueError(' More than one run selected in the argparser.')
-
     if siting_parameters['solution_method']['BB']['set']:
 
         custom_log(' BB chosen to solve the IP.')
@@ -139,8 +128,6 @@ if __name__ == '__main__':
         params = siting_parameters['solution_method']['LS']
 
         jl_dict = generate_jl_input(deployment_dict, site_coordinates)
-        path_to_sol = args['init_sol_folder'] + str(args['c']) + '_GRED_STGH_p' + str(args['p'])
-        path_to_init_sol_folder = join(data_path, 'output', path_to_sol)
 
         import julia
         j = julia.Julia(compiled_modules=False)
@@ -152,12 +139,11 @@ if __name__ == '__main__':
                                                              criticality_data, c, params['neighborhood'],
                                                              params['no_iterations'], params['no_epochs'],
                                                              params['initial_temp'], params['no_runs'],
-                                                             args['LS_init_algorithm'],
-                                                             args['p'], path_to_init_sol_folder)
+                                                             args['LS_init_algorithm'], args['p'])
         end = time.time()
         print(f"Average CPU time for c={c}: {round((end-start)/params['no_runs'], 1)} s")
 
-        output_folder = init_folder(model_parameters, c, suffix=f"_LS_{args['LS_init_algorithm']}_bis")
+        output_folder = init_folder(model_parameters, c, suffix=f"_LS_{args['LS_init_algorithm']}")
 
         with open(join(output_folder, 'config_model.yaml'), 'w') as outfile:
             yaml.dump(model_parameters, outfile, default_flow_style=False, sort_keys=False)