DoE objective evaluation in pytorch

bofire/data_models/constraints/interpoint.py

@@ -1,5 +1,5 @@
 import math
-from typing import Annotated, Literal, Optional
+from typing import Annotated, Dict, Literal, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -83,3 +83,6 @@ def __call__(self, experiments: pd.DataFrame) -> pd.Series:
 
     def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
         raise NotImplementedError("Method `jacobian` currently not implemented.")
+
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        raise NotImplementedError("Method `hessian` currently not implemented.")

bofire/data_models/constraints/linear.py

@@ -1,4 +1,4 @@
-from typing import Annotated, List, Literal, Tuple
+from typing import Annotated, Dict, List, Literal, Tuple, Union
 
 import numpy as np
 import pandas as pd
@@ -64,6 +64,9 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             columns=[f"dg/d{name}" for name in self.features],
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[int, str], float]:
+        return {idx: 0.0 for idx in range(experiments.shape[0])}
+
 
 class LinearEqualityConstraint(LinearConstraint, EqualityConstraint):
     """Linear equality constraint of the form `coefficients * x = rhs`.

bofire/data_models/constraints/nchoosek.py

@@ -1,4 +1,4 @@
-from typing import Literal
+from typing import Dict, Literal, Union
 
 import numpy as np
 import pandas as pd
@@ -119,3 +119,6 @@ def is_fulfilled(self, experiments: pd.DataFrame, tol: float = 1e-6) -> pd.Serie
 
     def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
         raise NotImplementedError("Jacobian not implemented for NChooseK constraints.")
+
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        raise NotImplementedError("Hessian not implemented for NChooseK constraints.")

bofire/data_models/constraints/nonlinear.py

@@ -1,6 +1,6 @@
 import inspect
 import warnings
-from typing import Callable, Literal, Optional, Union
+from typing import Callable, Dict, Literal, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -9,6 +9,7 @@
 
 try:
     import torch
+    from torch.autograd.functional import hessian as torch_hessian
     from torch.autograd.functional import jacobian as torch_jacobian
 
     torch_tensor = torch.tensor
@@ -21,6 +22,7 @@ def error_func(*args, **kwargs):
     torch_jacobian = error_func
     torch_tensor = error_func
     torch_diag = error_func
+    torch_hessian = error_func
 
 from bofire.data_models.constraints.constraint import (
     EqualityConstraint,
@@ -47,6 +49,9 @@ class NonlinearConstraint(IntrapointConstraint):
     jacobian_expression: Optional[Union[str, Callable]] = Field(
         default=None, validate_default=True
     )
+    hessian_expression: Optional[Union[str, Callable]] = Field(
+        default=None, validate_default=True
+    )
 
     def validate_inputs(self, inputs: Inputs):
         if self.features is not None:
@@ -100,12 +105,53 @@ def set_jacobian_expression(cls, jacobian_expression, info) -> Union[str, Callab
 
         return jacobian_expression
 
+    @field_validator("hessian_expression")
+    @classmethod
+    def set_hessian_expression(cls, hessian_expression, info) -> Union[str, Callable]:
+        if (
+            hessian_expression is None
+            and "features" in info.data.keys()
+            and "expression" in info.data.keys()
+        ):
+            try:
+                import sympy
+            except ImportError as e:
+                warnings.warn(e.msg)
+                warnings.warn("please run `pip install sympy` for this functionality.")
+                return hessian_expression
+            if info.data["features"] is not None:
+                if isinstance(info.data["expression"], str):
+                    return (
+                        "["
+                        + ", ".join(
+                            [
+                                "["
+                                + ", ".join(
+                                    [
+                                        str(
+                                            sympy.S(info.data["expression"])
+                                            .diff(key1)
+                                            .diff(key2)
+                                        )
+                                        for key1 in info.data["features"]
+                                    ]
+                                )
+                                + "]"
+                                for key2 in info.data["features"]
+                            ]
+                        )
+                        + "]"
+                    )
+
+        return hessian_expression
+
     def __call__(self, experiments: pd.DataFrame) -> pd.Series:
         if isinstance(self.expression, str):
             return experiments.eval(self.expression)
         elif isinstance(self.expression, Callable):
             func_input = {
-                col: torch_tensor(experiments[col]) for col in experiments.columns
+                col: torch_tensor(experiments[col], requires_grad=False)
+                for col in experiments.columns
             }
             return pd.Series(self.expression(**func_input).cpu().numpy())
 
@@ -129,7 +175,10 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             elif isinstance(self.jacobian_expression, Callable):
                 args = inspect.getfullargspec(self.jacobian_expression).args
 
-                func_input = {arg: torch_tensor(experiments[arg]) for arg in args}
+                func_input = {
+                    arg: torch_tensor(experiments[arg], requires_grad=False)
+                    for arg in args
+                }
                 result = self.jacobian_expression(**func_input)
 
                 return pd.DataFrame(
@@ -141,12 +190,14 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
                             for col in experiments.columns
                         ]
                     ),
-                    index=["dg/d" + name for name in args],
+                    index=["dg/d" + name for name in experiments.columns],
                 ).transpose()
         elif isinstance(self.expression, Callable):
             args = inspect.getfullargspec(self.expression).args
 
-            func_input = tuple([torch_tensor(experiments[arg]) for arg in args])
+            func_input = tuple(
+                [torch_tensor(experiments[arg], requires_grad=False) for arg in args]
+            )
 
             result = torch_jacobian(self.expression, func_input)
             result = [torch_diag(result[i]).cpu().numpy() for i in range(len(args))]
@@ -160,6 +211,77 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             "The jacobian of a nonlinear constraint cannot be evaluated if jacobian_expression is None and expression is not Callable.",
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> Dict[Union[str, int], pd.DataFrame]:
+        """
+        Computes a dict of dataframes where the key dimension is the index of the experiments dataframe
+        and the value is the hessian matrix of the constraint evaluated at the corresponding experiment.
+
+        Args:
+            experiments (pd.DataFrame): Dataframe to evaluate the constraint Hessian on.
+
+        Returns:
+            Dict[pd.DataFrame]: Dictionary of dataframes where the key is the index of the experiments dataframe
+            and the value is the Hessian matrix of the constraint evaluated at the corresponding experiment.
+        """
+        if self.hessian_expression is not None:
+            if isinstance(self.hessian_expression, str):
+                res = experiments.eval(self.hessian_expression)
+            else:
+                if not isinstance(self.hessian_expression, Callable):
+                    raise ValueError(
+                        "The hessian_expression must be a string or a callable.",
+                    )
+                args = inspect.getfullargspec(self.hessian_expression).args
+
+                func_input = {
+                    arg: torch_tensor(experiments[arg], requires_grad=False)
+                    for arg in args
+                }
+                res = self.hessian_expression(**func_input)
+            for i, _ in enumerate(res):
+                for j, entry in enumerate(res[i]):
+                    if not hasattr(entry, "__iter__"):
+                        res[i][j] = pd.Series(np.repeat(entry, experiments.shape[0]))
+            res = np.array(res)
+            names = self.features or [f"x{i}" for i in range(experiments.shape[1])]
+
+            return {
+                idx: pd.DataFrame(
+                    res[..., i],
+                    columns=[f"d/d{name}" for name in names],
+                    index=[f"d/d{name}" for name in names],
+                )
+                for i, idx in enumerate(experiments.index)
+            }
+
+        elif isinstance(self.expression, Callable):
+            args = inspect.getfullargspec(self.expression).args
+
+            func_input = {
+                idx: tuple(
+                    [
+                        torch_tensor(experiments[arg][idx], requires_grad=False)
+                        for arg in args
+                    ]
+                )
+                for idx in experiments.index
+            }
+
+            names = self.features or [f"x{i}" for i in range(experiments.shape[1])]
+            res = {
+                idx: pd.DataFrame(
+                    np.array(torch_hessian(self.expression, func_input[idx])),
+                    columns=[f"d/d{name}" for name in names],
+                    index=[f"d/d{name}" for name in names],
+                )
+                for idx in experiments.index
+            }
+            return res
+
+        raise ValueError(
+            "The hessian of a nonlinear constraint cannot be evaluated if hessian_expression is None and expression is not Callable.",
+        )
+
 
 class NonlinearEqualityConstraint(NonlinearConstraint, EqualityConstraint):
     """Nonlinear equality constraint of the form 'expression == 0'.

