diff --git a/src/tranquilo/acceptance_decision.py b/src/tranquilo/acceptance_decision.py index 56e5e69..c5c285c 100644 --- a/src/tranquilo/acceptance_decision.py +++ b/src/tranquilo/acceptance_decision.py @@ -7,6 +7,7 @@ from typing import NamedTuple import numpy as np +import pandas as pd from tranquilo.acceptance_sample_size import ( get_acceptance_sample_sizes, @@ -20,6 +21,7 @@ def get_acceptance_decider(acceptance_decider, acceptance_options): "classic": _accept_classic, "naive_noisy": accept_naive_noisy, "noisy": accept_noisy, + "classic_line_search": accept_classic_line_search, } out = get_component( @@ -85,6 +87,142 @@ def accept_naive_noisy( return out +def accept_classic_line_search( + subproblem_solution, + state, + history, + *, + wrapped_criterion, + min_improvement, + batch_size, + sample_points, + search_radius_factor, + rng, +): + # ================================================================================== + # Quick return if batch_size is 1 + + if batch_size == 1: + return _accept_classic( + subproblem_solution=subproblem_solution, + state=state, + history=history, + wrapped_criterion=wrapped_criterion, + min_improvement=min_improvement, + ) + + # ================================================================================== + # Add candidate to history + + candidate_x = subproblem_solution.x + candidate_index = history.add_xs(candidate_x) + + eval_info = {candidate_index: 1} + + # ================================================================================== + # Determine whether the candidate it sufficiently close to the border of the + # trustregion, in which case we perform a line search + + perform_line_search = _is_on_border(state.trustregion, x=candidate_x, rtol=1e-1) + + if perform_line_search: + alpha_grid = _generate_alpha_grid(batch_size) + + line_search_xs = _sample_on_line( + start_point=state.x, direction_point=candidate_x, alpha_grid=alpha_grid + ) + else: + line_search_xs = None + + # ================================================================================== + # Check whether there are any unallocated evaluations left, and if yes perform a + # speculative sampling + + n_evals_line_search = 0 if line_search_xs is None else len(line_search_xs) + n_unallocated_evals = batch_size - 1 - n_evals_line_search + + if n_unallocated_evals > 0: + speculative_xs = _generate_speculative_sample( + new_center=candidate_x, + search_radius_factor=search_radius_factor, + trustregion=state.trustregion, + sample_points=sample_points, + n_points=n_unallocated_evals, + history=history, + rng=rng, + ) + else: + speculative_xs = None + + # ================================================================================== + # Consolidate newly sampled points + + if line_search_xs is not None and speculative_xs is not None: + new_xs = np.vstack([line_search_xs, speculative_xs]) + elif line_search_xs is not None: + new_xs = line_search_xs + elif speculative_xs is not None: + new_xs = speculative_xs + + # ================================================================================== + # Add new points to history and evaluate criterion + + new_indices = history.add_xs(new_xs) + + for idx in new_indices: + eval_info[idx] = 1 + + wrapped_criterion(eval_info) + + # ================================================================================== + # Calculate rho + + candidate_fval = np.mean(history.get_fvals(candidate_index)) + + actual_improvement = -(candidate_fval - state.fval) + + rho = calculate_rho( + actual_improvement=actual_improvement, + expected_improvement=subproblem_solution.expected_improvement, + ) + + # ================================================================================== + # Check if there are any better points + + new_fvals = history.get_fvals(new_indices) + new_fvals = pd.Series({i: np.mean(fvals) for i, fvals in new_fvals.items()}) + new_fval_argmin = new_fvals.idxmin() + + found_better_candidate = new_fvals.loc[new_fval_argmin] < candidate_fval + + # If a better point was found, update the candidates + if found_better_candidate: + candidate_x = history.get_xs(new_fval_argmin) + candidate_fval = new_fvals.loc[new_fval_argmin] + candidate_index = new_fval_argmin + + # ================================================================================== + # Calculate the overall improvement using a potentially updated candidate and draw + # the acceptance conclusions based on that. + + overall_improvement = -(candidate_fval - state.fval) + is_accepted = overall_improvement >= min_improvement + + # ================================================================================== + # Return results + + res = _get_acceptance_result( + candidate_x=candidate_x, + candidate_fval=candidate_fval, + candidate_index=candidate_index, + rho=rho, + is_accepted=is_accepted, + old_state=state, + n_evals=1, + ) + return res + + def _accept_simple( subproblem_solution, state, @@ -247,3 +385,100 @@ def calculate_rho(actual_improvement, expected_improvement): else: rho = actual_improvement / expected_improvement return rho + + +# ====================================================================================== +# Helper functions for line search +# ====================================================================================== + + +def _generate_speculative_sample( + new_center, trustregion, sample_points, n_points, history, search_radius_factor, rng +): + """Generative a speculative sample. + + Args: + new_center (np.ndarray): New center of the trust region. + trustregion (Region): Current trust region. + sample_points (callable): Function to sample points. + n_points (int): Number of points to sample. + history (History): Tranquilo history. + search_radius_factor (float): Factor to multiply the trust region radius by to + get the search radius. + rng (np.random.Generator): Random number generator. + + Returns: + np.ndarray: Speculative sample. + + """ + search_region = trustregion._replace( + center=new_center, radius=search_radius_factor * trustregion.radius + ) + + old_indices = history.get_x_indices_in_region(search_region) + + old_xs = history.get_xs(old_indices) + + model_xs = old_xs + + new_xs = sample_points( + search_region, + n_points=n_points, + existing_xs=model_xs, + rng=rng, + ) + return new_xs + + +def _sample_on_line(start_point, direction_point, alpha_grid): + """Sample points on a line defined by startind and direction points. + + Args: + start_point (np.ndarray): Starting point of the line. + direction_point (np.ndarray): Direction point of the line. + alpha_grid (np.ndarray): Grid of alphas to sample points on the line. 0 + corresponds to the starting point and 1 corresponds to the direction point. + Points larger than 1 are beyond the direction point. + + Returns: + np.ndarray: Sampled points on the line. + + """ + xs = start_point + alpha_grid.reshape(-1, 1) * (direction_point - start_point) + return xs + + +def _is_on_border(trustregion, x, rtol): + """Check whether a point is sufficiently close to the border of a trust region. + + Args: + trustregion (Region): Trust region. + x (np.ndarray): Point to check. + rtol (float): Relative tolerance. + + Returns: + bool: True if the point is sufficiently close to the border of the trust region. + + """ + if trustregion.shape == "sphere": + candidate_on_border = _is_on_sphere_border(trustregion, x=x, rtol=rtol) + else: + candidate_on_border = _is_on_cube_border(trustregion, x=x, rtol=rtol) + return candidate_on_border + + +def _is_on_sphere_border(trustregion, x, rtol): + x_center_dist = np.linalg.norm(x - trustregion.center, ord=2) + return np.isclose(x_center_dist, trustregion.radius, rtol=rtol) + + +def _is_on_cube_border(trustregion, x, rtol): + cube_bounds = trustregion.cube_bounds + is_on_lower_border = np.isclose(x, cube_bounds.lower, rtol=rtol).any() + is_on_upper_border = np.isclose(x, cube_bounds.upper, rtol=rtol).any() + return is_on_lower_border or is_on_upper_border + + +def _generate_alpha_grid(batch_size): + n_points = min(batch_size, 4) - 1 + return 2 ** np.arange(1, n_points + 1, dtype=float) diff --git a/src/tranquilo/options.py b/src/tranquilo/options.py index f19d6b4..72bad1a 100644 --- a/src/tranquilo/options.py +++ b/src/tranquilo/options.py @@ -4,6 +4,25 @@ import numpy as np +def get_default_stagnation_options(noisy, batch_size): + if noisy: + out = StagnationOptions( + min_relative_step_keep=0.0, + drop=False, + ) + elif batch_size > 1: + out = StagnationOptions( + min_relative_step_keep=0.0, + drop=True, + ) + else: + out = StagnationOptions( + min_relative_step_keep=0.125, + drop=True, + ) + return out + + def get_default_radius_options(x): return RadiusOptions(initial_radius=0.1 * np.max(np.abs(x))) @@ -79,20 +98,6 @@ def get_default_n_evals_per_point(noisy, noise_adaptation_options): return noise_adaptation_options.min_n_evals if noisy else 1 -def get_default_stagnation_options(noisy): - if noisy: - out = StagnationOptions( - min_relative_step_keep=0.0, - drop=False, - ) - else: - out = StagnationOptions( - min_relative_step_keep=0.125, - drop=True, - ) - return out - - class StopOptions(NamedTuple): """Criteria for stopping without successful convergence.""" diff --git a/src/tranquilo/process_arguments.py b/src/tranquilo/process_arguments.py index 3c129c9..d0c0524 100644 --- a/src/tranquilo/process_arguments.py +++ b/src/tranquilo/process_arguments.py @@ -13,7 +13,6 @@ from tranquilo.history import History from tranquilo.options import ( ConvOptions, - get_default_stagnation_options, StopOptions, get_default_acceptance_decider, get_default_aggregator, @@ -26,6 +25,7 @@ get_default_radius_options, get_default_sample_size, get_default_search_radius_factor, + get_default_stagnation_options, update_option_bundle, NoiseAdaptationOptions, ) @@ -113,10 +113,6 @@ def process_arguments( x = _process_x(x) noisy = _process_noisy(noisy) n_cores = _process_n_cores(n_cores) - stagnation_options = update_option_bundle( - get_default_stagnation_options(noisy), - stagnation_options, - ) sampling_rng = _process_seed(seed) simulation_rng = _process_seed(seed + 1000) @@ -142,6 +138,9 @@ def process_arguments( acceptance_decider = _process_acceptance_decider(acceptance_decider, noisy) # process options that depend on arguments with dependent defaults + stagnation_options = update_option_bundle( + get_default_stagnation_options(noisy, batch_size=batch_size), stagnation_options + ) target_sample_size = _process_sample_size( sample_size=sample_size, model_type=model_type, @@ -332,6 +331,10 @@ def _process_n_evals_at_start(n_evals, noisy): return out +def next_multiple(n, base): + return int(np.ceil(n / base) * base) + + def _process_n_evals_per_point(n_evals, noisy, noise_adaptation_options): if n_evals is None: out = get_default_n_evals_per_point(noisy, noise_adaptation_options) diff --git a/src/tranquilo/tranquilo.py b/src/tranquilo/tranquilo.py index 2203676..ca1ee5e 100644 --- a/src/tranquilo/tranquilo.py +++ b/src/tranquilo/tranquilo.py @@ -13,7 +13,7 @@ ScalarModel, VectorModel, ) -from tranquilo.process_arguments import process_arguments +from tranquilo.process_arguments import process_arguments, next_multiple from tranquilo.region import Region from tranquilo.rho_noise import simulate_rho_noise from tranquilo.adjust_n_evals import adjust_n_evals @@ -52,7 +52,10 @@ def _internal_tranquilo( estimate_variance, accept_candidate, ): - eval_info = {0: n_evals_at_start} + if n_evals_at_start > 1: + eval_info = {0: next_multiple(n_evals_at_start, base=batch_size)} + else: + eval_info = {0: 1} evaluate_criterion(eval_info) @@ -125,9 +128,13 @@ def _internal_tranquilo( # ========================================================================== # sample points if necessary and do simple iteration # ========================================================================== + + n_new_points = max(0, target_sample_size - len(model_xs)) + n_new_points = next_multiple(n_new_points, base=batch_size) + new_xs = sample_points( trustregion=state.trustregion, - n_points=max(0, target_sample_size - len(model_xs)), + n_points=n_new_points, existing_xs=model_xs, rng=sampling_rng, ) @@ -211,17 +218,19 @@ def _internal_tranquilo( sample_counter = 0 while _relative_step_length < stagnation_options.min_relative_step: - if stagnation_options.drop: + n_to_drop = stagnation_options.sample_increment + + if stagnation_options.drop and len(model_xs) > n_to_drop: model_xs, model_indices = drop_worst_points( xs=model_xs, indices=model_indices, state=state, - n_to_drop=stagnation_options.sample_increment, + n_to_drop=n_to_drop, ) new_xs = sample_points( trustregion=state.trustregion, - n_points=stagnation_options.sample_increment, + n_points=n_to_drop, existing_xs=model_xs, rng=sampling_rng, ) @@ -291,6 +300,10 @@ def _internal_tranquilo( wrapped_criterion=evaluate_criterion, noise_variance=noise_variance, history=history, + search_radius_factor=search_radius_factor, + batch_size=batch_size, + sample_points=sample_points, + rng=sampling_rng, ) # ============================================================================== diff --git a/src/tranquilo/wrap_criterion.py b/src/tranquilo/wrap_criterion.py index 8485713..1a74962 100644 --- a/src/tranquilo/wrap_criterion.py +++ b/src/tranquilo/wrap_criterion.py @@ -8,12 +8,24 @@ def get_wrapped_criterion(criterion, batch_evaluator, n_cores, history): """Wrap the criterion function to do get parallelization and history handling. + Notes + ----- + The wrapped criterion function takes a dict mapping x_indices to required numbers of evaluations as only argument. It evaluates the criterion function in parallel and saves the resulting function evaluations in the history. The wrapped criterion function does not return anything. + Args: + criterion (function): The criterion function to wrap. + batch_evaluator (function): The batch evaluator to use. + n_cores (int): The number of cores to use. + history (History): The tranquilo history. + + Returns: + callable: The wrapped criterion function. + """ batch_evaluator = process_batch_evaluator(batch_evaluator) diff --git a/tests/test_acceptance_decision.py b/tests/test_acceptance_decision.py index 37a9f0e..062e772 100644 --- a/tests/test_acceptance_decision.py +++ b/tests/test_acceptance_decision.py @@ -2,13 +2,21 @@ import numpy as np import pytest +from tranquilo.sample_points import get_sampler from tranquilo.acceptance_decision import ( _accept_simple, _get_acceptance_result, calculate_rho, + _generate_alpha_grid, + _is_on_border, + _is_on_cube_border, + _is_on_sphere_border, + _sample_on_line, + _generate_speculative_sample, ) from tranquilo.history import History from tranquilo.region import Region +from tranquilo.bounds import Bounds from tranquilo.solve_subproblem import SubproblemResult from numpy.testing import assert_array_equal @@ -31,7 +39,7 @@ def subproblem_solution(): # ====================================================================================== -# Test accept_xxx +# Test accept_simple # ====================================================================================== @@ -140,3 +148,118 @@ def test_get_acceptance_result(): def test_calculate_rho(actual_improvement, expected_improvement, expected): rho = calculate_rho(actual_improvement, expected_improvement) assert rho == expected + + +# ====================================================================================== +# Test _generate_alpha_grid +# ====================================================================================== + +CASES = zip( + [1, 2, 4, 6], + [np.array([]), np.array([2]), np.array([2, 4, 8]), np.array([2, 4, 8])], +) + + +@pytest.mark.parametrize("batch_size, expected", CASES) +def test_generate_alpha_grid(batch_size, expected): + alpha_grid = _generate_alpha_grid(batch_size) + assert_array_equal(alpha_grid, expected) + + +# ====================================================================================== +# Test generating speculative sample +# ====================================================================================== + + +def test_generate_speculative_sample(): + trustregion = Region(center=np.zeros(2), radius=1.0) + + history = namedtuple("History", "get_x_indices_in_region, get_xs") + + history.get_x_indices_in_region = lambda _: None + history.get_xs = lambda _: np.atleast_2d(np.ones(2)) + + got = _generate_speculative_sample( + new_center=np.ones(2), + trustregion=trustregion, + sample_points=get_sampler("random_hull"), + n_points=3, + history=history, + search_radius_factor=1.0, + rng=np.random.default_rng(1234), + ) + + assert len(got) == 3 + for point in got: + assert _is_on_sphere_border( + trustregion._replace(center=np.ones(2)), point, rtol=0 + ) + + +# ====================================================================================== +# Test sampling on line +# ====================================================================================== + + +def test_sample_on_line(): + start_point = np.zeros(2) + direction_point = np.array([1, 2]) + alpha_grid = np.array([0, 1, 2, 3.5]) + got = _sample_on_line(start_point, direction_point, alpha_grid) + exp = np.array([[0, 0], [1, 2], [2, 4], [3.5, 7]]) + assert_array_equal(got, exp) + + +def test_sample_on_line_no_points(): + start_point = np.zeros(2) + direction_point = np.array([1, 2]) + alpha_grid = np.array([]) + got = _sample_on_line(start_point, direction_point, alpha_grid) + exp = np.array([]).reshape(0, 2) + assert_array_equal(got, exp) + + +# ====================================================================================== +# Test border check functions +# ====================================================================================== + +CASES = [ + (np.array([0, 0]), 1, True), + (np.array([0, 0]), 0.5, False), + (np.array([0, 1]), 0.0, True), + (np.array([0, 0.9]), 0.1, True), + (np.array([0, 0.9]), 0.09, False), +] + + +@pytest.mark.parametrize("x, rtol, expected", CASES) +def test_is_on_sphere_border(x, rtol, expected): + region = Region(center=np.zeros(2), radius=1.0) + assert _is_on_sphere_border(region, x, rtol) == expected + + +CASES = [ + (np.ones(2), 0, True), + (0.9 * np.ones(2), 0.1, True), + (0.8 * np.ones(2), 0.1, False), +] + + +@pytest.mark.parametrize("x, rtol, expected", CASES) +def test_is_on_cube_border(x, rtol, expected): + bounds = Bounds(lower=-np.ones(2), upper=np.ones(2)) + region = Region(center=np.zeros(2), radius=2.0, bounds=bounds) + assert _is_on_cube_border(region, x, rtol) == expected + + +def test_is_on_border_sphere(): + region = Region(center=np.zeros(2), radius=1.0) + assert _is_on_border(region, np.array([0, 0.9]), 0.1) + assert not _is_on_border(region, np.array([0, 0.9]), 0.09) + + +def test_is_on_border_cube(): + bounds = Bounds(lower=-np.ones(2), upper=np.ones(2)) + region = Region(center=np.zeros(2), radius=2.0, bounds=bounds) + assert _is_on_border(region, 0.9 * np.ones(2), 0.1) + assert not _is_on_border(region, 0.8 * np.ones(2), 0.1) diff --git a/tests/test_process_arguments.py b/tests/test_process_arguments.py index cd93edb..4e4af10 100644 --- a/tests/test_process_arguments.py +++ b/tests/test_process_arguments.py @@ -15,7 +15,7 @@ _process_acceptance_decider, _process_model_fitter, _process_residualize, - _process_n_evals_at_start, + next_multiple, ) @@ -126,12 +126,9 @@ def test_process_residualize_invalid(): _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 - assert _process_n_evals_at_start(n_evals=10, noisy=None) == 10 - - -def test_process_n_evals_at_start_negative(): - with pytest.raises(ValueError, match="n_initial_acceptance_evals must be"): - _process_n_evals_at_start(n_evals=-1, noisy=None) +def test_roundup_to_next_multiple_of_batch_size(): + assert next_multiple(1, 1) == 1 + assert next_multiple(456, 456) == 456 + assert next_multiple(123123, 1) == 123123 + assert next_multiple(1, 123123) == 123123 + assert next_multiple(4, 10) == 10 diff --git a/tests/test_tranquilo.py b/tests/test_tranquilo.py index 084534c..45f0a6b 100644 --- a/tests/test_tranquilo.py +++ b/tests/test_tranquilo.py @@ -176,7 +176,7 @@ def test_external_tranquilo_ls_sphere_defaults(): # ====================================================================================== -@pytest.mark.parametrize("algo", ["tranquilo", "tranquilo_ls"]) +@pytest.mark.parametrize("algo", [tranquilo, tranquilo_ls]) def test_tranquilo_with_noise_handling_and_deterministic_function(algo): def _f(x): return {"root_contributions": x, "value": x @ x} @@ -184,8 +184,8 @@ def _f(x): res = minimize( criterion=_f, params=np.arange(5), + algorithm=algo, algo_options={"noisy": True}, - algorithm=tranquilo, ) aaae(res.params, np.zeros(5), decimal=3)