Fix #291: Add constraint checking to initial designs

.gitignore

@@ -6,6 +6,9 @@
 /.mypy_cache/
 /.vscode/
 /.venv/
+/venv/
+env/
+virtualenv/
 
 __pycache__/
 *.py[cod]

emukit/core/initial_designs/base.py

@@ -1,14 +1,18 @@
-# Copyright 2020-2024 The Emukit Authors. All Rights Reserved.
+# Copyright 2020-2026 The Emukit Authors. All Rights Reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 # Copyright 2018-2020 Amazon.com, Inc. or its affiliates. All Rights Reserved.
 # SPDX-License-Identifier: Apache-2.0
 
 
+import logging
+
 import numpy as np
 
 from .. import ParameterSpace
 
+_log = logging.getLogger(__name__)
+
 
 class InitialDesignBase(object):
     """
@@ -18,15 +22,74 @@ class InitialDesignBase(object):
     def __init__(self, parameter_space: ParameterSpace):
         """
         :param parameter_space: The parameter space to generate design for.
-
         """
         self.parameter_space = parameter_space
 
-    def get_samples(self, point_count: int) -> np.ndarray:
+    def _generate_samples(self, point_count: int) -> np.ndarray:
         """
-        Generates requested amount of points.
+        Generate samples without constraint checking. Should be overridden by subclasses.
 
         :param point_count: Number of points required.
         :return: A numpy array of generated samples, shape (point_count x space_dim)
         """
         raise NotImplementedError("Subclasses should implement this method.")
+
+    def _check_constraints(self, samples: np.ndarray) -> np.ndarray:
+        """
+        Check which samples satisfy all constraints.
+
+        :param samples: Array of shape (n_points x n_dims)
+        :return: Boolean array of shape (n_points,) where True indicates the point satisfies all constraints
+        """
+        if not self.parameter_space.constraints:
+            return np.ones(samples.shape[0], dtype=bool)
+
+        # Start with all points being valid
+        valid = np.ones(samples.shape[0], dtype=bool)
+
+        # Check each constraint and keep only points that satisfy all
+        for constraint in self.parameter_space.constraints:
+            constraint_satisfaction = constraint.evaluate(samples)
+            # Ensure we're working with boolean arrays
+            constraint_satisfaction = np.asarray(constraint_satisfaction, dtype=bool)
+            valid = valid & constraint_satisfaction
+
+        return valid
+
+    def get_samples(self, point_count: int, max_retries: int = 100) -> np.ndarray:
+        """
+        Generates requested amount of points that satisfy all constraints.
+        Uses rejection sampling: if any constraints are present and violated,
+        the entire batch is regenerated.
+
+        :param point_count: Number of points required.
+        :param max_retries: Maximum number of retry attempts to generate valid samples when constraints are present.
+                           Default is 100.
+        :return: A numpy array of generated samples, shape (point_count x space_dim)
+        :raises RuntimeError: If unable to generate the required number of valid points after max_retries attempts.
+        """
+        # If there are no constraints, just generate and return
+        if not self.parameter_space.constraints:
+            return self._generate_samples(point_count)
+
+        # With constraints: use rejection sampling
+        for attempt in range(max_retries):
+            candidates = self._generate_samples(point_count)
+            valid_mask = self._check_constraints(candidates)
+
+            if np.all(valid_mask):
+                # All points are valid
+                return candidates
+            else:
+                valid_count = np.sum(valid_mask)
+                _log.debug(
+                    f"Initial design: {valid_count}/{point_count} points satisfy constraints. "
+                    f"Retrying (attempt {attempt + 1}/{max_retries})."
+                )
+
+        # Failed to generate valid samples after all retries
+        raise RuntimeError(
+            f"Could not generate {point_count} valid samples respecting all constraints "
+            f"after {max_retries} attempts. "
+            f"Consider relaxing constraints or increasing max_retries."
+        )

emukit/core/initial_designs/latin_design.py

@@ -26,9 +26,9 @@ def __init__(self, parameter_space: ParameterSpace) -> None:
         """
         super(LatinDesign, self).__init__(parameter_space)
 
