diff --git a/hypermapper/bo.py b/hypermapper/bo.py index 06b1339..40cb93d 100644 --- a/hypermapper/bo.py +++ b/hypermapper/bo.py @@ -120,6 +120,17 @@ def main(config, black_box_function=None, profiling=None): model_type = config["models"]["model"] optimization_method = config["optimization_method"] time_budget = config["time_budget"] + acquisition_function_optimizer = config["acquisition_function_optimizer"] + if ( + acquisition_function_optimizer == "cma_es" + and not param_space.is_space_continuous() + ): + print( + "Warning: CMA_ES can only be used with continuous search spaces (i.e. all parameters must be of type 'real')" + ) + print("Switching acquisition function optimizer to local search") + acquisition_function_optimizer = "local_search" + input_params = param_space.get_input_parameters() number_of_objectives = len(optimization_metrics) objective_limits = {} @@ -445,6 +456,7 @@ def main(config, black_box_function=None, profiling=None): objective_limits, classification_model, profiling, + acquisition_function_optimizer, ) else: @@ -459,7 +471,7 @@ def main(config, black_box_function=None, profiling=None): ) best_configuration = ( param_space.random_sample_configurations_without_repetitions( - tmp_fast_addressing_of_data_array, 1 + tmp_fast_addressing_of_data_array, 1, use_priors=False )[0] ) local_search_t1 = datetime.datetime.now() diff --git a/hypermapper/cma_es.py b/hypermapper/cma_es.py new file mode 100644 index 0000000..dbb3293 --- /dev/null +++ b/hypermapper/cma_es.py @@ -0,0 +1,170 @@ +import sys +import os +import numpy as np +import datetime +import cma +import contextlib +import copy +import random + +# ensure backward compatibility +try: + from hypermapper.local_search import get_min_configurations + from hypermapper import space + from hypermapper.utility_functions import * +except ImportError: + if os.getenv("HYPERMAPPER_HOME"): # noqa + warnings.warn( + "Found environment variable 'HYPERMAPPER_HOME', used to update the system path. Support might be discontinued in the future. Please make sure your installation is working without this environment variable, e.g., by installing with 'pip install hypermapper'.", + DeprecationWarning, + 2, + ) # noqa + sys.path.append(os.environ["HYPERMAPPER_HOME"]) # noqa + ppath = os.getenv("PYTHONPATH") + if ppath: + path_items = ppath.split(":") + + scripts_path = ["hypermapper/scripts", "hypermapper_dev/scripts"] + + if os.getenv("HYPERMAPPER_HOME"): + scripts_path.append(os.path.join(os.getenv("HYPERMAPPER_HOME"), "scripts")) + + truncated_items = [ + p for p in sys.path if len([q for q in scripts_path if q in p]) == 0 + ] + if len(truncated_items) < len(sys.path): + warnings.warn( + "Found hypermapper in PYTHONPATH. Usage is deprecated and might break things. " + "Please remove all hypermapper references from PYTHONPATH. Trying to import" + "without hypermapper in PYTHONPATH..." + ) + sys.path = truncated_items + + sys.path.append(".") # noqa + sys.path = list(OrderedDict.fromkeys(sys.path)) + + from hypermapper.local_search import get_min_configurations + from hypermapper.utility_functions import * + from hypermapper import space + + +def cma_es( + param_space, + data_array, + fast_addressing_of_data_array, + optimization_objective, + logfile, + optimization_function, + optimization_function_parameters, + cma_es_random_points=10000, + cma_es_starting_points=10, + sigma=0.001, +): + """ + Optimize the acquisition function using a mix of random and local search. + This algorithm random samples N points and then does a local search on the + best points from the random search and the best points from previous iterations (if any). + :param local_search_starting_points: an integer for the number of starting points for the local search. If 0, all points will be used. + :param local_search_random_points: number of random points to sample before the local search. + :param param_space: a space object containing the search space. + :param fast_addressing_of_data_array: A list containing the points that were already explored. + :param enable_feasible_predictor: whether to use constrained optimization. + :param optimization_function: the function that will be optimized by the local search. + :param optimization_function_parameters: a dictionary containing the parameters that will be passed to the optimization function. + :param optimization_objective: the name given to the scalarized values. + :param previous_points: previous points that have already been evaluated. + :return: all points evaluted and the best point found by the local search. + """ + t0 = datetime.datetime.now() + input_params = param_space.get_input_parameters() + input_param_objects = param_space.get_input_parameters_objects() + tmp_fast_addressing_of_data_array = copy.deepcopy(fast_addressing_of_data_array) + + # CMA-ES works better if optimizing between normalized ranges + param_min, param_max = {}, {} + normalizer = {} + unnormalizer = {} + for param in input_params: + param_min[param], param_max[param] = ( + input_param_objects[param].get_min(), + input_param_objects[param].get_max(), + ) + normalizer[param] = lambda x, input_param: (x - param_min[input_param]) / ( + param_max[input_param] - param_min[input_param] + ) + unnormalizer[param] = ( + lambda x, input_param: x * (param_max[input_param] - param_min[input_param]) + + param_min[input_param] + ) + + best_previous = get_min_configurations( + data_array, cma_es_starting_points, optimization_objective + ) + + concatenation_keys = input_params + cma_es_configurations = {} + cma_es_configurations = concatenate_data_dictionaries( + cma_es_configurations, best_previous, concatenation_keys + ) + # Get mode or random + param_mode = param_space.get_prior_mode() + cma_es_configurations = concatenate_data_dictionaries( + cma_es_configurations, param_mode, concatenation_keys + ) + + # Passing the dictionary with ** expands the key-value pairs into function parameters + # The acquisition functions return a tuple with two lists, cma wants only the first element of the first list + cmaes_black_box = lambda x: optimization_function( + configurations=[ + { + param: unnormalizer[param](x[param_idx], param) + for param_idx, param in enumerate(input_params) + } + ], + **optimization_function_parameters + )[0][0] + best_points = [] + best_point_vals = [] + + for configuration_idx in range( + len(cma_es_configurations[list(cma_es_configurations.keys())[0]]) + ): + x0 = [ + normalizer[param](cma_es_configurations[param][configuration_idx], param) + for param in input_params + ] + with open(logfile, "a") as f, contextlib.redirect_stdout(f): + es = cma.CMAEvolutionStrategy(x0, sigma, {"bounds": [0, 1]}) + es.optimize(cmaes_black_box) + + best_points.append(es.result.xbest) + best_point_vals.append(es.result.fbest) + + best_configuration_idx = np.argmin(best_point_vals) + best_configuration = { + param: unnormalizer[param]( + best_points[best_configuration_idx][param_idx], param + ) + for param_idx, param in enumerate(input_params) + } + configuration_string = param_space.get_unique_hash_string_from_values( + best_configuration + ) + # If the best configuration has already been evaluated before, remove it and get the next best configuration + if configuration_string in fast_addressing_of_data_array: + del best_point_vals[best_configuration_idx] + del best_points[best_configuration_idx] + best_configuration_idx = np.argmin(best_point_vals) + best_configuration = { + param: unnormalizer[param](best[param_idx], param) + for param_idx, param in enumerate(input_params) + } + configuration_string = param_space.get_unique_hash_string_from_values( + best_configuration + ) + + sys.stdout.write_to_logfile( + ("CMA-ES time %10.4f sec\n" % ((datetime.datetime.now() - t0).total_seconds())) + ) + + return best_configuration diff --git a/hypermapper/prior_optimization.py b/hypermapper/prior_optimization.py index 3751a7b..0e6349f 100644 --- a/hypermapper/prior_optimization.py +++ b/hypermapper/prior_optimization.py @@ -17,6 +17,7 @@ from hypermapper import models from hypermapper.local_search import local_search from hypermapper.utility_functions import dict_list_to_matrix + from hypermapper.cma_es import cma_es except ImportError: if os.getenv("HYPERMAPPER_HOME"): # noqa warnings.warn( @@ -51,6 +52,7 @@ from hypermapper import models from hypermapper.local_search import local_search from hypermapper.utility_functions import dict_list_to_matrix + from hypermapper.cma_es import cma_es def compute_probability_from_prior(configurations, param_space, objective_weights): @@ -178,6 +180,15 @@ def compute_EI_from_posteriors( :param posterior_normalization_limits: :param debug: whether to run in debug mode. """ + param_objects = param_space.get_input_parameters_objects() + for parameter in param_space.get_input_parameters(): + param_min, param_max = ( + param_objects[parameter].get_min(), + param_objects[parameter].get_max(), + ) + for configuration in configurations: + configuration[parameter] = min(configuration[parameter], param_max) + configuration[parameter] = max(configuration[parameter], param_min) user_prior_t0 = datetime.datetime.now() prior_good = compute_probability_from_prior( configurations, param_space, objective_weights @@ -336,6 +347,8 @@ def compute_EI_from_posteriors( good_bad_ratios = good_bad_ratios + feasibility_indicator good_bad_ratios = -1 * good_bad_ratios + good_bad_ratios[good_bad_ratios == float("inf")] = sys.maxsize + good_bad_ratios[good_bad_ratios == float("-inf")] = -1 * sys.maxsize good_bad_ratios = list(good_bad_ratios) sys.stdout.write_to_logfile( @@ -366,6 +379,7 @@ def prior_guided_optimization( objective_limits, classification_model=None, profiling=None, + acquisition_function_optimizer="local_search", ): """ Run a prior-guided bayesian optimization iteration. @@ -379,8 +393,6 @@ def prior_guided_optimization( :param objective_limits: estimated minimum and maximum limits for each objective. :param classification_model: feasibility classifier for constrained optimization. """ - local_search_starting_points = config["local_search_starting_points"] - local_search_random_points = config["local_search_random_points"] scalarization_key = config["scalarization_key"] number_of_cpus = config["number_of_cpus"] @@ -419,17 +431,46 @@ def prior_guided_optimization( float("-inf"), ] - _, best_configuration = local_search( - local_search_starting_points, - local_search_random_points, - param_space, - fast_addressing_of_data_array, - False, # set feasibility to false, we handle it inside the acquisition function - compute_EI_from_posteriors, - function_parameters, - scalarization_key, - number_of_cpus, - previous_points=data_array, - profiling=profiling, - ) + if acquisition_function_optimizer == "local_search": + local_search_starting_points = config["local_search_starting_points"] + local_search_random_points = config["local_search_random_points"] + _, best_configuration = local_search( + local_search_starting_points, + local_search_random_points, + param_space, + fast_addressing_of_data_array, + False, # set feasibility to false, we handle it inside the acquisition function + compute_EI_from_posteriors, + function_parameters, + scalarization_key, + number_of_cpus, + previous_points=data_array, + profiling=profiling, + ) + elif acquisition_function_optimizer == "cma_es": + logfile = deal_with_relative_and_absolute_path( + config["run_directory"], config["log_file"] + ) + sigma = config["cma_es_sigma"] + cma_es_starting_points = config["cma_es_starting_points"] + cma_es_random_points = config["cma_es_random_points"] + best_configuration = cma_es( + param_space, + data_array, + fast_addressing_of_data_array, + scalarization_key, + logfile, + compute_EI_from_posteriors, + function_parameters, + cma_es_random_points=cma_es_random_points, + cma_es_starting_points=cma_es_starting_points, + sigma=sigma, + ) + else: + print( + "Unrecognized acquisition function optimizer:", + acquisition_function_optimizer, + ) + raise SystemExit + return best_configuration diff --git a/hypermapper/random_scalarizations.py b/hypermapper/random_scalarizations.py index b449345..0d792cf 100644 --- a/hypermapper/random_scalarizations.py +++ b/hypermapper/random_scalarizations.py @@ -27,6 +27,7 @@ reciprocate_weights, compute_data_array_scalarization, ) + from hypermapper.cma_es import cma_es except ImportError: if os.getenv("HYPERMAPPER_HOME"): # noqa warnings.warn( @@ -66,6 +67,7 @@ reciprocate_weights, compute_data_array_scalarization, ) + from hypermapper.cma_es import cma_es def run_acquisition_function( @@ -514,6 +516,7 @@ def random_scalarizations( objective_limits, classification_model=None, profiling=None, + acquisition_function_optimizer="local_search", ): """ Run one iteration of bayesian optimization with random scalarizations. @@ -551,17 +554,44 @@ def random_scalarizations( ] optimization_function_parameters["number_of_cpus"] = config["number_of_cpus"] - _, best_configuration = local_search( - local_search_starting_points, - local_search_random_points, - param_space, - fast_addressing_of_data_array, - False, # we do not want the local search to consider feasibility constraints, only the acquisition functions - run_acquisition_function, - optimization_function_parameters, - scalarization_key, - number_of_cpus, - previous_points=data_array, - profiling=profiling, - ) + if acquisition_function_optimizer == "local_search": + _, best_configuration = local_search( + local_search_starting_points, + local_search_random_points, + param_space, + fast_addressing_of_data_array, + False, # we do not want the local search to consider feasibility constraints, only the acquisition functions + run_acquisition_function, + optimization_function_parameters, + scalarization_key, + number_of_cpus, + previous_points=data_array, + profiling=profiling, + ) + elif acquisition_function_optimizer == "cma_es": + logfile = deal_with_relative_and_absolute_path( + config["run_directory"], config["log_file"] + ) + sigma = config["cma_es_sigma"] + cma_es_starting_points = config["cma_es_starting_points"] + cma_es_random_points = config["cma_es_random_points"] + best_configuration = cma_es( + param_space, + data_array, + fast_addressing_of_data_array, + scalarization_key, + logfile, + compute_EI_from_posteriors, + function_parameters, + cma_es_random_points=cma_es_random_points, + cma_es_starting_points=cma_es_starting_points, + sigma=sigma, + ) + else: + print( + "Unrecognized acquisition function optimizer:", + acquisition_function_optimizer, + ) + raise SystemExit + return best_configuration diff --git a/hypermapper/schema.json b/hypermapper/schema.json index b99f5de..45ebdd8 100644 --- a/hypermapper/schema.json +++ b/hypermapper/schema.json @@ -178,6 +178,16 @@ ], "default" : "uniform" }, + "custom_gaussian_prior_means":{ + "type": ["array", "boolean"], + "description": "means for the custom gaussian prior for this parameter. Can only be used with real parameters for now.", + "default": false + }, + "custom_gaussian_prior_stds":{ + "type": ["array", "boolean"], + "description": "standard deviations for the custom gaussian prior for this parameter. Can only be used with real parameters for now.", + "default": false + }, "parameter_type": { "description": "The type of the parameter that is being defined.", "type": "string", @@ -354,8 +364,8 @@ }, "acquisition_function_optimizer":{ "type": "string", - "description": "which method to use to optimize the acquisition function.", - "enum":["local_search"], + "description": "which method to use to optimize the acquisition function. CMA_ES only works on continuous spaces (i.e. only parameters of type 'real').", + "enum":["local_search", "cma_es"], "default": "local_search" }, "evolution_population_size":{ diff --git a/hypermapper/space.py b/hypermapper/space.py index 001cdea..153b884 100644 --- a/hypermapper/space.py +++ b/hypermapper/space.py @@ -177,6 +177,24 @@ def randomly_select(self, size=1): # remove excess samples to get the proper amount samples = samples[0:size] + elif self.prior == "gaussian_mixture": + sample_idx = random.randint(0, len(self.custom_gaussian_prior_mean) - 1) + sample = norm.rvs( + loc=self.custom_gaussian_prior_mean[sample_idx], + scale=self.custom_gaussian_prior_std[sample_idx], + ) + sample_count = 1 + while sample > self.max_value or sample < self.min_value: + sample_count += 1 + if sample_count > 10000: + print( + "\n ####\n Warning: reached maximum number of invalid samples from the custom gaussian prior. \nThe prior sampling will stop now. Is the prior defined outside or too close to the interval edges?\n" + ) + raise SystemExit # too many random samples failed, the prior is probably too close to the edge + sample = norm.rvs( + loc=self.custom_gaussian_prior_mean[sample_idx], + scale=self.custom_gaussian_prior_std[sample_idx], + ) else: samples = beta(self.alpha, self.beta, size) samples = self.from_range_0_1_to_parameter_value(samples) @@ -260,6 +278,16 @@ def get_x_probability(self, x): loc=self.custom_gaussian_prior_mean, scale=self.custom_gaussian_prior_std, ) + elif self.prior == "gaussian_mixture": + prob = 0 + mixture_weights = 1 / len(self.custom_gaussian_prior_mean) + for idx in range(len(self.custom_gaussian_prior_mean)): + prob += mixture_weights * norm.pdf( + x, + loc=self.custom_gaussian_prior_mean[idx], + scale=self.custom_gaussian_prior_std[idx], + ) + return prob else: rescaled_x = self.from_parameter_value_to_0_1_range(x) pdf = stats.beta.pdf(rescaled_x, self.alpha, self.beta) @@ -290,6 +318,10 @@ def get_max(self): def set_custom_gaussian_prior_params(self, mean, std): if std == -1: std = (self.max_value - self.min_value) / 2 + if type(std) is list: + for idx in range(len(std)): + if std[idx] == -1: + std[idx] = (self.max_value - self.min_value) / 2 self.custom_gaussian_prior_mean = mean self.custom_gaussian_prior_std = std @@ -710,6 +742,7 @@ def __init__(self, config): self.all_input_parameters = OrderedDict() self.input_non_categorical_parameters = OrderedDict() self.input_categorical_parameters = OrderedDict() + self.is_continuous = True self.parse_input_parameters(config["input_parameters"]) if self.get_discrete_space_size() > max_number_of_predictions: @@ -777,7 +810,7 @@ def __init__(self, config): prior_estimation_data[self.optimization_metrics[0]], point_quantile ) indices = np.nonzero( - prior_estimation_data[self.optimization_metrics[0]] < threshold + prior_estimation_data[self.optimization_metrics[0]] <= threshold )[0] good_points = all_input_data[:, indices] self.pdf = gaussian_kde(good_points, bw_method=self.bw_selector) @@ -798,7 +831,7 @@ def __init__(self, config): point_quantile, ) indices_good = np.nonzero( - prior_estimation_data[self.optimization_metrics[0]] < threshold + prior_estimation_data[self.optimization_metrics[0]] <= threshold )[0] for input_param in self.get_input_parameters(): if self.parameters_type[input_param] == "real": @@ -810,34 +843,66 @@ def __init__(self, config): ) elif config["input_parameters"][param]["prior"] == "custom_gaussian": self.normalize_priors = True - if len(config["custom_gaussian_prior_means"]) == 1: - mean = config["custom_gaussian_prior_means"][0] - elif len(config["custom_gaussian_prior_means"]) == len( - config["input_parameters"] + if ( + config["input_parameters"][param]["custom_gaussian_prior_means"] + is False + ) or ( + config["input_parameters"][param]["custom_gaussian_prior_stds"] + is False ): - mean = config["custom_gaussian_prior_means"][param_idx] - else: - print( - "Error: the custom_gaussian prior means array must be either 1 or", - len(config["input_parameters"]), - "received", - len(config["custom_gaussian_prior_means"]), + if len(config["custom_gaussian_prior_means"]) == 1: + mean = config["input_parameters"][param][ + "custom_gaussian_prior_means" + ][0] + else: + mean = config["input_parameters"][param][ + "custom_gaussian_prior_means" + ] + if len(config["custom_gaussian_prior_stds"]) == 1: + std = config["input_parameters"][param][ + "custom_gaussian_prior_stds" + ][0] + else: + std = config["input_parameters"][param][ + "custom_gaussian_prior_stds" + ] + if len(std) != len(mean): + print( + "Invalid custom gaussian parameters. The number of means and standard deviations must be equal, but got:" + ) + print(f"{len(mean)} means and {len(std)} standard deviations") + self.all_input_parameters[param].set_custom_gaussian_prior_params( + mean, std ) - raise SystemExit - if len(config["custom_gaussian_prior_stds"]) == 1: - std = config["custom_gaussian_prior_stds"][0] - elif len(config["custom_gaussian_prior_stds"]) == len( - config["input_parameters"] - ): - std = config["custom_gaussian_prior_stds"][param_idx] else: - print( - "Error: the custom_gaussian prior stds array must be either 1 or", - len(config["input_parameters"]), - "received", - len(config["custom_gaussian_prior_stds"]), - ) - raise SystemExit + if len(config["custom_gaussian_prior_means"]) == 1: + mean = config["custom_gaussian_prior_means"][0] + elif len(config["custom_gaussian_prior_means"]) == len( + config["input_parameters"] + ): + mean = config["custom_gaussian_prior_means"][param_idx] + else: + print( + "Error: the custom_gaussian prior means array must be either 1 or", + len(config["input_parameters"]), + "received", + len(config["custom_gaussian_prior_means"]), + ) + raise SystemExit + if len(config["custom_gaussian_prior_stds"]) == 1: + std = config["custom_gaussian_prior_stds"][0] + elif len(config["custom_gaussian_prior_stds"]) == len( + config["input_parameters"] + ): + std = config["custom_gaussian_prior_stds"][param_idx] + else: + print( + "Error: the custom_gaussian prior stds array must be either 1 or", + len(config["input_parameters"]), + "received", + len(config["custom_gaussian_prior_stds"]), + ) + raise SystemExit self.all_input_parameters[param].set_custom_gaussian_prior_params( mean, std ) @@ -863,6 +928,7 @@ def parse_input_parameters(self, input_parameters_json): else: param_distribution = param["prior"] if param_type == "ordinal": + self.is_continuous = False param_values = param["values"] self.all_input_parameters[param_name] = OrdinalParameter( values=param_values, @@ -882,6 +948,7 @@ def parse_input_parameters(self, input_parameters_json): else: self.parameters_python_type[param_name] = "float" elif param_type == "categorical": + self.is_continuous = False param_values = param["values"] if isinstance(param_distribution, str): if param_distribution != "uniform": @@ -932,6 +999,7 @@ def parse_input_parameters(self, input_parameters_json): self.parameters_type[param_name] = "real" self.parameters_python_type[param_name] = "float" elif param_type == "integer": + self.is_continuous = False param_min, param_max = param["values"] if param_min > param_max: param_min, param_max = ( @@ -1108,6 +1176,12 @@ def get_estimate_prior_flags(self): """ return [self.estimate_prior_flag] + def is_space_continuous(self): + """ + :return: true if all parameters are real, otherwise returns false + """ + return self.is_continuous + def get_space_size(self): """ Get the the size of the Cartesian product of the parameters of input previously declared in the json file. diff --git a/setup.py b/setup.py index fc1c59d..1052d82 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setup( name="hypermapper", - version="2.2.5", + version="2.2.6", description="HyperMapper is a multi-objective black-box optimization tool based on Bayesian Optimization.", long_description=long_description, long_description_content_type="text/markdown", @@ -19,6 +19,7 @@ "Operating System :: OS Independent", ], install_requires=[ + "cma", "numpy", "matplotlib", "ply",