diff --git a/src/estimagic/optimization/pounders.py b/src/estimagic/optimization/pounders.py index d1410c7f5..afb98cb39 100644 --- a/src/estimagic/optimization/pounders.py +++ b/src/estimagic/optimization/pounders.py @@ -6,7 +6,7 @@ from estimagic.batch_evaluators import process_batch_evaluator from estimagic.config import DEFAULT_N_CORES from estimagic.decorators import mark_minimizer -from estimagic.optimization.history import LeastSquaresHistory +from estimagic.optimization.pounders_history import LeastSquaresHistory from estimagic.optimization.pounders_auxiliary import ( add_accepted_point_to_residual_model, add_geomtery_points_to_make_main_model_fully_linear, diff --git a/src/estimagic/optimization/history.py b/src/estimagic/optimization/pounders_history.py similarity index 100% rename from src/estimagic/optimization/history.py rename to src/estimagic/optimization/pounders_history.py diff --git a/src/estimagic/optimization/tranquilo/aggregate_models.py b/src/estimagic/optimization/tranquilo/aggregate_models.py index e8b3e7fb2..3de7c6408 100644 --- a/src/estimagic/optimization/tranquilo/aggregate_models.py +++ b/src/estimagic/optimization/tranquilo/aggregate_models.py @@ -50,7 +50,8 @@ def _aggregate_models_template(vector_model, aggregator): intercept=intercept, linear_terms=linear_terms, square_terms=square_terms, - region=vector_model.region, + shift=vector_model.shift, + scale=vector_model.scale, ) return scalar_model diff --git a/src/estimagic/optimization/tranquilo/estimate_variance.py b/src/estimagic/optimization/tranquilo/estimate_variance.py index cadd47562..c4503323d 100644 --- a/src/estimagic/optimization/tranquilo/estimate_variance.py +++ b/src/estimagic/optimization/tranquilo/estimate_variance.py @@ -4,7 +4,7 @@ import numpy as np from estimagic.optimization.tranquilo.get_component import get_component -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.options import VarianceEstimatorOptions diff --git a/src/estimagic/optimization/tranquilo/filter_points.py b/src/estimagic/optimization/tranquilo/filter_points.py index b1054af15..32263e01a 100644 --- a/src/estimagic/optimization/tranquilo/filter_points.py +++ b/src/estimagic/optimization/tranquilo/filter_points.py @@ -1,11 +1,8 @@ import numpy as np import scipy -from numba import njit -from scipy.linalg import qr_multiply from estimagic.optimization.tranquilo.clustering import cluster from estimagic.optimization.tranquilo.get_component import get_component -from estimagic.optimization.tranquilo.models import n_second_order_terms from estimagic.optimization.tranquilo.volume import get_radius_after_volume_scaling from estimagic.optimization.tranquilo.options import FilterOptions @@ -31,8 +28,6 @@ def get_sample_filter(sample_filter="keep_all", user_options=None): "discard_all": discard_all, "keep_all": keep_all, "clustering": keep_cluster_centers, - "keep_sphere": keep_sphere, - "drop_pounders": drop_collinear_pounders, "drop_excess": drop_excess, } @@ -55,22 +50,6 @@ def keep_all(xs, indices): return xs, indices -def keep_sphere(xs, indices, state): - dists = np.linalg.norm(xs - state.trustregion.center, axis=1) - keep = dists <= state.trustregion.radius - return xs[keep], indices[keep] - - -def drop_collinear_pounders(xs, indices, state): - """Drop collinear points using pounders filtering.""" - if xs.shape[0] <= xs.shape[1] + 1: - filtered_xs, filtered_indices = xs, indices - else: - filtered_xs, filtered_indices = _drop_collinear_pounders(xs, indices, state) - - return filtered_xs, filtered_indices - - def drop_excess(xs, indices, state, target_size): n_to_drop = max(0, len(xs) - target_size) @@ -100,8 +79,7 @@ def drop_worst_points(xs, indices, state, n_to_drop): n_dropped = 0 if n_dropped < n_to_drop: - order = 2 if state.trustregion.shape == "sphere" else np.inf - dists = np.linalg.norm(xs - state.trustregion.center, axis=1, ord=order) + dists = np.linalg.norm(xs - state.x, axis=1) while n_dropped < n_to_drop and (dists > state.trustregion.radius).any(): drop_index = np.argmax(dists) @@ -149,123 +127,3 @@ def keep_cluster_centers( # do I need to make sure trustregion center is in there? out = xs[centers], indices[centers] return out - - -def _drop_collinear_pounders(xs, indices, state): - theta2 = 1e-4 - n_samples, n_params = xs.shape - n_poly_terms = n_second_order_terms(n_params) - - indices_reverse = indices[::-1] - indexer_reverse = np.arange(n_samples)[::-1] - - radius = state.trustregion.radius - center = state.trustregion.center - index_center = int(np.where(indices_reverse == state.index)[0]) - centered_xs = (xs - center) / radius - - ( - linear_features, - square_features, - idx_list_n_plus_1, - index, - ) = _get_polynomial_feature_matrices( - centered_xs, - indexer_reverse, - index_center, - n_params, - n_samples, - n_poly_terms, - ) - - indexer_filtered = indexer_reverse[idx_list_n_plus_1].tolist() - _index_center = indexer_reverse[index_center] - - counter = n_params + 1 - - while (counter < n_samples) and (index >= 0): - if index == _index_center: - index -= 1 - continue - - linear_features[counter, 1:] = centered_xs[index] - square_features[counter, :] = _scaled_square_features( - linear_features[counter, 1:] - ) - - linear_features_pad = np.zeros((n_samples, n_samples)) - linear_features_pad[:n_samples, : n_params + 1] = linear_features - - n_z_mat, _ = qr_multiply( - linear_features_pad[: counter + 1, :], - square_features.T[:n_poly_terms, : counter + 1], - ) - beta = np.linalg.svd(n_z_mat.T[n_params + 1 :], compute_uv=False) - - if beta[min(counter - n_params, n_poly_terms) - 1] > theta2: - indexer_filtered += [index] - counter += 1 - - index -= 1 - - filtered_indices = indices[indexer_filtered] - filtered_xs = xs[indexer_filtered] - - return filtered_xs, filtered_indices - - -def _get_polynomial_feature_matrices( - centered_xs, indexer, index_center, n_params, n_samples, n_poly_terms -): - linear_features = np.zeros((n_samples, n_params + 1)) - square_features = np.zeros((n_samples, n_poly_terms)) - - linear_features[0, 1:] = centered_xs[indexer[index_center]] - square_features[0, :] = _scaled_square_features(linear_features[0, 1:]).flatten() - - _is_center_in_head = index_center < n_params - idx_list_n = [i for i in range(n_params + _is_center_in_head) if i != index_center] - idx_list_n_plus_1 = [index_center, *idx_list_n] - - linear_features[:, 0] = 1 - linear_features[: n_params + 1, 1:] = centered_xs[indexer[idx_list_n_plus_1]] - square_features[: n_params + 1, :] = _scaled_square_features( - linear_features[: n_params + 1, 1:] - ) - - idx = n_samples - _is_center_in_head - len(idx_list_n) - 1 - - return linear_features, square_features, idx_list_n_plus_1, idx - - -@njit -def _scaled_square_features(x): - """Construct scaled interaction and square terms. - - The interaction terms are scaled by 1 / sqrt{2} while the square terms are scaled - by 1 / 2. - - Args: - x (np.ndarray): Array of shape (n_samples, n_params). - - Returns: - np.ndarray: Scaled interaction and square terms. Has shape (n_samples, - n_params + (n_params - 1) * n_params / 1). - - """ - n_samples, n_params = np.atleast_2d(x).shape - 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): - poly_terms[idx] = xt[i] ** 2 / 2 - idx += 1 - - for j in range(i + 1, n_params): - poly_terms[idx] = xt[i] * xt[j] / np.sqrt(2) - idx += 1 - - return poly_terms.T diff --git a/src/estimagic/optimization/tranquilo/fit_models.py b/src/estimagic/optimization/tranquilo/fit_models.py index 9a419efc3..32f0b74ce 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -16,16 +16,23 @@ def get_fitter( - fitter, fitter_options=None, model_type=None, infinity_handling="relative" + fitter, + fitter_options=None, + model_type=None, + residualize=None, + infinity_handling=None, ): """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_type``. - - user_options (dict): Options for the fit method. The following are supported: + fitter (str or callable): Name of a fit method or a fit method. Arguments need + to be, in order, + - x (np.ndarray): Data points. + - y (np.ndarray): Corresponding function evaluations at data points. + - weighs (np.ndarray): Weights for the data points. + - model_type (str): Type of model to be fitted. + + fitter_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. @@ -34,12 +41,16 @@ def get_fitter( - "linear": Only linear effects and intercept. - "quadratic": Fully quadratic model. + residualize (bool): If True, the model is fitted to the residuals of the old + model. This introduces momentum when the coefficients are penalized. + + infinity_handling (str): How to handle infinite values in the data. Currently + supported: {"relative"}. See `handle_infinty.py`. + Returns: callable: The partialled fit method that only depends on x and y. """ - fitter_options = {} if fitter_options is None else fitter_options - built_in_fitters = { "ols": fit_ols, "ridge": fit_ridge, @@ -47,15 +58,13 @@ def get_fitter( "tranquilo": fit_tranquilo, } - default_options = FitterOptions(model_type=model_type) - mandatory_arguments = ["x", "y", "model_type"] _raw_fitter = get_component( name_or_func=fitter, component_name="fitter", func_dict=built_in_fitters, - default_options=default_options, + default_options=FitterOptions(), user_options=fitter_options, mandatory_signature=mandatory_arguments, ) @@ -67,7 +76,7 @@ def get_fitter( fitter=_raw_fitter, model_type=model_type, clip_infinite_values=clip_infinite_values, - residualize=fitter_options.get("residualize", False), + residualize=residualize, ) return fitter @@ -106,15 +115,13 @@ def _fitter_template( n_residuals = y.shape[1] y_clipped = clip_infinite_values(y) - x_centered = (x - region.center) / region.radius + x_unit = region.map_to_unit(x) if residualize: old_model_moved = move_model(old_model, region) - y_clipped = y_clipped - old_model_moved.predict(x_centered).reshape( - y_clipped.shape - ) + y_clipped = y_clipped - old_model_moved.predict(x_unit).reshape(y_clipped.shape) - coef = fitter(x=x_centered, y=y_clipped, weights=weights) + coef = fitter(x=x_unit, y=y_clipped, weights=weights, model_type=model_type) # results processing intercepts, linear_terms, square_terms = np.split(coef, (1, n_params + 1), axis=1) @@ -128,7 +135,13 @@ def _fitter_template( else: square_terms = None - results = VectorModel(intercepts, linear_terms, square_terms, region=region) + results = VectorModel( + intercepts, + linear_terms, + square_terms, + shift=region.effective_center, + scale=region.effective_radius, + ) if residualize: results = add_models(results, old_model_moved) diff --git a/src/estimagic/optimization/tranquilo/new_history.py b/src/estimagic/optimization/tranquilo/history.py similarity index 83% rename from src/estimagic/optimization/tranquilo/new_history.py rename to src/estimagic/optimization/tranquilo/history.py index 20c665041..88cd1018d 100644 --- a/src/estimagic/optimization/tranquilo/new_history.py +++ b/src/estimagic/optimization/tranquilo/history.py @@ -1,5 +1,4 @@ import numpy as np -from numba import njit class History: @@ -201,22 +200,9 @@ def get_x_indices_in_region(self, region): # early return if there are no entries if self.get_n_fun() == 0: return np.array([]) - - center = region.center - shape = region.shape - radius = region.radius - - if isinstance(center, int): - center = self.get_xs(x_indices=center) - xs = self.get_xs() - - out = _find_indices_in_trust_region(xs, center=center, radius=radius) - - if shape == "sphere": - mask = np.linalg.norm(xs[out] - region.center, axis=1) <= radius - out = out[mask] - + mask = np.linalg.norm(xs - region.center, axis=1) <= region.radius + out = np.arange(len(mask))[mask] return out def __repr__(self): @@ -243,40 +229,6 @@ def _add_entries_to_array(arr, new, position): return arr -@njit -def _find_indices_in_trust_region(xs, center, radius): - """Get the row indices of all parameter vectors in a trust region. - - This is for square trust regions, i.e. balls in term of an infinity norm. - - Args: - xs (np.ndarray): 2d numpy array where each row is a parameter vector. - center (np.ndarray): 1d numpy array that marks the center of the trust region. - radius (float): Radius of the trust region. - - Returns: - np.ndarray: The indices of parameters in the trust region. - - """ - n_obs, dim = xs.shape - out = np.zeros(n_obs).astype(np.int64) - success_counter = 0 - upper = center + radius - lower = center - radius - for i in range(n_obs): - success = True - for j in range(dim): - value = xs[i, j] - if not (lower[j] <= value <= upper[j]) or np.isnan(value): - success = False - continue - if success: - out[success_counter] = i - success_counter += 1 - - return out[:success_counter] - - def _extract_from_indices(arr, mapper, x_indices, n_xs): """Retrieve fvecs or fvals from history. diff --git a/src/estimagic/optimization/tranquilo/models.py b/src/estimagic/optimization/tranquilo/models.py index 9764061c2..bc2a4afcd 100644 --- a/src/estimagic/optimization/tranquilo/models.py +++ b/src/estimagic/optimization/tranquilo/models.py @@ -4,10 +4,8 @@ import numpy as np from numba import njit -from estimagic.optimization.tranquilo.region import Region - -@dataclass +@dataclass(frozen=True) class VectorModel: intercepts: np.ndarray # shape (n_residuals,) linear_terms: np.ndarray # shape (n_residuals, n_params) @@ -15,7 +13,10 @@ class VectorModel: np.ndarray, None ] = None # shape (n_residuals, n_params, n_params) - region: Region = None + # scale and shift correspond to effective_radius and effective_center of the region + # on which the model was fitted + scale: Union[float, np.ndarray] = None + shift: np.ndarray = None def predict(self, x: np.ndarray) -> np.ndarray: return _predict_vector(self, x) @@ -25,13 +26,16 @@ def _replace(self, **kwargs): return replace(self, **kwargs) -@dataclass +@dataclass(frozen=True) class ScalarModel: intercept: float linear_terms: np.ndarray # shape (n_params,) square_terms: Union[np.ndarray, None] = None # shape (n_params, n_params) - region: Region = None + # scale and shift correspond to effective_radius and effective_center of the region + # on which the model was fitted + scale: Union[float, np.ndarray] = None + shift: np.ndarray = None def predict(self, x: np.ndarray) -> np.ndarray: return _predict_scalar(self, x) @@ -41,8 +45,8 @@ def _replace(self, **kwargs): return replace(self, **kwargs) -def _predict_vector(model: VectorModel, centered_x: np.ndarray) -> np.ndarray: - """Evaluate a VectorModel at centered_x. +def _predict_vector(model: VectorModel, x_unit: np.ndarray) -> np.ndarray: + """Evaluate a VectorModel at x_unit. We utilize that a quadratic model can be written in the form: @@ -55,20 +59,20 @@ def _predict_vector(model: VectorModel, centered_x: np.ndarray) -> np.ndarray: f seperately. Args: - scalar_model (VectorModel): The aggregated model. Has entries: + model (VectorModel): The aggregated model. Has entries: - 'intercepts': corresponds to 'a' in the above equation - 'linear_terms': corresponds to 'g' in the above equation - 'square_terms': corresponds to 'H' in the above equation - x (np.ndarray): New data. Has shape (n_params,) or (n_samples, n_params). + x_unit (np.ndarray): New data. Has shape (n_params,) or (n_samples, n_params). Returns: np.ndarray: Model evaluations, has shape (n_samples, n_residuals) if x is 2d and (n_residuals,) if x is 1d. """ - is_flat_x = centered_x.ndim == 1 + is_flat_x = x_unit.ndim == 1 - x = np.atleast_2d(centered_x) + x = np.atleast_2d(x_unit) y = model.linear_terms @ x.T + model.intercepts.reshape(-1, 1) @@ -78,7 +82,7 @@ def _predict_vector(model: VectorModel, centered_x: np.ndarray) -> np.ndarray: if is_flat_x: out = y.flatten() else: - out = y.T.reshape(len(centered_x), -1) + out = y.T.reshape(len(x_unit), -1) return out @@ -97,11 +101,11 @@ def add_models(model1, model2): if type(model1) != type(model2): raise TypeError("Models must be of the same type.") - if not np.allclose(model1.region.center, model2.region.center): - raise ValueError("Models must have the same center.") + if not np.allclose(model1.shift, model2.shift): + raise ValueError("Models must have the same shift.") - if not np.allclose(model1.region.radius, model2.region.radius): - raise ValueError("Models must have the same radius.") + if not np.allclose(model1.scale, model2.scale): + raise ValueError("Models must have the same scale.") new = {} if isinstance(model1, ScalarModel): @@ -119,42 +123,6 @@ def add_models(model1, model2): return out -def subtract_models(model1, model2): - """Subtract two models. - - Args: - model1 (Union[ScalarModel, VectorModel]): The first model. - model2 (Union[ScalarModel, VectorModel]): The second model. - - Returns: - Union[ScalarModel, VectorModel]: The difference of the two models. - - """ - if type(model1) != type(model2): - raise TypeError("Models must be of the same type.") - - if not np.allclose(model1.region.center, model2.region.center): - raise ValueError("Models must have the same center.") - - if not np.allclose(model1.region.radius, model2.region.radius): - raise ValueError("Models must have the same radius.") - - new = {} - if isinstance(model1, ScalarModel): - new["intercept"] = model1.intercept - model2.intercept - else: - new["intercepts"] = model1.intercepts - model2.intercepts - - new["linear_terms"] = model1.linear_terms - model2.linear_terms - - if model1.square_terms is not None: - assert model2.square_terms is not None - new["square_terms"] = model1.square_terms - model2.square_terms - - out = replace(model1, **new) - return out - - def move_model(model, new_region): """Move a model to a new region. @@ -166,86 +134,78 @@ def move_model(model, new_region): Union[ScalarModel, VectorModel]: The moved model. """ - old_region = model.region - out = _scale_model(model, old_radius=old_region.radius, new_radius=1.0) + + # undo old scaling + out = _scale_model(model, factor=1 / model.scale) + + # shift the center + shift = new_region.effective_center - model.shift if isinstance(model, ScalarModel): - out = _shift_scalar_model( - out, old_center=old_region.center, new_center=new_region.center - ) + out = _shift_scalar_model(out, shift=shift) else: - out = _shift_vector_model( - out, old_center=old_region.center, new_center=new_region.center - ) - out = _scale_model(out, old_radius=1.0, new_radius=new_region.radius) + out = _shift_vector_model(out, shift=shift) + + # apply new scaling + new_scale = new_region.effective_radius + out = _scale_model(out, factor=new_scale) return out -def _scale_model(model, old_radius, new_radius): +def _scale_model(model, factor): """Scale a scalar or vector model to a new radius. Args: model (Union[ScalarModel, VectorModel]): The model to scale. - old_radius (float): The old radius. - new_radius (float): The new radius. + factor (Union[float, np.ndarray]): The scaling factor. Returns: Union[ScalarModel, VectorModel]: The scaled model. """ - factor = new_radius / old_radius - new_g = model.linear_terms * factor new_h = None if model.square_terms is None else model.square_terms * factor**2 out = model._replace( linear_terms=new_g, square_terms=new_h, - region=model.region._replace(radius=new_radius), + scale=model.scale * factor, ) return out -def _shift_scalar_model(model, old_center, new_center): +def _shift_scalar_model(model, shift): """Shift a scalar model to a new center. Args: model (ScalarModel): The model to shift. - old_center (np.ndarray): The old center. - new_center (np.ndarray): The new center. + shift (np.ndarray): The shift. Returns: ScalarModel: The shifted model. """ - shift = new_center - old_center - new_c = model.predict(shift) - - new_g = model.linear_terms + shift @ model.square_terms + new_g = model.linear_terms + model.square_terms @ shift out = model._replace( intercept=new_c, linear_terms=new_g, - region=model.region._replace(center=new_center), + shift=model.shift + shift, ) - return out -def _shift_vector_model(model, old_center, new_center): +def _shift_vector_model(model, shift): """Shift a vector model to a new center. Args: model (VectorModel): The model to shift. - old_center (np.ndarray): The old center. - new_center (np.ndarray): The new center. + shift (np.ndarray): The shift. Returns: VectorModel: The shifted model. """ - shift = new_center - old_center - new_c = model.predict(shift) new_g = model.linear_terms @@ -256,14 +216,13 @@ def _shift_vector_model(model, old_center, new_center): out = model._replace( intercepts=new_c, linear_terms=new_g, - region=model.region._replace(center=new_center), + shift=model.shift + shift, ) - return out -def _predict_scalar(model: ScalarModel, centered_x: np.ndarray) -> np.ndarray: - """Evaluate a ScalarModel at centered_x. +def _predict_scalar(model: ScalarModel, x_unit: np.ndarray) -> np.ndarray: + """Evaluate a ScalarModel at x_unit. We utilize that a quadratic model can be written in the form: @@ -274,11 +233,11 @@ def _predict_scalar(model: ScalarModel, centered_x: np.ndarray) -> np.ndarray: thought of as the gradient and Hessian. Args: - scalar_model (ScalarModel): The aggregated model. Has entries: + model (ScalarModel): The aggregated model. Has entries: - 'intercept': corresponds to 'a' in the above equation - 'linear_terms': corresponds to 'g' in the above equation - 'square_terms': corresponds to 'H' in the above equation - centered_x (np.ndarray): New data. Has shape (n_params,) or (n_samples, + x_unit (np.ndarray): New data. Has shape (n_params,) or (n_samples, n_params). Returns: @@ -286,9 +245,9 @@ def _predict_scalar(model: ScalarModel, centered_x: np.ndarray) -> np.ndarray: is 2d and a float otherwise. """ - is_flat_x = centered_x.ndim == 1 + is_flat_x = x_unit.ndim == 1 - x = np.atleast_2d(centered_x) + x = np.atleast_2d(x_unit) y = x @ model.linear_terms + model.intercept diff --git a/src/estimagic/optimization/tranquilo/options.py b/src/estimagic/optimization/tranquilo/options.py index 72f39957d..651f19100 100644 --- a/src/estimagic/optimization/tranquilo/options.py +++ b/src/estimagic/optimization/tranquilo/options.py @@ -27,13 +27,17 @@ def get_default_sample_size(model_type, x): def get_default_model_fitter(model_type, sample_size, x): n_params = n_free_params(dim=len(x), model_type=model_type) - if model_type == "linear": - fitter = "ridge" if sample_size < n_params else "ols" + if model_type == "linear" or sample_size >= n_params: + fitter = "ols" else: fitter = "tranquilo" return fitter +def get_default_residualize(model_fitter): + return model_fitter == "tranquilo" + + def get_default_subsolver(bounds, cube_subsolver, sphere_subsolver): return cube_subsolver if bounds.has_any else sphere_subsolver @@ -137,7 +141,6 @@ class SubsolverOptions(NamedTuple): class FitterOptions(NamedTuple): - model_type: str l2_penalty_linear: float = 0.0 l2_penalty_square: float = 0.1 p_intercept: float = 0.05 diff --git a/src/estimagic/optimization/tranquilo/process_arguments.py b/src/estimagic/optimization/tranquilo/process_arguments.py index 01a102be2..3836a292d 100644 --- a/src/estimagic/optimization/tranquilo/process_arguments.py +++ b/src/estimagic/optimization/tranquilo/process_arguments.py @@ -10,7 +10,7 @@ from estimagic.optimization.tranquilo.estimate_variance import get_variance_estimator from estimagic.optimization.tranquilo.filter_points import get_sample_filter from estimagic.optimization.tranquilo.fit_models import get_fitter -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History from estimagic.optimization.tranquilo.options import ( ConvOptions, StagnationOptions, @@ -19,6 +19,7 @@ get_default_aggregator, get_default_batch_size, get_default_model_fitter, + get_default_residualize, get_default_model_type, get_default_n_evals_at_start, get_default_radius_options, @@ -85,6 +86,7 @@ def process_arguments( variance_estimator="classic", variance_estimator_options=None, infinity_handler="relative", + residualize=None, ): # process convergence options conv_options = ConvOptions( @@ -133,6 +135,7 @@ def process_arguments( model_fitter = _process_model_fitter( model_fitter, model_type=model_type, sample_size=target_sample_size, x=x ) + residualize = _process_residualize(residualize, model_fitter=model_fitter) # initialize components history = History(functype=functype) @@ -172,6 +175,7 @@ def process_arguments( fitter_options=model_fitter_options, model_type=model_type, infinity_handling=infinity_handler, + residualize=residualize, ) aggregate_model = get_aggregator( @@ -292,6 +296,17 @@ def _process_model_fitter(model_fitter, model_type, sample_size, x): return out +def _process_residualize(residualize, model_fitter): + if residualize is None: + out = get_default_residualize(model_fitter) + else: + if not isinstance(residualize, bool): + raise ValueError("residualize must be a boolean.") + out = residualize + + return out + + def _process_n_evals_at_start(n_evals, noisy): if n_evals is None: out = get_default_n_evals_at_start(noisy) diff --git a/src/estimagic/optimization/tranquilo/region.py b/src/estimagic/optimization/tranquilo/region.py index e77a9a961..48abe2bcb 100644 --- a/src/estimagic/optimization/tranquilo/region.py +++ b/src/estimagic/optimization/tranquilo/region.py @@ -16,12 +16,19 @@ class Region: 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" + def __post_init__(self): + shape = _get_shape(self.center, self.radius, self.bounds) + cube_bounds = _get_cube_bounds(self.center, self.radius, self.bounds, shape) + cube_center = _get_cube_center(cube_bounds) + effective_center = _get_effective_center(shape, self.center, cube_center) + effective_radius = _get_effective_radius(shape, self.radius, cube_bounds) + + # cannot use standard __setattr__ because it is frozen + super().__setattr__("shape", shape) + super().__setattr__("_cube_bounds", cube_bounds) + super().__setattr__("_cube_center", cube_center) + super().__setattr__("effective_center", effective_center) + super().__setattr__("effective_radius", effective_radius) @property def cube_bounds(self) -> Bounds: @@ -29,9 +36,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)) - bounds = _get_cube_bounds(center=self.center, radius=radius, bounds=self.bounds) - return bounds + return self._cube_bounds @property def cube_center(self) -> np.ndarray: @@ -39,8 +44,7 @@ def cube_center(self) -> np.ndarray: raise AttributeError( "The trustregion is a sphere, and thus has no cube center." ) - center = _get_cube_center(bounds=self.cube_bounds) - return center + return self._cube_center def map_to_unit(self, x: np.ndarray) -> np.ndarray: """Map points from the trustregion to the unit sphere or cube.""" @@ -90,7 +94,39 @@ def _map_from_unit_sphere(x, center, radius): return out -def _get_cube_bounds(center, radius, bounds): +def _get_shape(center, radius, bounds): + any_bounds_binding = _any_bounds_binding( + bounds=bounds, center=center, radius=radius + ) + return "cube" if any_bounds_binding else "sphere" + + +def _get_cube_bounds(center, radius, bounds, shape): + if shape == "cube": + radius = get_radius_of_cube_with_volume_of_sphere(radius, len(center)) + cube_bounds = _create_cube_bounds(center=center, radius=radius, bounds=bounds) + return cube_bounds + + +def _get_cube_center(cube_bounds): + cube_center = (cube_bounds.lower + cube_bounds.upper) / 2 + return cube_center + + +def _get_effective_center(shape, center, cube_center): + effective_center = center if shape == "sphere" else cube_center + return effective_center + + +def _get_effective_radius(shape, radius, cube_bounds): + if shape == "sphere": + effective_radius = radius + else: + effective_radius = (cube_bounds.upper - cube_bounds.lower) / 2 + return effective_radius + + +def _create_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 @@ -104,12 +140,6 @@ 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 @@ -118,5 +148,5 @@ def _any_bounds_binding(bounds, center, radius): 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 + out = np.any(lower_binding) or np.any(upper_binding) return out diff --git a/src/estimagic/optimization/tranquilo/rho_noise.py b/src/estimagic/optimization/tranquilo/rho_noise.py index 8bfc8d387..9df77335e 100644 --- a/src/estimagic/optimization/tranquilo/rho_noise.py +++ b/src/estimagic/optimization/tranquilo/rho_noise.py @@ -46,9 +46,9 @@ def simulate_rho_noise( n_samples, n_params = xs.shape n_residuals = len(noise_cov) - centered_xs = (xs - trustregion.center) / trustregion.radius + x_unit = trustregion.map_to_unit(xs) - true_fvecs = vector_model.predict(centered_xs) + true_fvecs = vector_model.predict(x_unit) true_scalar_model = model_aggregator(vector_model=vector_model) @@ -74,7 +74,7 @@ def simulate_rho_noise( sim_scalar_model = model_aggregator(vector_model=sim_vector_model) sim_sub_sol = subsolver(sim_scalar_model, trustregion) - sim_candidate_fval = true_scalar_model.predict(sim_sub_sol.centered_x) + sim_candidate_fval = true_scalar_model.predict(sim_sub_sol.x_unit) sim_actual_improvement = -(sim_candidate_fval - true_current_fval) sim_rho = calculate_rho( diff --git a/src/estimagic/optimization/tranquilo/solve_subproblem.py b/src/estimagic/optimization/tranquilo/solve_subproblem.py index f9a0a9fe1..3ddeb3f15 100644 --- a/src/estimagic/optimization/tranquilo/solve_subproblem.py +++ b/src/estimagic/optimization/tranquilo/solve_subproblem.py @@ -176,7 +176,7 @@ def _solve_subproblem_template( raw_result = solver(model, **_bounds, **options) - x = _uncenter_and_unscale(raw_result["x"], trustregion) + x = trustregion.map_from_unit(raw_result["x"]) if "lower_bounds" in bounds: x = np.clip(x, bounds["lower_bounds"], np.inf) @@ -186,17 +186,18 @@ def _solve_subproblem_template( # make sure expected improvement is calculated accurately in case of clipping and # does not depend on whether the subsolver ignores intercepts or not. - fval_at_center = model.predict(np.zeros_like(x)) + fval_old = model.predict(trustregion.map_to_unit(trustregion.center)) fval_candidate = model.predict(raw_result["x"]) - expected_improvement = -(fval_candidate - fval_at_center) + expected_improvement = -(fval_candidate - fval_old) result = SubproblemResult( x=x, expected_improvement=expected_improvement, n_iterations=raw_result["n_iterations"], success=raw_result["success"], - centered_x=raw_result["x"], + x_unit=raw_result["x"], + shape=trustregion.shape, ) return result @@ -209,7 +210,7 @@ def _get_centered_and_scaled_bounds(bounds, trustregion): if bounds["lower_bounds"] is None: lower_bounds = np.full(n_params, -1) else: - lower_bounds = _center_and_scale(bounds["lower_bounds"], trustregion) + lower_bounds = trustregion.map_to_unit(bounds["lower_bounds"]) lower_bounds = np.nan_to_num(lower_bounds, nan=-1, neginf=-1) out["lower_bounds"] = lower_bounds @@ -217,20 +218,12 @@ def _get_centered_and_scaled_bounds(bounds, trustregion): if bounds["upper_bounds"] is None: upper_bounds = np.ones(n_params) else: - upper_bounds = _center_and_scale(bounds["upper_bounds"], trustregion) + upper_bounds = trustregion.map_to_unit(bounds["upper_bounds"]) upper_bounds = np.nan_to_num(upper_bounds, nan=1, posinf=1) out["upper_bounds"] = upper_bounds return out -def _center_and_scale(vec, trustregion): - return (vec - trustregion.center) / trustregion.radius - - -def _uncenter_and_unscale(vec, trustregion): - return vec * trustregion.radius + trustregion.center - - class SubproblemResult(NamedTuple): """Result of the subproblem solver.""" @@ -238,4 +231,5 @@ class SubproblemResult(NamedTuple): expected_improvement: float n_iterations: int success: bool - centered_x: np.ndarray + x_unit: np.ndarray + shape: str diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 7ea1a4b3e..b502c1203 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -58,7 +58,8 @@ def _internal_tranquilo( intercepts=_init_fvec, linear_terms=np.zeros((len(_init_fvec), len(x))), square_terms=np.zeros((len(_init_fvec), len(x), len(x))), - region=trustregion, + shift=trustregion.center, + scale=trustregion.radius, ) _init_model = aggregate_model(_init_vector_model) diff --git a/src/estimagic/optimization/tranquilo/tranquilo_history.py b/src/estimagic/optimization/tranquilo/tranquilo_history.py deleted file mode 100644 index 9dc8546e6..000000000 --- a/src/estimagic/optimization/tranquilo/tranquilo_history.py +++ /dev/null @@ -1,218 +0,0 @@ -import numpy as np -from numba import njit - - -class History: - """Container to save and retrieve history entries. - - These entries are: xs, fvecs and fvals. - - fvals don't need to be added explicitly, as they are computed internally - whenever new entries are added. - - """ - - def __init__(self, functype): - self.xs = None - self.fvecs = None - self.fvals = None - self.n_fun = 0 - - self.functype = functype - - if functype == "scalar": - self.aggregate = lambda x: x.flatten() - elif functype == "likelihood": - self.aggregate = lambda x: x.sum(axis=-1) - elif functype == "least_squares": - self.aggregate = lambda x: (x**2).sum(axis=-1) - else: - raise ValueError( - "funtype must be 'scalar', 'likelihood' or 'least_squares'." - ) - - def add_entries(self, xs, fvecs): - """Add new parameter vectors and fvecs to the history. - - Args: - xs (np.ndarray or list): 1d or 2d array or list of 1d arrays with - parameter vectors. - fvecs (np.ndarray or list): 1d or 2d array or list of 1d arrays with - least square fvecs. - - """ - xs = np.atleast_2d(xs) - - n_new_points = len(xs) if xs.size != 0 else 0 - - if n_new_points == 0: - return - - if self.functype == "scalar": - fvecs = np.reshape(fvecs, (-1, 1)) - else: - fvecs = np.atleast_2d(fvecs) - - fvals = np.atleast_1d(self.aggregate(fvecs)) - - if n_new_points != len(fvecs): - raise ValueError() - - self.xs = _add_entries_to_array(self.xs, xs, self.n_fun) - self.fvecs = _add_entries_to_array(self.fvecs, fvecs, self.n_fun) - self.fvals = _add_entries_to_array(self.fvals, fvals, self.n_fun) - - self.n_fun += len(xs) - - def get_entries(self, index=None): - """Retrieve xs, fvecs and fvals from the history. - - Args: - index (None, int or np.ndarray): Specifies the subset of rows that will - be returned. - - Returns: - np.ndarray: 1d or 2d array with parameter vectors. - np.ndarray: 1d or 2d array with fvecs. - np.ndarray: Float or 1d array with criterion values. - - """ - names = ["xs", "fvecs", "fvals"] - - out = (getattr(self, name)[: self.n_fun] for name in names) - - # Reducing arrays to length n_fun ensures that invalid indices raise IndexError - if index is not None: - out = [arr[index] for arr in out] - - return tuple(out) - - def get_xs(self, index=None): - """Retrieve xs from history. - - Args: - index (None, int or np.ndarray): Specifies the subset of rows that will - be returned. - - Returns: - np.ndarray: 1d or 2d array with parameter vectors - - """ - out = self.xs[: self.n_fun] - out = out[index] if index is not None else out - - return out - - def get_fvecs(self, index=None): - """Retrieve fvecs from history. - - Args: - index (None, int or np.ndarray): Specifies the subset of rows that will - be returned. - - Returns: - np.ndarray: 1d or 2d array with fvecs. - - """ - out = self.fvecs[: self.n_fun] - out = out[index] if index is not None else out - - return out - - def get_fvals(self, index=None): - """Retrieve fvals from history. - - Args: - index (None, int or np.ndarray): Specifies the subset of rows that will - be returned. - - Returns: - np.ndarray: Float or 1d array with criterion values. - - """ - out = self.fvals[: self.n_fun] - out = out[index] if index is not None else out - - return out - - def get_n_fun(self): - return self.n_fun - - def get_indices_in_region(self, region): - # early return if there are no entries - if self.get_n_fun() == 0: - return np.array([]) - - center = region.center - shape = region.shape - radius = region.radius - - if isinstance(center, int): - center = self.get_xs(index=center) - - xs = self.get_xs() - - out = _find_indices_in_trust_region(xs, center=center, radius=radius) - - if shape == "sphere": - mask = np.linalg.norm(xs[out] - region.center, axis=1) <= radius - out = out[mask] - - return out - - def __repr__(self): - return f"History for {self.functype} function with {self.n_fun} entries." - - -def _add_entries_to_array(arr, new, position): - if arr is None: - shape = 100_000 if new.ndim == 1 else (100_000, new.shape[1]) - arr = np.full(shape, np.nan) - - n_new_points = len(new) if new.size != 0 else 0 - - if len(arr) - position - n_new_points < 0: - n_extend = max(len(arr), n_new_points) - if arr.ndim == 2: - extension_shape = (n_extend, arr.shape[1]) - arr = np.vstack([arr, np.full(extension_shape, np.nan)]) - else: - arr = np.hstack([arr, np.full(n_extend, np.nan)]) - - arr[position : position + n_new_points] = new - - return arr - - -@njit -def _find_indices_in_trust_region(xs, center, radius): - """Get the row indices of all parameter vectors in a trust region. - - This is for square trust regions, i.e. balls in term of an infinity norm. - - Args: - xs (np.ndarray): 2d numpy array where each row is a parameter vector. - center (np.ndarray): 1d numpy array that marks the center of the trust region. - radius (float): Radius of the trust region. - - Returns: - np.ndarray: The indices of parameters in the trust region. - - """ - n_obs, dim = xs.shape - out = np.zeros(n_obs).astype(np.int64) - success_counter = 0 - upper = center + radius - lower = center - radius - for i in range(n_obs): - success = True - for j in range(dim): - value = xs[i, j] - if not (lower[j] <= value <= upper[j]) or np.isnan(value): - success = False - continue - if success: - out[success_counter] = i - success_counter += 1 - - return out[:success_counter] diff --git a/tests/optimization/test_history.py b/tests/optimization/test_pounders_history.py similarity index 98% rename from tests/optimization/test_history.py rename to tests/optimization/test_pounders_history.py index 6c7572c64..fc38905a9 100644 --- a/tests/optimization/test_history.py +++ b/tests/optimization/test_pounders_history.py @@ -1,7 +1,7 @@ """Test the history class for least-squares optimizers.""" import numpy as np import pytest -from estimagic.optimization.history import LeastSquaresHistory +from estimagic.optimization.pounders_history import LeastSquaresHistory from numpy.testing import assert_array_almost_equal as aaae ENTRIES = [ diff --git a/tests/optimization/test_pounders_unit.py b/tests/optimization/test_pounders_unit.py index db7fbc1f1..1c6e51e8d 100644 --- a/tests/optimization/test_pounders_unit.py +++ b/tests/optimization/test_pounders_unit.py @@ -8,7 +8,7 @@ import yaml from estimagic.batch_evaluators import joblib_batch_evaluator from estimagic.config import TEST_FIXTURES_DIR -from estimagic.optimization.history import LeastSquaresHistory +from estimagic.optimization.pounders_history import LeastSquaresHistory from estimagic.optimization.pounders_auxiliary import ( add_geomtery_points_to_make_main_model_fully_linear, create_initial_residual_model, diff --git a/tests/optimization/tranquilo/test_acceptance_decision.py b/tests/optimization/tranquilo/test_acceptance_decision.py index 88fb5169e..57457be2f 100644 --- a/tests/optimization/tranquilo/test_acceptance_decision.py +++ b/tests/optimization/tranquilo/test_acceptance_decision.py @@ -7,7 +7,7 @@ _get_acceptance_result, calculate_rho, ) -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.solve_subproblem import SubproblemResult from numpy.testing import assert_array_equal @@ -24,7 +24,8 @@ def subproblem_solution(): expected_improvement=1.0, n_iterations=1, success=True, - centered_x=None, + x_unit=None, + shape=None, ) return res diff --git a/tests/optimization/tranquilo/test_estimate_variance.py b/tests/optimization/tranquilo/test_estimate_variance.py index 067f5d5d5..8e7cd1793 100644 --- a/tests/optimization/tranquilo/test_estimate_variance.py +++ b/tests/optimization/tranquilo/test_estimate_variance.py @@ -3,7 +3,7 @@ from estimagic.optimization.tranquilo.estimate_variance import ( _estimate_variance_classic, ) -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History from estimagic.optimization.tranquilo.tranquilo import Region from numpy.testing import assert_array_almost_equal as aaae diff --git a/tests/optimization/tranquilo/test_filter_points.py b/tests/optimization/tranquilo/test_filter_points.py index 1c32783c0..97ded5a11 100644 --- a/tests/optimization/tranquilo/test_filter_points.py +++ b/tests/optimization/tranquilo/test_filter_points.py @@ -1,229 +1,48 @@ -import numpy as np -import pytest -from estimagic.optimization.tranquilo.filter_points import ( - _scaled_square_features, - drop_collinear_pounders, -) -from estimagic.optimization.tranquilo.region import Region +from estimagic.optimization.tranquilo.filter_points import get_sample_filter 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 -from numpy.testing import assert_equal +from estimagic.optimization.tranquilo.region import Region +from numpy.testing import assert_array_equal as aae +import pytest +import numpy as np @pytest.fixture() -def basic_case(): - x_accepted = np.array([0.16004745, 0.00572722, 0.01158929]) - radius = 0.0125 - - xs = np.array( - [ - [0.15, 0.008, 0.01], - [0.25, 0.008, 0.01], - [0.15, 0.108, 0.01], - [0.15, 0.008, 0.11], - [0.15961778, -0.07539625, 0.08766385], - [0.2, 0.00851182, -0.00302887], - [0.15049526, -0.04199751, 0.00993654], - [0.13739276, 0.00793654, -0.03838443], - [0.15046527, 0.00796766, 0.01269269], - [0.14986784, 0.00809919, 0.00927703], - [0.12530518, 0.00613383, 0.01349762], - [0.14987566, 0.0081864, 0.00937541], - [0.15076496, 0.00570962, 0.01295807], - [0.15074537, 0.00526659, 0.01240602], - [0.15069081, 0.00552219, 0.0121367], - [0.15067245, 0.00559504, 0.01191949], - [0.15141789, 0.0056498, 0.01210095], - [0.16317245, 0.00558118, 0.01208116], - [0.15692245, 0.00559149, 0.01197266], - [0.15379745, 0.00562833, 0.01182871], - [0.16004745, 0.00572722, 0.01158929], # 20 - ] - ) - indices = np.arange(len(xs)) - - trustregion = Region(center=x_accepted, radius=radius) - - state = State( - trustregion=trustregion, +def state(): + out = State( + trustregion=Region(center=np.ones(2), radius=0.3), model_indices=None, model=None, - index=20, - x=x_accepted, - fval=0, + vector_model=None, + candidate_index=5, + candidate_x=np.array([1.1, 1.2]), + index=2, + x=np.ones(2), + fval=15, rho=None, accepted=True, - new_indices=None, - old_indices_discarded=None, old_indices_used=None, - candidate_index=None, - candidate_x=None, - vector_model=None, - ) - - expected_indices = np.array([20, 19, 18, 17, 16, 15, 13, 12, 8, 5, 4, 3, 2, 1, 0]) - expected_xs = np.array( - [ - [0.16004745, 0.00572722, 0.01158929], - [0.15379745, 0.00562833, 0.01182871], - [0.15692245, 0.00559149, 0.01197266], - [0.16317245, 0.00558118, 0.01208116], - [0.15141789, 0.0056498, 0.01210095], - [0.15067245, 0.00559504, 0.01191949], - [0.15074537, 0.00526659, 0.01240602], - [0.15076496, 0.00570962, 0.01295807], - [0.15046527, 0.00796766, 0.01269269], - [0.2, 0.00851182, -0.00302887], - [0.15961778, -0.07539625, 0.08766385], - [0.15, 0.008, 0.11], - [0.15, 0.108, 0.01], - [0.25, 0.008, 0.01], - [0.15, 0.008, 0.01], - ] - ) - - return xs, indices, state, expected_xs, expected_indices - - -@pytest.fixture() -def reordered_case(basic_case): - _, indices, state, expected_xs, expected_indices = basic_case - state = state._replace(index=13) - - xs = np.array( - [ - [0.15, 0.008, 0.01], - [0.25, 0.008, 0.01], - [0.15, 0.108, 0.01], - [0.15, 0.008, 0.11], - [0.15961778, -0.07539625, 0.08766385], - [0.2, 0.00851182, -0.00302887], - [0.15049526, -0.04199751, 0.00993654], - [0.13739276, 0.00793654, -0.03838443], - [0.15046527, 0.00796766, 0.01269269], - [0.14986784, 0.00809919, 0.00927703], - [0.12530518, 0.00613383, 0.01349762], - [0.14987566, 0.0081864, 0.00937541], - [0.15076496, 0.00570962, 0.01295807], - [0.16004745, 0.00572722, 0.01158929], - [0.15074537, 0.00526659, 0.01240602], - [0.15069081, 0.00552219, 0.0121367], - [0.15067245, 0.00559504, 0.01191949], - [0.15141789, 0.0056498, 0.01210095], - [0.16317245, 0.00558118, 0.01208116], - [0.15692245, 0.00559149, 0.01197266], - [0.15379745, 0.00562833, 0.01182871], - ] - ) - - expected_indices = np.array([13, 20, 19, 18, 17, 16, 14, 12, 8, 5, 4, 3, 2, 1, 0]) - - return xs, indices, state, expected_xs, expected_indices - - -@pytest.fixture() -def truncated_case(reordered_case): - _, _, state, expected_xs, _ = reordered_case - - xs = np.array( - [ - [0.15049526, -0.04199751, 0.00993654], - [0.13739276, 0.00793654, -0.03838443], - [0.15046527, 0.00796766, 0.01269269], - [0.14986784, 0.00809919, 0.00927703], - [0.12530518, 0.00613383, 0.01349762], - [0.14987566, 0.0081864, 0.00937541], - [0.15076496, 0.00570962, 0.01295807], - [0.16004745, 0.00572722, 0.01158929], - [0.15074537, 0.00526659, 0.01240602], - [0.15067245, 0.00559504, 0.01191949], - [0.15141789, 0.0056498, 0.01210095], - [0.16317245, 0.00558118, 0.01208116], - [0.15692245, 0.00559149, 0.01197266], - [0.15379745, 0.00562833, 0.01182871], - ] - ) - indices = np.array([6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20]) - - expected_indices = np.array([13, 20, 19, 18, 17, 16, 14, 12, 8]) - expected_xs = expected_xs[:-6] - - return xs, indices, state, expected_xs, expected_indices - - -@pytest.fixture() -def sparse_case(reordered_case): - _, _, state, *_ = reordered_case - - xs = np.array( - [ - [0.15046527, 0.00796766, 0.01269269], - [0.14986784, 0.00809919, 0.00927703], - [0.12530518, 0.00613383, 0.01349762], - [0.14987566, 0.0081864, 0.00937541], - [0.15076496, 0.00570962, 0.01295807], - [0.16004745, 0.00572722, 0.01158929], - [0.15074537, 0.00526659, 0.01240602], - [0.15067245, 0.00559504, 0.01191949], - [0.15141789, 0.0056498, 0.01210095], - [0.15379745, 0.00562833, 0.01182871], - ] - ) - indices = np.array([8, 9, 10, 11, 12, 13, 14, 16, 17, 20]) - - expected_indices = np.array([13, 20, 17, 16, 14, 12, 11, 10, 9, 8]) - expected_xs = np.array( - [ - [0.16004745, 0.00572722, 0.01158929], - [0.15379745, 0.00562833, 0.01182871], - [0.15141789, 0.0056498, 0.01210095], - [0.15067245, 0.00559504, 0.01191949], - [0.15074537, 0.00526659, 0.01240602], - [0.15076496, 0.00570962, 0.01295807], - [0.14987566, 0.0081864, 0.00937541], - [0.12530518, 0.00613383, 0.01349762], - [0.14986784, 0.00809919, 0.00927703], - [0.15046527, 0.00796766, 0.01269269], - ] - ) - - return xs, indices, state, expected_xs, expected_indices - - -def test_indices_in_trust_region(basic_case): - functype = "scalar" - history = History(functype=functype) - - xs, *_ = basic_case - x_accepted = np.array([0.16004745, 0.00572722, 0.01158929]) - radius = 0.0125 - - 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) - - expected_indices = np.array([0, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]) - assert_equal(indices_in_tr, expected_indices) - - -@pytest.mark.parametrize( - "test_case", ["basic_case", "reordered_case", "truncated_case", "sparse_case"] -) -def test_drop_collinear_pounders(test_case, request): - old_xs, old_indices, state, expected_xs, expected_indices = request.getfixturevalue( - test_case + old_indices_discarded=None, + new_indices=None, + step_length=0.1, + relative_step_length=0.1 / 0.3, ) + return out - filtered_xs, filtered_indices = drop_collinear_pounders(old_xs, old_indices, state) - assert_equal(filtered_indices, expected_indices) - aaae(filtered_xs, expected_xs) +def test_discard_all(state): + filter = get_sample_filter("discard_all") + xs = np.arange(10).reshape(5, 2) + indices = np.arange(5) + got_xs, got_idxs = filter(xs=xs, indices=indices, state=state) + expected_xs = np.ones((1, 2)) + aae(got_xs, expected_xs) + aae(got_idxs, np.array([2])) -def test_scaled_square_features(): - x = np.arange(4).reshape(2, 2) - got = _scaled_square_features(x) - expected = np.array([[0, 0, 1 / 2], [2, 6 / np.sqrt(2), 9 / 2]]) - assert_equal(got, expected) +def test_keep_all(): + filter = get_sample_filter("keep_all") + xs = np.arange(10).reshape(5, 2) + indices = np.arange(5) + got_xs, got_idxs = filter(xs=xs, indices=indices, state=None) + aae(got_xs, xs) + aae(got_idxs, indices) diff --git a/tests/optimization/tranquilo/test_fit_models.py b/tests/optimization/tranquilo/test_fit_models.py index a4c6be5f1..de5b9f065 100644 --- a/tests/optimization/tranquilo/test_fit_models.py +++ b/tests/optimization/tranquilo/test_fit_models.py @@ -6,13 +6,13 @@ from numpy.testing import assert_array_almost_equal, assert_array_equal -def aaae(x, y, case=None): +def aaae(x, y, decimal=None, case=None): tolerance = { - None: 3, "hessian": 2, "gradient": 3, } - assert_array_almost_equal(x, y, decimal=tolerance[case]) + decimal = decimal or tolerance.get(case, None) + assert_array_almost_equal(x, y, decimal=decimal) # ====================================================================================== @@ -62,41 +62,45 @@ def func(x): # ====================================================================================== -def test_fit_ols_against_truth(quadratic_case): - fit_ols = get_fitter("ols", model_type="quadratic") - got = fit_ols( - x=quadratic_case["x"], - y=quadratic_case["y"], - region=Region(center=np.zeros(4), radius=1.0), - old_model=None, +@pytest.mark.parametrize("fitter", ["ols", "ridge", "powell", "tranquilo"]) +def test_fit_against_truth_quadratic(fitter, quadratic_case): + options = {"l2_penalty_square": 0} + fit_pounders = get_fitter( + fitter, + options, + model_type="quadratic", + residualize=False, + infinity_handling="relative", ) - - n_p = got.linear_terms.size - aaae(got.linear_terms.flatten(), quadratic_case["linear_terms_expected"]) - aaae(got.square_terms.reshape(n_p, n_p), quadratic_case["square_terms_expected"]) - - -def test_fit_powell_against_truth_quadratic(quadratic_case): - fit_pounders = get_fitter("powell", model_type="quadratic") got = fit_pounders( quadratic_case["x"], quadratic_case["y"], region=Region(center=np.zeros(4), radius=1.0), old_model=None, ) - - aaae(got.linear_terms.flatten(), quadratic_case["linear_terms_expected"]) - aaae(got.square_terms.reshape((4, 4)), quadratic_case["square_terms_expected"]) + decimal = 3 if fitter != "ridge" else 1 + aaae( + got.linear_terms.flatten(), + quadratic_case["linear_terms_expected"], + decimal=decimal, + ) + aaae( + got.square_terms.reshape((4, 4)), + quadratic_case["square_terms_expected"], + decimal=decimal, + ) -@pytest.mark.parametrize("model", ["ols", "ridge"]) +@pytest.mark.parametrize("model", ["ols", "ridge", "tranquilo"]) def test_fit_ols_against_gradient(model, quadratic_case): - if model == "ridge": - options = {"l2_penalty_square": 0} - else: - options = None - - fit_ols = get_fitter(model, options, model_type="quadratic") + options = {"l2_penalty_square": 0} + fit_ols = get_fitter( + model, + options, + model_type="quadratic", + residualize=False, + infinity_handling="relative", + ) got = fit_ols( quadratic_case["x"], quadratic_case["y"], @@ -112,12 +116,16 @@ def test_fit_ols_against_gradient(model, quadratic_case): aaae(gradient["derivative"], grad, case="gradient") -@pytest.mark.parametrize( - "model, options", - [("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, model_type="quadratic") +@pytest.mark.parametrize("model", ("ols", "ridge", "tranquilo", "powell")) +def test_fit_ols_against_hessian(model, quadratic_case): + options = {"l2_penalty_square": 0} + fit_ols = get_fitter( + model, + options, + model_type="quadratic", + residualize=False, + infinity_handling="relative", + ) got = fit_ols( quadratic_case["x"], quadratic_case["y"], diff --git a/tests/optimization/tranquilo/test_new_history.py b/tests/optimization/tranquilo/test_history.py similarity index 96% rename from tests/optimization/tranquilo/test_new_history.py rename to tests/optimization/tranquilo/test_history.py index b36a5f13b..bf2410df5 100644 --- a/tests/optimization/tranquilo/test_new_history.py +++ b/tests/optimization/tranquilo/test_history.py @@ -1,11 +1,11 @@ """Test the history class for least-squares optimizers.""" -from collections import namedtuple - import numpy as np import pytest -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History +from estimagic.optimization.tranquilo.region import Region from numpy.testing import assert_array_almost_equal as aaae + XS = [ np.arange(3), np.arange(3).tolist(), @@ -117,10 +117,9 @@ def test_get_indices_in_trustregion(): indices = history.add_xs(xs) history.add_evals(indices, fvecs) - trustregion = namedtuple("Region", ["center", "radius", "shape"])( + trustregion = Region( center=np.ones(2), radius=0.3, - shape="sphere", ) indices = history.get_x_indices_in_region(trustregion) diff --git a/tests/optimization/tranquilo/test_models.py b/tests/optimization/tranquilo/test_models.py index 12899c958..cffffa93e 100644 --- a/tests/optimization/tranquilo/test_models.py +++ b/tests/optimization/tranquilo/test_models.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.models import ( ScalarModel, VectorModel, @@ -11,9 +12,7 @@ n_free_params, n_interactions, n_second_order_terms, - subtract_models, ) -from estimagic.optimization.tranquilo.region import Region from numpy.testing import assert_array_almost_equal as aaae from numpy.testing import assert_array_equal @@ -23,7 +22,6 @@ def test_predict_scalar(): intercept=1.0, linear_terms=np.arange(2), square_terms=(np.arange(4) + 1).reshape(2, 2), - region=None, ) x = np.array([[0, 0], [0, 1], [1, 0], [1, 2]]) exp = np.array([1, 4, 1.5, 16.5]) @@ -36,7 +34,6 @@ def test_predict_vector(): intercepts=1 + np.arange(3), linear_terms=np.arange(6).reshape(3, 2), square_terms=(np.arange(3 * 2 * 2) + 1).reshape(3, 2, 2), - region=None, ) x = np.array([[0, 0], [0, 1], [1, 0], [1, 2]], dtype=float) exp = np.array( @@ -112,7 +109,8 @@ 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), + shift=np.array([0.2, 0.3]), + scale=0.6, ) return out @@ -129,48 +127,49 @@ def vector_model(): [[0.8, 0.2], [0.2, 0.7]], ] ), - region=Region(center=np.array([0.2, 0.3]), radius=0.6), + shift=np.array([0.2, 0.3]), + scale=0.6, ) return out def test_move_scalar_model(scalar_model): - old_region = scalar_model.region + old_region = Region(center=scalar_model.shift, radius=scalar_model.scale) 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]]) - x_old = (x_unscaled - old_region.center) / old_region.radius - x_new = (x_unscaled - new_region.center) / new_region.radius + x_old = old_region.map_to_unit(x_unscaled) + x_new = new_region.map_to_unit(x_unscaled) new_model = move_model(old_model, new_region) old_prediction = old_model.predict(x_old) new_prediction = new_model.predict(x_new) - assert new_model.region.radius == new_region.radius - aaae(new_model.region.center, new_region.center) + assert new_model.scale == new_region.radius + aaae(new_model.shift, new_region.center) assert np.allclose(old_prediction, new_prediction) def test_move_vector_model(vector_model): - old_region = vector_model.region + old_region = Region(center=vector_model.shift, radius=vector_model.scale) new_region = Region(center=np.array([-0.1, 0.1]), radius=0.45) old_model = vector_model x_unscaled = np.array([[0.5, 0.5]]) - x_old = (x_unscaled - old_region.center) / old_region.radius - x_new = (x_unscaled - new_region.center) / new_region.radius + x_old = old_region.map_to_unit(x_unscaled) + x_new = new_region.map_to_unit(x_unscaled) new_model = move_model(old_model, new_region) old_prediction = old_model.predict(x_old) new_prediction = new_model.predict(x_new) - assert new_model.region.radius == new_region.radius - aaae(new_model.region.center, new_region.center) + assert new_model.scale == new_region.radius + aaae(new_model.shift, new_region.center) assert np.allclose(old_prediction, new_prediction) @@ -189,19 +188,3 @@ def test_add_vector_models(vector_model): assert np.allclose(got.intercepts, vector_model.intercepts * 2) aaae(got.linear_terms, vector_model.linear_terms * 2) aaae(got.square_terms, vector_model.square_terms * 2) - - -def test_subtract_scalar_model(scalar_model): - got = subtract_models(scalar_model, scalar_model) - - assert got.intercept == 0.0 - aaae(got.linear_terms, np.zeros_like(scalar_model.linear_terms)) - aaae(got.square_terms, np.zeros_like(scalar_model.square_terms)) - - -def test_subtract_vector_model(vector_model): - got = subtract_models(vector_model, vector_model) - - assert np.allclose(got.intercepts, np.zeros_like(vector_model.intercepts)) - aaae(got.linear_terms, np.zeros_like(vector_model.linear_terms)) - aaae(got.square_terms, np.zeros_like(vector_model.square_terms)) diff --git a/tests/optimization/tranquilo/test_process_arguments.py b/tests/optimization/tranquilo/test_process_arguments.py index 52786e963..8532e65ea 100644 --- a/tests/optimization/tranquilo/test_process_arguments.py +++ b/tests/optimization/tranquilo/test_process_arguments.py @@ -14,6 +14,7 @@ _process_search_radius_factor, _process_acceptance_decider, _process_model_fitter, + _process_residualize, _process_n_evals_at_start, ) @@ -106,12 +107,6 @@ def test_process_model_fitter(): ) == "ols" ) - assert ( - _process_model_fitter( - model_fitter=None, model_type="linear", sample_size=2, x=np.arange(3) - ) - == "ridge" - ) assert ( _process_model_fitter( model_fitter="xyz", model_type=None, sample_size=None, x=None @@ -120,6 +115,17 @@ def test_process_model_fitter(): ) +def test_process_residualize(): + assert _process_residualize(residualize=None, model_fitter="tranquilo") is True + assert _process_residualize(residualize=None, model_fitter="ols") is False + assert _process_residualize(residualize=False, model_fitter="custom") is False + + +def test_process_residualize_invalid(): + with pytest.raises(ValueError, match="residualize must be a boolean."): + _process_residualize(residualize="invalid", model_fitter=None) + + def test_process_n_evals_at_start(): assert _process_n_evals_at_start(n_evals=None, noisy=True) == 5 assert _process_n_evals_at_start(n_evals=None, noisy=False) == 1 diff --git a/tests/optimization/tranquilo/test_region.py b/tests/optimization/tranquilo/test_region.py index 496036475..f4643a005 100644 --- a/tests/optimization/tranquilo/test_region.py +++ b/tests/optimization/tranquilo/test_region.py @@ -1,17 +1,20 @@ import numpy as np -import pytest from estimagic.optimization.tranquilo.bounds import Bounds from estimagic.optimization.tranquilo.region import ( Region, _any_bounds_binding, + _get_shape, _get_cube_bounds, _get_cube_center, + _get_effective_radius, + _get_effective_center, _map_from_unit_cube, _map_from_unit_sphere, _map_to_unit_cube, _map_to_unit_sphere, ) from numpy.testing import assert_array_equal +import pytest def test_map_to_unit_sphere(): @@ -48,41 +51,69 @@ def test_any_bounds_binding_false(): assert not out +def test_get_shape_sphere(): + out = _get_shape(center=np.zeros(2), radius=1, bounds=None) + assert out == "sphere" + + +def test_get_shape_cube(): + bounds = Bounds(lower=np.zeros(2), upper=np.ones(2)) + out = _get_shape(center=np.zeros(2), radius=1, bounds=bounds) + assert out == "cube" + + 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) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds, shape="sphere") 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) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds, shape="sphere") 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)) + out = _get_cube_bounds(center=np.zeros(2), radius=1, bounds=bounds, shape="cube") + np.all(out.lower > -np.ones(2)) + np.all(out.lower < np.zeros(2)) + np.all(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) + out = _get_cube_center(cube_bounds=bounds) assert_array_equal(out, np.array([0.5, 5.25])) +def test_get_effective_radius(): + bounds = Bounds(lower=np.array([0, 0.5]), upper=np.array([1, 10])) + out = _get_effective_radius(shape="cube", radius=None, cube_bounds=bounds) + assert_array_equal(out, np.array([0.5, 4.75])) + + +def test_get_effective_center_sphere(): + out = _get_effective_center(shape="sphere", center=np.ones(2), cube_center=None) + assert_array_equal(out, np.ones(2)) + + +def test_get_effective_center_cube(): + out = _get_effective_center(shape="cube", center=None, cube_center=np.zeros(2)) + assert_array_equal(out, np.zeros(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): + with pytest.raises(AttributeError, match="The trustregion is a sphere"): region.cube_bounds - with pytest.raises(AttributeError): + with pytest.raises(AttributeError, match="The trustregion is a sphere"): region.cube_center diff --git a/tests/optimization/tranquilo/test_rho_noise.py b/tests/optimization/tranquilo/test_rho_noise.py index bf2db04a8..4e2a145ff 100644 --- a/tests/optimization/tranquilo/test_rho_noise.py +++ b/tests/optimization/tranquilo/test_rho_noise.py @@ -46,7 +46,12 @@ 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.0) - model_fitter = get_fitter(fitter="ols", model_type=model_type) + model_fitter = get_fitter( + fitter="ols", + model_type=model_type, + residualize=False, + infinity_handling="relative", + ) vector_model = model_fitter( xs, fvecs, weights=None, region=trustregion, old_model=None diff --git a/tests/optimization/tranquilo/test_solve_subproblem.py b/tests/optimization/tranquilo/test_solve_subproblem.py index 64b692344..1acda4bd0 100644 --- a/tests/optimization/tranquilo/test_solve_subproblem.py +++ b/tests/optimization/tranquilo/test_solve_subproblem.py @@ -1,8 +1,7 @@ -from collections import namedtuple - import numpy as np import pytest from estimagic.optimization.tranquilo.models import ScalarModel +from estimagic.optimization.tranquilo.region import Region from estimagic.optimization.tranquilo.solve_subproblem import get_subsolver from numpy.testing import assert_array_almost_equal as aaae @@ -33,7 +32,7 @@ def test_without_bounds(solver_name): intercept=0, linear_terms=linear_terms, square_terms=quadratic_terms ) - trustregion = namedtuple("Trustregion", ["center", "radius"])( + trustregion = Region( center=np.zeros(3), radius=1, ) diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index 4f720a6b6..e628a3d1d 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -129,9 +129,9 @@ def test_external_tranquilo_scalar_sphere_defaults(): "model_fitter": ["ols"], "model_type": ["linear"], }, - "pounders_filtering": { - "sample_filter": ["drop_pounders"], - "model_fitter": ["ols"], + "tranquilo": { + "sample_filter": ["keep_all", "discard_all"], + "model_fitter": ["tranquilo"], "model_type": ["linear"], }, } @@ -207,3 +207,28 @@ def _f(x): ) aaae(res.params, np.zeros(3), decimal=1) + + +# ====================================================================================== +# Bounded case +# ====================================================================================== + + +def sum_of_squares(x): + contribs = x**2 + return {"value": contribs.sum(), "contributions": contribs, "root_contributions": x} + + +@pytest.mark.parametrize("algorithm", ["tranquilo", "tranquilo_ls"]) +def test_tranquilo_with_binding_bounds(algorithm): + res = minimize( + criterion=sum_of_squares, + params=np.array([3, 2, -3]), + lower_bounds=np.array([1, -np.inf, -np.inf]), + upper_bounds=np.array([np.inf, np.inf, -1]), + algorithm=algorithm, + collect_history=True, + skip_checks=True, + ) + assert res.success in [True, None] + aaae(res.params, np.array([1, 0, -1]), decimal=3) diff --git a/tests/optimization/tranquilo/test_tranquilo_history.py b/tests/optimization/tranquilo/test_tranquilo_history.py deleted file mode 100644 index 536dbe45a..000000000 --- a/tests/optimization/tranquilo/test_tranquilo_history.py +++ /dev/null @@ -1,112 +0,0 @@ -"""Test the history class for least-squares optimizers.""" -from collections import namedtuple - -import numpy as np -import pytest -from estimagic.optimization.tranquilo.tranquilo_history import History -from numpy.testing import assert_array_almost_equal as aaae - -ENTRIES = [ - (np.arange(3), [np.arange(5)]), - ([np.arange(3)], list(range(5))), - (np.arange(3).reshape(1, 3), np.arange(5).reshape(1, 5)), -] - -TEST_CASES = [] -for entries in ENTRIES: - # leave it in in case centered stuff comes back - for is_center in [False]: - TEST_CASES.append((entries, is_center)) - - -@pytest.mark.parametrize("entries, is_center", TEST_CASES) -def test_add_entries_not_initialized(entries, is_center): - history = History(functype="least_squares") - - if is_center: - c_info = {"x": np.zeros(3), "fvecs": np.zeros(5), "radius": 1} - history.add_centered_entries(*entries, c_info) - else: - history.add_entries(*entries) - - xs, fvecs, fvals = history.get_entries() - xs_sinlge = history.get_xs() - fvecs_sinlge = history.get_fvecs() - fvals_sinlge = history.get_fvals() - - for entry in xs, fvecs, fvals: - assert isinstance(entry, np.ndarray) - - aaae(xs, np.arange(3).reshape(1, 3)) - aaae(xs_sinlge, np.arange(3).reshape(1, 3)) - aaae(fvecs, np.arange(5).reshape(1, 5)) - aaae(fvecs_sinlge, np.arange(5).reshape(1, 5)) - aaae(fvals, np.array([30.0])) - aaae(fvals_sinlge, np.array([30.0])) - - -@pytest.mark.parametrize("entries, is_center", TEST_CASES) -def test_add_entries_initialized_with_space(entries, is_center): - history = History(functype="least_squares") - history.add_entries(np.ones((4, 3)), np.zeros((4, 5))) - - if is_center: - c_info = {"x": np.zeros(3), "fvecs": np.zeros(5), "radius": 1} - history.add_centered_entries(*entries, c_info) - else: - history.add_entries(*entries) - - xs, fvecs, fvals = history.get_entries(index=-1) - xs_sinlge = history.get_xs(index=-1) - fvecs_sinlge = history.get_fvecs(index=-1) - fvals_sinlge = history.get_fvals(index=-1) - - for entry in xs, fvecs: - assert isinstance(entry, np.ndarray) - - aaae(xs, np.arange(3)) - aaae(xs_sinlge, np.arange(3)) - aaae(fvecs, np.arange(5)) - aaae(fvecs_sinlge, np.arange(5)) - assert fvals == 30 - assert fvals_sinlge == 30 - - -def test_add_entries_initialized_extension_needed(): - history = History(functype="least_squares") - history.add_entries(np.ones((4, 3)), np.zeros((4, 5))) - history.xs = history.xs[:5] - history.fvecs = history.fvecs[:5] - history.fvals = history.fvals[:5] - - history.add_entries(np.arange(12).reshape(4, 3), np.arange(20).reshape(4, 5)) - - assert len(history.xs) == 10 - assert len(history.fvecs) == 10 - assert len(history.fvals) == 10 - - xs, fvecs, _ = history.get_entries(index=-1) - xs_sinlge = history.get_xs(index=-1) - fvecs_sinlge = history.get_fvecs(index=-1) - - for entry in xs, xs_sinlge, fvecs, fvecs_sinlge: - assert isinstance(entry, np.ndarray) - - assert history.get_n_fun() == 8 - - -def test_get_indices_in_trustregion(): - history = History(functype="least_squares") - xs = [[1, 1], [1.1, 1.2], [1.5, 1], [0.9, 0.9]] - fvecs = np.zeros((4, 3)) - history.add_entries(xs, fvecs) - - trustregion = namedtuple("TrustRegion", ["center", "radius", "shape"])( - center=np.ones(2), - radius=0.3, - shape="sphere", - ) - - indices = history.get_indices_in_region(trustregion) - - aaae(indices, np.array([0, 1, 3])) diff --git a/tests/optimization/tranquilo/test_wrap_criterion.py b/tests/optimization/tranquilo/test_wrap_criterion.py index e9e4b04f0..3886a36cd 100644 --- a/tests/optimization/tranquilo/test_wrap_criterion.py +++ b/tests/optimization/tranquilo/test_wrap_criterion.py @@ -2,7 +2,7 @@ import numpy as np import pytest -from estimagic.optimization.tranquilo.new_history import History +from estimagic.optimization.tranquilo.history import History from estimagic.optimization.tranquilo.wrap_criterion import get_wrapped_criterion from numpy.testing import assert_array_almost_equal as aaae diff --git a/tests/visualization/test_visualize_tranquilo.py b/tests/visualization/test_visualize_tranquilo.py index ed4e3b7eb..9d2d2dfd3 100644 --- a/tests/visualization/test_visualize_tranquilo.py +++ b/tests/visualization/test_visualize_tranquilo.py @@ -7,13 +7,13 @@ "random_hull": { "sampler": "random_hull", "subsolver": "gqtpar_fast", - "sample_filter": "drop_pounders", + "sample_filter": "keep_all", "stopping.max_iterations": 10, }, "optimal_hull": { "sampler": "optimal_hull", "subsolver": "gqtpar_fast", - "sample_filter": "drop_pounders", + "sample_filter": "keep_all", "stopping.max_iterations": 10, }, }