update

hvd/delta_p.py

@@ -55,7 +55,7 @@ def compute(self, X: np.ndarray = None, Y: np.ndarray = None) -> float:
         return np.mean(self.D[np.arange(len(Y)), self.indices] ** self.p) ** (1 / self.p)
 
     def compute_derivatives(
-        self, X: np.ndarray, Y: np.ndarray = None, compute_hessian: bool = True
+        self, X: np.ndarray, Y: np.ndarray = None, compute_hessian: bool = True, Jacobian: np.ndarray = None, **kwargs,
     ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
         """compute the derivatives of the generational distance^p
 
@@ -76,7 +76,12 @@ def compute_derivatives(
             Y = np.array([self.func(x) for x in X])
 
         self._compute_indices(Y)
-        J = np.array([self.jac(x) for x in X])  # (N, n_objective, dim)
+        # Jacobian of the objective function
+        if Jacobian is None:
+            J = np.array([self.jac(x) for x in X])  # (N, n_obj, dim)
+        else:
+            J = Jacobian
+        # J = np.array([self.jac(x) for x in X])  # (N, n_objective, dim)
         diff = Y - self.ref[self.indices]  # (N, n_objective)
         diff_norm = np.sqrt(np.sum(diff**2, axis=1)).reshape(-1, 1)  # (N, 1)
         grad_ = np.einsum("ijk,ij->ik", J, diff)  # (N, dim)
@@ -139,6 +144,7 @@ def match(self, Y: Union[dict, np.ndarray], Y_idx: Union[List, None] = None) ->
             out[Y_idx[v]] = medoids
             for i, j in enumerate(Y_idx[v]):
                 self._medoids_idx[j] = (k, i)
+        # TODO: make `out` a property
         return out
 
     def set_medoids(self, medroid: np.ndarray, k: int):
@@ -249,6 +255,7 @@ def _match(self, Y: np.ndarray, Y_label: np.ndarray = None):
         Y_ = Y
         Y_idx = None
         # if the reference set is clustered, then also try to cluster the approximation set
+        # TODO: implement clustering of `Y`
         if Y_label is not None:
             n_cluster = len(np.unique(Y_label))
             Y_idx = [np.nonzero(Y_label == i)[0] for i in range(n_cluster)]

hvd/newton.py

@@ -12,15 +12,8 @@
 from .delta_p import GenerationalDistance, InvertedGenerationalDistance
 from .hypervolume import hypervolume
 from .hypervolume_derivatives import HypervolumeDerivatives
-from .utils import (
-    compute_chim,
-    get_logger,
-    get_non_dominated,
-    merge_lists,
-    non_domin_sort,
-    precondition_hessian,
-    set_bounds,
-)
+from .utils import (compute_chim, get_logger, get_non_dominated, merge_lists,
+                    non_domin_sort, precondition_hessian, set_bounds)
 
 np.seterr(divide="ignore", invalid="ignore")
 
@@ -822,7 +815,7 @@ def _compute_R(
         primal_vars, dual_vars = state.primal, state.dual
         if grad is None:
             grad = self.active_indicator.compute_derivatives(
-                primal_vars, state.Y, self.Y_label, compute_hessian=False, Jacobian=state.J
+                X=primal_vars, Y=state.Y, Y_label=self.Y_label, compute_hessian=False, Jacobian=state.J
             )
         R = grad  # the unconstrained case
         dH, idx = None, None
@@ -839,7 +832,7 @@ def _compute_netwon_step(self) -> Tuple[np.ndarray, np.ndarray]:
         R = np.zeros((self.N, self.dim))  # the root-finding problem
         # gradient and Hessian of the incumbent indicator
         grad, Hessian = self.active_indicator.compute_derivatives(
-            primal_vars, self.state.Y, self.Y_label, Jacobian=self.state.J
+            X=primal_vars, Y=self.state.Y, Jacobian=self.state.J, Y_label=self.Y_label, 
         )
         # the root-finding problem and the gradient of the active constraints
         R_list, dH, active_indices = self._compute_R(self.state, grad=grad)
@@ -872,7 +865,7 @@ def _shift_reference_set(self):
         1. always shift the first reference set (`self.iter_count == 0`); Otherwise, weird matching can happen
         2. if at least one approximation point is close to its matched target and the Newton step is not zero.
         """
-        distance = np.linalg.norm(self.state.Y - self.active_indicator._medoids, axis=1)
+        distance = np.linalg.norm(self.state.Y - self._igd._medoids, axis=1)
         masks = (
             np.array([True] * self.N)
             if self.iter_count == 0
@@ -897,15 +890,17 @@ def _shift_reference_set(self):
         for i, k in enumerate(indices):
             n = self._eta[self.Y_label[k]]
             # NOTE: initial shift CF1: 0.6, CF2/3: 0.2
+            # DTLZ4: 0.08 seems to work a bit better
+            # TODO: create a configuration class to set those hyperparameter of this method, e.g., shift amount
             v = 0.05 * n if self.iter_count > 0 else 0.2 * n  # the initial shift is a bit larger
-            self.active_indicator.shift_medoids(v, k)
+            self._igd.shift_medoids(v, k)
 
         if self.iter_count == 0:  # record the initial medoids
-            self.history_medoids = [[m.copy()] for m in self.active_indicator._medoids]
+            self.history_medoids = [[m.copy()] for m in self._igd._medoids]
         else:
             # log the updated medoids
             for i, k in enumerate(indices):
-                self.history_medoids[k].append(self.active_indicator._medoids[indices][i])
+                self.history_medoids[k].append(self._igd._medoids[indices][i])
         self.logger.info(f"{len(indices)} target points are shifted")
 
     def _backtracking_line_search(

hvd/problems/base.py

@@ -4,6 +4,7 @@
 # enable double-precision of JAX
 os.environ["JAX_ENABLE_X64"] = "True"
 
+import jax
 import jax.numpy as jnp
 import numpy as np
 from jax import jacfwd, jacrev, jit
@@ -26,10 +27,9 @@ def func(x):
 
 class MOOAnalytical:
     def __init__(self):
-        obj_func = partial(self.__class__._objective, self)
-        self._obj_func = jit(obj_func)
-        self._objective_jacobian = jit(jacrev(obj_func))
-        self._objective_hessian = jit(hessian(obj_func))
+        self._obj_func = jit(partial(self.__class__._objective, self))
+        self._objective_jacobian = jit(jacrev(self._obj_func))
+        self._objective_hessian = jit(hessian(self._obj_func))
         self.CPU_time: int = 0  # in nanoseconds
 
     def objective(self, x: np.ndarray) -> np.ndarray:
@@ -154,6 +154,11 @@ def get_pareto_front(self, N: int = 1000) -> np.ndarray:
 
 
 class CONV4(MOOAnalytical):
+    pass
+
+
+class CONV42F(MOOAnalytical):
+    """Convex Problem 4 with 2 disconnected Pareto fronts"""
     def __init__(self):
         self.n_obj = 4
         self.n_var = 4
@@ -168,15 +173,15 @@ def _objective(self, x: jnp.ndarray) -> jnp.ndarray:
         fa4 = jnp.array([2, 2, 2, 0])
         fa1 = jnp.array([0, 2, 2, 2])
         deltay = fa4 - fa1
-
-        if jnp.all(x < 0):
-            z = x + deltaa
-            y = jnp.array([jnp.sum((z - a[i]) ** 2) - 1.1 * deltay[i] for i in range(4)])
-        else:
-            y = jnp.array([jnp.sum((x - a[i]) ** 2) for i in range(4)])
+        z = x + deltaa
+        y = jax.lax.select(jnp.all(x < 0), jnp.array([jnp.sum((z - a[i]) ** 2) - 1.1 * deltay[i] for i in range(4)]), jnp.array([jnp.sum((x - a[i]) ** 2) for i in range(4)]))
+        # if jnp.all(x < 0):
+        #     z = x + deltaa
+        #     y = jnp.array([jnp.sum((z - a[i]) ** 2) - 1.1 * deltay[i] for i in range(4)])
+        # else:
+        #     y = jnp.array([jnp.sum((x - a[i]) ** 2) for i in range(4)])
         return y
 
-
 class UF7(MOOAnalytical):
     def __init__(self, n_var: int = 30) -> None:
         self.n_obj = 2

hvd/problems/cf.py

@@ -365,7 +365,6 @@ def get_pareto_front(self, N: int = 1000) -> np.ndarray:
 
 class CF9(CF8, ConstrainedMOOAnalytical):
     def __init__(self, n_var: int = 10, boundry_constraints: bool = False) -> None:
-        self.n_var = n_var
         super().__init__(n_var, boundry_constraints)
         self.xl = np.r_[0, 0, np.zeros(self.n_var - 2) - 2]
         self.xu = np.r_[1, 1, np.zeros(self.n_var - 2) + 2]

hvd/problems/eqdtlz.py

@@ -4,15 +4,16 @@
 from ..utils import timeit
 from .base import ConstrainedMOOAnalytical, MOOAnalytical
 
+# TODO: unify the DTLZ classes here with the ones in `dtlz.py`
 
 class _DTLZ(MOOAnalytical):
-    def __init__(self, n_obj: int = 3, n_var: int = 11):
+    def __init__(self, n_obj: int = 3, n_var: int = 11, boundry_constraints: bool = False):
         self.n_obj = n_obj
         # the default decision space of 11-dimensional
         self.n_var = n_var if n_var is not None else self.n_obj + 8
         self.xl = np.zeros(self.n_var)
         self.xu = np.ones(self.n_var)
-        super().__init__()
+        super().__init__(boundry_constraints=boundry_constraints)
 
     def get_pareto_set(self, N: int = 1000, kind="uniform") -> np.ndarray:
         M = self.n_obj
@@ -124,7 +125,7 @@ def _objective(self, x: jnp.ndarray) -> jnp.ndarray:
             [1], jnp.sin(x_[::-1] * jnp.pi / 2)
         ]
 
-
+# TODO: somehow simplify the way of defining the following classes since the constraints are the same
 class Eq1DTLZ1(DTLZ1):
     n_eq_constr = 1