-    def get_samples(self, point_count: int) -> np.ndarray:
+    def _generate_samples(self, point_count: int) -> np.ndarray:
         """
-        Generates requested amount of points.
+        Generates requested amount of points (without constraint checking).
 
         :param point_count: Number of points required.
         :return: A numpy array of generated samples, shape (point_count x space_dim)

emukit/core/initial_designs/random_design.py

@@ -23,9 +23,9 @@ def __init__(self, parameter_space: ParameterSpace) -> None:
         """
         super(RandomDesign, self).__init__(parameter_space)
 
-    def get_samples(self, point_count: int) -> np.ndarray:
+    def _generate_samples(self, point_count: int) -> np.ndarray:
         """
-        Generates requested amount of points.
+        Generates requested amount of points (without constraint checking).
 
         :param point_count: Number of points required.
         :return: A numpy array of generated samples, shape (point_count x space_dim)

emukit/core/initial_designs/sobol_design.py

@@ -21,13 +21,13 @@ class SobolDesign(InitialDesignBase):
 
     def __init__(self, parameter_space: ParameterSpace) -> None:
         """
-        param parameter_space: The parameter space to generate design for.
+        :param parameter_space: The parameter space to generate design for.
         """
         super(SobolDesign, self).__init__(parameter_space)
 
-    def get_samples(self, point_count: int) -> np.ndarray:
+    def _generate_samples(self, point_count: int) -> np.ndarray:
         """
-        Generates requested amount of points.
+        Generates requested amount of points (without constraint checking).
 
         :param point_count: Number of points required.
         :return: A numpy array of generated samples, shape (point_count x space_dim)

tests/emukit/core/test_initial_designs.py

