batch constraints for DoE and __call__ for interpointEqualityConstrai…

…nt (experimental-design#322)

bofire/data_models/constraints/interpoint.py

@@ -43,7 +43,31 @@ def is_fulfilled(
         return pd.Series([True])
 
     def __call__(self, experiments: pd.DataFrame) -> pd.Series:
-        raise NotImplementedError("Method `__call__` currently not implemented.")
+        """Numerically evaluates the constraint. Returns the distance to the constraint fulfillment
+        for each batch of size batch_size.
+
+        Args:
+            experiments (pd.DataFrame): Dataframe to evaluate the constraint on.
+
+        Returns:
+            pd.Series: Distance to reach constraint fulfillment.
+        """
+        multiplicity = self.multiplicity or len(experiments)
+        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]
+        feature_values = feature_values.reshape(n_batches, multiplicity).T
+
+        batchwise_constraint_matrix = np.zeros(shape=(multiplicity - 1, multiplicity))
+        batchwise_constraint_matrix[:, 0] = 1.0
+        batchwise_constraint_matrix[:, 1:] = -np.eye(multiplicity - 1)
+
+        return pd.Series(
+            np.linalg.norm(batchwise_constraint_matrix @ feature_values, axis=0, ord=2)
+            ** 2,
+            index=[f"batch_{i}" for i in range(n_batches)],
+        )
 
     def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
         raise NotImplementedError("Method `jacobian` currently not implemented.")

bofire/data_models/constraints/linear.py

@@ -6,12 +6,12 @@
 
 from bofire.data_models.constraints.constraint import (
     Coefficients,
-    Constraint,
     FeatureKeys,
+    IntrapointConstraint,
 )
 
 
-class LinearConstraint(Constraint):
+class LinearConstraint(IntrapointConstraint):
     """Abstract base class for linear equality and inequality constraints.
 
     Attributes:

bofire/data_models/constraints/nchoosek.py

@@ -4,14 +4,14 @@
 import pandas as pd
 from pydantic import root_validator, validator
 
-from bofire.data_models.constraints.constraint import Constraint, FeatureKeys
+from bofire.data_models.constraints.constraint import FeatureKeys, IntrapointConstraint
 
 
 def narrow_gaussian(x, ell=1e-3):
     return np.exp(-0.5 * (x / ell) ** 2)
 
 
-class NChooseKConstraint(Constraint):
+class NChooseKConstraint(IntrapointConstraint):
     """NChooseK constraint that defines how many ingredients are allowed in a formulation.
 
     Attributes:

bofire/data_models/constraints/nonlinear.py

@@ -5,10 +5,10 @@
 import pandas as pd
 from pydantic import validator
 
-from bofire.data_models.constraints.constraint import Constraint, FeatureKeys
+from bofire.data_models.constraints.constraint import FeatureKeys, IntrapointConstraint
 
 
-class NonlinearConstraint(Constraint):
+class NonlinearConstraint(IntrapointConstraint):
     """Base class for nonlinear equality and inequality constraints.
 
     Attributes:

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,22 @@ 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, LinearEqualityConstraint) or isinstance(
+            c, 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, NonlinearEqualityConstraint) or isinstance(
+            c, 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 +224,130 @@ 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, LinearEqualityConstraint) or isinstance(
+        c, 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, NonlinearEqualityConstraint) or isinstance(
+        c, 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
+        multiplicity = c.multiplicity or len(domain.inputs)
+        n_batches = int(np.ceil(n_experiments / multiplicity))
+        lb = np.zeros(n_batches * (multiplicity - 1))
+        ub = np.zeros(n_batches * (multiplicity - 1))
+
+        # write constraint as matrix
+        feature_idx = 0
+        if c.feature not in domain.inputs.get_keys():
+            raise ValueError(f"Feature {c.feature} is not part of the domain {domain}.")
+        for i, name in enumerate(domain.inputs.get_keys()):
+            if name == c.feature:
+                feature_idx = i
+
+        A = np.zeros(shape=(n_batches * (multiplicity - 1), D * n_experiments))
+        for batch in range(n_batches):
+            for i in range(multiplicity - 1):
+                if batch * multiplicity + i + 2 <= n_experiments:
+                    A[
+                        batch * (multiplicity - 1) + i,
+                        batch * multiplicity * D + feature_idx,
+                    ] = 1.0
+                    A[
+                        batch * (multiplicity - 1) + i,
+                        (batch * multiplicity + i + 1) * D + feature_idx,
+                    ] = -1.0
+
+        # remove overflow in last batch
+        if (n_experiments % multiplicity) != 0:
+            n_overflow = multiplicity - (n_experiments % 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

@@ -6,6 +6,9 @@
 import tests.bofire.data_models.specs.api as specs
 from bofire.data_models.constraints.api import (
     Constraint,
+    InterpointConstraint,
+    InterpointEqualityConstraint,
+    IntrapointConstraint,
     LinearConstraint,
     LinearEqualityConstraint,
     LinearInequalityConstraint,
@@ -78,6 +81,7 @@ def test_as_smaller_equal():
 c6 = LinearInequalityConstraint.from_smaller_equal(
     features=["f1", "f2", "f3"], coefficients=[1, 1, 1], rhs=100.0
 )
+c7 = InterpointEqualityConstraint(feature="f2", multiplicity=2)
 
 if1 = ContinuousInput(key="f1", bounds=(0, 2))
 if2 = ContinuousInput(key="f2", bounds=(0, 4))
@@ -89,6 +93,8 @@ def test_as_smaller_equal():
 constraints2 = Constraints(constraints=[c4, c5])
 constraints3 = Constraints(constraints=[c6])
 constraints4 = Constraints(constraints=[c3])
+constraints5 = Constraints(constraints=[c7])
+constraints6 = Constraints(constraints=[c1, c7])
 
 
 @pytest.mark.parametrize(
@@ -133,12 +139,27 @@ def test_constraints_plus():
     [
         (constraints2, 5),
         (constraints4, 5),
+        (constraints5, 5),
+        (constraints6, 5),
     ],
 )
 def test_constraints_call(constraints, num_candidates):
     candidates = inputs.sample(num_candidates, SamplingMethodEnum.UNIFORM)
     returned = constraints(candidates)
-    assert returned.shape == (num_candidates, len(constraints))
+
+    num_rows = 0
+    if np.any([isinstance(c, IntrapointConstraint) for c in constraints]):
+        num_rows += num_candidates
+
+    max_num_batches = 0
+    for c in constraints:
+        if isinstance(c, InterpointConstraint):
+            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))
 
 
 @pytest.mark.parametrize(