diff --git a/examples/examples_bbo/plot_black_box_optimization_noisy.py b/examples/examples_bbo/plot_black_box_optimization_noisy.py
index 9d2cc24a..65a818c9 100644
--- a/examples/examples_bbo/plot_black_box_optimization_noisy.py
+++ b/examples/examples_bbo/plot_black_box_optimization_noisy.py
@@ -7,7 +7,7 @@
 In this tutorial, we show you how to manage **noisy** `black-box optimization (Wikipedia) <https://en.wikipedia.org/wiki/Derivative-free_optimization>`_ (a.k.a., derivative-free optimization) with DeepHyper.
 
 Black-box optimization is a field of optimization research where an objective function :math:`f(x) = y \in \mathbb{R}` is optimized only based on input-output observations :math:`\{ (x_1,y_1), \ldots, (x_n, y_n) \}`.
- 
+
 Let's start by installing DeepHyper!
 """
 
@@ -59,6 +59,7 @@ def f(job):
     obs = np.random.binomial(n=1, p=p)
     return obs
 
+
 # %%
 # Search Space of Input Variables
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -79,7 +80,7 @@ def f(job):
 # %%
 # Evaluator Interface
 # -------------------
-# 
+#
 # DeepHyper uses an API called :class:`deephyper.evaluator.Evaluator` to distribute the computation of black-box functions and adapt to different backends (e.g., threads, processes, MPI, Ray). An ``Evaluator`` object wraps the black-box function ``f`` that we want to optimize. Then a ``method`` parameter is used to select the backend and ``method_kwargs`` defines some available options of this backend.
 #
 #
@@ -95,19 +96,18 @@ def f(job):
 evaluator = Evaluator.create(
     f,
     method="thread",
-    method_kwargs={
-        "num_workers": 1,
-        "callbacks": [TqdmCallback()]
-    },
+    method_kwargs={"num_workers": 1, "callbacks": [TqdmCallback()]},
 )
 
-print(f"Evaluator has {evaluator.num_workers} available worker{'' if evaluator.num_workers == 1 else 's'}")
+print(
+    f"Evaluator has {evaluator.num_workers} available worker{'' if evaluator.num_workers == 1 else 's'}"
+)
 
 # %%
 # Search Algorithm
 # ----------------
-# 
-# The next step is to define the search algorithm that we want to use. Here, we choose :class:`deephyper.hpo.CBO` (Centralized Bayesian Optimization) which is a sampling based Bayesian optimization strategy. 
+#
+# The next step is to define the search algorithm that we want to use. Here, we choose :class:`deephyper.hpo.CBO` (Centralized Bayesian Optimization) which is a sampling based Bayesian optimization strategy.
 # This algorithm has the advantage of being asynchronous which is crutial to keep a good utilization of the resources when the number of available workers increases.
 # We also choose, how to optimize the acquisition function of the Bayesian optimization with ``"ga"`` (i.e., continuous Genetic Algorithm).
 #
@@ -130,13 +130,14 @@ def create_search():
     )
     return search
 
+
 max_evals = 300
 search = create_search()
 results = search.search(evaluator, max_evals)
 
 # %%
 # Finally, let us visualize the results. The ``search(...)`` returns a DataFrame also saved locally under ``results.csv`` (in case of crash we don't want to lose the possibly expensive evaluations already performed).
-# 
+#
 # The DataFrame contains the usual columns:
 #
 # 1. the optimized hyperparameters: such as :math:`x` with name ``p:x``.
@@ -159,7 +160,9 @@ def create_search():
 from deephyper.analysis.hpo import parameters_at_max
 
 
-parameters, objective = parameters_at_max(results, column="sol.objective", prefix="sol.p:", n_last=20)
+parameters, objective = parameters_at_max(
+    results, column="sol.objective", prefix="sol.p:", n_last=20
+)
 print("\nEstimated Optimum values")
 print("x:", parameters["x"])
 print("objective:", objective)
@@ -178,4 +181,4 @@ def create_search():
 plot_search_trajectory_single_objective_hpo(results, column="sol.p:x", mode="max", ax=ax)
 _ = ax.set_ylabel(r"Estimated solution $x$")
 _ = ax.set_ylim(-10, 10)
-_ = plt.title("Search Trajectory")
\ No newline at end of file
+_ = plt.title("Search Trajectory")
diff --git a/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler.py b/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler.py
index e9fd29ab..54f797e4 100644
--- a/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler.py
+++ b/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler.py
@@ -70,7 +70,7 @@
     parameters_at_max,
     filter_failed_objectives,
 )
-from deephyper.hpo import HpProblem, CBO
+from deephyper.hpo import HpProblem, CBO, RandomSearch, RegularizedEvolution
 
 # %%
 # Custom Sampler
@@ -85,10 +85,15 @@
 # - all generated samples satisfy :math:`x_i < x_{i+1}` by construction;
 # - the sampler focuses on the feasible region, avoiding wasted evaluations.
 
-n = 10
-m = 32
+m = 80
+n = 40
 
-print("optimum:", sum([m - i - 1 for i in range(n)]))
+MAX_SUM = False
+if MAX_SUM:
+    optimum_value = sum([m - i - 1 for i in range(n)])
+else:
+    optimum_value = -sum([i for i in range(n)])
+print("optimum:", optimum_value)
 
 pb = HpProblem()
 for i in range(n):
@@ -105,8 +110,8 @@ def sampling_fn(size: int) -> list[dict]:
 
     def sample_one():
         vals = np.sort(rng.choice(indexes, size=n, replace=False)).tolist()
-        return {k: v for k, v in zip(pb.hyperparameter_names, vals)}
-    
+        return {f"x{i}": v for i, v in zip(range(n), vals)}
+
     return [sample_one() for _ in range(size)]
 
 
@@ -124,22 +129,32 @@ def sample_one():
 # in the objective function so that CBO learns to avoid them.
 
 
-def constraint_fn(df: pd.DataFrame) -> pd.Series:
-    accept = pd.Series(np.ones((len(df)), dtype=bool))
-    for i in range(n - 1):
-        accept = accept & (df[f"x{i}"] < df[f"x{i + 1}"])
-    return accept
+def constraint_fn1(df: pd.DataFrame) -> pd.Series:
+    """Returns booleans for equality constraint."""
+    x = df[[f"x{i}" for i in range(n)]].to_numpy()
+    violations = (x[:, :-1] < x[:, 1:]).all(axis=1)
+    return pd.Series(violations, index=df.index)
+
+
+def constraint_fn2(df: pd.DataFrame) -> pd.Series:
+    x = df[[f"x{i}" for i in range(n)]].to_numpy()
+    diff = x[:, :-1] - x[:, 1:] + 1  # if diff >= 0 → violation
+    violations = diff.sum(axis=1)
+    return pd.Series(violations, index=df.index)
 
 
-pb.set_constraint_fn(constraint_fn)
+pb.set_constraint_fn(constraint_fn1)
 
 
 def f(job):
     """Objective function: maximize sum(x_i)."""
     df = pd.DataFrame([job.parameters])
-    accept = constraint_fn(df)
+    accept = constraint_fn1(df)
     if all(accept):
-        return sum(job.parameters.values())
+        if MAX_SUM:
+            return sum([job.parameters[f"x{i}"] for i in range(n)])
+        else:
+            return -sum([job.parameters[f"x{i}"] for i in range(n)])
     else:
         return "F_constraint"
 
@@ -159,29 +174,58 @@ def f(job):
 #
 # The search runs for ``max_evals=300`` iterations.
 
-search = CBO(
-    pb,
-    surrogate_model="ET",
-    surrogate_model_kwargs={"max_features": "sqrt"},
-    acq_optimizer="mixedga",
-    acq_optimizer_kwargs={
-        "n_points": 1_000,
-        "acq_optimizer_freq": 2,
-        "filter_failures": "mean",
-    },
-    acq_func_kwargs={
-        # Exploration/Exploitation mechanism
-        "kappa": 10.0,
-        "scheduler": {
-            "type": "periodic-exp-decay",
-            "period": 20,
-            "kappa_final": 0.1,
-        },
-    },
-    verbose=1,
-)
-results = search.search(f, max_evals=300)
 
+def make_search():
+    search = CBO(
+        pb,
+        surrogate_model="ET",
+        surrogate_model_kwargs={
+            "max_features": "sqrt",
+            # "min_samples_split": 2,
+            # "min_samples_leaf": 1,
+            "bootstrap": True,
+        },
+        acq_optimizer="mixedga",
+        acq_optimizer_kwargs={
+            "n_points": 10_000,
+            "acq_optimizer_freq": 2,
+            "filter_failures": "mean",
+            # "ga_pop_size": 100,
+        },
+        acq_func_kwargs={
+            # Exploration/Exploitation mechanism
+            "kappa": 1.96,
+            # "scheduler": {
+            #     "type": "periodic-exp-decay",
+            #     "period": 20,
+            #     "kappa_final": 0.01,
+            # },
+        },
+        n_initial_points=10,
+        initial_points=[pb.default_configuration],
+        objective_scaler="identity",
+        verbose=1,
+    )
+    # search = RandomSearch(pb, verbose=1)
+    # search = RegularizedEvolution(pb, verbose=1)
+    return search
+
+
+# n_repetitions = 10
+# results_total = []
+# for i in range(n_repetitions):
+#     search = make_search()
+#     results = search.search(f, max_evals=50)
+#     if i > 0:
+#         search.fit_surrogate(pd.concat(results_total + [results]))
+#     results = search.search(f, max_evals=50)
+#     results_total.append(results)
+#     print(len(results_total))
+
+# results = pd.concat(results_total)
+
+search = make_search()
+results = search.search(f, max_evals=500)
 # %%
 results
 
@@ -216,8 +260,15 @@ def f(job):
 HEIGHT_PLOTS = WIDTH_PLOTS / 1.618
 
 fig, ax = plt.subplots(figsize=(WIDTH_PLOTS, HEIGHT_PLOTS))
-plot_search_trajectory_single_objective_hpo(results, mode="max", ax=ax)
+plot_search_trajectory_single_objective_hpo(results, mode="max" if MAX_SUM else "min", ax=ax)
+y_max = sum([m - i - 1 for i in range(n)])
+y_min= sum([i for i in range(n)])
+ax.set_ylim(y_min-10, y_max + 10)
+yhline = optimum_value if MAX_SUM else -optimum_value
+ax.axhline(yhline, color="black", linestyle="--", linewidth=2, label="Optimum")
+ax.legend()
 _ = plt.title("Search Trajectory")
+plt.show()
 
 # %%
 # Visualizing the Feasible Region and Evaluations
@@ -231,7 +282,7 @@ def f(job):
 
 results, _ = filter_failed_objectives(results)
 
-p_columns = [col for col in results.columns if col.startswith("p:")]
+p_columns = [f"p:x{i}" for i in range(n)]
 
 # Create a normalizer over the objective range
 obj_vals = results["objective"]
diff --git a/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler_reformulation.py b/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler_reformulation.py
new file mode 100644
index 00000000..28c00e0c
--- /dev/null
+++ b/examples/examples_bbo/plot_constrained_black_box_optimization_chained_sampler_reformulation.py
@@ -0,0 +1,337 @@
+r"""
+Constrained Black-Box Optimization with Custom Sampler
+======================================================
+
+**Author(s)**: Romain Egele.
+
+This tutorial demonstrates how to solve *constrained*
+`black-box optimization <https://en.wikipedia.org/wiki/Derivative-free_optimization>`_
+using **DeepHyper**, focusing on how to encode structural constraints directly
+in the **sampling strategy** of the search algorithm.
+
+Black-box optimization aims to optimize an unknown function
+:math:`f(x) = y \in \mathbb{R}` using only input–output evaluations
+:math:`\{(x_1, y_1), \ldots, (x_n, y_n)\}`.
+No analytical gradients or structural properties of :math:`f` are required.
+
+In *constrained* settings, the search is further restricted to parameters that
+satisfy one or more feasibility rules. Constraints can significantly reshape
+the search space and modify the behavior of the optimizer.
+
+Problem Setting
+---------------
+In this example, we consider a *discrete*, *ordered* search space of
+dimension :math:`N`.
+Each variable :math:`x_i` must satisfy the **monotonicity constraint**
+
+.. math::
+
+    x_0 < x_1 < \cdots < x_{N-1}.
+
+Each :math:`x_i` is bounded between :math:`i` and :math:`m - N + i`.
+This could tipically represent the layer indexes to drop in Depth pruning of Large language models.
+The objective is to **maximize** the sum
+
+.. math::
+
+   f(x) = \sum_{i=0}^{N-1} x_i.
+
+Since the optimal strategy is to push every variable as high as possible while
+respecting monotonicity, the theoretical optimum is:
+
+.. math::
+
+   \sum_{i=0}^{N-1} (m - i).
+
+DeepHyper offers several ways to incorporate constraints:
+
+#. **Custom sampler** *(this tutorial)*: constraints are enforced
+   directly when generating new candidate points.
+
+#. **Rejection sampling**:
+   see :ref:`Constrained Black-Box Optimization with Rejection Sampling <sphx_glr_examples_examples_bbo_plot_constrained_black_box_optimization.py>`.
+
+#. **Learning to avoid failures** (CBO auto-handles failed evaluations):
+   see this tutorial and also
+   :ref:`Learn to Avoid Failures with Bayesian Optimization <sphx_glr_examples_examples_bbo_plot_notify_failures_hpo.py>`.
+
+#. **Multi-objective approach** where the constraint becomes an additional
+   objective (tutorial forthcoming).
+"""
+
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import matplotlib.colors as colors
+import numpy as np
+import pandas as pd
+from numpy.random import Generator
+from deephyper.analysis.hpo import (
+    plot_search_trajectory_single_objective_hpo,
+    parameters_at_max,
+    filter_failed_objectives,
+)
+from deephyper.hpo import HpProblem, CBO, RandomSearch, RegularizedEvolution
+
+# %%
+# Custom Sampler
+# --------------
+#
+# Because every :math:`x_i` must be strictly larger than :math:`x_{i-1}`, the
+# usual independent sampling over each variable would frequently violate the
+# constraint.
+#
+# Instead, we implement a custom sampler to ensure that:
+#
+# - all generated samples satisfy :math:`x_i < x_{i+1}` by construction;
+# - the sampler focuses on the feasible region, avoiding wasted evaluations.
+
+m = 80
+n = 40
+
+MAX_SUM = False
+if MAX_SUM:
+    optimum_value = sum([m - i - 1 for i in range(n)])
+else:
+    optimum_value = -sum([i for i in range(n)])
+print("optimum:", optimum_value)
+
+pb = HpProblem()
+pb.add((0, m-n), "g0")
+for i in range(1, n):
+    pb.add((1, m - n + 1), f"g{i}")
+print(pb)
+
+
+def sampling_fn(size: int) -> list[dict]:
+    if n > m:
+        raise ValueError(f"Cannot sample {n} items from {m} elements.")
+
+    rng = np.random.default_rng(None)
+    indexes = np.arange(m)
+
+    def sample_one():
+        # Sample in x_i space then map to g_i space
+        vals = np.sort(rng.choice(indexes, size=n, replace=False))
+        vals[1:] = vals[1:] - vals[:-1]
+        return {f"g{i}": v for i, v in zip(range(n), vals.tolist())}
+
+    return [sample_one() for _ in range(size)]
+
+
+pb.set_sampling_fn(sampling_fn)
+
+# %%
+# Constraint Function
+# -------------------
+#
+# Although the sampler already generates feasible points, we explicitly define a
+# ``constraint_fn``. This allows DeepHyper to properly handle *failed* trials
+# (e.g., from manually constructed parameter sets or mutation-based acquisition
+# optimizers).
+# Not only that, this will help report non-feasible points using ``"F_constraint"``
+# in the objective function so that CBO learns to avoid them.
+
+
+def constraint_fn1(df: pd.DataFrame) -> pd.Series:
+    """Returns booleans for equality constraint."""
+    g = df[[f"g{i}" for i in range(n)]].to_numpy()
+    x = np.cumsum(g, axis=1)
+    violations = (x < m).all(axis=1) & (x[:, :-1] < x[:, 1:]).all(axis=1)
+    return pd.Series(violations, index=df.index)
+
+
+def constraint_fn2(df: pd.DataFrame) -> pd.Series:
+    x = df[[f"x{i}" for i in range(n)]].to_numpy()
+    diff = x[:, :-1] - x[:, 1:] + 1  # if diff >= 0 → violation
+    violations = diff.sum(axis=1)
+    return pd.Series(violations, index=df.index)
+
+
+pb.set_constraint_fn(constraint_fn1)
+
+def repair_fn(df: pd.DataFrame | dict) -> pd.DataFrame | dict:
+    g_columns = [f"g{i}" for i in range(n)]
+    g = df[g_columns].to_numpy()
+    x = np.cumsum(g, axis=1)
+    x = np.clip(x, np.arange(n), m-n+np.arange(n))
+    x[:, 1:] = x[:, 1:] - x[:, :-1]
+    g = x
+    df[g_columns] = g
+    return df
+
+pb.set_repair_fn(repair_fn)
+
+
+def f(job):
+    """Objective function: maximize sum(x_i)."""
+    df = pd.DataFrame([job.parameters])
+    accept = constraint_fn1(df)
+    if all(accept):
+        if MAX_SUM:
+            return sum(np.cumsum([job.parameters[f"g{i}"] for i in range(n)]))
+        else:
+            return -sum(np.cumsum([job.parameters[f"g{i}"] for i in range(n)]))
+    else:
+        return "F_constraint"
+
+
+# %%
+# Bayesian Optimization with Mixed-GA Acquisition Optimization
+# ------------------------------------------------------------
+#
+# We run a **Centralized Bayesian Optimization (CBO)** search using:
+#
+# - Ensemble of Trees surrogate model (``"ET"``).
+# - A **mixed genetic algorithm** (``"mixedga"``) to optimize the acquisition
+#   function.
+# - A **periodically decaying scheduler** on the exploration parameter ``kappa``.
+#
+# This setup is well suited for discrete, irregularly constrained spaces.
+#
+# The search runs for ``max_evals=300`` iterations.
+
+
+def make_search():
+    search = CBO(
+        pb,
+        surrogate_model="ET",
+        surrogate_model_kwargs={
+            # "max_features": "sqrt",
+            # "min_samples_split": 2,
+            # "min_samples_leaf": 1,
+            # "bootstrap": True,
+        },
+        acq_optimizer="mixedga",
+        acq_optimizer_kwargs={
+            "n_points": 10_000,
+            "acq_optimizer_freq": 2,
+            "filter_failures": "mean",
+            # "ga_pop_size": 100,
+        },
+        acq_func_kwargs={
+            # Exploration/Exploitation mechanism
+            "kappa": 10.96,
+            "scheduler": {
+                "type": "periodic-exp-decay",
+                "period": 20,
+                "kappa_final": 0.01,
+            },
+        },
+        n_initial_points=10,
+        objective_scaler="identity",
+        verbose=1,
+    )
+    # search = RandomSearch(pb, verbose=1)
+    # search = RegularizedEvolution(pb, verbose=1)
+    return search
+
+
+# n_repetitions = 10
+# results_total = []
+# for i in range(n_repetitions):
+#     search = make_search()
+#     results = search.search(f, max_evals=50)
+#     if i > 0:
+#         search.fit_surrogate(pd.concat(results_total + [results]))
+#     results = search.search(f, max_evals=50)
+#     results_total.append(results)
+#     print(len(results_total))
+
+# results = pd.concat(results_total)
+
+search = make_search()
+results = search.search(f, max_evals=500)
+# %%
+results
+
+# %%
+# Extracting the Best Parameters
+# ------------------------------
+# To recover the parameters corresponding to the best observed objective value,
+# we can use :func:`deephyper.analysis.hpo.parameters_at_max`.
+
+g_columns = [f"p:g{i}" for i in range(n)]
+x_columns = [f"p:x{i}" for i in range(n)]
+
+g = results[g_columns].to_numpy()
+results[x_columns] = np.cumsum(g, axis=1)
+results.drop(columns=g_columns, inplace=True)
+
+parameters, objective = parameters_at_max(results)
+print("\nOptimum values")
+for i in range(n):
+    print(f"x{i}: {parameters[f'x{i}']:.3f}")
+print("objective:", objective)
+
+# %%
+# Visualization
+# ---------------
+# We conclude with:
+#
+# - a **search trajectomakery plot** showing the best objective value over time,
+#   where the periodic exploration schedule is clearly visible;
+#
+# - a **feasible-space evaluation plot** showing all sampled curves
+#   :math:`i \mapsto x_i` (each curve is one evaluation), colored by objective
+#   value.
+#
+# These visualizations confirm that the optimizer progressively learns the
+# structure of the monotonic constraint and approaches the theoretical optimum.
+
+WIDTH_PLOTS = 8
+HEIGHT_PLOTS = WIDTH_PLOTS / 1.618
+
+fig, ax = plt.subplots(figsize=(WIDTH_PLOTS, HEIGHT_PLOTS))
+plot_search_trajectory_single_objective_hpo(results, mode="max" if MAX_SUM else "min", ax=ax)
+y_max = sum([m - i - 1 for i in range(n)])
+y_min= sum([i for i in range(n)])
+ax.set_ylim(y_min-10, y_max + 10)
+yhline = optimum_value if MAX_SUM else -optimum_value
+ax.axhline(yhline, color="black", linestyle="--", linewidth=2, label="Optimum")
+ax.legend()
+_ = plt.title("Search Trajectory")
+plt.show()
+
+# %%
+# Visualizing the Feasible Region and Evaluations
+# -----------------------------------------------
+# We now plot all evaluated points in the (x, y) plane, color-coded by
+# objective value, along with the constraint boundary ``x + y = 10``.
+
+
+# sphinx_gallery_thumbnail_number = 2
+
+
+results, _ = filter_failed_objectives(results)
+
+p_columns = [f"p:x{i}" for i in range(n)]
+
+# Create a normalizer over the objective range
+obj_vals = results["objective"]
+norm = colors.Normalize(vmin=obj_vals.min(), vmax=obj_vals.max())
+
+# Choose a colormap (viridis is a good default)
+cmap = plt.get_cmap("viridis")
+
+fig, ax = plt.subplots(figsize=(WIDTH_PLOTS, HEIGHT_PLOTS))
+
+for i, row in results.iterrows():
+    x_values = row[p_columns].values
+    y_values = np.arange(n)
+    obj_value = row["objective"]
+
+    color = cmap(norm(obj_value))  # map objective → color
+    ax.plot(x_values, y_values, color=color, alpha=0.9)
+
+# Optionally add a colorbar
+sm = cm.ScalarMappable(norm=norm, cmap=cmap)
+cbar = fig.colorbar(sm, ax=ax)
+cbar.set_label("Objective value")
+ax.grid()
+ax.set_ylim(0, n - 1)
+ax.set_xlim(0, m)
+ax.set_ylabel(r"$i$")
+ax.set_xlabel(r"$x_i$")
+ax.set_yticks(list(range(n)), [str(i) for i in range(n)])
+ax.set_xticks(list(range(0, m, 2)), [str(i) for i in range(0, m, 2)])
+plt.show()
diff --git a/src/deephyper/hpo/_problem.py b/src/deephyper/hpo/_problem.py
index 27d9e14e..e33644b1 100644
--- a/src/deephyper/hpo/_problem.py
+++ b/src/deephyper/hpo/_problem.py
@@ -1,10 +1,14 @@
 import copy
+from typing import Callable
 import warnings
 
+from numbers import Number
+
 import ConfigSpace as cs
 import ConfigSpace.hyperparameters as csh
 import numpy as np
 import pandas as pd
+
 from sklearn.utils import check_random_state
 
 import deephyper.skopt
@@ -224,6 +228,7 @@ def __init__(self, config_space=None, seed: int | None = None):
 
         self.constraint_fn = None
         self.sampling_fn = None
+        self.repair_fn = None
 
     def __str__(self):
         return repr(self)
@@ -359,7 +364,7 @@ def add(self, value, name=None, default_value=None) -> None:
     def sample(
         self,
         size: int = 1,
-        strict: bool = False,
+        strict: bool = True,
         max_trials: int = 5,
         n_jobs: int = 1,
     ) -> list[dict]:
@@ -375,6 +380,8 @@ def sample(
         Returns:
             list[dict]: the list of sampled configurations.
         """
