update

examples/HVN/convex_example.py

@@ -52,7 +52,7 @@ def g_Hessian(x):
 ref = np.array([1, 1])
 max_iters = 10
 
-N = 10
+N = 15
 x1 = np.random.rand(N)
 x2 = x1**2 - 2 * x1 + 1
 x0 = np.c_[x1, x2]
@@ -75,6 +75,7 @@ def g_Hessian(x):
     minimization=True,
     max_iters=max_iters,
     verbose=True,
+    preconditioning=True,
 )
 X, Y, stop = opt.run()
 
@@ -112,8 +113,8 @@ def g_Hessian(x):
 ax0.legend([r"$Y_{\text{final}}$", r"$Y_0$", "Pareto front"])
 
 ax12 = ax1.twinx()
-ax1.plot(range(1, len(opt.history_HV) + 1), opt.history_HV, "b-")
-ax12.semilogy(range(1, len(opt.history_HV) + 1), opt.hist_R_norm, "g--")
+ax1.plot(range(1, len(opt.history_indicator_value) + 1), opt.history_indicator_value, "b-")
+ax12.semilogy(range(1, len(opt.history_R_norm) + 1), opt.history_R_norm, "g--")
 ax1.set_ylabel("HV", color="b")
 ax12.set_ylabel(r"$||G(\mathbf{X})||$", color="g")
 ax1.set_title("Performance")

examples/HV_Hessian/2D_example_Hessian.py

