diff --git a/CHANGELOG.md b/CHANGELOG.md index d2c4ef4f..3997cde5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,11 +4,17 @@ ### Changes -- +- Added the get_Brownie_Bee_Pareto() function to support Pareto front plotting in the Brownie Bee user interface. +- Added the get_Pareto_front_compromise() function to facilitate a default highlighted compromise in the Pareto plot. +- Added the get_Brownie_Bee_1d_plot() function to support custom 1D dependency plots in the Brownie Bee user interface. ### Bugfixes -- +- opt.ask() is now reproducible for multiobjective optimization through seeding of the NSGAII algorithm and + a small adjustment to the way optimizer._tell() chooses between Steinerberger and NSGAII sampling. +- Pareto front calculations are now reproducible. +- Fixes small error in legend generation of plot_objective_1d for the case where the user specifies a + specific set of settings to create the plot at ## Version 1.1.1 [published] diff --git a/ProcessOptimizer/optimizer/optimizer.py b/ProcessOptimizer/optimizer/optimizer.py index e44cce8e..9dd8e5ce 100644 --- a/ProcessOptimizer/optimizer/optimizer.py +++ b/ProcessOptimizer/optimizer/optimizer.py @@ -753,7 +753,7 @@ def _tell(self, x, y, fit=True): prob_stbr = 0.25 # Simulate a random number - random_uniform_number = np.random.uniform() + random_uniform_number = self.rng.uniform() # The random number decides what strategy to use for the next point if random_uniform_number < prob_stbr: @@ -1309,12 +1309,18 @@ def __MinimalDistance(X, Y): # This function calls NSGAII to estimate the Pareto Front def NSGAII(self, MU=40): from ._NSGA2 import NSGAII - + + if MU % 4 != 0: + raise ValueError( + "Number of simulated points must be divisible by 4 for the NSGAII algorithm" + ) + pop, logbook, front = NSGAII( self.n_objectives, self.__ObjectiveGP, np.array(self.space.transformed_bounds), MU=MU, + seed=42, ) return pop, logbook, front diff --git a/ProcessOptimizer/plots.py b/ProcessOptimizer/plots.py index dd60acb5..20a8d053 100644 --- a/ProcessOptimizer/plots.py +++ b/ProcessOptimizer/plots.py @@ -11,6 +11,7 @@ from scipy.stats.mstats import mquantiles from scipy.stats import norm from scipy.ndimage import gaussian_filter1d +from scipy.signal import savgol_filter from warnings import warn from ProcessOptimizer import expected_minimum, expected_minimum_random_sampling from .space import Categorical, Integer @@ -1453,9 +1454,9 @@ def plot_objective_1d( highlight_label = "Expected minimum" elif pars == "expected_minimum_random": highlight_label = "Simulated minimum" - elif isinstance(pars, list): - # The case where the user specifies [x[0], x[1], ...] - highlight_label = "Point: " + str(pars) + elif isinstance(pars, list): + # The case where the user specifies [x[0], x[1], ...] + highlight_label = "Point: " + str(pars) # Legend icon for the highlighted value legend_hl = mpl.lines.Line2D( [], @@ -1514,6 +1515,102 @@ def plot_objective_1d( ax, space, ylabel=ylabel, dim_labels=dimensions ) + +def get_Brownie_Bee_1d_plot( + result, + x_eval=None, + n_points=60, + n_samples=250, +): + """Returns the information needed to produce single factor dependence plots + of the model in `result`. This utility function is mainly intended for use + with the Brownie Bee user interface. + + The data returned allows visualization of how Y depends on dimension `i` + when all other factor values are locked to those provided in x_eval. + + Parameters + ---------- + * `result` [`OptimizeResult`] + The result for which to create the plotting data. + + * `x_eval` [list or np.ndarray of floats, ints and/or strings, default=None] + [x[0], x[1], ..., x[n]] - Factor settings to create the dependency plot + at, defined by the values in this list. Depending on the system, this + list can contain a mixture of floats, ints and strings. If no settings + are provided the function defaults to using expected_minimum on the + model. + + * `n_points` [int, default=60] + Number of points at which to evaluate the partial dependence + along each dimension. + + * `n_samples` [int, default=250] + Number of random samples to use for averaging the model function + at each of the `n_points`. + + + Returns + ------- + * `plot_list`: [`list of lists`]: + A list of lists that provide the necessary data for creating dependency + plots along each dimension of the space. Each list contains three lists + and a float of the x-axis value to highlight in the same plot: + [[x-axis], [y_low], [y_high], x_highlight]. + + The last entry in plot_list contains information on the mean and std of + the model predictions at x_eval. This data allows the user to build a + histogram of expected results at these settings. + """ + space = result.space + model = result.models[-1] + + if x_eval is None: + # Identify the location of the expected minimum, and its mean and std + x_eval, [res_mean, res_std] = expected_minimum( + result, + n_random_starts=20, + random_state=42, + return_std=True, + ) + elif isinstance(x_eval, list) or isinstance(x_eval, np.ndarray): + assert len(x_eval) == len(space), "Input settings must have same length as number of features" + res_mean, res_std = model.predict(np.array(x_eval).reshape(1, -1), return_std=True) + else: + raise TypeError("x_eval must be a list of settings, or None") + + rvs_transformed = space.transform(space.rvs(n_samples=n_samples)) + # Map the settings to highlight so we automatically take care of + # categorical factors that use 1-hot encoding + _, highlight, _ = _map_categories(space, result.x_iters, x_eval) + + # Gather all data relevant for plotting + plot_list = [] + + # Generate 1D plot information for each dimension in space + for i in range(space.n_dims): + xi, yi, stddevs = dependence( + space, + model, + i, + j=None, + sample_points=rvs_transformed, + n_points=n_points, + x_eval=x_eval, + ) + plot_list.append([ + xi.tolist(), + (yi-1.96*stddevs).tolist(), + (yi+1.96*stddevs).tolist(), + highlight[i], + ]) + + # Add information about the expected results at x_eval + plot_list.append([res_mean, res_std]) + + return plot_list + + def plot_brownie_bee_frontend( result, n_points=60, @@ -2571,3 +2668,76 @@ def plot_Pareto_bokeh( json_item = bh_embed.json_item(p) return json_item + +def get_Brownie_Bee_Pareto(optimizer, n_points=100): + """Calculate Pareto front in two dimensions and return its points, as well + as uncertainty band around each objective along the front. This function is + mainly intended for use with the Brownie Bee user interface. + + Parameters + ---------- + * `optimizer` [`Optimizer`] + The optimizer containing data and the multiobjective model + * `n_points` [int, default=100] + The number of points to simulate for the Pareto front. Must be a + multiple of 4 for the NSGAII algorithm. + + Returns + ------- + * `front_x`: [numpy.ndarray]: + Pareto front locations in optimizer X-space + * `front_y`: [numpy.ndarray]: + Pareto front locations in optimizer Y-space + * `objective1_error`: [numpy.ndarray]: + Uncertainty (1.96*std) of objective 1 values at the front_y locations, + including observational noise + * `objective2_error`: [numpy.ndarray]: + Uncertainty (1.96*std) of objective 2 values at the front_y locations, + including observational noise + """ + + if optimizer.models == []: + raise ValueError("No models have been fitted yet") + + if optimizer.n_objectives == 1: + raise ValueError( + "get_Pareto_points is not possible with single objective optimization" + ) + + if optimizer.n_objectives > 2: + raise ValueError("get_Pareto_points is not possible with >2 objectives") + + if n_points % 4 != 0: + raise ValueError( + "Number of simulated points must be divisible by 4 for the NSGAII algorithm" + ) + + # Estimate the Pareto front of the models in the optimizer object + front_x, logbook, front_y = optimizer.NSGAII(MU=n_points) + + front_x = np.asarray(front_x) + front_x = np.asarray( + optimizer.space.inverse_transform( + front_x.reshape(len(front_x), optimizer.space.transformed_n_dims) + ) + ) + + # Sort the points in ascending order on objective 1 + idx = np.argsort(front_y[:, 0]) + front_x = front_x[idx, :] + front_y = front_y[idx, :] + + # Extract estimates of the objective functions along the Pareto front + output = optimizer.estimate(front_x) + # Grab objective errors and smooth them slightly along the front + objective1_error = [1.96*point.Y1.std for point in output] + objective1_error = savgol_filter(objective1_error, 11, 1, mode="interp") + objective2_error = [1.96*point.Y2.std for point in output] + objective2_error = savgol_filter(objective2_error, 11, 1, mode="interp") + + return ( + front_x, + front_y, + objective1_error, + objective2_error, + ) diff --git a/ProcessOptimizer/tests/test_multiobjective.py b/ProcessOptimizer/tests/test_multiobjective.py index 8f863782..366ea848 100644 --- a/ProcessOptimizer/tests/test_multiobjective.py +++ b/ProcessOptimizer/tests/test_multiobjective.py @@ -6,6 +6,7 @@ from ProcessOptimizer import Optimizer +from ProcessOptimizer.model_systems import get_model_system @pytest.mark.fast_test @@ -46,3 +47,47 @@ def test_Pareto_in_space(): # Assert that Pareto points are in space for x in pop: assert_equal(opt.space.__contains__(x), True) + + +@pytest.mark.fast_test +def test_Pareto_reproducible(): + gold_model_system = get_model_system('gold_map') + distance_model_system = get_model_system('distance_map', camp_coordinates=(4,10)) + + opt = Optimizer(gold_model_system.space, n_initial_points=4, n_objectives=2) + + gold_model_system.noise_model.set_seed(40) + distance_model_system.noise_model.set_seed(40) + + for i in range(10): + new_dig_site = opt.ask() + gold_found = gold_model_system.get_score(new_dig_site) + distance = distance_model_system.get_score(new_dig_site) + opt.tell(new_dig_site, [gold_found, distance]) + + # Calculate Pareto front twice and compare + pop1, logbook, front1 = opt.NSGAII() + pop2, logbook, front2 = opt.NSGAII() + assert_equal(pop1, pop2) + assert_equal(front1, front2) + + +@pytest.mark.fast_test +def test_multiobjective_ask_reproducible(): + gold_model_system = get_model_system('gold_map') + distance_model_system = get_model_system('distance_map', camp_coordinates=(4,10)) + + opt1 = Optimizer(gold_model_system.space, n_initial_points=4, n_objectives=2) + opt2 = Optimizer(gold_model_system.space, n_initial_points=4, n_objectives=2) + + gold_model_system.noise_model.set_seed(40) + distance_model_system.noise_model.set_seed(40) + # Check that the two optimizers stay in sync + for i in range(40): + new_dig_site_1 = opt1.ask() + new_dig_site_2 = opt2.ask() + assert_equal(new_dig_site_1, new_dig_site_2) + gold_found = gold_model_system.get_score(new_dig_site_1) + distance = distance_model_system.get_score(new_dig_site_1) + opt1.tell(new_dig_site_1, [gold_found, distance]) + opt2.tell(new_dig_site_2, [gold_found, distance]) diff --git a/ProcessOptimizer/utils/utils.py b/ProcessOptimizer/utils/utils.py index 4eb60675..c64bd30c 100644 --- a/ProcessOptimizer/utils/utils.py +++ b/ProcessOptimizer/utils/utils.py @@ -583,3 +583,43 @@ def y_coverage(res, return_plot=False, random_state=None, horizontal=False): plt.show() return (observed_min, observed_max), (expected_min, expected_max) + + +def get_Pareto_front_compromise(front_y): + """Support function to identify a default compromise between two objectives + in a Pareto front. This function is mainly intended for use with the + Brownie Bee user interface. + + Parameters + ---------- + * `front_y` [numpy.ndarray] + Pareto front locations in optimizer Y-space + + + Returns + ------- + * `best_idx`: [int]: + The index of the point in front_y that is closest to the normalized + ideal solution for both objectives. + """ + if not isinstance(front_y, np.ndarray): + raise TypeError("front_y must be a numpy.ndarray, got (%s)." % type(front_y)) + + if front_y.shape[1] != 2: + raise ValueError("front_y must have two columns, got (%s)" % front_y.shape[1]) + + # Get the best predicted point for each objective + best_obj1 = np.min(front_y[:, 0]) + best_obj2 = np.min(front_y[:, 1]) + # Get the span of each objective for normalization + span1 = np.ptp(front_y[:, 0]) + span2 = np.ptp(front_y[:, 1]) + + # Calculate distance from ideal solution to points along the front + dist = np.sqrt( + ((front_y[:, 0] - best_obj1)/span1)**2 + + ((front_y[:, 1] - best_obj2)/span2)**2 + ) + best_idx = np.argmin(dist) + + return best_idx diff --git a/pyproject.toml b/pyproject.toml index ea3f676c..2e94aa42 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ browniebee = [ "scipy==1.13.0", "six==1.16.0", "patsy==1.0.1", + "torch==2.8.0", ] release = [ "build==1.2.2.post1"