+        if size < 1:
+            raise ValueError(f"{size=} should be > 0")
 
         def _sample_dimension(dim, i, n_samples, random_state, out):
             """Wrapper to sample dimension for joblib parallelization."""
@@ -387,6 +394,8 @@ def sample_fn(size: int) -> list[dict]:
                 )
                 if sample_with_config_space:
                     samples = self._space.sample_configuration(size=size)
+                    if size == 1:
+                        samples = [samples]
                     samples = [dict(s) for s in samples]
                 else:
                     # Regular sampling without transfer learning from flat search space
@@ -432,9 +441,13 @@ def sample_fn(size: int) -> list[dict]:
             # Convert batch into DataFrame only once
             df = pd.DataFrame(batch)
 
-            # Apply constraint ---
+            # Apply constraint
             accept_mask = self.constraint_fn(df)
 
+            # If constraint value is defined numerically to leverage "constraint domination"
+            if isinstance(accept_mask[0], Number):
+                accept_mask = ~(accept_mask > 0)
+
             df = df[accept_mask]
             accepted.extend(df.to_dict(orient="records"))
 
@@ -451,13 +464,25 @@ def sample_fn(size: int) -> list[dict]:
             else:
                 batch_size = int((size - len(accepted)) / ratio_accept + 0.5)
 