@@ -21,8 +21,8 @@ def MOP1_Hessian(x):
 
 ref = np.array([20, 20])
 hvh = HypervolumeDerivatives(
-    n_decision_var=2,
-    n_objective=2,
+    n_var=2,
+    n_obj=2,
     ref=ref,
     func=MOP1,
     jac=MOP1_Jacobian,
@@ -33,7 +33,7 @@ def MOP1_Hessian(x):
     X = np.random.rand(10, 2)
     Y = np.array([MOP1(_) for _ in X])
     idx = get_non_dominated(-1 * Y, return_index=True)
-    out = hvh.compute(X[idx, :])
+    out = hvh._compute_hessian(X[idx, :])
     HdY2, HdX2 = out["HVdY2"], out["HVdX2"]
     HdX = out["HVdX"]
     w, v = np.linalg.eigh(HdX2)

examples/HV_Hessian/3D_example_Hessian.py

@@ -16,14 +16,14 @@
     ref = np.array([2, 2, 2])
 
     hvh = HypervolumeDerivatives(
-        n_decision_var=dim,
-        n_objective=3,
+        n_var=dim,
+        n_obj=3,
         ref=ref,
         func=f.objective,
         jac=f.objective_jacobian,
         hessian=f.objective_hessian,
     )
-    out = hvh.compute(X)
+    out = hvh._compute_hessian(X)
     # automatic differentiation is very slow;
     # used for verifying the result of the analytical computation
     AD = hvh.compute_automatic_differentiation(X)

examples/HV_Hessian/4D_example_Hessian.py

@@ -10,7 +10,7 @@
 ref = np.array([17, 35, 7, 10])
 hvh = HypervolumeDerivatives(4, 4, ref, minimization=True)
 X = np.array([(16, 23, 1, 8), (14, 32, 2, 5), (12, 27, 3, 1), (10, 21, 4, 9), (8, 33, 5, 3)])
-out = hvh.compute(X)
+out = hvh._compute_hessian(X)
 
 assert np.all(
     out["HVdY2"]

examples/MMD/example_ZDT.py

@@ -74,6 +74,7 @@
     max_iters=max_iters,
     verbose=True,
     metrics=metrics,
+    preconditioning=True,
 )
 X, Y, _ = opt.run()
 

hvd/base.py

@@ -1,6 +1,7 @@
-import numpy as np
+from copy import deepcopy
+from typing import Callable, Self, Tuple
 
-from .utils import merge_lists
+import numpy as np
 
 
 class State:
@@ -9,24 +10,31 @@ def __init__(
         n_var: int,
         n_eq: int,
         n_ieq: int,
-        func: callable,
-        jac: callable,
-        h: callable = None,
-        h_jac: callable = None,
-        g: callable = None,
-        g_jac: callable = None,
+        func: Callable,
+        jac: Callable,
+        h: Callable = None,
+        h_jac: Callable = None,
+        h_hess: Callable = None,
+        g: Callable = None,
+        g_jac: Callable = None,
+        g_hess: Callable = None,
     ) -> None:
         self.n_var: int = n_var  # the number of the primal variables
         self.n_eq: int = n_eq
         self.n_ieq: int = n_ieq
-        self.func: callable = func
-        self._jac: callable = jac
-        self.h: callable = h
-        self.h_jac: callable = h_jac
-        self.g: callable = g
-        self.g_jac: callable = g_jac
+        self.n_cstr: int = self.n_eq + self.n_ieq
+        self.func: Callable = func
+        self._jac: Callable = jac
+        self.h: Callable = h
+        self.h_jac: Callable = h_jac
+        self.h_hess: Callable = h_hess
+        self.g: Callable = g
+        self.g_jac: Callable = g_jac
+        self.g_hess: Callable = g_hess
         self._constrained: bool = self.g is not None or self.h is not None
         self.n_jac_evals: int = 0
+        self.n_cstr_jac_evals: int = 0
+        self.n_cstr_hess_evals: int = 0
 
     @property
     def N(self):
@@ -37,32 +45,67 @@ def jac(self, x: np.ndarray) -> np.ndarray:
         self.n_jac_evals += 1
         return self._jac(x)
 
+    def evaluate_cstr(self, x: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
+        H = self.h(x) if self.h is not None else []
+        G = self.g(x) if self.g is not None else []
+        if isinstance(H, (int, float)):
+            H = np.array([H])
+        if isinstance(G, (int, float)):
+            G = np.array([G])
+        active_indices = [True] * self.n_eq
+        if self.g is not None:
+            active_indices += (G >= -1e-4).tolist()
+        return np.r_[H, G], np.array(active_indices)
+
+    def evaluate_cstr_jac(self, x: np.ndarray) -> np.ndarray:
+        if self.g_jac is not None or self.h_jac is not None:
+            self.n_cstr_jac_evals += 1
+        # TODO: only evaluate active inequality constraints
+        dH = self.h_jac(x).tolist() if self.h_jac is not None else []
+        dG = self.g_jac(x).tolist() if self.g_jac is not None else []
+        return np.array(dH + dG).reshape(self.n_cstr, self.n_var)
+
+    def evaluate_cstr_hess(self, x: np.ndarray) -> np.ndarray:
+        if self.g_hess is not None or self.h_hess is not None:
+            self.n_cstr_hess_evals += 1
+        # TODO: only evaluate active inequality constraints
+        ddH = self.h_hess(x).tolist() if self.h_hess is not None else []
+        ddG = self.g_hess(x).tolist() if self.g_hess is not None else []
+        return np.array(ddH + ddG).reshape(self.n_cstr, self.n_var, self.n_var)
+
     def update(self, X: np.ndarray, compute_gradient: bool = True):
         self.X = X
         primal_vars = self.primal
         self.Y = np.array([self.func(x) for x in primal_vars])
         self.J = np.array([self.jac(x) for x in primal_vars]) if compute_gradient else None
-        if self._constrained:
-            eq = self._evaluate_constraints(primal_vars, type="eq", compute_gradient=compute_gradient)
-            ieq = self._evaluate_constraints(primal_vars, type="ieq", compute_gradient=compute_gradient)
-            cstr_value, active_indices, dH = merge_lists(eq, ieq)
-            self.cstr_value = cstr_value
-            self.active_indices = active_indices
-            self.dH = dH
+        cstr_value, active_indices = list(zip(*[self.evaluate_cstr(x) for x in primal_vars]))
+        self.cstr_value = np.array(cstr_value)
+        self.active_indices = np.array(active_indices)
+        self.cstr_grad = np.array([self.evaluate_cstr_jac(x) for x in primal_vars])
+        self.cstr_hess = np.array([self.evaluate_cstr_hess(x) for x in primal_vars])
 
     def update_one(self, x: np.ndarray, k: int):
+        primal_vars = x[: self.n_var]
         self.X[k] = x
-        x = np.atleast_2d(x)
-        primal_vars = x[:, : self.n_var]
-        self.Y[k] = self.func(primal_vars[0])
-        self.J[k] = self.jac(primal_vars[0])
-        if self._constrained:
-            eq = self._evaluate_constraints(primal_vars, type="eq", compute_gradient=True)
-            ieq = self._evaluate_constraints(primal_vars, type="ieq", compute_gradient=True)
-            cstr_value, active_indices, dH = merge_lists(eq, ieq)
+        self.Y[k] = self.func(primal_vars)
+        self.J[k] = self.jac(primal_vars)
+        cstr_value, active_indices = self.evaluate_cstr(primal_vars)
+        cstr_grad = self.evaluate_cstr_jac(primal_vars)
+        cstr_hess = self.evaluate_cstr_hess(primal_vars)
+        if cstr_value is not None:
             self.cstr_value[k] = cstr_value
             self.active_indices[k] = active_indices
-            self.dH[k] = dH
+        if cstr_grad is not None:
+            self.cstr_grad[k] = cstr_grad
+        if cstr_hess is not None:
+            self.cstr_hess[k] = cstr_hess
+
+    def is_feasible(self) -> np.ndarray:
+        # TODO: should the active ieq cstr be considered infeasible?
+        ieq_active_idx = self.active_indices[:, self.n_eq :]
+        eq_feasible = np.all(np.isclose(self.H, 0, atol=1e-4, rtol=0), axis=1)
+        ieq_feasible = ~np.any(ieq_active_idx, axis=1)
+        return np.bitwise_and(eq_feasible, ieq_feasible)
 
     @property
     def primal(self) -> np.ndarray:
@@ -84,18 +127,19 @@ def G(self) -> np.ndarray:
         """Inequality constraint values"""
         return self.cstr_value[:, self.n_eq :]
 
-    def _evaluate_constraints(
-        self,
-        primal_vars: np.ndarray,
-        type: str = "eq",
-        compute_gradient: bool = True,
-    ):
-        N = len(primal_vars)
-        func = self.h if type == "eq" else self.g
-        jac = self.h_jac if type == "eq" else self.g_jac
-        value = np.array([func(x) for x in primal_vars]).reshape(N, -1)
-        active_indices = [[True] * self.n_eq] * N if type == "eq" else [v >= -1e-4 for v in value]
-        active_indices = np.array(active_indices)
-        if compute_gradient:
-            dH = np.array([jac(x).reshape(-1, self.n_var) for x in primal_vars])
-        return (value, active_indices, dH) if compute_gradient else (value, active_indices)
+    def __getitem__(self, indices: np.ndarray) -> Self:
+        """Slicing the state class
+
+        Args:
+            indices (np.ndarray): indices to select
+        """
+        obj = deepcopy(self)
+        obj.n_jac_evals = 0
+        obj.X = obj.X[indices]
+        obj.Y = obj.Y[indices]
+        obj.J = obj.J[indices]
+        if self._constrained:
+            obj.cstr_value = obj.cstr_value[indices]
+            obj.active_indices = obj.active_indices[indices]
+            obj.cstr_grad = obj.cstr_grad[indices]
+        return obj