bofire/data_models/constraints/product.py

@@ -83,6 +83,11 @@ def jacobian(self, experiments: pd.DataFrame) -> pd.DataFrame:
             "Jacobian for product constraints is not yet implemented.",
         )
 
+    def hessian(self, experiments: pd.DataFrame) -> List[pd.DataFrame]:
+        raise NotImplementedError(
+            "Hessian for product constraints is not yet implemented.",
+        )
+
 
 class ProductEqualityConstraint(ProductConstraint, EqualityConstraint):
     """Represents a product constraint of the form `sign * x1**e1 * x2**e2 * ... * xn**en == rhs`.

bofire/data_models/strategies/doe.py

@@ -112,6 +112,7 @@ class DoEStrategy(Strategy):
 
     verbose: bool = False  # get rid of this at a later stage
     ipopt_options: Optional[Dict] = None
+    use_hessian: bool = False
 
     @classmethod
     def is_constraint_implemented(cls, my_type: Type[Constraint]) -> bool:

bofire/strategies/doe/branch_and_bound.py

@@ -235,6 +235,7 @@ def find_local_max_ipopt_BaB(
     categorical_groups: Optional[List[List[ContinuousInput]]] = None,
     discrete_variables: Optional[Dict[str, Tuple[ContinuousInput, List[float]]]] = None,
     verbose: bool = False,
+    use_hessian: bool = False,
 ) -> pd.DataFrame:
     """Function computing a d-optimal design" for a given domain and model.
     It allows for the problem to have categorical values which is solved by Branch-and-Bound
