Merge pull request #3 from luinardi/prior_guided_optimization

Added prior-guided optimization approach.

example_outputs/synthetic_functions/prior_currinexp/output_prior_currinexp.txt

@@ -0,0 +1,92 @@
+Design of experiment phase, number of doe samples = 3 .......
+x1,x2,Value,Timestamp
+0.0711680251749421,0.918860499459033,4.02779694878494,6950
+0.08374201043118659,0.08391745617300503,10.42748160722958,6950
+0.13139443915285315,0.5365545892167386,7.692592659010549,6950
+
+
+End of doe/resume phase, the number of configuration runs is: 3
+
+Starting optimization iteration 1
+x1,x2,Value,Timestamp
+0.4509563903201488,0.7069863095345909,6.117712074031937,6965
+
+Starting optimization iteration 2
+x1,x2,Value,Timestamp
+0.36392434241928656,1.0,5.035876129222752,15253
+
+Starting optimization iteration 3
+x1,x2,Value,Timestamp
+4.062834327761441e-06,1.0,1.1805742753340724,23483
+
+Starting optimization iteration 4
+x1,x2,Value,Timestamp
+3.05653254848465e-07,1.0,1.1804205284412628,30936
+
+Starting optimization iteration 5
+x1,x2,Value,Timestamp
+1.0265688688718155e-05,1.0,1.1808281011384816,38248
+
+Starting optimization iteration 6
+x1,x2,Value,Timestamp
+1.2548908632819513e-06,1.0,1.1804593720139331,46522
+
+Starting optimization iteration 7
+x1,x2,Value,Timestamp
+1.778651187016753e-05,1.0,1.1811358592545058,54803
+
+Starting optimization iteration 8
+x1,x2,Value,Timestamp
+0.3398424549396122,0.362590602351231,9.737453904476649,54820
+
+Starting optimization iteration 9
+x1,x2,Value,Timestamp
+0.0024892400780651457,1.0,1.282252625178526,63329
+
+Starting optimization iteration 10
+x1,x2,Value,Timestamp
+0.0006171945147863797,1.0,1.2056637700336026,71662
+
+Starting optimization iteration 11
+x1,x2,Value,Timestamp
+1.3516466735036529e-05,1.0,1.180961125588085,79518
+
+Starting optimization iteration 12
+x1,x2,Value,Timestamp
+5.50726757179193e-05,1.0,1.1826616383188775,86852
+
+Starting optimization iteration 13
+x1,x2,Value,Timestamp
+6.715085215939413e-07,1.0,1.1804354995355162,94167
+
+Starting optimization iteration 14
+x1,x2,Value,Timestamp
+2.746991426485545e-06,1.0,1.180520429977757,101582
+
+Starting optimization iteration 15
+x1,x2,Value,Timestamp
+9.88987177950969e-05,1.0,1.1844550326094196,108838
+
+Starting optimization iteration 16
+x1,x2,Value,Timestamp
+6.517528293725511e-07,1.0,1.1804346911165693,116122
+
+Starting optimization iteration 17
+x1,x2,Value,Timestamp
+0.014024969093860385,1.0,1.7515590044049345,123734
+
+Starting optimization iteration 18
+x1,x2,Value,Timestamp
+1.134778693051246e-06,1.0,1.1804544569265605,131117
+
+Starting optimization iteration 19
+x1,x2,Value,Timestamp
+0.003972875152055387,1.0,1.3429146725521282,138537
+
+Starting optimization iteration 20
+x1,x2,Value,Timestamp
+2.5976861786358505e-05,1.0,1.181471014894109,145806
+
+End of Prior Optimization
+### End of the hypermapper script.
+End of CurrinExp.

example_outputs/synthetic_functions/prior_currinexp/prior_currinexp_output_samples.csv

