I optimality criterion

bofire/data_models/constraints/api.py

@@ -4,7 +4,7 @@
  Constraint,
  ConstraintError,
  ConstraintNotFulfilledError,
- EqalityConstraint,
+ EqualityConstraint,
  InequalityConstraint,
  IntrapointConstraint,
 )
@@ -37,7 +37,7 @@
  InterpointConstraint,
  ProductConstraint,
  InequalityConstraint,
- EqalityConstraint,
+ EqualityConstraint,
 ]
 
 AnyConstraint = Union[

bofire/data_models/constraints/constraint.py

@@ -69,7 +69,7 @@ class IntrapointConstraint(Constraint):
  type: str
 
 
-class EqalityConstraint(IntrapointConstraint):
+class EqualityConstraint(IntrapointConstraint):
  type: str
 
  def is_fulfilled(self, experiments: pd.DataFrame, tol: float = 1e-6) -> pd.Series:

bofire/data_models/constraints/linear.py

@@ -5,7 +5,7 @@
 from pydantic import Field, model_validator
 
 from bofire.data_models.constraints.constraint import (
- EqalityConstraint,
+ EqualityConstraint,
  InequalityConstraint,
  IntrapointConstraint,
 )
@@ -64,7 +64,7 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
  )
 
 
-class LinearEqualityConstraint(LinearConstraint, EqalityConstraint):
+class LinearEqualityConstraint(LinearConstraint, EqualityConstraint):
  """Linear equality constraint of the form `coefficients * x = rhs`.
 
  Attributes:

bofire/data_models/constraints/nonlinear.py

@@ -6,7 +6,7 @@
 from pydantic import Field, field_validator
 
 from bofire.data_models.constraints.constraint import (
- EqalityConstraint,
+ EqualityConstraint,
  InequalityConstraint,
  IntrapointConstraint,
 )
@@ -89,7 +89,7 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
  )
 
 
-class NonlinearEqualityConstraint(NonlinearConstraint, EqalityConstraint):
+class NonlinearEqualityConstraint(NonlinearConstraint, EqualityConstraint):
  """Nonlinear equality constraint of the form 'expression == 0'.
 
  Attributes:

bofire/data_models/constraints/product.py

@@ -5,7 +5,7 @@
 from pydantic import Field, model_validator
 
 from bofire.data_models.constraints.constraint import (
- EqalityConstraint,
+ EqualityConstraint,
  InequalityConstraint,
  IntrapointConstraint,
 )
@@ -84,7 +84,7 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
  )
 
 
