feat: added option to use fixed epsilon (#163)

* feat: fixed epsilon implemented in core

* chore: added tests for constant epsilon

Author: kjappelbaum

docs/background.rst

@@ -7,6 +7,8 @@ This implementation has the following features:
 
 - We ensure that the sampling is scale-invariant and that the algorithm can deal with positive and negative objective values.
 
+- In contrast to the `original implementation by Zuluaga et al. <https://jmlr.org/papers/v17/15-047.html>`_ we do not assume that the range of the objectives is known a priori. In their implementation it is used to calculate fixed tolerance values :math:`\epsilon_i \cdot r_i` (where :math:`r_i` is the range of objective :math:`i`). We instead use by default :math:`\epsilon_i \cdot |\mu_i|`.
+
 - Instead of using the predicted :math:`\hat{\mu}` and :math:`\hat{\sigma}` also for the sampled points we use the measured :math:`\mu` and :math:`\sigma`.
 
 - This implementation is directly scalable to :math:`n`-dimensional problems.

pyepal/pal/core.py

@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint:disable=anomalous-backslash-in-string
 # Copyright 2020 PyePAL authors
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
@@ -147,15 +148,16 @@ def _union_one_dim(
     return np.array(out_lows), np.array(out_ups)
 
 
-def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
+def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals, too-many-branches
     pareto_optimal_0: np.array,
     not_pareto_optimal_0: np.array,
     unclassified_0: np.array,
     rectangle_lows: np.array,
     rectangle_ups: np.array,
     epsilon: np.array,
+    is_fixed_epsilon: bool = False,
 ) -> Tuple[np.array, np.array, np.array]:
-    """Performs the classification part of the algorithm
+    """Performs the classification of the algorithm
     (p. 4 of the PAL paper, see algorithm 1/2 of the epsilon-PAL paper)
 
     One core concept is that once a point is classified,
@@ -172,6 +174,13 @@ def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
         rectangle_lows (np.array): lower uncertainty boundaries
         rectangle_ups (np.array): upper uncertainty boundaries
         epsilon (np.array): granularity parameter (one per dimension)
+        is_fixed_epsilon (bool): If true it assumes that epsilon contains *absolute*
+            tolerance values for every objective. These would typically be calculated as
+            :math:`\epsilon_i = \varepsilon_i \cdot r_i`, where :math:`r_i` is the range
+            of objective :math:`i`. By default this is False. This is, we will use
+            :math:`\epsilon_i = \varepsilon_i \cdot y_i` to compute the objectives,
+            which hence avoids the need of knowing the range of a objective before
+            using the algorithm.
 
     Returns:
         Tuple[list, list, list]: binary encoded list of Pareto optimal,
@@ -186,11 +195,17 @@ def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
     if sum(pareto_optimal_0) > 0:
         pareto_indices = np.where(pareto_optimal_0)[0]
         pareto_pessimistic_lows = rectangle_lows[pareto_indices]  # p_pess(P)
+
+        if is_fixed_epsilon:
+            tolerances_0 = np.tile(epsilon, (len(pareto_pessimistic_lows), 1))
+        else:
+            tolerances_0 = np.abs(epsilon * pareto_pessimistic_lows)
+
         for i in range(0, len(unclassified_0)):
             if unclassified_t[i] == 1:
                 # discard if any lower-bound epsilon dominates the upper bound
                 if dominance_check_jitted_2(
-                    pareto_pessimistic_lows + np.abs(epsilon * pareto_pessimistic_lows),
+                    pareto_pessimistic_lows + tolerances_0,
                     rectangle_ups[i],
                 ):
                     not_pareto_optimal_t[i] = True
@@ -207,14 +222,21 @@ def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
     pareto_unclassified_pessimistic_points = pareto_unclassified_lows[
         pareto_unclassified_pessimistic_mask
     ]
+
+    if is_fixed_epsilon:
+        tolerances_1 = np.tile(
+            epsilon, (len(pareto_unclassified_pessimistic_points), 1)
+        )
+    else:
+        tolerances_1 = epsilon * np.abs(pareto_unclassified_pessimistic_points)
+
     for i in range(0, len(unclassified_t)):  # pylint:disable=consider-using-enumerate
         # We can only discard points that are unclassified so far
         # We cannot discard points that are part of p_pess(P \cup U)
         if unclassified_t[i] and (i not in original_indices):
             # discard if any lower-bound epsilon dominates the upper bound
             if dominance_check_jitted_2(
-                epsilon * np.abs(pareto_unclassified_pessimistic_points)
-                + pareto_unclassified_pessimistic_points,
+                tolerances_1 + pareto_unclassified_pessimistic_points,
                 rectangle_ups[i],
             ):
                 not_pareto_optimal_t[i] = True
@@ -228,6 +250,11 @@ def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
 
     index_map = {index: i for i, index in enumerate(unclassified_indices)}
 
+    if is_fixed_epsilon:
+        tolerances_2 = np.tile(epsilon, (len(rectangle_lows), 1))
+    else:
+        tolerances_2 = epsilon * np.abs(rectangle_lows)
+
     # The index map helps us to mask the current point from the unclassified_ups list
     for i in range(0, len(unclassified_t)):  # pylint:disable=consider-using-enumerate
         # again, we only care about unclassified points
@@ -237,7 +264,7 @@ def _pareto_classify(  # pylint:disable=too-many-arguments, too-many-locals
             # the current point is epsilon-accurate Pareto optimal
             if not dominance_check_jitted_3(
                 unclassified_ups,
-                rectangle_lows[i] + epsilon * np.abs(rectangle_lows[i]),
+                rectangle_lows[i] + tolerances_2[i],
                 index_map[i],
             ):
                 pareto_optimal_t[i] = True

pyepal/pal/pal_base.py

@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint:disable=anomalous-backslash-in-string
 # Copyright 2020 PyePAL authors
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
@@ -40,6 +41,7 @@
     validate_epsilon,
     validate_goals,
     validate_ndim,
+    validate_ranges,
 )
 
 PAL_LOGGER = logging.getLogger("PALLogger")
@@ -53,7 +55,7 @@
 __all__ = ["PALBase", "PAL_LOGGER"]
 
 
-class PALBase:  # pylint:disable=too-many-instance-attributes
+class PALBase:  # pylint:disable=too-many-instance-attributes, too-many-public-methods
     """PAL base class"""
 
     def __init__(  # pylint:disable=too-many-arguments
@@ -66,6 +68,7 @@ def __init__(  # pylint:disable=too-many-arguments
         beta_scale: float = 1 / 9,
         goals: List[str] = None,
         coef_var_threshold: float = 3,
+        ranges: Union[np.ndarray, None] = None,
     ):
         """Initialize the PAL instance
 
@@ -87,6 +90,11 @@ def __init__(  # pylint:disable=too-many-arguments
             coef_var_threshold (float, optional): Use only points with
                 a coefficient of variation below this threshold
                 in the classification step. Defaults to 3.
+            ranges (np.ndarray, optional): Numpy array of length ndmin,
+                where each element contains the value range of given objective.
+                If this is provided, we will use :math:`\epsilon \cdot ranges`
+                to computer the uncertainties of the hyperrectangles instead
+                of the default behavior :math:`\epsilon \cdot |\mu|`
 
         """
         self.cross_val_points = 10  # maybe we make it an argument at some point
@@ -112,6 +120,7 @@ def __init__(  # pylint:disable=too-many-arguments
         self.design_space = X_design
         self.beta = None
         self.goals = validate_goals(goals, ndim)
+        self.ranges = validate_ranges(ranges, ndim)
 
         # self.y is what needs to be used for train/predict
         # as there the data has been turned into maximization
@@ -128,6 +137,16 @@ def __repr__(self):
         {self.number_discarded_points} discarded points, \
         {self.number_unclassified_points} unclassified points."
 
+    def _uses_fixed_epsilon(self):
+        if self.ranges is not None:
+            return True
+        return False
+
+    @property
+    def uses_fixed_epsilon(self):
+        """True if it uses the fixed epsilon :math:`\epsilon \cdot ranges`"""
+        return self._uses_fixed_epsilon()
+
     def _reset(self):
         self.pareto_optimal = np.array([False] * self.number_design_points)
         self.discarded = np.array([False] * self.number_design_points)
@@ -373,15 +392,26 @@ def _update_coef_var_mask(self):
 
     def _classify(self):
         self._update_coef_var_mask()
-        pareto_optimal, discarded, unclassified = _pareto_classify(
-            self.pareto_optimal[self.coef_var_mask],
-            self.discarded[self.coef_var_mask],
-            self.unclassified[self.coef_var_mask],
-            self.rectangle_lows[self.coef_var_mask],
-            self.rectangle_ups[self.coef_var_mask],
-            self.epsilon,
-        )
-
+        if self.uses_fixed_epsilon:
+            pareto_optimal, discarded, unclassified = _pareto_classify(
+                self.pareto_optimal[self.coef_var_mask],
+                self.discarded[self.coef_var_mask],
+                self.unclassified[self.coef_var_mask],
+                self.rectangle_lows[self.coef_var_mask],
+                self.rectangle_ups[self.coef_var_mask],
+                self.epsilon * self.ranges,
+                is_fixed_epsilon=True,
+            )
+        else:
+            pareto_optimal, discarded, unclassified = _pareto_classify(
+                self.pareto_optimal[self.coef_var_mask],
+                self.discarded[self.coef_var_mask],
+                self.unclassified[self.coef_var_mask],
+                self.rectangle_lows[self.coef_var_mask],
+                self.rectangle_ups[self.coef_var_mask],
+                self.epsilon,
+                is_fixed_epsilon=False,
+            )
         self.pareto_optimal[self.coef_var_mask] = pareto_optimal
         self.discarded[self.coef_var_mask] = discarded
         self.unclassified[self.coef_var_mask] = unclassified

pyepal/pal/validate_inputs.py

@@ -18,7 +18,7 @@
 import collections
 import warnings
 from copy import deepcopy
-from typing import Any, Iterable, List, Sequence
+from typing import Any, Iterable, List, Sequence, Union
 
 import numpy as np
 from sklearn.gaussian_process import GaussianProcessRegressor
@@ -446,3 +446,20 @@ def validate_positive_integer_list(
             raise ValueError("{} must be a positive integer".format(parameter_name))
 
     return seq
+
+
+def validate_ranges(ranges: Any, ndim: int) -> Union[None, np.ndarray]:
+    """Make sure that it has the correct numnber of elements and that all
+    elements are positive."""
+    if not isinstance(ranges, (np.ndarray, list)):
+        return None
+
+    if not len(ranges) == ndim:
+        raise ValueError(
+            "The number of elements in ranges must match the number of objectives."
+        )
+    for elem in ranges:
+        if not elem > 0:
+            raise ValueError("Ranges must be positive.")
+
+    return np.array(ranges)

tests/test_pal_base.py

@@ -52,6 +52,10 @@ def test_pal_base(make_random_dataset):
         palinstance.sample()
 
     assert palinstance.y.shape == (100, 3)
+    assert not palinstance.uses_fixed_epsilon
+
+    palinstance = PALBase(make_random_dataset[0], ["model"], 3, ranges=[1, 1, 1])
+    assert palinstance.uses_fixed_epsilon
 
 
 def test_update_train_set(make_random_dataset):

tests/test_pal_core.py

@@ -273,6 +273,28 @@ def test_pareto_classify(binh_korn_points):  # pylint:disable=too-many-locals
         == np.array([False, False, False, False, False, False, False, True])
     ).all()
 
+    pareto_optimal_t, discarded_t, unclassified_t = _pareto_classify(
+        is_pareto_optimal,
+        is_discarded,
+        is_unclassified,
+        rectangle_lows,
+        rectangle_ups,
+        np.array([0.1, 0.1]),
+        is_fixed_epsilon=True,
+    )
+
+    assert (
+        pareto_optimal_t
+        == np.array([True, True, True, False, True, False, False, False])
+    ).all()
+    assert (
+        discarded_t == np.array([False, False, False, True, False, True, True, False])
+    ).all()
+    assert (
+        unclassified_t
+        == np.array([False, False, False, False, False, False, False, True])
+    ).all()
+
     # 3D arrays, but 3rd dimenension alsways 0
 
     pareto_optimal_points = np.array([[0.5, 2, 0], [3, 1, 0], [4, 0.5, 0]])

tests/test_pal_sklearn.py

@@ -277,6 +277,34 @@ def test_orchestration_run_one_step_batch(  # pylint:disable=too-many-statements
     for model in palinstance.models:
         assert check_is_fitted(model) is None
 
+    # test using the "fixed" epsilon
+    gpr_0 = GaussianProcessRegressor(
+        RBF(), normalize_y=True, n_restarts_optimizer=6, random_state=10
+    )
+    gpr_1 = GaussianProcessRegressor(
+        RBF(), normalize_y=True, n_restarts_optimizer=6, random_state=10
+    )
+    palinstance = PALSklearn(
+        X_binh_korn,
+        [gpr_0, gpr_1],
+        2,
+        beta_scale=1 / 9,
+        ranges=np.ptp(y_binh_korn, axis=0),
+    )
+    assert palinstance.uses_fixed_epsilon
+    palinstance.cross_val_points = 0
+    sample_idx = np.array([1, 10, 20, 40, 70, 90])
+    palinstance.update_train_set(sample_idx, y_binh_korn[sample_idx])
+    idx = palinstance.run_one_step(batch_size=1)
+    for index in idx:
+        assert index not in [1, 10, 20, 40, 70, 90]
+    assert palinstance.number_sampled_points > 0
+    assert sum(palinstance.unclassified) > 0
+    assert sum(palinstance.discarded) == 0
+
+    for model in palinstance.models:
+        assert check_is_fitted(model) is None
+
 
 def test_orchestration_run_one_step_parallel(binh_korn_points):
     """Test the parallel processing"""

tests/test_validate_inputs.py

@@ -41,6 +41,7 @@
     validate_number_models,
     validate_optimizers,
     validate_positive_integer_list,
+    validate_ranges,
 )
 
 
@@ -316,3 +317,18 @@ def test_validate_positive_integer_list():
         validate_positive_integer_list(-1, 2)
 
     assert validate_positive_integer_list(1, 2) == [1, 1]
+
+
+def test_validate_ranges():
+    """Check that the range validation works"""
+    arr = np.array([1, 1, 1])
+    assert (arr == validate_ranges(arr, 3)).all()
+
+    with pytest.raises(ValueError):
+        validate_ranges(arr, 2)
+
+    with pytest.raises(ValueError):
+        validate_ranges(np.array([-0.1, 0.1, 1]), 2)
+
+    assert validate_ranges(None, 3) is None
+    assert (validate_ranges([1, 1, 1], 3) == np.array([1, 1, 1])).all()