@@ -0,0 +1,24 @@
+x1,x2,Value,Timestamp
+0.0711680251749421,0.918860499459033,4.02779694878494,6950
+0.08374201043118659,0.08391745617300503,10.42748160722958,6950
+0.13139443915285315,0.5365545892167386,7.692592659010549,6950
+0.4509563903201488,0.7069863095345909,6.117712074031937,6965
+0.36392434241928656,1.0,5.035876129222752,15253
+4.062834327761441e-06,1.0,1.1805742753340724,23483
+3.05653254848465e-07,1.0,1.1804205284412628,30936
+1.0265688688718155e-05,1.0,1.1808281011384816,38248
+1.2548908632819513e-06,1.0,1.1804593720139331,46522
+1.778651187016753e-05,1.0,1.1811358592545058,54803
+0.3398424549396122,0.362590602351231,9.737453904476649,54820
+0.0024892400780651457,1.0,1.282252625178526,63329
+0.0006171945147863797,1.0,1.2056637700336026,71662
+1.3516466735036529e-05,1.0,1.180961125588085,79518
+5.50726757179193e-05,1.0,1.1826616383188775,86852
+6.715085215939413e-07,1.0,1.1804354995355162,94167
+2.746991426485545e-06,1.0,1.180520429977757,101582
+9.88987177950969e-05,1.0,1.1844550326094196,108838
+6.517528293725511e-07,1.0,1.1804346911165693,116122
+0.014024969093860385,1.0,1.7515590044049345,123734
+1.134778693051246e-06,1.0,1.1804544569265605,131117
+0.003972875152055387,1.0,1.3429146725521282,138537
+2.5976861786358505e-05,1.0,1.181471014894109,145806

example_scenarios/synthetic/prior_currinexp/prior_currinexp.py

@@ -0,0 +1,29 @@
+#!/usr/bin/python
+import math, sys
+from subprocess import Popen, PIPE
+sys.path.append('scripts')
+import hypermapper
+
+def currin_exp_function(X):
+    """
+    Compute the CurrinExp function.
+    The value is computed as defined in https://www.sfu.ca/~ssurjano/curretal88exp.html
+    :param X: dictionary containing the input points.
+    :return: the value of the CurrinExp function.
+    """
+    x1 = X['x1']
+    x2 = X['x2']
+    factor1 = 1 - math.exp(-1/(2*x2))
+    factor2 = 2300*x1*x1*x1 + 1900*x1*x1 + 2092*x1 + 60
+    factor3 = 100*x1*x1*x1 + 500*x1*x1 + 4*x1 + 20
+    y_value = factor1*factor2/factor3
+
+    return y_value
+
+def main():
+    parameters_file = "example_scenarios/synthetic/prior_currinexp/prior_currinexp_scenario.json"
+    hypermapper.optimize(parameters_file, currin_exp_function)
+    print("End of CurrinExp.")
+
+if __name__ == "__main__":
+    main()

example_scenarios/synthetic/prior_currinexp/prior_currinexp_scenario.json

@@ -0,0 +1,22 @@
+{
+    "application_name": "prior_currinexp",
+    "optimization_objectives": ["Value"],
+    "optimization_iterations": 20,
+    "optimization_method": "prior_guided_optimization",
+    "design_of_experiment": {
+        "number_of_samples": 3
+    },
+    "input_parameters" : {
+        "x1": {
+            "parameter_type" : "real",
+            "values" : [0, 1],
+            "prior": "decay"
+        },
+        "x2": {
+            "parameter_type" : "real",
+            "values" : [0.00001, 1],
+            "prior": "exponential"
+
+        }
+    }
+}

scripts/hypermapper.py

@@ -1,6 +1,7 @@
 import sys
 import random_scalarizations
 import local_search
+import prior_optimization
 import compute_pareto
 import plot_dse
 import json
@@ -54,10 +55,12 @@ def optimize(parameters_file, black_box_function=None):
 
     optimization_method = config["optimization_method"]
 
-    if optimization_method == "random_scalarizations":
+    if (optimization_method == "random_scalarizations") or (optimization_method == "bayesian_optimization"):
         random_scalarizations.main(config, black_box_function=black_box_function)
     elif optimization_method == "local_search":
         local_search.main(config, black_box_function=black_box_function)