+        # Apply repair if available
+        if self.repair_fn:
+            if len(accepted) < size:
+                # Sample a batch
+                batch = sample_fn(size - len(accepted))
+
+                # Convert batch into DataFrame only once
+                df = pd.DataFrame(batch)
+
+                df = self.repair_fn(df)
+                accepted.extend(df.to_dict(orient="records"))
+
         # If constraints are too strict, return what we have (or raise)
         # You can choose to raise if you need strictly size samples
         if strict:
             accepted = accepted[:size]
 
             if len(accepted) < size:
-                return RuntimeError(f"The number of samples is less than {size=}!")
+                raise RuntimeError(f"The number of samples is less than {size=}!")
 
         return accepted
 
@@ -556,3 +581,41 @@ def constraint_fn(df: pd.DataFrame) -> pd.Series:
     def set_sampling_fn(self, fn: callable):
         """Set the sampling function."""
         self.sampling_fn = fn
+
+    def is_feasible(self, x: pd.DataFrame | dict) -> pd.Series | bool:
+        """Check if a configuration or set of configurations are feasible w.r.t. constraints."""
+        if isinstance(x, dict):
+            df = pd.DataFrame([x])
+        else:
+            df = x
+
+        # If a constraint function was defined it has priority
+        if self.constraint_fn is not None:
+            accept = self.constraint_fn(df)
+
+        # If there are forbidden or conditions defined through ConfigSpace
+        elif len(self._space.forbidden_clauses) > 0 or len(self._space.conditions) > 0:
+            accept = []
+            for conf in df.to_dict(orient="records"):
+                cs_conf = cs.Configuration(self._space, conf, allow_inactive_with_values=True)
+                try:
+                    cs_conf.check_valid_configuration()
+                except ValueError:
+                    accept.append(False)
+                else:
+                    accept.append(True)
+            accept = pd.Series(accept)
+
+        # Finally case if there is no constraints
+        else:
+            accept = pd.Series(np.ones(len(df)))
+
+        if isinstance(x, dict):
+            return accept.values.tolist()[0]
+        else:
+            return accept
+
+    def set_repair_fn(self, fn: Callable):
+        """Set the repair function."""
+        # x: pd.DataFrame | dict -> pd.DataFrame | dict
+        self.repair_fn = fn
diff --git a/src/deephyper/hpo/_random.py b/src/deephyper/hpo/_random.py
index 95a58cf7..b2703435 100644
--- a/src/deephyper/hpo/_random.py
+++ b/src/deephyper/hpo/_random.py
@@ -2,9 +2,11 @@
 
 import numpy as np
 
