update the benchmark scripts

hvd/newton.py

@@ -600,7 +600,6 @@ def __init__(
         g: Callable = None,
         g_jac: Callable = None,
         x0: np.ndarray = None,
-        y0: np.ndarray = None,
         max_iters: Union[int, str] = np.inf,
         xtol: float = 0,
         verbose: bool = True,
@@ -650,7 +649,6 @@ def __init__(
         self._initialize(x0)
         self._set_indicator(ref, func, jac, hessian)
         self._set_logging(verbose)
-        self.y0 = y0
         # parameters controlling stop criteria
         self.xtol = xtol
         self.max_iters: int = self.N * 10 if max_iters is None else max_iters
@@ -763,11 +761,7 @@ def run(self) -> Tuple[np.ndarray, np.ndarray, Dict]:
 
     def newton_iteration(self):
         # compute the initial indicator values. The first clustering and matching is executed here.
-        if self.iter_count == 0:
-            # TODO: use `self.state.Y` instead
-            self._compute_indicator_value(self.y0)
-        else:
-            self._compute_indicator_value(self.state.Y)
+        self._compute_indicator_value(self.state.Y)
         # shift the reference set if needed
         self._shift_reference_set()
         # compute the Newton step
@@ -887,19 +881,6 @@ def _shift_reference_set(self):
                 np.isclose(np.linalg.norm(self.step[:, : self.dim_p], axis=1), 0),
             )
         )
-        if 11 < 2:
-            import matplotlib.pyplot as plt
-
-            fig, ax = plt.subplots(1, 1, figsize=(10, 6.5))
-            x = self.active_indicator._medoids
-            y = self.state.Y
-            ax.plot(x[:, 0], x[:, 1], "r+")
-            ax.plot(y[:, 0], y[:, 1], "k+")
-            for i in range(len(x)):
-                ax.plot((x[i, 0], y[i, 0]), (x[i, 1], y[i, 1]), "k--")
-            plt.savefig(f"{self.iter_count}.pdf")
-            plt.close(fig)
-
         indices = np.nonzero(masks)[0]
         if len(indices) == 0:
             return

hvd/utils.py

@@ -43,6 +43,20 @@
 #     return step, max_step_size
 
 
+# if 11 < 2:
+#     import matplotlib.pyplot as plt
+
+#     fig, ax = plt.subplots(1, 1, figsize=(10, 6.5))
+#     x = self.active_indicator._medoids
+#     y = self.state.Y
+#     ax.plot(x[:, 0], x[:, 1], "r+")
+#     ax.plot(y[:, 0], y[:, 1], "k+")
+#     for i in range(len(x)):
+#         ax.plot((x[i, 0], y[i, 0]), (x[i, 1], y[i, 1]), "k--")
+#     plt.savefig(f"{self.iter_count}.pdf")
+#     plt.close(fig)
+
+
 def merge_lists(x, y):
     if x is None:
         return y

scripts/benchmark_CF.py

@@ -40,7 +40,7 @@
 pareto_front = problem.get_pareto_front(1000)
 
 path = "CF_gen_300"
-emoa = "NSGA-II"
+emoa = "NSGA-III"
 gen = 300
 
 
@@ -158,7 +158,6 @@ def execute(run: int):
     x0 = x0[idx]
     y0 = y0[idx]
     Y_label = Y_label[idx]
-    y0 = np.array([problem.objective(_) for _ in x0])
 
     # TODO: this is an ad-hoc solution. Maybe fix this special case in the `ReferenceSet` class
     # if the minimal number of points in the `ref` clusters is smaller than
@@ -183,7 +182,6 @@ def execute(run: int):
         g_jac=problem.ieq_jacobian,
         N=len(x0),
         x0=x0,
-        y0=y0,
         xl=problem.xl,
         xu=problem.xu,
         max_iters=max_iters,

scripts/benchmark_DTLZ.py

@@ -15,6 +15,7 @@
 from hvd.hypervolume import hypervolume
 from hvd.newton import DpN
 from hvd.problems import DTLZ1, DTLZ2, DTLZ3, DTLZ4, DTLZ5, DTLZ6, DTLZ7, PymooProblemWithAD
+from hvd.utils import get_non_dominated
 
 plt.style.use("ggplot")
 rcParams["font.size"] = 17
@@ -33,16 +34,16 @@
 
 max_iters = 6
 n_jobs = 30
-ref_point = np.array([11, 11])
 problem_name = sys.argv[1]
 print(problem_name)
 f = locals()[problem_name]()
 problem = PymooProblemWithAD(f)
 pareto_front = problem.get_pareto_front()
 
