diff --git a/examples/MMD/example_DTLZ1.py b/examples/MMD/example_DTLZ1.py
index 94abf1f..385a155 100644
--- a/examples/MMD/example_DTLZ1.py
+++ b/examples/MMD/example_DTLZ1.py
@@ -30,8 +30,7 @@
 rcParams["ytick.major.width"] = 1
 
 np.random.seed(66)
-max_iters = 15
-
+max_iters = 17
 f = DTLZ1(boundry_constraints=True)
 problem = PymooProblemWithAD(f)
 pareto_front = problem.get_pareto_front()
@@ -63,8 +62,6 @@
 ref = {0: ref1_, 1: ref2_}
 eta = {0: -1 * np.array([1 / np.sqrt(3)] * 3), 1: -1 * np.array([1 / np.sqrt(3)] * 3)}
 Y_idx = [list(range(0, len(Y_component1))), list(range(len(Y_component1), len(Y0)))]
-Y_label = np.array([0] * len(X_component1) + [1] * len(X_component2))
-
 
 # X0 = X_component2
 # Y0 = Y_component2
@@ -140,7 +137,7 @@
     preconditioning=False,
 )
 if 1 < 2:
-    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)
@@ -163,7 +160,7 @@
     verbose=True,
     type="igd",
     eta=eta,
-    Y_label=Y_label,
+    Y_label=None,
     pareto_front=pareto_front,
 )
 X_DpN, Y_DpN, _ = opt_dpn.run()
@@ -179,9 +176,7 @@
 ax0.set_box_aspect((1, 1, 1))
 ax0.view_init(45, 45)
 ax0.plot(Y0[:, 0], Y0[:, 1], Y0[:, 2], "r+", ms=8, alpha=0.6)
-# ax0.plot(ref[:, 0], ref[:, 1], ref[:, 2], "g.", ms=6, alpha=0.6)
-ax0.plot(ref1_[:, 0], ref1_[:, 1], ref1_[:, 2], "g.", ms=6, alpha=0.6)
-ax0.plot(ref2_[:, 0], ref2_[:, 1], ref2_[:, 2], "g.", 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.", mec="none", ms=5, alpha=0.4)
 
 ax0.set_title("Initialization")
diff --git a/hvd/mmd_newton.py b/hvd/mmd_newton.py
index bc3eee1..565d62a 100644
--- a/hvd/mmd_newton.py
+++ b/hvd/mmd_newton.py
@@ -396,9 +396,8 @@ def bootstrap_reference_set(optimizer, problem, interval: int = 5) -> Tuple[np.n
             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]
-            # ref_archive.append(Y[idx].copy())
-            # ref = np.concatenate(ref_archive, axis=0)
+            ref_archive.append(Y[idx].copy())
+            ref = np.concatenate(ref_archive, axis=0)
             # km = KMedoids(
             #     n_clusters=optimizer.state.N, method="alternate", random_state=0, init="k-medoids++"
             # ).fit(ref)
@@ -414,7 +413,7 @@ def bootstrap_reference_set(optimizer, problem, interval: int = 5) -> Tuple[np.n
             optimizer.indicator.ref = optimizer.ref = ref
             optimizer.indicator.compute(Y=optimizer.state.Y)  # to compute the medoids
 
-            if 1 < 2:
+            if 11 < 2:
                 import matplotlib.pyplot as plt
 
                 pareto_front = problem.get_pareto_front()
@@ -432,12 +431,10 @@ def bootstrap_reference_set(optimizer, problem, interval: int = 5) -> Tuple[np.n
                     ms=6,
                     alpha=0.6,
                 )
-                rs = ref.reference_set
-                ax0.plot(rs[:, 0], rs[:, 1], rs[:, 2], "g.", 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.3)
                 ax0.plot(m[:, 0], m[:, 1], m[:, 2], "r^", ms=6, alpha=0.6)
                 plt.show()
-                breakpoint()
 
         optimizer.newton_iteration()
         optimizer.log()