+from deephyper.hpo._problem import HpProblem
 from deephyper.hpo._search import Search
 from deephyper.hpo._solution import SolutionSelection
 from deephyper.hpo.utils import get_inactive_value_of_hyperparameter
+from deephyper.stopper._stopper import Stopper
 
 __all__ = ["RandomSearch"]
 
@@ -52,15 +54,13 @@ class RandomSearch(Search):
 
     def __init__(
         self,
-        problem,
-        random_state=None,
-        log_dir=".",
-        verbose=0,
-        stopper=None,
+        problem: HpProblem,
+        random_state: int | np.random.RandomState | None = None,
+        log_dir: str = ".",
+        verbose: int = 0,
+        stopper: Stopper | None = None,
         checkpoint_history_to_csv: bool = True,
-        solution_selection: Optional[
-            Literal["argmax_obs", "argmax_est"] | SolutionSelection
-        ] = None,
+        solution_selection: Literal["argmax_obs", "argmax_est"] | SolutionSelection | None = None,
     ):
         super().__init__(
             problem,
@@ -86,7 +86,7 @@ def _ask(self, n: int = 1) -> list[dict[str, Optional[str | int | float]]]:
 
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", category=DeprecationWarning)
-            new_samples = self._problem.space.sample_configuration(size=n)
+            new_samples = self._problem.sample(size=n)
 
         if not (isinstance(new_samples, list)):
             new_samples = [new_samples]
diff --git a/src/deephyper/hpo/_regevo.py b/src/deephyper/hpo/_regevo.py
index 4c0766a8..52c9d5b1 100644
--- a/src/deephyper/hpo/_regevo.py
+++ b/src/deephyper/hpo/_regevo.py
@@ -1,13 +1,15 @@
 from collections import deque
-from typing import Dict, List, Literal, Optional
+from typing import Any, Literal, Optional
 
 import numpy as np
-
+import pandas as pd
 from ConfigSpace.util import deactivate_inactive_hyperparameters
 
+from deephyper.hpo._problem import HpProblem
 from deephyper.hpo._search import Search
 from deephyper.hpo._solution import SolutionSelection
 from deephyper.hpo.utils import get_inactive_value_of_hyperparameter
+from deephyper.stopper._stopper import Stopper
 
 __all__ = ["RegularizedEvolution"]
 
@@ -59,21 +61,25 @@ class RegularizedEvolution(Search):
 
         sample_size (int, optional):
             The number of samples to draw from the population. Defaults to ``10``.
+
+        max_trials_rejection_sampling (int):
+            The maximum number of trials for rejection sampling with resolving contraints. Defaults
+            to ``-1``.
     """
 
     def __init__(
         self,
-        problem,
-        random_state=None,
-        log_dir=".",
-        verbose=0,
-        stopper=None,
+        problem: HpProblem,
+        random_state: int | np.random.RandomState | None = None,
+        log_dir: str = ".",
+        verbose: int = 0,
+        stopper: Stopper | None = None,
         checkpoint_history_to_csv: bool = True,
-        solution_selection: Optional[
-            Literal["argmax_obs", "argmax_est"] | SolutionSelection
-        ] = None,
+        solution_selection: Literal["argmax_obs", "argmax_est"] | SolutionSelection | None = None,
         population_size: int = 100,
         sample_size: int = 10,
+        max_trials_rejection_sampling: int = -1,
+        init_population: list[tuple[dict, Any]] | None = None,
     ):
         super().__init__(
             problem,
@@ -88,9 +94,14 @@ def __init__(
         assert population_size > sample_size, "population_size must be greater than sample_size"
         self.population_size = population_size
         self.sample_size = sample_size
-        self._population = deque(maxlen=self.population_size)
+        if init_population is None:
+            init_population = []
+        self._population: deque[tuple[dict, Any]] = deque(
+            init_population, maxlen=self.population_size
+        )
+        self._max_trials_rejection_sampling = max_trials_rejection_sampling
 
-    def _ask(self, n: int = 1) -> List[Dict]:
+    def _ask(self, n: int = 1) -> list[dict[str, Any]]:
         """Ask the search for new configurations to evaluate.
 
         Args:
@@ -99,8 +110,6 @@ def _ask(self, n: int = 1) -> List[Dict]:
         Returns:
             List[Dict]: a list of hyperparameter configurations to evaluate.
         """
-        space = self._problem.space
-
         # Random sampling
         if len(self._population) < self.population_size:
             import warnings
@@ -108,22 +117,15 @@ def _ask(self, n: int = 1) -> List[Dict]:
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore", category=DeprecationWarning)
 
-                new_samples = space.sample_configuration(size=n)
+                new_samples = self._problem.sample(size=n)
 
             if not (isinstance(new_samples, list)):
                 new_samples = [new_samples]
 
             for i, sample in enumerate(new_samples):
                 sample = dict(sample)
-                for hp_name in self._problem.hyperparameter_names:
-                    # If the parameter is inactive due to some conditions then we attribute the
-                    # lower bound value to break symmetries and enforce the same representation.
-                    if hp_name not in sample:
-                        sample[hp_name] = get_inactive_value_of_hyperparameter(space[hp_name])
 
-                    # Make sure to have JSON serializable values
-                    if type(sample[hp_name]).__module__ == np.__name__:
-                        sample[hp_name] = sample[hp_name].tolist()
+                sample = self._set_inactive(sample)
 
                 new_samples[i] = sample
 
@@ -131,40 +133,41 @@ def _ask(self, n: int = 1) -> List[Dict]:
         else:
             new_samples = []
             for i in range(n):
+                # Get a sample of parents from the population
                 samples_idxs = self._random_state.choice(
                     self.population_size, size=self.sample_size, replace=False
                 )
 
                 samples = [self._population[i] for i in samples_idxs]
 
+                # Select the parent
                 parent_sample = max(samples, key=lambda x: x[1])[0]
 
-                child_sample = parent_sample.copy()
-                active_hyperparameter_names = list(
-                    space.get_active_hyperparameters(
-                        deactivate_inactive_hyperparameters(child_sample, space)
-                    )
-                )
-                hp_name = self._random_state.choice(active_hyperparameter_names)
-                hp = space[hp_name]
-                hp_value = hp.rvs(size=None, random_state=space.random)
-
-                child_sample[hp_name] = hp_value
-                child_sample = dict(deactivate_inactive_hyperparameters(child_sample, space))
-
-                for hp_name in self._problem.hyperparameter_names:
-                    # If the parameter is inactive due to some conditions then we attribute the
-                    # lower bound value to break symmetries and enforce the same representation.
-                    if hp_name not in child_sample:
-                        child_sample[hp_name] = get_inactive_value_of_hyperparameter(
-                            self._problem.space[hp_name]
-                        )
+                # Produce the child
+                n_trials = 0
+                child_sample = self._mutate(parent_sample)
 
-                    # Make sure to have JSON serializable values
-                    if type(child_sample[hp_name]).__module__ == np.__name__:
-                        child_sample[hp_name] = child_sample[hp_name].tolist()
+                def is_not_max_trials():
+                    return (
+                        self._max_trials_rejection_sampling < 0
+                        or n_trials < self._max_trials_rejection_sampling
+                    )
 
-                new_samples.append(child_sample)
+                while is_not_max_trials() and not self._problem.is_feasible(child_sample):
+                    child_sample = self._mutate(parent_sample)
+                    n_trials += 1
+
+                # If we can't produce a feasible child for _max_trials_rejection_sampling
+                # Then we resample a new fresh child
+                if is_not_max_trials():
+                    new_samples.append(child_sample)
+                else:
+                    try:
+                        child_sample = self._problem.sample(size=1, strict=True)[0]
+                    except RuntimeError:
+                        raise RuntimeError(
+                            "Could not resolve constraints through rejection sampling!"
+                        )
 
         return new_samples
 
@@ -182,3 +185,47 @@ def _tell(
             if isinstance(obj, str):
                 continue
             self._population.append((config, obj))
+
+    def _mutate(self, parent_sample: dict) -> dict:
+        space = self._problem.space
+        child_sample = parent_sample.copy()
+        active_hyperparameter_names = list(
+            space.get_active_hyperparameters(
+                deactivate_inactive_hyperparameters(child_sample, space)
+            )
+        )
+        hp_name = self._random_state.choice(active_hyperparameter_names)
+        hp = space[hp_name]
+        hp_value = hp.rvs(size=None, random_state=space.random)
+
+        child_sample[hp_name] = hp_value
+        child_sample = dict(deactivate_inactive_hyperparameters(child_sample, space))
+
+        child_sample = self._set_inactive(child_sample)
+
+        child_sample = self._repair(child_sample)
+
+        return child_sample
+
+    def _set_inactive(self, sample: dict):
+        sample = sample.copy()
+        space = self._problem.space
+        for hp_name in self._problem.hyperparameter_names:
+            # If the parameter is inactive due to some conditions then we attribute the
+            # lower bound value to break symmetries and enforce the same repsresentation.
+            if hp_name not in sample:
+                sample[hp_name] = get_inactive_value_of_hyperparameter(space[hp_name])
+
+            # Make sure to have JSON serializable values
+            if type(sample[hp_name]).__module__ == np.__name__:
+                sample[hp_name] = sample[hp_name].tolist()
+        return sample
+
+    def _repair(self, sample: dict):
+        sample = sample.copy()
+        if self._problem.repair_fn is None:
+            return sample
+
+        df = pd.DataFrame([sample])
+        df = self._problem.repair_fn(df)
+        return df.to_dict(orient="records")[0]
diff --git a/src/deephyper/hpo/_search.py b/src/deephyper/hpo/_search.py
index d0d4b766..9665b108 100644
--- a/src/deephyper/hpo/_search.py
+++ b/src/deephyper/hpo/_search.py
@@ -6,8 +6,8 @@
 import os
 import pathlib
 import time
-from typing import Any, Dict, List, Literal, Optional
 from inspect import iscoroutinefunction
+from typing import Any, Dict, List, Literal, Optional
 
 import numpy as np
 import pandas as pd
@@ -16,7 +16,7 @@
     get_mask_of_rows_without_failures,
     read_results_from_csv,
 )
-from deephyper.evaluator import Evaluator, HPOJob, MaximumJobsSpawnReached, JobStatus
+from deephyper.evaluator import Evaluator, HPOJob, JobStatus, MaximumJobsSpawnReached
 from deephyper.evaluator.callback import TqdmCallback
 from deephyper.hpo._problem import HpProblem
 from deephyper.hpo._solution import (
@@ -25,8 +25,8 @@
     Solution,
     SolutionSelection,
 )
-from deephyper.stopper import Stopper
 from deephyper.skopt.moo import non_dominated_set
+from deephyper.stopper import Stopper
 
 __all__ = ["Search", "SearchHistory"]
 
@@ -239,6 +239,37 @@ def compute_pareto_efficiency(self):
             for job, pf in zip(self.jobs, self.pareto_efficient):
                 job.pareto_efficient = pf
 
+    def check_objective_convergence(self, objective_tol: float | None, window: int = 20):
+        """Check if an optimization has converged based on objective tolerance.
+
+        It will evaluate the average positive improvement over a window of most recent objectives.
+
+        Args:
+            objectives (list[float]):
+                Sequence of objective values.
+
+            objective_tol (float, optional):
+                Minimum significant improvement. If the absolute improvement
+                stays below this threshold over the last `window` steps, we
+                declare convergence.
+
+            window (int):
+                Number of recent improvements to examine.
+
+        Return:
+            bool: True if converged, False otherwise.
+        """
+        if objective_tol is None or len(self.jobs) <= window:
+            return False
+
+        recent = [j.objective for j in self.jobs[-window:] if not isinstance(j.objective, str)]
+        if len(recent) < window:
+            return False
+
+        improvements = [max(recent[i] - recent[i - 1], 0) for i in range(1, window)]
+        has_converged = sum(improvements) / (window - 1) < objective_tol
+        return has_converged
+
 
 class Search(abc.ABC):
     """Base class search/optimization algorithms.
@@ -274,15 +305,13 @@ class Search(abc.ABC):
 
     def __init__(
         self,
-        problem,
-        random_state=None,
+        problem: HpProblem,
+        random_state: int | np.random.RandomState | None = None,
         log_dir: str = ".",
         verbose: int = 0,
-        stopper: Optional[Stopper] = None,
+        stopper: Stopper | None = None,
         checkpoint_history_to_csv: bool = True,
-        solution_selection: Optional[
-            Literal["argmax_obs", "argmax_est"] | SolutionSelection
-        ] = None,
+        solution_selection: Literal["argmax_obs", "argmax_est"] | SolutionSelection | None = None,
         checkpoint_restart: bool = False,
     ):
         # get the __init__ parameters
@@ -441,6 +470,7 @@ def search(
         max_evals: int = -1,
         timeout: Optional[int | float] = None,
         max_evals_strict: bool = False,
+        objective_tol: float | None = None,
     ) -> pd.DataFrame:
         """Execute the search algorithm.
 
@@ -457,6 +487,8 @@ def search(
             max_evals_strict (bool, optional): If ``True`` the search will not spawn more than
                 ``max_evals`` jobs. Defaults to ``False``.
 
+            objective_tol (float, optional): ...
+
         Returns:
             pd.DataFrame: A pandas DataFrame containing the evaluations performed or ``None`` if the
                 search could not evaluate any configuration.
@@ -473,6 +505,8 @@ def search(
         """
         logger.info(f"Starting search with {type(self).__name__}")
 
+        self._objective_tol = objective_tol
+
         # Configure evaluator
         # if a callable is directly passed wrap it around the serial evaluator
         self.check_evaluator(evaluator)
@@ -605,9 +639,19 @@ def num_evals():
         # Update the number of evals in case the `search.search(...)` was previously called
         max_evals = max_evals if max_evals < 0 else max_evals + num_evals()
 
+        # Stopping criteria
+        objective_tol = self._objective_tol
+
+        def stopping_criteria() -> bool:
+            return (
+                not self.stopped
+                and (max_evals < 0 or num_evals() < max_evals)
+                and not self.history.check_objective_convergence(objective_tol, window=20)
+            )
+
         n_ask = self._evaluator.num_workers
 
-        while not self.stopped and (max_evals < 0 or num_evals() < max_evals):
+        while stopping_criteria():
             new_batch = self.ask(n_ask)
 
             logger.info(f"Submitting {len(new_batch)} configurations...")
diff --git a/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_ga.py b/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_ga.py
index dfed64ae..ec64660d 100644
--- a/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_ga.py
+++ b/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_ga.py
@@ -35,6 +35,7 @@ def _evaluate(self, x, out, *args, **kwargs):
         out["F"] = y
 
         if self.constraint_fn is not None:
+            # "H" are equality constraints
             out["H"] = self.constraint_fn(x)
 
 
diff --git a/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_mixedga.py b/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_mixedga.py
index 71de6e5a..49f15c08 100644
--- a/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_mixedga.py
+++ b/src/deephyper/skopt/optimizer/acq_optimizer/pymoo_mixedga.py
@@ -1,6 +1,7 @@
 """Mixed-Integer genetic-algorithm optimization for the acquisition function."""
 
 from collections import OrderedDict
+from typing import Callable
 
 import numpy as np
 from ConfigSpace.forbidden import ForbiddenClause, ForbiddenConjunction, ForbiddenRelation
@@ -22,8 +23,8 @@
 from pymoo.termination.max_gen import MaximumGenerationTermination
 from pymoo.termination.robust import RobustTermination
 from sklearn.utils import check_random_state
-
 import deephyper.skopt.space as skopt_space
+from deephyper.skopt.space.space import Space
 
 Config.warnings["not_compiled"] = False
 
@@ -31,10 +32,7 @@
 # https://pymoo.org/interface/problem.html
 
 
-Config.warnings["not_compiled"] = False
-
-
-def convert_space_to_pymoo_mixed(space):
+def convert_space_to_pymoo_mixed(space: Space) -> OrderedDict[str, Real | Integer | Choice]:
     """Convert a DeepHyper space to a pymoo space.
 
     Optimizing in the source input space.
@@ -43,7 +41,7 @@ def convert_space_to_pymoo_mixed(space):
         space (Space): from deephyper.skopt.space.
 
     Returns:
-        dict: a pymoo space.
+        OrderedDict[str, Real | Integer | Choice]: a pymoo space.
     """
     pymoo_space = OrderedDict()
     for dim in space.dimensions:
@@ -69,7 +67,8 @@ def __init__(self, space, acq_func=None, constraint_fn=None, **kwargs):
         super().__init__(
             vars=convert_space_to_pymoo_mixed(space),
             n_obj=1,
-            n_eq_constr=int(constraint_fn is not None),
+            n_eq_constr=int(constraint_fn is not None), # H
+            # n_ieq_constr=int(constraint_fn is not None), # G
             **kwargs,
         )
         self.space = space
@@ -92,8 +91,13 @@ def _evaluate(self, x, out, *args, **kwargs):
         out["F"] = y
 
         if self.constraint_fn is not None:
+            # Equality Constraints
             out["H"] = self.constraint_fn(x)
 
+            # Inequality Constraints
+            # out["G"] = self.constraint_fn(x)
+            # print(np.asarray(x[0]).tolist(), out["F"][0], out["G"][0])
+
 
 class PyMOOMixedElementWiseProblem(ElementwiseProblem):
     """Pymoo mixed-integer problem definition (element-wise)."""
@@ -140,10 +144,11 @@ def __init__(self, ftol=1e-6, period=30, n_max_gen=1000, **kwargs) -> None:
 class ConfigSpaceRepair(Repair):
     """Pymoo repair operator for ConfigSpace conditions/forbiddens."""
 
-    def __init__(self, space):
+    def __init__(self, space: Space, repair_fn: Callable | None=None):
         super().__init__()
         self.space = space
         self.config_space = self.space.config_space
+        self.repair_fn = repair_fn
 
     def _do(self, problem, x, **kwargs):
         def deactivate_inactive_dimensions(x: dict):
@@ -192,6 +197,10 @@ def deactivate_inactive_dimensions(x: dict):
                 )
             )
 
