From cc2b17ae08e389cd6951f3858e065347fe7c54ee Mon Sep 17 00:00:00 2001 From: timmens Date: Tue, 7 Mar 2023 13:57:33 +0100 Subject: [PATCH 1/8] Add attribute to Bounds --- src/estimagic/optimization/tranquilo/options.py | 16 +++++++++++++++- .../optimization/tranquilo/solve_subproblem.py | 11 +++++------ .../optimization/tranquilo/tranquilo.py | 14 ++------------ 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/options.py b/src/estimagic/optimization/tranquilo/options.py index d4eb2b4b0..4fef04fcc 100644 --- a/src/estimagic/optimization/tranquilo/options.py +++ b/src/estimagic/optimization/tranquilo/options.py @@ -1,3 +1,4 @@ +from dataclasses import dataclass from typing import NamedTuple import numpy as np @@ -8,12 +9,16 @@ ) -class Bounds(NamedTuple): +@dataclass +class Bounds: """Stopping criteria.""" lower: np.ndarray upper: np.ndarray + def __post_init__(self): + self.has_any = _check_if_there_are_bounds(self.lower, self.upper) + class StopOptions(NamedTuple): """Criteria for stopping without successful convergence.""" @@ -76,3 +81,12 @@ class StagnationOptions(NamedTuple): sample_increment: int = 1 max_trials: int = 1 drop: bool = True + + +def _check_if_there_are_bounds(lb, ub): + out = False + if lb is not None and np.isfinite(lb).any(): + out = True + if ub is not None and np.isfinite(ub).any(): + out = True + return out diff --git a/src/estimagic/optimization/tranquilo/solve_subproblem.py b/src/estimagic/optimization/tranquilo/solve_subproblem.py index 576e8e617..360bfdb09 100644 --- a/src/estimagic/optimization/tranquilo/solve_subproblem.py +++ b/src/estimagic/optimization/tranquilo/solve_subproblem.py @@ -56,7 +56,7 @@ def get_subsolver(solver, user_options=None, bounds=None): subproblem ("gqtpar"). - k_hard (float): Stopping criterion for the "hard" case in the trust-region subproblem ("gqtpar"). - bounds (NamedTuple or None): + bounds (Bounds): Returns: callable: The subsolver. @@ -108,12 +108,11 @@ def get_subsolver(solver, user_options=None, bounds=None): valid_bounds = {"lower_bounds", "upper_bounds"}.intersection(args) bounds_dict = {"lower_bounds": None, "upper_bounds": None} - if bounds is not None: + if bounds is not None and bounds.has_any: for type_ in ["lower", "upper"]: - if hasattr(bounds, type_): - candidate = getattr(bounds, type_) - if candidate is not None and np.isfinite(candidate).any(): - bounds_dict[f"{type_}_bounds"] = candidate + candidate = getattr(bounds, type_) + if candidate is not None and np.isfinite(candidate).any(): + bounds_dict[f"{type_}_bounds"] = candidate for name, value in bounds_dict.items(): if name not in valid_bounds and value is not None: diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 936137205..f9207f268 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -158,10 +158,10 @@ def _tranquilo( functype=functype, ) - _has_bounds = _check_if_there_are_bounds(lower_bounds, upper_bounds) + bounds = Bounds(lower=lower_bounds, upper=upper_bounds) if trustregion_shape is None: - trustregion_shape = "sphere" if not _has_bounds else "cube" + trustregion_shape = "sphere" if not bounds.has_any else "cube" if sampler is None: sampler = f"optimal_{trustregion_shape}" @@ -223,7 +223,6 @@ def _tranquilo( history=history, ) - bounds = Bounds(lower=lower_bounds, upper=upper_bounds) sample_points = get_sampler( sampler, bounds=bounds, @@ -734,15 +733,6 @@ def _process_sample_size(user_sample_size, model_info, x, sample_size_factor): ) -def _check_if_there_are_bounds(lb, ub): - out = False - if lb is not None and np.isfinite(lb).any(): - out = True - if ub is not None and np.isfinite(ub).any(): - out = True - return out - - def _get_search_region(trustregion, search_options): shape = trustregion.shape dim = len(trustregion.center) From fc15618bbfc2e7cb9370ba9135a23f4fa5fe6855 Mon Sep 17 00:00:00 2001 From: timmens Date: Tue, 7 Mar 2023 14:23:41 +0100 Subject: [PATCH 2/8] Reduce search options to search_radius_factor --- .../optimization/tranquilo/options.py | 5 --- .../optimization/tranquilo/tranquilo.py | 33 ++++--------------- 2 files changed, 7 insertions(+), 31 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/options.py b/src/estimagic/optimization/tranquilo/options.py index 4fef04fcc..09150241c 100644 --- a/src/estimagic/optimization/tranquilo/options.py +++ b/src/estimagic/optimization/tranquilo/options.py @@ -60,11 +60,6 @@ class Region(NamedTuple): shape: str -class HistorySearchOptions(NamedTuple): - radius_type: str = "inscribed" - radius_factor: float = 5 - - class AcceptanceOptions(NamedTuple): confidence_level: float = 0.8 power_level: float = 0.8 diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index f9207f268..3f999eb91 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -26,7 +26,6 @@ AcceptanceOptions, Bounds, ConvOptions, - HistorySearchOptions, RadiusOptions, Region, StagnationOptions, @@ -62,7 +61,7 @@ def _tranquilo( n_cores=1, silence_experimental_warning=False, infinity_handling="relative", - history_search_options=None, + search_radius_factor=None, noisy=False, sample_size_factor=None, acceptance_decider=None, @@ -104,7 +103,7 @@ def _tranquilo( - "quadratic: 0.5 * n * (n + 1) + n + 1 surrogate_model (str): Type of surrogate model to fit. Both a "linear" and "quadratic" surrogate model are supported. - radius_options (NemdTuple or NoneType): Options for trust-region radius + radius_options (NamedTuple or NoneType): Options for trust-region radius management. sampler_options (dict or NoneType): Additional keyword arguments passed to the sampler function. @@ -169,6 +168,9 @@ def _tranquilo( if subsolver is None: subsolver = "bntr_fast" if trustregion_shape == "cube" else "gqtpar_fast" + if search_radius_factor is None: + search_radius_factor = 4.25 if functype == "scalar" else 5.0 + target_sample_size = _process_sample_size( user_sample_size=sample_size, model_info=model_info, @@ -194,12 +196,6 @@ def _tranquilo( if conv_options is None: conv_options = ConvOptions() - if history_search_options is None: - if functype == "scalar": - history_search_options = HistorySearchOptions(radius_factor=4.25) - else: - history_search_options = HistorySearchOptions() - if acceptance_decider is None: acceptance_decider = "noisy" if noisy else "classic" @@ -309,9 +305,8 @@ def _tranquilo( # find, filter and count points # ============================================================================== - search_region = _get_search_region( - trustregion=state.trustregion, - search_options=history_search_options, + search_region = state.trustregion._replace( + radius=search_radius_factor * state.trustregion.radius ) old_indices = history.get_x_indices_in_region(search_region) @@ -733,20 +728,6 @@ def _process_sample_size(user_sample_size, model_info, x, sample_size_factor): ) -def _get_search_region(trustregion, search_options): - shape = trustregion.shape - dim = len(trustregion.center) - - if shape == "sphere" and search_options.radius_type == "circumscribed": - radius_factor = np.sqrt(dim) * search_options.radius_factor - else: - radius_factor = search_options.radius_factor - - search_radius = radius_factor * trustregion.radius - - return trustregion._replace(radius=search_radius) - - def _concatenate_indices(first, second): first = np.atleast_1d(first).astype(int) second = np.atleast_1d(second).astype(int) From 97813bcf9fe0a306bfc898c44588acdea2a217a9 Mon Sep 17 00:00:00 2001 From: timmens Date: Tue, 7 Mar 2023 17:01:53 +0100 Subject: [PATCH 3/8] Replace ModelInfo class by model_type string --- .../tranquilo/aggregate_models.py | 47 ++--- .../optimization/tranquilo/fit_models.py | 138 ++++++--------- .../optimization/tranquilo/models.py | 34 +--- .../optimization/tranquilo/tranquilo.py | 40 ++--- .../optimization/tranquilo/test_fit_models.py | 164 ++---------------- tests/optimization/tranquilo/test_models.py | 39 ++--- .../optimization/tranquilo/test_rho_noise.py | 13 +- .../optimization/tranquilo/test_tranquilo.py | 39 ++--- 8 files changed, 145 insertions(+), 369 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/aggregate_models.py b/src/estimagic/optimization/tranquilo/aggregate_models.py index 963d16340..b95efa42d 100644 --- a/src/estimagic/optimization/tranquilo/aggregate_models.py +++ b/src/estimagic/optimization/tranquilo/aggregate_models.py @@ -5,7 +5,7 @@ from estimagic.optimization.tranquilo.models import ScalarModel -def get_aggregator(aggregator, functype, model_info): +def get_aggregator(aggregator, functype, model_type): """Get a function that aggregates a VectorModel into a ScalarModel. Args: @@ -13,8 +13,9 @@ def get_aggregator(aggregator, functype, model_info): The function must take as argument: - vector_model (VectorModel): A fitted vector model. functype (str): One of "scalar", "least_squares" and "likelihood". - model_info (ModelInfo): Information that describes the functional form of - the model. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: callable: The partialled aggregator that only depends on vector_model. @@ -37,11 +38,11 @@ def get_aggregator(aggregator, functype, model_info): _using_built_in_aggregator = False else: raise ValueError( - "Invalid aggregator: {aggregator}. Must be one of " + "Invalid aggregator: {aggregator}. Must be one of " f"{list(built_in_aggregators)} or a callable." ) - # determine if aggregator is compatible with functype and model_info + # determine if aggregator is compatible with functype and model_type aggregator_compatible_with_functype = { "scalar": ("identity", "sum"), "least_squares": ("least_squares_linear",), @@ -51,17 +52,9 @@ def get_aggregator(aggregator, functype, model_info): ), } - aggregator_compatible_with_model_info = { - # keys are names of aggregators and values are functions of model_info that - # return False in case of incompatibility - "identity": lambda x: True, # noqa: ARG005 - "sum": _is_second_order_model, - "information_equality_linear": lambda model_info: not _is_second_order_model( - model_info - ), - "least_squares_linear": lambda model_info: not _is_second_order_model( - model_info - ), + aggregator_compatible_with_model_type = { + "linear": {"information_equality_linear", "least_squares_linear"}, + "quadratic": {"identity", "sum"}, } if _using_built_in_aggregator: @@ -71,12 +64,12 @@ def get_aggregator(aggregator, functype, model_info): f"Aggregator {_aggregator_name} is not compatible with functype " f"{functype}. It would not produce a quadratic main model." ) - if not aggregator_compatible_with_model_info[_aggregator_name](model_info): + if _aggregator_name not in aggregator_compatible_with_model_type[model_type]: raise ValueError( - f"ModelInfo {model_info} is not compatible with aggregator " - f"{_aggregator_name}. Depending on the aggregator this may be because " - "it would not produce a quadratic main model or that the aggregator " - "requires a different residual model for theoretical reasons." + f"Aggregator {_aggregator_name} is not compatible with model_type " + f"{model_type}. This is because the combination would not produce a " + "quadratic main model or that the aggregator requires a different " + "residual model." ) # create aggregator @@ -114,7 +107,7 @@ def aggregator_identity(vector_model): Assumptions ----------- 1. functype: scalar - 2. ModelInfo: has squares or interactions + 2. model_type: quadratic """ intercept = float(vector_model.intercepts) @@ -137,7 +130,7 @@ def aggregator_sum(vector_model): Assumptions ----------- 1. functype: likelihood - 2. ModelInfo: has squares or interactions + 2. model_type: quadratic """ vm_intercepts = vector_model.intercepts @@ -157,7 +150,7 @@ def aggregator_least_squares_linear(vector_model): Assumptions ----------- 1. functype: least_squares - 2. ModelInfo: has intercept but no squares and no interaction + 2. model_type: linear References ---------- @@ -185,7 +178,7 @@ def aggregator_information_equality_linear(vector_model): Assumptions ----------- 1. functype: likelihood - 2. ModelInfo: has no squares and no interaction + 2. model_type: linear """ vm_linear_terms = vector_model.linear_terms @@ -198,7 +191,3 @@ def aggregator_information_equality_linear(vector_model): square_terms = -fisher_information / 2 return intercept, linear_terms, square_terms - - -def _is_second_order_model(model_info): - return model_info.has_squares or model_info.has_interactions diff --git a/src/estimagic/optimization/tranquilo/fit_models.py b/src/estimagic/optimization/tranquilo/fit_models.py index 04dda5ec5..c4cb15615 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -7,44 +7,36 @@ from estimagic.optimization.tranquilo.get_component import get_component from estimagic.optimization.tranquilo.handle_infinity import get_infinity_handler from estimagic.optimization.tranquilo.models import ( - ModelInfo, VectorModel, add_models, move_model, - n_interactions, n_second_order_terms, ) def get_fitter( - fitter, fitter_options=None, model_info=None, infinity_handling="relative" + fitter, fitter_options=None, model_type=None, infinity_handling="relative" ): """Get a fit-function with partialled options. Args: fitter (str or callable): Name of a fit method or a fit method. The first argument of any fit method needs to be ``x``, second ``y`` and third - ``model_info``. + ``model_type``. user_options (dict): Options for the fit method. The following are supported: - l2_penalty_linear (float): Penalty that is applied to all linear terms. - l2_penalty_square (float): Penalty that is applied to all square terms, that is the quadratic and interaction terms. - model_info (ModelInfo): Information that describes the functional form of - the model. Has entries: - - has_squares (bool): Whether to use quadratic terms as features in the - regression. - - has_interactions (bool): Whether to use interaction terms as features - in the regression. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: callable: The partialled fit method that only depends on x and y. """ - if model_info is None: - model_info = ModelInfo() - fitter_options = {} if fitter_options is None else fitter_options built_in_fitters = { @@ -57,13 +49,13 @@ def get_fitter( default_options = { "l2_penalty_linear": 0, "l2_penalty_square": 0.1, - "model_info": model_info, + "model_type": model_type, "p_intercept": 0.05, "p_linear": 0.4, "p_square": 1.0, } - mandatory_arguments = ["x", "y", "model_info"] + mandatory_arguments = ["x", "y", "model_type"] _raw_fitter = get_component( name_or_func=fitter, @@ -79,7 +71,7 @@ def get_fitter( fitter = partial( _fitter_template, fitter=_raw_fitter, - model_info=model_info, + model_type=model_type, clip_infinite_values=clip_infinite_values, residualize=fitter_options.get("residualize", False), ) @@ -94,7 +86,7 @@ def _fitter_template( old_model, weights=None, fitter=None, - model_info=None, + model_type=None, clip_infinite_values=None, residualize=False, ): @@ -107,9 +99,10 @@ def _fitter_template( evaluations that have been centered around the function value at the trust region center. fitter (callable): Fit method. The first argument of any fit method needs to be - ``x``, second ``y`` and third ``model_info``. - model_info (ModelInfo): Information that describes the functional form of - the model. + ``x``, second ``y`` and third ``model_type``. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: VectorModel or ScalarModel: Results container. @@ -134,12 +127,10 @@ def _fitter_template( intercepts = intercepts.flatten() # construct final square terms - if model_info.has_interactions: + if model_type == "quadratic": square_terms = _reshape_square_terms_to_hess( - square_terms, n_params, n_residuals, model_info.has_squares + square_terms, n_params, n_residuals ) - elif model_info.has_squares: - square_terms = 2 * np.stack([np.diag(a) for a in square_terms]) else: square_terms = None @@ -151,7 +142,7 @@ def _fitter_template( return results -def fit_ols(x, y, weights, model_info): +def fit_ols(x, y, weights, model_type): """Fit a linear model using ordinary least squares. Args: @@ -160,14 +151,15 @@ def fit_ols(x, y, weights, model_info): y (np.ndarray): Array of shape (n_samples, n_residuals) with function evaluations that have been centered around the function value at the trust region center. - model_info (ModelInfo): Information that describes the functional form of the - model. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: np.ndarray: The model coefficients. """ - features = _build_feature_matrix(x, model_info) + features = _build_feature_matrix(x, model_type) features_w, y_w = _add_weighting(features, y, weights) coef = _fit_ols(features_w, y_w) @@ -191,7 +183,7 @@ def _fit_ols(x, y): return coef -def fit_tranquilo(x, y, weights, model_info, p_intercept, p_linear, p_square): +def fit_tranquilo(x, y, weights, model_type, p_intercept, p_linear, p_square): """Fit a linear model using ordinary least squares. The difference to fit_ols is that the linear terms are penalized less strongly @@ -203,14 +195,15 @@ def fit_tranquilo(x, y, weights, model_info, p_intercept, p_linear, p_square): y (np.ndarray): Array of shape (n_samples, n_residuals) with function evaluations that have been centered around the function value at the trust region center. - model_info (ModelInfo): Information that describes the functional form of the - model. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: np.ndarray: The model coefficients. """ - features = _build_feature_matrix(x, model_info) + features = _build_feature_matrix(x, model_type) features_w, y_w = _add_weighting(features, y, weights) n_params = x.shape[1] @@ -232,7 +225,7 @@ def fit_ridge( x, y, weights, - model_info, + model_type, l2_penalty_linear, l2_penalty_square, ): @@ -244,8 +237,9 @@ def fit_ridge( y (np.ndarray): Array of shape (n_samples, n_residuals) with function evaluations that have been centered around the function value at the trust region center. - model_info (ModelInfo): Information that describes the functional form of the - model. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. l2_penalty_linear (float): Penalty that is applied to all linear terms. l2_penalty_square (float): Penalty that is applied to all square terms, that is the quadratic and interaction terms. @@ -254,7 +248,7 @@ def fit_ridge( np.ndarray: The model coefficients. """ - features = _build_feature_matrix(x, model_info) + features = _build_feature_matrix(x, model_type) features_w, y_w = _add_weighting(features, y, weights) @@ -293,7 +287,7 @@ def _fit_ridge(x, y, penalty): return coef -def fit_powell(x, y, model_info): +def fit_powell(x, y, model_type): """Fit a model, switching between penalized and unpenalized fitting. For: @@ -309,8 +303,9 @@ def fit_powell(x, y, model_info): y (np.ndarray): Array of shape (n_samples, n_residuals) with function evaluations that have been centered around the function value at the trust region center. - model_info (ModelInfo): Information that describes the functional form of the - model. + model_type (str): Type of the model that is fitted. The following are supported: + - "linear": Only linear effects and intercept. + - "quadratic": Fully quadratic model. Returns: np.ndarray: The model coefficients. @@ -321,26 +316,23 @@ def fit_powell(x, y, model_info): _switch_to_linear = n_samples <= n_params + 1 _n_just_identified = n_params + 1 - if model_info.has_squares: - _n_just_identified += n_params - if model_info.has_interactions: - _n_just_identified += int(0.5 * n_params * (n_params - 1)) + if model_type == "quadratic": + _n_just_identified += n_second_order_terms(n_params) if _switch_to_linear: - model_info = model_info._replace(has_squares=False, has_interactions=False) - coef = fit_ols(x, y, weights=None, model_info=model_info) + coef = fit_ols(x, y, weights=None, model_type="linear") n_resid, n_present = coef.shape padding = np.zeros((n_resid, _n_just_identified - n_present)) coef = np.hstack([coef, padding]) elif n_samples >= _n_just_identified: - coef = fit_ols(x, y, weights=None, model_info=model_info) + coef = fit_ols(x, y, weights=None, model_type=model_type) else: - coef = _fit_minimal_frobenius_norm_of_hessian(x, y, model_info) + coef = _fit_minimal_frobenius_norm_of_hessian(x, y, model_type) return coef -def _fit_minimal_frobenius_norm_of_hessian(x, y, model_info): +def _fit_minimal_frobenius_norm_of_hessian(x, y): """Fit a quadraitc model using the powell fitting method. The solution represents the quadratic whose Hessian matrix is of @@ -360,8 +352,6 @@ def _fit_minimal_frobenius_norm_of_hessian(x, y, model_info): y (np.ndarray): Array of shape (n_samples, n_residuals) with function evaluations that have been centered around the function value at the trust region center. - model_info (ModelInfo): Information that describes the functional form of the - model. Returns: np.ndarray: The model coefficients. @@ -377,19 +367,14 @@ def _fit_minimal_frobenius_norm_of_hessian(x, y, model_info): if n_samples >= _n_too_many: raise ValueError("Too may points for minimum frobenius fitting") - has_squares = model_info.has_squares - - if has_squares: - n_poly_features = n_params * (n_params + 1) // 2 - else: - n_poly_features = n_params * (n_params - 1) // 2 + n_poly_features = n_second_order_terms(n_params) ( m_mat, n_mat, z_mat, n_z_mat, - ) = _get_feature_matrices_minimal_frobenius_norm_of_hessian(x, model_info) + ) = _get_feature_matrices_minimal_frobenius_norm_of_hessian(x) coef = _get_current_fit_minimal_frobenius_norm_of_hessian( y=y, @@ -441,11 +426,11 @@ def _get_current_fit_minimal_frobenius_norm_of_hessian( return np.atleast_2d(coef) -def _get_feature_matrices_minimal_frobenius_norm_of_hessian(x, model_info): +def _get_feature_matrices_minimal_frobenius_norm_of_hessian(x): n_samples, n_params = x.shape - has_squares = model_info.has_squares - features = _polynomial_features(x, has_squares) + intercept = np.ones((n_samples, 1)) + features = np.concatenate((intercept, _quadratic_features(x)), axis=1) m_mat, n_mat = np.split(features, (n_params + 1,), axis=1) m_mat_pad = np.zeros((n_samples, n_samples)) @@ -469,20 +454,15 @@ def _get_feature_matrices_minimal_frobenius_norm_of_hessian(x, model_info): ) -def _build_feature_matrix(x, model_info): - if model_info.has_interactions: - features = _polynomial_features(x, model_info.has_squares) - else: - data = (np.ones(len(x)), x) - data = (*data, x**2) if model_info.has_squares else data - features = np.column_stack(data) - +def _build_feature_matrix(x, model_type): + raw = x if model_type == "linear" else _quadratic_features(x) + intercept = np.ones((len(x), 1)) + features = np.concatenate((intercept, raw), axis=1) return features -def _reshape_square_terms_to_hess(square_terms, n_params, n_residuals, has_squares): - offset = 0 if has_squares else 1 - idx1, idx2 = np.triu_indices(n_params, k=offset) +def _reshape_square_terms_to_hess(square_terms, n_params, n_residuals): + idx1, idx2 = np.triu_indices(n_params) hess = np.zeros((n_residuals, n_params, n_params), dtype=np.float64) hess[:, idx1, idx2] = square_terms hess = hess + np.triu(hess).transpose(0, 2, 1) @@ -491,27 +471,21 @@ def _reshape_square_terms_to_hess(square_terms, n_params, n_residuals, has_squar @njit -def _polynomial_features(x, has_squares): +def _quadratic_features(x): + # Create fully quadratic features without intercept n_samples, n_params = x.shape - - if has_squares: - n_poly_terms = n_second_order_terms(n_params) - else: - n_poly_terms = n_interactions(n_params) + n_poly_terms = n_second_order_terms(n_params) poly_terms = np.empty((n_poly_terms, n_samples), np.float64) xt = x.T idx = 0 for i in range(n_params): - j_start = i if has_squares else i + 1 + j_start = i for j in range(j_start, n_params): poly_terms[idx] = xt[i] * xt[j] idx += 1 - - intercept = np.ones((1, n_samples), x.dtype) - out = np.concatenate((intercept, xt, poly_terms), axis=0) - + out = np.concatenate((xt, poly_terms), axis=0) return out.T diff --git a/src/estimagic/optimization/tranquilo/models.py b/src/estimagic/optimization/tranquilo/models.py index a9a6dbe81..c9a4d03d2 100644 --- a/src/estimagic/optimization/tranquilo/models.py +++ b/src/estimagic/optimization/tranquilo/models.py @@ -1,5 +1,5 @@ from dataclasses import dataclass, replace -from typing import NamedTuple, Union +from typing import Union import numpy as np from numba import njit @@ -41,11 +41,6 @@ def _replace(self, **kwargs): return replace(self, **kwargs) -class ModelInfo(NamedTuple): - has_squares: bool = True - has_interactions: bool = True - - def _predict_vector(model: VectorModel, centered_x: np.ndarray) -> np.ndarray: """Evaluate a VectorModel at centered_x. @@ -308,25 +303,12 @@ def _predict_scalar(model: ScalarModel, centered_x: np.ndarray) -> np.ndarray: return out -def n_free_params(dim, info_or_name): +def n_free_params(dim, model_type): """Number of free parameters in a model specified by name or model_info.""" out = dim + 1 - if isinstance(info_or_name, ModelInfo): - info = info_or_name - if info.has_squares: - out += dim - if info.has_interactions: - out += n_interactions(dim) - elif isinstance(info_or_name, str) and info_or_name in ( - "linear", - "quadratic", - "diagonal", - ): - name = info_or_name - if name == "quadratic": + if model_type in ("linear", "quadratic"): + if model_type == "quadratic": out += n_second_order_terms(dim) - elif name == "diagonal": - out += dim else: raise ValueError() return out @@ -346,12 +328,10 @@ def n_interactions(dim): def is_second_order_model(model_or_info): """Check if a model has any second order terms.""" - if isinstance(model_or_info, ModelInfo): - model_info = model_or_info - out = model_info.has_interactions or model_info.has_squares + if isinstance(model_or_info, str): + out = model_or_info == "quadratic" elif isinstance(model_or_info, (ScalarModel, VectorModel)): - model = model_or_info - out = model.square_terms is not None + out = model_or_info.square_terms is not None else: raise TypeError() return out diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 3f999eb91..4dfacf7ca 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -16,7 +16,6 @@ ) from estimagic.optimization.tranquilo.fit_models import get_fitter from estimagic.optimization.tranquilo.models import ( - ModelInfo, ScalarModel, VectorModel, n_free_params, @@ -152,7 +151,7 @@ def _tranquilo( if solver_options is None: solver_options = {} - model_info = _process_surrogate_model( + model_type = _process_surrogate_model( surrogate_model=surrogate_model, functype=functype, ) @@ -173,7 +172,7 @@ def _tranquilo( target_sample_size = _process_sample_size( user_sample_size=sample_size, - model_info=model_info, + model_type=model_type, x=x, sample_size_factor=sample_size_factor, ) @@ -222,7 +221,7 @@ def _tranquilo( sample_points = get_sampler( sampler, bounds=bounds, - model_info=model_info, + model_info=model_type, user_options=sampler_options, ) @@ -231,13 +230,13 @@ def _tranquilo( aggregate_vector_model = get_aggregator( aggregator=aggregator, functype=functype, - model_info=model_info, + model_type=model_type, ) fit_model = get_fitter( fitter=fitter, fitter_options=fit_options, - model_info=model_info, + model_type=model_type, infinity_handling=infinity_handling, ) @@ -663,26 +662,21 @@ def _process_surrogate_model(surrogate_model, functype): else: surrogate_model = "linear" - if isinstance(surrogate_model, ModelInfo): - out = surrogate_model - elif isinstance(surrogate_model, str): - if surrogate_model == "linear": - out = ModelInfo(has_squares=False, has_interactions=False) - elif surrogate_model == "diagonal": - out = ModelInfo(has_squares=True, has_interactions=False) - elif surrogate_model == "quadratic": - out = ModelInfo(has_squares=True, has_interactions=True) - else: - raise ValueError(f"Invalid surrogate model: {surrogate_model}") - + if isinstance(surrogate_model, str): + if surrogate_model not in ("linear", "quadratic"): + raise ValueError( + f"Invalid surrogate model: {surrogate_model} must be in ('linear', " + "'quadratic')" + ) else: raise TypeError(f"Invalid surrogate model: {surrogate_model}") - return out + + return surrogate_model -def _process_sample_size(user_sample_size, model_info, x, sample_size_factor): +def _process_sample_size(user_sample_size, model_type, x, sample_size_factor): if user_sample_size is None: - if model_info.has_squares or model_info.has_interactions: + if model_type == "quadratic": out = 2 * len(x) + 1 else: out = len(x) + 1 @@ -690,11 +684,11 @@ def _process_sample_size(user_sample_size, model_info, x, sample_size_factor): elif isinstance(user_sample_size, str): user_sample_size = user_sample_size.replace(" ", "") if user_sample_size in ["linear", "n+1"]: - out = n_free_params(dim=len(x), info_or_name="linear") + out = n_free_params(dim=len(x), model_type="linear") elif user_sample_size in ["powell", "2n+1", "2*n+1"]: out = 2 * len(x) + 1 elif user_sample_size == "quadratic": - out = n_free_params(dim=len(x), info_or_name="quadratic") + out = n_free_params(dim=len(x), model_type="quadratic") else: raise ValueError(f"Invalid sample size: {user_sample_size}") diff --git a/tests/optimization/tranquilo/test_fit_models.py b/tests/optimization/tranquilo/test_fit_models.py index 9ba9c23bb..9386a18a4 100644 --- a/tests/optimization/tranquilo/test_fit_models.py +++ b/tests/optimization/tranquilo/test_fit_models.py @@ -1,10 +1,7 @@ import numpy as np import pytest -import yaml from estimagic import first_derivative, second_derivative -from estimagic.config import TEST_FIXTURES_DIR -from estimagic.optimization.tranquilo.fit_models import _polynomial_features, get_fitter -from estimagic.optimization.tranquilo.models import ModelInfo +from estimagic.optimization.tranquilo.fit_models import _quadratic_features, get_fitter from estimagic.optimization.tranquilo.options import Region from numpy.testing import assert_array_almost_equal, assert_array_equal @@ -18,13 +15,6 @@ def aaae(x, y, case=None): assert_array_almost_equal(x, y, decimal=tolerance[case]) -def read_yaml(path): - with open(rf"{path}") as file: - data = yaml.full_load(file) - - return data - - # ====================================================================================== # Fixtures # ====================================================================================== @@ -38,7 +28,7 @@ def quadratic_case(): """ n_params = 4 - n_samples = 2_000 + n_samples = 15 # theoretical terms linear_terms = 1 + np.arange(n_params) @@ -52,46 +42,8 @@ def func(x): x0 = np.ones(n_params) # random data - x = np.array( - [x0 + np.random.uniform(-0.01 * x0, 0.01 * x0) for _ in range(n_samples)] - ) - y = np.array([func(_x) for _x in list(x)]).reshape(-1, 1) - - out = { - "func": func, - "x0": x0, - "x": x, - "y": y, - "linear_terms_expected": linear_terms, - "square_terms_expected": square_terms, - } - return out - - -@pytest.fixture() -def just_identified_case(): - """Test scenario with true quadratic function and n + 1 points. - - We return true function, and function evaluations and data on random points. - - """ - n_params = 4 - n_samples = n_params + 1 - - # theoretical terms - linear_terms = 1 + np.arange(n_params) - square_terms = np.zeros((n_params, n_params)) - - def func(x): - y = -10 + linear_terms @ x + 0.5 * x.T @ square_terms @ x - return y - - x0 = np.ones(n_params) - - # random data - x = np.array( - [x0 + np.random.uniform(-0.01 * x0, 0.01 * x0) for _ in range(n_samples)] - ) + rng = np.random.default_rng(56789) + x = np.array([x0 + rng.uniform(-0.01 * x0, 0.01 * x0) for _ in range(n_samples)]) y = np.array([func(_x) for _x in list(x)]).reshape(-1, 1) out = { @@ -105,82 +57,13 @@ def func(x): return out -@pytest.fixture() -def data_fit_pounders(): - """Test data from Tao Pounders.""" - test_data = read_yaml(TEST_FIXTURES_DIR / "get_coefficients_residual_model.yaml") - - n_params = 3 - n_samples = 2 * n_params + 1 - n_poly_features = n_params * (n_params + 1) // 2 - _is_just_identified = False - - inputs_dict = { - "y": np.array(test_data["f_interpolated"]), - "m_mat": np.array(test_data["x_sample_monomial_basis"])[ - : n_params + 1, : n_params + 1 - ], - "n_mat": np.array(test_data["monomial_basis"])[:n_samples], - "z_mat": np.array(test_data["basis_null_space"]), - "n_z_mat": np.array(test_data["lower_triangular"])[:, n_params + 1 : n_samples], - "n_params": n_params, - "n_poly_features": n_poly_features, - "_is_just_identified": _is_just_identified, - } - - expected = { - "linear_terms": np.array(test_data["linear_terms_expected"]), - "square_terms": np.array(test_data["square_terms_expected"]), - } - - return inputs_dict, expected - - -@pytest.fixture() -def data_get_feature_matrices_pounders(): - test_data = read_yaml( - TEST_FIXTURES_DIR / "get_interpolation_matrices_residual_model.yaml" - ) - - n_params = 3 - n_samples = 2 * n_params + 1 - n_poly_features = n_params * (n_params + 1) // 2 - center = np.array(test_data["x_accepted"]) - radius = test_data["delta"] - - model_indices = np.array([13, 12, 11, 10, 9, 8, 6]) - history_x = np.array(test_data["history_x"]) - x = (history_x[model_indices] - center) / radius - - inputs_dict = { - "x": x, - "model_indices": np.array(test_data["model_indices"]), - "n_params": n_params, - "n_samples": n_samples, - "n_poly_features": n_poly_features, - } - - expected = { - "m_mat": np.array(test_data["x_sample_monomial_basis_expected"])[ - : n_params + 1, : n_params + 1 - ], - "n_mat": np.array(test_data["monomial_basis_expected"]), - "z_mat": np.array(test_data["basis_null_space_expected"]), - "n_z_mat": np.array(test_data["lower_triangular_expected"])[ - :, n_params + 1 : n_samples - ], - } - - return inputs_dict, expected - - # ====================================================================================== # Tests # ====================================================================================== def test_fit_ols_against_truth(quadratic_case): - fit_ols = get_fitter("ols") + fit_ols = get_fitter("ols", model_type="quadratic") got = fit_ols( x=quadratic_case["x"], y=quadratic_case["y"], @@ -192,21 +75,17 @@ def test_fit_ols_against_truth(quadratic_case): aaae(got.square_terms.squeeze(), quadratic_case["square_terms_expected"]) -@pytest.mark.parametrize("scenario", ["just_identified_case", "quadratic_case"]) -def test_fit_powell_against_truth(scenario, request): - test_case = request.getfixturevalue(scenario) - - model_info = ModelInfo(has_squares=True, has_interactions=True) - fit_pounders = get_fitter("powell", model_info=model_info) +def test_fit_powell_against_truth_quadratic(quadratic_case): + fit_pounders = get_fitter("powell", model_type="quadratic") got = fit_pounders( - test_case["x"], - test_case["y"], + quadratic_case["x"], + quadratic_case["y"], region=Region(center=np.zeros(4), radius=1.0, shape="sphere"), old_model=None, ) - aaae(got.linear_terms.squeeze(), test_case["linear_terms_expected"]) - aaae(got.square_terms.squeeze(), test_case["square_terms_expected"]) + aaae(got.linear_terms.squeeze(), quadratic_case["linear_terms_expected"]) + aaae(got.square_terms.squeeze(), quadratic_case["square_terms_expected"]) @pytest.mark.parametrize("model", ["ols", "ridge"]) @@ -216,7 +95,7 @@ def test_fit_ols_against_gradient(model, quadratic_case): else: options = None - fit_ols = get_fitter(model, options) + fit_ols = get_fitter(model, options, model_type="quadratic") got = fit_ols( quadratic_case["x"], quadratic_case["y"], @@ -237,7 +116,7 @@ def test_fit_ols_against_gradient(model, quadratic_case): [("ols", None), ("ridge", {"l2_penalty_linear": 0, "l2_penalty_square": 0})], ) def test_fit_ols_against_hessian(model, options, quadratic_case): - fit_ols = get_fitter(model, options) + fit_ols = get_fitter(model, options, model_type="quadratic") got = fit_ols( quadratic_case["x"], quadratic_case["y"], @@ -249,18 +128,9 @@ def test_fit_ols_against_hessian(model, options, quadratic_case): aaae(hessian["derivative"], hess, case="hessian") -@pytest.mark.parametrize("has_squares", [True, False]) -def test_polynomial_features(has_squares): +def test_quadratic_features(): x = np.array([[0, 1, 2], [3, 4, 5]]) - expected = { - # has_squares: expected value, - True: np.array( - [[1, 0, 1, 2, 0, 0, 0, 1, 2, 4], [1, 3, 4, 5, 9, 12, 15, 16, 20, 25]] - ), - False: np.array([[1, 0, 1, 2, 0, 0, 2], [1, 3, 4, 5, 12, 15, 20]]), - } - - got = _polynomial_features(x, has_squares=has_squares) - - assert_array_equal(got, expected[has_squares]) + expected = np.array([[0, 1, 2, 0, 0, 0, 1, 2, 4], [3, 4, 5, 9, 12, 15, 16, 20, 25]]) + got = _quadratic_features(x) + assert_array_equal(got, expected) diff --git a/tests/optimization/tranquilo/test_models.py b/tests/optimization/tranquilo/test_models.py index 451a0ed87..daf255b85 100644 --- a/tests/optimization/tranquilo/test_models.py +++ b/tests/optimization/tranquilo/test_models.py @@ -1,7 +1,6 @@ import numpy as np import pytest from estimagic.optimization.tranquilo.models import ( - ModelInfo, ScalarModel, VectorModel, _predict_scalar, @@ -52,44 +51,32 @@ def test_predict_vector(): def test_n_free_params_name_quadratic(): - assert n_free_params(dim=2, info_or_name="quadratic") == 1 + 2 + 3 - assert n_free_params(dim=3, info_or_name="quadratic") == 1 + 3 + 6 - assert n_free_params(dim=9, info_or_name="quadratic") == 1 + 9 + 45 - - -def test_n_free_params_name_diagonal(): - assert n_free_params(dim=2, info_or_name="diagonal") == 1 + 2 + 2 - assert n_free_params(dim=3, info_or_name="diagonal") == 1 + 3 + 3 - assert n_free_params(dim=9, info_or_name="diagonal") == 1 + 9 + 9 + assert n_free_params(dim=2, model_type="quadratic") == 1 + 2 + 3 + assert n_free_params(dim=3, model_type="quadratic") == 1 + 3 + 6 + assert n_free_params(dim=9, model_type="quadratic") == 1 + 9 + 45 def test_n_free_params_name_invalid(): with pytest.raises(ValueError): - assert n_free_params(dim=3, info_or_name="invalid") + assert n_free_params(dim=3, model_type="invalid") @pytest.mark.parametrize("dim", [2, 3, 9]) def test_n_free_params_info_linear(dim): - info = ModelInfo(has_squares=False, has_interactions=False) - assert n_free_params(dim, info) == 1 + dim - - -@pytest.mark.parametrize("dim", [2, 3, 9]) -def test_n_free_params_info_diagonal(dim): - info = ModelInfo(has_squares=True, has_interactions=False) - assert n_free_params(dim, info) == 1 + dim + dim + assert n_free_params(dim, model_type="linear") == 1 + dim @pytest.mark.parametrize("dim", [2, 3, 9]) def test_n_free_params_info_quadratic(dim): - info = ModelInfo(has_squares=True, has_interactions=True) - assert n_free_params(dim, info) == 1 + dim + dim + (dim * (dim - 1) // 2) + assert n_free_params(dim, model_type="quadratic") == 1 + dim + n_second_order_terms( + dim + ) def test_n_free_params_invalid(): model = ScalarModel(intercept=1.0, linear_terms=np.ones(1), square_terms=np.ones(1)) with pytest.raises(ValueError): - n_free_params(dim=1, info_or_name=model) + n_free_params(dim=1, model_type=model) def test_n_second_order_terms(): @@ -100,11 +87,9 @@ def test_n_interactions(): assert n_interactions(3) == 3 -@pytest.mark.parametrize("has_squares", [True, False]) -@pytest.mark.parametrize("has_interactions", [True, False]) -def test_is_second_order_model_info(has_squares, has_interactions): - model_info = ModelInfo(has_squares=has_squares, has_interactions=has_interactions) - assert is_second_order_model(model_info) == has_squares or has_interactions +@pytest.mark.parametrize("model_type", ("linear", "quadratic")) +def test_is_second_order_model_type(model_type): + assert is_second_order_model(model_type) == (model_type == "quadratic") def test_is_second_order_model_model(): diff --git a/tests/optimization/tranquilo/test_rho_noise.py b/tests/optimization/tranquilo/test_rho_noise.py index 15612fbbd..872a1ce8b 100644 --- a/tests/optimization/tranquilo/test_rho_noise.py +++ b/tests/optimization/tranquilo/test_rho_noise.py @@ -2,7 +2,6 @@ import pytest from estimagic.optimization.tranquilo.aggregate_models import get_aggregator from estimagic.optimization.tranquilo.fit_models import get_fitter -from estimagic.optimization.tranquilo.models import ModelInfo from estimagic.optimization.tranquilo.options import Region from estimagic.optimization.tranquilo.rho_noise import simulate_rho_noise from estimagic.optimization.tranquilo.solve_subproblem import get_subsolver @@ -10,7 +9,7 @@ @pytest.mark.parametrize("functype", ["scalar", "least_squares"]) -def test_convergence_to_one_if_nois_is_tiny(functype): +def test_convergence_to_one_if_noise_is_tiny(functype): """Test simulate_rho_noise. For the test, the "true" model is a standard sphere function. @@ -31,27 +30,27 @@ def test_convergence_to_one_if_nois_is_tiny(functype): if functype == "least_squares": fvecs = xs.copy() - model_info = ModelInfo(False, False) + model_type = "linear" model_aggregator = get_aggregator( aggregator="least_squares_linear", functype="least_squares", - model_info=model_info, + model_type=model_type, ) n_residuals = 2 else: fvecs = (xs**2).sum(axis=1).reshape(-1, 1) - model_info = ModelInfo(True, True) + model_type = "quadratic" model_aggregator = get_aggregator( aggregator="identity", functype="scalar", - model_info=model_info, + model_type=model_type, ) n_residuals = 1 noise_cov = np.eye(n_residuals) * 1e-12 trustregion = Region(center=np.ones(2) * 0.5, radius=1, shape="sphere") - model_fitter = get_fitter(fitter="ols", model_info=model_info) + model_fitter = get_fitter(fitter="ols", model_type=model_type) vector_model = model_fitter( xs, fvecs, weights=None, region=trustregion, old_model=None diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index 1e23898cf..17940ec8b 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -3,7 +3,6 @@ import numpy as np import pytest from estimagic.optimization.optimize import minimize -from estimagic.optimization.tranquilo.models import ModelInfo from estimagic.optimization.tranquilo.tranquilo import ( _process_sample_size, _process_surrogate_model, @@ -192,41 +191,29 @@ def test_external_tranquilo_ls_sphere_defaults(): def test_process_surrogate_model_none_scalar(): got = _process_surrogate_model(None, functype="scalar") - assert got.has_interactions is True - assert got.has_squares is True + assert got == "quadratic" @pytest.mark.parametrize("functype", ["least_squares", "likelihood"]) def test_process_surrogate_model_none_not_scalar(functype): got = _process_surrogate_model(None, functype=functype) - assert got.has_interactions is False - assert got.has_squares is False + assert got == "linear" -@pytest.mark.parametrize("has_interactions", [True, False]) -@pytest.mark.parametrize("has_squares", [True, False]) -def test_process_surrogate_model_info(has_interactions, has_squares): - model_info = ModelInfo(has_squares=has_squares, has_interactions=has_interactions) - got = _process_surrogate_model(model_info, functype="whatever") - assert got == model_info +@pytest.mark.parametrize("model_type", ("linear", "quadratic")) +def test_process_surrogate_model_info(model_type): + got = _process_surrogate_model(model_type, functype="whatever") + assert got == model_type def test_process_surrogate_model_str_linear(): got = _process_surrogate_model("linear", functype="scalar") - assert got.has_interactions is False - assert got.has_squares is False - - -def test_process_surrogate_model_str_diagonal(): - got = _process_surrogate_model("diagonal", functype="least_squares") - assert got.has_interactions is False - assert got.has_squares is True + assert got == "linear" def test_process_surrogate_model_str_quadratic(): got = _process_surrogate_model("quadratic", functype="likelihood") - assert got.has_interactions is True - assert got.has_squares is True + assert got == "quadratic" def test_process_surrogate_model_str_invalid(): @@ -241,15 +228,13 @@ def test_process_surrogate_model_invalid(functype): _process_surrogate_model(surrogate_model, functype=functype) -@pytest.mark.parametrize("has_interactions", [True, False]) -@pytest.mark.parametrize("has_squares", [True, False]) -def test_process_sample_size_none_linear(has_interactions, has_squares): - model_info = ModelInfo(has_interactions=has_interactions, has_squares=has_squares) +@pytest.mark.parametrize("model_type", ("linear", "quadratic")) +def test_process_sample_size_none_linear(model_type): x = np.ones((3, 2)) got = _process_sample_size( - None, model_info=model_info, x=x, sample_size_factor=None + None, model_type=model_type, x=x, sample_size_factor=None ) - if has_interactions or has_squares: + if model_type == "quadratic": assert got == 7 else: assert got == 4 From aa779b85c1e42234b149196f974050d9144f8c76 Mon Sep 17 00:00:00 2001 From: timmens Date: Tue, 7 Mar 2023 17:44:22 +0100 Subject: [PATCH 4/8] Fix bug when calling _fit_minimal_frobenius_norm_of_hessian --- src/estimagic/optimization/tranquilo/fit_models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/estimagic/optimization/tranquilo/fit_models.py b/src/estimagic/optimization/tranquilo/fit_models.py index c4cb15615..89926921c 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -327,7 +327,7 @@ def fit_powell(x, y, model_type): elif n_samples >= _n_just_identified: coef = fit_ols(x, y, weights=None, model_type=model_type) else: - coef = _fit_minimal_frobenius_norm_of_hessian(x, y, model_type) + coef = _fit_minimal_frobenius_norm_of_hessian(x, y) return coef From 0dcbc52c3f7f941c7a39a62668f861e85d797c87 Mon Sep 17 00:00:00 2001 From: Tim Mensinger Date: Mon, 13 Mar 2023 17:32:14 +0100 Subject: [PATCH 5/8] Add functionality to Region class (#444) --- .../optimization/tranquilo/bounds.py | 27 ++++ .../tranquilo/estimate_variance.py | 2 +- .../optimization/tranquilo/geometry.py | 28 ++--- .../optimization/tranquilo/models.py | 2 +- .../optimization/tranquilo/options.py | 27 ---- .../optimization/tranquilo/region.py | 116 ++++++++++++++++++ .../optimization/tranquilo/sample_points.py | 98 ++++----------- .../optimization/tranquilo/tranquilo.py | 18 ++- .../optimization/tranquilo/weighting.py | 2 +- src/estimagic/visualization/deviation_plot.py | 4 +- .../visualization/visualize_tranquilo.py | 3 - .../tranquilo/test_acceptance_decision.py | 10 +- tests/optimization/tranquilo/test_bounds.py | 38 ++++++ .../tranquilo/test_estimate_variance.py | 2 +- .../tranquilo/test_filter_points.py | 6 +- .../optimization/tranquilo/test_fit_models.py | 10 +- tests/optimization/tranquilo/test_geometry.py | 24 ++-- tests/optimization/tranquilo/test_models.py | 10 +- tests/optimization/tranquilo/test_region.py | 90 ++++++++++++++ .../optimization/tranquilo/test_rho_noise.py | 4 +- .../tranquilo/test_sample_points.py | 43 +++---- 21 files changed, 369 insertions(+), 195 deletions(-) create mode 100644 src/estimagic/optimization/tranquilo/bounds.py create mode 100644 src/estimagic/optimization/tranquilo/region.py create mode 100644 tests/optimization/tranquilo/test_bounds.py create mode 100644 tests/optimization/tranquilo/test_region.py diff --git a/src/estimagic/optimization/tranquilo/bounds.py b/src/estimagic/optimization/tranquilo/bounds.py new file mode 100644 index 000000000..2b683e5d1 --- /dev/null +++ b/src/estimagic/optimization/tranquilo/bounds.py @@ -0,0 +1,27 @@ +from dataclasses import dataclass, replace + +import numpy as np + + +@dataclass +class Bounds: + """Parameter bounds.""" + + lower: np.ndarray + upper: np.ndarray + + def __post_init__(self): + self.has_any = _any_finite(self.lower, self.upper) + + # make it behave like a NamedTuple + def _replace(self, **kwargs): + return replace(self, **kwargs) + + +def _any_finite(lb, ub): + out = False + if lb is not None and np.isfinite(lb).any(): + out = True + if ub is not None and np.isfinite(ub).any(): + out = True + return out diff --git a/src/estimagic/optimization/tranquilo/estimate_variance.py b/src/estimagic/optimization/tranquilo/estimate_variance.py index bbc1faa42..4d5177610 100644 --- a/src/estimagic/optimization/tranquilo/estimate_variance.py +++ b/src/estimagic/optimization/tranquilo/estimate_variance.py @@ -5,7 +5,7 @@ from estimagic.optimization.tranquilo.get_component import get_component from estimagic.optimization.tranquilo.new_history import History -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region def get_variance_estimator(fitter, user_options): diff --git a/src/estimagic/optimization/tranquilo/geometry.py b/src/estimagic/optimization/tranquilo/geometry.py index 3b28ccb14..9edaf2062 100644 --- a/src/estimagic/optimization/tranquilo/geometry.py +++ b/src/estimagic/optimization/tranquilo/geometry.py @@ -2,12 +2,8 @@ import numpy as np -from estimagic.optimization.tranquilo.options import Region -from estimagic.optimization.tranquilo.sample_points import ( - _get_effective_bounds, - _map_from_feasible_trustregion, - get_sampler, -) +from estimagic.optimization.tranquilo.region import Region +from estimagic.optimization.tranquilo.sample_points import get_sampler def get_geometry_checker_pair( @@ -26,7 +22,7 @@ def get_geometry_checker_pair( samples drawn inside a box or a ball, respectively. n_params (int): Number of parameters. n_simulations (int): Number of simulations for the mean calculation. - bounds (Bounds): The parameter bounds. + bounds (Bounds): The parameter bounds. See module bounds.py. Returns: callable: The sample quality calculator. @@ -45,7 +41,7 @@ def get_geometry_checker_pair( _checker = built_in_checker[checker] - quality_calculator = partial(_checker["quality_calculator"], bounds=bounds) + quality_calculator = _checker["quality_calculator"] cutoff_simulator = partial( _checker["cutoff_simulator"], reference_sampler=reference_sampler, @@ -66,7 +62,7 @@ def log_d_cutoff_simulator( rng (np.random.Generator): The random number generator. reference_sampler (str): Either "box" or "ball", corresponding to comparison samples drawn inside a box or a ball, respectively. - bounds (Bounds): The parameter bounds. + bounds (Bounds): The parameter bounds. See module bounds.py. n_params (int): Dimensionality of the sample. n_simulations (int): Number of simulations for the mean calculation. @@ -74,18 +70,18 @@ def log_d_cutoff_simulator( float: The simulated mean logarithm of the d-optimality criterion. """ - _sampler = get_sampler(reference_sampler, bounds) - trustregion = Region(center=np.zeros(n_params), radius=1, shape=None) + _sampler = get_sampler(reference_sampler) + trustregion = Region(center=np.zeros(n_params), radius=1.0, bounds=bounds) sampler = partial(_sampler, trustregion=trustregion) raw = [] for _ in range(n_simulations): x = sampler(n_points=n_samples, rng=rng) - raw.append(log_d_quality_calculator(x, trustregion, bounds)) + raw.append(log_d_quality_calculator(x, trustregion)) out = np.nanmean(raw) return out -def log_d_quality_calculator(sample, trustregion, bounds): +def log_d_quality_calculator(sample, trustregion): """Logarithm of the d-optimality criterion. For a data sample x the log_d_criterion is defined as log(det(x.T @ x)). If the @@ -94,15 +90,13 @@ def log_d_quality_calculator(sample, trustregion, bounds): Args: sample (np.ndarray): The data sample, shape = (n, p). - trustregion (TrustRegion): NamedTuple with attributes center and radius. - bounds (Bounds): The parameter bounds. + trustregion (Region): Trustregion. See module region.py. Returns: np.ndarray: The criterion values, shape = (n, ). """ - effective_bounds = _get_effective_bounds(trustregion, bounds) - points = _map_from_feasible_trustregion(sample, effective_bounds) + points = trustregion.map_to_unit(sample) n_samples, n_params = points.shape xtx = points.T @ points det = np.linalg.det(xtx / n_samples) diff --git a/src/estimagic/optimization/tranquilo/models.py b/src/estimagic/optimization/tranquilo/models.py index c9a4d03d2..9764061c2 100644 --- a/src/estimagic/optimization/tranquilo/models.py +++ b/src/estimagic/optimization/tranquilo/models.py @@ -4,7 +4,7 @@ import numpy as np from numba import njit -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region @dataclass diff --git a/src/estimagic/optimization/tranquilo/options.py b/src/estimagic/optimization/tranquilo/options.py index 09150241c..4c3cac90d 100644 --- a/src/estimagic/optimization/tranquilo/options.py +++ b/src/estimagic/optimization/tranquilo/options.py @@ -1,4 +1,3 @@ -from dataclasses import dataclass from typing import NamedTuple import numpy as np @@ -9,17 +8,6 @@ ) -@dataclass -class Bounds: - """Stopping criteria.""" - - lower: np.ndarray - upper: np.ndarray - - def __post_init__(self): - self.has_any = _check_if_there_are_bounds(self.lower, self.upper) - - class StopOptions(NamedTuple): """Criteria for stopping without successful convergence.""" @@ -54,12 +42,6 @@ class RadiusOptions(NamedTuple): max_radius_to_step_ratio: float = np.inf -class Region(NamedTuple): - center: np.ndarray - radius: float - shape: str - - class AcceptanceOptions(NamedTuple): confidence_level: float = 0.8 power_level: float = 0.8 @@ -76,12 +58,3 @@ class StagnationOptions(NamedTuple): sample_increment: int = 1 max_trials: int = 1 drop: bool = True - - -def _check_if_there_are_bounds(lb, ub): - out = False - if lb is not None and np.isfinite(lb).any(): - out = True - if ub is not None and np.isfinite(ub).any(): - out = True - return out diff --git a/src/estimagic/optimization/tranquilo/region.py b/src/estimagic/optimization/tranquilo/region.py new file mode 100644 index 000000000..bdaeb5443 --- /dev/null +++ b/src/estimagic/optimization/tranquilo/region.py @@ -0,0 +1,116 @@ +from dataclasses import dataclass, replace + +import numpy as np + +from estimagic.optimization.tranquilo.bounds import Bounds +from estimagic.optimization.tranquilo.volume import ( + get_radius_of_cube_with_volume_of_sphere, +) + + +@dataclass +class Region: + """Trust region.""" + + center: np.ndarray + radius: float + bounds: Bounds = None + + @property + def shape(self) -> str: + any_bounds_binding = _any_bounds_binding( + bounds=self.bounds, center=self.center, radius=self.radius + ) + return "cube" if any_bounds_binding else "sphere" + + @property + def cube_bounds(self) -> Bounds: + if self.shape == "sphere": + raise AttributeError( + "The trustregion is a sphere, and thus has no cube bounds." + ) + radius = get_radius_of_cube_with_volume_of_sphere( + self.radius, len(self.center), scaling_factor=1.0 + ) + bounds = _get_cube_bounds(center=self.center, radius=radius, bounds=self.bounds) + return bounds + + @property + def cube_center(self) -> np.ndarray: + if self.shape == "sphere": + raise AttributeError( + "The trustregion is a sphere, and thus has no cube center." + ) + cube_bounds = self.cube_bounds + center = (cube_bounds.lower + cube_bounds.upper) / 2 + return center + + def map_to_unit(self, x: np.ndarray) -> np.ndarray: + """Map points from the trustregion to the unit sphere or cube.""" + if self.shape == "sphere": + out = _map_to_unit_sphere(x, center=self.center, radius=self.radius) + else: + out = _map_to_unit_cube(x, cube_bounds=self.cube_bounds) + return out + + def map_from_unit(self, x: np.ndarray) -> np.ndarray: + """Map points from the unit sphere or cube to the trustregion.""" + if self.shape == "sphere": + out = _map_from_unit_sphere(x, center=self.center, radius=self.radius) + else: + out = _map_from_unit_cube(x, cube_bounds=self.cube_bounds) + return out + + # make it behave like a NamedTuple + def _replace(self, **kwargs): + return replace(self, **kwargs) + + +def _map_to_unit_cube(x, cube_bounds): + """Map points from the trustregion to the unit cube.""" + out = 2 * (x - cube_bounds.lower) / (cube_bounds.upper - cube_bounds.lower) - 1 + return out + + +def _map_to_unit_sphere(x, center, radius): + """Map points from the trustregion to the unit sphere.""" + out = (x - center) / radius + return out + + +def _map_from_unit_cube(x, cube_bounds): + """Map points from the unit cube to the trustregion.""" + out = (cube_bounds.upper - cube_bounds.lower) * (x + 1) / 2 + cube_bounds.lower + return out + + +def _map_from_unit_sphere(x, center, radius): + """Map points from the unit sphere to the trustregion.""" + out = x * radius + center + return out + + +def _get_cube_bounds(center, radius, bounds): + """Get new bounds that define the intersection of the trustregion and the bounds.""" + lower_bounds = center - radius + upper_bounds = center + radius + + if bounds is not None and bounds.lower is not None: + lower_bounds = np.clip(lower_bounds, bounds.lower, np.inf) + + if bounds is not None and bounds.upper is not None: + upper_bounds = np.clip(upper_bounds, -np.inf, bounds.upper) + + return Bounds(lower=lower_bounds, upper=upper_bounds) + + +def _any_bounds_binding(bounds, center, radius): + """Check if any bound is binding, i.e. inside the trustregion.""" + out = False + if bounds is not None and bounds.has_any: + if bounds.lower is not None: + lower_binding = np.min(center - bounds.lower) <= radius + if bounds.upper is not None: + upper_binding = np.min(bounds.upper - center) <= radius + out = lower_binding or upper_binding + return out diff --git a/src/estimagic/optimization/tranquilo/sample_points.py b/src/estimagic/optimization/tranquilo/sample_points.py index 209d8e108..8f1e05930 100644 --- a/src/estimagic/optimization/tranquilo/sample_points.py +++ b/src/estimagic/optimization/tranquilo/sample_points.py @@ -6,10 +6,9 @@ import estimagic as em from estimagic.optimization.tranquilo.get_component import get_component -from estimagic.optimization.tranquilo.options import Bounds -def get_sampler(sampler, bounds, model_info=None, user_options=None): +def get_sampler(sampler, model_info=None, user_options=None): """Get sampling function partialled options. Args: @@ -19,7 +18,6 @@ def get_sampler(sampler, bounds, model_info=None, user_options=None): 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. If sampler is 'hull_sampler' or 'optimal_hull_sampler' the user options must contain the @@ -56,12 +54,10 @@ def get_sampler(sampler, bounds, model_info=None, user_options=None): raise ValueError(msg) default_options = { - "bounds": bounds, "model_info": model_info, } mandatory_args = [ - "bounds", "trustregion", "n_points", "existing_xs", @@ -85,7 +81,6 @@ def _box_sampler( n_points, rng, existing_xs=None, # noqa: ARG001 - bounds=None, ): """Naive random generation of trustregion points inside a box. @@ -100,28 +95,31 @@ def _box_sampler( outside of the sampler. Args: - trustregion (TrustRegion): NamedTuple with attributes center and radius. + trustregion (Region): Trustregion. See module region.py. n_points (int): how many new points to sample rng (numpy.random.Generator): Random number generator. 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_params = len(trustregion.center) - effective_bounds = _get_effective_bounds(trustregion, bounds=bounds) + + bounds = trustregion.cube_bounds points = rng.uniform( - low=effective_bounds.lower, - high=effective_bounds.upper, + low=bounds.lower, + high=bounds.upper, size=(n_points, n_params), ) return points def _ball_sampler( - trustregion, n_points, rng, existing_xs=None, bounds=None # noqa: ARG001 + trustregion, + n_points, + rng, + existing_xs=None, # noqa: ARG001 ): """Naive random generation of trustregion points inside a ball. @@ -131,24 +129,21 @@ def _ball_sampler( Code is adapted from https://tinyurl.com/y3p2dz6b. Args: - trustregion (TrustRegion): NamedTuple with attributes center and radius. + trustregion (Region): Trustregion. See module region.py. n_points (int): how many new points to sample rng (numpy.random.Generator): Random number generator. 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_params = len(trustregion.center) - effective_bounds = _get_effective_bounds(trustregion, bounds=bounds) raw = rng.normal(size=(n_points, n_params)) norm = np.linalg.norm(raw, axis=1, ord=2) scale = gammainc(n_params / 2, norm**2 / 2) ** (1 / n_params) / norm points = raw * scale.reshape(-1, 1) - - out = _map_into_feasible_trustregion(points, bounds=effective_bounds) + out = trustregion.map_from_unit(points) return out @@ -159,7 +154,6 @@ def _hull_sampler( order, distribution=None, existing_xs=None, # noqa: ARG001 - bounds=None, ): """Random generation of trustregion points on the hull of general sphere / cube. @@ -168,7 +162,7 @@ def _hull_sampler( defined by the intersection of the trustregion and the bounds. Args: - trustregion (TrustRegion): NamedTuple with attributes center and radius. + trustregion (Region): Trustregion. See module region.py. n_points (int): how many new points to sample rng (numpy.random.Generator): Random number generator. order (int): Type of norm to use when scaling the sampled points. For 2 it will @@ -178,18 +172,16 @@ def _hull_sampler( 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_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 + raw = _draw_from_distribution(distribution, rng=rng, size=(n_points, n_params)) + points = _project_onto_unit_hull(raw, order=order) + out = trustregion.map_from_unit(points) + return out def _optimal_hull_sampler( @@ -200,7 +192,6 @@ def _optimal_hull_sampler( distribution=None, hardness=1, existing_xs=None, - bounds=None, algorithm="scipy_lbfgsb", algo_options=None, criterion=None, @@ -218,7 +209,7 @@ def _optimal_hull_sampler( seek: https://tinyurl.com/mrythbk4. Args: - trustregion (TrustRegion): NamedTuple with attributes center and radius. + trustregion (Region): Trustregion. See module region.py. n_points (int): how many new points to sample rng (numpy.random.Generator): Random number generator. order (int): Type of norm to use when scaling the sampled points. For 2 it will @@ -231,7 +222,6 @@ def _optimal_hull_sampler( 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 @@ -273,11 +263,10 @@ def _optimal_hull_sampler( if "stopping_max_iterations" not in algo_options: algo_options["stopping_max_iterations"] = 2 * n_params + 5 - 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) + + existing_xs_unit = trustregion.map_to_unit(existing_xs) if criterion == "distance": dist_to_center = np.linalg.norm(existing_xs_unit, axis=1) @@ -348,7 +337,7 @@ def _optimal_hull_sampler( opt_params = x0 points = _project_onto_unit_hull(opt_params.reshape(-1, n_params), order=order) - points = _map_into_feasible_trustregion(points, bounds=effective_bounds) + points = trustregion.map_from_unit(points) # Collect additional information. Mostly used for testing. info = { @@ -473,38 +462,6 @@ def _draw_from_distribution(distribution, rng, size): 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 bounds 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. @@ -520,16 +477,3 @@ def _project_onto_unit_hull(x, order): 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 - - if bounds is not None and bounds.lower is not None: - lower_bounds = np.clip(lower_bounds, bounds.lower, np.inf) - - if bounds is not None and bounds.upper is not None: - upper_bounds = np.clip(upper_bounds, -np.inf, bounds.upper) - - return Bounds(lower=lower_bounds, upper=upper_bounds) diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 4dfacf7ca..8b755a6bd 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -9,6 +9,7 @@ from estimagic.optimization.tranquilo.acceptance_decision import get_acceptance_decider from estimagic.optimization.tranquilo.adjust_radius import adjust_radius from estimagic.optimization.tranquilo.aggregate_models import get_aggregator +from estimagic.optimization.tranquilo.bounds import Bounds from estimagic.optimization.tranquilo.estimate_variance import get_variance_estimator from estimagic.optimization.tranquilo.filter_points import ( drop_worst_points, @@ -23,12 +24,11 @@ from estimagic.optimization.tranquilo.new_history import History from estimagic.optimization.tranquilo.options import ( AcceptanceOptions, - Bounds, ConvOptions, RadiusOptions, - Region, StagnationOptions, ) +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.sample_points import get_sampler from estimagic.optimization.tranquilo.solve_subproblem import get_subsolver from estimagic.optimization.tranquilo.wrap_criterion import get_wrapped_criterion @@ -67,7 +67,6 @@ def _tranquilo( acceptance_options=None, variance_estimator="classic", variance_estimation_options=None, - trustregion_shape=None, stagnation_options=None, n_evals_per_point=1, ): @@ -158,14 +157,14 @@ def _tranquilo( bounds = Bounds(lower=lower_bounds, upper=upper_bounds) - if trustregion_shape is None: - trustregion_shape = "sphere" if not bounds.has_any else "cube" - if sampler is None: - sampler = f"optimal_{trustregion_shape}" + sampler = "optimal_cube" if bounds.has_any else "optimal_sphere" if subsolver is None: - subsolver = "bntr_fast" if trustregion_shape == "cube" else "gqtpar_fast" + if bounds.has_any: + subsolver = "bntr_fast" + else: + subsolver = "gqtpar_fast" if search_radius_factor is None: search_radius_factor = 4.25 if functype == "scalar" else 5.0 @@ -220,7 +219,6 @@ def _tranquilo( sample_points = get_sampler( sampler, - bounds=bounds, model_info=model_type, user_options=sampler_options, ) @@ -265,7 +263,7 @@ def _tranquilo( _init_fvec = history.get_fvecs(0).mean(axis=0) _init_radius = radius_options.initial_radius * np.max(np.abs(x)) - _init_region = Region(center=x, radius=_init_radius, shape=trustregion_shape) + _init_region = Region(center=x, radius=_init_radius, bounds=bounds) _init_vector_model = VectorModel( intercepts=_init_fvec, diff --git a/src/estimagic/optimization/tranquilo/weighting.py b/src/estimagic/optimization/tranquilo/weighting.py index 6fb7ecc4e..1f655ff83 100644 --- a/src/estimagic/optimization/tranquilo/weighting.py +++ b/src/estimagic/optimization/tranquilo/weighting.py @@ -6,7 +6,7 @@ def get_sample_weighter(weighter, bounds): The resulting function takes the following arguments: - xs (np.ndarray): A 2d numpy array containing a sample. - - trustregion (TrustRegion): The current trustregion. + - trustregion (Region): Trustregion. See module region.py. Args: weighter (str) diff --git a/src/estimagic/visualization/deviation_plot.py b/src/estimagic/visualization/deviation_plot.py index 62f5d78fe..8f8280803 100644 --- a/src/estimagic/visualization/deviation_plot.py +++ b/src/estimagic/visualization/deviation_plot.py @@ -70,7 +70,9 @@ def deviation_plot( .reset_index() ) average_deviations = ( - deviations.groupby(["algorithm", "n_evaluations"]).mean()[outcome].reset_index() + deviations.groupby(["algorithm", "n_evaluations"]) + .mean(numeric_only=True)[outcome] + .reset_index() ) fig = px.line(average_deviations, x="n_evaluations", y=outcome, color="algorithm") diff --git a/src/estimagic/visualization/visualize_tranquilo.py b/src/estimagic/visualization/visualize_tranquilo.py index 41a8dbd97..9b1a8f1c7 100644 --- a/src/estimagic/visualization/visualize_tranquilo.py +++ b/src/estimagic/visualization/visualize_tranquilo.py @@ -11,7 +11,6 @@ from estimagic.optimization.optimize_result import OptimizeResult from estimagic.optimization.tranquilo.clustering import cluster from estimagic.optimization.tranquilo.geometry import log_d_quality_calculator -from estimagic.optimization.tranquilo.options import Bounds from estimagic.optimization.tranquilo.volume import get_radius_after_volume_scaling @@ -424,13 +423,11 @@ def _plot_distances_from_center(history, state, fig, col, rows): def _get_fekete_criterion(res): states = res.algorithm_output["states"][1:] history = np.array(res.history["params"]) - bounds = Bounds(lower=-np.inf, upper=np.inf) out = [np.nan] + [ log_d_quality_calculator( sample=history[state.model_indices], trustregion=state.trustregion, - bounds=bounds, ) for state in states ] diff --git a/tests/optimization/tranquilo/test_acceptance_decision.py b/tests/optimization/tranquilo/test_acceptance_decision.py index ea2e0cf23..df68a8fc3 100644 --- a/tests/optimization/tranquilo/test_acceptance_decision.py +++ b/tests/optimization/tranquilo/test_acceptance_decision.py @@ -8,10 +8,8 @@ calculate_rho, ) from estimagic.optimization.tranquilo.new_history import History -from estimagic.optimization.tranquilo.options import ( - AcceptanceOptions, - Region, -) +from estimagic.optimization.tranquilo.options import AcceptanceOptions +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.solve_subproblem import SubproblemResult from numpy.testing import assert_array_equal @@ -42,7 +40,7 @@ def acceptance_options(): # ====================================================================================== -trustregion = Region(center=np.zeros(2), radius=2.0, shape="sphere") +trustregion = Region(center=np.zeros(2), radius=2.0) State = namedtuple("State", "x trustregion fval index") states = [ # we will parametrize over `states` State(np.arange(2.0), trustregion, 0.25, 0), # better than candidate @@ -92,7 +90,7 @@ def test_get_acceptance_result(): candidate_fval = 0 candidate_index = 0 rho = 1 - tr = Region(center=np.zeros(2), radius=2.0, shape="sphere") + tr = Region(center=np.zeros(2), radius=2.0) old_state = namedtuple("State", "x fval index trustregion")(np.arange(2), 1, 1, tr) ar_when_accepted = _get_acceptance_result( diff --git a/tests/optimization/tranquilo/test_bounds.py b/tests/optimization/tranquilo/test_bounds.py new file mode 100644 index 000000000..fec2dbbf6 --- /dev/null +++ b/tests/optimization/tranquilo/test_bounds.py @@ -0,0 +1,38 @@ +import numpy as np +import pytest +from estimagic.optimization.tranquilo.bounds import Bounds, _any_finite + +CASES = [ + (np.array([1, 2]), np.array([5, 6]), True), + (np.array([1, 2]), None, True), + (None, np.array([5, 6]), True), + (None, None, False), + (np.array([np.inf, np.inf]), np.array([np.inf, np.inf]), False), + (np.array([-np.inf, -np.inf]), np.array([np.inf, np.inf]), False), + (np.array([1, 2]), np.array([np.inf, np.inf]), True), +] + + +@pytest.mark.parametrize("lb, ub, exp", CASES) +def test_any_finite_true(lb, ub, exp): + out = _any_finite(lb, ub) + assert out is exp + + +def test_bounds_none(): + bounds = Bounds(lower=None, upper=None) + assert bounds.has_any is False + + +def test_bounds_inifinite(): + lb = np.array([np.inf, np.inf]) + ub = np.array([np.inf, np.inf]) + bounds = Bounds(lower=lb, upper=ub) + assert bounds.has_any is False + + +def test_bounds_finite(): + lb = np.array([1, 2]) + ub = np.array([5, 6]) + bounds = Bounds(lower=lb, upper=ub) + assert bounds.has_any is True diff --git a/tests/optimization/tranquilo/test_estimate_variance.py b/tests/optimization/tranquilo/test_estimate_variance.py index 34a27a963..067f5d5d5 100644 --- a/tests/optimization/tranquilo/test_estimate_variance.py +++ b/tests/optimization/tranquilo/test_estimate_variance.py @@ -29,7 +29,7 @@ def test_estimate_variance_classic(model_type): history.add_evals(idxs.repeat(repetitions), evals) got = _estimate_variance_classic( - trustregion=Region(center=np.array([0.0, 0.0]), radius=1.0, shape="sphere"), + trustregion=Region(center=np.array([0.0, 0.0]), radius=1.0), history=history, model_type=model_type, max_distance_factor=1.0, diff --git a/tests/optimization/tranquilo/test_filter_points.py b/tests/optimization/tranquilo/test_filter_points.py index ccb545487..1c32783c0 100644 --- a/tests/optimization/tranquilo/test_filter_points.py +++ b/tests/optimization/tranquilo/test_filter_points.py @@ -4,7 +4,7 @@ _scaled_square_features, drop_collinear_pounders, ) -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.tranquilo import State from estimagic.optimization.tranquilo.tranquilo_history import History from numpy.testing import assert_array_almost_equal as aaae @@ -43,7 +43,7 @@ def basic_case(): ) indices = np.arange(len(xs)) - trustregion = Region(radius=radius, center=x_accepted, shape="sphere") + trustregion = Region(center=x_accepted, radius=radius) state = State( trustregion=trustregion, @@ -199,7 +199,7 @@ def test_indices_in_trust_region(basic_case): x_accepted = np.array([0.16004745, 0.00572722, 0.01158929]) radius = 0.0125 - trustregion = Region(center=x_accepted, radius=radius, shape="sphere") + trustregion = Region(center=x_accepted, radius=radius) history.add_entries(xs, np.zeros(xs.shape[0])) indices_in_tr = history.get_indices_in_region(trustregion) diff --git a/tests/optimization/tranquilo/test_fit_models.py b/tests/optimization/tranquilo/test_fit_models.py index 9386a18a4..942004ef5 100644 --- a/tests/optimization/tranquilo/test_fit_models.py +++ b/tests/optimization/tranquilo/test_fit_models.py @@ -2,7 +2,7 @@ import pytest from estimagic import first_derivative, second_derivative from estimagic.optimization.tranquilo.fit_models import _quadratic_features, get_fitter -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region from numpy.testing import assert_array_almost_equal, assert_array_equal @@ -67,7 +67,7 @@ def test_fit_ols_against_truth(quadratic_case): got = fit_ols( x=quadratic_case["x"], y=quadratic_case["y"], - region=Region(center=np.zeros(4), radius=1.0, shape="sphere"), + region=Region(center=np.zeros(4), radius=1.0), old_model=None, ) @@ -80,7 +80,7 @@ def test_fit_powell_against_truth_quadratic(quadratic_case): got = fit_pounders( quadratic_case["x"], quadratic_case["y"], - region=Region(center=np.zeros(4), radius=1.0, shape="sphere"), + region=Region(center=np.zeros(4), radius=1.0), old_model=None, ) @@ -99,7 +99,7 @@ def test_fit_ols_against_gradient(model, quadratic_case): got = fit_ols( quadratic_case["x"], quadratic_case["y"], - region=Region(center=np.zeros(4), radius=1.0, shape="sphere"), + region=Region(center=np.zeros(4), radius=1.0), old_model=None, ) @@ -120,7 +120,7 @@ def test_fit_ols_against_hessian(model, options, quadratic_case): got = fit_ols( quadratic_case["x"], quadratic_case["y"], - region=Region(center=np.zeros(4), radius=1.0, shape="sphere"), + region=Region(center=np.zeros(4), radius=1.0), old_model=None, ) hessian = second_derivative(quadratic_case["func"], quadratic_case["x0"]) diff --git a/tests/optimization/tranquilo/test_geometry.py b/tests/optimization/tranquilo/test_geometry.py index e391d32ad..b3bf78f88 100644 --- a/tests/optimization/tranquilo/test_geometry.py +++ b/tests/optimization/tranquilo/test_geometry.py @@ -1,13 +1,13 @@ import numpy as np from estimagic.optimization.tranquilo.geometry import get_geometry_checker_pair -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.sample_points import get_sampler def test_geometry_checker(): - rng = np.random.default_rng() - sampler = get_sampler("sphere", bounds=None) - trustregion = Region(center=np.zeros(2), radius=1, shape="sphere") + rng = np.random.default_rng(12345) + sampler = get_sampler("sphere") + trustregion = Region(center=np.zeros(2), radius=1.0) x = sampler(trustregion, n_points=10, rng=rng) x_scaled = x * 0.5 @@ -16,8 +16,8 @@ def test_geometry_checker(): "d_optimality", reference_sampler="ball", n_params=2, bounds=None ) - x_quality = quality_calculator(x, trustregion, bounds=None) - x_scaled_quality = quality_calculator(x_scaled, trustregion, bounds=None) + x_quality = quality_calculator(x, trustregion) + x_scaled_quality = quality_calculator(x_scaled, trustregion) cutoff = cutoff_simulator(n_samples=10, rng=rng, n_simulations=1_000) @@ -26,11 +26,11 @@ def test_geometry_checker(): def test_geometry_checker_scale_invariance(): - rng = np.random.default_rng() - sampler = get_sampler("sphere", bounds=None) + rng = np.random.default_rng(12345) + sampler = get_sampler("sphere") - trustregion = Region(center=np.zeros(2), radius=1, shape="sphere") - trustregion_scaled = Region(center=np.ones(2), radius=2, shape="sphere") + trustregion = Region(center=np.zeros(2), radius=1.0) + trustregion_scaled = Region(center=np.ones(2), radius=2.0) x = sampler(trustregion, n_points=10, rng=rng) x_scaled = 1 + 2 * x @@ -39,7 +39,7 @@ def test_geometry_checker_scale_invariance(): "d_optimality", reference_sampler="ball", n_params=2, bounds=None ) - x_quality = quality_calculator(x, trustregion, bounds=None) - x_scaled_quality = quality_calculator(x_scaled, trustregion_scaled, bounds=None) + x_quality = quality_calculator(x, trustregion) + x_scaled_quality = quality_calculator(x_scaled, trustregion_scaled) assert x_quality == x_scaled_quality diff --git a/tests/optimization/tranquilo/test_models.py b/tests/optimization/tranquilo/test_models.py index daf255b85..12899c958 100644 --- a/tests/optimization/tranquilo/test_models.py +++ b/tests/optimization/tranquilo/test_models.py @@ -13,7 +13,7 @@ n_second_order_terms, subtract_models, ) -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region from numpy.testing import assert_array_almost_equal as aaae from numpy.testing import assert_array_equal @@ -112,7 +112,7 @@ def scalar_model(): intercept=0.5, linear_terms=np.array([-0.3, 0.3]), square_terms=np.array([[0.8, 0.2], [0.2, 0.7]]), - region=Region(center=np.array([0.2, 0.3]), radius=0.6, shape="sphere"), + region=Region(center=np.array([0.2, 0.3]), radius=0.6), ) return out @@ -129,14 +129,14 @@ def vector_model(): [[0.8, 0.2], [0.2, 0.7]], ] ), - region=Region(center=np.array([0.2, 0.3]), radius=0.6, shape="sphere"), + region=Region(center=np.array([0.2, 0.3]), radius=0.6), ) return out def test_move_scalar_model(scalar_model): old_region = scalar_model.region - new_region = Region(center=np.array([-0.1, 0.1]), radius=0.45, shape="sphere") + new_region = Region(center=np.array([-0.1, 0.1]), radius=0.45) old_model = scalar_model x_unscaled = np.array([[0.5, 0.5]]) @@ -156,7 +156,7 @@ def test_move_scalar_model(scalar_model): def test_move_vector_model(vector_model): old_region = vector_model.region - new_region = Region(center=np.array([-0.1, 0.1]), radius=0.45, shape="sphere") + new_region = Region(center=np.array([-0.1, 0.1]), radius=0.45) old_model = vector_model diff --git a/tests/optimization/tranquilo/test_region.py b/tests/optimization/tranquilo/test_region.py new file mode 100644 index 000000000..208ffcef5 --- /dev/null +++ b/tests/optimization/tranquilo/test_region.py @@ -0,0 +1,90 @@ +import numpy as np +import pytest +from estimagic.optimization.tranquilo.bounds import Bounds +from estimagic.optimization.tranquilo.region import ( + Region, + _any_bounds_binding, + _get_cube_bounds, + _map_from_unit_cube, + _map_from_unit_sphere, + _map_to_unit_cube, + _map_to_unit_sphere, +) +from numpy.testing import assert_array_equal + + +def test_map_to_unit_sphere(): + got = _map_to_unit_sphere(np.ones(2), center=2 * np.ones(1), radius=2) + assert_array_equal(got, -0.5 * np.ones(2)) + + +def test_map_to_unit_cube(): + bounds = Bounds(lower=np.zeros(2), upper=2 * np.ones(2)) + got = _map_to_unit_cube(np.ones(2), cube_bounds=bounds) + assert_array_equal(got, np.zeros(2)) + + +def test_map_from_unit_sphere(): + got = _map_from_unit_sphere(-0.5 * np.ones(2), center=2 * np.ones(1), radius=2) + assert_array_equal(got, np.ones(2)) + + +def test_map_from_unit_cube(): + bounds = Bounds(lower=np.zeros(2), upper=2 * np.ones(2)) + got = _map_from_unit_cube(np.zeros(2), cube_bounds=bounds) + assert_array_equal(got, np.ones(2)) + + +def test_any_bounds_binding_true(): + bounds = Bounds(lower=-np.ones(2), upper=np.ones(2)) + out = _any_bounds_binding(bounds, center=np.zeros(2), radius=2) + assert out + + +def test_any_bounds_binding_false(): + bounds = Bounds(lower=-np.ones(2), upper=np.ones(2)) + out = _any_bounds_binding(bounds, center=np.zeros(2), radius=0.5) + assert not out + + +def test_get_cube_bounds(): + bounds = Bounds(lower=-np.ones(2), upper=np.ones(2)) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds) + assert_array_equal(out.lower, bounds.lower) + assert_array_equal(out.upper, bounds.upper) + + +def test_get_cube_bounds_no_bounds(): + bounds = Bounds(lower=None, upper=None) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds) + assert_array_equal(out.lower, -np.ones(2)) + assert_array_equal(out.upper, np.ones(2)) + + +def test_get_cube_bounds_updated_upper_bounds(): + bounds = Bounds(lower=-2 * np.ones(2), upper=0.5 * np.ones(2)) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds) + assert_array_equal(out.lower, -np.ones(2)) + assert_array_equal(out.upper, 0.5 * np.ones(2)) + + +def test_region_non_binding_bounds(): + region = Region(center=np.zeros(2), radius=1) + assert region.shape == "sphere" + assert region.radius == 1 + assert region.bounds is None + with pytest.raises(AttributeError): + region.cube_bounds + with pytest.raises(AttributeError): + region.cube_center + + +def test_region_binding_bounds(): + bounds = Bounds(lower=-np.ones(2), upper=0.5 * np.ones(2)) + region = Region(center=np.zeros(2), radius=1, bounds=bounds) + assert region.shape == "cube" + assert region.radius == 1 + assert region.bounds is bounds + # shrinkage because cube radius is smaller than (spherical) radius + assert np.all(bounds.lower - region.cube_bounds.lower < 0) + assert_array_equal(region.cube_bounds.upper, bounds.upper) diff --git a/tests/optimization/tranquilo/test_rho_noise.py b/tests/optimization/tranquilo/test_rho_noise.py index 872a1ce8b..c59a92ae5 100644 --- a/tests/optimization/tranquilo/test_rho_noise.py +++ b/tests/optimization/tranquilo/test_rho_noise.py @@ -2,7 +2,7 @@ import pytest from estimagic.optimization.tranquilo.aggregate_models import get_aggregator from estimagic.optimization.tranquilo.fit_models import get_fitter -from estimagic.optimization.tranquilo.options import Region +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.rho_noise import simulate_rho_noise from estimagic.optimization.tranquilo.solve_subproblem import get_subsolver from numpy.testing import assert_array_almost_equal as aaae @@ -49,7 +49,7 @@ def test_convergence_to_one_if_noise_is_tiny(functype): noise_cov = np.eye(n_residuals) * 1e-12 - trustregion = Region(center=np.ones(2) * 0.5, radius=1, shape="sphere") + trustregion = Region(center=np.ones(2) * 0.5, radius=1.0) model_fitter = get_fitter(fitter="ols", model_type=model_type) vector_model = model_fitter( diff --git a/tests/optimization/tranquilo/test_sample_points.py b/tests/optimization/tranquilo/test_sample_points.py index 75b9432cb..a2093f416 100644 --- a/tests/optimization/tranquilo/test_sample_points.py +++ b/tests/optimization/tranquilo/test_sample_points.py @@ -1,6 +1,7 @@ import numpy as np import pytest -from estimagic.optimization.tranquilo.options import Bounds, Region +from estimagic.optimization.tranquilo.bounds import Bounds +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.sample_points import ( _draw_from_distribution, _minimal_pairwise_distance_on_hull, @@ -15,10 +16,11 @@ @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=Region(center=np.zeros(2), radius=1, shape="sphere"), + bounds = Bounds(lower=np.array([-2.0, -2.0]), upper=np.array([0.25, 0.5])) + _sampler = get_sampler(sampler) + trustregion = Region(center=np.array([0.0, 0]), radius=1.5, bounds=bounds) + sample = _sampler( + trustregion=trustregion, n_points=5, rng=np.random.default_rng(1234), ) @@ -30,10 +32,11 @@ def test_bounds_are_satisfied(sampler): @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}) + bounds = Bounds(lower=np.array([-2.0, -2]), upper=np.array([0.25, 0.5])) + sampler = get_sampler("hull_sampler", user_options={"order": order}) + trustregion = Region(center=np.array([0.0, 0]), radius=1.5, bounds=bounds) sample = sampler( - trustregion=Region(center=np.zeros(2), radius=1, shape="sphere"), + trustregion=trustregion, n_points=5, rng=np.random.default_rng(1234), ) @@ -46,12 +49,10 @@ def test_bounds_are_satisfied_general_hull_sampler(order): @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=sampler, - bounds=Bounds(lower=-np.ones(3), upper=np.ones(3)), - ) + sampler = get_sampler(sampler=sampler) + bounds = Bounds(lower=-np.ones(3), upper=np.ones(3)) calculated = sampler( - trustregion=Region(center=np.zeros(3), radius=1, shape="sphere"), + trustregion=Region(center=np.zeros(3), radius=1, bounds=bounds), n_points=0, existing_xs=np.empty((5, 3)), rng=np.random.default_rng(1234), @@ -65,11 +66,10 @@ 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, ) + bounds = Bounds(lower=-np.ones(3), upper=np.ones(3)) calculated = sampler( - trustregion=Region(center=np.zeros(3), radius=1, shape="sphere"), + trustregion=Region(center=np.zeros(3), radius=1, bounds=bounds), n_points=3, existing_xs=np.ones((2, 3)), # same point implies min distance of zero always rng=np.random.default_rng(1234), @@ -83,17 +83,15 @@ 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)), ) - + bounds = Bounds(lower=-np.ones(3), upper=np.ones(3)) distances = [] for sampler in [standard_sampler, optimal_sampler]: sample = sampler( - trustregion=Region(center=np.zeros(3), radius=1, shape="sphere"), + trustregion=Region(center=np.zeros(3), radius=1, bounds=bounds), n_points=5, rng=np.random.default_rng(1234), ) @@ -111,12 +109,11 @@ def test_randomsearch(sampler, n_points_randomsearch): _sampler = get_sampler( sampler="optimal_" + sampler, - bounds=bounds, ) # optimal sampling without randomsearch _, info = _sampler( - trustregion=Region(center=np.zeros(3), radius=1, shape=sampler), + trustregion=Region(center=np.zeros(3), radius=1, bounds=bounds), n_points=5, rng=np.random.default_rng(0), return_info=True, @@ -124,7 +121,7 @@ def test_randomsearch(sampler, n_points_randomsearch): # optimal sampling with randomsearch _, info_randomsearch = _sampler( - trustregion=Region(center=np.zeros(3), radius=1, shape=sampler), + trustregion=Region(center=np.zeros(3), radius=1, bounds=bounds), n_points=5, rng=np.random.default_rng(0), n_points_randomsearch=n_points_randomsearch, From bb52773b89f4f01b563b2a33a7cab7b8b41f2afe Mon Sep 17 00:00:00 2001 From: timmens Date: Wed, 15 Mar 2023 11:43:39 +0100 Subject: [PATCH 6/8] Change dtype in test_tranquilo.py to np.float64 --- tests/optimization/tranquilo/test_tranquilo.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index 17940ec8b..a0686a53a 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -64,7 +64,7 @@ def test_internal_tranquilo_scalar_sphere_defaults( ): res = tranquilo( criterion=lambda x: x @ x, - x=np.arange(4), + x=np.arange(4, dtype=np.float64), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -104,7 +104,7 @@ def test_internal_tranquilo_scalar_sphere_imprecise_defaults( ): res = tranquilo( criterion=lambda x: x @ x, - x=np.arange(4), + x=np.arange(4, dtype=np.float64), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -121,7 +121,7 @@ def test_internal_tranquilo_scalar_sphere_imprecise_defaults( def test_external_tranquilo_scalar_sphere_defaults(): res = minimize( criterion=lambda x: x @ x, - params=np.arange(4), + params=np.arange(4, dtype=np.float64), algorithm="tranquilo", ) @@ -160,7 +160,7 @@ def test_internal_tranquilo_ls_sphere_defaults( ): res = tranquilo_ls( criterion=lambda x: x, - x=np.arange(5), + x=np.arange(5, dtype=np.float64), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -177,7 +177,7 @@ def test_internal_tranquilo_ls_sphere_defaults( def test_external_tranquilo_ls_sphere_defaults(): res = minimize( criterion=lambda x: x, - params=np.arange(5), + params=np.arange(5, dtype=np.float64), algorithm="tranquilo_ls", ) From 9a2cc4a8f7e2803d191de4394b30a2da5e782955 Mon Sep 17 00:00:00 2001 From: timmens Date: Wed, 15 Mar 2023 12:04:27 +0100 Subject: [PATCH 7/8] Revert dtype and lower testing accuracy --- tests/optimization/tranquilo/test_tranquilo.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index a0686a53a..b58de75ca 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -64,7 +64,7 @@ def test_internal_tranquilo_scalar_sphere_defaults( ): res = tranquilo( criterion=lambda x: x @ x, - x=np.arange(4, dtype=np.float64), + x=np.arange(4), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -104,7 +104,7 @@ def test_internal_tranquilo_scalar_sphere_imprecise_defaults( ): res = tranquilo( criterion=lambda x: x @ x, - x=np.arange(4, dtype=np.float64), + x=np.arange(4), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -121,7 +121,7 @@ def test_internal_tranquilo_scalar_sphere_imprecise_defaults( def test_external_tranquilo_scalar_sphere_defaults(): res = minimize( criterion=lambda x: x @ x, - params=np.arange(4, dtype=np.float64), + params=np.arange(4), algorithm="tranquilo", ) @@ -160,7 +160,7 @@ def test_internal_tranquilo_ls_sphere_defaults( ): res = tranquilo_ls( criterion=lambda x: x, - x=np.arange(5, dtype=np.float64), + x=np.arange(5), sample_filter=sample_filter, fitter=fitter, surrogate_model=surrogate_model, @@ -177,7 +177,7 @@ def test_internal_tranquilo_ls_sphere_defaults( def test_external_tranquilo_ls_sphere_defaults(): res = minimize( criterion=lambda x: x, - params=np.arange(5, dtype=np.float64), + params=np.arange(5), algorithm="tranquilo_ls", ) @@ -289,7 +289,7 @@ def _f(x): algo_options={"noisy": True}, ) - aaae(res.params, np.zeros(5), decimal=4) + aaae(res.params, np.zeros(5), decimal=3) @pytest.mark.slow() From 25cb4eb838a56afe4584f72f27e4502abe473a6b Mon Sep 17 00:00:00 2001 From: timmens Date: Wed, 15 Mar 2023 13:20:37 +0100 Subject: [PATCH 8/8] Incorporate comments from review --- src/estimagic/optimization/tranquilo/bounds.py | 5 +++-- src/estimagic/optimization/tranquilo/region.py | 15 +++++++++------ src/estimagic/optimization/tranquilo/volume.py | 12 +++++++----- tests/optimization/tranquilo/test_region.py | 7 +++++++ tests/optimization/tranquilo/test_volume.py | 12 ++++++++++++ 5 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/bounds.py b/src/estimagic/optimization/tranquilo/bounds.py index 2b683e5d1..f73275421 100644 --- a/src/estimagic/optimization/tranquilo/bounds.py +++ b/src/estimagic/optimization/tranquilo/bounds.py @@ -3,7 +3,7 @@ import numpy as np -@dataclass +@dataclass(frozen=True) class Bounds: """Parameter bounds.""" @@ -11,7 +11,8 @@ class Bounds: upper: np.ndarray def __post_init__(self): - self.has_any = _any_finite(self.lower, self.upper) + # cannot use standard __setattr__ because it is frozen + super().__setattr__("has_any", _any_finite(self.lower, self.upper)) # make it behave like a NamedTuple def _replace(self, **kwargs): diff --git a/src/estimagic/optimization/tranquilo/region.py b/src/estimagic/optimization/tranquilo/region.py index bdaeb5443..410d44dd6 100644 --- a/src/estimagic/optimization/tranquilo/region.py +++ b/src/estimagic/optimization/tranquilo/region.py @@ -8,7 +8,7 @@ ) -@dataclass +@dataclass(frozen=True) class Region: """Trust region.""" @@ -29,9 +29,7 @@ def cube_bounds(self) -> Bounds: raise AttributeError( "The trustregion is a sphere, and thus has no cube bounds." ) - radius = get_radius_of_cube_with_volume_of_sphere( - self.radius, len(self.center), scaling_factor=1.0 - ) + radius = get_radius_of_cube_with_volume_of_sphere(self.radius, len(self.center)) bounds = _get_cube_bounds(center=self.center, radius=radius, bounds=self.bounds) return bounds @@ -41,8 +39,7 @@ def cube_center(self) -> np.ndarray: raise AttributeError( "The trustregion is a sphere, and thus has no cube center." ) - cube_bounds = self.cube_bounds - center = (cube_bounds.lower + cube_bounds.upper) / 2 + center = _get_cube_center(bounds=self.cube_bounds) return center def map_to_unit(self, x: np.ndarray) -> np.ndarray: @@ -104,6 +101,12 @@ def _get_cube_bounds(center, radius, bounds): return Bounds(lower=lower_bounds, upper=upper_bounds) +def _get_cube_center(bounds): + """Get center of region defined by bounds.""" + center = (bounds.lower + bounds.upper) / 2 + return center + + def _any_bounds_binding(bounds, center, radius): """Check if any bound is binding, i.e. inside the trustregion.""" out = False diff --git a/src/estimagic/optimization/tranquilo/volume.py b/src/estimagic/optimization/tranquilo/volume.py index 09a2f1b0c..1c092f84a 100644 --- a/src/estimagic/optimization/tranquilo/volume.py +++ b/src/estimagic/optimization/tranquilo/volume.py @@ -15,26 +15,28 @@ def get_radius_after_volume_scaling(radius, dim, scaling_factor): return out -def get_radius_of_sphere_with_volume_of_cube(cube_radius, dim, scaling_factor): +def get_radius_of_sphere_with_volume_of_cube(cube_radius, dim, scaling_factor=None): log_radius = ( loggamma(dim / 2 + 1) / dim - np.log(np.pi) / 2 + np.log(2) + np.log(cube_radius) - + np.log(scaling_factor) / dim ) + if scaling_factor is not None: + log_radius += np.log(scaling_factor) / dim out = np.exp(log_radius) return out -def get_radius_of_cube_with_volume_of_sphere(sphere_radius, dim, scaling_factor): +def get_radius_of_cube_with_volume_of_sphere(sphere_radius, dim, scaling_factor=None): log_radius = ( - np.log(scaling_factor) / dim - + np.log(np.pi) / 2 + np.log(np.pi) / 2 + np.log(sphere_radius) - np.log(2) - loggamma(dim / 2 + 1) / dim ) + if scaling_factor is not None: + log_radius += np.log(scaling_factor) / dim out = np.exp(log_radius) return out diff --git a/tests/optimization/tranquilo/test_region.py b/tests/optimization/tranquilo/test_region.py index 208ffcef5..496036475 100644 --- a/tests/optimization/tranquilo/test_region.py +++ b/tests/optimization/tranquilo/test_region.py @@ -5,6 +5,7 @@ Region, _any_bounds_binding, _get_cube_bounds, + _get_cube_center, _map_from_unit_cube, _map_from_unit_sphere, _map_to_unit_cube, @@ -68,6 +69,12 @@ def test_get_cube_bounds_updated_upper_bounds(): assert_array_equal(out.upper, 0.5 * np.ones(2)) +def test_get_cube_center(): + bounds = Bounds(lower=np.array([0, 0.5]), upper=np.array([1, 10])) + out = _get_cube_center(bounds) + assert_array_equal(out, np.array([0.5, 5.25])) + + def test_region_non_binding_bounds(): region = Region(center=np.zeros(2), radius=1) assert region.shape == "sphere" diff --git a/tests/optimization/tranquilo/test_volume.py b/tests/optimization/tranquilo/test_volume.py index b63bfafb6..e09c500fa 100644 --- a/tests/optimization/tranquilo/test_volume.py +++ b/tests/optimization/tranquilo/test_volume.py @@ -44,6 +44,18 @@ def test_get_radius_of_cube_with_volume_of_sphere(dim): assert np.allclose(got, expected) +def test_get_radius_of_sphere_with_volume_of_cube_no_scaling(): + v1 = get_radius_of_sphere_with_volume_of_cube(2.0, 2, None) + v2 = get_radius_of_sphere_with_volume_of_cube(2.0, 2, 1.0) + assert v1 == v2 + + +def test_get_radius_of_cube_with_volume_of_sphere_no_scaling(): + v1 = get_radius_of_cube_with_volume_of_sphere(2.0, 2, None) + v2 = get_radius_of_cube_with_volume_of_sphere(2.0, 2, 1.0) + assert v1 == v2 + + @pytest.mark.parametrize("dim", dims) def test_radius_after_volume_rescaling_scaling_factor_sphere(dim): radius = 0.6