From d3c49717a2f7405acbaf9d041e39838665ff75b3 Mon Sep 17 00:00:00 2001 From: Janos Gabler Date: Tue, 1 Nov 2022 14:25:29 +0100 Subject: [PATCH 01/11] Add some helper functions to models.py. --- .../tranquilo/aggregate_models.py | 2 +- .../optimization/tranquilo/models.py | 42 ++++++++++++++++ tests/optimization/tranquilo/test_models.py | 50 +++++++++++++++++++ 3 files changed, 93 insertions(+), 1 deletion(-) create mode 100644 tests/optimization/tranquilo/test_models.py diff --git a/src/estimagic/optimization/tranquilo/aggregate_models.py b/src/estimagic/optimization/tranquilo/aggregate_models.py index 27cc5fd50..6341524a3 100644 --- a/src/estimagic/optimization/tranquilo/aggregate_models.py +++ b/src/estimagic/optimization/tranquilo/aggregate_models.py @@ -46,7 +46,7 @@ def get_aggregator(aggregator, functype, model_info): # determine if aggregator is compatible with functype and model_info aggregator_compatible_with_functype = { - "scalar": ("identity", "identity_old", "sum"), + "scalar": ("identity", "sum"), "least_squares": ("least_squares_linear",), "likelihood": ( "sum", diff --git a/src/estimagic/optimization/tranquilo/models.py b/src/estimagic/optimization/tranquilo/models.py index 7da5ccf33..43f211b38 100644 --- a/src/estimagic/optimization/tranquilo/models.py +++ b/src/estimagic/optimization/tranquilo/models.py @@ -53,3 +53,45 @@ def evaluate_model(scalar_model, centered_x): y += x.T @ scalar_model.square_terms @ x / 2 return y + + +def n_free_params(dim, info_or_name): + """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): + name = info_or_name + if name == "quadratic": + out += n_second_order_terms(dim) + elif name == "diagonal": + out += dim + else: + raise ValueError() + + return out + + +def n_second_order_terms(dim): + """Number of free second order terms in a quadratic model.""" + return int(dim * (dim + 1) * 0.5) + + +def n_interactions(dim): + """Number of free interaction terms in a quadratic model.""" + return int(dim * (dim - 1) * 0.5) + + +def is_second_order_model(model_or_info): + """Check if a model has any second order terms.""" + if isinstance(model_or_info, ModelInfo): + out = model_or_info.has_interactions or model_or_info.has_squares + elif model_or_info.square_terms is not None: + out = True + else: + raise TypeError() + return out diff --git a/tests/optimization/tranquilo/test_models.py b/tests/optimization/tranquilo/test_models.py new file mode 100644 index 000000000..4e3479ae7 --- /dev/null +++ b/tests/optimization/tranquilo/test_models.py @@ -0,0 +1,50 @@ +from estimagic.optimization.tranquilo.models import n_interactions +from estimagic.optimization.tranquilo.models import n_second_order_terms + + +def test_n_free_params_name_quadratic(): + pass + + +def test_n_free_params_name_diagonal(): + pass + + +def test_n_free_params_name_invalid(): + pass + + +def test_n_free_params_info_linear(): + pass + + +def test_n_free_params_info_diagonal(): + pass + + +def test_n_free_params_info_quadratic(): + pass + + +def test_n_free_params_invalid(): + pass + + +def test_n_second_order_terms(): + assert n_second_order_terms(3) == 6 + + +def test_n_interactions(): + assert n_interactions(3) == 3 + + +def test_is_second_order_model_info(): + pass + + +def test_is_second_order_model_model(): + pass + + +def test_is_second_order_model_invalid(): + pass From 532f30b714989632293549a055ba3f894f8998f3 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 10:52:56 +0100 Subject: [PATCH 02/11] Add functions for calculation of model degrees of freedom --- .../optimization/tranquilo/filter_points.py | 5 +- .../optimization/tranquilo/models.py | 22 +++++--- .../optimization/tranquilo/tranquilo.py | 6 +- tests/optimization/tranquilo/test_models.py | 56 ++++++++++++++----- 4 files changed, 63 insertions(+), 26 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/filter_points.py b/src/estimagic/optimization/tranquilo/filter_points.py index 7ea1ccc0d..921c6e27a 100644 --- a/src/estimagic/optimization/tranquilo/filter_points.py +++ b/src/estimagic/optimization/tranquilo/filter_points.py @@ -1,4 +1,5 @@ import numpy as np +from estimagic.optimization.tranquilo.models import n_second_order_terms from numba import njit from scipy.linalg import qr_multiply @@ -62,7 +63,7 @@ def drop_collinear_pounders(xs, indices, state): def _drop_collinear_pounders(xs, indices, state): theta2 = 1e-4 n_samples, n_params = xs.shape - n_poly_terms = n_params * (n_params + 1) // 2 + n_poly_terms = n_second_order_terms(n_params) indices_reverse = indices[::-1] indexer_reverse = np.arange(n_samples)[::-1] @@ -147,7 +148,7 @@ def _get_polynomial_feature_matrices( @njit def _square_features(x): n_samples, n_params = np.atleast_2d(x).shape - n_poly_terms = n_params * (n_params + 1) // 2 + n_poly_terms = n_second_order_terms(n_params) poly_terms = np.empty((n_poly_terms, n_samples), x.dtype) xt = x.T diff --git a/src/estimagic/optimization/tranquilo/models.py b/src/estimagic/optimization/tranquilo/models.py index 43f211b38..0440d272d 100644 --- a/src/estimagic/optimization/tranquilo/models.py +++ b/src/estimagic/optimization/tranquilo/models.py @@ -2,6 +2,7 @@ from typing import Union import numpy as np +from numba import njit class VectorModel(NamedTuple): @@ -64,7 +65,11 @@ def n_free_params(dim, info_or_name): out += dim if info.has_interactions: out += n_interactions(dim) - elif isinstance(info_or_name, str): + elif isinstance(info_or_name, str) and info_or_name in ( + "linear", + "quadratic", + "diagonal", + ): name = info_or_name if name == "quadratic": out += n_second_order_terms(dim) @@ -72,26 +77,29 @@ def n_free_params(dim, info_or_name): out += dim else: raise ValueError() - return out +@njit def n_second_order_terms(dim): """Number of free second order terms in a quadratic model.""" - return int(dim * (dim + 1) * 0.5) + return dim * (dim + 1) // 2 +@njit def n_interactions(dim): """Number of free interaction terms in a quadratic model.""" - return int(dim * (dim - 1) * 0.5) + return dim * (dim - 1) // 2 def is_second_order_model(model_or_info): """Check if a model has any second order terms.""" if isinstance(model_or_info, ModelInfo): - out = model_or_info.has_interactions or model_or_info.has_squares - elif model_or_info.square_terms is not None: - out = True + model_info = model_or_info + out = model_info.has_interactions or model_info.has_squares + elif isinstance(model_or_info, (ScalarModel, VectorModel)): + model = model_or_info + out = model.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 65a0cf6a5..ca9dc79b7 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -10,6 +10,7 @@ from estimagic.optimization.tranquilo.filter_points import get_sample_filter from estimagic.optimization.tranquilo.fit_models import get_fitter from estimagic.optimization.tranquilo.models import ModelInfo +from estimagic.optimization.tranquilo.models import n_free_params from estimagic.optimization.tranquilo.models import ScalarModel from estimagic.optimization.tranquilo.options import Bounds from estimagic.optimization.tranquilo.options import ConvOptions @@ -292,7 +293,6 @@ def _calculate_rho(actual_improvement, expected_improvement): rho = -np.inf else: rho = actual_improvement / expected_improvement - return rho @@ -371,11 +371,11 @@ def _process_sample_size(user_sample_size, model_info, x): elif isinstance(user_sample_size, str): user_sample_size = user_sample_size.replace(" ", "") if user_sample_size in ["linear", "n+1"]: - out = len(x) + 1 + out = n_free_params(dim=len(x), info_or_name="linear") elif user_sample_size in ["powell", "2n+1", "2*n+1"]: out = 2 * len(x) + 1 elif user_sample_size == "quadratic": - out = len(x) * (len(x) + 1) // 2 + len(x) + 1 + out = n_free_params(dim=len(x), info_or_name="quadratic") else: raise ValueError(f"Invalid sample size: {user_sample_size}") diff --git a/tests/optimization/tranquilo/test_models.py b/tests/optimization/tranquilo/test_models.py index 4e3479ae7..4b99aaf63 100644 --- a/tests/optimization/tranquilo/test_models.py +++ b/tests/optimization/tranquilo/test_models.py @@ -1,33 +1,52 @@ +import numpy as np +import pytest +from estimagic.optimization.tranquilo.models import is_second_order_model +from estimagic.optimization.tranquilo.models import ModelInfo +from estimagic.optimization.tranquilo.models import n_free_params from estimagic.optimization.tranquilo.models import n_interactions from estimagic.optimization.tranquilo.models import n_second_order_terms +from estimagic.optimization.tranquilo.models import ScalarModel def test_n_free_params_name_quadratic(): - pass + 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(): - pass + 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 def test_n_free_params_name_invalid(): - pass + with pytest.raises(ValueError): + assert n_free_params(dim=3, info_or_name="invalid") -def test_n_free_params_info_linear(): - pass +@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 -def test_n_free_params_info_diagonal(): - pass +@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 -def test_n_free_params_info_quadratic(): - pass +@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) def test_n_free_params_invalid(): - pass + 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) def test_n_second_order_terms(): @@ -38,13 +57,22 @@ def test_n_interactions(): assert n_interactions(3) == 3 -def test_is_second_order_model_info(): - pass +@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 def test_is_second_order_model_model(): - pass + model = ScalarModel(intercept=1.0, linear_terms=np.ones(1)) + assert is_second_order_model(model) is False + + model = ScalarModel(intercept=1.0, linear_terms=np.ones(1), square_terms=np.ones(1)) + assert is_second_order_model(model) is True def test_is_second_order_model_invalid(): - pass + model = np.linalg.lstsq + with pytest.raises(TypeError): + is_second_order_model(model) From 5859cc61d3a6a65a5ace971105cd171586f66057 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 12:28:28 +0100 Subject: [PATCH 03/11] Remove unused arguments in aggregator interface --- .../tranquilo/aggregate_models.py | 44 +++-------- .../optimization/tranquilo/tranquilo.py | 1 - .../tranquilo/test_aggregate_models.py | 76 ++++++++----------- 3 files changed, 44 insertions(+), 77 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/aggregate_models.py b/src/estimagic/optimization/tranquilo/aggregate_models.py index 6341524a3..05943029d 100644 --- a/src/estimagic/optimization/tranquilo/aggregate_models.py +++ b/src/estimagic/optimization/tranquilo/aggregate_models.py @@ -9,18 +9,14 @@ def get_aggregator(aggregator, functype, model_info): Args: aggregator (str or callable): Name of an aggregator or aggregator function. - The function must take as arguments (in that order): + The function must take as argument: - vector_model (VectorModel): A fitted vector model. - - fvec_center (np.ndarray): A 1d array of the residuals at the center of the - trust-region. In the noisy case, this may be an average. - - model_info (ModelInfo): The model information. functype (str): One of "scalar", "least_squares" and "likelihood". model_info (ModelInfo): Information that describes the functional form of the model. Returns: - callable: The partialled aggregator that only depends on vector_model and - fvec_center. + callable: The partialled aggregator that only depends on vector_model. """ built_in_aggregators = { @@ -83,45 +79,29 @@ def get_aggregator(aggregator, functype, model_info): ) # create aggregator - out = partial( - _aggregate_models_template, aggregator=_aggregator, model_info=model_info - ) + out = partial(_aggregate_models_template, aggregator=_aggregator) return out -def _aggregate_models_template(vector_model, fvec_center, aggregator, model_info): +def _aggregate_models_template(vector_model, aggregator): """Aggregate a VectorModel into a ScalarModel. - Note on fvec_center: - -------------------- - Let x0 be the x-value at which the x-sample is centered. If there is little noise - and the criterion function f is evaluated at x0, then fvec_center = f(x0). If, - however, the criterion function is very noisy or only evaluated in a neighborhood - around x0, then fvec_center is constructed as an average over evaluations of f - with x close to x0. - Args: vector_model (VectorModel): The VectorModel to aggregate. - fvec_center (np.ndarray): A 1d array of the residuals at the center of the - trust-region. In the noisy case, this may be an average. aggregator (callable): The function that does the actual aggregation. - model_info (ModelInfo): Information that describes the functional form of - the model. Returns: ScalarModel: The aggregated model """ - intercept, linear_terms, square_terms = aggregator( - vector_model, fvec_center, model_info - ) + intercept, linear_terms, square_terms = aggregator(vector_model) scalar_model = ScalarModel( intercept=intercept, linear_terms=linear_terms, square_terms=square_terms ) return scalar_model -def aggregator_identity(vector_model, fvec_center, model_info): +def aggregator_identity(vector_model): """Aggregate quadratic VectorModel using identity function. This aggregation is useful if the underlying maximization problem is a scalar @@ -135,14 +115,14 @@ def aggregator_identity(vector_model, fvec_center, model_info): """ intercept = float(vector_model.intercepts) linear_terms = np.squeeze(vector_model.linear_terms) - if model_info.has_squares or model_info.has_interactions: - square_terms = np.squeeze(vector_model.square_terms) - else: + if vector_model.square_terms is None: square_terms = np.zeros((len(linear_terms), len(linear_terms))) + else: + square_terms = np.squeeze(vector_model.square_terms) return intercept, linear_terms, square_terms -def aggregator_sum(vector_model, fvec_center, model_info): +def aggregator_sum(vector_model): """Aggregate quadratic VectorModel using sum function. This aggregation is useful if the underlying maximization problem is a likelihood @@ -163,7 +143,7 @@ def aggregator_sum(vector_model, fvec_center, model_info): return intercept, linear_terms, square_terms -def aggregator_least_squares_linear(vector_model, fvec_center, model_info): +def aggregator_least_squares_linear(vector_model): """Aggregate linear VectorModel assuming a least_squares functype. This aggregation is useful if the underlying maximization problem is a least-squares @@ -190,7 +170,7 @@ def aggregator_least_squares_linear(vector_model, fvec_center, model_info): return intercept, linear_terms, square_terms -def aggregator_information_equality_linear(vector_model, fvec_center, model_info): +def aggregator_information_equality_linear(vector_model): """Aggregate linear VectorModel using the Fisher information equality. This aggregation is useful if the underlying maximization problem is a likelihood diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index ca9dc79b7..397a7cbdc 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -221,7 +221,6 @@ def _tranquilo( scalar_model = aggregate_vector_model( vector_model=vector_model, - fvec_center=state.fvec, ) sub_sol = solve_subproblem(model=scalar_model, trustregion=trustregion) diff --git a/tests/optimization/tranquilo/test_aggregate_models.py b/tests/optimization/tranquilo/test_aggregate_models.py index 5ef0197b8..057ca9657 100644 --- a/tests/optimization/tranquilo/test_aggregate_models.py +++ b/tests/optimization/tranquilo/test_aggregate_models.py @@ -8,54 +8,45 @@ aggregator_least_squares_linear, ) from estimagic.optimization.tranquilo.aggregate_models import aggregator_sum -from estimagic.optimization.tranquilo.models import ModelInfo +from estimagic.optimization.tranquilo.models import ScalarModel from estimagic.optimization.tranquilo.models import VectorModel from numpy.testing import assert_array_equal -@pytest.mark.parametrize( - "model, expected_square_terms", - [ - ("quadratic", np.arange(9).reshape(3, 3)), - ("linear", np.zeros(9).reshape(3, 3)), - ], -) -def test_aggregator_identity(model, expected_square_terms): - if model == "linear": - model_info = ModelInfo(has_squares=False, has_interactions=False) - else: - model_info = ModelInfo() - - fvec_center = np.array([2.0]) - +@pytest.mark.parametrize("square_terms", [np.arange(9).reshape(1, 3, 3), None]) +def test_aggregator_identity(square_terms): vector_model = VectorModel( - intercepts=fvec_center, + intercepts=np.array([2.0]), linear_terms=np.arange(3).reshape(1, 3), - square_terms=np.arange(9).reshape(1, 3, 3), + square_terms=square_terms, ) - got = tuple(aggregator_identity(vector_model, fvec_center, model_info)) + if square_terms is None: + expected_square_terms = np.zeros((3, 3)) + else: + expected_square_terms = np.arange(9).reshape(3, 3) - assert got[0] == 2.0 - assert_array_equal(got[1], np.arange(3)) - assert_array_equal(got[2], expected_square_terms) + got = ScalarModel(*aggregator_identity(vector_model)) + assert_array_equal(got.intercept, 2.0) + assert_array_equal(got.linear_terms, np.arange(3)) + assert_array_equal(got.square_terms, expected_square_terms) -def test_aggregator_sum(): - model_info = ModelInfo() - fvec_center = np.array([1.0, 2.0]) +def test_aggregator_sum(): vector_model = VectorModel( - intercepts=fvec_center, + intercepts=np.array([1.0, 2.0]), linear_terms=np.arange(6).reshape(2, 3), square_terms=np.arange(18).reshape(2, 3, 3), ) - got = tuple(aggregator_sum(vector_model, fvec_center, model_info)) + got = ScalarModel(*aggregator_sum(vector_model)) - assert got[0] == 3.0 - assert_array_equal(got[1], np.array([3, 5, 7])) - assert_array_equal(got[2], np.array([[9, 11, 13], [15, 17, 19], [21, 23, 25]])) + assert_array_equal(got.intercept, 3.0) + assert_array_equal(got.linear_terms, np.array([3, 5, 7])) + assert_array_equal( + got.square_terms, np.array([[9, 11, 13], [15, 17, 19], [21, 23, 25]]) + ) def test_aggregator_least_squares_linear(): @@ -65,30 +56,27 @@ def test_aggregator_least_squares_linear(): square_terms=np.arange(18).reshape(2, 3, 3), # should not be used by aggregator ) - got = tuple(aggregator_least_squares_linear(vector_model, None, None)) + got = ScalarModel(*aggregator_least_squares_linear(vector_model)) - assert got[0] == 4.0 - assert_array_equal(got[1], np.array([12, 16, 20])) - assert_array_equal(got[2], np.array([[18, 24, 30], [24, 34, 44], [30, 44, 58]])) + assert_array_equal(got.intercept, 4.0) + assert_array_equal(got.linear_terms, np.array([12, 16, 20])) + assert_array_equal( + got.square_terms, np.array([[18, 24, 30], [24, 34, 44], [30, 44, 58]]) + ) def test_aggregator_information_equality_linear(): - model_info = ModelInfo() - fvec_center = np.array([1.0, 2.0]) - vector_model = VectorModel( - intercepts=fvec_center, + intercepts=np.array([1.0, 2.0]), linear_terms=np.arange(6).reshape(2, 3), square_terms=np.arange(18).reshape(2, 3, 3), # should not be used by aggregator ) - got = tuple( - aggregator_information_equality_linear(vector_model, fvec_center, model_info) - ) + got = ScalarModel(*aggregator_information_equality_linear(vector_model)) - assert got[0] == 3.0 - assert_array_equal(got[1], np.array([3, 5, 7])) + assert_array_equal(got.intercept, 3.0) + assert_array_equal(got.linear_terms, np.array([3, 5, 7])) assert_array_equal( - got[2], + got.square_terms, np.array([[-4.5, -6.0, -7.5], [-6.0, -8.5, -11.0], [-7.5, -11.0, -14.5]]), ) From 1e557d9d06ea616eceea111df5b78ebb0339819d Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 13:52:13 +0100 Subject: [PATCH 04/11] Improve construction of test cases in test_tranquilo.py --- .../optimization/tranquilo/test_tranquilo.py | 167 +++++++++--------- 1 file changed, 82 insertions(+), 85 deletions(-) diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index 25f33a2b2..b5dcac10e 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -1,4 +1,4 @@ -from itertools import product +import itertools import numpy as np import pytest @@ -8,63 +8,56 @@ from numpy.testing import assert_array_almost_equal as aaae +def _product(sample_filter, fitter, surrogate_model, sample_size): + # is used to create products of test cases + return list(itertools.product(sample_filter, fitter, surrogate_model, sample_size)) + + # ====================================================================================== # Scalar Tranquilo # ====================================================================================== -_sample_filter = ["discard_all", "keep_all"] -_fitter = ["ols"] -_surrogate_model = ["quadratic"] -_sample_size = ["quadratic"] -ols = list(product(_sample_filter, _fitter, _surrogate_model, _sample_size)) - -_sample_filter = ["keep_all"] -_fitter = ["ols"] -_surrogate_model = ["quadratic"] -_sample_size = ["powell"] -ols_keep_all = list(product(_sample_filter, _fitter, _surrogate_model, _sample_size)) - -_sample_filter = ["discard_all"] -_fitter = ["powell"] -_surrogate_model = ["quadratic"] -_sample_size = ["quadratic"] -pounders_discard_all = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) - -_sample_filter = ["keep_all"] -_fitter = ["powell"] -_surrogate_model = ["quadratic"] -_sample_size = ["linear", "powell", "quadratic"] -pounders_keep_all = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) - -_sample_filter = ["drop_pounders"] -_fitter = ["ols"] -_surrogate_model = ["quadratic"] -_sample_size = ["powell", "quadratic"] -ols_pounders_filtering = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) - -_sample_filter = ["drop_pounders"] -_fitter = ["powell"] -_surrogate_model = ["quadratic"] -_sample_size = ["linear", "powell", "quadratic"] -pounders_filtering = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) - - -TEST_CASES = ( - ols - + ols_keep_all - + pounders_discard_all - + pounders_keep_all - + ols_pounders_filtering - + pounders_filtering -) +TEST_CASES = { + "ols": { + "sample_filter": ["discard_all", "keep_all"], + "fitter": ["ols"], + "surrogate_model": ["quadratic"], + "sample_size": ["quadratic"], + }, + "ols_keep_all": { + "sample_filter": ["keep_all"], + "fitter": ["ols"], + "surrogate_model": ["quadratic"], + "sample_size": ["powell"], + }, + "pounders_discard_all": { + "sample_filter": ["discard_all"], + "fitter": ["powell"], + "surrogate_model": ["quadratic"], + "sample_size": ["quadratic"], + }, + "pounders_keep_all": { + "sample_filter": ["keep_all"], + "fitter": ["powell"], + "surrogate_model": ["quadratic"], + "sample_size": ["linear", "powell", "quadratic"], + }, + "ols_pounders_filtering": { + "sample_filter": ["drop_pounders"], + "fitter": ["ols"], + "surrogate_model": ["quadratic"], + "sample_size": ["powell", "quadratic"], + }, + "pounders_filtering": { + "sample_filter": ["drop_pounders"], + "fitter": ["powell"], + "surrogate_model": ["quadratic"], + "sample_size": ["linear", "powell", "quadratic"], + }, +} + +TEST_CASES = [_product(**kwargs) for kwargs in TEST_CASES.values()] +TEST_CASES = itertools.chain.from_iterable(TEST_CASES) @pytest.mark.parametrize( @@ -88,25 +81,27 @@ def test_internal_tranquilo_scalar_sphere_defaults( # Imprecise options for scalar tranquilo # ====================================================================================== -_sample_filter = ["keep_all"] -_fitter = ["ols"] -_surrogate_model = ["quadratic"] -_sample_size = ["linear"] -ols_keep_all = list(product(_sample_filter, _fitter, _surrogate_model, _sample_size)) - -_sample_filter = ["discard_all"] -_fitter = ["powell"] -_surrogate_model = ["quadratic"] -_sample_size = ["powell"] -pounders_discard_all = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) - -TEST_CASES_IMPRECISE = ols_keep_all + pounders_discard_all +TEST_CASES = { + "ls_keep": { + "sample_filter": ["keep_all"], + "fitter": ["ols"], + "surrogate_model": ["quadratic"], + "sample_size": ["linear"], + }, + "pounders_discard_all": { + "sample_filter": ["discard_all"], + "fitter": ["powell"], + "surrogate_model": ["quadratic"], + "sample_size": ["powell"], + }, +} + +TEST_CASES = [_product(**kwargs) for kwargs in TEST_CASES.values()] +TEST_CASES = itertools.chain.from_iterable(TEST_CASES) @pytest.mark.parametrize( - "sample_filter, fitter, surrogate_model, sample_size", TEST_CASES_IMPRECISE + "sample_filter, fitter, surrogate_model, sample_size", TEST_CASES ) def test_internal_tranquilo_scalar_sphere_imprecise_defaults( sample_filter, fitter, surrogate_model, sample_size @@ -141,22 +136,24 @@ def test_external_tranquilo_scalar_sphere_defaults(): # Least-squares Tranquilo # ====================================================================================== -_sample_filter = ["keep_all", "discard_all"] -_fitter = ["ols"] -_surrogate_model = ["linear"] -_sample_size = ["linear"] -ols = list(product(_sample_filter, _fitter, _surrogate_model, _sample_size)) - - -_sample_filter = ["drop_pounders"] -_fitter = ["ols"] -_surrogate_model = ["linear"] -_sample_size = ["linear"] -pounders_filtering = list( - product(_sample_filter, _fitter, _surrogate_model, _sample_size) -) -TEST_CASES = ols + pounders_filtering +TEST_CASES = { + "ols": { + "sample_filter": ["keep_all", "discard_all"], + "fitter": ["ols"], + "surrogate_model": ["linear"], + "sample_size": ["linear"], + }, + "pounders_filtering": { + "sample_filter": ["drop_pounders"], + "fitter": ["ols"], + "surrogate_model": ["linear"], + "sample_size": ["linear"], + }, +} + +TEST_CASES = [_product(**kwargs) for kwargs in TEST_CASES.values()] +TEST_CASES = itertools.chain.from_iterable(TEST_CASES) @pytest.mark.parametrize( From 63d505267e8bda8c70297ec7b66cb8fe77da3a85 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 14:57:08 +0100 Subject: [PATCH 05/11] Test input processing functions in tranquilo.py --- .../optimization/tranquilo/test_tranquilo.py | 113 ++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/tests/optimization/tranquilo/test_tranquilo.py b/tests/optimization/tranquilo/test_tranquilo.py index b5dcac10e..c82a3cf4b 100644 --- a/tests/optimization/tranquilo/test_tranquilo.py +++ b/tests/optimization/tranquilo/test_tranquilo.py @@ -3,11 +3,19 @@ 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 +from estimagic.optimization.tranquilo.tranquilo import _process_surrogate_model from estimagic.optimization.tranquilo.tranquilo import tranquilo from estimagic.optimization.tranquilo.tranquilo import tranquilo_ls from numpy.testing import assert_array_almost_equal as aaae +# ====================================================================================== +# Test tranquilo end-to-end +# ====================================================================================== + + def _product(sample_filter, fitter, surrogate_model, sample_size): # is used to create products of test cases return list(itertools.product(sample_filter, fitter, surrogate_model, sample_size)) @@ -186,3 +194,108 @@ def test_external_tranquilo_ls_sphere_defaults(): ) aaae(res.params, np.zeros(5), decimal=5) + + +# ====================================================================================== +# Test input processing functions +# ====================================================================================== + + +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 + + +@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 + + +@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 + + +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 + + +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 + + +def test_process_surrogate_model_str_invalid(): + with pytest.raises(ValueError): + _process_surrogate_model("whatever", None) + + +@pytest.mark.parametrize("functype", ["scalar", "least_squares"]) +def test_process_surrogate_model_invalid(functype): + surrogate_model = np.linalg.lstsq + with pytest.raises(ValueError): + _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) + x = np.ones((3, 2)) + got = _process_sample_size(None, model_info=model_info, x=x) + if has_interactions or has_squares: + assert got == 7 + else: + assert got == 4 + + +STR_TEST_CASES = [ # assume len(x) = 3 + # (user_sample_size, expected) # noqa: E800 + ("linear", 3 + 1), + ("n+1", 3 + 1), + ("n + 1", 3 + 1), + ("powell", 2 * 3 + 1), + ("2n+1", 2 * 3 + 1), + ("2 n + 1", 2 * 3 + 1), + ("2*n + 1", 2 * 3 + 1), + ("quadratic", 6 + 3 + 1), +] + + +@pytest.mark.parametrize("user_sample_size, expected", STR_TEST_CASES) +def test_process_sample_size_str(user_sample_size, expected): + x = np.ones((3, 2)) + got = _process_sample_size(user_sample_size, None, x=x) + assert got == expected + + +def test_process_sample_size_str_invalid(): + with pytest.raises(ValueError): + _process_sample_size("n**2", None, None) + + +@pytest.mark.parametrize("user_sample_size", [1, 10, -100, 10.5]) +def test_process_sample_size_number(user_sample_size): + got = _process_sample_size(user_sample_size, None, None) + assert got == int(user_sample_size) + + +def test_process_sample_size_invalid(): + x = np.ones((3, 2)) + with pytest.raises(ValueError): + _process_sample_size(np.zeros_like(x), None, x) From b677a4291adf8e5c0782d2f65877eb832ccb4588 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 15:40:54 +0100 Subject: [PATCH 06/11] Refactoring --- src/estimagic/optimization/tranquilo/fit_models.py | 6 ++++-- src/estimagic/optimization/tranquilo/tranquilo.py | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/fit_models.py b/src/estimagic/optimization/tranquilo/fit_models.py index 3d46a99e5..48d872be2 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -4,6 +4,8 @@ import numpy as np from estimagic.optimization.tranquilo.models import ModelInfo +from estimagic.optimization.tranquilo.models import n_interactions +from estimagic.optimization.tranquilo.models import n_second_order_terms from estimagic.optimization.tranquilo.models import VectorModel from numba import njit from scipy.linalg import qr_multiply @@ -442,9 +444,9 @@ def _polynomial_features(x, has_squares): n_samples, n_params = x.shape if has_squares: - n_poly_terms = n_params * (n_params + 1) // 2 + n_poly_terms = n_second_order_terms(n_params) else: - n_poly_terms = n_params * (n_params - 1) // 2 + n_poly_terms = n_interactions(n_params) poly_terms = np.empty((n_poly_terms, n_samples), x.dtype) xt = x.T diff --git a/src/estimagic/optimization/tranquilo/tranquilo.py b/src/estimagic/optimization/tranquilo/tranquilo.py index 397a7cbdc..9c9b6e3cf 100644 --- a/src/estimagic/optimization/tranquilo/tranquilo.py +++ b/src/estimagic/optimization/tranquilo/tranquilo.py @@ -348,10 +348,10 @@ def _process_surrogate_model(surrogate_model, functype): elif isinstance(surrogate_model, str): if surrogate_model == "linear": out = ModelInfo(has_squares=False, has_interactions=False) - elif surrogate_model == "quadratic": - out = ModelInfo(has_squares=True, has_interactions=True) 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}") From a0f454b2258b840fb485283308205d8f8bb27eb6 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 16:13:54 +0100 Subject: [PATCH 07/11] Add tests for _scaled_square_features --- .../optimization/tranquilo/filter_points.py | 27 ++++++++++++++----- .../tranquilo/test_filter_points.py | 8 ++++++ 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/filter_points.py b/src/estimagic/optimization/tranquilo/filter_points.py index 921c6e27a..3a31ea661 100644 --- a/src/estimagic/optimization/tranquilo/filter_points.py +++ b/src/estimagic/optimization/tranquilo/filter_points.py @@ -98,7 +98,9 @@ def _drop_collinear_pounders(xs, indices, state): continue linear_features[counter, 1:] = centered_xs[index] - square_features[counter, :] = _square_features(linear_features[counter, 1:]) + 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 @@ -128,7 +130,7 @@ def _get_polynomial_feature_matrices( square_features = np.zeros((n_samples, n_poly_terms)) linear_features[0, 1:] = centered_xs[indexer[index_center]] - square_features[0, :] = _square_features(linear_features[0, 1:]).flatten() + 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] @@ -136,7 +138,7 @@ def _get_polynomial_feature_matrices( linear_features[:, 0] = 1 linear_features[: n_params + 1, 1:] = centered_xs[indexer[idx_list_n_plus_1]] - square_features[: n_params + 1, :] = _square_features( + square_features[: n_params + 1, :] = _scaled_square_features( linear_features[: n_params + 1, 1:] ) @@ -146,16 +148,29 @@ def _get_polynomial_feature_matrices( @njit -def _square_features(x): +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), x.dtype) + poly_terms = np.empty((n_poly_terms, n_samples), float) xt = x.T idx = 0 for i in range(n_params): - poly_terms[idx] = 0.5 * xt[i] ** 2 + poly_terms[idx] = xt[i] ** 2 / 2 idx += 1 for j in range(i + 1, n_params): diff --git a/tests/optimization/tranquilo/test_filter_points.py b/tests/optimization/tranquilo/test_filter_points.py index 6407e8d15..f99ebb56b 100644 --- a/tests/optimization/tranquilo/test_filter_points.py +++ b/tests/optimization/tranquilo/test_filter_points.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from estimagic.optimization.tranquilo.filter_points import _scaled_square_features from estimagic.optimization.tranquilo.filter_points import drop_collinear_pounders from estimagic.optimization.tranquilo.options import TrustRegion from estimagic.optimization.tranquilo.tranquilo import State @@ -208,3 +209,10 @@ def test_drop_collinear_pounders(test_case, request): assert_equal(filtered_indices, expected_indices) aaae(filtered_xs, expected_xs) + + +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) From ccc04a686a122defec76106ddd7930d73ce8a907 Mon Sep 17 00:00:00 2001 From: Tim Date: Wed, 2 Nov 2022 16:41:55 +0100 Subject: [PATCH 08/11] Change dtype to float in _polynomial_features --- 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 48d872be2..89c0bc368 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -448,7 +448,7 @@ def _polynomial_features(x, has_squares): else: n_poly_terms = n_interactions(n_params) - poly_terms = np.empty((n_poly_terms, n_samples), x.dtype) + poly_terms = np.empty((n_poly_terms, n_samples), float) xt = x.T idx = 0 From 30077834dafe9e90aa8bf59f32051d171708722d Mon Sep 17 00:00:00 2001 From: Tim Date: Thu, 3 Nov 2022 17:08:59 +0100 Subject: [PATCH 09/11] Replace solve by lstsq in fit_models --- src/estimagic/optimization/tranquilo/filter_points.py | 2 +- src/estimagic/optimization/tranquilo/fit_models.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/estimagic/optimization/tranquilo/filter_points.py b/src/estimagic/optimization/tranquilo/filter_points.py index 3a31ea661..115662a9f 100644 --- a/src/estimagic/optimization/tranquilo/filter_points.py +++ b/src/estimagic/optimization/tranquilo/filter_points.py @@ -165,7 +165,7 @@ def _scaled_square_features(x): 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), float) + poly_terms = np.empty((n_poly_terms, n_samples), np.float64) xt = x.T idx = 0 diff --git a/src/estimagic/optimization/tranquilo/fit_models.py b/src/estimagic/optimization/tranquilo/fit_models.py index 89c0bc368..861ccd826 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -371,16 +371,15 @@ def _get_current_fit_minimal_frobenius_norm_of_hessian( for k in range(n_residuals): z_y_vec = np.dot(z_mat.T, y[:, k]) - coeffs_first_stage = np.linalg.solve( - np.atleast_2d(n_z_mat_square), - np.atleast_1d(z_y_vec), + coeffs_first_stage, *_ = np.linalg.lstsq( + np.atleast_2d(n_z_mat_square), np.atleast_1d(z_y_vec), rcond=None ) coeffs_second_stage = np.atleast_2d(n_z_mat) @ coeffs_first_stage rhs = y[:, k] - n_mat @ coeffs_second_stage - alpha = np.linalg.solve(m_mat, rhs[: n_params + 1]) + alpha, *_ = np.linalg.lstsq(m_mat, rhs[: n_params + 1], rcond=None) coeffs_linear[k, :] = alpha[offset : (n_params + 1)] coeffs_square[k] = coeffs_second_stage From 7ed4c76dc6ae1d3604fc37d977d26fe6e82d8eff Mon Sep 17 00:00:00 2001 From: Sebastian Gsell Date: Fri, 4 Nov 2022 12:10:30 +0100 Subject: [PATCH 10/11] Change dtype to np.float64 in _polynomial_features. --- 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 861ccd826..6c487c492 100644 --- a/src/estimagic/optimization/tranquilo/fit_models.py +++ b/src/estimagic/optimization/tranquilo/fit_models.py @@ -447,7 +447,7 @@ def _polynomial_features(x, has_squares): else: n_poly_terms = n_interactions(n_params) - poly_terms = np.empty((n_poly_terms, n_samples), float) + poly_terms = np.empty((n_poly_terms, n_samples), np.float64) xt = x.T idx = 0 From e00d3026685998277327d40e8c2e3e0fb7f2bb3c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 4 Nov 2022 14:51:24 +0100 Subject: [PATCH 11/11] [pre-commit.ci] pre-commit autoupdate (#394) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ad3e35f11..3806ee0bc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,7 +6,7 @@ repos: - id: debug-statements - id: end-of-file-fixer - repo: https://github.com/asottile/reorder_python_imports - rev: v3.8.5 + rev: v3.9.0 hooks: - id: reorder-python-imports types: [python]