update

examples/HVN/2D_example.py

@@ -93,7 +93,7 @@ def h_Hessian(x):
 )
 X, Y, stop = opt.run()
 
-fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(22, 5))
+fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(22, 8))
 plt.subplots_adjust(right=0.93, left=0.05)
 ciricle = plt.Circle((0, 0), 1, color="r", fill=False, ls="--", lw=1.5)
 

examples/MMD/example_DTLZ1.py

@@ -29,7 +29,7 @@
 rcParams["ytick.major.width"] = 1
 
 np.random.seed(66)
-max_iters = 10
+max_iters = 13
 f = DTLZ1(boundry_constraints=True)
 problem = PymooProblemWithAD(f)
 pareto_front = problem.get_pareto_front()
@@ -67,8 +67,8 @@
     metrics=metrics,
     preconditioning=False,
 )
-if 11 < 2:
-    X, Y, _ = bootstrap_reference_set(opt, problem, ref_, 5)
+if 1 < 2:
+    X, Y, _ = bootstrap_reference_set(opt, problem, 5)
 else:
     X, Y, _ = opt.run()
 Y = get_non_dominated(Y)

hvd/base.py

@@ -38,6 +38,7 @@ def __init__(
 
     @property
     def N(self):
+        # TODO: maybe `N` -> `__len__`
         return len(self.X)
 
     def jac(self, x: np.ndarray) -> np.ndarray:
@@ -63,15 +64,17 @@ def evaluate_cstr_jac(self, x: np.ndarray) -> np.ndarray:
         # 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)
+        out = np.array(dH + dG)
+        return out if len(out) == 0 else out.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)
+        out = np.array(ddH + ddG)
+        return out if len(out) == 0 else out.reshape(self.n_cstr, self.n_var, self.n_var)
 
     def update(self, X: np.ndarray, compute_gradient: bool = True):
         self.X = X

hvd/mmd_newton.py

@@ -370,9 +370,7 @@ def _handle_box_constraint(self, step: np.ndarray) -> Tuple[np.ndarray, np.ndarr
         return step, max_step_size
 
 
-def bootstrap_reference_set(
-    optimizer, problem, init_ref: np.ndarray, interval: int = 5
-) -> Tuple[np.ndarray, np.ndarray, Dict]:
+def bootstrap_reference_set(optimizer, problem, interval: int = 5) -> Tuple[np.ndarray, np.ndarray, Dict]:
     """Bootstrap the reference set with the intermediate population of an MOO algorithm
 
     Args:
@@ -385,18 +383,26 @@ def bootstrap_reference_set(
         Tuple[np.ndarray, np.ndarray, Dict]: (efficient set, Pareto front approximation,
             the stopping criteria)
     """
+    alpha = 0.05
     for i in range(optimizer.max_iters):
+        # TODO: generalize the condition to trigger the boostrap: maybe if the R norm for most points
+        # are near zero?
         if i % interval == 0 and i > 0:
-            # TODO: find a better way to take out outliersin `state.Y`
-            ref_ = np.r_[optimizer.state.Y, init_ref]
-            ref_ = get_non_dominated(ref_)
+            # relax dominance resistent points (improper dominance)
+            func = lambda y: (1 - alpha) * y + alpha * y.mean(axis=1).reshape(len(y), -1)
+            Y = func(optimizer.state.Y)
+            idx = get_non_dominated(Y, return_index=True)
+            ref_ = Y[idx]
             eta = compute_chim(ref_)
             ref_ += 0.08 * eta
             Y_idx = None
             ref = ClusteredReferenceSet(ref=ref_, eta={0: eta}, Y_idx=Y_idx)
-            optimizer.ref = ref
-            optimizer.indicator.ref = ref
-            optimizer.indicator.compute(Y=optimizer.state.Y)  # To compute the medoids
+            # only keep the non-dominated points
+            optimizer.N = len(idx)
+            optimizer.step = optimizer.step[idx]
+            optimizer.state = optimizer.state[idx]
+            optimizer.indicator.ref = optimizer.ref = ref
+            optimizer.indicator.compute(Y=optimizer.state.Y)  # to compute the medoids
 
             if 11 < 2:
                 import matplotlib.pyplot as plt
@@ -408,9 +414,17 @@ def bootstrap_reference_set(
                 ax0 = fig.add_subplot(1, 1, 1, projection="3d")
                 ax0.set_box_aspect((1, 1, 1))
                 ax0.view_init(45, 45)
+                ax0.plot(
+                    optimizer.state.Y[:, 0],
+                    optimizer.state.Y[:, 1],
+                    optimizer.state.Y[:, 2],
+                    "r+",
+                    ms=6,
+                    alpha=0.6,
+                )
                 ax0.plot(ref_[:, 0], ref_[:, 1], ref_[:, 2], "g.", ms=6, alpha=0.6)
-                ax0.plot(pareto_front[:, 0], pareto_front[:, 1], pareto_front[:, 2], "k.", ms=6, alpha=0.6)
-                ax0.plot(m[:, 0], m[:, 1], m[:, 2], "r+", ms=6, alpha=0.6)
+                ax0.plot(pareto_front[:, 0], pareto_front[:, 1], pareto_front[:, 2], "k.", ms=6, alpha=0.3)
+                ax0.plot(m[:, 0], m[:, 1], m[:, 2], "r^", ms=6, alpha=0.6)
                 plt.show()
 
         optimizer.newton_iteration()