add DENT example

examples/HVN/concave_example.py

@@ -4,6 +4,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 from matplotlib import rcParams
 
 from hvd.newton import HVN
@@ -51,39 +52,39 @@ def g_Hessian(x):
 
 ref = np.array([1, 1])
 max_iters = 10
-N = 10
-x1 = np.random.rand(N)
-x2 = 1 - x1**2
+x1 = pd.read_csv("convex_concave/concave_x.csv", header=None, index_col=None).values
+x2 = pd.read_csv("convex_concave/concave_y.csv", header=None, index_col=None).values
 x0 = np.c_[x1, x2]
 y0 = x0.copy()
-
+N = len(x0)
 opt = HVN(
     n_var=2,
     n_obj=2,
     ref=ref,
     func=concave,
     jac=concave_Jacobian,
     hessian=concave_Hessian,
-    h=g,
-    h_jac=g_Jacobian,
-    h_hessian=g_Hessian,
+    g=g,
+    g_jac=g_Jacobian,
+    g_hessian=g_Hessian,
     N=len(x0),
     X0=x0,
     xl=0,
     xu=1,
     max_iters=max_iters,
     verbose=True,
+    preconditioning=True,
 )
 X, Y, stop = opt.run()
 
 fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(13, 6.5))
 plt.subplots_adjust(right=0.93, left=0.05)
 ax0.set_aspect("equal")
-ax0.plot(Y[:, 0], Y[:, 1], "r*")
+ax0.plot(Y[:, 0], Y[:, 1], "r*", ms=8)
 ax0.plot(y0[:, 0], y0[:, 1], "k+", ms=8)
-trajectory = np.array([y0] + opt.history_Y)
 
 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(
@@ -103,18 +104,19 @@ def g_Hessian(x):
 
 x_vals = np.linspace(0, 1, 100)
 y_vals = 1 - x_vals**2
-ax0.plot(x_vals, y_vals, "r--")
+ax0.plot(x_vals, y_vals, "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_{\text{final}}$", r"$Y_0$", "Pareto front"])
 
-ax12 = ax1.twinx()
+ax22 = ax1.twinx()
 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--")
+ax22.semilogy(range(1, len(opt.history_R_norm) + 1), opt.history_R_norm, "g--")
 ax1.set_ylabel("HV", color="b")
-ax12.set_ylabel(r"$||R(\mathbf{X})||$", color="g")
+ax22.set_ylabel(r"$||R(\mathbf{X})||$", color="g")
 ax1.set_title("Performance")
 ax1.set_xlabel("iteration")
 ax1.set_xticks(range(1, max_iters + 1))
+
 plt.savefig(f"concave-example-{N}.pdf", dpi=1000)

examples/HVN/convex_example.py

@@ -3,6 +3,7 @@
 sys.path.insert(0, "./")
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 from matplotlib import rcParams
 
 from hvd.newton import HVN
@@ -52,12 +53,11 @@ def g_Hessian(x):
 ref = np.array([1, 1])
 max_iters = 10
 
-N = 15
-x1 = np.random.rand(N)
-x2 = x1**2 - 2 * x1 + 1
-x0 = np.c_[x1, x2]
-y0 = x0.copy()
-
+x1 = pd.read_csv("convex_concave/convex_x.csv", header=None, index_col=None).values
+x2 = pd.read_csv("convex_concave/convex_y.csv", header=None, index_col=None).values
+X0 = np.c_[x1, x2]
+y0 = X0.copy()
+N = len(X0)
 opt = HVN(
     n_var=2,
     n_obj=2,
@@ -68,8 +68,8 @@ def g_Hessian(x):
     g=g,
     g_jac=g_Jacobian,
     g_hessian=g_Hessian,
-    N=len(x0),
-    X0=x0,
+    N=len(X0),
+    X0=X0,
     xl=0,
     xu=1,
     max_iters=max_iters,
@@ -81,11 +81,11 @@ def g_Hessian(x):
 fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(13, 6.5))
 plt.subplots_adjust(right=0.93, left=0.05)
 ax0.set_aspect("equal")
-ax0.plot(Y[:, 0], Y[:, 1], "r*")
+ax0.plot(Y[:, 0], Y[:, 1], "r*", ms=6)
 ax0.plot(y0[:, 0], y0[:, 1], "k+", ms=8)