-gen = 100
-path = f"./DTLZ-gen{gen}/"
-emoa = "NSGA-III"
+gen = 300
+# path = f"./DTLZ-gen{gen}/"
+path = "./DTLZ_refs/"
+emoa = "NSGA-II"
 
 
 def plot(y0, Y, ref, hist_Y, history_medoids, hist_IGD, hist_R_norm, fig_name):
@@ -164,12 +165,17 @@ def execute(run: int):
     x0 = x0[idx]
     y0 = y0[idx]
     Y_label = Y_label[idx]
+
+    # TODO: this is an ad-hoc solution. Maybe fix this special case in the `ReferenceSet` class
+    # if the minimal number of points in the `ref` clusters is smaller than
+    # the maximal number of points in `y0` clusters, then we merge all clusters
+    min_point_ref_cluster = np.min([len(r) for r in ref.values()])
+    max_point_y_cluster = np.max(np.unique(Y_label, return_counts=True)[1])
     # if the number of clusters of `Y` is more than that of the reference set
-    if len(np.unique(Y_label)) > len(ref):
+    if (len(np.unique(Y_label)) > len(ref)) or (max_point_y_cluster > min_point_ref_cluster):
         ref = np.vstack([r for r in ref.values()])
         Y_label = np.zeros(len(y0))
         eta = None
-
     N = len(x0)
     opt = DpN(
         dim=problem.n_var,
@@ -182,6 +188,7 @@ def execute(run: int):
         g_jac=problem.ieq_jacobian,
         N=len(x0),
         x0=x0,
+        # y0=y0,
         xl=problem.xl,
         xu=problem.xu,
         max_iters=max_iters,
@@ -192,29 +199,31 @@ def execute(run: int):
         Y_label=Y_label,
     )
     X, Y, _ = opt.run()
-    data = np.concatenate([np.c_[[0] * N, y0], np.c_[[max_iters] * N, opt.hist_Y[-1]]], axis=0)
-    df = pd.DataFrame(data, columns=["iteration", "f1", "f2", "f3"])
-    df.to_csv(f"./tmp/{problem_name}_DpN_{emoa}_run{run}_{gen}.csv")
+    # remove the dominated solution in Y
+    Y = get_non_dominated(Y)
+    # data = np.concatenate([np.c_[[0] * N, y0], np.c_[[max_iters] * N, opt.hist_Y[-1]]], axis=0)
+    # df = pd.DataFrame(data, columns=["iteration", "f1", "f2", "f3"])
+    # df.to_csv(f"./tmp/{problem_name}_DpN_{emoa}_run{run}_{gen}.csv")
     # fig_name = f"./figure/{problem_name}_DpN_{emoa}_run{run}_{gen}.pdf"
     # plot(y0, Y, all_ref, opt.hist_Y, opt.history_medoids, opt.hist_IGD, opt.hist_R_norm, fig_name)
     gd_value = GenerationalDistance(pareto_front).compute(Y=Y)
     igd_value = InvertedGenerationalDistance(pareto_front).compute(Y=Y)
-    return np.array([igd_value, gd_value])
+    return np.array([igd_value, gd_value, opt.state.n_jac_evals])
 
 
 # get all run IDs
 run_id = [
     int(re.findall(r"run_(\d+)_", s)[0])
     for s in glob(f"{path}/{problem_name}_{emoa}_run_*_lastpopu_x_gen{gen}.csv")
 ]
-if gen == 110 and problem_name == "DTLZ4":
-    run_id = list(set(run_id) - set([14]))
-
-if 1 < 2:
+print(run_id)
+if 11 < 2:
+    data = []
     for i in run_id:
-        execute(i)
+        print(i)
+        data.append(execute(i))
 else:
     data = Parallel(n_jobs=n_jobs)(delayed(execute)(run=i) for i in run_id)
-    # df = pd.DataFrame(np.array(data), columns=["IGD", "GD", "HV"])
-    df = pd.DataFrame(np.array(data), columns=["IGD", "GD"])
-    df.to_csv(f"{problem_name}-DpN-{emoa}-{gen}.csv", index=False)
+
+df = pd.DataFrame(np.array(data), columns=["IGD", "GD", "Jac_calls"])
+df.to_csv(f"results/{problem_name}-DpN-{emoa}-{gen}.csv", index=False)

scripts/benchmark_EA.py

@@ -28,6 +28,8 @@
 
 
 class ProblemWrapper(ElementwiseProblem):