@@ -261,6 +262,7 @@ def find_local_max_ipopt_BaB(
             verbose (bool): if true, print information during the optimization process
             transform_range (Optional[Bounds]): range to which the input variables are transformed.
                 If None is provided, the features will not be scaled. Defaults to None.
+            use_hessian (bool): if true, the hessian will be used in the optimization process. Defaults to False.
         Returns:
             A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
             local optimum.
@@ -324,6 +326,7 @@ def find_local_max_ipopt_BaB(
         sampling,
         None,
         partially_fixed_experiments=initial_branch,
+        use_hessian=use_hessian,
     )
     initial_value = objective_function.evaluate(
         initial_design.to_numpy().flatten(),
@@ -367,6 +370,7 @@ def find_local_max_ipopt_exhaustive(
     categorical_groups: Optional[List[List[ContinuousInput]]] = None,
     discrete_variables: Optional[Dict[str, Tuple[ContinuousInput, List[float]]]] = None,
     verbose: bool = False,
+    use_hessian: bool = False,
 ) -> pd.DataFrame:
     """Function computing a d-optimal design" for a given domain and model.
     It allows for the problem to have categorical values which is solved by exhaustive search
@@ -392,6 +396,7 @@ def find_local_max_ipopt_exhaustive(
                 with key:(relaxed variable, valid values). Defaults to None
             verbose (bool): if true, print information during the optimization process
             transform_range (Optional[Bounds]): range to which the input variables are transformed.
+            use_hessian (bool): if true, the hessian will be used in the optimization process. Defaults to False.
         Returns:
             A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
             local optimum.
@@ -503,6 +508,7 @@ def find_local_max_ipopt_exhaustive(
                 sampling,
                 None,
                 one_set_of_experiments,
+                use_hessian=use_hessian,
             )
             domain.validate_candidates(
                 candidates=current_design.apply(lambda x: np.round(x, 8)),

bofire/strategies/doe/design.py

@@ -4,7 +4,6 @@
 import numpy as np
 import pandas as pd
 from formulaic import Formula
-from scipy.optimize._minimize import standardize_constraints
 
 from bofire.data_models.constraints.api import (
     ConstraintNotFulfilledError,
@@ -17,6 +16,7 @@
 from bofire.data_models.strategies.doe import AnyOptimalityCriterion
 from bofire.strategies.doe.objective import get_objective_function
 from bofire.strategies.doe.utils import (
+    _minimize,
     constraints_as_scipy_constraints,
     nchoosek_constraints_as_bounds,
 )
@@ -31,6 +31,7 @@ def find_local_max_ipopt(
     sampling: Optional[pd.DataFrame] = None,
     fixed_experiments: Optional[pd.DataFrame] = None,
     partially_fixed_experiments: Optional[pd.DataFrame] = None,
+    use_hessian: bool = False,
 ) -> pd.DataFrame:
     """Function computing an optimal design for a given domain and model.
 
@@ -48,6 +49,7 @@ def find_local_max_ipopt(
             Variables can be fixed to one value or can be set to a range by setting a tuple with lower and upper bound
             Non-fixed variables have to be set to None or nan.
         criterion (OptimalityCriterion): OptimalityCriterion object indicating which criterion function to use.
+        use_hessian (bool): If True, the hessian of the objective function is used. Default is False.
 
     Returns:
         A pd.DataFrame object containing the best found input for the experiments. In general, this is only a
@@ -58,14 +60,6 @@ def find_local_max_ipopt(
     # Checks and preparation steps
     #
 
-    # warn user if IPOPT scipy interface is not available
-    try:
-        from cyipopt import minimize_ipopt  # type: ignore
-    except ImportError:
-        raise ImportError(
-            "cyipopt is not installed. Install it via `conda install -c conda-forge cyipopt`"
-        )
-
     objective_function = get_objective_function(
         criterion, domain=domain, n_experiments=n_experiments
     )
@@ -166,28 +160,26 @@ def find_local_max_ipopt(
     # set ipopt options
     if ipopt_options is None:
         ipopt_options = {}
-    _ipopt_options = {"maxiter": 500, "disp": 0}
+    _ipopt_options = {"max_iter": 500, "print_level": 0}
     for key in ipopt_options.keys():
         _ipopt_options[key] = ipopt_options[key]
-    if _ipopt_options["disp"] > 12:
-        _ipopt_options["disp"] = 0
+    if _ipopt_options["print_level"] > 12:
+        _ipopt_options["print_level"] = 0
 
     #
     # Do the optimization
     #
-
-    result = minimize_ipopt(
-        objective_function.evaluate,
+    x = _minimize(
+        objective_function=objective_function,
         x0=x0,
         bounds=bounds,
-        # "SLSQP" has no deeper meaning here and just ensures correct constraint standardization
-        constraints=standardize_constraints(constraints, x0, "SLSQP"),
-        options=_ipopt_options,
-        jac=objective_function.evaluate_jacobian,
+        constraints=constraints,
+        use_hessian=use_hessian,
+        ipopt_options=_ipopt_options,
     )
 
     design = pd.DataFrame(
-        result["x"].reshape(n_experiments, len(domain.inputs)),
+        x.reshape(n_experiments, len(domain.inputs)),
         columns=domain.inputs.get_keys(),
         index=[f"exp{i}" for i in range(n_experiments)],
     )