+        # Custom repair function defined in the HpProblem
+        if self.repair_fn:
+            x = self.repair_fn(x)
+
         return x
 
 
@@ -200,13 +209,14 @@ class MixedGAPymooAcqOptimizer:
 
     def __init__(
         self,
-        space,
-        x_init,
-        y_init,
+        space: Space,
+        x_init: np.ndarray | list[list],
+        y_init: np.ndarray | list,
         pop_size: int = 100,
-        random_state=None,
-        termination_kwargs=None,
-        constraint_fn=None,
+        random_state: int | np.random.RandomState | None=None,
+        termination_kwargs: dict | None=None,
+        constraint_fn: Callable | None=None,
+        repair_fn = Callable | None,
     ):
         self.space = space
         self.x_init = np.array(x_init)
@@ -224,6 +234,7 @@ def __init__(
         default_termination_kwargs.update(termination_kwargs)
         self.termination_kwargs = default_termination_kwargs
         self.constraint_fn = constraint_fn
+        self.repair_fn = repair_fn
 
     def minimize(self, acq_func):
         """Minimize the acquisition function."""
@@ -240,7 +251,7 @@ def minimize(self, acq_func):
             self.y_init,
         )
 
-        repair = ConfigSpaceRepair(self.space)
+        repair = ConfigSpaceRepair(self.space, self.repair_fn)
         eliminate_duplicates = MixedVariableDuplicateElimination()
         algorithm = MixedVariableGA(
             pop_size=self.pop_size,
diff --git a/src/deephyper/skopt/optimizer/optimizer.py b/src/deephyper/skopt/optimizer/optimizer.py
index b596b1a5..93186988 100644
--- a/src/deephyper/skopt/optimizer/optimizer.py
+++ b/src/deephyper/skopt/optimizer/optimizer.py
@@ -1044,7 +1044,7 @@ def _tell(self, x, y, fit=True):
                     )
 
                     pop_size = self._pymoo_pop_size