+    """Wrap a `MOOAnalytical` problem into a problem in `Pymoo"""
+
     def __init__(self, problem: MOOAnalytical) -> None:
         self._problem = problem
         super().__init__(
@@ -44,10 +46,10 @@ def pareto_front(self, N: int) -> np.ndarray:
 
     def _evaluate(self, x: np.ndarray, out: dict, *args, **kwargs) -> None:
         out["F"] = self._problem.objective(x)  # objective value
-        if self.n_eq_constr > 0:
-            out["H"] = self._problem.eq_constraint(x)  # equality constraint value
-        if self.n_ieq_constr > 0:
-            out["G"] = self._problem.ieq_constraint(x)  # inequality constraint value
+        if self._problem.n_eq_constr > 0:
+            out["H"] = np.array([self._problem.eq_constraint(_) for _ in x])  # equality constraint value
+        if self._problem.n_ieq_constr > 0:
+            out["G"] = np.array([self._problem.ieq_constraint(_) for _ in x])  # inequality constraint value
 
 
 class ModifiedObjective(Problem):
@@ -107,18 +109,20 @@ def minimize(
 
     # store the deep copied algorithm in the result object
     res.algorithm = algorithm
-    if problem.n_obj == 3:
-        if problem_name not in ["DTLZ4", "DTLZ6", "DTLZ7"]:
+    if problem.name().startswith("DTLZ"):
+        if problem_name in ["DTLZ5", "DTLZ6", "DTLZ7"]:
+            pareto_front = problem.pareto_front()
+        else:
             ref_dirs = UniformReferenceDirectionFactory(3, n_partitions=30).do()
             pareto_front = problem.pareto_front(ref_dirs)
-        else:
-            pareto_front = problem.pareto_front()
     else:
         pareto_front = problem.pareto_front(1000)
-    gd_value = GenerationalDistance(pareto_front).compute(Y=res.F)
-    igd_value = InvertedGenerationalDistance(pareto_front).compute(Y=res.F)
-    # return np.array([igd_value, gd_value, hypervolume(res.F, ref_point)])
-    # return np.array([igd_value, gd_value, np.sum(np.bitwise_or(res.X < 0, res.X > 1))])
+    # TODO: ad-hoc solution for `res.F` being `None`. Figure out why..
+    if res.F is None:
+        igd_value, gd_value = np.nan, np.nan
+    else:
+        gd_value = GenerationalDistance(pareto_front).compute(Y=res.F)
+        igd_value = InvertedGenerationalDistance(pareto_front).compute(Y=res.F)
     return np.array([igd_value, gd_value])
 
 
@@ -156,8 +160,9 @@ def get_Jacobian_calls(path, problem_name, algorithm_name, gen):
     return int(np.median(pd.read_csv(f"{path}/{problem_name}-DpN-{algorithm_name}-{gen}.csv").Jac_calls))
 
 
-n_iter_newton = 8
+n_iter_newton = 6
 gen = 300
+# NOTE: the following running budget is estimated with upper bounds of AD's theory
 # gen_func = lambda n_var, scale: 4 * scale + 10 * n_var
 # NOTE: 1.836 is obtained on ZDTs
 gen_func = lambda n_var, scale: int(1.836 * scale + 3)
@@ -166,8 +171,8 @@ def get_Jacobian_calls(path, problem_name, algorithm_name, gen):
 
 for problem_name in [problem]:
     print(problem_name)
-    # problem = get_problem(problem_name)
-    problem = ProblemWrapper(locals()[problem_name]())
+    problem = get_problem(problem_name)
+    # problem = ProblemWrapper(locals()[problem_name]())
     pop_size = 100 if problem.n_obj == 2 else 300
     constrained = problem.n_eq_constr > 0 or problem.n_ieq_constr > 0
 
@@ -183,7 +188,9 @@ def get_Jacobian_calls(path, problem_name, algorithm_name, gen):
             for i in range(N)
         )
         # df = pd.DataFrame(np.array(data), columns=["IGD", "GD", "HV"])
-        df = pd.DataFrame(np.array(data), columns=["IGD", "GD"])
+        data = np.array(data)
+        data = data[~np.any(np.isnan(data), axis=1)]
+        df = pd.DataFrame(data, columns=["IGD", "GD"])
         df.to_csv(f"results/{problem_name}-{algorithm_name}-{gen}.csv", index=False)
         # data = pd.concat(data, axis=0)
         # data.to_csv(f"./data/{problem_name.upper()}_{algorithm_name}.csv", index=False)

scripts/benchmark_ZDT.py

@@ -177,7 +177,6 @@ def execute(run: int):
         g_jac=problem.ieq_jacobian,
         N=len(x0),
         x0=x0,
-        y0=y0,
         xl=problem.xl,
         xu=problem.xu,
         max_iters=max_iters,