Adaptive noise handling (#5)

src/tranquilo/acceptance_decision.py

@@ -80,7 +80,7 @@ def accept_naive_noisy(
         history=history,
         wrapped_criterion=wrapped_criterion,
         min_improvement=min_improvement,
-        n_evals=5,
+        n_evals=10,
     )
     return out
 
@@ -131,6 +131,7 @@ def _accept_simple(
         rho=rho,
         is_accepted=is_accepted,
         old_state=state,
+        n_evals=n_evals,
     )
 
     return res
@@ -189,6 +190,7 @@ def accept_noisy(
         rho=rho,
         is_accepted=is_accepted,
         old_state=state,
+        n_evals=n_2,
     )
 
     return res
@@ -204,6 +206,7 @@ class AcceptanceResult(NamedTuple):
     relative_step_length: float
     candidate_index: int
     candidate_x: np.ndarray
+    n_evals_acceptance: int
 
 
 def _get_acceptance_result(
@@ -213,6 +216,7 @@ def _get_acceptance_result(
     rho,
     is_accepted,
     old_state,
+    n_evals,
 ):
     x = candidate_x if is_accepted else old_state.x
     fval = candidate_fval if is_accepted else old_state.fval
@@ -230,6 +234,7 @@ def _get_acceptance_result(
         relative_step_length=relative_step_length,
         candidate_index=candidate_index,
         candidate_x=candidate_x,
+        n_evals_acceptance=n_evals,
     )
     return out
 

src/tranquilo/adjust_n_evals.py

@@ -0,0 +1,29 @@
+def adjust_n_evals(n_evals, rho, rho_noise, options):
+    """Adjust the number of evaluations based on the noise adaptation options.
+
+    Args:
+        n_evals (int): The current number of evaluations.
+        rho_noise (np.ndarray): The simulated rho_noises.
+        options (NoiseAdaptationOptions): Options for noise adaptation.
+
+    Returns:
+        int: The updated number of evaluations.
+        bool: Whether the number of evaluations was increased.
+
+    """
+    # most rhos are very high -> decrease
+    if (rho_noise > options.high_rho).mean() > options.min_share_high_rho:
+        new_n_evals = max(n_evals - 1, options.min_n_evals)
+
+    # most rhos are above rho low -> keep constant
+    elif (
+        rho_noise > options.low_rho
+    ).mean() > options.min_share_low_rho or rho >= options.good_rho_threshold:
+        new_n_evals = n_evals
+
+    # many rhos are below rho low -> increase
+    else:
+        new_n_evals = min(n_evals + 1, options.max_n_evals)
+
+    is_increased = new_n_evals > n_evals
+    return new_n_evals, is_increased

src/tranquilo/adjust_radius.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 
-def adjust_radius(radius, rho, step_length, options):
+def adjust_radius(radius, rho, step_length, options, n_evals_is_increased):
     """Adjust the trustregion radius based on relative improvement and stepsize.
 
     This is just a slight generalization of the pounders radius adjustment. With default
@@ -16,6 +16,7 @@ def adjust_radius(radius, rho, step_length, options):
             parameter vectors.
         step (np.ndarray): The step between the last two accepted parameter vectors.
         options (NamedTuple): Options for radius management.
+        n_evals_is_increased (bool): Whether the number of evaluations was increased.
 
     Returns:
         float: The updated radius.
@@ -25,7 +26,7 @@ def adjust_radius(radius, rho, step_length, options):
 
     if rho >= options.rho_increase and is_large_step:
         new_radius = radius * options.expansion_factor
-    elif rho >= options.rho_decrease:
+    elif rho >= options.rho_decrease or n_evals_is_increased:
         new_radius = radius
     else:
         new_radius = radius * options.shrinking_factor

src/tranquilo/history.py

@@ -173,6 +173,11 @@ def get_fvals(self, x_indices):
         )
         return out
 
+    def get_n_evals(self, x_indices):
+        fvals = self.get_fvals(x_indices)
+        n_evals = {k: len(v) for k, v in fvals.items()}
+        return n_evals
+
     def get_model_data(self, x_indices, average=True):
         if np.isscalar(x_indices):
             x_indices = [x_indices]

src/tranquilo/options.py

@@ -75,6 +75,24 @@ def get_default_n_evals_at_start(noisy):
     return 5 if noisy else 1
 
 
+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."""
 
@@ -111,20 +129,20 @@ class RadiusOptions(NamedTuple):
 
 
 class AcceptanceOptions(NamedTuple):
-    confidence_level: float = 0.8
-    power_level: float = 0.8
+    confidence_level: float = 0.95
+    power_level: float = 0.95
     n_initial: int = 5
-    n_min: int = 5
-    n_max: int = 100
+    n_min: int = 4
+    n_max: int = 50
     min_improvement: float = 0.0
 
 
 class StagnationOptions(NamedTuple):
-    min_relative_step_keep: float = 0.125
+    min_relative_step_keep: float
+    drop: bool
     min_relative_step: float = 0.05
     sample_increment: int = 1
     max_trials: int = 1
-    drop: bool = True
 
 
 class SubsolverOptions(NamedTuple):
@@ -168,6 +186,18 @@ class SamplerOptions(NamedTuple):
     return_info: bool = False
 
 
+class NoiseAdaptationOptions(NamedTuple):
+    rho_noise_n_draws: int = 100
+    high_rho: float = 0.6
+    low_rho: float = 0.1
+    ignore_corelation: bool = True
+    min_share_high_rho: float = 0.7
+    min_share_low_rho: float = 0.9
+    min_n_evals: int = 1
+    max_n_evals: int = 30
+    good_rho_threshold: float = 0.1
+
+
 def update_option_bundle(default_options, user_options=None):
     """Update default options with user options.
 

src/tranquilo/process_arguments.py

@@ -13,7 +13,7 @@
 from tranquilo.history import History
 from tranquilo.options import (
     ConvOptions,
-    StagnationOptions,
+    get_default_stagnation_options,
     StopOptions,
     get_default_acceptance_decider,
     get_default_aggregator,
@@ -22,10 +22,12 @@
     get_default_residualize,
     get_default_model_type,
     get_default_n_evals_at_start,
+    get_default_n_evals_per_point,
     get_default_radius_options,
     get_default_sample_size,
     get_default_search_radius_factor,
     update_option_bundle,
+    NoiseAdaptationOptions,
 )
 from tranquilo.region import Region
 from tranquilo.sample_points import get_sampler
@@ -64,12 +66,13 @@ def process_arguments(
     sample_size=None,
     model_type=None,
     search_radius_factor=None,
-    n_evals_per_point=1,
+    n_evals_per_point=None,
     n_evals_at_start=None,
     seed=925408,
     # bundled advanced options
     radius_options=None,
     stagnation_options=None,
+    noise_adaptation_options=None,
     # component names and related options
     sampler="optimal_hull",
     sampler_options=None,
@@ -110,13 +113,26 @@ def process_arguments(
     x = _process_x(x)
     noisy = _process_noisy(noisy)
     n_cores = _process_n_cores(n_cores)
-    stagnation_options = update_option_bundle(StagnationOptions(), stagnation_options)
-    n_evals_per_point = int(n_evals_per_point)
+    stagnation_options = update_option_bundle(
+        get_default_stagnation_options(noisy),
+        stagnation_options,
+    )
     sampling_rng = _process_seed(seed)
+    simulation_rng = _process_seed(seed + 1000)
+
     n_evals_at_start = _process_n_evals_at_start(
         n_evals_at_start,
         noisy,
     )
+    noise_adaptation_options = update_option_bundle(
+        NoiseAdaptationOptions(), noise_adaptation_options
+    )
+
+    n_evals_per_point = _process_n_evals_per_point(
+        n_evals=n_evals_per_point,
+        noisy=noisy,
+        noise_adaptation_options=noise_adaptation_options,
+    )
 
     # process options that depend on arguments with static defaults
     search_radius_factor = _process_search_radius_factor(search_radius_factor, functype)
@@ -198,11 +214,13 @@ def process_arguments(
         "batch_size": batch_size,
         "target_sample_size": target_sample_size,
         "stagnation_options": stagnation_options,
+        "noise_adaptation_options": noise_adaptation_options,
         "search_radius_factor": search_radius_factor,
         "n_evals_per_point": n_evals_per_point,
         "n_evals_at_start": n_evals_at_start,
         "trustregion": trustregion,
         "sampling_rng": sampling_rng,
+        "simulation_rng": simulation_rng,
         "history": history,
         "sample_points": sample_points,
         "solve_subproblem": solve_subproblem,
@@ -312,3 +330,15 @@ def _process_n_evals_at_start(n_evals, noisy):
         raise ValueError("n_initial_acceptance_evals must be non-negative.")
 
     return out
+
+
+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)
+    else:
+        out = int(n_evals)
+
+    if out < 1:
+        raise ValueError("n_evals_per_point must be non-negative.")
+
+    return out

src/tranquilo/rho_noise.py

@@ -6,14 +6,14 @@
 def simulate_rho_noise(
     xs,
     vector_model,
+    old_vector_model,
     trustregion,
     noise_cov,
     model_fitter,
     model_aggregator,
     subsolver,
     rng,
-    n_draws=100,
-    ignore_corelation=True,
+    options,
 ):
     """Simulate a rho that would obtain on average if there is no approximation error.
 
@@ -29,6 +29,8 @@ def simulate_rho_noise(
         vector_model (VectorModel): A vector surrogate model that is taken as true model
             for the simulation. In many cases this model was fitted on xs but this is
             not a requirement.
+        old_vector_model (VectorModel): A vector surrogate model that is potentially
+            used to fit the new model (if residualize=True).
         trustregion (Region): The trustregion in which the optimization is performed.
         noise_cov(np.ndarray): Covariance matrix of the noise. The noise is assumed to
             be drawn from a multivariate normal distribution with mean zero and this
@@ -38,11 +40,12 @@ def simulate_rho_noise(
             scalar model.
         subsolver (callable): A function that solves the subproblem.
         rng (np.random.Generator): Random number generator.
-        n_draws (int): Number of draws used to estimate the rho noise.
-        ignore_corelation (bool): If True, the noise is assumed to be uncorrelated and
-            only the diagonal entries of the covariance matrix are used.
+        options (NoiseAdaptationOptions): Options for the noise adaptation.
 
     """
+    n_draws = options.rho_noise_n_draws
+    ignore_corelation = options.ignore_corelation
+
     n_samples, n_params = xs.shape
     n_residuals = len(noise_cov)
 
@@ -52,7 +55,9 @@ def simulate_rho_noise(
 
     true_scalar_model = model_aggregator(vector_model=vector_model)
 
-    true_current_fval = true_scalar_model.predict(np.zeros(n_params))
+    true_current_fval = true_scalar_model.predict(
+        trustregion.map_to_unit(trustregion.center)
+    )
 
     if ignore_corelation:
         noise_cov = np.diag(np.diag(noise_cov))
@@ -69,7 +74,7 @@ def simulate_rho_noise(
             sim_fvecs,
             weights=None,
             region=trustregion,
-            old_model=None,
+            old_model=old_vector_model,
         )
         sim_scalar_model = model_aggregator(vector_model=sim_vector_model)
         sim_sub_sol = subsolver(sim_scalar_model, trustregion)