+    elif optimization_method == 'prior_guided_optimization':
+        prior_optimization.main(config, black_box_function)
     else:
         print("Unrecognized optimization method:", optimization_method)
         raise SystemExit

scripts/local_search.py

@@ -286,39 +286,62 @@ def local_search(
     if str_data not in fast_addressing_of_data_array:
         tmp_fast_addressing_of_data_array[str_data] = 1
         if local_search_random_points - 1 > 0:
-            configurations = [default_configuration] + param_space.random_sample_configurations_without_repetitions(tmp_fast_addressing_of_data_array, local_search_random_points-1)
+            uniform_configurations = [default_configuration] + param_space.random_sample_configurations_without_repetitions(tmp_fast_addressing_of_data_array, local_search_random_points-1, use_priors=False)
     else:
-        configurations = param_space.random_sample_configurations_without_repetitions(tmp_fast_addressing_of_data_array, local_search_random_points)
+        uniform_configurations = [default_configuration] + param_space.random_sample_configurations_without_repetitions(tmp_fast_addressing_of_data_array, local_search_random_points-1, use_priors=False)
+
+    prior_configurations = param_space.random_sample_configurations_without_repetitions(tmp_fast_addressing_of_data_array, local_search_random_points, use_priors=True) # will be uniform random if no prior
 
     # Passing the dictionary with ** expands the key-value pairs into function parameters
-    function_values, feasibility_indicators = optimization_function(configurations=configurations, **optimization_function_parameters)
+    function_values_uniform, feasibility_indicators_uniform = optimization_function(configurations=uniform_configurations, **optimization_function_parameters)
+    function_values_prior, feasibility_indicators_prior = optimization_function(configurations=prior_configurations, **optimization_function_parameters)
 
     # This will concatenate the entire neighbors array if all configurations were evaluated
     # but only the evaluated configurations if we reached the budget and did not evaluate all
-    function_values_size = len(function_values)
-    new_data_array = concatenate_list_of_dictionaries(configurations[:function_values_size])
-    new_data_array[scalarization_key] = function_values
+    function_values_uniform_size = len(function_values_uniform)
+    new_data_array_uniform = concatenate_list_of_dictionaries(uniform_configurations[:function_values_uniform_size])
+    new_data_array_uniform[scalarization_key] = function_values_uniform
+
+    function_values_prior_size = len(function_values_prior)
+    new_data_array_prior   = concatenate_list_of_dictionaries(prior_configurations[:function_values_prior_size])
+    new_data_array_prior[scalarization_key] = function_values_prior
+
     if enable_feasible_predictor:
-        new_data_array[feasible_parameter] = feasibility_indicators
+        new_data_array_uniform[feasible_parameter] = feasibility_indicators_uniform
+        new_data_array_prior[feasible_parameter] = feasibility_indicators_prior
 
+    new_data_array = concatenate_data_dictionaries(new_data_array_uniform, new_data_array_prior)
     data_array = concatenate_data_dictionaries(data_array, new_data_array)
 
     # If some configurations were not evaluated, we reached the budget and must stop
-    if function_values_size < len(configurations):
+    if (function_values_uniform_size < len(uniform_configurations)) or (function_values_prior_size < len(prior_configurations)):
         sys.stdout.write_to_logfile("Out of budget, not all configurations were evaluated, stopping local search\n")
         end_of_search = True
 
     if enable_feasible_predictor:
-        local_search_configurations = get_min_feasible_configurations(
-                                                                    data_array,
-                                                                    local_search_starting_points,
-                                                                    scalarization_key,
-                                                                    feasible_parameter)
+        local_search_configurations_uniform = get_min_feasible_configurations(
+                                                                            new_data_array_uniform,
+                                                                            local_search_starting_points,
+                                                                            scalarization_key,
+                                                                            feasible_parameter)
+        local_search_configurations_prior = get_min_feasible_configurations(
+                                                                            new_data_array_prior,
+                                                                            local_search_starting_points,
+                                                                            scalarization_key,
+                                                                            feasible_parameter)
     else:
-        local_search_configurations = get_min_configurations(
-                                                            data_array,
-                                                            local_search_starting_points,
-                                                            scalarization_key)
+        local_search_configurations_uniform = get_min_configurations(
+                                                                    new_data_array_uniform,
+                                                                    local_search_starting_points,
+                                                                    scalarization_key)
+        local_search_configurations_prior = get_min_configurations(
+                                                                new_data_array_prior,
+                                                                local_search_starting_points,
+                                                                scalarization_key)
+
+    local_search_configurations = concatenate_data_dictionaries(
+                                                            local_search_configurations_uniform,
+                                                            local_search_configurations_prior)
 
     if previous_points is not None:
         concatenation_keys = input_params + [scalarization_key]
@@ -335,6 +358,7 @@ def local_search(
                                                     local_search_starting_points,
                                                     scalarization_key)
 
+
         local_search_configurations = concatenate_data_dictionaries(local_search_configurations, best_previous, concatenation_keys)
         data_array = concatenate_data_dictionaries(data_array, previous_points, concatenation_keys)
 

scripts/models.py

@@ -55,7 +55,7 @@ def generate_multi_output_regression_model(data_array
     if len (Y_test) == 0:
         Y_test = Y[:]
 
-    regressor = customRegressor(n_estimators=n_estimators, bootstrap=False, min_samples_split=5, max_features = max_features, n_jobs=1)
+    regressor = customRegressor(n_estimators=n_estimators, max_features = max_features, n_jobs=1, bootstrap=False, min_samples_split=5) 
     regressor.fit(X_train, y_train)
 
     if print_importances:
@@ -71,7 +71,7 @@ def generate_mono_output_regression_models(data_array
                                , Ycols
                                , learn_ratio
                                , debug=False
-                               , model_type="gaussian_process"
+                               , model_type="random_forest"
                                , number_of_cpus=0
                                , print_importances=False
                                , **model_kwargs):
@@ -316,7 +316,6 @@ def compute_parameter_importance(model, input_params, param_space):
             feature_idx += 1
     return parameter_importances
 
-
 def parallel_model_prediction(model, bufferx, param_space, debug=False, number_of_cpus=0):
     """
     This function explicitly parallelize the prediction of the Random Forest model.
@@ -651,6 +650,36 @@ def transform_rf_using_uniform_splits(regression_models, data_array, param_space
                 tree.tree_.threshold[node_idx] = new_split
     return regression_models
 
+def compute_probability_from_model(
+                                model_means,
+                                model_stds,
+                                param_space,
+                                objective_weights,
+                                threshold,
+                                compute_bad=True):
+    """
+    Compute the probability of a configuration being good or bad according to the model.
+    :param model_means: predicted means of the model for each configuration.
+    :param model_means: predicted std of the model for each configuration.
+    :param param_space: Space object for the optimization problem.
+    :param objective_weights: objective weights for multi-objective optimization. Not implemented yet.
+    :param threshold: threshold on objective values separating good points and bad points.
+    :param compute_bad: whether to compute the probability of being good or bad.
+    """
+    optimization_parameters = param_space.get_optimization_parameters()
+    probabilities = np.ones(len(model_means[optimization_parameters[0]]))
+
+    for parameter in optimization_parameters:
+        parameter_means = model_means[parameter]
+        parameter_stds = model_stds[parameter]
+        if compute_bad:
+            p = 1 - stats.norm.cdf((threshold[parameter] - parameter_means)/parameter_stds)
+        else:
+            p = stats.norm.cdf((threshold[parameter] - parameter_means)/parameter_stds)
+        probabilities *= p**objective_weights[parameter]
+
+    return probabilities
+
 def compute_gp_prediction_mean_and_std(bufferx, model, param_space):
     """
     Compute the mean and std prediction of a GP model for a list of points.
@@ -688,4 +717,4 @@ def sample_gp_posterior(bufferx, model, param_space):
     for objective in model:
         gp_samples[objective] = model[objective].posterior_samples_f(normalized_bufferx, size=1)
 
-    return gp_samples
+    return gp_samples