From a7c02320462be5a99ff0a591735d1654f26d693a Mon Sep 17 00:00:00 2001 From: Tim Mensinger Date: Sun, 6 Nov 2022 14:04:58 +0100 Subject: [PATCH] Improve trust region sampling in tranquilo (#388) --- .../optimization/tranquilo/options.py | 7 + .../optimization/tranquilo/sample_points.py | 331 ++++++++++++++++-- .../optimization/tranquilo/tranquilo.py | 24 +- .../tranquilo/tranquilo_history.py | 15 +- .../tranquilo/test_filter_points.py | 1 + .../tranquilo/test_sample_points.py | 142 +++++++- 6 files changed, 480 insertions(+), 40 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/options.py b/src/estimagic/optimization/tranquilo/options.py index a45143c4b..17009bbf6 100644 --- a/src/estimagic/optimization/tranquilo/options.py +++ b/src/estimagic/optimization/tranquilo/options.py @@ -56,3 +56,10 @@ class RadiusOptions(NamedTuple): class TrustRegion(NamedTuple): center: np.ndarray radius: float + + +class RadiusFactors(NamedTuple): + accepatance: float = 0.02 + centric: float = 0.1 + outer: float = 0.6 + neighborhood: float = 1.5 diff --git a/src/estimagic/optimization/tranquilo/sample_points.py b/src/estimagic/optimization/tranquilo/sample_points.py index a945b6dac..2fdc4f556 100644 --- a/src/estimagic/optimization/tranquilo/sample_points.py +++ b/src/estimagic/optimization/tranquilo/sample_points.py @@ -2,34 +2,46 @@ import warnings from functools import partial +import estimagic as em import numpy as np from estimagic.optimization.tranquilo.options import Bounds +from scipy.spatial.distance import pdist +from scipy.special import logsumexp -def get_sampler(sampler, bounds, user_options=None): +def get_sampler( + sampler, bounds, model_info=None, radius_factors=None, user_options=None +): """Get sampling function partialled options. Args: sampler (str or callable): Name of a sampling method or sampling function. - The arguments of sampling functions need to be: ``lower_bounds``, - ``upper_bounds``, ``target_size``, ``existing_xs``, ``existing_fvals``. + The arguments of sampling functions need to be: ``trustregion``, + ``target_size``, ``rng``, ``existing_xs`` and ``bounds``. Sampling functions need to return a dictionary with the entry "points" (and arbitrary additional information). See ``reference_sampler`` for details. bounds (Bounds): A NamedTuple with attributes ``lower`` and ``upper`` user_options (dict): Additional keyword arguments for the sampler. Options that - are not used by the sampler are ignored with a warning. + are not used by the sampler are ignored with a warning. If sampler is + 'hull_sampler' or 'optimal_hull_sampler' the user options must contain the + argument 'order', which is a positive integer. Returns: callable: Function that depends on trustregion, target_size, existing_xs and - existing_fvals and returns a new sample. + existing_fvals, model_info and and returns a new sample. """ user_options = {} if user_options is None else user_options built_in_samplers = { "naive": _naive_sampler, - "sphere": _sphere_sampler, + "hull_sampler": _hull_sampler, + "optimal_hull_sampler": _optimal_hull_sampler, + "cube": partial(_hull_sampler, order=np.inf), + "sphere": partial(_hull_sampler, order=2), + "optimal_cube": partial(_optimal_hull_sampler, order=np.inf), + "optimal_sphere": partial(_optimal_hull_sampler, order=2), } if isinstance(sampler, str) and sampler in built_in_samplers: @@ -40,10 +52,19 @@ def get_sampler(sampler, bounds, user_options=None): _sampler_name = getattr(sampler, "__name__", "your sampler") else: raise ValueError( - "Invalid sampler: {sampler}. Must be one of {list(built_in_samplers)} " + f"Invalid sampler: {sampler}. Must be one of {list(built_in_samplers)} " "or a callable." ) + if "hull_sampler" in _sampler_name and "order" not in user_options: + msg = ( + "The hull_sampler and optimal_hull_sampler require the argument 'order' to " + "be prespecfied in the user_options dictionary. Order is a positive " + "integer. For order = 2 the hull_sampler equals the sphere_sampler, and " + "for order = np.inf it equals the cube_sampler." + ) + raise ValueError(msg) + args = set(inspect.signature(_sampler).parameters) mandatory_args = { @@ -54,6 +75,13 @@ def get_sampler(sampler, bounds, user_options=None): "rng", } + optional_kwargs = { + "model_info": model_info, + "radius_factors": radius_factors, + } + + optional_kwargs = {k: v for k, v in optional_kwargs.items() if k in args} + problematic = mandatory_args - args if problematic: raise ValueError( @@ -77,6 +105,7 @@ def get_sampler(sampler, bounds, user_options=None): out = partial( _sampler, bounds=bounds, + **optional_kwargs, **reduced, ) @@ -93,7 +122,8 @@ def _naive_sampler( """Naive random generation of trustregion points. This is just a reference implementation to illustrate the interface of trustregion - samplers. + samplers. Mathematically it samples uniformaly from inside the cube defined by the + intersection of the trustregion and the bounds. All arguments but seed are mandatory, even if not used. @@ -110,47 +140,300 @@ def _naive_sampler( x vector at which the criterion function has already been evaluated, that satisfies lower_bounds <= existing_xs <= upper_bounds. rng (numpy.random.Generator): Random number generator. - bounds (Bounds or None): NamedTuple + bounds (Bounds or None): NamedTuple. """ - n_points = _get_effective_n_points(target_size, existing_xs) + n_points = _get_effective_n_points(target_size, existing_xs=existing_xs) n_params = len(trustregion.center) - region_bounds = _get_effective_bounds(trustregion, bounds) + effective_bounds = _get_effective_bounds(trustregion, bounds=bounds) points = rng.uniform( - low=region_bounds.lower, - high=region_bounds.upper, + low=effective_bounds.lower, + high=effective_bounds.upper, size=(n_points, n_params), ) + return points + + +def _hull_sampler( + trustregion, + target_size, + rng, + order, + distribution=None, + existing_xs=None, + bounds=None, +): + """Random generation of trustregion points on the hull of general sphere / cube. + + Points are sampled randomly on a hull (of a sphere for order=2 and of a cube for + order=np.inf). These points are then mapped into the feasible region, which is + defined by the intersection of the trustregion and the bounds. + + Args: + trustregion (TrustRegion): NamedTuple with attributes center and radius. + target_size (int): Target number of points in the combined sample of existing_xs + and newly sampled points. The sampler does not have to guarantee that this + number will actually be reached. + rng (numpy.random.Generator): Random number generator. + order (int): Type of norm to use when scaling the sampled points. For 2 it will + result in sphere sampling, for np.inf in cube sampling. + distribution (str): Distribution to use for initial sample before points are + projected onto unit hull. Must be in {'normal', 'uniform'}. + existing_xs (np.ndarray or None): 2d numpy array in which each row is an + x vector at which the criterion function has already been evaluated, that + satisfies lower_bounds <= existing_xs <= upper_bounds. + bounds (Bounds or None): NamedTuple. + """ + n_points = _get_effective_n_points(target_size, existing_xs=existing_xs) + n_params = len(trustregion.center) + effective_bounds = _get_effective_bounds(trustregion, bounds=bounds) + + if distribution is None: + distribution = "normal" if order <= 3 else "uniform" + points = _draw_from_distribution(distribution, rng=rng, size=(n_points, n_params)) + points = _project_onto_unit_hull(points, order=order) + points = _map_into_feasible_trustregion(points, bounds=effective_bounds) return points -def _sphere_sampler( +def _optimal_hull_sampler( trustregion, target_size, rng, + model_info, + radius_factors, + order, + distribution=None, + hardness=1, existing_xs=None, bounds=None, + algorithm="scipy_lbfgsb", + algo_options=None, ): - n_points = _get_effective_n_points(target_size, existing_xs) + """Optimal generation of trustregion points on the hull of general sphere / cube. + + Points are sampled optimally on a hull (of a sphere for order=2 and of a cube for + order=np.inf), where the criterion that is maximized is the minimum distance of all + pairs of points, except for pairs of existing points. These points are then mapped + into the feasible region, which is defined by the intersection of the trustregion + and the bounds. Instead of using a hard minimum we return the soft minimum, whose + accuracy we govern by the hardness factor. For more information on the soft-minimum, + seek: https://tinyurl.com/mrythbk4. + + Args: + trustregion (TrustRegion): NamedTuple with attributes center and radius. + target_size (int): Target number of points in the combined sample of existing_xs + and newly sampled points. The sampler does not have to guarantee that this + number will actually be reached. + rng (numpy.random.Generator): Random number generator. + order (int): Type of norm to use when scaling the sampled points. For 2 it will + result in sphere sampling, for np.inf in cube sampling. + distribution (str): Distribution to use for initial sample before points are + projected onto unit hull. Must be in {'normal', 'uniform'}. + hardness (float): Positive scaling factor. As hardness tends to infinity the + soft minimum (logsumexp) approaches the hard minimum. Default is 1. A + detailed explanation is given in the docstring. + existing_xs (np.ndarray or None): 2d numpy array in which each row is an + x vector at which the criterion function has already been evaluated, that + satisfies lower_bounds <= existing_xs <= upper_bounds. + bounds (Bounds or None): NamedTuple. + algorithm (str): Optimization algorithm. + algo_options (dict): Algorithm specific configuration of the optimization. See + :ref:`list_of_algorithms` for supported options of each algorithm. Default + sets ``stopping_max_iterations=n_params``. + + Returns: + np.ndarray: Generated points. Has shape (target_size, len(trustregion.center)). + + """ + n_points = _get_effective_n_points(target_size, existing_xs=existing_xs) n_params = len(trustregion.center) - raw = rng.normal(size=(n_points, n_params)) - denom = np.linalg.norm(raw, axis=1).reshape(-1, 1) - points = trustregion.radius * raw / denom + trustregion.center + if n_points <= 0: + return np.array([]) + + algo_options = {} if algo_options is None else algo_options + if "stopping_max_iterations" not in algo_options: + algo_options["stopping_max_iterations"] = 2 * n_params + + effective_bounds = _get_effective_bounds(trustregion, bounds=bounds) + + if existing_xs is not None: + # map existing points into unit space for easier optimization + existing_xs_unit = _map_from_feasible_trustregion(existing_xs, effective_bounds) + + dist_to_center = np.linalg.norm(existing_xs_unit, axis=1) + not_centric = dist_to_center >= radius_factors.centric - if bounds is not None and (bounds.lower is not None or bounds.upper is not None): - bounds = _get_effective_bounds(trustregion, bounds) - points = np.clip( - points, - a_min=bounds.lower, - a_max=bounds.upper, + if not_centric.any(): + existing_xs_unit = existing_xs_unit[not_centric] + else: + existing_xs_unit = None + + else: + existing_xs_unit = None + + # start params + if distribution is None: + distribution = "normal" if order <= 3 else "uniform" + x0 = _draw_from_distribution(distribution, rng=rng, size=(n_points, n_params)) + x0 = _project_onto_unit_hull(x0, order=order) + x0 = x0.flatten() # flatten so that em.maximize uses fast path + + # This would raise an error because there are zero pairs to calculate the + # pairwise distance + if existing_xs_unit is None and n_points == 1: + opt_params = x0 + else: + res = em.maximize( + criterion=_minimal_pairwise_distance_on_hull, + params=x0, + algorithm=algorithm, + criterion_kwargs={ + "existing_xs": existing_xs_unit, + "order": order, + "hardness": hardness, + "n_params": n_params, + }, + lower_bounds=-np.ones_like(x0), + upper_bounds=np.ones_like(x0), + algo_options=algo_options, ) + opt_params = res.params + + points = _project_onto_unit_hull(opt_params.reshape(-1, n_params), order=order) + points = _map_into_feasible_trustregion(points, bounds=effective_bounds) return points +# ====================================================================================== +# Helper functions +# ====================================================================================== + + +def _minimal_pairwise_distance_on_hull(x, existing_xs, order, hardness, n_params): + """Compute minimal pairwise distance of new and existing points. + + Instead of optimizing the distance of points in the feasible trustregion, this + criterion function leads to the maximization of the minimum distance of the points + in the unit space. These can then be mapped into the feasible trustregion. We do not + consider the distances between existing points. Instead of using a hard minimum we + return the soft minimum, whose accuracy we govern by the hardness factor. For more + information on the soft-minimum, seek: https://tinyurl.com/mrythbk4. + + Args: + x (np.ndarray): Flattened 1d array of internal points. Each value is in [-1, 1]. + existing_xs (np.ndarray or None): 2d numpy array in which each row is an + x vector at which the criterion function has already been evaluated, that + satisfies -1 <= existing_xs <= 1. + order (int): Type of norm to use when scaling the sampled points. For 2 we + project onto the hull of a sphere, for np.inf onto the hull of a cube. + hardness (float): Positive scaling factor. As hardness tends to infinity the + soft minimum (logsumexp) approaches the hard minimum. Default is 1. A + detailed explanation is given in the docstring. + n_params (int): Dimensionality of the problem. + + Returns: + float: The criterion value. + + """ + x = x.reshape(-1, n_params) + x = _project_onto_unit_hull(x, order=order) + + if existing_xs is not None: + sample = np.row_stack([x, existing_xs]) + n_existing_pairs = len(existing_xs) * (len(existing_xs) - 1) // 2 + slc = slice(0, -n_existing_pairs) if n_existing_pairs else slice(None) + else: + sample = x + slc = slice(None) + + dist = pdist(sample) ** 2 + + # drop distances between existing points. They could introduce flat spots. + dist = dist[slc] + + # soft minimum + crit_value = -logsumexp(-hardness * dist) + return crit_value + + +def _draw_from_distribution(distribution, rng, size): + """Draw points from distribution. + + Args: + distribution (str): Distribution to use for initial sample before points are + projected onto unit hull. Must be in {'normal', 'uniform'}. + rng (np.random.Generator): Random number generator. + size (Union[int, tuple[int]]): Output shape. + + Returns: + np.ndarray: Randomly drawn points. + + """ + if distribution == "normal": + draw = rng.normal(size=size) + elif distribution == "uniform": + draw = rng.uniform(-1, 1, size=size) + else: + raise ValueError( + f"distribution is {distribution}, but needs to be in ('normal', 'uniform')." + ) + return draw + + +def _map_into_feasible_trustregion(points, bounds): + """Map points from the unit space into trustregion defined by bounds. + + Args: + points (np.ndarray): 2d array of points to be mapped. Each value is in [-1, 1]. + bounds (Bounds): A NamedTuple with attributes ``lower`` and ``upper``, where + lower and upper define the rectangle that is the feasible trustregion. + + Returns: + np.ndarray: Points in trustregion. + + """ + out = (bounds.upper - bounds.lower) * (points + 1) / 2 + bounds.lower + return out + + +def _map_from_feasible_trustregion(points, bounds): + """Map points from a feasible trustregion definde by boudns into unit space. + + Args: + points (np.ndarray): 2d array of points to be mapped. Each value is in [-1, 1]. + bounds (Bounds): A NamedTuple with attributes ``lower`` and ``upper``, where + lower and upper define the rectangle that is the feasible trustregion. + + Returns: + np.ndarray: Points in unit space. + + """ + out = 2 * (points - bounds.lower) / (bounds.upper - bounds.lower) - 1 + return out + + +def _project_onto_unit_hull(x, order): + """Project points from the unit space onto the hull of a geometric figure. + + Args: + x (np.ndarray): 2d array of points to be projects. Each value is in [-1, 1]. + order (int): Type of norm to use when scaling the sampled points. For 2 we + project onto the hull of a sphere, for np.inf onto the hull of a cube. + + Returns: + np.ndarray: The projected points. + + """ + norm = np.linalg.norm(x, axis=1, ord=order).reshape(-1, 1) + projected = x / norm + return projected + + def _get_effective_bounds(trustregion, bounds): lower_bounds = trustregion.center - trustregion.radius upper_bounds = trustregion.center + trustregion.radius diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 9c9b6e3cf..129293ce7 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -14,6 +14,7 @@ from estimagic.optimization.tranquilo.models import ScalarModel from estimagic.optimization.tranquilo.options import Bounds from estimagic.optimization.tranquilo.options import ConvOptions +from estimagic.optimization.tranquilo.options import RadiusFactors from estimagic.optimization.tranquilo.options import RadiusOptions from estimagic.optimization.tranquilo.options import TrustRegion from estimagic.optimization.tranquilo.sample_points import get_sampler @@ -32,11 +33,12 @@ def _tranquilo( random_seed=925408, sampler="sphere", sample_filter="keep_all", - fitter="ols", + fitter=None, subsolver="bntr", sample_size=None, surrogate_model=None, radius_options=None, + radius_factors=None, sampler_options=None, fit_options=None, solver_options=None, @@ -110,6 +112,8 @@ def _tranquilo( if radius_options is None: radius_options = RadiusOptions() + if radius_factors is None: + radius_factors = RadiusFactors() if sampler_options is None: sampler_options = {} if fit_options is None: @@ -128,6 +132,12 @@ def _tranquilo( x=x, ) + if fitter is None: + if functype == "scalar": + fitter = "powell" + else: + fitter = "ols" + if functype == "scalar": aggregator = "identity" elif functype == "likelihood": @@ -147,7 +157,13 @@ def _tranquilo( history = History(functype=functype) bounds = Bounds(lower=lower_bounds, upper=upper_bounds) trustregion = TrustRegion(center=x, radius=radius_options.initial_radius) - sample_points = get_sampler(sampler, bounds=bounds, user_options=sampler_options) + sample_points = get_sampler( + sampler, + bounds=bounds, + model_info=model_info, + radius_factors=radius_factors, + user_options=sampler_options, + ) filter_points = get_sample_filter(sample_filter) aggregate_vector_model = get_aggregator( @@ -179,6 +195,7 @@ def _tranquilo( x=history.get_xs(0), fvec=history.get_fvecs(0), fval=history.get_fvals(0), + model_indices=np.array([0]), ) states = [state] @@ -255,6 +272,7 @@ def _tranquilo( x=candidate_x, fvec=history.get_fvecs(candidate_index), fval=history.get_fvals(candidate_index), + model_indices=model_indices, ) else: @@ -263,6 +281,7 @@ def _tranquilo( model=scalar_model, rho=rho, radius=new_radius, + model_indices=model_indices, ) states.append(state) @@ -303,6 +322,7 @@ class State(NamedTuple): x: np.ndarray fvec: np.ndarray fval: np.ndarray + model_indices: np.ndarray def _is_converged(states, options): diff --git a/src/estimagic/optimization/tranquilo/tranquilo_history.py b/src/estimagic/optimization/tranquilo/tranquilo_history.py index 51e870a9d..cdb8cc784 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo_history.py +++ b/src/estimagic/optimization/tranquilo/tranquilo_history.py @@ -41,7 +41,10 @@ def add_entries(self, xs, fvecs): least square fvecs. """ xs = np.atleast_2d(xs) - if len(xs) == 0: + + n_new_points = len(xs) if xs.size != 0 else 0 + + if n_new_points == 0: return if self.functype == "scalar": @@ -50,7 +53,7 @@ def add_entries(self, xs, fvecs): fvecs = np.atleast_2d(fvecs) fvals = np.atleast_1d(self.aggregate(fvecs)) - if len(xs) != len(fvecs): + if n_new_points != len(fvecs): raise ValueError() self.xs = _add_entries_to_array(self.xs, xs, self.n_fun) @@ -160,15 +163,17 @@ def _add_entries_to_array(arr, new, position): shape = 100_000 if new.ndim == 1 else (100_000, new.shape[1]) arr = np.full(shape, np.nan) - if len(arr) - position - len(new) < 0: - n_extend = max(len(arr), len(new)) + n_new_points = len(new) if new.size != 0 else 0 + + if len(arr) - position - n_new_points < 0: + n_extend = max(len(arr), n_new_points) if arr.ndim == 2: extension_shape = (n_extend, arr.shape[1]) arr = np.vstack([arr, np.full(extension_shape, np.nan)]) else: arr = np.hstack([arr, np.full(n_extend, np.nan)]) - arr[position : position + len(new)] = new + arr[position : position + n_new_points] = new return arr diff --git a/tests/optimization/tranquilo/test_filter_points.py b/tests/optimization/tranquilo/test_filter_points.py index f99ebb56b..cd7336dbb 100644 --- a/tests/optimization/tranquilo/test_filter_points.py +++ b/tests/optimization/tranquilo/test_filter_points.py @@ -49,6 +49,7 @@ def basic_case(): x=x_accepted, fvec=0, fval=0, + model_indices=None, ) expected_indices = np.array([20, 19, 18, 17, 16, 15, 13, 12, 8, 5, 4, 3, 2, 1, 0]) diff --git a/tests/optimization/tranquilo/test_sample_points.py b/tests/optimization/tranquilo/test_sample_points.py index 28cab7e21..b6b471ecc 100644 --- a/tests/optimization/tranquilo/test_sample_points.py +++ b/tests/optimization/tranquilo/test_sample_points.py @@ -1,22 +1,146 @@ -from collections import namedtuple - import numpy as np +import pytest +from estimagic.optimization.tranquilo.options import Bounds +from estimagic.optimization.tranquilo.options import RadiusFactors from estimagic.optimization.tranquilo.options import TrustRegion +from estimagic.optimization.tranquilo.sample_points import _draw_from_distribution +from estimagic.optimization.tranquilo.sample_points import ( + _minimal_pairwise_distance_on_hull, +) +from estimagic.optimization.tranquilo.sample_points import _project_onto_unit_hull from estimagic.optimization.tranquilo.sample_points import get_sampler +from numpy.testing import assert_array_almost_equal as aaae +from scipy.spatial.distance import pdist + + +SAMPLERS = ["naive", "cube", "sphere", "optimal_cube", "optimal_sphere"] + + +@pytest.mark.parametrize("sampler", SAMPLERS) +def test_bounds_are_satisfied(sampler): + bounds = Bounds(lower=-2 * np.ones(2), upper=np.array([0.25, 0.5])) + sampler = get_sampler(sampler, bounds) + sample = sampler( + trustregion=TrustRegion(center=np.zeros(2), radius=1), + target_size=5, + rng=np.random.default_rng(1234), + ) + lower = np.full_like(sample, bounds.lower) + upper = np.full_like(sample, bounds.upper) + assert np.all(lower <= sample) + assert np.all(sample <= upper) + + +@pytest.mark.parametrize("order", [3, 10, 100]) +def test_bounds_are_satisfied_general_hull_sampler(order): + bounds = Bounds(lower=-2 * np.ones(2), upper=np.array([0.25, 0.5])) + sampler = get_sampler("hull_sampler", bounds, user_options={"order": order}) + sample = sampler( + trustregion=TrustRegion(center=np.zeros(2), radius=1), + target_size=5, + rng=np.random.default_rng(1234), + ) + lower = np.full_like(sample, bounds.lower) + upper = np.full_like(sample, bounds.upper) + assert np.all(lower <= sample) + assert np.all(sample <= upper) -def test_integration_of_get_sampler_and_refercen_sampler(): +@pytest.mark.parametrize("sampler", SAMPLERS) +def test_enough_existing_points(sampler): + # test that if enough existing points exist an empty array is returned sampler = get_sampler( - sampler="naive", - bounds=namedtuple("Bounds", ["lower", "upper"])(-np.ones(3), np.ones(3)), + sampler=sampler, + bounds=Bounds(lower=-np.ones(3), upper=np.ones(3)), + ) + calculated = sampler( + trustregion=TrustRegion(center=np.zeros(3), radius=1), + target_size=5, + existing_xs=np.empty((5, 3)), + rng=np.random.default_rng(1234), ) + assert calculated.size == 0 + + +@pytest.mark.parametrize("sampler", ["optimal_cube", "optimal_sphere"]) +def test_optimization_ignores_existing_points(sampler): + # test that existing points behave as constants in the optimal sampling + sampler = get_sampler( + sampler=sampler, + bounds=Bounds(lower=-np.ones(3), upper=np.ones(3)), + model_info=None, + radius_factors=RadiusFactors(), + ) calculated = sampler( - trustregion=TrustRegion(center=0.5 * np.ones(3), radius=1), + trustregion=TrustRegion(center=np.zeros(3), radius=1), target_size=5, + existing_xs=np.ones((2, 3)), # same point implies min distance of zero always rng=np.random.default_rng(1234), ) - assert calculated.shape == (5, 3) - assert (calculated <= 1).all() - assert (calculated >= -1).all() + assert pdist(calculated).min() > 0 + + +@pytest.mark.parametrize("sampler", ["sphere", "cube"]) +def test_optimality(sampler): + # test that optimal versions of hull samplers produce better criterion value + standard_sampler = get_sampler( + sampler=sampler, + bounds=Bounds(lower=-np.ones(3), upper=np.ones(3)), + ) + optimal_sampler = get_sampler( + sampler="optimal_" + sampler, + bounds=Bounds(lower=-np.ones(3), upper=np.ones(3)), + ) + + distances = [] + for sampler in [standard_sampler, optimal_sampler]: + sample = sampler( + trustregion=TrustRegion(center=np.zeros(3), radius=1), + target_size=5, + rng=np.random.default_rng(1234), + ) + distances.append(pdist(sample).min()) + + assert distances[1] > distances[0] + + +@pytest.mark.parametrize("order", [2, np.inf]) +def test_pairwise_distance_on_hull(order): + + # equal points imply zero distance + value = _minimal_pairwise_distance_on_hull( + x=np.ones(4), existing_xs=None, hardness=1, order=order, n_params=2 + ) + assert value == 0 + + # non-equal points imply positive distance + value = _minimal_pairwise_distance_on_hull( + x=np.arange(4), existing_xs=None, hardness=1, order=order, n_params=2 + ) + assert value > 0 + + +@pytest.mark.parametrize("order", [2, np.inf]) +def test_project_onto_unit_hull(order): + + rng = np.random.default_rng(1234) + old = rng.uniform(-1, 1, size=10).reshape(5, 2) + new = _project_onto_unit_hull(old, order) + + norm = np.linalg.norm(old, axis=1, ord=order) + with pytest.raises(AssertionError): + aaae(1, norm) + + norm = np.linalg.norm(new, axis=1, ord=order) + aaae(1, norm) + + +@pytest.mark.parametrize("distribution", ["normal", "uniform"]) +def test_draw_from_distribution(distribution): + rng = np.random.default_rng() + draw = _draw_from_distribution(distribution, rng=rng, size=(3, 2)) + if distribution == "uniform": + assert (-1 <= draw).all() and (draw <= 1).all() + assert draw.shape == (3, 2)