-trajectory = np.array([y0] + opt.history_Y)
 
 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(
@@ -105,18 +105,19 @@ def g_Hessian(x):
 
 x_vals = np.linspace(0, 1, 100)
 y_vals = x_vals**2 - 2 * x_vals + 1
-ax0.plot(x_vals, y_vals, "r--")
+ax0.plot(x_vals, y_vals, "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_{\text{final}}$", r"$Y_0$", "Pareto front"])
 
-ax12 = ax1.twinx()
+ax22 = ax1.twinx()
 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--")
+ax22.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")
+ax22.set_ylabel(r"$||G(\mathbf{X})||$", color="g")
 ax1.set_title("Performance")
 ax1.set_xlabel("iteration")
 ax1.set_xticks(range(1, max_iters + 1))
+
 plt.savefig(f"convex-example-{N}.pdf", dpi=1000)

examples/HVN/dent_example.py

@@ -0,0 +1,94 @@
+import sys
+
+sys.path.insert(0, "./")
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from matplotlib import rcParams
+
+from hvd.newton import HVN
+from hvd.problems.misc import DENT
+
+plt.style.use("ggplot")
+plt.rc("text.latex", preamble=r"\usepackage{amsmath}")
+rcParams["font.size"] = 13
+rcParams["xtick.direction"] = "out"
+rcParams["ytick.direction"] = "out"
+rcParams["text.usetex"] = True
+rcParams["legend.numpoints"] = 1
+rcParams["xtick.labelsize"] = 13
+rcParams["ytick.labelsize"] = 13
+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 = 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
+problem = DENT()
+pareto_front = problem.get_pareto_front(1000)
+N = len(X0)
+
+opt = HVN(
+    n_var=2,
+    n_obj=2,
+    ref=ref,
+    func=problem.objective,
+    jac=problem.objective_jacobian,
+    hessian=problem.objective_hessian,
+    N=len(X0),
+    X0=X0,
+    xl=-2,
+    xu=2,
+    max_iters=max_iters,
+    verbose=True,
+    preconditioning=False,
+)
+X, Y, stop = opt.run()
+
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(13, 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,
+        )
+
+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_{\text{final}}$", r"$Y_0$", "Pareto front"])
+
+ax22 = ax1.twinx()
+ax1.plot(range(1, len(opt.history_indicator_value) + 1), opt.history_indicator_value, "b-")
+ax22.semilogy(range(1, len(opt.history_R_norm) + 1), opt.history_R_norm, "g--")
+ax1.set_ylabel("HV", color="b")
+ax22.set_ylabel(r"$||R(\mathbf{X})||$", color="g")
+ax1.set_title("Performance")
+ax1.set_xlabel("iteration")
+ax1.set_xticks(range(1, max_iters + 1))
+
+plt.savefig(f"dent-example-{N}.pdf", dpi=1000)

hvd/problems/misc.py

@@ -11,6 +11,30 @@
 from .base import ConstrainedMOOAnalytical, MOOAnalytical
 
 
+class DENT(MOOAnalytical):
+    def __init__(self, **kwargs):
+        self.n_obj = 2
+        self.n_var = 2
+        self.xl = -2 * np.ones(self.n_var)
+        self.xu = 2 * np.ones(self.n_var)
+        super().__init__(**kwargs)
+
+    def _objective(self, x: jnp.ndarray) -> jnp.ndarray:
+        term1 = jnp.sqrt(1 + (x[0] + x[1]) ** 2) + jnp.sqrt(1 + (x[0] - x[1]) ** 2)
+        term2 = 0.85 * jnp.exp(-((x[0] - x[1]) ** 2))
+        y1 = 0.5 * (term1 + x[0] - x[1]) + term2
+        y2 = 0.5 * (term1 - x[0] + x[1]) + term2
+        return jnp.array([y1, y2])
+
+    def get_pareto_set(self, n: int = 100) -> np.ndarray:
+        x = np.linspace(-2, 2, n)
+        return np.c_[x, -x]
+
+    def get_pareto_front(self, n: int = 100) -> np.ndarray:
+        X = self.get_pareto_set(n)
+        return np.array([self.objective(x) for x in X])
+
+
 class CONV3(ConstrainedMOOAnalytical):
     def __init__(self, **kwargs):
         self.n_obj = 3