-
+                    
                     idx_sorted = np.argsort(values)
                     x_init = [Xsample[i] for i in idx_sorted[:pop_size]]
                     x_init = list(
@@ -1057,21 +1057,29 @@ def _tell(self, x, y, fit=True):
 
                     # Constraint handling
                     constraint_fn = None
-                    if (
-                        self.space.custom_sampler is not None
-                        and hasattr(self.space.custom_sampler, "constraint_fn")
-                        and self.space.custom_sampler.constraint_fn
-                    ):
-                        inner_fn = self.space.custom_sampler.constraint_fn
-                        dim_names = self.space.dimension_names
-
-                        def constraint_fn(x):
-                            df = pd.DataFrame(x, columns=dim_names)
-                            accept = inner_fn(df)
-                            G = (~accept.values).astype(float).reshape(-1)
-                            return G
-
-                        constraint_fn = constraint_fn
+                    repair_fn = None
+                    if self.space.custom_sampler is not None:
+                        if hasattr(self.space.custom_sampler, "constraint_fn") and self.space.custom_sampler.constraint_fn:
+                            inner_fn = self.space.custom_sampler.constraint_fn
+                            dim_names = self.space.dimension_names
+
+                            def constraint_fn(x):
+                                df = pd.DataFrame(x, columns=dim_names)
+                                accept = inner_fn(df)
+                                if isinstance(accept[0], numbers.Number):
+                                    G = accept.values.astype(float).reshape(-1)
+                                else:
+                                    G = (~accept.values).astype(float).reshape(-1)
+                                return G
+
+                            constraint_fn = constraint_fn
+
+                        if hasattr(self.space.custom_sampler, "repair_fn") and self.space.custom_sampler.repair_fn:
+                            def repair_fn(x):
+                                df = pd.DataFrame(map(lambda xi: xi, x), columns=dim_names)
+                                df = self.space.custom_sampler.repair_fn(df)
+                                df = np.asarray(df.to_dict(orient="records"))
+                                return df
 
                     acq_opt = MixedGAPymooAcqOptimizer(
                         space=self.space,
@@ -1082,6 +1090,7 @@ def constraint_fn(x):
                         random_state=self.rng.randint(0, np.iinfo(np.int32).max),
                         termination_kwargs=self._pymoo_termination_kwargs,
                         constraint_fn=constraint_fn,
+                        repair_fn=repair_fn,
                     )
 
                     args = (est, np.min(yi), cand_acq_func, False, self.acq_func_kwargs)
@@ -1109,7 +1118,14 @@ def constraint_fn(x):
                         def constraint_fn(x):
                             df = pd.DataFrame(x, columns=dim_names)
                             accept = inner_fn(df)
-                            G = (~accept.values).astype(float).reshape(-1)
+                            # Check if we directly have constraints function values or just booleans
+                            if isinstance(accept[0], numbers.Number):
+                                G = accept.values.astype(float).reshape(-1)
+                            else:
+                                # The user defined a function returning a boolean
+                                # True: The constraint is valid
+                                # False: The constraint is invalid
+                                G = (~accept.values).astype(float).reshape(-1)
                             return G
 
                         constraint_fn = constraint_fn
diff --git a/tests/hpo/test_hp_problem.py b/tests/hpo/test_hp_problem.py
index cd53eff2..17f80f52 100644
--- a/tests/hpo/test_hp_problem.py
+++ b/tests/hpo/test_hp_problem.py
@@ -173,7 +173,7 @@ def constraint_fn(df: pd.DataFrame) -> pd.Series:
 
     pb.set_constraint_fn(constraint_fn)
 
-    samples = pb.sample(size=100)
+    samples = pb.sample(size=100, strict=False)
     df = pd.DataFrame(samples)
     assert all(df["x"] >= 9)