add support for batch constraints in DoE

bofire/data_models/constraints/interpoint.py

@@ -53,7 +53,7 @@ def __call__(self, experiments: pd.DataFrame) -> pd.Series:
             pd.Series: Distance to reach constraint fulfillment.
         """
         multiplicity = self.multiplicity or len(experiments)
-        n_batches = (experiments.shape[0] + self.multiplicity - 1) // multiplicity
+        n_batches = int(np.ceil((experiments.shape[0] / multiplicity)))
         feature_values = np.zeros(n_batches * multiplicity)
         feature_values[: experiments.shape[0]] = experiments[self.feature].values
         feature_values[experiments.shape[0] :] = feature_values[-multiplicity]

bofire/strategies/doe/utils.py

@@ -8,12 +8,14 @@
 from scipy.optimize import LinearConstraint, NonlinearConstraint
 
 from bofire.data_models.constraints.api import (
+    Constraint,
+    InterpointEqualityConstraint,
     LinearEqualityConstraint,
     LinearInequalityConstraint,
     NChooseKConstraint,
     NonlinearEqualityConstraint,
+    NonlinearInequalityConstraint,
 )
-from bofire.data_models.constraints.nonlinear import NonlinearInequalityConstraint
 from bofire.data_models.domain.domain import Domain
 from bofire.data_models.features.continuous import ContinuousInput
 from bofire.data_models.strategies.api import (
@@ -187,7 +189,7 @@ def constraints_as_scipy_constraints(
     n_experiments: int,
     ignore_nchoosek: bool = True,
 ) -> List:
-    """Formulates opti constraints as scipy constraints.
+    """Formulates bofire constraints as scipy constraints.
 
     Args:
         domain (Domain): Domain whose constraints should be formulated as scipy constraints.
@@ -197,80 +199,20 @@ def constraints_as_scipy_constraints(
     Returns:
         A list of scipy constraints corresponding to the constraints of the given opti problem.
     """
-    D = len(domain.inputs)
 
     # reformulate constraints
     constraints = []
     if len(domain.constraints) == 0:
         return constraints
     for c in domain.constraints:
-        if isinstance(c, LinearEqualityConstraint):
-            # write lower/upper bound as vector
-            lb = np.ones(n_experiments) * (c.rhs / np.linalg.norm(c.coefficients))
-            ub = np.ones(n_experiments) * (c.rhs / np.linalg.norm(c.coefficients))
-
-            # write constraint as matrix
-            lhs = {
-                c.features[i]: c.coefficients[i] / np.linalg.norm(c.coefficients)
-                for i in range(len(c.features))
-            }
-            row = np.zeros(D)
-            for i, name in enumerate(domain.inputs.get_keys()):
-                if name in lhs.keys():
-                    row[i] = lhs[name]
-
-            A = np.zeros(shape=(n_experiments, D * n_experiments))
-            for i in range(n_experiments):
-                A[i, i * D : (i + 1) * D] = row
-
-            constraints.append(LinearConstraint(A, lb, ub))  # type: ignore
-
-        elif isinstance(c, LinearInequalityConstraint):
-            # write upper/lowe bound as vector
-            lb = -np.inf * np.ones(n_experiments)
-            ub = np.ones(n_experiments) * c.rhs / np.linalg.norm(c.coefficients)
-
-            # write constraint as matrix
-            lhs = {
-                c.features[i]: c.coefficients[i] / np.linalg.norm(c.coefficients)
-                for i in range(len(c.features))
-            }
-            row = np.zeros(D)
-            for i, name in enumerate(domain.inputs.get_keys()):
-                if name in lhs.keys():
-                    row[i] = lhs[name]
-
-            A = np.zeros(shape=(n_experiments, D * n_experiments))
-            for i in range(n_experiments):
-                A[i, i * D : (i + 1) * D] = row
-
+        if isinstance(c, Union[LinearEqualityConstraint, LinearInequalityConstraint]):
+            A, lb, ub = get_constraint_function_and_bounds(c, domain, n_experiments)
             constraints.append(LinearConstraint(A, lb, ub))  # type: ignore
 
