From 87819b33d81a8d4b3de0788471253f9757ca5fa0 Mon Sep 17 00:00:00 2001 From: wangronin Date: Fri, 12 Jul 2024 08:11:06 +0200 Subject: [PATCH] update --- examples/HVN/2D_example.py | 2 +- examples/MMD/example_DTLZ1.py | 6 +++--- hvd/base.py | 7 +++++-- hvd/mmd_newton.py | 36 ++++++++++++++++++++++++----------- 4 files changed, 34 insertions(+), 17 deletions(-) diff --git a/examples/HVN/2D_example.py b/examples/HVN/2D_example.py index 0c9a05d..d3b1e0e 100644 --- a/examples/HVN/2D_example.py +++ b/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) diff --git a/examples/MMD/example_DTLZ1.py b/examples/MMD/example_DTLZ1.py index e8b6223..f79073d 100644 --- a/examples/MMD/example_DTLZ1.py +++ b/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) diff --git a/hvd/base.py b/hvd/base.py index 370626e..2addb43 100644 --- a/hvd/base.py +++ b/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,7 +64,8 @@ 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: @@ -71,7 +73,8 @@ def evaluate_cstr_hess(self, x: np.ndarray) -> np.ndarray: # 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 diff --git a/hvd/mmd_newton.py b/hvd/mmd_newton.py index f816193..e715245 100644 --- a/hvd/mmd_newton.py +++ b/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()