update

.gitignore

@@ -15,6 +15,7 @@
 data/
 model/
 hv/
+RSG/
 *gen*/
 
 **/*.pyc

examples/HVN/dent_example.py

@@ -28,9 +28,11 @@
 
 ref = np.array([5, 5])
 max_iters = 10
-X0 = pd.read_csv("dent/dent_x.csv", header=None, index_col=None).values
-Y0 = pd.read_csv("dent/dent_y.csv", header=None, index_col=None).values
+# X0 = pd.read_csv("dent/dent_x.csv", header=None, index_col=None).values
+# Y0 = pd.read_csv("dent/dent_y.csv", header=None, index_col=None).values
 problem = DENT()
+X0 = problem.get_pareto_set(50)
+Y0 = np.array([problem.objective(x) for x in X0])
 pareto_front = problem.get_pareto_front(1000)
 N = len(X0)
 

examples/MMD/example_DTLZ1.py

@@ -7,10 +7,11 @@
 import pandas as pd
 from matplotlib import rcParams
 
+from hvd.bootstrap import bootstrap_reference_set
 from hvd.delta_p import GenerationalDistance, InvertedGenerationalDistance
-from hvd.mmd_newton import MMDNewton, bootstrap_reference_set
+from hvd.mmd_newton import MMDNewton
 from hvd.newton import DpN
-from hvd.problems import DTLZ1, PymooProblemWithAD
+from hvd.problems import DTLZ1
 from hvd.reference_set import ClusteredReferenceSet
 from hvd.utils import get_non_dominated
 
@@ -30,15 +31,14 @@
 np.random.seed(66)
 max_iters = 16
 
-f = DTLZ1(boundry_constraints=True)
-problem = PymooProblemWithAD(f)
+problem = DTLZ1(n_var=7, boundry_constraints=True)
 pareto_front = problem.get_pareto_front()
 
 # read the reference set data
-ref_ = pd.read_csv("./DTLZ1/DTLZ1_RANDOM_run_1_ref_1_gen0.csv", header=None).values
-X0 = pd.read_csv("./DTLZ1/DTLZ1_RANDOM_run_1_lastpopu_x_gen0.csv", header=None).values
-Y0 = pd.read_csv("./DTLZ1/DTLZ1_RANDOM_run_1_lastpopu_y_gen0.csv", header=None).values
-eta = {0: pd.read_csv("./DTLZ1/DTLZ1_RANDOM_run_1_eta_1_gen0.csv", header=None).values.ravel()}
+ref_ = pd.read_csv("./data/DTLZ1/DTLZ1_RANDOM_run_1_ref_1_gen0.csv", header=None).values
+X0 = pd.read_csv("./data/DTLZ1/DTLZ1_RANDOM_run_1_lastpopu_x_gen0.csv", header=None).values
+Y0 = pd.read_csv("./data/DTLZ1/DTLZ1_RANDOM_run_1_lastpopu_y_gen0.csv", header=None).values
+eta = {0: pd.read_csv("./data/DTLZ1/DTLZ1_RANDOM_run_1_eta_1_gen0.csv", header=None).values.ravel()}
 Y_idx = None
 Y_label = np.array([0] * len(ref_))
 ref = {0: ref_}
@@ -67,11 +67,11 @@
 if 1 < 2:
     opt.indicator.beta = 0.25
     # TODO: figure out how to determine when to bootstrap automatically
-    X, Y, _ = bootstrap_reference_set(opt, problem, 7)
+    X, Y, _ = bootstrap_reference_set(opt, problem, 5)
 else:
     X, Y, _ = opt.run()
 
-Y = get_non_dominated(Y)
+# Y = get_non_dominated(Y)
 igd_mmd = igd.compute(Y=Y)
 
 opt_dpn = DpN(
@@ -145,5 +145,4 @@
 ax1.set_ylabel(r"$f_3$")
 
 plt.tight_layout()
-# plt.show()
-plt.savefig(f"MMD-{f.__class__.__name__}.pdf", dpi=1000)
+plt.savefig(f"MMD-{problem.__class__.__name__}.pdf", dpi=1000)

examples/subset_selection/Dent/example.py

@@ -0,0 +1,139 @@
+import sys
+
+sys.path.insert(0, "./")
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from matplotlib import rcParams
+
+from examples.subset_selection.Dent.spline import pareto_front_approx
+from hvd.newton import HVN
+
+plt.style.use("ggplot")
+plt.rc("text.latex", preamble=r"\usepackage{amsmath}")
+rcParams["font.size"] = 11
+rcParams["xtick.direction"] = "out"
+rcParams["ytick.direction"] = "out"
+rcParams["text.usetex"] = True
+rcParams["legend.numpoints"] = 1
+rcParams["xtick.labelsize"] = 11
+rcParams["ytick.labelsize"] = 11
+rcParams["xtick.major.size"] = 7
+rcParams["xtick.major.width"] = 1
+rcParams["ytick.major.size"] = 7
+rcParams["ytick.major.width"] = 1
+
+np.random.seed(66)
+
+ref = np.array([5, 5])
+max_iters = 25
+X0 = pd.read_csv("examples/subset_selection/Dent/points.csv", header=None, index_col=None).values
+X0 = X0[X0[:, 1].argsort()]
+Y0 = X0 = X0[1:-1, :]
+N = len(X0)
+
+
+def objective(x: np.ndarray) -> np.ndarray:
+    return x
+
+
+def jacobian(_: np.ndarray) -> np.ndarray:
+    return np.eye(2)
+
+
+def hessian(_) -> np.ndarray:
+    return np.array([np.zeros((2, 2))] * 2)
+
+
+def h(x: np.ndarray) -> float:
+    return pareto_front_approx(x[0])[0] - x[1]
+
+
+def h_Jacobian(x: np.ndarray) -> np.ndarray:
+    return np.array([pareto_front_approx(x[0])[1], -1])
+
+
+def h_Hessian(x: np.ndarray) -> np.ndarray:
+    return np.array([[pareto_front_approx(x[0])[2], 0], [0, 0]])
+
+
+opt = HVN(
+    n_var=2,
+    n_obj=2,
+    ref=ref,
+    func=objective,
+    jac=jacobian,
+    hessian=hessian,
+    h=h,
+    h_jac=h_Jacobian,
+    h_hessian=h_Hessian,
+    N=len(X0),
+    X0=X0,
+    xl=[0.561606693680458, 0.5616],
+    xu=[4.55710916559211, 4.5582],
+    max_iters=max_iters,
+    verbose=True,
+    preconditioning=False,
+)
+X, Y, stop = opt.run()
+
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(14, 6.5))
+plt.subplots_adjust(right=0.93, left=0.05)
+ax0.set_aspect("equal")
+# ax0.plot(Y[:, 0], Y[:, 1], "r+", ms=8)
+ax0.plot(Y0[:, 0], Y0[:, 1], "k+", ms=8)
+
+if 11 < 2:
+    trajectory = np.array([Y0] + opt.history_Y)
+    for i in range(N):
+        x, y = trajectory[:, i, 0], trajectory[:, i, 1]
+        ax0.quiver(
+            x[:-1],
+            y[:-1],
+            x[1:] - x[:-1],
+            y[1:] - y[:-1],
+            scale_units="xy",
+            angles="xy",
+            scale=1,
+            color="k",
+            width=0.005,
+            alpha=0.5,
+            headlength=4.7,
+            headwidth=2.7,
+        )
+
+x = np.linspace(0.561606693680458, 4.55710916559211, 1000)
+y = np.array([pareto_front_approx(_)[0] for _ in x])
+pareto_front = np.c_[x, y]
+ax0.plot(pareto_front[:, 0], pareto_front[:, 1], "k--", alpha=0.5)
+ax0.set_title("Objective space")
+ax0.set_xlabel(r"$f_1$")
+ax0.set_ylabel(r"$f_2$")
+ax0.legend([r"$Y_0$", "Approximated Pareto front"])
+
+
+ax1.plot(Y[:, 0], Y[:, 1], "r+", ms=8)
+ax1.plot(pareto_front[:, 0], pareto_front[:, 1], "k--", alpha=0.5)
+ax1.set_title("Objective space")
+ax1.set_xlabel(r"$f_1$")
+ax1.set_ylabel(r"$f_2$")
+ax1.legend([r"$Y_{\text{final}}$", "Approximated Pareto front"])
+
+# ax22 = ax1.twinx()
+# HV0 = opt.history_indicator_value[0]
+# impr = (np.array(opt.history_indicator_value) - HV0) / HV0
+# lns = ax1.plot(range(1, len(opt.history_indicator_value) + 1), opt.history_indicator_value, "b-")
+# lns += ax22.semilogy(range(1, len(opt.history_R_norm) + 1), opt.history_R_norm, "g--")
+# ax1.legend(lns, ["HV", r"$||R(\mathbf{X})||$"], loc=1)
+# ax1.set_ylabel("HV", color="b")
+# ax22.set_ylabel(r"$||R(\mathbf{X})||$", color="g")
+# ax1.set_title(f"Performance with ref: {ref}")
+# ax1.set_xlabel("iteration")
+# ax1.set_xticks(range(1, max_iters + 1))
+plt.tight_layout()
+
+plt.savefig(f"dent-example-{N}.pdf", dpi=1000)
+
+data = pd.DataFrame(np.c_[Y0, Y], columns=["initial y1", "initial y2", "final y1", "final y2"])
+data.to_csv(f"dent-example-{N}.csv")

examples/subset_selection/Dent/spline.py

@@ -0,0 +1,26 @@
+from typing import Tuple
+
+import numpy as np
+import pandas as pd
+
+coefs = pd.read_csv("examples/subset_selection/Dent/coeff.csv", header=None).values
+breaks = pd.read_csv("examples/subset_selection/Dent/knots.csv", header=None).values.ravel()
+
+
+def pareto_front_approx(x: float) -> Tuple[float, float, float]:
+    interval = np.where(breaks <= x)[0][-1]
+    interval_min = breaks[interval]
+    interval_coefs = coefs[interval, :]
+    y = (
+        interval_coefs[0] * (x - interval_min) ** 3
+        + interval_coefs[1] * (x - interval_min) ** 2
+        + interval_coefs[2] * (x - interval_min)
+        + interval_coefs[3]
+    )
+    dy = (
+        3 * interval_coefs[0] * (x - interval_min) ** 2
+        + 2 * interval_coefs[1] * (x - interval_min)
+        + interval_coefs[2]
+    )
+    d2y = 6 * interval_coefs[0] * (x - interval_min) + 2 * interval_coefs[1]
+    return y, dy, d2y