-        elif isinstance(c, NonlinearEqualityConstraint):
-            # write upper/lower bound as vector
-            lb = np.zeros(n_experiments)
-            ub = np.zeros(n_experiments)
-
-            # define constraint evaluation (and gradient if provided)
-            fun = ConstraintWrapper(
-                constraint=c, domain=domain, n_experiments=n_experiments
-            )
-
-            if c.jacobian_expression is not None:
-                constraints.append(NonlinearConstraint(fun, lb, ub, jac=fun.jacobian))
-            else:
-                constraints.append(NonlinearConstraint(fun, lb, ub))
-
-        elif isinstance(c, NonlinearInequalityConstraint):
-            # write upper/lower bound as vector
-            lb = -np.inf * np.ones(n_experiments)
-            ub = np.zeros(n_experiments)
-
-            # define constraint evaluation (and gradient if provided)
-            fun = ConstraintWrapper(
-                constraint=c, domain=domain, n_experiments=n_experiments
-            )
-
+        elif isinstance(
+            c, Union[NonlinearEqualityConstraint, NonlinearInequalityConstraint]
+        ):
+            fun, lb, ub = get_constraint_function_and_bounds(c, domain, n_experiments)
             if c.jacobian_expression is not None:
                 constraints.append(NonlinearConstraint(fun, lb, ub, jac=fun.jacobian))
             else:
@@ -280,23 +222,124 @@ def constraints_as_scipy_constraints(
             if ignore_nchoosek:
                 pass
             else:
-                # write upper/lower bound as vector
-                lb = -np.inf * np.ones(n_experiments)
-                ub = np.zeros(n_experiments)
-
-                # define constraint evaluation (and gradient if provided)
-                fun = ConstraintWrapper(
-                    constraint=c, domain=domain, n_experiments=n_experiments
+                fun, lb, ub = get_constraint_function_and_bounds(
+                    c, domain, n_experiments
                 )
-
                 constraints.append(NonlinearConstraint(fun, lb, ub, jac=fun.jacobian))
 
+        elif isinstance(c, InterpointEqualityConstraint):
+            A, lb, ub = get_constraint_function_and_bounds(c, domain, n_experiments)
+            constraints.append(LinearConstraint(A, lb, ub))  # type: ignore
+
         else:
             raise NotImplementedError(f"No implementation for this constraint: {c}")
 
     return constraints
 
 
+def get_constraint_function_and_bounds(
+    c: Constraint, domain: Domain, n_experiments: int
+) -> List:
+    """Returns the function definition and bounds for a given constraint and domain.
+
+    Args:
+        c (Constraint): Constraint for which the constraint matrix should be determined.
+        domain (Domain): Domain for which the constraint matrix should be determined.
+        n_experiments (int): Number of experiments for which the constraint matrix should be determined.
+
+    Returns:
+        A list containing the constraint defining function and the lower and upper bounds.
+    """
+    D = len(domain.inputs)
+
+    if isinstance(c, Union[LinearEqualityConstraint, LinearInequalityConstraint]):
+        # write constraint as matrix
+        lhs = {
+            c.features[i]: c.coefficients[i] / np.linalg.norm(c.coefficients)
+            for i in range(len(c.features))
+        }
+        row = np.zeros(D)
+        for i, name in enumerate(domain.inputs.get_keys()):
+            if name in lhs.keys():
+                row[i] = lhs[name]
+
+        A = np.zeros(shape=(n_experiments, D * n_experiments))
+        for i in range(n_experiments):
+            A[i, i * D : (i + 1) * D] = row
+
+        # write upper/lower bound as vector
+        lb = -np.inf * np.ones(n_experiments)
+        ub = np.ones(n_experiments) * (c.rhs / np.linalg.norm(c.coefficients))
+        if isinstance(c, LinearEqualityConstraint):
+            lb = np.ones(n_experiments) * (c.rhs / np.linalg.norm(c.coefficients))
+
+        return [A, lb, ub]
+
+    elif isinstance(
+        c, Union[NonlinearEqualityConstraint, NonlinearInequalityConstraint]
+    ):
+        # define constraint evaluation (and gradient if provided)
+        fun = ConstraintWrapper(
+            constraint=c, domain=domain, n_experiments=n_experiments
+        )
+
+        # write upper/lower bound as vector
+        lb = -np.inf * np.ones(n_experiments)
+        ub = np.zeros(n_experiments)
+        if isinstance(c, NonlinearEqualityConstraint):
+            lb = np.zeros(n_experiments)
+
+        return [fun, lb, ub]
+
+    elif isinstance(c, NChooseKConstraint):
+        # define constraint evaluation (and gradient if provided)
+        fun = ConstraintWrapper(
+            constraint=c, domain=domain, n_experiments=n_experiments
+        )
+
+        # write upper/lower bound as vector
+        lb = -np.inf * np.ones(n_experiments)
+        ub = np.zeros(n_experiments)
+
+        return [fun, lb, ub]
+
+    elif isinstance(c, InterpointEqualityConstraint):
+        # write lower/upper bound as vector
+        n_batches = int(np.ceil(n_experiments / c.multiplicity))
+        lb = np.zeros(n_batches * (c.multiplicity - 1))
+        ub = np.zeros(n_batches * (c.multiplicity - 1))
+
+        # write constraint as matrix
+        for i, name in enumerate(domain.inputs.get_keys()):
+            if name == c.feature:
+                feature_idx = i
+
+        A = np.zeros(shape=(n_batches * (c.multiplicity - 1), D * n_experiments))
+        for batch in range(n_batches):
+            for i in range(c.multiplicity - 1):
+                if batch * c.multiplicity + i + 2 <= n_experiments:
+                    A[
+                        batch * (c.multiplicity - 1) + i,
+                        batch * c.multiplicity * D + feature_idx,
+                    ] = 1.0
+                    A[
+                        batch * (c.multiplicity - 1) + i,
+                        (batch * c.multiplicity + i + 1) * D + feature_idx,
+                    ] = -1.0
+
+        # remove overflow in last batch
+        if (n_experiments % c.multiplicity) != 0:
+            n_overflow = c.multiplicity - (n_experiments % c.multiplicity)
+            A = A[:-n_overflow, :]
+            lb = lb[:-n_overflow]
+            ub = ub[:-n_overflow]
+
+        return [A, lb, ub]
+
+    else:
+        raise NotImplementedError(f"No implementation for this constraint: {c}")
+
+
 class ConstraintWrapper:
     """Wrapper for nonlinear constraints."""
 

tests/bofire/data_models/test_constraints.py

@@ -154,7 +154,9 @@ def test_constraints_call(constraints, num_candidates):
     max_num_batches = 0
     for c in constraints:
         if isinstance(c, InterpointConstraint):
-            max_num_batches = max(max_num_batches, num_candidates // c.multiplicity + 1)
+            max_num_batches = max(
+                max_num_batches, int(np.ceil(num_candidates / c.multiplicity))
+            )
     num_rows += max_num_batches
 
     assert returned.shape == (num_rows, len(constraints))

tests/bofire/strategies/doe/test_utils.py

@@ -5,6 +5,7 @@
 from scipy.optimize import LinearConstraint, NonlinearConstraint
 
 from bofire.data_models.constraints.api import (
+    InterpointEqualityConstraint,
     LinearEqualityConstraint,
     LinearInequalityConstraint,
     NChooseKConstraint,
@@ -326,6 +327,31 @@ def test_constraints_as_scipy_constraints():
         constraints[0].fun(np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0])), [2, 0, 0]
     )
 
+    # domain with batch constraint
+    domain = Domain.from_lists(
+        inputs=[ContinuousInput(key=f"x{i}", bounds=(0, 1)) for i in range(3)],
+        outputs=[ContinuousOutput(key="y")],
+        constraints=[
+            InterpointEqualityConstraint(feature="x0", multiplicity=3),
+        ],
+    )
+    n_experiments = 5
+
+    constraints = constraints_as_scipy_constraints(domain, n_experiments)
+    assert len(constraints) == 1
+    assert isinstance(constraints[0], LinearConstraint)
+    A = np.array(
+        [
+            [1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0],
+        ],
+        dtype=float,
+    )
+    assert np.allclose(constraints[0].A, A)
+    assert np.allclose(constraints[0].lb, np.zeros(3))
+    assert np.allclose(constraints[0].ub, np.zeros(3))
+
 
 def test_ConstraintWrapper():
     # define domain with all types of constraints