@@ -0,0 +1,182 @@
+# Copyright 2020-2026 The Emukit Authors. All Rights Reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+# Copyright 2018-2020 Amazon.com, Inc. or its affiliates. All Rights Reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+import numpy as np
+import pytest
+
+from emukit.core import CategoricalParameter, ContinuousParameter, DiscreteParameter, ParameterSpace
+from emukit.core.constraints import LinearInequalityConstraint, NonlinearInequalityConstraint
+from emukit.core.initial_designs import RandomDesign
+from emukit.core.initial_designs.latin_design import LatinDesign
+from emukit.core.initial_designs.sobol_design import SobolDesign
+
+
+def create_initial_designs(space: ParameterSpace):
+    return [RandomDesign(space), LatinDesign(space), SobolDesign(space)]
+
+
+def test_design_returns_correct_number_of_points():
+    p = ContinuousParameter("c", 1.0, 5.0)
+    space = ParameterSpace([p])
+    points_count = 5
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        assert points_count == len(points)
+        assert all([len(p) == 1 for p in points])
+
+
+def test_design_returns_points_within_bounds():
+    p1 = ContinuousParameter("p1", 0.01, 0.05)
+    p2 = ContinuousParameter("p2", -100.0, -90.0)
+    space = ParameterSpace([p1, p2])
+    points_count = 5
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        for i, p in enumerate(space.parameters):
+            assert np.all(p.min <= points[:, i])
+            assert np.all(points[:, i] <= p.max)
+
+
+def test_design_with_mixed_domain(encoding):
+    p1 = ContinuousParameter("p1", 1.0, 5.0)
+    p2 = CategoricalParameter("p2", encoding)
+    p3 = DiscreteParameter("p3", [1, 2, 5, 6])
+    space = ParameterSpace([p1, p2, p3])
+    points_count = 5
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        assert points_count == len(points)
+        # columns count is 1 for continuous plus 1 for discrete plus number of categories
+        columns_count = 1 + 1 + len(encoding.categories)
+        assert all([len(p) == columns_count for p in points])
+
+
+# Tests for constraint-respecting designs
+
+
+def test_designs_respect_linear_inequality_constraints():
+    """Test that designs respect linear inequality constraints."""
+    p1 = ContinuousParameter("p1", 0.0, 10.0)
+    p2 = ContinuousParameter("p2", 0.0, 10.0)
+
+    # Constraint: p1 + p2 <= 18 (loose enough to be achievable)
+    constraint = LinearInequalityConstraint(
+        constraint_matrix=np.array([[1.0, 1.0]]), lower_bound=np.array([-np.inf]), upper_bound=np.array([18.0])
+    )
+
+    space = ParameterSpace([p1, p2], constraints=[constraint])
+    points_count = 10
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        # Verify all points satisfy the constraint
+        assert points.shape == (points_count, 2)
+        constraint_values = points[:, 0] + points[:, 1]
+        assert np.all(constraint_values <= 18.0 + 1e-6)  # Small tolerance for numerical errors
+
+
+def test_designs_respect_nonlinear_constraints():
+    """Test that designs respect nonlinear constraints."""
+    p1 = ContinuousParameter("p1", 0.0, 5.0)
+    p2 = ContinuousParameter("p2", 0.0, 5.0)
+
+    # Constraint: p1^2 + p2^2 <= 22 (circle of radius ~4.69, covers ~75% of space)
+    # Note: constraint function receives a 1-d array (single point), not 2-d
+    def circle_constraint(x):
+        return x[0] ** 2 + x[1] ** 2
+
+    constraint = NonlinearInequalityConstraint(
+        constraint_function=circle_constraint, lower_bound=np.array([-np.inf]), upper_bound=np.array([22.0])
+    )
+
+    space = ParameterSpace([p1, p2], constraints=[constraint])
+    points_count = 5
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        # Verify all points satisfy the constraint
+        assert points.shape == (points_count, 2)
+        constraint_values = np.array([circle_constraint(p) for p in points])
+        assert np.all(constraint_values <= 22.0 + 1e-6)  # Small tolerance for numerical errors
+
+
+def test_designs_with_multiple_constraints():
+    """Test that designs respect multiple constraints simultaneously."""
+    p1 = ContinuousParameter("p1", 0.0, 10.0)
+    p2 = ContinuousParameter("p2", 0.0, 10.0)
+
+    # Constraint 1: p1 >= 0.5 (loose constraint, 95% of space)
+    constraint1 = LinearInequalityConstraint(
+        constraint_matrix=np.array([[1.0, 0.0]]), lower_bound=np.array([0.5]), upper_bound=np.array([np.inf])
+    )
+
+    # Constraint 2: p2 <= 9.5 (loose constraint, 95% of space)
+    constraint2 = LinearInequalityConstraint(
+        constraint_matrix=np.array([[0.0, 1.0]]), lower_bound=np.array([-np.inf]), upper_bound=np.array([9.5])
+    )
+
+    space = ParameterSpace([p1, p2], constraints=[constraint1, constraint2])
+    points_count = 5
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        points = design.get_samples(points_count)
+
+        # Verify all points satisfy both constraints
+        assert points.shape == (points_count, 2)
+        assert np.all(points[:, 0] >= 0.5 - 1e-6)
+        assert np.all(points[:, 1] <= 9.5 + 1e-6)
+
+
+def test_design_fails_with_impossible_constraints():
+    """Test that design raises error when constraints are impossible to satisfy."""
+    p1 = ContinuousParameter("p1", 0.0, 5.0)
+
+    # Constraint: p1 > 10 (impossible given bounds)
+    constraint = LinearInequalityConstraint(
+        constraint_matrix=np.array([[1.0]]), lower_bound=np.array([10.0]), upper_bound=np.array([np.inf])
+    )
+
+    space = ParameterSpace([p1], constraints=[constraint])
+
+    designs = create_initial_designs(space)
+    for design in designs:
+        with pytest.raises(RuntimeError, match="Could not generate"):
+            design.get_samples(10)
+
+
+def test_design_respects_max_retries():
+    """Test that max_retries parameter controls retry behavior."""
+    p1 = ContinuousParameter("p1", 0.0, 10.0)
+    p2 = ContinuousParameter("p2", 0.0, 10.0)
+
+    # Very restrictive constraint that's hard to satisfy
+    constraint = LinearInequalityConstraint(
+        constraint_matrix=np.array([[1.0, 1.0]]),
+        lower_bound=np.array([19.5]),  # Very close to maximum
+        upper_bound=np.array([20.0]),
+    )
+
+    space = ParameterSpace([p1, p2], constraints=[constraint])
+
+    # Test with all design types
+    designs = create_initial_designs(space)
+    for design in designs:
+        with pytest.raises(RuntimeError, match="Could not generate"):
+            design.get_samples(5, max_retries=10)