-class ProductEqualityConstraint(ProductConstraint, EqalityConstraint):
+class ProductEqualityConstraint(ProductConstraint, EqualityConstraint):
  """
  Represents a product constraint of the form `sign * x1**e1 * x2**e2 * ... * xn**en == rhs`.
 

bofire/strategies/doe/design.py

@@ -19,7 +19,7 @@
 from bofire.data_models.features.api import ContinuousInput, Input
 from bofire.data_models.strategies.api import RandomStrategy as RandomStrategyDataModel
 from bofire.data_models.types import Bounds
-from bofire.strategies.doe.objective import get_objective_class
+from bofire.strategies.doe.objective import Objective, get_objective_class
 from bofire.strategies.doe.utils import (
  constraints_as_scipy_constraints,
  get_formula_from_string,
@@ -482,14 +482,17 @@ def find_local_max_ipopt(
  )
 
  # get objective function and its jacobian
- objective_class = get_objective_class(objective)
- objective_function = objective_class(
- domain=domain,
- model=model_formula,
- n_experiments=n_experiments,
- delta=delta,
- transform_range=transform_range,
- )
+ if isinstance(objective, Objective):
+ objective_function = objective
+ else:
+ objective_class = get_objective_class(objective)
+ objective_function = objective_class(
+ domain=domain,
+ model=model_formula,
+ n_experiments=n_experiments,
+ delta=delta,
+ transform_range=transform_range,
+ )
 
  # write constraints as scipy constraints
  constraints = constraints_as_scipy_constraints(

bofire/strategies/doe/objective.py

@@ -1,16 +1,30 @@
+import warnings
 from abc import abstractmethod
 from copy import deepcopy
+from itertools import product
 from typing import Optional, Type
 
 import numpy as np
 import pandas as pd
 import torch
+from cyipopt import minimize_ipopt # type: ignore
 from formulaic import Formula
+from scipy.optimize._minimize import standardize_constraints
 from torch import Tensor
 
+from bofire.data_models.constraints.api import (
+ ConstraintNotFulfilledError,
+ EqualityConstraint,
+)
 from bofire.data_models.domain.api import Domain
+from bofire.data_models.enum import SamplingMethodEnum
 from bofire.data_models.types import Bounds
 from bofire.strategies.doe.transform import IndentityTransform, MinMaxTransform
+from bofire.strategies.doe.utils import (
+ constraints_as_scipy_constraints,
+ get_formula_from_string,
+ nchoosek_constraints_as_bounds,
+)
 from bofire.strategies.enum import OptimalityCriterionEnum
 from bofire.utils.torch_tools import tkwargs
 
@@ -474,6 +488,205 @@ def _convert_input_to_tensor(
  return torch.tensor(X.values, requires_grad=requires_grad, **tkwargs)
 
 
+class IOptimality(Objective):
+ """A class implementing the evaluation of I-criterion and its jacobian w.r.t. the inputs."""
+
+ def __init__(
+ self,
+ domain: Domain,
+ model: Formula,
+ n_experiments: int,
+ delta: float = 1e-6,
+ transform_range: Optional[Bounds] = None,
+ n_space: Optional[int] = None,
+ ipopt_options: Optional[dict] = None,
+ ) -> None:
+ """
+ Args:
+ domain (Domain): A domain defining the DoE domain together with model_type.
+ model_type (str or Formula): A formula containing all model terms.
+ n_experiments (int): Number of experiments
+ delta (float): A regularization parameter for the information matrix. Default value is 1e-3.
+ transform_range (Bounds, optional): range to which the input variables are transformed before applying the objective function. Default is None.
+ n_space (int, optional): Number of space filling points. If none is provided,
+ n_space = n_experiments is assumed. Only relevant if SpaceFilling is used
+ to fill the feasible space, i.e. in presence of equality constraints.
+ Otherwise a uniform grid is generated. If None is provided, n_space = 10 * n_experiments is assumed.
+ Default is None.
+ ipopt_options (dict, optional): Options for the Ipopt solver to generate space filling point.
+ If None is provided, the default options (maxiter = 500) are used.
+ """
+
+ if transform_range is not None:
+ raise ValueError(
+ "IOptimality does not support transformations of the input variables."
+ )
+
+ super().__init__(
+ domain=domain,
+ model=model,
+ n_experiments=n_experiments,
+ delta=delta,
+ transform_range=transform_range,
+ )
+
+ # uniformly fill the design space
+ if np.any([isinstance(obj, EqualityConstraint) for obj in domain.constraints]):
+ warnings.warn(
+ "Equality constraints were detected. No equidistant grid of points can be generated. The design space will be filled via SpaceFilling.",
+ UserWarning,
+ )
+ if n_space is None:
+ n_space = n_experiments * 10
+
+ print(
+ f"Filling the design space for the I-optimality criterion with {n_space} points..."
+ )
+ x0 = (
+ domain.inputs.sample(n=n_space, method=SamplingMethodEnum.UNIFORM)
+ .to_numpy()
+ .flatten()
+ )
+ objective = SpaceFilling(
+ domain, model, n_space, delta, transform_range=None
+ )
+ constraints = constraints_as_scipy_constraints(
+ domain, n_space, ignore_nchoosek=True
+ )
+ bounds = nchoosek_constraints_as_bounds(domain, n_space)
+ if ipopt_options is None:
+ ipopt_options = {}
+ _ipopt_options = {"maxiter": 500, "disp": 0}
+ for key in ipopt_options.keys():
+ _ipopt_options[key] = ipopt_options[key]
+ if _ipopt_options["disp"] > 12:
+ _ipopt_options["disp"] = 0
+
+ result = minimize_ipopt(
+ objective.evaluate,
+ x0=x0,
+ bounds=bounds,
+ constraints=standardize_constraints(constraints, x0, "SLSQP"),
+ options=_ipopt_options,
+ jac=objective.evaluate_jacobian,
+ )
+
+ design = pd.DataFrame(
+ result["x"].reshape(n_space, len(domain.inputs)),
+ columns=domain.inputs.get_keys(),
+ index=[f"exp{i}" for i in range(n_space)],
+ )
+ else:
+ low, high = domain.inputs.get_bounds(specs={})
+ points = [
+ list(
+ np.linspace(
+ low[i],
+ high[i],
+ int(100 * (high[i] - low[i])),
+ )
+ )
+ for i in range(len(low))
+ ]
+ points = np.array(list(product(*points)))
+ points = pd.DataFrame(points, columns=domain.inputs.get_keys())
+ if len(domain.constraints) > 0:
+ fulfilled = domain.constraints(experiments=points)
+ fulfilled = np.array(
+ [
+ np.array(fulfilled.iloc[:, i]) <= 0.0
+ for i in range(fulfilled.shape[1])
+ ]
+ )
+ fulfilled = np.array(np.prod(fulfilled, axis=0), dtype=bool)
+ design = points[fulfilled]
+ else:
+ design = points
+ n_space = len(design)
+ print(
+ f"Filling the design space with {len(design)} equally spaced grid points."
+ )
+
+ try:
+ domain.validate_candidates(
+ candidates=design.apply(lambda x: np.round(x, 8)),
+ only_inputs=True,
+ tol=1e-4,
+ )
+ except (ValueError, ConstraintNotFulfilledError):
+ warnings.warn(
+ "Some points do not lie inside the domain or violate constraints. Please check if the \
+ results lie within your tolerance.",
+ UserWarning,
+ )
+
+ model_formula = get_formula_from_string(
+ model_type=model, rhs_only=True, domain=domain
+ )
+ X = model_formula.get_model_matrix(design).to_numpy()
+ self.space_filling_design = design
+
+ self.YtY = torch.from_numpy(X.T @ X) / n_space
+ self.YtY.requires_grad = False
+
+ print("Done!")
+
+ def _evaluate(self, x: np.ndarray) -> float:
+ """Computes trace((Y.T@Y) / nY @ inv(X.T@X + delta)).
+ Where X is the model matrix corresponding to x, Y is the model matrix of points
+ uniformly filling up the feasible space and nY is the number of such points.
+
+ Args:
+ x (np.ndarray): values of design variables a 1d array.
+
+ Returns:
+ trace((Y.T@Y) / nY @ inv(X.T@X + delta))
+
+ """
+ X = self._convert_input_to_model_tensor(x, requires_grad=False)
+ return float(
+ torch.trace(
+ self.YtY.detach()
+ @ torch.linalg.inv(
+ X.detach().T @ X.detach()
+ + self.delta * torch.eye(self.n_model_terms)
+ )
+ )
+ )
+
+ def _evaluate_jacobian(self, x: np.ndarray) -> np.ndarray:
+ """Computes the jacobian of trace((Y.T@Y) / nY @ inv(X.T@X + delta)).
+
+ Args:
+ x (np.ndarray): values of design variables a 1d array.
+
+ Returns:
+ The jacobian of trace((Y.T@Y) / nY @ inv(X.T@X + delta))
+
+ """
+ # get model matrix X
+ X = self._convert_input_to_model_tensor(x, requires_grad=True)
+
+ # first part of jacobian
+ torch.trace(
+ self.YtY.detach()
+ @ torch.linalg.inv(X.T @ X + self.delta * torch.eye(self.n_model_terms))
+ ).backward()
+ J1 = X.grad.detach().numpy() # type: ignore
+ J1 = np.repeat(J1, self.n_vars, axis=0).reshape(
+ self.n_experiments, self.n_vars, self.n_model_terms
+ )
+
+ # second part of jacobian
+ J2 = self._model_jacobian_t(x)
+
+ # combine both parts
+ J = J1 * J2
+ J = np.sum(J, axis=-1)
+
+ return J.flatten()
+
+
 def get_objective_class(objective: OptimalityCriterionEnum) -> Type:
  objective = OptimalityCriterionEnum(objective)
 
@@ -489,3 +702,5 @@ def get_objective_class(objective: OptimalityCriterionEnum) -> Type:
  return KOptimality
  elif objective == OptimalityCriterionEnum.SPACE_FILLING:
  return SpaceFilling
+ elif objective == OptimalityCriterionEnum.I_OPTIMALITY:
+ return IOptimality

bofire/strategies/enum.py

@@ -8,3 +8,4 @@ class OptimalityCriterionEnum(Enum):
  G_OPTIMALITY = "G_OPTIMALITY"
  K_OPTIMALITY = "K_OPTIMALITY"
  SPACE_FILLING = "SPACE_FILLING"
+ I_OPTIMALITY = "I_OPTIMALITY"