diff --git a/imbens/datasets/_openml.py b/imbens/datasets/_openml.py
index 845b08f..761cd83 100644
--- a/imbens/datasets/_openml.py
+++ b/imbens/datasets/_openml.py
@@ -21,6 +21,95 @@
 import tqdm
 
 
+def _get_data_path(data_home):
+    if data_home is None:
+        data_home = Path(user_cache_dir("imbens")) / "datasets"
+        data_home.mkdir(parents=True, exist_ok=True)
+        return data_home
+    if not isinstance(data_home, Path):
+        data_home = Path(data_home)
+    if not data_home.is_dir():
+        raise ValueError(f"data_home {data_home} is not a directory")
+    data_home.mkdir(parents=True, exist_ok=True)
+    return data_home
+
+
+def _save_data(X, y, data_home, file_name):
+    np.savez_compressed(data_home / file_name, X=X, y=y)
+
+
+def _load_data(data_home, file_name):
+    data = np.load(data_home / file_name, allow_pickle=True)
+    return data["X"], data["y"]
+
+
+def _validate_openml_params(openml_id, cat_preprocess, data_home):
+    assert isinstance(openml_id, int), "openml_id must be an integer"
+    assert cat_preprocess in [
+        "drop",
+        "onehot",
+        "ordinal",
+    ], "cat_preprocess must be one of ['drop', 'onehot', 'ordinal']"
+    assert isinstance(
+        data_home, (str, Path, type(None))
+    ), "data_home must be a string, Path or None"
+
+
+def _encode_target(y):
+    y_map = y.value_counts().sort_values()[::-1].index
+    y_map = {v: i for i, v in enumerate(y_map)}
+    return y.map(y_map)
+
+
+def _standardize_numeric(X, feats_type):
+    scaler = StandardScaler()
+    X[feats_type["num"]] = scaler.fit_transform(X[feats_type["num"]])
+    return X
+
+
+def _preprocess_feature(feat, X, feats_type):
+    if is_any_real_numeric_dtype(X[feat]):
+        feats_type["num"].append(feat)
+        if X[feat].isnull().sum() > 0:
+            X[feat] = X[feat].fillna(X[feat].mean())
+    else:
+        if X[feat].isnull().sum() > 0:
+            X[feat] = X[feat].fillna(X[feat].mode().iloc[0])
+        n_unique = len(X[feat].unique())
+        if n_unique > 2:
+            try:
+                X[feat] = pd.to_numeric(X[feat])
+            except ValueError:
+                if n_unique <= 50:
+                    feats_type["multi_cat"].append(feat)
+                else:
+                    feats_type["drop"].append(feat)
+        else:
+            feats_type["bin_cat"].append(feat)
+
+
+def _cat_drop(X, feats_type, ord_encoder):
+    return X.drop(columns=feats_type["multi_cat"] + feats_type["bin_cat"])
+
+
+def _cat_onehot(X, feats_type, ord_encoder):
+    X[feats_type["bin_cat"]] = ord_encoder.fit_transform(X[feats_type["bin_cat"]])
+    return pd.get_dummies(X, columns=feats_type["multi_cat"])
+
+
+def _cat_ordinal(X, feats_type, ord_encoder):
+    X[feats_type["bin_cat"]] = ord_encoder.fit_transform(X[feats_type["bin_cat"]])
+    X[feats_type["multi_cat"]] = ord_encoder.fit_transform(X[feats_type["multi_cat"]])
+    return X
+
+
+_CAT_PREPROCESSORS = {
+    "drop": _cat_drop,
+    "onehot": _cat_onehot,
+    "ordinal": _cat_ordinal,
+}
+
+
 def _fetch_openml_data(openml_id, cat_preprocess="onehot", data_home=None):
     """
     Fetches a dataset from OpenML, preprocesses it, and caches it locally.
@@ -52,112 +141,35 @@ def _fetch_openml_data(openml_id, cat_preprocess="onehot", data_home=None):
     If the dataset is already cached locally, it will be loaded from the cache.
     Otherwise, the dataset will be downloaded, preprocessed, and saved locally.
     """
+    _validate_openml_params(openml_id, cat_preprocess, data_home)
 
-    def get_data_path(data_home):
-        if data_home is None:
-            data_home = Path(user_cache_dir("imbens")) / "datasets"
-            data_home.mkdir(parents=True, exist_ok=True)
-            return data_home
-        else:
-            # check if store_path is a valid path
-            if not isinstance(data_home, Path):
-                data_home = Path(data_home)
-            if not data_home.is_dir():
-                raise ValueError(f"data_home {data_home} is not a directory")
-            data_home.mkdir(parents=True, exist_ok=True)
-            return data_home
-
-    def save_data(X, y, data_home, file_name):
-        np.savez_compressed(data_home / file_name, X=X, y=y)
-        # print(f"Data saved to {data_home / file_name}")
-
-    def load_data(data_home, file_name):
-        data = np.load(data_home / file_name, allow_pickle=True)
-        X = data["X"]
-        y = data["y"]
-        return X, y
-
-    assert isinstance(openml_id, int), "openml_id must be an integer"
-    assert cat_preprocess in [
-        "drop",
-        "onehot",
-        "ordinal",
-    ], "cat_preprocess must be one of ['drop', 'onehot', 'ordinal']"
-    assert isinstance(
-        data_home, (str, Path, type(None))
-    ), "data_home must be a string, Path or None"
-
-    data_home = get_data_path(data_home)
+    data_home = _get_data_path(data_home)
     dataset = openml.datasets.get_dataset(openml_id)
     file_name = f"{dataset.id}_{cat_preprocess}_{dataset.name}.npz"
 
-    # check if data is already cached
     try:
-        X, y = load_data(data_home, file_name)
-        # print(f"Data loaded from {data_home / file_name}")
+        X, y = _load_data(data_home, file_name)
         return X, y
     except FileNotFoundError:
-        # print(f"Data not cached in {data_home}, processing from scratch.")
         pass
 
     X, y, cat_ind, feat_names = dataset.get_data(
         target=dataset.default_target_attribute
     )
 
-    # target encoding, smaller classes aer assigned larger values
-    y_map = y.value_counts().sort_values()[::-1].index
-    y_map = {v: i for i, v in enumerate(y_map)}
-    y = y.map(y_map)
+    y = _encode_target(y)
 
     feats_type = {"num": [], "bin_cat": [], "multi_cat": [], "drop": []}
-    # preprocessing
     for feat in X.columns:
-        if is_any_real_numeric_dtype(X[feat]):
-            feats_type["num"].append(feat)
-            if X[feat].isnull().sum() > 0:
-                # for numerical columns, fill nan with mean
-                X[feat] = X[feat].fillna(X[feat].mean())
-        else:  # categorical column
-            if X[feat].isnull().sum() > 0:
-                # for categorical columns, fill nan with most frequent value
-                X[feat] = X[feat].fillna(X[feat].mode().iloc[0])
-            n_unique = len(X[feat].unique())
-            if n_unique > 2:
-                # try to convert to numeric
-                try:
-                    X[feat] = pd.to_numeric(X[feat])
-                except ValueError:
-                    if n_unique <= 50:
-                        feats_type["multi_cat"].append(feat)
-                    else:
-                        feats_type["drop"].append(feat)
-            else:
-                feats_type["bin_cat"].append(feat)
+        _preprocess_feature(feat, X, feats_type)
 
     ord_encoder = OrdinalEncoder()
     X = X.drop(columns=feats_type["drop"])
-    # encode categorical columns
-    if cat_preprocess == "drop":
-        X = X.drop(columns=feats_type["multi_cat"])
-        X = X.drop(columns=feats_type["bin_cat"])
-    elif cat_preprocess == "onehot":
-        X[feats_type["bin_cat"]] = ord_encoder.fit_transform(X[feats_type["bin_cat"]])
-        X = pd.get_dummies(X, columns=feats_type["multi_cat"])
-    elif cat_preprocess == "ordinal":
-        # ordinal encoding for multi categorical columns
-        X[feats_type["bin_cat"]] = ord_encoder.fit_transform(X[feats_type["bin_cat"]])
-        X[feats_type["multi_cat"]] = ord_encoder.fit_transform(
-            X[feats_type["multi_cat"]]
-        )
-    else:
-        raise ValueError(f"Unknown cat_preprocess: {cat_preprocess}")
+    X = _CAT_PREPROCESSORS[cat_preprocess](X, feats_type, ord_encoder)
 
-    # standardize numerical columns
-    scaler = StandardScaler()
-    X[feats_type["num"]] = scaler.fit_transform(X[feats_type["num"]])
+    X = _standardize_numeric(X, feats_type)
 
-    # save data
-    save_data(X, y, data_home, file_name)
+    _save_data(X, y, data_home, file_name)
 
     return X, y
 
diff --git a/imbens/datasets/_zenodo.py b/imbens/datasets/_zenodo.py
index 1bac6a5..cf4ee6f 100644
--- a/imbens/datasets/_zenodo.py
+++ b/imbens/datasets/_zenodo.py
@@ -100,6 +100,57 @@
     MAP_ID_NAME[v + 1] = k
 
 
+def _resolve_filter_data(filter_data):
+    """Validate and resolve filter_data into a list of dataset names."""
+    if filter_data is None:
+        return list(MAP_NAME_ID.keys())
+
+    list_data = MAP_NAME_ID.keys()
+    filter_data_ = []
+    for it in filter_data:
+        if isinstance(it, str):
+            if it not in list_data:
+                raise ValueError(
+                    f"{it} is not a dataset available. "
+                    f"The available datasets are {list_data}"
+                )
+            filter_data_.append(it)
+        elif isinstance(it, int):
+            if it < 1 or it > 27:
+                raise ValueError(
+                    f"The dataset with the ID={it} is not an "
+                    f"available dataset. The IDs are "
+                    f"{range(1, 28)}"
+                )
+            filter_data_.append(MAP_ID_NAME[it])
+        else:
+            raise ValueError(
+                f"The value in the tuple should be str or int."
+                f" Got {type(it)} instead."
+            )
+    return filter_data_
+
+
+def _load_dataset(dataset_name, zenodo_dir, download_if_missing, verbose):
+    """Load a single dataset from disk, downloading if necessary."""
+    filename = PRE_FILENAME + str(MAP_NAME_ID[dataset_name]) + POST_FILENAME
+    filepath = join(zenodo_dir, filename)
+    available = isfile(filepath)
+
+    if download_if_missing and not available:
+        makedirs(zenodo_dir, exist_ok=True)
+        if verbose:
+            print("Downloading %s" % URL)
+        f = BytesIO(urlopen(URL).read())
+        tar = tarfile.open(fileobj=f)
+        tar.extractall(path=zenodo_dir)
+    elif not download_if_missing and not available:
+        raise IOError("Data not found and `download_if_missing` is False")
+
+    data = np.load(filepath)
+    return data["data"], data["label"]
+
+
 @_deprecate_positional_args
 def fetch_zenodo_datasets(
     *,
@@ -165,11 +216,11 @@ def fetch_zenodo_datasets(
     +--+--------------+-------------------------------+-------+---------+-----+
     |2 |optical_digits| UCI, target: 8                | 9.1:1 | 5,620   | 64  |
     +--+--------------+-------------------------------+-------+---------+-----+
-    |3 |satimage      | UCI, target: 4                | 9.3:1 | 6,435   | 36  |
+    |3 |satimage      | UCI, target: 4                | 9.4:1 | 6,435   | 36  |
     +--+--------------+-------------------------------+-------+---------+-----+
     |4 |pen_digits    | UCI, target: 5                | 9.4:1 | 10,992  | 16  |
     +--+--------------+-------------------------------+-------+---------+-----+
-    |5 |abalone       | UCI, target: 7                | 9.7:1 | 4,177   | 10  |
+    |5 |abalone       | UCI, target: 7                | 9.4:1 | 4,177   | 10  |
     +--+--------------+-------------------------------+-------+---------+-----+
     |6 |sick_euthyroid| UCI, target: sick euthyroid   | 9.8:1 | 3,163   | 42  |
     +--+--------------+-------------------------------+-------+---------+-----+
@@ -227,55 +278,10 @@ def fetch_zenodo_datasets(
     zenodo_dir = join(data_home, "zenodo")
     datasets = OrderedDict()
 
-    if filter_data is None:
-        filter_data_ = MAP_NAME_ID.keys()
-    else:
-        list_data = MAP_NAME_ID.keys()
-        filter_data_ = []
-        for it in filter_data:
-            if isinstance(it, str):
-                if it not in list_data:
-                    raise ValueError(
-                        f"{it} is not a dataset available. "
-                        f"The available datasets are {list_data}"
-                    )
-                else:
-                    filter_data_.append(it)
-            elif isinstance(it, int):
-                if it < 1 or it > 27:
-                    raise ValueError(
-                        f"The dataset with the ID={it} is not an "
-                        f"available dataset. The IDs are "
-                        f"{range(1, 28)}"
-                    )
-                else:
-                    # The index start at one, then we need to remove one
-                    # to not have issue with the indexing.
-                    filter_data_.append(MAP_ID_NAME[it])
-            else:
-                raise ValueError(
-                    f"The value in the tuple should be str or int."
-                    f" Got {type(it)} instead."
-                )
+    filter_data_ = _resolve_filter_data(filter_data)
 
-    # go through the list and check if the data are available
     for it in filter_data_:
-        filename = PRE_FILENAME + str(MAP_NAME_ID[it]) + POST_FILENAME
-        filename = join(zenodo_dir, filename)
-        available = isfile(filename)
-
-        if download_if_missing and not available:
-            makedirs(zenodo_dir, exist_ok=True)
-            if verbose:
-                print("Downloading %s" % URL)
-            f = BytesIO(urlopen(URL).read())
-            tar = tarfile.open(fileobj=f)
-            tar.extractall(path=zenodo_dir)
-        elif not download_if_missing and not available:
-            raise IOError("Data not found and `download_if_missing` is False")
-
-        data = np.load(filename)
-        X, y = data["data"], data["label"]
+        X, y = _load_dataset(it, zenodo_dir, download_if_missing, verbose)
 
         if shuffle:
             ind = np.arange(X.shape[0])
diff --git a/imbens/ensemble/_bagging.py b/imbens/ensemble/_bagging.py
index 0bacb6d..615bf44 100644
--- a/imbens/ensemble/_bagging.py
+++ b/imbens/ensemble/_bagging.py
@@ -55,6 +55,7 @@
 import itertools
 import numbers
 from abc import ABCMeta, abstractmethod
+from types import SimpleNamespace
 from warnings import warn
 
 import numpy as np
@@ -200,6 +201,69 @@ class ResampleBaggingClassifier(
         "training_type": _training_type,
     }
 
+    # --- Properties for backward compatibility ---
+    # (store data internally in _sampler_cfg / _train_state)
+
+    @property
+    def sampling_strategy(self):
+        return self._sampler_cfg.sampling_strategy
+
+    @sampling_strategy.setter
+    def sampling_strategy(self, value):
+        self._sampler_cfg.sampling_strategy = value
+
+    @property
+    def _sampling_type(self):
+        return self._sampler_cfg.sampling_type
+
+    @_sampling_type.setter
+    def _sampling_type(self, value):
+        self._sampler_cfg.sampling_type = value
+
+    @property
+    def _max_samples(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.max_samples
+        return None
+
+    @property
+    def _max_features(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.max_features
+        return None
+
+    @property
+    def _n_samples(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.n_samples
+        return None
+
+    @property
+    def _seeds(self):
+        return self._train_state.seeds
+
+    @_seeds.setter
+    def _seeds(self, value):
+        self._train_state.seeds = value
+
+    @property
+    def eval_datasets_(self):
+        if not hasattr(self, '_train_state'):
+            return {}
+        return self._train_state.eval_datasets
+
+    @property
+    def eval_metrics_(self):
+        if not hasattr(self, '_train_state'):
+            return {}
+        return self._train_state.eval_metrics
+
+    @property
+    def train_verbose_(self):
+        if not hasattr(self, '_train_state'):
+            return None
+        return self._train_state.train_verbose
+
     @_deprecate_positional_args
     def __init__(
         self,
@@ -220,9 +284,11 @@ def __init__(
         verbose=0,
     ):
 
-        self.sampling_strategy = sampling_strategy
-        self._sampling_type = sampling_type
-        self.sampler = sampler
+        self._sampler_cfg = SimpleNamespace(
+            sampler=sampler,
+            sampling_type=sampling_type,
+            sampling_strategy=sampling_strategy,
+        )
 
         super().__init__(
             estimator=estimator,
@@ -242,19 +308,19 @@ def _validate_y(self, y):
         """Validate the label vector."""
         y_encoded = super()._validate_y(y)
         if (
-            isinstance(self.sampling_strategy, dict)
-            and self.sampler_._sampling_type != "bypass"
+            isinstance(self._sampler_cfg.sampling_strategy, dict)
+            and self._sampler_cfg.sampler_._sampling_type != "bypass"
         ):
-            self._sampling_strategy = {
+            self._sampler_cfg.internal_sampling_strategy = {
                 np.where(self.classes_ == key)[0][0]: value
                 for key, value in check_sampling_strategy(
-                    self.sampling_strategy,
+                    self._sampler_cfg.sampling_strategy,
                     y,
-                    self.sampler_._sampling_type,
+                    self._sampler_cfg.sampler_._sampling_type,
                 ).items()
             }
         else:
-            self._sampling_strategy = self.sampling_strategy
+            self._sampler_cfg.internal_sampling_strategy = self._sampler_cfg.sampling_strategy
         return y_encoded
 
     def _validate_estimator(self, default=DecisionTreeClassifier()):
@@ -276,9 +342,9 @@ def _validate_estimator(self, default=DecisionTreeClassifier()):
             estimator = clone(default)
 
         # validate sampler and sampler_kwargs
-        # validated sampler stored in self.sampler_
+        # validated sampler stored in _sampler_cfg.sampler_
         try:
-            self.sampler_ = clone(self.sampler)
+            self._sampler_cfg.sampler_ = clone(self._sampler_cfg.sampler)
         except Exception as e:
             e_args = list(e.args)
             e_args[0] = (
@@ -287,12 +353,14 @@ def _validate_estimator(self, default=DecisionTreeClassifier()):
             e.args = tuple(e_args)
             raise e
 
-        if self.sampler_._sampling_type != "bypass":
-            self.sampler_.set_params(sampling_strategy=self._sampling_strategy)
-            self.sampler_.set_params(**self.sampler_kwargs_)
+        if self._sampler_cfg.sampler_._sampling_type != "bypass":
+            self._sampler_cfg.sampler_.set_params(
+                sampling_strategy=self._sampler_cfg.internal_sampling_strategy
+            )
+            self._sampler_cfg.sampler_.set_params(**self._sampler_cfg.sampler_kwargs)
 
         self._estimator = Pipeline(
-            [("sampler", self.sampler_), ("classifier", estimator)]
+            [("sampler", self._sampler_cfg.sampler_), ("classifier", estimator)]
         )
         try:
             # scikit-learn < 1.2
@@ -300,67 +368,18 @@ def _validate_estimator(self, default=DecisionTreeClassifier()):
         except AttributeError:
             pass
 
-    @_deprecate_positional_args
-    @FuncSubstitution(
-        eval_datasets=_get_parameter_docstring("eval_datasets"),
-        eval_metrics=_get_parameter_docstring("eval_metrics"),
-        train_verbose=_get_parameter_docstring("train_verbose", **_properties),
-    )
-    def _fit(
-        self,
-        X,
-        y,
-        *,
-        sample_weight=None,
-        sampler_kwargs: dict = {},
-        max_samples=None,
-        eval_datasets: dict = None,
-        eval_metrics: dict = None,
-        train_verbose: bool or int or dict,
+    def _validate_input_data(
+        self, X, y, sample_weight, sampler_kwargs, eval_datasets, eval_metrics, train_verbose
     ):
-        """Build a Bagging ensemble of estimators from the training set (X, y).
-
-        Parameters
-        ----------
-        X : {array-like, sparse matrix} of shape (n_samples, n_features)
-            The training input samples. Sparse matrices are accepted only if
-            they are supported by the base estimator.
-
-        y : array-like of shape (n_samples,)
-            The target values (class labels in classification, real numbers in
-            regression).
-
-        sample_weight : array-like of shape (n_samples,), default=None
-            Sample weights. If None, then samples are equally weighted.
-            Note that this is supported only if the base estimator supports
-            sample weighting.
-
-        sampler_kwargs : dict, default={}
-            The kwargs to use as additional parameters when instantiating a
-            new sampler. If none are given, default parameters are used.
-
-        max_samples : int or float, default=None
-            Argument to use instead of self.max_samples.
-
-        %(eval_datasets)s
-
-        %(eval_metrics)s
-
-        %(train_verbose)s
-
-        Returns
-        -------
-        self : object
-        """
-
-        # Check data, sampler_kwargs and random_state
         check_target_type(y)
 
-        self.sampler_kwargs_ = check_type(sampler_kwargs, "sampler_kwargs", dict)
+        self._train_state = SimpleNamespace()
+        self._sampler_cfg.sampler_kwargs = check_type(
+            sampler_kwargs, "sampler_kwargs", dict
+        )
 
         random_state = check_random_state(self.random_state)
 
-        # Convert data (X is required to be 2d and indexable)
         check_x_y_args = {
             "accept_sparse": ["csr", "csc"],
             "dtype": None,
@@ -369,14 +388,13 @@ def _fit(
         }
         X, y = validate_data(self, X, y, **check_x_y_args)
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        self._train_state.eval_datasets = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
 
-        # Check evaluation metrics
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        self._train_state.eval_metrics = check_eval_metrics(eval_metrics)
 
-        # Check verbose
-        self.train_verbose_ = check_train_verbose(
+        self._train_state.train_verbose = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
         self._init_training_log_format()
@@ -384,15 +402,15 @@ def _fit(
         if sample_weight is not None:
             sample_weight = _check_sample_weight(sample_weight, X, dtype=None)
 
-        # Remap output
-        n_samples, self.n_features_in_ = X.shape
-        self._n_samples = n_samples
+        n_samples = X.shape[0]
+        self._train_state.n_samples = n_samples
         y = self._validate_y(y)
 
-        # Check parameters
         self._validate_estimator()
 
-        # Validate max_samples
+        return X, y, random_state, sample_weight
+
+    def _validate_params(self, X, max_samples):
         if max_samples is None:
             max_samples = self.max_samples
         if not isinstance(max_samples, numbers.Integral):
@@ -401,10 +419,8 @@ def _fit(
         if not (0 < max_samples <= X.shape[0]):
             raise ValueError("max_samples must be in (0, n_samples]")
 
-        # Store validated integer row sampling value
-        self._max_samples = max_samples
+        self._train_state.max_samples = max_samples
 
-        # Validate max_features
         if isinstance(self.max_features, numbers.Integral):
             max_features = self.max_features
         elif isinstance(self.max_features, float):
@@ -416,11 +432,9 @@ def _fit(
             raise ValueError("max_features must be in (0, n_features]")
 
         max_features = max(1, int(max_features))
+        self._train_state.max_features = max_features
 
-        # Store validated integer feature sampling value
-        self._max_features = max_features
-
-        # Other checks
+    def _check_bootstrap_oob_warm_start(self):
         if not self.bootstrap and self.oob_score:
             raise ValueError(
                 "Out of bag estimation only available" " if bootstrap=True"
@@ -435,7 +449,6 @@ def _fit(
             del self.oob_score_
 
         if not self.warm_start or not hasattr(self, "estimators_"):
-            # Free allocated memory, if any
             self.estimators_ = []
             self.estimators_features_ = []
             self.estimators_n_training_samples_ = []
@@ -449,26 +462,21 @@ def _fit(
                 % (self.n_estimators, len(self.estimators_))
             )
 
-        elif n_more_estimators == 0:
-            warn(
-                "Warm-start fitting without increasing n_estimators does not "
-                "fit new trees."
-            )
-            return self
+        return n_more_estimators
 
-        # Parallel loop
+    def _parallel_fit_estimators(
+        self, X, y, sample_weight, n_more_estimators, random_state
+    ):
         n_jobs, n_estimators, starts = _partition_estimators(
             n_more_estimators, self.n_jobs
         )
         total_n_estimators = sum(n_estimators)
 
-        # Advance random state to state after training
-        # the first n_estimators
         if self.warm_start and len(self.estimators_) > 0:
             random_state.randint(MAX_INT, size=len(self.estimators_))
 
         seeds = random_state.randint(MAX_INT, size=n_more_estimators)
-        self._seeds = seeds
+        self._train_state.seeds = seeds
 
         all_results = Parallel(
             n_jobs=n_jobs, verbose=self.verbose, **self._parallel_args()
@@ -486,7 +494,6 @@ def _fit(
             for i in range(n_jobs)
         )
 
-        # Reduce
         self.estimators_ += list(
             itertools.chain.from_iterable(t[0] for t in all_results)
         )
@@ -497,10 +504,78 @@ def _fit(
             itertools.chain.from_iterable(t[2] for t in all_results)
         )
 
+    @_deprecate_positional_args
+    @FuncSubstitution(
+        eval_datasets=_get_parameter_docstring("eval_datasets"),
+        eval_metrics=_get_parameter_docstring("eval_metrics"),
+        train_verbose=_get_parameter_docstring("train_verbose", **_properties),
+    )
+    def _fit(
+        self,
+        X,
+        y,
+        *,
+        sample_weight=None,
+        sampler_kwargs: dict = {},
+        max_samples=None,
+        eval_datasets: dict = None,
+        eval_metrics: dict = None,
+        train_verbose: bool or int or dict,
+    ):
+        """Build a Bagging ensemble of estimators from the training set (X, y).
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix} of shape (n_samples, n_features)
+            The training input samples. Sparse matrices are accepted only if
+            they are supported by the base estimator.
+
+        y : array-like of shape (n_samples,)
+            The target values (class labels in classification, real numbers in
+            regression).
+
+        sample_weight : array-like of shape (n_samples,), default=None
+            Sample weights. If None, then samples are equally weighted.
+            Note that this is supported only if the base estimator supports
+            sample weighting.
+
+        sampler_kwargs : dict, default={}
+            The kwargs to use as additional parameters when instantiating a
+            new sampler. If none are given, default parameters are used.
+
+        max_samples : int or float, default=None
+            Argument to use instead of self.max_samples.
+
+        %(eval_datasets)s
+
+        %(eval_metrics)s
+
+        %(train_verbose)s
+
+        Returns
+        -------
+        self : object
+        """
+
+        X, y, random_state, sample_weight = self._validate_input_data(
+            X, y, sample_weight, sampler_kwargs, eval_datasets, eval_metrics,
+            train_verbose
+        )
+        self._validate_params(X, max_samples)
+        n_more_estimators = self._check_bootstrap_oob_warm_start()
+        if n_more_estimators == 0:
+            warn(
+                "Warm-start fitting without increasing n_estimators does not "
+                "fit new trees."
+            )
+            return self
+        self._parallel_fit_estimators(
+            X, y, sample_weight, n_more_estimators, random_state
+        )
+
         if self.oob_score:
             self._set_oob_score(X, y)
 
-        # Print training infomation to console.
         self._training_log_to_console()
 
         return self
diff --git a/imbens/ensemble/_boost.py b/imbens/ensemble/_boost.py
index 4a86c15..74b13cd 100644
--- a/imbens/ensemble/_boost.py
+++ b/imbens/ensemble/_boost.py
@@ -24,12 +24,12 @@
         check_train_verbose,
         check_type,
     )
-    from .base import MAX_INT, ImbalancedEnsembleClassifierMixin
+    from .base import MAX_INT, ImbalancedEnsembleClassifierMixin, _TrainingState
 else:  # pragma: no cover
     import sys  # For local test
 
     sys.path.append("..")
-    from ensemble.base import ImbalancedEnsembleClassifierMixin, MAX_INT
+    from ensemble.base import ImbalancedEnsembleClassifierMixin, MAX_INT, _TrainingState
     from utils._docstring import FuncGlossarySubstitution
     from utils._validation import _deprecate_positional_args
     from utils._validation_data import check_eval_datasets
@@ -58,7 +58,7 @@
 _super = AdaBoostClassifier
 
 
-class ResampleBoostClassifier(
+class ResampleBoostClassifier(  # pylint: disable=too-many-instance-attributes
     ImbalancedEnsembleClassifierMixin, AdaBoostClassifier, metaclass=ABCMeta
 ):
     """Base class for all resampling + boosting imbalanced ensemble classifier.
@@ -131,10 +131,10 @@ def _make_sampler(self, append=True, random_state=None, **overwrite_kwargs):
         """
 
         sampler = clone(self.sampler_)
-        sampler.set_params(**self.sampler_kwargs_)
+        sampler.set_params(**self._train_state_.sampler_kwargs_)
 
         # Arguments passed to _make_sampler function have higher priority,
-        # they will overwrite the self.sampler_kwargs_
+        # they will overwrite the self._train_state_.sampler_kwargs_
         sampler.set_params(**overwrite_kwargs)
 
         if random_state is not None:
@@ -346,51 +346,44 @@ def _boost_discrete(
 
         return sample_weight, estimator_weight, estimator_error
 
-    @_deprecate_positional_args
-    def _fit(
-        self,
-        X,
-        y,
-        *,
-        sample_weight,
-        sampler_kwargs: dict,
-        target_label: int,
-        n_target_samples: int or dict,
-        balancing_schedule: str or function = "uniform",
-        update_x_y_after_resample: bool = False,
-        eval_datasets: dict,
-        eval_metrics: dict,
-        train_verbose: bool or int or dict,
-    ):
-
-        update_x_y_after_resample = check_type(
-            update_x_y_after_resample, "update_x_y_after_resample", bool
-        )
-
-        self.sampler_kwargs_ = check_type(sampler_kwargs, "sampler_kwargs", dict)
+    def _check_early_stop(self, iboost, sample_weight, estimator_error, early_termination_):
+        if sample_weight is None and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight is None)."
+            )
+            return True
 
-        early_termination_ = check_type(
-            self.early_termination, "early_termination", bool
-        )
+        if estimator_error == 0 and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (training error is 0)."
+            )
+            return True
 
-        # Check that algorithm is supported.
-        if self.algorithm not in ("SAMME", "SAMME.R"):
-            raise ValueError("algorithm %s is not supported" % self.algorithm)
+        if early_termination_ and np.sum(sample_weight) <= 0:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight_sum <= 0)."
+            )
+            return True
 
-        # Check parameters.
-        if self.learning_rate <= 0:
-            raise ValueError("learning_rate must be greater than zero")
+        return False
 
-        if self.estimator == None or isinstance(
+    def _get_data_type_config(self):
+        if self.estimator is None or isinstance(
             self.estimator, (BaseDecisionTree, BaseForest)
         ):
-            DTYPE = np.float64  # from fast_dict.pxd
-            dtype = DTYPE
-            accept_sparse = "csc"
-        else:
-            dtype = None
-            accept_sparse = ["csr", "csc"]
+            return np.float64, "csc"
+        return None, ["csr", "csc"]
 
+    def _validate_and_prepare_data(
+        self, X, y, dtype, accept_sparse, eval_datasets, eval_metrics, train_verbose,
+        target_label, n_target_samples, balancing_schedule,
+    ):
         check_x_y_args = {
             "accept_sparse": accept_sparse,
             "ensure_2d": True,
@@ -400,67 +393,113 @@ def _fit(
         }
         X, y = validate_data(self, X, y, **check_x_y_args)
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        self._train_state_.eval_datasets_ = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
 
-        self.classes_, self._y_encoded = np.unique(y, return_inverse=True)
+        self.classes_, self._train_state_._y_encoded = np.unique(
+            y, return_inverse=True
+        )
         self.n_classes_ = len(self.classes_)
 
-        # Store original class distribution
-        self.origin_distr_ = dict(Counter(y))
+        self._train_state_.origin_distr_ = dict(Counter(y))
         (
             self.target_label_,
-            self.target_distr_,
+            self._train_state_.target_distr_,
         ) = check_target_label_and_n_target_samples(
             y, target_label, n_target_samples, self._sampling_type
         )
 
-        self.balancing_schedule_ = check_balancing_schedule(balancing_schedule)
+        self._train_state_.balancing_schedule_ = check_balancing_schedule(
+            balancing_schedule
+        )
 
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        self._train_state_.eval_metrics_ = check_eval_metrics(eval_metrics)
 
-        self.train_verbose_ = check_train_verbose(
+        self._train_state_.train_verbose_ = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
 
         self._init_training_log_format()
 
-        # Check sample weight
+        return X, y
+
+    def _init_sample_weight(self, sample_weight, X):
         sample_weight = _check_sample_weight(sample_weight, X, np.float64)
         sample_weight /= sample_weight.sum()
         if np.any(sample_weight < 0):
             raise ValueError("sample_weight cannot contain negative weights")
 
-        self.raw_sample_weight_ = sample_weight
+        self._train_state_.raw_sample_weight_ = sample_weight
 
-        sample_weight = copy(self.raw_sample_weight_)
+        return copy(sample_weight)
 
+    def _init_ensemble(self, sample_weight):
         self._validate_estimator()
 
-        # Check random state
         random_state = check_random_state(self.random_state)
 
-        # Clear any previous fit results.
         self.estimators_ = []
         self.estimator_weights_ = np.zeros(self.n_estimators, dtype=np.float64)
         self.estimator_errors_ = np.ones(self.n_estimators, dtype=np.float64)
-        self.estimators_n_training_samples_ = np.zeros(self.n_estimators, dtype=int)
+        self._train_state_.estimators_n_training_samples_ = np.zeros(
+            self.n_estimators, dtype=int
+        )
 
         self.samplers_ = []
 
-        # Genrate random seeds array
         seeds = random_state.randint(MAX_INT, size=self.n_estimators)
-        self._seeds = seeds
+        self._train_state_._seeds = seeds
 
         epsilon = np.finfo(sample_weight.dtype).eps
-        zero_weight_mask = sample_weight == 0.0
+        self._train_state_._epsilon = epsilon
+
+        return random_state, seeds
+
+    @_deprecate_positional_args
+    def _fit(
+        self,
+        X,
+        y,
+        *,
+        sample_weight,
+        sampler_kwargs: dict,
+        target_label: int,
+        n_target_samples: int or dict,
+        balancing_schedule: str or function = "uniform",
+        update_x_y_after_resample: bool = False,
+        eval_datasets: dict,
+        eval_metrics: dict,
+        train_verbose: bool or int or dict,
+    ):
+
+        update_x_y_after_resample = check_type(
+            update_x_y_after_resample, "update_x_y_after_resample", bool
+        )
+
+        self._train_state_ = _TrainingState()
+        self._train_state_.sampler_kwargs_ = check_type(
+            sampler_kwargs, "sampler_kwargs", dict
+        )
+
+        early_termination_ = check_type(
+            self.early_termination, "early_termination", bool
+        )
+
+        dtype, accept_sparse = self._get_data_type_config()
+        X, y = self._validate_and_prepare_data(
+            X, y, dtype, accept_sparse, eval_datasets, eval_metrics, train_verbose,
+            target_label, n_target_samples, balancing_schedule,
+        )
+        sample_weight = self._init_sample_weight(sample_weight, X)
+        random_state, seeds = self._init_ensemble(sample_weight)
         sampler_ = self.sampler_
 
         for iboost in range(self.n_estimators):
 
-            current_iter_distr = self.balancing_schedule_(
-                origin_distr=self.origin_distr_,
-                target_distr=self.target_distr_,
+            current_iter_distr = self._train_state_.balancing_schedule_(
+                origin_distr=self._train_state_.origin_distr_,
+                target_distr=self._train_state_.target_distr_,
                 i_estimator=iboost,
                 total_estimator=self.n_estimators,
             )
@@ -471,17 +510,14 @@ def _fit(
                 sampling_strategy=current_iter_distr,
             )
 
-            # Perform re-sampling
             X_resampled, y_resampled, sample_weight_resampled = sampler.fit_resample(
                 X, y, sample_weight=sample_weight
             )
 
-            # Update X, y, sample_weight if update_x_y_after_resample is True
             if update_x_y_after_resample:
                 X, y = X_resampled, y_resampled
                 sample_weight = sample_weight_resampled
 
-            # Boosting step.
             sample_weight, estimator_weight, estimator_error = self._boost(
                 iboost,
                 X_resampled,
@@ -495,46 +531,16 @@ def _fit(
 
             self.estimator_weights_[iboost] = estimator_weight
             self.estimator_errors_[iboost] = estimator_error
-            self.estimators_n_training_samples_[iboost] = y_resampled.shape[0]
+            self._train_state_.estimators_n_training_samples_[iboost] = y_resampled.shape[0]
 
-            # Print training infomation to console.
             self._training_log_to_console(iboost, y_resampled)
 
-            # Early termination.
-            if sample_weight is None and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight is None)."
-                )
-                break
-
-            # Stop if error is zero.
-            if estimator_error == 0 and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (training error is 0)."
-                )
+            if self._check_early_stop(iboost, sample_weight, estimator_error, early_termination_):
                 break
 
             sample_weight_sum = np.sum(sample_weight)
 
-            # Stop if the sum of sample weights has become non-positive.
-            if (
-                sample_weight is not None
-                and sample_weight_sum <= 0
-                and early_termination_
-            ):
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight_sum <= 0)."
-                )
-                break
-
             if iboost < self.n_estimators - 1:
-                # Normalize.
                 sample_weight /= sample_weight_sum
 
         return self
@@ -642,8 +648,6 @@ def __init__(
         random_state=None,
     ):
 
-        self.early_termination = early_termination
-
         super(ReweightBoostClassifier, self).__init__(
             estimator=estimator,
             n_estimators=n_estimators,
@@ -652,6 +656,8 @@ def __init__(
             random_state=random_state,
         )
 
+        self._train_state_ = _TrainingState(early_termination=early_termination)
+
     def _compute_mult_in_exp_weights_array(self, y_true, y_pred):
         """
         Compute the additional weights that need to be multiplied
@@ -683,7 +689,10 @@ def _preprocess_sample_weight(self, sample_weight, y):
 
     def _set_cost_matrix(self, how: str = "inverse"):
         """Set the cost matrix according to the 'how' parameter."""
-        classes, origin_distr = self._encode_map.values(), self.origin_distr_
+        classes, origin_distr = (
+            self._train_state_._encode_map.values(),
+            self._train_state_.origin_distr_,
+        )
         cost_matrix = []
         for c_pred in classes:
             cost_c = [
@@ -726,7 +735,9 @@ def _boost_real(self, iboost, X, y, sample_weight, random_state):
         y_predict_proba = estimator.predict_proba(X)
 
         y_predict = self.classes_.take(np.argmax(y_predict_proba, axis=1), axis=0)
-        y_predict = np.array(list(map(lambda x: self._encode_map[x], y_predict)))
+        y_predict = np.array(
+            list(map(lambda x: self._train_state_._encode_map[x], y_predict))
+        )
 
         # Instances incorrectly classified
         incorrect = y_predict != y
@@ -766,11 +777,11 @@ def _boost_real(self, iboost, X, y, sample_weight, random_state):
 
         # Compute additional weights for multiplication
         mult_in_exp_weight = self._compute_mult_in_exp_weights_array(
-            y_true=self._y_encoded, y_pred=y_predict
+            y_true=self._train_state_._y_encoded, y_pred=y_predict
         )
         mult_in_exp_weight /= mult_in_exp_weight.max()
         mult_out_exp_weight = self._compute_mult_out_exp_weights_array(
-            y_true=self._y_encoded, y_pred=y_predict
+            y_true=self._train_state_._y_encoded, y_pred=y_predict
         )
 
         # Only boost the weights if it will fit again
@@ -806,7 +817,9 @@ def _boost_discrete(self, iboost, X, y, sample_weight, random_state):
         # Instances incorrectly classified
         incorrect = y_predict != y
 
-        y_predict = np.array(list(map(lambda x: self._encode_map[x], y_predict)))
+        y_predict = np.array(
+            list(map(lambda x: self._train_state_._encode_map[x], y_predict))
+        )
 
         # Error fraction
         estimator_error = np.mean(np.average(incorrect, weights=sample_weight, axis=0))
@@ -833,10 +846,10 @@ def _boost_discrete(self, iboost, X, y, sample_weight, random_state):
 
         # Compute additional weights for multiplication
         mult_in_exp_weight = self._compute_mult_in_exp_weights_array(
-            y_true=self._y_encoded, y_pred=y_predict
+            y_true=self._train_state_._y_encoded, y_pred=y_predict
         )
         mult_out_exp_weight = self._compute_mult_out_exp_weights_array(
-            y_true=self._y_encoded, y_pred=y_predict
+            y_true=self._train_state_._y_encoded, y_pred=y_predict
         )
 
         # Only boost the weights if I will fit again
@@ -854,41 +867,22 @@ def _boost_discrete(self, iboost, X, y, sample_weight, random_state):
 
         return sample_weight, estimator_weight, estimator_error
 
-    @_deprecate_positional_args
-    def _fit(
-        self,
-        X,
-        y,
-        *,
-        sample_weight,
-        cost_matrix,
-        eval_datasets: dict,
-        eval_metrics: dict,
-        train_verbose: bool or int or dict,
-    ):
-
-        early_termination_ = check_type(
-            self.early_termination, "early_termination", bool
-        )
-
-        # Check that algorithm is supported.
+    def _validate_input_params(self):
         if self.algorithm not in ("SAMME", "SAMME.R"):
             raise ValueError("algorithm %s is not supported" % self.algorithm)
-
-        # Check parameters.
         if self.learning_rate <= 0:
             raise ValueError("learning_rate must be greater than zero")
 
-        if self.estimator == None or isinstance(
+    def _get_data_type_config(self):
+        if self.estimator is None or isinstance(
             self.estimator, (BaseDecisionTree, BaseForest)
         ):
-            DTYPE = np.float64  # from fast_dict.pxd
-            dtype = DTYPE
-            accept_sparse = "csc"
-        else:
-            dtype = None
-            accept_sparse = ["csr", "csc"]
+            return np.float64, "csc"
+        return None, ["csr", "csc"]
 
+    def _validate_and_prepare_data(
+        self, X, y, dtype, accept_sparse, eval_datasets, eval_metrics, train_verbose
+    ):
         check_x_y_args = {
             "accept_sparse": accept_sparse,
             "ensure_2d": True,
@@ -898,106 +892,146 @@ def _fit(
         }
         X, y = validate_data(self, X, y, **check_x_y_args)
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        self._train_state_.eval_datasets_ = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
 
         raw_y = y.copy()
-        self.classes_, self._y_encoded = np.unique(y, return_inverse=True)
-        self._encode_map = {
+        self.classes_, self._train_state_._y_encoded = np.unique(
+            y, return_inverse=True
+        )
+        self._train_state_._encode_map = {
             c: np.where(self.classes_ == c)[0][0] for c in self.classes_
         }
         self.n_classes_ = len(self.classes_)
 
-        # Store original class distribution
-        self.origin_distr_ = dict(Counter(self._y_encoded))
-        self.target_distr_ = dict(Counter(self._y_encoded))
+        self._train_state_.origin_distr_ = dict(
+            Counter(self._train_state_._y_encoded)
+        )
+        self._train_state_.target_distr_ = dict(
+            Counter(self._train_state_._y_encoded)
+        )
 
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        self._train_state_.eval_metrics_ = check_eval_metrics(eval_metrics)
 
-        self.train_verbose_ = check_train_verbose(
+        self._train_state_.train_verbose_ = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
 
         self._init_training_log_format()
 
-        # Check sample weight
+        return X, y, raw_y
+
+    def _init_sample_weight(self, sample_weight, X, y_encoded):
         sample_weight = _check_sample_weight(sample_weight, X, np.float64)
-        sample_weight = self._preprocess_sample_weight(sample_weight, self._y_encoded)
+        sample_weight = self._preprocess_sample_weight(sample_weight, y_encoded)
         sample_weight /= sample_weight.sum()
         if np.any(sample_weight < 0):
             raise ValueError("sample_weight cannot contain negative weights")
 
-        self.raw_sample_weight_ = sample_weight
+        self._train_state_.raw_sample_weight_ = sample_weight
 
-        sample_weight = copy(self.raw_sample_weight_)
+        return copy(sample_weight)
 
-        # Initialize & validate cost matrix
+    def _init_cost_matrix(self, cost_matrix):
         if cost_matrix is None:
             cost_matrix = self._set_cost_matrix()
         elif isinstance(cost_matrix, str):
             cost_matrix = self._set_cost_matrix(how=cost_matrix)
         cost_matrix = self._validate_cost_matrix(cost_matrix, self.n_classes_)
-        self.cost_matrix_ = cost_matrix
+        self._train_state_.cost_matrix_ = cost_matrix
 
+    def _init_ensemble(self):
         self._validate_estimator()
 
-        # Check random state
         random_state = check_random_state(self.random_state)
 
-        # Clear any previous fit results.
         self.estimators_ = []
         self.estimator_weights_ = np.zeros(self.n_estimators, dtype=np.float64)
         self.estimator_errors_ = np.ones(self.n_estimators, dtype=np.float64)
-        self.estimators_n_training_samples_ = np.zeros(self.n_estimators, dtype=int)
+        self._train_state_.estimators_n_training_samples_ = np.zeros(
+            self.n_estimators, dtype=int
+        )
 
-        # Genrate random seeds array
         seeds = random_state.randint(MAX_INT, size=self.n_estimators)
-        self._seeds = seeds
+        self._train_state_._seeds = seeds
+
+        return random_state
+
+    def _check_early_stop(self, iboost, sample_weight, estimator_error, early_termination_):
+        if sample_weight is None and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight is None)."
+            )
+            return True
+
+        if estimator_error == 0 and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (training error is 0)."
+            )
+            return True
+
+        if early_termination_ and np.sum(sample_weight) <= 0:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight_sum <= 0)."
+            )
+            return True
+
+        return False
+
+    @_deprecate_positional_args
+    def _fit(
+        self,
+        X,
+        y,
+        *,
+        sample_weight,
+        cost_matrix,
+        eval_datasets: dict,
+        eval_metrics: dict,
+        train_verbose: bool or int or dict,
+    ):
+
+        early_termination_ = check_type(
+            self._train_state_.early_termination, "early_termination", bool
+        )
+
+        self._validate_input_params()
+        dtype, accept_sparse = self._get_data_type_config()
+        X, y, raw_y = self._validate_and_prepare_data(
+            X, y, dtype, accept_sparse, eval_datasets, eval_metrics, train_verbose
+        )
+        sample_weight = self._init_sample_weight(
+            sample_weight, X, self._train_state_._y_encoded
+        )
+        self._init_cost_matrix(cost_matrix)
+        random_state = self._init_ensemble()
 
         for iboost in range(self.n_estimators):
-            # Boosting step
             sample_weight, estimator_weight, estimator_error = self._boost(
                 iboost, X, raw_y, sample_weight, random_state
             )
 
             self.estimator_weights_[iboost] = estimator_weight
             self.estimator_errors_[iboost] = estimator_error
-            self.estimators_n_training_samples_[iboost] = y.shape[0]
+            self._train_state_.estimators_n_training_samples_[iboost] = y.shape[0]
 
-            # Print training infomation to console.
             self._training_log_to_console(iboost, y)
 
-            # Early termination.
-            if sample_weight is None and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight is None)."
-                )
-                break
-
-            # Stop if error is zero.
-            if estimator_error == 0 and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (training error is 0)."
-                )
+            if self._check_early_stop(
+                iboost, sample_weight, estimator_error, early_termination_
+            ):
                 break
 
             sample_weight_sum = np.sum(sample_weight)
 
-            # Stop if the sum of sample weights has become non-positive.
-            if sample_weight_sum <= 0 and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight_sum <= 0)."
-                )
-                break
-
             if iboost < self.n_estimators - 1:
-                # Normalize.
                 sample_weight /= sample_weight_sum
 
         return self
diff --git a/imbens/ensemble/_compatible/adaboost_compatible.py b/imbens/ensemble/_compatible/adaboost_compatible.py
index 9c56148..46e64fe 100644
--- a/imbens/ensemble/_compatible/adaboost_compatible.py
+++ b/imbens/ensemble/_compatible/adaboost_compatible.py
@@ -22,12 +22,12 @@
         check_train_verbose,
         check_type,
     )
-    from ..base import MAX_INT, ImbalancedEnsembleClassifierMixin
+    from ..base import _TrainingState, MAX_INT, ImbalancedEnsembleClassifierMixin
 else:  # pragma: no cover
     import sys  # For local test
 
     sys.path.append("../..")
-    from ensemble.base import ImbalancedEnsembleClassifierMixin, MAX_INT
+    from ensemble.base import _TrainingState, ImbalancedEnsembleClassifierMixin, MAX_INT
     from utils._validation_data import check_eval_datasets
     from utils._validation_param import check_train_verbose, check_eval_metrics
     from utils._validation import _deprecate_positional_args
@@ -158,6 +158,9 @@ class CompatibleAdaBoostClassifier(
     {example}
     """
 
+    __name__ = _method_name
+    _properties = _properties
+
     def __init__(
         self,
         estimator=None,
@@ -178,8 +181,7 @@ def __init__(
             random_state=random_state,
         )
 
-        self.__name__ = _method_name
-        self._properties = _properties
+        self._train_state_ = _TrainingState()
 
     @_deprecate_positional_args
     @FuncSubstitution(
@@ -223,23 +225,47 @@ def fit(
         self : object
         """
 
+        early_termination_, check_x_y_args = self._validate_parameters()
+        X, y = validate_data(self, X, y, **check_x_y_args)
+
+        self._init_training_state(
+            X, y, sample_weight, eval_datasets, eval_metrics, train_verbose, check_x_y_args
+        )
+        sample_weight = copy(self._train_state_.raw_sample_weight_)
+
+        for iboost in range(self.n_estimators):
+            sample_weight, estimator_weight, estimator_error = self._boost(
+                iboost, X, y, sample_weight, self._train_state_._random_state
+            )
+
+            self.estimator_weights_[iboost] = estimator_weight
+            self.estimator_errors_[iboost] = estimator_error
+            self._train_state_.estimators_n_training_samples_[iboost] = y.shape[0]
+
+            self._training_log_to_console(iboost, y)
+
+            if self._check_early_termination(
+                iboost, sample_weight, estimator_error, early_termination_
+            ):
+                break
+
+        return self
+
+    def _validate_parameters(self):
         early_termination_ = check_type(
             self.early_termination, "early_termination", bool
         )
 
-        # Check that algorithm is supported.
         if self.algorithm not in ("SAMME", "SAMME.R"):
             raise ValueError("algorithm %s is not supported" % self.algorithm)
 
-        # Check parameters.
         if self.learning_rate <= 0:
             raise ValueError("learning_rate must be greater than zero")
 
         if self.estimator is None or isinstance(
             self.estimator, (BaseDecisionTree, BaseForest)
         ):
-            DTYPE = np.float64  # from fast_dict.pxd
-            dtype = DTYPE
+            dtype = np.float64
             accept_sparse = "csc"
         else:
             dtype = None
@@ -252,97 +278,82 @@ def fit(
             "dtype": dtype,
             "y_numeric": False,
         }
-        X, y = validate_data(self, X, y, **check_x_y_args)
+        return early_termination_, check_x_y_args
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+    def _init_training_state(
+        self, X, y, sample_weight, eval_datasets, eval_metrics, train_verbose, check_x_y_args
+    ):
+        self._train_state_.eval_datasets_ = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
 
         self.classes_, _ = np.unique(y, return_inverse=True)
         self.n_classes_ = len(self.classes_)
 
-        self.origin_distr_ = dict(Counter(y))
-        self.target_distr_ = dict(Counter(y))
+        self._train_state_.origin_distr_ = dict(Counter(y))
+        self._train_state_.target_distr_ = dict(Counter(y))
 
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        self._train_state_.eval_metrics_ = check_eval_metrics(eval_metrics)
 
-        self.train_verbose_ = check_train_verbose(
+        self._train_state_.train_verbose_ = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
 
         self._init_training_log_format()
 
-        # Check sample weight
         sample_weight = _check_sample_weight(sample_weight, X, np.float64)
         sample_weight /= sample_weight.sum()
         if np.any(sample_weight < 0):
             raise ValueError("sample_weight cannot contain negative weights")
 
-        self.raw_sample_weight_ = sample_weight
-
-        sample_weight = copy(self.raw_sample_weight_)
+        self._train_state_.raw_sample_weight_ = sample_weight
 
         self._validate_estimator()
 
-        # Check random state
         random_state = check_random_state(self.random_state)
 
-        # Clear any previous fit results.
         self.estimators_ = []
         self.estimator_weights_ = np.zeros(self.n_estimators, dtype=np.float64)
         self.estimator_errors_ = np.ones(self.n_estimators, dtype=np.float64)
-        self.estimators_n_training_samples_ = np.zeros(self.n_estimators, dtype=int)
+        self._train_state_.estimators_n_training_samples_ = np.zeros(
+            self.n_estimators, dtype=int
+        )
 
-        # Genrate random seeds array
         seeds = random_state.randint(MAX_INT, size=self.n_estimators)
-        self._seeds = seeds
-
-        for iboost in range(self.n_estimators):
-            # Boosting step
-            sample_weight, estimator_weight, estimator_error = self._boost(
-                iboost, X, y, sample_weight, random_state
+        self._train_state_._seeds = seeds
+        self._train_state_._random_state = random_state
+
+    def _check_early_termination(self, iboost, sample_weight, estimator_error, early_termination_):
+        if sample_weight is None and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight is None)."
             )
+            return True
 
-            self.estimator_weights_[iboost] = estimator_weight
-            self.estimator_errors_[iboost] = estimator_error
-            self.estimators_n_training_samples_[iboost] = y.shape[0]
-
-            # Print training infomation to console.
-            self._training_log_to_console(iboost, y)
-
-            # Early termination.
-            if sample_weight is None and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight is None)."
-                )
-                break
-
-            # Stop if error is zero.
-            if estimator_error == 0 and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (training error is 0)."
-                )
-                break
+        if estimator_error == 0 and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (training error is 0)."
+            )
+            return True
 
-            sample_weight_sum = np.sum(sample_weight)
+        sample_weight_sum = np.sum(sample_weight)
 
-            # Stop if the sum of sample weights has become non-positive.
-            if sample_weight_sum <= 0 and early_termination_:
-                print(
-                    f"Training early-stop at iteration"
-                    f" {iboost+1}/{self.n_estimators}"
-                    f" (sample_weight_sum <= 0)."
-                )
-                break
+        if sample_weight_sum <= 0 and early_termination_:
+            print(
+                f"Training early-stop at iteration"
+                f" {iboost+1}/{self.n_estimators}"
+                f" (sample_weight_sum <= 0)."
+            )
+            return True
 
-            if iboost < self.n_estimators - 1:
-                # Normalize.
-                sample_weight /= sample_weight_sum
+        if iboost < self.n_estimators - 1:
+            sample_weight /= sample_weight_sum
 
-        return self
+        return False
 
     @FuncGlossarySubstitution(_super.decision_function, "classes_")
     def decision_function(self, X):
diff --git a/imbens/ensemble/_compatible/bagging_compatible.py b/imbens/ensemble/_compatible/bagging_compatible.py
index d8ef509..f5038cf 100644
--- a/imbens/ensemble/_compatible/bagging_compatible.py
+++ b/imbens/ensemble/_compatible/bagging_compatible.py
@@ -19,13 +19,13 @@
     from ...utils._validation_data import check_eval_datasets
     from ...utils._validation_param import check_eval_metrics, check_train_verbose
     from .._bagging import _parallel_build_estimators
-    from ..base import MAX_INT, ImbalancedEnsembleClassifierMixin
+    from ..base import _TrainingState, MAX_INT, ImbalancedEnsembleClassifierMixin
 else:  # pragma: no cover
     import sys  # For local test
 
     sys.path.append("../..")
-    from ensemble.base import ImbalancedEnsembleClassifierMixin, MAX_INT
     from ensemble._bagging import _parallel_build_estimators
+    from ensemble.base import _TrainingState, ImbalancedEnsembleClassifierMixin, MAX_INT
     from utils._validation_data import check_eval_datasets
     from utils._validation_param import check_train_verbose, check_eval_metrics
     from utils._validation import _deprecate_positional_args, check_target_type
@@ -194,6 +194,37 @@ class CompatibleBaggingClassifier(ImbalancedEnsembleClassifierMixin, BaggingClas
     {example}
     """
 
+    __name__ = _method_name
+    _properties = _properties
+
+    @property
+    def _max_samples(self):
+        if hasattr(self, '_train_state_'):
+            return self._train_state_._max_samples
+        return self.max_samples
+
+    @property
+    def _max_features(self):
+        if hasattr(self, '_train_state_'):
+            return self._train_state_._max_features
+        return self.max_features
+
+    @property
+    def _n_samples(self):
+        if hasattr(self, '_train_state_'):
+            return self._train_state_._n_samples
+        return 0
+
+    @property
+    def _seeds(self):
+        if hasattr(self, '_train_state_'):
+            return self._train_state_._seeds
+        return None
+
+    @_seeds.setter
+    def _seeds(self, value):
+        self._train_state_._seeds = value
+
     @_deprecate_positional_args
     def __init__(
         self,
@@ -225,8 +256,7 @@ def __init__(
             verbose=verbose,
         )
 
-        self.__name__ = _method_name
-        self._properties = _properties
+        self._train_state_ = _TrainingState()
 
     @_deprecate_positional_args
     @FuncSubstitution(
@@ -279,7 +309,6 @@ def fit(
 
         random_state = check_random_state(self.random_state)
 
-        # Convert data (X is required to be 2d and indexable)
         check_x_y_args = {
             "accept_sparse": ["csr", "csc"],
             "dtype": None,
@@ -288,30 +317,93 @@ def fit(
         }
         X, y = validate_data(self, X, y, **check_x_y_args)
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
-
-        # Check evaluation metrics
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
-
-        # Check verbose
-        self.train_verbose_ = check_train_verbose(
-            train_verbose, self.n_estimators, **self._properties
+        self._pre_fit_setup(
+            eval_datasets, eval_metrics, train_verbose, check_x_y_args, X, y
         )
-        self._init_training_log_format()
 
         if sample_weight is not None:
             sample_weight = _check_sample_weight(sample_weight, X, dtype=None)
 
-        # Remap output
         n_samples, self.n_features_in_ = X.shape
-        self._n_samples = n_samples
+        self._train_state_._n_samples = n_samples
         y = self._validate_y(y)
 
-        # Check parameters
         self._validate_estimator(self._get_estimator())
 
-        # Validate max_samples
+        max_samples = self._validate_max_samples(max_samples, X)
+        max_features = self._validate_max_features()
+
+        self._validate_common_checks()
+
+        n_more_estimators = self._initialize_warm_start()
+        if n_more_estimators == 0:
+            return self
+
+        self._parallel_fit(n_more_estimators, X, y, sample_weight, random_state)
+
+        self._post_fit_process(X, y)
+
+        return self
+
+    def _pre_fit_setup(
+        self, eval_datasets, eval_metrics, train_verbose, check_x_y_args, X, y
+    ):
+        self._train_state_.eval_datasets_ = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
+        self._train_state_.eval_metrics_ = check_eval_metrics(eval_metrics)
+        self._train_state_.train_verbose_ = check_train_verbose(
+            train_verbose, self.n_estimators, **self._properties
+        )
+        self._init_training_log_format()
+
+    def _parallel_fit(
+        self, n_more_estimators, X, y, sample_weight, random_state
+    ):
+        n_jobs, n_estimators, starts = _partition_estimators(
+            n_more_estimators, self.n_jobs
+        )
+        total_n_estimators = sum(n_estimators)
+
+        if self.warm_start and len(self.estimators_) > 0:
+            random_state.randint(MAX_INT, size=len(self.estimators_))
+
+        seeds = random_state.randint(MAX_INT, size=n_more_estimators)
+        self._train_state_._seeds = seeds
+
+        all_results = Parallel(
+            n_jobs=n_jobs, verbose=self.verbose, **self._parallel_args()
+        )(
+            delayed(_parallel_build_estimators)(
+                n_estimators[i],
+                self,
+                X,
+                y,
+                sample_weight,
+                seeds[starts[i] : starts[i + 1]],
+                total_n_estimators,
+                verbose=self.verbose,
+            )
+            for i in range(n_jobs)
+        )
+
+        self.estimators_ += list(
+            itertools.chain.from_iterable(t[0] for t in all_results)
+        )
+        self.estimators_features_ += list(
+            itertools.chain.from_iterable(t[1] for t in all_results)
+        )
+        self.estimators_n_training_samples_ += list(
+            itertools.chain.from_iterable(t[2] for t in all_results)
+        )
+
+    def _post_fit_process(self, X, y):
+        if self.oob_score:
+            self._set_oob_score(X, y)
+
+        self._training_log_to_console()
+
+    def _validate_max_samples(self, max_samples, X):
         if max_samples is None:
             max_samples = self.max_samples
         if not isinstance(max_samples, numbers.Integral):
@@ -320,10 +412,10 @@ def fit(
         if not (0 < max_samples <= X.shape[0]):
             raise ValueError("max_samples must be in (0, n_samples]")
 
-        # Store validated integer row sampling value
-        self._max_samples = max_samples
+        self._train_state_._max_samples = max_samples
+        return max_samples
 
-        # Validate max_features
+    def _validate_max_features(self):
         if isinstance(self.max_features, numbers.Integral):
             max_features = self.max_features
         elif isinstance(self.max_features, float):
@@ -336,10 +428,10 @@ def fit(
 
         max_features = max(1, int(max_features))
 
-        # Store validated integer feature sampling value
-        self._max_features = max_features
+        self._train_state_._max_features = max_features
+        return max_features
 
-        # Other checks
+    def _validate_common_checks(self):
         if not self.bootstrap and self.oob_score:
             raise ValueError(
                 "Out of bag estimation only available" " if bootstrap=True"
@@ -353,8 +445,8 @@ def fit(
         if hasattr(self, "oob_score_") and self.warm_start:
             del self.oob_score_
 
+    def _initialize_warm_start(self):
         if not self.warm_start or not hasattr(self, "estimators_"):
-            # Free allocated memory, if any
             self.estimators_ = []
             self.estimators_features_ = []
             self.estimators_n_training_samples_ = []
@@ -373,56 +465,8 @@ def fit(
                 "Warm-start fitting without increasing n_estimators does not "
                 "fit new trees."
             )
-            return self
-
-        # Parallel loop
-        n_jobs, n_estimators, starts = _partition_estimators(
-            n_more_estimators, self.n_jobs
-        )
-        total_n_estimators = sum(n_estimators)
-
-        # Advance random state to state after training
-        # the first n_estimators
-        if self.warm_start and len(self.estimators_) > 0:
-            random_state.randint(MAX_INT, size=len(self.estimators_))
-
-        seeds = random_state.randint(MAX_INT, size=n_more_estimators)
-        self._seeds = seeds
-
-        all_results = Parallel(
-            n_jobs=n_jobs, verbose=self.verbose, **self._parallel_args()
-        )(
-            delayed(_parallel_build_estimators)(
-                n_estimators[i],
-                self,
-                X,
-                y,
-                sample_weight,
-                seeds[starts[i] : starts[i + 1]],
-                total_n_estimators,
-                verbose=self.verbose,
-            )
-            for i in range(n_jobs)
-        )
-
-        # Reduce
-        self.estimators_ += list(
-            itertools.chain.from_iterable(t[0] for t in all_results)
-        )
-        self.estimators_features_ += list(
-            itertools.chain.from_iterable(t[1] for t in all_results)
-        )
-        self.estimators_n_training_samples_ += list(
-            itertools.chain.from_iterable(t[2] for t in all_results)
-        )
 
-        if self.oob_score:
-            self._set_oob_score(X, y)
-
-        # Print training infomation to console.
-        self._training_log_to_console()
-
-        return self
+        return n_more_estimators
 
     @FuncGlossarySubstitution(_super.predict_proba, "classes_")
     def predict_proba(self, X):
diff --git a/imbens/ensemble/_over_sampling/kmeans_smote_boost.py b/imbens/ensemble/_over_sampling/kmeans_smote_boost.py
index 9d689ef..46da07b 100644
--- a/imbens/ensemble/_over_sampling/kmeans_smote_boost.py
+++ b/imbens/ensemble/_over_sampling/kmeans_smote_boost.py
@@ -211,10 +211,12 @@ def __init__(
 
         self.k_neighbors = k_neighbors
         self.k_neighbors_ = check_type(k_neighbors, "k_neighbors", numbers.Integral)
-        self.n_jobs_sampler = n_jobs_sampler
-        self.kmeans_estimator = kmeans_estimator
-        self.cluster_balance_threshold = cluster_balance_threshold
-        self.density_exponent = density_exponent
+        self.kmeans_smote_kwargs = {
+            "n_jobs": n_jobs_sampler,
+            "kmeans_estimator": kmeans_estimator,
+            "cluster_balance_threshold": cluster_balance_threshold,
+            "density_exponent": density_exponent,
+        }
         if self.k_neighbors_ < 1:
             raise ValueError(
                 f"The 'k_neighbors' parameter of NearestNeighbors must be"
@@ -277,10 +279,7 @@ def fit(
 
         kmeans_smote_sampler_kwargs = {
             "k_neighbors": self.k_neighbors_,
-            "n_jobs": self.n_jobs_sampler,
-            "kmeans_estimator": self.kmeans_estimator,
-            "cluster_balance_threshold": self.cluster_balance_threshold,
-            "density_exponent": self.density_exponent,
+            **self.kmeans_smote_kwargs,
         }
         update_x_y_after_resample = True
 
diff --git a/imbens/ensemble/_reweighting/adauboost.py b/imbens/ensemble/_reweighting/adauboost.py
index e11db41..3d61f08 100644
--- a/imbens/ensemble/_reweighting/adauboost.py
+++ b/imbens/ensemble/_reweighting/adauboost.py
@@ -234,7 +234,10 @@ def _check_beta(self, beta: str or dict) -> dict:
     def _set_beta(self, how: str = "inverse") -> dict:
         """Set the self.beta_ by 'how'."""
 
-        classes, origin_distr = self._encode_map.values(), self.origin_distr_
+        classes, origin_distr = (
+            self._train_state_._encode_map.values(),
+            self._train_state_.origin_distr_,
+        )
         c_maj = max(origin_distr.keys(), key=(lambda x: origin_distr[x]))
         beta = [origin_distr[c_maj] / origin_distr[c_min] for c_min in classes]
         if how == "uniform":
diff --git a/imbens/ensemble/_under_sampling/balance_cascade.py b/imbens/ensemble/_under_sampling/balance_cascade.py
index 1901a64..61cb3c2 100644
--- a/imbens/ensemble/_under_sampling/balance_cascade.py
+++ b/imbens/ensemble/_under_sampling/balance_cascade.py
@@ -24,12 +24,12 @@
         check_target_label_and_n_target_samples,
         check_train_verbose,
     )
-    from ..base import MAX_INT, BaseImbalancedEnsemble
+    from ..base import MAX_INT, BaseImbalancedEnsemble, _TrainingState
 else:  # pragma: no cover
     import sys  # For local test
 
     sys.path.append("../..")
-    from ensemble.base import BaseImbalancedEnsemble, MAX_INT
+    from ensemble.base import BaseImbalancedEnsemble, MAX_INT, _TrainingState
     from sampler._under_sampling import BalanceCascadeUnderSampler
     from utils._validation_data import check_eval_datasets
     from utils._validation_param import (
@@ -50,6 +50,15 @@
 
 import numpy as np
 
+# Attributes stored in _TrainingState for backward-compatible delegation
+_TRAIN_STATE_ATTRS = frozenset({
+    'eval_datasets_', 'origin_distr_', 'target_label_', 'target_distr_',
+    'balancing_schedule_', 'eval_metrics_', 'train_verbose_',
+    'train_verbose_format_', 'keep_ratios_', '_seeds',
+    'y_pred_proba_latest', 'estimators_n_training_samples_',
+    'sample_weights_',
+})
+
 # Properties
 _method_name = 'BalanceCascadeClassifier'
 _sampler_class = BalanceCascadeUnderSampler
@@ -186,6 +195,13 @@ def __init__(
 
         self.replacement = replacement
 
+    def __getattr__(self, name):
+        if '_train_state_' in self.__dict__ and name in _TRAIN_STATE_ATTRS:
+            return getattr(self._train_state_, name)
+        raise AttributeError(
+            f"'{type(self).__name__}' object has no attribute '{name}'"
+        )
+
     @_deprecate_positional_args
     @FuncSubstitution(
         target_label=_get_parameter_docstring('target_label', **_properties),
@@ -254,16 +270,20 @@ def _fit(
             self.replacement,
         )
 
+        # Initialize training state container
+        self._train_state_ = _TrainingState()
+
         # Check evaluation data
         check_x_y_args = self.check_x_y_args
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        self._train_state_.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
 
         # Check target sample strategy
         origin_distr_ = dict(Counter(y))
         target_label_, target_distr_ = check_target_label_and_n_target_samples(
             y, target_label, n_target_samples, self._sampling_type
         )
-        self.origin_distr_, self.target_label_, self.target_distr_ = (
+        ts = self._train_state_
+        ts.origin_distr_, ts.target_label_, ts.target_distr_ = (
             origin_distr_,
             target_label_,
             target_distr_,
@@ -271,13 +291,13 @@ def _fit(
 
         # Check balancing schedule
         balancing_schedule_ = check_balancing_schedule(balancing_schedule)
-        self.balancing_schedule_ = balancing_schedule_
+        ts.balancing_schedule_ = balancing_schedule_
 
         # Check evaluation metrics
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        ts.eval_metrics_ = check_eval_metrics(eval_metrics)
 
         # Check training train_verbose format
-        self.train_verbose_ = check_train_verbose(
+        ts.train_verbose_ = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
 
@@ -287,13 +307,13 @@ def _fit(
         # Clear any previous fit results.
         self.estimators_ = []
         self.estimators_features_ = []
-        self.estimators_n_training_samples_ = np.zeros(n_estimators, dtype=int)
+        ts.estimators_n_training_samples_ = np.zeros(n_estimators, dtype=int)
         self.samplers_ = []
-        self.sample_weights_ = []
+        ts.sample_weights_ = []
 
         # Genrate random seeds array
         seeds = random_state.randint(MAX_INT, size=n_estimators)
-        self._seeds = seeds
+        ts._seeds = seeds
 
         # Initialize the keep_ratios and dropped_index
         keep_ratios = {
@@ -302,7 +322,7 @@ def _fit(
             )
             for label in classes_
         }
-        self.keep_ratios_ = keep_ratios
+        ts.keep_ratios_ = keep_ratios
         dropped_index = np.full_like(y, fill_value=False, dtype=bool)
 
         # Check if sample_weight is specified
@@ -348,11 +368,11 @@ def _fit(
             resample_out = sampler.fit_resample(
                 X,
                 y,
-                y_pred_proba=self.y_pred_proba_latest,
+                y_pred_proba=ts.y_pred_proba_latest,
                 dropped_index=dropped_index,
                 keep_populations=keep_populations,
                 classes_=classes_,
-                encode_map=self._encode_map,
+                encode_map=self._internal_state_['encode_map'],
                 sample_weight=sample_weight,
             )
 
@@ -373,8 +393,8 @@ def _fit(
                 (X_resampled, y_resampled, dropped_index) = resample_out
                 estimator.fit(X_resampled, y_resampled)
 
-            self.estimators_features_.append(self.features_)
-            self.estimators_n_training_samples_[i_iter] = y_resampled.shape[0]
+            self.estimators_features_.append(self._internal_state_['features'])
+            ts.estimators_n_training_samples_[i_iter] = y_resampled.shape[0]
 
             # Print training infomation to console.
             self._training_log_to_console(i_iter, y_resampled)
@@ -386,14 +406,15 @@ def _update_cached_prediction_probabilities(self, i_iter, X):
         data during ensemble training. Must be called in each iteration before fit the
         estimator."""
 
+        ts = self._train_state_
         if i_iter == 0:
-            self.y_pred_proba_latest = np.zeros(
-                (self._n_samples, self.n_classes_), dtype=np.float64
+            ts.y_pred_proba_latest = np.zeros(
+                (self._internal_state_['n_samples'], self.n_classes_), dtype=np.float64
             )
         else:
-            y_pred_proba_latest = self.y_pred_proba_latest
+            y_pred_proba_latest = ts.y_pred_proba_latest
             y_pred_proba_new = self.estimators_[-1].predict_proba(X)
-            self.y_pred_proba_latest = (
+            ts.y_pred_proba_latest = (
                 y_pred_proba_latest * i_iter + y_pred_proba_new
             ) / (i_iter + 1)
         return
diff --git a/imbens/ensemble/_under_sampling/balanced_random_forest.py b/imbens/ensemble/_under_sampling/balanced_random_forest.py
index a67b87f..3eafd59 100644
--- a/imbens/ensemble/_under_sampling/balanced_random_forest.py
+++ b/imbens/ensemble/_under_sampling/balanced_random_forest.py
@@ -43,6 +43,7 @@
 import numbers
 from collections import Counter
 from copy import deepcopy
+from types import SimpleNamespace
 from warnings import warn
 
 import numpy as np
@@ -73,7 +74,7 @@
 _ensemble_type = "random-forest"
 _training_type = "parallel"
 
-_properties = {
+_CLASS_PROPERTIES = {
     "sampling_type": _sampling_type,
     "solution_type": _solution_type,
     "ensemble_type": _ensemble_type,
@@ -119,11 +120,11 @@ def _local_parallel_build_trees(
 
 
 @Substitution(
-    random_state=_get_parameter_docstring("random_state", **_properties),
-    n_jobs=_get_parameter_docstring("n_jobs", **_properties),
+    random_state=_get_parameter_docstring("random_state", **_CLASS_PROPERTIES),
+    n_jobs=_get_parameter_docstring("n_jobs", **_CLASS_PROPERTIES),
     example=_get_example_docstring(_method_name),
 )
-class BalancedRandomForestClassifier(
+class BalancedRandomForestClassifier(  # pylint: disable=too-many-instance-attributes
     ImbalancedEnsembleClassifierMixin, RandomForestClassifier
 ):
     """A balanced random forest classifier.
@@ -323,6 +324,59 @@ class labels (multi-output problem).
     {example}
     """
 
+    # --- Properties for backward compatibility ---
+    # (store data internally in _config)
+
+    @property
+    def __name__(self):
+        return self._config['name']
+
+    @property
+    def _sampling_type(self):
+        return self._config['sampling_type']
+
+    @property
+    def _sampler_class(self):
+        return _sampler_class
+
+    @property
+    def _properties(self):
+        return self._config['properties']
+
+    @property
+    def sampling_strategy(self):
+        return self._config['sampling_strategy']
+
+    @sampling_strategy.setter
+    def sampling_strategy(self, value):
+        self._config['sampling_strategy'] = value
+
+    @property
+    def replacement(self):
+        return self._config['replacement']
+
+    @replacement.setter
+    def replacement(self, value):
+        self._config['replacement'] = value
+
+    @property
+    def eval_datasets_(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.eval_datasets
+        return {}
+
+    @property
+    def eval_metrics_(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.eval_metrics
+        return {}
+
+    @property
+    def train_verbose_(self):
+        if hasattr(self, '_train_state'):
+            return self._train_state.train_verbose
+        return None
+
     @_deprecate_positional_args
     def __init__(
         self,
@@ -369,13 +423,13 @@ def __init__(
             max_samples=max_samples,
         )
 
-        self.__name__ = _method_name
-        self._sampling_type = _sampling_type
-        self._sampler_class = _sampler_class
-        self._properties = _properties
-
-        self.sampling_strategy = sampling_strategy
-        self.replacement = replacement
+        self._config = {
+            'name': _method_name,
+            'sampling_type': _sampling_type,
+            'properties': _CLASS_PROPERTIES,
+            'sampling_strategy': sampling_strategy,
+            'replacement': replacement,
+        }
 
     def _validate_estimator(self, default=DecisionTreeClassifier()):
         """Check the estimator and the n_estimator attribute, set the
@@ -396,7 +450,7 @@ def _validate_estimator(self, default=DecisionTreeClassifier()):
             self.estimator_ = clone(default)
 
         self.sampler_ = RandomUnderSampler(
-            sampling_strategy=self._sampling_strategy,
+            sampling_strategy=self._train_state.sampling_strategy,
             replacement=self.replacement,
         )
 
@@ -419,7 +473,7 @@ def _make_sampler_estimator(self, random_state=None):
     @FuncSubstitution(
         eval_datasets=_get_parameter_docstring("eval_datasets"),
         eval_metrics=_get_parameter_docstring("eval_metrics"),
-        train_verbose=_get_parameter_docstring("train_verbose", **_properties),
+        train_verbose=_get_parameter_docstring("train_verbose", **_CLASS_PROPERTIES),
     )
     def fit(
         self,
@@ -463,7 +517,30 @@ def fit(
             The fitted instance.
         """
 
-        # Validate or convert input data
+        X, y_encoded, sample_weight, n_samples_bootstrap = self._validate_and_prepare_data(
+            X, y, sample_weight, eval_datasets, eval_metrics, train_verbose
+        )
+
+        random_state = check_random_state(self.random_state)
+        n_more_estimators = self._get_n_more_estimators(random_state)
+
+        if n_more_estimators > 0:
+            self._build_estimators(
+                X, y_encoded, sample_weight, n_samples_bootstrap,
+                n_more_estimators, random_state
+            )
+
+        self._compute_oob_score_if_needed(X, y_encoded)
+
+        if hasattr(self, "classes_") and self.n_outputs_ == 1:
+            self.n_classes_ = self.n_classes_[0]
+            self.classes_ = self.classes_[0]
+
+        self._training_log_to_console()
+
+        return self
+
+    def _validate_input_data(self, X, y):
         if issparse(y):
             raise ValueError("sparse multilabel-indicator for y is not supported.")
 
@@ -474,29 +551,26 @@ def fit(
         }
         X, y = validate_data(self, X, y, **check_x_y_args)
 
-        # Check evaluation data
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        if issparse(X):
+            X.sort_indices()
+
+        _, self.n_features_in_ = X.shape
+
+        return X, y, check_x_y_args
 
-        # Check evaluation metrics
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+    def _init_training_state(self, eval_datasets, eval_metrics, train_verbose, X, y, check_x_y_args):
+        self._train_state = SimpleNamespace()
 
-        # Check verbose
-        self.train_verbose_ = check_train_verbose(
+        self._train_state.eval_datasets = check_eval_datasets(
+            eval_datasets, X, y, **check_x_y_args
+        )
+        self._train_state.eval_metrics = check_eval_metrics(eval_metrics)
+        self._train_state.train_verbose = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
         self._init_training_log_format()
 
-        if sample_weight is not None:
-            sample_weight = _check_sample_weight(sample_weight, X)
-
-        if issparse(X):
-            # Pre-sort indices to avoid that each individual tree of the
-            # ensemble sorts the indices.
-            X.sort_indices()
-
-        # Remap output
-        _, self.n_features_in_ = X.shape
-
+    def _prepare_target(self, y):
         y = np.atleast_1d(y)
         if y.ndim == 2 and y.shape[1] == 1:
             warn(
@@ -508,8 +582,6 @@ def fit(
             )
 
         if y.ndim == 1:
-            # reshape is necessary to preserve the data contiguity against vs
-            # [:, np.newaxis] that does not.
             y = np.reshape(y, (-1, 1))
 
         self.n_outputs_ = y.shape[1]
@@ -519,8 +591,11 @@ def fit(
         if getattr(y, "dtype", None) != DOUBLE or not y.flags.contiguous:
             y_encoded = np.ascontiguousarray(y_encoded, dtype=DOUBLE)
 
+        return y_encoded, expanded_class_weight
+
+    def _resolve_sampling_strategy(self, y):
         if isinstance(self.sampling_strategy, dict):
-            self._sampling_strategy = {
+            self._train_state.sampling_strategy = {
                 np.where(self.classes_[0] == key)[0][0]: value
                 for key, value in check_sampling_strategy(
                     self.sampling_strategy,
@@ -529,29 +604,46 @@ def fit(
                 ).items()
             }
         else:
-            self._sampling_strategy = self.sampling_strategy
+            self._train_state.sampling_strategy = self.sampling_strategy
 
+    def _apply_expanded_weights(self, sample_weight, expanded_class_weight):
         if expanded_class_weight is not None:
             if sample_weight is not None:
                 sample_weight = sample_weight * expanded_class_weight
             else:
                 sample_weight = expanded_class_weight
+        return sample_weight
+
+    def _validate_and_prepare_data(self, X, y, sample_weight,
+                                    eval_datasets, eval_metrics, train_verbose):
+        X, y, check_x_y_args = self._validate_input_data(X, y)
+
+        self._init_training_state(
+            eval_datasets, eval_metrics, train_verbose, X, y, check_x_y_args
+        )
+
+        if sample_weight is not None:
+            sample_weight = _check_sample_weight(sample_weight, X)
+
+        y_encoded, expanded_class_weight = self._prepare_target(y)
+
+        self._resolve_sampling_strategy(y)
+
+        sample_weight = self._apply_expanded_weights(sample_weight, expanded_class_weight)
 
-        # Get bootstrap sample size
         n_samples_bootstrap = _get_n_samples_bootstrap(
             n_samples=X.shape[0], max_samples=self.max_samples
         )
 
-        # Check parameters
         self._validate_estimator()
 
         if not self.bootstrap and self.oob_score:
             raise ValueError("Out of bag estimation only available if bootstrap=True")
 
-        random_state = check_random_state(self.random_state)
+        return X, y_encoded, sample_weight, n_samples_bootstrap
 
+    def _get_n_more_estimators(self, random_state):
         if not self.warm_start or not hasattr(self, "estimators_"):
-            # Free allocated memory, if any
             self.estimators_ = []
             self.estimators_n_training_samples_ = []
             self.samplers_ = []
@@ -573,64 +665,48 @@ def fit(
             )
         else:
             if self.warm_start and len(self.estimators_) > 0:
-                # We draw from the random state to get the random state we
-                # would have got if we hadn't used a warm_start.
                 random_state.randint(MAX_INT, size=len(self.estimators_))
 
-            trees = []
-            samplers = []
-            for _ in range(n_more_estimators):
-                tree, sampler = self._make_sampler_estimator(random_state=random_state)
-                trees.append(tree)
-                samplers.append(sampler)
-
-            # Parallel loop: we prefer the threading backend as the Cython code
-            # for fitting the trees is internally releasing the Python GIL
-            # making threading more efficient than multiprocessing in
-            # that case. However, we respect any parallel_backend contexts set
-            # at a higher level, since correctness does not rely on using
-            # threads.
-            samplers_trees = Parallel(
-                n_jobs=self.n_jobs, verbose=self.verbose, prefer="threads"
-            )(
-                delayed(_local_parallel_build_trees)(
-                    s,
-                    t,
-                    self.bootstrap,
-                    X,
-                    y_encoded,
-                    sample_weight,
-                    i,
-                    len(trees),
-                    verbose=self.verbose,
-                    class_weight=self.class_weight,
-                    n_samples_bootstrap=n_samples_bootstrap,
-                    forest=self,
-                )
-                for i, (s, t) in enumerate(zip(samplers, trees))
-            )
-            samplers, trees, n_training_samples = zip(*samplers_trees)
-
-            # Collect newly grown trees
-            self.estimators_.extend(trees)
-            self.samplers_.extend(samplers)
-            self.estimators_n_training_samples_.extend(n_training_samples)
-
-            # Create pipeline with the fitted samplers and trees
-            self.pipelines_.extend(
-                [
-                    make_pipeline(deepcopy(s), deepcopy(t))
-                    for s, t in zip(samplers, trees)
-                ]
+        return n_more_estimators
+
+    def _build_estimators(self, X, y_encoded, sample_weight,
+                          n_samples_bootstrap, n_more_estimators, random_state):
+        trees = []
+        samplers = []
+        for _ in range(n_more_estimators):
+            tree, sampler = self._make_sampler_estimator(random_state=random_state)
+            trees.append(tree)
+            samplers.append(sampler)
+
+        samplers_trees = Parallel(
+            n_jobs=self.n_jobs, verbose=self.verbose, prefer="threads"
+        )(
+            delayed(_local_parallel_build_trees)(
+                s, t, self.bootstrap, X, y_encoded, sample_weight,
+                i, len(trees), verbose=self.verbose,
+                class_weight=self.class_weight,
+                n_samples_bootstrap=n_samples_bootstrap,
+                forest=self,
             )
+            for i, (s, t) in enumerate(zip(samplers, trees))
+        )
+        samplers, trees, n_training_samples = zip(*samplers_trees)
+
+        self.estimators_.extend(trees)
+        self.samplers_.extend(samplers)
+        self.estimators_n_training_samples_.extend(n_training_samples)
+
+        self.pipelines_.extend(
+            [
+                make_pipeline(deepcopy(s), deepcopy(t))
+                for s, t in zip(samplers, trees)
+            ]
+        )
 
+    def _compute_oob_score_if_needed(self, X, y_encoded):
         if self.oob_score:
-            y_type = type_of_target(y)
+            y_type = type_of_target(y_encoded)
             if y_type in ("multiclass-multioutput", "unknown"):
-                # FIXME: we could consider to support multiclass-multioutput if
-                # we introduce or reuse a constructor parameter (e.g.
-                # oob_score) allowing our user to pass a callable defining the
-                # scoring strategy on OOB sample.
                 raise ValueError(
                     "The type of target cannot be used to compute OOB "
                     f"estimates. Got {y_type} while only the following are "
@@ -639,16 +715,6 @@ def fit(
                 )
             self._set_oob_score_and_attributes(X, y_encoded)
 
-        # Decapsulate classes_ attributes
-        if hasattr(self, "classes_") and self.n_outputs_ == 1:
-            self.n_classes_ = self.n_classes_[0]
-            self.classes_ = self.classes_[0]
-
-        # Print training infomation to console.
-        self._training_log_to_console()
-
-        return self
-
     def _set_oob_score_and_attributes(self, X, y):
         """Compute and set the OOB score and attributes.
         Parameters
diff --git a/imbens/ensemble/_under_sampling/self_paced_ensemble.py b/imbens/ensemble/_under_sampling/self_paced_ensemble.py
index 79dae42..10b2aac 100644
--- a/imbens/ensemble/_under_sampling/self_paced_ensemble.py
+++ b/imbens/ensemble/_under_sampling/self_paced_ensemble.py
@@ -24,12 +24,12 @@
         check_target_label_and_n_target_samples,
         check_train_verbose,
     )
-    from ..base import MAX_INT, BaseImbalancedEnsemble
+    from ..base import MAX_INT, BaseImbalancedEnsemble, _TrainingState
 else:  # pragma: no cover
     import sys  # For local test
 
     sys.path.append("../..")
-    from ensemble.base import BaseImbalancedEnsemble, MAX_INT
+    from ensemble.base import BaseImbalancedEnsemble, MAX_INT, _TrainingState
     from sampler._under_sampling import SelfPacedUnderSampler
     from utils._validation_data import check_eval_datasets
     from utils._validation_param import (
@@ -66,13 +66,21 @@
     'training_type': _training_type,
 }
 
+# Attributes stored in _TrainingState for backward-compatible delegation
+_TRAIN_STATE_ATTRS = frozenset({
+    'eval_datasets_', 'origin_distr_', 'target_label_', 'target_distr_',
+    'balancing_schedule_', 'eval_metrics_', 'train_verbose_',
+    'train_verbose_format_', '_seeds',
+    'y_pred_proba_latest', 'estimators_n_training_samples_',
+})
+
 
 @Substitution(
     random_state=_get_parameter_docstring('random_state', **_properties),
     n_jobs=_get_parameter_docstring('n_jobs', **_properties),
     example=_get_example_docstring(_method_name),
 )
-class SelfPacedEnsembleClassifier(BaseImbalancedEnsemble):
+class SelfPacedEnsembleClassifier(BaseImbalancedEnsemble):  # pylint: disable=too-many-instance-attributes
     """A self-paced ensemble (SPE) Classifier for class-imbalanced learning.
 
     Self-paced Ensemble (SPE) [1]_ is an ensemble learning framework for massive highly
@@ -201,6 +209,13 @@ def __init__(
         self.soft_resample_flag = soft_resample_flag
         self.replacement = replacement
 
+    def __getattr__(self, name):
+        if '_train_state_' in self.__dict__ and name in _TRAIN_STATE_ATTRS:
+            return getattr(self._train_state_, name)
+        raise AttributeError(
+            f"'{type(self).__name__}' object has no attribute '{name}'"
+        )
+
     @_deprecate_positional_args
     @FuncSubstitution(
         target_label=_get_parameter_docstring('target_label', **_properties),
@@ -277,16 +292,20 @@ def _fit(
             self.classes_,
         )
 
+        # Initialize training state container
+        self._train_state_ = _TrainingState()
+        ts = self._train_state_
+
         # Check evaluation data
         check_x_y_args = self.check_x_y_args
-        self.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
+        ts.eval_datasets_ = check_eval_datasets(eval_datasets, X, y, **check_x_y_args)
 
         # Check target sample strategy
         origin_distr_ = dict(Counter(y))
         target_label_, target_distr_ = check_target_label_and_n_target_samples(
             y, target_label, n_target_samples, self._sampling_type
         )
-        self.origin_distr_, self.target_label_, self.target_distr_ = (
+        ts.origin_distr_, ts.target_label_, ts.target_distr_ = (
             origin_distr_,
             target_label_,
             target_distr_,
@@ -294,13 +313,13 @@ def _fit(
 
         # Check balancing schedule
         balancing_schedule_ = check_balancing_schedule(balancing_schedule)
-        self.balancing_schedule_ = balancing_schedule_
+        ts.balancing_schedule_ = balancing_schedule_
 
         # Check evaluation metrics
-        self.eval_metrics_ = check_eval_metrics(eval_metrics)
+        ts.eval_metrics_ = check_eval_metrics(eval_metrics)
 
         # Check training train_verbose format
-        self.train_verbose_ = check_train_verbose(
+        ts.train_verbose_ = check_train_verbose(
             train_verbose, self.n_estimators, **self._properties
         )
 
@@ -310,12 +329,12 @@ def _fit(
         # Clear any previous fit results.
         self.estimators_ = []
         self.estimators_features_ = []
-        self.estimators_n_training_samples_ = np.zeros(n_estimators, dtype=int)
+        ts.estimators_n_training_samples_ = np.zeros(n_estimators, dtype=int)
         self.samplers_ = []
 
         # Genrate random seeds array
         seeds = random_state.randint(MAX_INT, size=n_estimators)
-        self._seeds = seeds
+        ts._seeds = seeds
 
         # Check if sample_weight is specified
         specified_sample_weight = sample_weight is not None
@@ -348,10 +367,10 @@ def _fit(
             resample_out = sampler.fit_resample(
                 X,
                 y,
-                y_pred_proba=self.y_pred_proba_latest,
+                y_pred_proba=ts.y_pred_proba_latest,
                 alpha=alpha,
                 classes_=classes_,
-                encode_map=self._encode_map,
+                encode_map=self._internal_state_['encode_map'],
                 sample_weight=sample_weight,
             )
 
@@ -367,8 +386,8 @@ def _fit(
                 (X_resampled, y_resampled) = resample_out
                 estimator.fit(X_resampled, y_resampled)
 
-            self.estimators_features_.append(self.features_)
-            self.estimators_n_training_samples_[i_iter] = y_resampled.shape[0]
+            self.estimators_features_.append(self._internal_state_['features'])
+            ts.estimators_n_training_samples_[i_iter] = y_resampled.shape[0]
 
             # Print training infomation to console.
             self._training_log_to_console(i_iter, y_resampled)
@@ -380,14 +399,15 @@ def _update_cached_prediction_probabilities(self, i_iter, X):
         data during ensemble training. Must be called in each iteration before fit the
         estimator."""
 
+        ts = self._train_state_
         if i_iter == 0:
-            self.y_pred_proba_latest = np.zeros(
-                (self._n_samples, self.n_classes_), dtype=np.float64
+            ts.y_pred_proba_latest = np.zeros(
+                (self._internal_state_['n_samples'], self.n_classes_), dtype=np.float64
             )
         else:
-            y_pred_proba_latest = self.y_pred_proba_latest
+            y_pred_proba_latest = ts.y_pred_proba_latest
             y_pred_proba_new = self.estimators_[-1].predict_proba(X)
-            self.y_pred_proba_latest = (
+            ts.y_pred_proba_latest = (
                 y_pred_proba_latest * i_iter + y_pred_proba_new
             ) / (i_iter + 1)
         return
diff --git a/imbens/ensemble/base.py b/imbens/ensemble/base.py
index 5e67ae0..0cd030f 100644
--- a/imbens/ensemble/base.py
+++ b/imbens/ensemble/base.py
@@ -55,6 +55,87 @@
 MAX_INT = np.iinfo(np.int32).max
 
 
+class _TrainingState:  # pylint: disable=too-many-instance-attributes
+    """Container for internal training state attributes.
+
+    Groups fit-time internal attributes to reduce the number of
+    top-level instance attributes on the ensemble classes.
+    """
+
+    __slots__ = (
+        '_y_encoded',
+        '_seeds',
+        '_n_samples',
+        '_max_samples',
+        '_max_features',
+        'raw_sample_weight_',
+        'sampler_kwargs_',
+        'balancing_schedule_',
+        '_encode_map',
+        'origin_distr_',
+        'target_distr_',
+        'target_label_',
+        'eval_datasets_',
+        'eval_metrics_',
+        'train_verbose_',
+        'train_verbose_format_',
+        'cost_matrix_',
+        'estimators_n_training_samples_',
+        'sample_weights_',
+        'keep_ratios_',
+        'y_pred_proba_latest',
+        'early_termination',
+    )
+
+    def __init__(
+        self,
+        _y_encoded=None,
+        _seeds=None,
+        _n_samples=None,
+        _max_samples=None,
+        _max_features=None,
+        raw_sample_weight_=None,
+        sampler_kwargs_=None,
+        balancing_schedule_=None,
+        _encode_map=None,
+        origin_distr_=None,
+        target_distr_=None,
+        target_label_=None,
+        eval_datasets_=None,
+        eval_metrics_=None,
+        train_verbose_=None,
+        train_verbose_format_=None,
+        cost_matrix_=None,
+        estimators_n_training_samples_=None,
+        sample_weights_=None,
+        keep_ratios_=None,
+        y_pred_proba_latest=None,
+        early_termination=None,
+    ):
+        self._y_encoded = _y_encoded
+        self._seeds = _seeds
+        self._n_samples = _n_samples
+        self._max_samples = _max_samples
+        self._max_features = _max_features
+        self.raw_sample_weight_ = raw_sample_weight_
+        self.sampler_kwargs_ = sampler_kwargs_
+        self.balancing_schedule_ = balancing_schedule_
+        self._encode_map = _encode_map
+        self.origin_distr_ = origin_distr_
+        self.target_distr_ = target_distr_
+        self.target_label_ = target_label_
+        self.eval_datasets_ = eval_datasets_
+        self.eval_metrics_ = eval_metrics_
+        self.train_verbose_ = train_verbose_
+        self.train_verbose_format_ = train_verbose_format_
+        self.cost_matrix_ = cost_matrix_
+        self.estimators_n_training_samples_ = estimators_n_training_samples_
+        self.sample_weights_ = sample_weights_
+        self.keep_ratios_ = keep_ratios_
+        self.y_pred_proba_latest = y_pred_proba_latest
+        self.early_termination = early_termination
+
+
 def sort_dict_by_key(d):
     """Sort a dict by key, return sorted dict."""
     return dict(sorted(d.items(), key=lambda k: k[0]))
@@ -81,6 +162,28 @@ class ImbalancedEnsembleClassifierMixin(ClassifierMixin):
 
     _estimator_ensemble_type = "imbens_classifier"
 
+    def _get_ts(self):
+        """Return the training state container if available, else self.
+
+        Allows the mixin methods to work with classes that use
+        _TrainingState (ReweightBoostClassifier hierarchy) and those
+        that store attributes directly on self.
+        """
+        return self._train_state_ if hasattr(self, '_train_state_') else self
+
+    def _compute_metrics(self, y_eval, y_predict_proba, eval_metrics, classes_):
+        scores = {}
+        for metric_name, (metric_func, kwargs, ac_proba, ac_labels) in eval_metrics.items():
+            if ac_labels:
+                kwargs["labels"] = classes_
+            if ac_proba:
+                score = metric_func(y_eval, y_predict_proba, **kwargs)
+            else:
+                y_predict = classes_.take(np.argmax(y_predict_proba, axis=1), axis=0)
+                score = metric_func(y_eval, y_predict, **kwargs)
+            scores[metric_name] = score
+        return scores
+
     def _evaluate(
         self,
         dataset_name: str,
@@ -92,48 +195,25 @@ def _evaluate(
         ensemble training process.
         """
 
-        eval_datasets_ = self.eval_datasets_
-        classes_ = self.classes_
-        verbose_format_ = self.train_verbose_format_
-
-        # Temporarily disable verbose
+        ts = self._get_ts()
         support_verbose = hasattr(self, "verbose")
         if support_verbose:
             verbose, self.verbose = self.verbose, 0
 
-        # If no eval_metrics is given, use self.eval_metrics_
-        if eval_metrics == None:
-            eval_metrics = self.eval_metrics_
+        if eval_metrics is None:
+            eval_metrics = ts.eval_metrics_
 
-        # If return numerical results
-        if return_value_dict == True:
+        if return_value_dict:
             value_dict = {}
-            for data_name, (X_eval, y_eval) in eval_datasets_.items():
+            for data_name, (X_eval, y_eval) in ts.eval_datasets_.items():
                 y_predict_proba = self.predict_proba(X_eval)
-                data_value_dict = {}
-                for metric_name, (
-                    metric_func,
-                    kwargs,
-                    ac_proba,
-                    ac_labels,
-                ) in eval_metrics.items():
-                    if ac_labels:
-                        kwargs["labels"] = classes_
-                    if ac_proba:  # If the metric take predict probabilities
-                        score = metric_func(y_eval, y_predict_proba, **kwargs)
-                    else:  # If the metric do not take predict probabilities
-                        y_predict = classes_.take(
-                            np.argmax(y_predict_proba, axis=1), axis=0
-                        )
-                        score = metric_func(y_eval, y_predict, **kwargs)
-                    data_value_dict[metric_name] = score
-                value_dict[data_name] = data_value_dict
+                value_dict[data_name] = self._compute_metrics(
+                    y_eval, y_predict_proba, eval_metrics, self.classes_
+                )
             out = value_dict
-
-        # If return string
         else:
             eval_info = ""
-            if return_header == True:
+            if return_header:
                 for metric_name in eval_metrics.keys():
                     eval_info = self._training_log_add_block(
                         eval_info,
@@ -141,39 +221,27 @@ def _evaluate(
                         "",
                         "",
                         " ",
-                        verbose_format_["len_metrics"][metric_name],
+                        ts.train_verbose_format_["len_metrics"][metric_name],
                         strip=False,
                     )
             else:
-                (X_eval, y_eval) = eval_datasets_[dataset_name]
+                X_eval, y_eval = ts.eval_datasets_[dataset_name]
                 y_predict_proba = self.predict_proba(X_eval)
-                for metric_name, (
-                    metric_func,
-                    kwargs,
-                    ac_proba,
-                    ac_labels,
-                ) in eval_metrics.items():
-                    if ac_labels:
-                        kwargs["labels"] = classes_
-                    if ac_proba:  # If the metric take predict probabilities
-                        score = metric_func(y_eval, y_predict_proba, **kwargs)
-                    else:  # If the metric do not take predict probabilities
-                        y_predict = classes_.take(
-                            np.argmax(y_predict_proba, axis=1), axis=0
-                        )
-                        score = metric_func(y_eval, y_predict, **kwargs)
+                metrics = self._compute_metrics(
+                    y_eval, y_predict_proba, eval_metrics, self.classes_
+                )
+                for metric_name in eval_metrics.keys():
                     eval_info = self._training_log_add_block(
                         eval_info,
-                        "{:.3f}".format(score),
+                        "{:.3f}".format(metrics[metric_name]),
                         "",
                         "",
                         " ",
-                        verbose_format_["len_metrics"][metric_name],
+                        ts.train_verbose_format_["len_metrics"][metric_name],
                         strip=False,
                     )
             out = eval_info[:-1]
 
-        # Recover verbose state
         if support_verbose:
             self.verbose = verbose
 
@@ -182,16 +250,17 @@ def _evaluate(
     def _init_training_log_format(self):
         """Private function for initialization of the training verbose format"""
 
-        if self.train_verbose_:
+        ts = self._get_ts()
+        if ts.train_verbose_:
             len_iter = (
                 max(len(str(self.n_estimators)), len(TRAINING_LOG_HEAD_TITLES["iter"]))
                 + 2
             )
-            if self.train_verbose_["print_distribution"]:
+            if ts.train_verbose_["print_distribution"]:
                 len_class_distr = (
                     max(
-                        len(str(self.target_distr_)),
-                        len(str(self.origin_distr_)),
+                        len(str(ts.target_distr_)),
+                        len(str(ts.origin_distr_)),
                         len(TRAINING_LOG_HEAD_TITLES["class_distr"]),
                     )
                     + 2
@@ -200,16 +269,16 @@ def _init_training_log_format(self):
                 len_class_distr = 0
             len_metrics = {
                 metric_name: max(len(metric_name), 5) + 2
-                for metric_name in self.eval_metrics_.keys()
+                for metric_name in ts.eval_metrics_.keys()
             }
             metrics_total_length = sum(len_metrics.values()) + len(len_metrics) - 1
             len_datasets = {
                 dataset_name: max(
                     metrics_total_length, len("Data: " + dataset_name) + 2
                 )
-                for dataset_name in self.eval_datasets_.keys()
+                for dataset_name in ts.eval_datasets_.keys()
             }
-            self.train_verbose_format_ = {
+            ts.train_verbose_format_ = {
                 "len_iter": len_iter,
                 "len_class_distr": len_class_distr,
                 "len_metrics": len_metrics,
@@ -235,21 +304,26 @@ def _training_log_add_line(
     ):
         """Private function for adding a line to training log."""
 
+        ts = self._get_ts()
         if texts == None:
-            texts = ("", "", tuple("" for _ in self.eval_datasets_.keys()))
+            texts = (
+                "",
+                "",
+                tuple("" for _ in ts.eval_datasets_.keys()),
+            )
         if tabs == None:
             tabs = ("┃", "┃", "┃", " ")
         if widths == None:
             widths = (
-                self.train_verbose_format_["len_iter"],
-                self.train_verbose_format_["len_class_distr"],
-                tuple(self.train_verbose_format_["len_datasets"].values()),
+                ts.train_verbose_format_["len_iter"],
+                ts.train_verbose_format_["len_class_distr"],
+                tuple(ts.train_verbose_format_["len_datasets"].values()),
             )
         if flags == None:
             flags = (
                 True,
-                self.train_verbose_["print_distribution"],
-                self.train_verbose_["print_metrics"],
+                ts.train_verbose_["print_distribution"],
+                ts.train_verbose_["print_metrics"],
             )
         (sta_char, mid_char, end_char, fill_char) = tabs
         (flag_iter, flag_distr, flag_metric) = flags
@@ -274,6 +348,7 @@ def _training_log_add_line(
     def _training_log_to_console_head(self):
         """Private function for printing a table header."""
 
+        ts = self._get_ts()
         # line 1
         info = (
             self._training_log_add_line(
@@ -289,7 +364,8 @@ def _training_log_to_console_head(self):
                     "",
                     "",
                     tuple(
-                        "Data: " + data_name for data_name in self.eval_datasets_.keys()
+                        "Data: " + data_name
+                        for data_name in ts.eval_datasets_.keys()
                     ),
                 ),
             )
@@ -302,7 +378,9 @@ def _training_log_to_console_head(self):
                 texts=(
                     TRAINING_LOG_HEAD_TITLES["iter"],
                     TRAINING_LOG_HEAD_TITLES["class_distr"],
-                    tuple("Metric" for data_name in self.eval_datasets_.keys()),
+                    tuple(
+                        "Metric" for data_name in ts.eval_datasets_.keys()
+                    ),
                 ),
             )
             + "\n"
@@ -316,7 +394,7 @@ def _training_log_to_console_head(self):
                     "",
                     tuple(
                         self._evaluate("", return_header=True)
-                        for data_name in self.eval_datasets_.keys()
+                        for data_name in ts.eval_datasets_.keys()
                     ),
                 ),
             )
@@ -330,7 +408,7 @@ def _training_log_to_console_head(self):
     def _training_log_to_console(self, i_iter=None, y=None):
         """Private function for printing training log to sys.stdout."""
 
-        if self.train_verbose_:
+        if self._get_ts().train_verbose_:
 
             if not hasattr(self, "_properties"):
                 raise AttributeError(
@@ -361,9 +439,10 @@ def _training_log_to_console_iterative(self, i_iter, y_resampled):
         if i_iter == 0:
             print(self._training_log_to_console_head())
 
-        eval_data_names = self.eval_datasets_.keys()
+        ts = self._get_ts()
+        eval_data_names = ts.eval_datasets_.keys()
 
-        if (i_iter + 1) % self.train_verbose_["granularity"] == 0 or i_iter == 0:
+        if (i_iter + 1) % ts.train_verbose_["granularity"] == 0 or i_iter == 0:
             print(
                 self._training_log_add_line(
                     texts=(
@@ -395,7 +474,7 @@ def _training_log_to_console_parallel(self):
         """Private function for printing training log to sys.stdout.
         (for ensemble classifiers that train in a parallel manner)"""
 
-        eval_data_names = self.eval_datasets_.keys()
+        eval_data_names = self._get_ts().eval_datasets_.keys()
         print(self._training_log_to_console_head())
         print(
             self._training_log_add_line(
@@ -443,7 +522,7 @@ def _parallel_decision_function(estimators, estimators_features, X):
     random_state=_get_parameter_docstring("random_state"),
     n_jobs=_get_parameter_docstring("n_jobs", **_properties),
 )
-class BaseImbalancedEnsemble(
+class BaseImbalancedEnsemble(  # pylint: disable=too-many-instance-attributes
     ImbalancedEnsembleClassifierMixin, BaseEnsemble, metaclass=ABCMeta
 ):
     """Base class for all imbalanced-ensemble classes that are
@@ -481,6 +560,38 @@ class BaseImbalancedEnsemble(
         The collection of fitted base estimators.
     """
 
+    @property
+    def random_state(self):
+        return self._config['random_state']
+
+    @random_state.setter
+    def random_state(self, value):
+        self._config['random_state'] = value
+
+    @property
+    def n_jobs(self):
+        return self._config['n_jobs']
+
+    @n_jobs.setter
+    def n_jobs(self, value):
+        self._config['n_jobs'] = value
+
+    @property
+    def verbose(self):
+        return self._config['verbose']
+
+    @verbose.setter
+    def verbose(self, value):
+        self._config['verbose'] = value
+
+    @property
+    def check_x_y_args(self):
+        return self._config['check_x_y_args']
+
+    @check_x_y_args.setter
+    def check_x_y_args(self, value):
+        self._config['check_x_y_args'] = value
+
     def __init__(
         self,
         estimator,
@@ -491,13 +602,15 @@ def __init__(
         verbose=0,
     ):
 
-        self.random_state = random_state
-        self.n_jobs = n_jobs
-        self.verbose = verbose
-        self.check_x_y_args = {
-            "accept_sparse": ["csr", "csc"],
-            "ensure_all_finite": False,
-            "dtype": None,
+        self._config = {
+            'random_state': random_state,
+            'n_jobs': n_jobs,
+            'verbose': verbose,
+            'check_x_y_args': {
+                "accept_sparse": ["csr", "csc"],
+                "ensure_all_finite": False,
+                "dtype": None,
+            },
         }
 
         super(BaseImbalancedEnsemble, self).__init__(
@@ -581,10 +694,12 @@ def fit(self, X, y, *, sample_weight=None, **kwargs):
 
         # Remap output
         n_samples, self.n_features_in_ = X.shape
-        self.features_ = np.arange(self.n_features_in_)
-        self._n_samples = n_samples
+        self._internal_state_ = {
+            'features': np.arange(self.n_features_in_),
+            'n_samples': n_samples,
+        }
         y = self._validate_y(y)
-        self._encode_map = {
+        self._internal_state_['encode_map'] = {
             c: np.where(self.classes_ == c)[0][0] for c in self.classes_
         }
 
diff --git a/imbens/metrics/_classification.py b/imbens/metrics/_classification.py
index c8a222d..6e240b9 100644
--- a/imbens/metrics/_classification.py
+++ b/imbens/metrics/_classification.py
@@ -37,145 +37,29 @@
 from sklearn.utils.multiclass import unique_labels
 from sklearn.utils.validation import check_consistent_length, column_or_1d
 
+_AVERAGE_OPTIONS = (None, "micro", "macro", "weighted", "samples")
 
-@_deprecate_positional_args
-def sensitivity_specificity_support(
-    y_true,
-    y_pred,
-    *,
-    labels=None,
-    pos_label=1,
-    average=None,
-    warn_for=("sensitivity", "specificity"),
-    sample_weight=None,
-):
-    """Compute sensitivity, specificity, and support for each class
-
-    The sensitivity is the ratio ``tp / (tp + fn)`` where ``tp`` is the number
-    of true positives and ``fn`` the number of false negatives. The sensitivity
-    quantifies the ability to avoid false negatives_[1].
-
-    The specificity is the ratio ``tn / (tn + fp)`` where ``tn`` is the number
-    of true negatives and ``fn`` the number of false negatives. The specificity
-    quantifies the ability to avoid false positives_[1].
-
-    The support is the number of occurrences of each class in ``y_true``.
-
-    If ``pos_label is None`` and in binary classification, this function
-    returns the average sensitivity and specificity if ``average``
-    is one of ``'weighted'``.
-
-    Read more in the `User Guide <https://imbalanced-learn.org/stable/metrics.html#sensitivity-specificity>`_.
 
-    Parameters
-    ----------
-    y_true : ndarray of shape (n_samples,)
-        Ground truth (correct) target values.
-
-    y_pred : ndarray of shape (n_samples,)
-        Estimated targets as returned by a classifier.
-
-    labels : list, default=None
-        The set of labels to include when ``average != 'binary'``, and their
-        order if ``average is None``. Labels present in the data can be
-        excluded, for example to calculate a multiclass average ignoring a
-        majority negative class, while labels not present in the data will
-        result in 0 components in a macro average. For multilabel targets,
-        labels are column indices. By default, all labels in ``y_true`` and
-        ``y_pred`` are used in sorted order.
-
-    pos_label : str or int, default=1
-        The class to report if ``average='binary'`` and the data is binary.
-        If the data are multiclass, this will be ignored;
-        setting ``labels=[pos_label]`` and ``average != 'binary'`` will report
-        scores for that label only.
-
-    average : str, default=None
-        If ``None``, the scores for each class are returned. Otherwise, this
-        determines the type of averaging performed on the data:
-
-        ``'binary'``:
-            Only report results for the class specified by ``pos_label``.
-            This is applicable only if targets (``y_{true,pred}``) are binary.
-        ``'micro'``:
-            Calculate metrics globally by counting the total true positives,
-            false negatives and false positives.
-        ``'macro'``:
-            Calculate metrics for each label, and find their unweighted
-            mean.  This does not take label imbalance into account.
-        ``'weighted'``:
-            Calculate metrics for each label, and find their average, weighted
-            by support (the number of true instances for each label). This
-            alters 'macro' to account for label imbalance; it can result in an
-            F-score that is not between precision and recall.
-        ``'samples'``:
-            Calculate metrics for each instance, and find their average (only
-            meaningful for multilabel classification where this differs from
-            :func:`accuracy_score`).
-
-    warn_for : tuple or set of {{"sensitivity", "specificity"}}, for internal use
-        This determines which warnings will be made in the case that this
-        function is being used to return only one of its metrics.
-
-    sample_weight : ndarray of shape (n_samples,), default=None
-        Sample weights.
-
-    Returns
-    -------
-    sensitivity : float (if `average is None`) or ndarray of \
-            shape (n_unique_labels,)
-        The sensitivity metric.
-
-    specificity : float (if `average is None`) or ndarray of \
-            shape (n_unique_labels,)
-        The specificity metric.
-
-    support : int (if `average is None`) or ndarray of \
-            shape (n_unique_labels,)
-        The number of occurrences of each label in ``y_true``.
-
-    References
-    ----------
-    .. [1] `Wikipedia entry for the Sensitivity and specificity
-           <https://en.wikipedia.org/wiki/Sensitivity_and_specificity>`_
-
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from imbens.metrics import sensitivity_specificity_support
-    >>> y_true = np.array(['cat', 'dog', 'pig', 'cat', 'dog', 'pig'])
-    >>> y_pred = np.array(['cat', 'pig', 'dog', 'cat', 'cat', 'dog'])
-    >>> sensitivity_specificity_support(y_true, y_pred, average='macro')
-    (0.33333333333333331, 0.66666666666666663, None)
-    >>> sensitivity_specificity_support(y_true, y_pred, average='micro')
-    (0.33333333333333331, 0.66666666666666663, None)
-    >>> sensitivity_specificity_support(y_true, y_pred, average='weighted')
-    (0.33333333333333331, 0.66666666666666663, None)
-    """
-    average_options = (None, "micro", "macro", "weighted", "samples")
+def _check_average_validity(average, average_options):
     if average not in average_options and average != "binary":
         raise ValueError("average has to be one of " + str(average_options))
 
-    y_type, y_true, y_pred = _check_targets(y_true, y_pred)
-    present_labels = unique_labels(y_true, y_pred)
 
+def _handle_binary_average(average, y_type, pos_label, present_labels, labels):
     if average == "binary":
         if y_type == "binary":
             if pos_label not in present_labels:
                 if len(present_labels) < 2:
-                    # Only negative labels
-                    return (0.0, 0.0, 0)
-                else:
-                    raise ValueError(
-                        "pos_label=%r is not a valid label: %r"
-                        % (pos_label, present_labels)
-                    )
-            labels = [pos_label]
-        else:
-            raise ValueError(
-                "Target is %s but average='binary'. Please "
-                "choose another average setting." % y_type
-            )
+                    return True, (0.0, 0.0, 0), None
+                raise ValueError(
+                    "pos_label=%r is not a valid label: %r"
+                    % (pos_label, present_labels)
+                )
+            return False, None, [pos_label]
+        raise ValueError(
+            "Target is %s but average='binary'. Please "
+            "choose another average setting." % y_type
+        )
     elif pos_label not in (None, 1):
         warnings.warn(
             "Note that pos_label (set to %r) is ignored when "
@@ -184,18 +68,22 @@ def sensitivity_specificity_support(
             % (pos_label, average),
             UserWarning,
         )
+    return False, None, labels
+
 
+def _resolve_labels(labels, present_labels):
     if labels is None:
-        labels = present_labels
-        n_labels = None
-    else:
-        n_labels = len(labels)
-        labels = np.hstack(
-            [labels, np.setdiff1d(present_labels, labels, assume_unique=True)]
-        )
+        return present_labels, None
+    n_labels = len(labels)
+    labels = np.hstack(
+        [labels, np.setdiff1d(present_labels, labels, assume_unique=True)]
+    )
+    return labels, n_labels
 
-    # Calculate tp_sum, pred_sum, true_sum ###
 
+def _compute_confusion_statistics(
+    y_true, y_pred, labels, n_labels, y_type, average, sample_weight
+):
     if y_type.startswith("multilabel"):
         raise ValueError("imblearn does not support multilabel")
     elif average == "samples":
@@ -204,56 +92,74 @@ def sensitivity_specificity_support(
             "not meaningful outside multilabel "
             "classification. See the accuracy_score instead."
         )
+
+    le = LabelEncoder()
+    le.fit(labels)
+    y_true = le.transform(y_true)
+    y_pred = le.transform(y_pred)
+    sorted_labels = le.classes_
+
+    tp = y_true == y_pred
+    tp_bins = y_true[tp]
+    if sample_weight is not None:
+        tp_bins_weights = np.asarray(sample_weight)[tp]
     else:
-        le = LabelEncoder()
-        le.fit(labels)
-        y_true = le.transform(y_true)
-        y_pred = le.transform(y_pred)
-        sorted_labels = le.classes_
+        tp_bins_weights = None
 
-        # labels are now from 0 to len(labels) - 1 -> use bincount
-        tp = y_true == y_pred
-        tp_bins = y_true[tp]
-        if sample_weight is not None:
-            tp_bins_weights = np.asarray(sample_weight)[tp]
-        else:
-            tp_bins_weights = None
+    if len(tp_bins):
+        tp_sum = np.bincount(
+            tp_bins, weights=tp_bins_weights, minlength=len(labels)
+        )
+    else:
+        true_sum = pred_sum = tp_sum = np.zeros(len(labels))
+    if len(y_pred):
+        pred_sum = np.bincount(y_pred, weights=sample_weight, minlength=len(labels))
+    if len(y_true):
+        true_sum = np.bincount(y_true, weights=sample_weight, minlength=len(labels))
 
-        if len(tp_bins):
-            tp_sum = np.bincount(
-                tp_bins, weights=tp_bins_weights, minlength=len(labels)
-            )
-        else:
-            # Pathological case
-            true_sum = pred_sum = tp_sum = np.zeros(len(labels))
-        if len(y_pred):
-            pred_sum = np.bincount(y_pred, weights=sample_weight, minlength=len(labels))
-        if len(y_true):
-            true_sum = np.bincount(y_true, weights=sample_weight, minlength=len(labels))
+    tn_sum = y_true.size - (pred_sum + true_sum - tp_sum)
 
-        # Compute the true negative
-        tn_sum = y_true.size - (pred_sum + true_sum - tp_sum)
+    indices = np.searchsorted(sorted_labels, labels[:n_labels])
+    tp_sum = tp_sum[indices]
+    true_sum = true_sum[indices]
+    pred_sum = pred_sum[indices]
+    tn_sum = tn_sum[indices]
+
+    return tp_sum, pred_sum, true_sum, tn_sum
 
-        # Retain only selected labels
-        indices = np.searchsorted(sorted_labels, labels[:n_labels])
-        tp_sum = tp_sum[indices]
-        true_sum = true_sum[indices]
-        pred_sum = pred_sum[indices]
-        tn_sum = tn_sum[indices]
 
+def _apply_micro_average(tp_sum, pred_sum, true_sum, tn_sum, average):
     if average == "micro":
         tp_sum = np.array([tp_sum.sum()])
         pred_sum = np.array([pred_sum.sum()])
         true_sum = np.array([true_sum.sum()])
         tn_sum = np.array([tn_sum.sum()])
+    return tp_sum, pred_sum, true_sum, tn_sum
 
-    # Finally, we have all our sufficient statistics. Divide! #
 
-    with np.errstate(divide="ignore", invalid="ignore"):
-        # Divide, and on zero-division, set scores to 0 and warn:
+def _compute_final_average(
+    sensitivity, specificity, average, true_sum, sample_weight
+):
+    if average == "weighted":
+        weights = true_sum
+        if weights.sum() == 0:
+            return 0, 0, None
+    elif average == "samples":
+        weights = sample_weight
+    else:
+        weights = None
 
-        # Oddly, we may get an "invalid" rather than a "divide" error
-        # here.
+    if average is not None:
+        assert average != "binary" or len(specificity) == 1
+        specificity = np.average(specificity, weights=weights)
+        sensitivity = np.average(sensitivity, weights=weights)
+        true_sum = None
+
+    return sensitivity, specificity, true_sum
+
+
+def _compute_sensitivity_specificity(tn_sum, pred_sum, tp_sum, true_sum, average, warn_for):
+    with np.errstate(divide="ignore", invalid="ignore"):
         specificity = _prf_divide(
             tn_sum,
             tn_sum + pred_sum - tp_sum,
@@ -265,23 +171,82 @@ def sensitivity_specificity_support(
         sensitivity = _prf_divide(
             tp_sum, true_sum, "sensitivity", "true", average, warn_for
         )
+    return sensitivity, specificity
 
-    # Average the results
 
-    if average == "weighted":
-        weights = true_sum
-        if weights.sum() == 0:
-            return 0, 0, None
-    elif average == "samples":
-        weights = sample_weight
-    else:
-        weights = None
+@_deprecate_positional_args
+def sensitivity_specificity_support(
+    y_true,
+    y_pred,
+    *,
+    labels=None,
+    pos_label=1,
+    average=None,
+    warn_for=("sensitivity", "specificity"),
+    sample_weight=None,
+):
+    """Compute sensitivity, specificity, and support for each class.
 
-    if average is not None:
-        assert average != "binary" or len(specificity) == 1
-        specificity = np.average(specificity, weights=weights)
-        sensitivity = np.average(sensitivity, weights=weights)
-        true_sum = None  # return no support
+    The sensitivity is the ratio ``tp / (tp + fn)`` and the specificity
+    is the ratio ``tn / (tn + fp)``.
+
+    Parameters
+    ----------
+    y_true : ndarray of shape (n_samples,)
+        Ground truth target values.
+    y_pred : ndarray of shape (n_samples,)
+        Estimated targets as returned by a classifier.
+    labels : list, default=None
+        The set of labels to include when ``average != 'binary'``.
+    pos_label : str or int, default=1
+        The class to report if ``average='binary'`` and the data is binary.
+    average : str, default=None
+        Averaging method: ``'binary'``, ``'micro'``, ``'macro'``,
+        ``'weighted'``, ``'samples'``, or ``None``.
+    warn_for : tuple, default=("sensitivity", "specificity")
+        Determines which warnings will be made.
+    sample_weight : ndarray of shape (n_samples,), default=None
+        Sample weights.
+
+    Returns
+    -------
+    sensitivity : float or ndarray
+    specificity : float or ndarray
+    support : int or ndarray
+
+    References
+    ----------
+    .. [1] Wikipedia entry for Sensitivity and specificity.
+    .. [2] https://imbalanced-learn.org/stable/metrics.html
+    """
+    _check_average_validity(average, _AVERAGE_OPTIONS)
+
+    y_type, y_true, y_pred = _check_targets(y_true, y_pred)
+    present_labels = unique_labels(y_true, y_pred)
+
+    early_return, result, labels = _handle_binary_average(
+        average, y_type, pos_label, present_labels, labels
+    )
+    if early_return:
+        return result
+
+    labels, n_labels = _resolve_labels(labels, present_labels)
+
+    tp_sum, pred_sum, true_sum, tn_sum = _compute_confusion_statistics(
+        y_true, y_pred, labels, n_labels, y_type, average, sample_weight
+    )
+
+    tp_sum, pred_sum, true_sum, tn_sum = _apply_micro_average(
+        tp_sum, pred_sum, true_sum, tn_sum, average
+    )
+
+    sensitivity, specificity = _compute_sensitivity_specificity(
+        tn_sum, pred_sum, tp_sum, true_sum, average, warn_for
+    )
+
+    sensitivity, specificity, true_sum = _compute_final_average(
+        sensitivity, specificity, average, true_sum, sample_weight
+    )
 
     return sensitivity, specificity, true_sum
 
diff --git a/imbens/pipeline.py b/imbens/pipeline.py
index de73f2a..8cd8943 100644
--- a/imbens/pipeline.py
+++ b/imbens/pipeline.py
@@ -1,4 +1,4 @@
-"""
+"""
 The :mod:`imbens.pipeline` module implements utilities to build a
 composite estimator, as a chain of transforms, samples and estimators.
 """
diff --git a/imbens/sampler/_over_sampling/_adasyn.py b/imbens/sampler/_over_sampling/_adasyn.py
index 7cd52a5..05f20e8 100644
--- a/imbens/sampler/_over_sampling/_adasyn.py
+++ b/imbens/sampler/_over_sampling/_adasyn.py
@@ -1,4 +1,4 @@
-"""Class to perform over-sampling using ADASYN."""
+"""Class to perform over-sampling using ADASYN."""
 
 # Adapted from imbalanced-learn
 
diff --git a/imbens/sampler/_over_sampling/_random_over_sampler.py b/imbens/sampler/_over_sampling/_random_over_sampler.py
index 557316c..41d250b 100644
--- a/imbens/sampler/_over_sampling/_random_over_sampler.py
+++ b/imbens/sampler/_over_sampling/_random_over_sampler.py
@@ -1,4 +1,4 @@
-"""Class to perform random over-sampling."""
+"""Class to perform random over-sampling."""
 
 # Adapted from imbalanced-learn
 
@@ -153,9 +153,7 @@ def _check_X_y(self, X, y):
         )
         return X, y, binarize_y
 
-    def _fit_resample(self, X, y, sample_weight=None):
-        random_state = check_random_state(self.random_state)
-
+    def _validate_shrinkage(self, X):
         if isinstance(self.shrinkage, Real):
             self.shrinkage_ = {
                 klass: self.shrinkage for klass in self.sampling_strategy_
@@ -169,34 +167,86 @@ def _fit_resample(self, X, y, sample_weight=None):
                 f"Got {repr(self.shrinkage)} instead."
             )
 
-        if self.shrinkage_ is not None:
-            missing_shrinkage_keys = (
-                self.sampling_strategy_.keys() - self.shrinkage_.keys()
+        if self.shrinkage_ is None:
+            return X
+
+        missing_shrinkage_keys = (
+            self.sampling_strategy_.keys() - self.shrinkage_.keys()
+        )
+        if missing_shrinkage_keys:
+            raise ValueError(
+                f"`shrinkage` should contain a shrinkage factor for "
+                f"each class that will be resampled. The missing "
+                f"classes are: {repr(missing_shrinkage_keys)}"
             )
-            if missing_shrinkage_keys:
+
+        for klass, shrink_factor in self.shrinkage_.items():
+            if shrink_factor < 0:
                 raise ValueError(
-                    f"`shrinkage` should contain a shrinkage factor for "
-                    f"each class that will be resampled. The missing "
-                    f"classes are: {repr(missing_shrinkage_keys)}"
+                    f"The shrinkage factor needs to be >= 0. "
+                    f"Got {shrink_factor} for class {klass}."
                 )
 
-            for klass, shrink_factor in self.shrinkage_.items():
-                if shrink_factor < 0:
-                    raise ValueError(
-                        f"The shrinkage factor needs to be >= 0. "
-                        f"Got {shrink_factor} for class {klass}."
-                    )
+        try:
+            X = check_array(X, accept_sparse=["csr", "csc"], dtype="numeric")
+        except ValueError as exc:
+            raise ValueError(
+                "When shrinkage is not None, X needs to contain only "
+                "numerical data to later generate a smoothed bootstrap "
+                "sample."
+            ) from exc
+        return X
+
+    def _make_smoothed_bootstrap(self, X, target_class_indices, bootstrap_indices, num_samples, class_sample, random_state):
+        n_samples, n_features = X.shape
+        smoothing_constant = (4 / ((n_features + 2) * n_samples)) ** (
+            1 / (n_features + 4)
+        )
+        if sparse.issparse(X):
+            _, X_class_variance = mean_variance_axis(
+                X[target_class_indices, :],
+                axis=0,
+            )
+            X_class_scale = np.sqrt(X_class_variance, out=X_class_variance)
+        else:
+            X_class_scale = np.std(X[target_class_indices, :], axis=0)
+        smoothing_matrix = np.diagflat(
+            self.shrinkage_[class_sample] * smoothing_constant * X_class_scale
+        )
+        X_new = random_state.randn(num_samples, n_features)
+        X_new = X_new.dot(smoothing_matrix) + X[bootstrap_indices, :]
+        if sparse.issparse(X):
+            X_new = sparse.csr_matrix(X_new, dtype=X.dtype)
+        return X_new
 
-            # smoothed bootstrap imposes to make numerical operation; we need
-            # to be sure to have only numerical data in X
-            try:
-                X = check_array(X, accept_sparse=["csr", "csc"], dtype="numeric")
-            except ValueError as exc:
-                raise ValueError(
-                    "When shrinkage is not None, X needs to contain only "
-                    "numerical data to later generate a smoothed bootstrap "
-                    "sample."
-                ) from exc
+    def _init_resampling(self, X):
+        random_state = check_random_state(self.random_state)
+        X = self._validate_shrinkage(X)
+        return random_state, X
+
+    def _concatenate_resampled(self, X_resampled, y_resampled, X):
+        if sparse.issparse(X):
+            X_resampled = sparse.vstack(X_resampled, format=X.format)
+        else:
+            X_resampled = np.vstack(X_resampled)
+        y_resampled = np.hstack(y_resampled)
+        return X_resampled, y_resampled
+
+    def _build_sample_weight_result(self, sample_weight, y_resampled, y):
+        sample_weight_new = np.empty(
+            y_resampled.shape[0] - y.shape[0], dtype=np.float64
+        )
+        sample_weight_new[:] = np.mean(sample_weight)
+        sample_weight_resampled = np.hstack(
+            [sample_weight, sample_weight_new]
+        ).reshape(-1, 1)
+        sample_weight_resampled = np.squeeze(
+            normalize(sample_weight_resampled, axis=0, norm="l1")
+        )
+        return sample_weight_resampled
+
+    def _fit_resample(self, X, y, sample_weight=None):
+        random_state, X = self._init_resampling(X)
 
         X_resampled = [X.copy()]
         y_resampled = [y.copy()]
@@ -216,29 +266,13 @@ def _fit_resample(self, X, y, sample_weight=None):
             )
             sample_indices = np.append(sample_indices, bootstrap_indices)
             if self.shrinkage_ is not None:
-                # generate a smoothed bootstrap with a perturbation
-                n_samples, n_features = X.shape
-                smoothing_constant = (4 / ((n_features + 2) * n_samples)) ** (
-                    1 / (n_features + 4)
-                )
-                if sparse.issparse(X):
-                    _, X_class_variance = mean_variance_axis(
-                        X[target_class_indices, :],
-                        axis=0,
+                X_resampled.append(
+                    self._make_smoothed_bootstrap(
+                        X, target_class_indices, bootstrap_indices,
+                        num_samples, class_sample, random_state,
                     )
-                    X_class_scale = np.sqrt(X_class_variance, out=X_class_variance)
-                else:
-                    X_class_scale = np.std(X[target_class_indices, :], axis=0)
-                smoothing_matrix = np.diagflat(
-                    self.shrinkage_[class_sample] * smoothing_constant * X_class_scale
                 )
-                X_new = random_state.randn(num_samples, n_features)
-                X_new = X_new.dot(smoothing_matrix) + X[bootstrap_indices, :]
-                if sparse.issparse(X):
-                    X_new = sparse.csr_matrix(X_new, dtype=X.dtype)
-                X_resampled.append(X_new)
             else:
-                # generate a bootstrap
                 X_resampled.append(_safe_indexing(X, bootstrap_indices))
 
             y_resampled.append(_safe_indexing(y, bootstrap_indices))
@@ -249,29 +283,16 @@ def _fit_resample(self, X, y, sample_weight=None):
                 )
 
         self.sample_indices_ = np.array(sample_indices)
+        X_resampled, y_resampled = self._concatenate_resampled(
+            X_resampled, y_resampled, X
+        )
 
-        if sparse.issparse(X):
-            X_resampled = sparse.vstack(X_resampled, format=X.format)
-        else:
-            X_resampled = np.vstack(X_resampled)
-        y_resampled = np.hstack(y_resampled)
-
-        # If given sample_weight
         if sample_weight_flag:
-            # sample_weight is already validated in self.fit_resample()
-            sample_weight_new = np.empty(
-                y_resampled.shape[0] - y.shape[0], dtype=np.float64
-            )
-            sample_weight_new[:] = np.mean(sample_weight)
-            sample_weight_resampled = np.hstack(
-                [sample_weight, sample_weight_new]
-            ).reshape(-1, 1)
-            sample_weight_resampled = np.squeeze(
-                normalize(sample_weight_resampled, axis=0, norm="l1")
+            sample_weight_resampled = self._build_sample_weight_result(
+                sample_weight, y_resampled, y
             )
             return X_resampled, y_resampled, sample_weight_resampled
-        else:
-            return X_resampled, y_resampled
+        return X_resampled, y_resampled
 
     def _more_tags(self):  # pragma: no cover
         return {
diff --git a/imbens/sampler/_over_sampling/_smote/cluster.py b/imbens/sampler/_over_sampling/_smote/cluster.py
index 15a5e6d..074d62a 100644
--- a/imbens/sampler/_over_sampling/_smote/cluster.py
+++ b/imbens/sampler/_over_sampling/_smote/cluster.py
@@ -211,6 +211,87 @@ def _find_cluster_sparsity(self, X):
         )
         return (mean_distance ** exponent) / X.shape[0]
 
+    def _get_valid_clusters_and_sparsities(
+        self, X, y, X_clusters, class_sample, n_samples, total_inp_samples
+    ):
+        valid_clusters = []
+        cluster_sparsities = []
+
+        for cluster_idx in range(self.kmeans_estimator_.n_clusters):
+            cluster_mask = np.flatnonzero(X_clusters == cluster_idx)
+            X_cluster = _safe_indexing(X, cluster_mask)
+            y_cluster = _safe_indexing(y, cluster_mask)
+
+            if cluster_mask.sum() == 0:
+                continue
+
+            cluster_class_mean = (y_cluster == class_sample).mean()
+            balance_threshold = (
+                n_samples / total_inp_samples / 2
+                if self.cluster_balance_threshold_ == "auto"
+                else self.cluster_balance_threshold_
+            )
+
+            if cluster_class_mean < balance_threshold:
+                continue
+
+            anticipated_samples = cluster_class_mean * X_cluster.shape[0]
+            if anticipated_samples < self.nn_k_.n_neighbors:
+                continue
+
+            X_cluster_class = _safe_indexing(
+                X_cluster, np.flatnonzero(y_cluster == class_sample)
+            )
+            valid_clusters.append(cluster_mask)
+            cluster_sparsities.append(
+                self._find_cluster_sparsity(X_cluster_class)
+            )
+
+        return valid_clusters, np.array(cluster_sparsities)
+
+    def _sample_from_cluster(
+        self, X, y, valid_cluster_idx, valid_cluster, class_sample,
+        n_samples, cluster_n_samples_list, cluster_weights,
+        X_resampled, y_resampled
+    ):
+        X_cluster = _safe_indexing(X, valid_cluster)
+        y_cluster = _safe_indexing(y, valid_cluster)
+
+        X_cluster_class = _safe_indexing(
+            X_cluster, np.flatnonzero(y_cluster == class_sample)
+        )
+
+        self.nn_k_.fit(X_cluster_class)
+        nns = self.nn_k_.kneighbors(X_cluster_class, return_distance=False)[
+            :, 1:
+        ]
+
+        cluster_n_samples = (
+            int(n_samples - sum(cluster_n_samples_list))
+            if valid_cluster_idx == self.kmeans_estimator_.n_clusters - 1
+            else math.floor(
+                n_samples * cluster_weights[valid_cluster_idx]
+            )
+        )
+
+        cluster_n_samples_list[valid_cluster_idx] = cluster_n_samples
+
+        X_new, y_new = self._make_samples(
+            X_cluster_class,
+            y.dtype,
+            class_sample,
+            X_cluster_class,
+            nns,
+            cluster_n_samples,
+            1.0,
+        )
+
+        stack = [np.vstack, sparse.vstack][int(sparse.issparse(X_new))]
+        X_resampled = stack((X_resampled, X_new))
+        y_resampled = np.hstack((y_resampled, y_new))
+
+        return X_resampled, y_resampled
+
     def _fit_resample(self, X, y, sample_weight=None):
         self._validate_estimator()
         X_resampled = X.copy()
@@ -222,44 +303,12 @@ def _fit_resample(self, X, y, sample_weight=None):
                 continue
 
             X_clusters = self.kmeans_estimator_.fit_predict(X)
-            valid_clusters = []
-            cluster_sparsities = []
-
-            # identify cluster which are answering the requirements
-            for cluster_idx in range(self.kmeans_estimator_.n_clusters):
-
-                cluster_mask = np.flatnonzero(X_clusters == cluster_idx)
-                X_cluster = _safe_indexing(X, cluster_mask)
-                y_cluster = _safe_indexing(y, cluster_mask)
-                
-                # empty cluster
-                if cluster_mask.sum() == 0:
-                    continue
-
-                cluster_class_mean = (y_cluster == class_sample).mean()
-
-                if self.cluster_balance_threshold_ == "auto":
-                    balance_threshold = n_samples / total_inp_samples / 2
-                else:
-                    balance_threshold = self.cluster_balance_threshold_
-
-                # the cluster is already considered balanced
-                if cluster_class_mean < balance_threshold:
-                    continue
-
-                # not enough samples to apply SMOTE
-                anticipated_samples = cluster_class_mean * X_cluster.shape[0]
-                if anticipated_samples < self.nn_k_.n_neighbors:
-                    continue
-
-                X_cluster_class = _safe_indexing(
-                    X_cluster, np.flatnonzero(y_cluster == class_sample)
+            valid_clusters, cluster_sparsities = (
+                self._get_valid_clusters_and_sparsities(
+                    X, y, X_clusters, class_sample, n_samples, total_inp_samples
                 )
+            )
 
-                valid_clusters.append(cluster_mask)
-                cluster_sparsities.append(self._find_cluster_sparsity(X_cluster_class))
-
-            cluster_sparsities = np.array(cluster_sparsities)
             cluster_weights = cluster_sparsities / cluster_sparsities.sum()
             cluster_n_samples_list = np.zeros_like(cluster_weights)
 
@@ -272,51 +321,29 @@ def _fit_resample(self, X, y, sample_weight=None):
                 )
 
             for valid_cluster_idx, valid_cluster in enumerate(valid_clusters):
-                X_cluster = _safe_indexing(X, valid_cluster)
-                y_cluster = _safe_indexing(y, valid_cluster)
-
-                X_cluster_class = _safe_indexing(
-                    X_cluster, np.flatnonzero(y_cluster == class_sample)
+                X_resampled, y_resampled = self._sample_from_cluster(
+                    X, y, valid_cluster_idx, valid_cluster, class_sample,
+                    n_samples, cluster_n_samples_list, cluster_weights,
+                    X_resampled, y_resampled
                 )
 
-                self.nn_k_.fit(X_cluster_class)
-                nns = self.nn_k_.kneighbors(X_cluster_class, return_distance=False)[
-                    :, 1:
-                ]
-
-                if valid_cluster_idx == self.kmeans_estimator_.n_clusters - 1:
-                    cluster_n_samples = int(n_samples - sum(cluster_n_samples_list))
-                else:
-                    cluster_n_samples = math.floor(
-                        n_samples * cluster_weights[valid_cluster_idx]
-                    )
-                cluster_n_samples_list[valid_cluster_idx] = cluster_n_samples
-
-                X_new, y_new = self._make_samples(
-                    X_cluster_class,
-                    y.dtype,
-                    class_sample,
-                    X_cluster_class,
-                    nns,
-                    cluster_n_samples,
-                    1.0,
-                )
+        return self._handle_sample_weight(
+            sample_weight, X_resampled, y_resampled, y
+        )
 
-                stack = [np.vstack, sparse.vstack][int(sparse.issparse(X_new))]
-                X_resampled = stack((X_resampled, X_new))
-                y_resampled = np.hstack((y_resampled, y_new))
-        
-        # If given sample_weight
+    def _handle_sample_weight(self, sample_weight, X_resampled, y_resampled, y):
         if sample_weight is not None:
-            # sample_weight is already validated in self.fit_resample()
             sample_weight_new = \
                 np.empty(y_resampled.shape[0] - y.shape[0], dtype=np.float64)
             sample_weight_new[:] = np.mean(sample_weight)
-            sample_weight_resampled = np.hstack([sample_weight, sample_weight_new]).reshape(-1, 1)
+            sample_weight_resampled = np.hstack(
+                [sample_weight, sample_weight_new]
+            ).reshape(-1, 1)
             sample_weight_resampled = \
                 np.squeeze(normalize(sample_weight_resampled, axis=0, norm='l1'))
             return X_resampled, y_resampled, sample_weight_resampled
-        else: return X_resampled, y_resampled
+
+        return X_resampled, y_resampled
 
 
 # %%
diff --git a/imbens/sampler/_over_sampling/_smote/filter.py b/imbens/sampler/_over_sampling/_smote/filter.py
index b86862b..390ddf5 100644
--- a/imbens/sampler/_over_sampling/_smote/filter.py
+++ b/imbens/sampler/_over_sampling/_smote/filter.py
@@ -1,4 +1,4 @@
-"""SMOTE variant applying some filtering before the generation process."""
+"""SMOTE variant applying some filtering before the generation process."""
 # Adapted from imbalanced-learn
 
 # Authors: Guillaume Lemaitre
@@ -397,13 +397,16 @@ def _fit_resample(self, X, y, sample_weight=None):
             self.nn_k_.fit(X_class)
             fractions = random_state.beta(10, 10)
             n_generated_samples = int(fractions * (n_samples + 1))
+            X_new_list = []
+            y_new_list = []
+
             if np.count_nonzero(danger_bool) > 0:
                 nns = self.nn_k_.kneighbors(
                     _safe_indexing(support_vector, np.flatnonzero(danger_bool)),
                     return_distance=False,
                 )[:, 1:]
 
-                X_new_1, y_new_1 = self._make_samples(
+                X_new, y_new = self._make_samples(
                     _safe_indexing(support_vector, np.flatnonzero(danger_bool)),
                     y.dtype,
                     class_sample,
@@ -412,6 +415,8 @@ def _fit_resample(self, X, y, sample_weight=None):
                     n_generated_samples,
                     step_size=1.0,
                 )
+                X_new_list.append(X_new)
+                y_new_list.append(y_new)
 
             if np.count_nonzero(safety_bool) > 0:
                 nns = self.nn_k_.kneighbors(
@@ -419,7 +424,7 @@ def _fit_resample(self, X, y, sample_weight=None):
                     return_distance=False,
                 )[:, 1:]
 
-                X_new_2, y_new_2 = self._make_samples(
+                X_new, y_new = self._make_samples(
                     _safe_indexing(support_vector, np.flatnonzero(safety_bool)),
                     y.dtype,
                     class_sample,
@@ -428,25 +433,15 @@ def _fit_resample(self, X, y, sample_weight=None):
                     n_samples - n_generated_samples,
                     step_size=-self.out_step,
                 )
+                X_new_list.append(X_new)
+                y_new_list.append(y_new)
 
-            if np.count_nonzero(danger_bool) > 0 and np.count_nonzero(safety_bool) > 0:
-                if sparse.issparse(X_resampled):
-                    X_resampled = sparse.vstack([X_resampled, X_new_1, X_new_2])
-                else:
-                    X_resampled = np.vstack((X_resampled, X_new_1, X_new_2))
-                y_resampled = np.concatenate((y_resampled, y_new_1, y_new_2), axis=0)
-            elif np.count_nonzero(danger_bool) == 0:
-                if sparse.issparse(X_resampled):
-                    X_resampled = sparse.vstack([X_resampled, X_new_2])
-                else:
-                    X_resampled = np.vstack((X_resampled, X_new_2))
-                y_resampled = np.concatenate((y_resampled, y_new_2), axis=0)
-            elif np.count_nonzero(safety_bool) == 0:
+            if X_new_list:
                 if sparse.issparse(X_resampled):
-                    X_resampled = sparse.vstack([X_resampled, X_new_1])
+                    X_resampled = sparse.vstack([X_resampled] + X_new_list)
                 else:
-                    X_resampled = np.vstack((X_resampled, X_new_1))
-                y_resampled = np.concatenate((y_resampled, y_new_1), axis=0)
+                    X_resampled = np.vstack([X_resampled] + X_new_list)
+                y_resampled = np.concatenate([y_resampled] + y_new_list, axis=0)
 
         # If given sample_weight
         if sample_weight is not None:
diff --git a/imbens/sampler/_under_sampling/_prototype_selection/_edited_nearest_neighbours.py b/imbens/sampler/_under_sampling/_prototype_selection/_edited_nearest_neighbours.py
index 74e723d..5d0bc16 100644
--- a/imbens/sampler/_under_sampling/_prototype_selection/_edited_nearest_neighbours.py
+++ b/imbens/sampler/_under_sampling/_prototype_selection/_edited_nearest_neighbours.py
@@ -298,13 +298,13 @@ def _validate_estimator(self):
                 f" Got {type(self.max_iter)} instead."
             )
 
-        self.nn_ = check_neighbors_object(
+        nn = check_neighbors_object(
             "n_neighbors", self.n_neighbors, additional_neighbor=1
         )
 
         self.enn_ = EditedNearestNeighbours(
             sampling_strategy=self.sampling_strategy,
-            n_neighbors=self.nn_,
+            n_neighbors=nn,
             kind_sel=self.kind_sel,
             n_jobs=self.n_jobs,
         )
diff --git a/imbens/sampler/_under_sampling/_prototype_selection/_random_under_sampler.py b/imbens/sampler/_under_sampling/_prototype_selection/_random_under_sampler.py
index d58ccf5..06a062f 100644
--- a/imbens/sampler/_under_sampling/_prototype_selection/_random_under_sampler.py
+++ b/imbens/sampler/_under_sampling/_prototype_selection/_random_under_sampler.py
@@ -102,33 +102,34 @@ def _check_X_y(self, X, y):
         )
         return X, y, binarize_y
 
-    def _fit_resample(self, X, y, sample_weight=None, sample_proba=None):
-        random_state = check_random_state(self.random_state)
-
+    def _validate_sample_proba(self, sample_proba, y):
         if sample_proba is None:
-            pass
-        elif not isinstance(sample_proba, (np.ndarray, list)):
+            return None
+        if not isinstance(sample_proba, (np.ndarray, list)):
             raise TypeError(
                 f"`sample_proba` should be an array-like of shape (n_samples,),"
                 f" got {type(sample_proba)} instead."
             )
-        else:
-            sample_proba = np.asarray(sample_proba)
-            if sample_proba.shape != y.shape:
-                raise ValueError(
-                    f"`sample_proba` should be of shape {y.shape}, got {sample_proba.shape}."
-                )
-            else:
-                try:
-                    sample_proba = sample_proba.astype(float)
-                except Exception as e:
-                    e_args = list(e.args)
-                    e_args[0] += (
-                        f"\n`sample_proba` should be an array-like with dtype == float,"
-                        + f" please check your usage."
-                    )
-                    e.args = tuple(e_args)
-                    raise e
+        sample_proba = np.asarray(sample_proba)
+        if sample_proba.shape != y.shape:
+            raise ValueError(
+                f"`sample_proba` should be of shape {y.shape}, got {sample_proba.shape}."
+            )
+        try:
+            sample_proba = sample_proba.astype(float)
+        except Exception as e:
+            e_args = list(e.args)
+            e_args[0] += (
+                f"\n`sample_proba` should be an array-like with dtype == float,"
+                + f" please check your usage."
+            )
+            e.args = tuple(e_args)
+            raise e
+        return sample_proba
+
+    def _fit_resample(self, X, y, sample_weight=None, sample_proba=None):
+        random_state = check_random_state(self.random_state)
+        sample_proba = self._validate_sample_proba(sample_proba, y)
 
         idx_under = np.empty((0,), dtype=int)
 
diff --git a/imbens/sampler/_under_sampling/_prototype_selection/_self_paced_under_sampler.py b/imbens/sampler/_under_sampling/_prototype_selection/_self_paced_under_sampler.py
index 79e4522..ab0cddb 100644
--- a/imbens/sampler/_under_sampling/_prototype_selection/_self_paced_under_sampler.py
+++ b/imbens/sampler/_under_sampling/_prototype_selection/_self_paced_under_sampler.py
@@ -252,24 +252,9 @@ def _fit_resample(
         else:
             return _safe_indexing(X, index_spu), _safe_indexing(y, index_spu)
 
-    def _undersample_single_class(
-        self, hardness_c, n_target_samples_c, index_c, alpha, random_state
-    ):
-        """Perform self-paced under-sampling in a single class"""
+    def _compute_bin_info(self, hardness_c, n_target_samples_c, alpha):
         k_bins = self.k_bins
-        soft_resample_flag = self.soft_resample_flag
-        replacement = self.replacement
-        n_samples_c = hardness_c.shape[0]
-
-        # if hardness is not distinguishable or no sample will be dropped
-        if hardness_c.max() == hardness_c.min() or n_target_samples_c == n_samples_c:
-            # perform random under-sampling
-            return random_state.choice(
-                index_c, size=n_target_samples_c, replace=replacement
-            )
-
         with np.errstate(divide="ignore", invalid="ignore"):
-            # compute population & hardness contribution of each bin
             populations, edges = np.histogram(hardness_c, bins=k_bins)
             contributions = np.zeros(k_bins)
             index_bins = []
@@ -283,11 +268,11 @@ def _undersample_single_class(
                 if populations[i_bin] > 0:
                     contributions[i_bin] = hardness_c[index_bin].mean()
 
-            # compute the expected number of samples to be sampled from each bin
             bin_weights = 1.0 / (contributions + alpha)
             bin_weights[np.isnan(bin_weights) | np.isinf(bin_weights)] = 0
-            n_target_samples_bins = n_target_samples_c * bin_weights / bin_weights.sum()
-            # check whether exists empty bins
+            n_target_samples_bins = (
+                n_target_samples_c * bin_weights / bin_weights.sum()
+            )
             n_invalid_samples = sum(n_target_samples_bins[populations == 0])
             if n_invalid_samples > 0:
                 n_valid_samples = n_target_samples_c - n_invalid_samples
@@ -296,62 +281,108 @@ def _undersample_single_class(
             n_target_samples_bins = n_target_samples_bins.astype(int)
             while True:
                 for i in np.flip(np.argsort(populations)):
-                    n_target_diff = n_target_samples_c - n_target_samples_bins.sum()
+                    n_target_diff = (
+                        n_target_samples_c - n_target_samples_bins.sum()
+                    )
                     if n_target_diff == 0:
                         break
                     elif populations[i] > 0:
                         n_target_samples_bins[i] += 1
                 if n_target_diff == 0:
                     break
-            assert n_target_samples_c == n_target_samples_bins.sum()
-
-        if soft_resample_flag:
-            with np.errstate(divide="ignore", invalid="ignore"):
-                # perform soft (weighted) self-paced under-sampling
-                soft_spu_bin_weights = n_target_samples_bins / populations
-                soft_spu_bin_weights[~np.isfinite(soft_spu_bin_weights)] = 0
-                # print ('soft_spu_bin_weights: ', soft_spu_bin_weights)
-            # compute sampling probabilities
-            soft_spu_sample_proba = np.zeros_like(hardness_c)
-            for i_bin in range(k_bins):
-                soft_spu_sample_proba[index_bins[i_bin]] = soft_spu_bin_weights[i_bin]
-            soft_spu_sample_proba /= soft_spu_sample_proba.sum()
-            # sample with respect to the sampling probabilities
+
+        return index_bins, populations, n_target_samples_bins
+
+    def _soft_resample_single_class(
+        self,
+        hardness_c,
+        index_bins,
+        n_target_samples_bins,
+        populations,
+        index_c,
+        n_target_samples_c,
+        random_state,
+    ):
+        k_bins = self.k_bins
+        with np.errstate(divide="ignore", invalid="ignore"):
+            soft_spu_bin_weights = n_target_samples_bins / populations
+            soft_spu_bin_weights[~np.isfinite(soft_spu_bin_weights)] = 0
+        soft_spu_sample_proba = np.zeros_like(hardness_c)
+        for i_bin in range(k_bins):
+            soft_spu_sample_proba[index_bins[i_bin]] = soft_spu_bin_weights[i_bin]
+        soft_spu_sample_proba /= soft_spu_sample_proba.sum()
+        return random_state.choice(
+            index_c,
+            size=n_target_samples_c,
+            replace=self.replacement,
+            p=soft_spu_sample_proba,
+        )
+
+    def _hard_resample_single_class(
+        self,
+        index_bins,
+        n_target_samples_bins,
+        populations,
+        index_c,
+        random_state,
+    ):
+        k_bins = self.k_bins
+        index_c_results = []
+        for i_bin in range(k_bins):
+            if (
+                populations[i_bin] < n_target_samples_bins[i_bin]
+                and not self.replacement
+            ):
+                raise RuntimeError(
+                    f"Met {i_bin}-th bin with insufficient number of data samples"
+                    f" ({populations[i_bin]}, expected"
+                    f" >= {n_target_samples_bins[i_bin]})."
+                    f" Set 'soft_resample_flag' or 'replacement' to `True` to."
+                    f" avoid this issue."
+                )
+            index_c_bin = index_c[index_bins[i_bin]]
+            if len(index_c_bin) > 0:
+                index_c_results.append(
+                    random_state.choice(
+                        index_c_bin,
+                        size=n_target_samples_bins[i_bin],
+                        replace=self.replacement,
+                    )
+                )
+        return np.hstack(index_c_results)
+
+    def _undersample_single_class(
+        self, hardness_c, n_target_samples_c, index_c, alpha, random_state
+    ):
+        n_samples_c = hardness_c.shape[0]
+
+        if hardness_c.max() == hardness_c.min() or n_target_samples_c == n_samples_c:
             return random_state.choice(
+                index_c, size=n_target_samples_c, replace=self.replacement
+            )
+
+        index_bins, populations, n_target_samples_bins = self._compute_bin_info(
+            hardness_c, n_target_samples_c, alpha
+        )
+
+        if self.soft_resample_flag:
+            return self._soft_resample_single_class(
+                hardness_c,
+                index_bins,
+                n_target_samples_bins,
+                populations,
                 index_c,
-                size=n_target_samples_c,
-                replace=replacement,
-                p=soft_spu_sample_proba,
+                n_target_samples_c,
+                random_state,
             )
         else:
-            # perform hard self-paced under-sampling
-            index_c_results = []
-            for i_bin in range(k_bins):
-                # if no sufficient data in bin for resampling, raise an RuntimeError
-                if (
-                    populations[i_bin] < n_target_samples_bins[i_bin]
-                    and not replacement
-                ):
-                    raise RuntimeError(
-                        f"Met {i_bin}-th bin with insufficient number of data samples"
-                        f" ({populations[i_bin]}, expected"
-                        f" >= {n_target_samples_bins[i_bin]})."
-                        f" Set 'soft_resample_flag' or 'replacement' to `True` to."
-                        f" avoid this issue."
-                    )
-                index_c_bin = index_c[index_bins[i_bin]]
-                # random sample from each bin
-                if len(index_c_bin) > 0:
-                    index_c_results.append(
-                        random_state.choice(
-                            index_c_bin,
-                            size=n_target_samples_bins[i_bin],
-                            replace=replacement,
-                        )
-                    )
-            # concatenate and return the result
-            index_c_result = np.hstack(index_c_results)
-            return index_c_result
+            return self._hard_resample_single_class(
+                index_bins,
+                n_target_samples_bins,
+                populations,
+                index_c,
+                random_state,
+            )
 
     def _more_tags(self):  # pragma: no cover
         return {
diff --git a/imbens/sampler/_under_sampling/_prototype_selection/tests/test_balance_cascade_under_sampler.py b/imbens/sampler/_under_sampling/_prototype_selection/tests/test_balance_cascade_under_sampler.py
index e643e40..55abe3e 100644
--- a/imbens/sampler/_under_sampling/_prototype_selection/tests/test_balance_cascade_under_sampler.py
+++ b/imbens/sampler/_under_sampling/_prototype_selection/tests/test_balance_cascade_under_sampler.py
@@ -95,6 +95,20 @@
 sampling_strategy_normal = {2: 6, 1: 4, 0: 2}
 sampling_strategy_error = {2: 200, 1: 100, 0: 90}
 
+_SAMPLING_STRATEGY_MAP = {
+    'org': sampling_strategy_org,
+    'auto': 'auto',
+    'default': sampling_strategy_default,
+    'normal': sampling_strategy_normal,
+    'error': sampling_strategy_error,
+}
+
+_KEEP_POP_MAP = {
+    'org': sampling_strategy_org,
+    'default': sampling_strategy_default,
+    'normal': sampling_strategy_normal,
+}
+
 
 @pytest.mark.parametrize("as_frame", [True, False], ids=["dataframe", "array"])
 @pytest.mark.parametrize("proba_det", [True, False], ids=["det", "normal"])
@@ -113,70 +127,41 @@ def test_spus_fit_resample(as_frame, proba_det, sampling, keep_pop):
     else:
         X_ = X
 
-    # prepare y_pred_proba
-    if proba_det:
-        y_pred_proba = pred_proba_det
-    else:
-        y_pred_proba = pred_proba_normal
-
-    if sampling == 'org':
-        sampling_strategy = sampling_strategy_org
-    elif sampling == 'auto':
-        sampling_strategy = "auto"
-    elif sampling == 'default':
-        sampling_strategy = sampling_strategy_default
-    elif sampling == 'normal':
-        sampling_strategy = sampling_strategy_normal
-    elif sampling == 'error':
-        sampling_strategy = sampling_strategy_error
-
-    if keep_pop == 'org':
-        keep_populations = sampling_strategy_org
-    elif keep_pop == 'default':
-        keep_populations = sampling_strategy_default
-    elif keep_pop == 'normal':
-        keep_populations = sampling_strategy_normal
+    y_pred_proba = pred_proba_det if proba_det else pred_proba_normal
+    sampling_strategy = _SAMPLING_STRATEGY_MAP[sampling]
+    keep_populations = _KEEP_POP_MAP[keep_pop]
 
-    # init resampler
     bcus = BalanceCascadeUnderSampler(
         sampling_strategy=sampling_strategy,
         replacement=False,
         random_state=RND_SEED,
     )
 
-    # resampling
     if sampling == 'error':
         with pytest.raises(ValueError, match="With under-sampling methods"):
-            X_resampled, y_resampled, new_dropped_index = bcus.fit_resample(
-                X_,
-                y,
-                y_pred_proba=y_pred_proba,
-                classes_=classes_,
-                encode_map=encode_map,
-                dropped_index=dropped_index,
+            bcus.fit_resample(
+                X_, y, y_pred_proba=y_pred_proba, classes_=classes_,
+                encode_map=encode_map, dropped_index=dropped_index,
                 keep_populations=keep_populations,
             )
-    else:
-        X_resampled, y_resampled, new_dropped_index = bcus.fit_resample(
-            X_,
-            y,
-            y_pred_proba=y_pred_proba,
-            classes_=classes_,
-            encode_map=encode_map,
-            dropped_index=dropped_index,
-            keep_populations=keep_populations,
-        )
+        return
+
+    X_resampled, y_resampled, new_dropped_index = bcus.fit_resample(
+        X_, y, y_pred_proba=y_pred_proba, classes_=classes_,
+        encode_map=encode_map, dropped_index=dropped_index,
+        keep_populations=keep_populations,
+    )
+
+    if as_frame and hasattr(X_resampled, "loc"):
+        X_resampled = X_resampled.to_numpy()
 
-        if as_frame:
-            if hasattr(X_resampled, "loc"):
-                X_resampled = X_resampled.to_numpy()
+    if sampling == 'auto':
+        sampling_strategy = sampling_strategy_default
 
-        if sampling == 'auto':
-            sampling_strategy = sampling_strategy_default
-        assert X_resampled.shape == (sum(sampling_strategy.values()), 2)
-        assert y_resampled.shape == (sum(sampling_strategy.values()),)
-        assert dict(Counter(y_resampled)) == sampling_strategy
-        assert dict(Counter(y[~new_dropped_index])) == keep_populations
+    assert X_resampled.shape == (sum(sampling_strategy.values()), 2)
+    assert y_resampled.shape == (sum(sampling_strategy.values()),)
+    assert dict(Counter(y_resampled)) == sampling_strategy
+    assert dict(Counter(y[~new_dropped_index])) == keep_populations
 
 
 def test_bcus_no_sufficient_data():
diff --git a/imbens/sampler/base.py b/imbens/sampler/base.py
index 41fae54..78e618b 100644
--- a/imbens/sampler/base.py
+++ b/imbens/sampler/base.py
@@ -1,4 +1,4 @@
-"""Base class for sampling
+"""Base class for sampling
 """
 
 # Adapted from imbalanced-learn
diff --git a/imbens/utils/_validation_param.py b/imbens/utils/_validation_param.py
index 53c4821..f39e054 100644
--- a/imbens/utils/_validation_param.py
+++ b/imbens/utils/_validation_param.py
@@ -99,9 +99,64 @@ def _check_n_target_samples_int(y, n_target_samples, sampling_type):
         raise SamplingKindError
 
 
+def _check_n_target_samples_dict_under(target_stats, n_target_samples):
+    target_distr = copy(target_stats)
+    for class_label, n_target_sample in n_target_samples.items():
+        n_origin_sample = target_stats[class_label]
+        if n_target_sample > n_origin_sample:
+            raise ValueError(
+                f" The target number of samples of class {class_label}"
+                f" should be < {n_origin_sample} (number of samples"
+                f" in class {class_label}) to perform under-sampling,"
+                f" got {n_target_sample}."
+            )
+        else:
+            target_distr[class_label] = n_target_sample
+    return target_distr
+
+
+def _check_n_target_samples_dict_over(target_stats, n_target_samples):
+    target_distr = copy(target_stats)
+    for class_label, n_target_sample in n_target_samples.items():
+        n_origin_sample = target_stats[class_label]
+        if n_target_sample < n_origin_sample:
+            raise ValueError(
+                f" The target number of samples of class {class_label}"
+                f" should be > {n_origin_sample} (number of samples"
+                f" in class {class_label}) to perform over-sampling,"
+                f" got {n_target_sample}."
+            )
+        else:
+            target_distr[class_label] = n_target_sample
+    return target_distr
+
+
+def _check_n_target_samples_dict_hybrid(target_stats, n_target_samples):
+    target_distr = copy(target_stats)
+    if all(
+        n_target_samples[label] <= target_stats[label]
+        for label in n_target_samples.keys()
+    ):
+        raise Warning(
+            f"The target number of samples is smaller than the number"
+            f" of original samples for all classes. ONLY under-sampling"
+            f" will be carried out."
+        )
+    elif all(
+        n_target_samples[label] >= target_stats[label]
+        for label in n_target_samples.keys()
+    ):
+        raise Warning(
+            f"The target number of samples is greater than the number"
+            f" of original samples for all classes. ONLY over-sampling"
+            f" will be carried out."
+        )
+    target_distr.update(n_target_samples)
+    return target_distr
+
+
 def _check_n_target_samples_dict(y, n_target_samples, sampling_type):
     target_stats = dict(Counter(y))
-    # check that all keys in n_target_samples are also in y
     set_diff_sampling_strategy_target = set(n_target_samples.keys()) - set(
         target_stats.keys()
     )
@@ -110,7 +165,6 @@ def _check_n_target_samples_dict(y, n_target_samples, sampling_type):
             f"The {set_diff_sampling_strategy_target} target class is/are not "
             f"present in the data."
         )
-    # check that there is no negative number
     if any(n_samples <= 0 for n_samples in n_target_samples.values()):
         raise ValueError(
             f"The number of samples in a class must > 0. "
@@ -118,58 +172,11 @@ def _check_n_target_samples_dict(y, n_target_samples, sampling_type):
         )
 
     if sampling_type == 'under-sampling':
-        target_distr = copy(target_stats)
-        for class_label, n_target_sample in n_target_samples.items():
-            n_origin_sample = target_stats[class_label]
-            if n_target_sample > n_origin_sample:
-                raise ValueError(
-                    f" The target number of samples of class {class_label}"
-                    f" should be < {n_origin_sample} (number of samples"
-                    f" in class {class_label}) to perform under-sampling,"
-                    f" got {n_target_sample}."
-                )
-            else:
-                target_distr[class_label] = n_target_sample
-        return target_distr
-
+        return _check_n_target_samples_dict_under(target_stats, n_target_samples)
     elif sampling_type == 'over-sampling':
-        target_distr = copy(target_stats)
-        for class_label, n_target_sample in n_target_samples.items():
-            n_origin_sample = target_stats[class_label]
-            if n_target_sample < n_origin_sample:
-                raise ValueError(
-                    f" The target number of samples of class {class_label}"
-                    f" should be > {n_origin_sample} (number of samples"
-                    f" in class {class_label}) to perform over-sampling,"
-                    f" got {n_target_sample}."
-                )
-            else:
-                target_distr[class_label] = n_target_sample
-        return target_distr
-
+        return _check_n_target_samples_dict_over(target_stats, n_target_samples)
     elif sampling_type == "multi-class-hybrid-sampling":
-        target_distr = copy(target_stats)
-        if all(
-            n_target_samples[label] <= target_stats[label]
-            for label in n_target_samples.keys()
-        ):
-            raise Warning(
-                f"The target number of samples is smaller than the number"
-                f" of original samples for all classes. ONLY under-sampling"
-                f" will be carried out."
-            )
-        elif all(
-            n_target_samples[label] >= target_stats[label]
-            for label in n_target_samples.keys()
-        ):
-            raise Warning(
-                f"The target number of samples is greater than the number"
-                f" of original samples for all classes. ONLY over-sampling"
-                f" will be carried out."
-            )
-        target_distr.update(n_target_samples)
-        return target_distr
-
+        return _check_n_target_samples_dict_hybrid(target_stats, n_target_samples)
     else:
         raise SamplingKindError
 
@@ -490,76 +497,71 @@ def check_eval_metrics(eval_metrics):
 )
 
 
-def check_train_verbose(
-    train_verbose: bool or numbers.Integral or dict,
-    n_estimators_ensemble: int,
-    training_type: str,
-    **ignored_properties,
-):
-    # n_estimators_ensemble:int,):
-
-    train_verbose_ = copy(TRAIN_VERBOSE_DEFAULT)
-    train_verbose_.update({'granularity': max(1, int(n_estimators_ensemble / 10))})
+def _check_train_verbose_parallel(train_verbose, n_estimators_ensemble):
+    if isinstance(train_verbose, bool):
+        if train_verbose:
+            train_verbose_ = copy(TRAIN_VERBOSE_DEFAULT)
+            train_verbose_.update({'granularity': max(1, int(n_estimators_ensemble / 10))})
+            train_verbose_['print_distribution'] = False
+            return train_verbose_
+        return False
+    raise TypeError(
+        f"'train_verbose' can only be of type `bool`"
+        f" for ensemble classifiers trained in parallel,"
+        f" gor {type(train_verbose)}."
+    )
 
-    if training_type == 'parallel':
-        # For ensemble classifiers trained in parallel
-        # train_verbose can only be of type bool
-        if isinstance(train_verbose, bool):
-            if train_verbose == True:
-                train_verbose_['print_distribution'] = False
-                return train_verbose_
-            if train_verbose == False:
-                return False
-        else:
-            raise TypeError(
-                f"'train_verbose' can only be of type `bool`"
-                f" for ensemble classifiers trained in parallel,"
-                f" gor {type(train_verbose)}."
-            )
 
-    elif training_type == 'iterative':
-        # For ensemble classifiers trained in iterative manner
-        # train_verbose can be of type bool / int / dict
-        if isinstance(train_verbose, bool):
-            if train_verbose == True:
-                return train_verbose_
-            if train_verbose == False:
-                return False
-
-        if isinstance(train_verbose, numbers.Integral):
-            train_verbose_.update({'granularity': train_verbose})
+def _check_train_verbose_iterative(train_verbose, train_verbose_):
+    if isinstance(train_verbose, bool):
+        if train_verbose:
             return train_verbose_
+        return False
+
+    if isinstance(train_verbose, numbers.Integral):
+        train_verbose_.update({'granularity': train_verbose})
+        return train_verbose_
 
-        if isinstance(train_verbose, dict):
-            # check key value type
-            set_diff_verbose_keys = set(train_verbose.keys()) - set(
-                TRAIN_VERBOSE_TYPE.keys()
+    if isinstance(train_verbose, dict):
+        set_diff_verbose_keys = set(train_verbose.keys()) - set(TRAIN_VERBOSE_TYPE.keys())
+        if set_diff_verbose_keys:
+            raise ValueError(
+                f"'train_verbose' keys {set_diff_verbose_keys} are not supported."
+                + TRAIN_VERBOSE_DICT_INFO
             )
-            if len(set_diff_verbose_keys) > 0:
-                raise ValueError(
-                    f"'train_verbose' keys {set_diff_verbose_keys} are not supported."
+        for key, value in train_verbose.items():
+            if not isinstance(value, TRAIN_VERBOSE_TYPE[key]):
+                raise TypeError(
+                    f"train_verbose['{key}'] has wrong data type, should be {TRAIN_VERBOSE_TYPE[key]}."
                     + TRAIN_VERBOSE_DICT_INFO
                 )
-            for key, value in train_verbose.items():
-                if not isinstance(value, TRAIN_VERBOSE_TYPE[key]):
-                    raise TypeError(
-                        f"train_verbose['{key}'] has wrong data type, should be {TRAIN_VERBOSE_TYPE[key]}."
-                        + TRAIN_VERBOSE_DICT_INFO
-                    )
-            train_verbose_.update(train_verbose)
-            return train_verbose_
+        train_verbose_.update(train_verbose)
+        return train_verbose_
 
-        else:
-            raise TypeError(
-                f"'train_verbose' should be of type `bool`, `int`, or `dict`, got {type(train_verbose)} instead."
-                + TRAIN_VERBOSE_DICT_INFO
-            )
+    raise TypeError(
+        f"'train_verbose' should be of type `bool`, `int`, or `dict`, got {type(train_verbose)} instead."
+        + TRAIN_VERBOSE_DICT_INFO
+    )
 
-    else:
-        raise NotImplementedError(
-            f"'check_train_verbose' for 'training_type' = {training_type}"
-            f" needs to be implemented."
-        )
+
+def check_train_verbose(
+    train_verbose: bool or numbers.Integral or dict,
+    n_estimators_ensemble: int,
+    training_type: str,
+    **ignored_properties,
+):
+    if training_type == 'parallel':
+        return _check_train_verbose_parallel(train_verbose, n_estimators_ensemble)
+
+    if training_type == 'iterative':
+        train_verbose_ = copy(TRAIN_VERBOSE_DEFAULT)
+        train_verbose_.update({'granularity': max(1, int(n_estimators_ensemble / 10))})
+        return _check_train_verbose_iterative(train_verbose, train_verbose_)
+
+    raise NotImplementedError(
+        f"'check_train_verbose' for 'training_type' = {training_type}"
+        f" needs to be implemented."
+    )
 
 
 VISUALIZER_ENSEMBLES_EXAMPLE_INFO = " Example: {..., ensemble_name: ensemble, ...}"
diff --git a/imbens/visualizer/visualizer.py b/imbens/visualizer/visualizer.py
index 4b2bb2a..a94fd8a 100644
--- a/imbens/visualizer/visualizer.py
+++ b/imbens/visualizer/visualizer.py
@@ -98,6 +98,19 @@ def set_ax_border(ax, border_color="black", border_width=2):
     return ax
 
 
+class _FitData:
+    """Container for data produced during ImbalancedEnsembleVisualizer.fit()."""
+
+    def __init__(self):
+        self.ensembles_ = None
+        self.vis_format_ = None
+        self.ensembles_n_training_samples_ = None
+        self.ensembles_has_n_training_samples_ = None
+        self.granularity_ = None
+        self.perf_dataframe_ = None
+        self.conf_matrices_ = None
+
+
 class ImbalancedEnsembleVisualizer:
     """A visualizer, providing several utilities to visualize:
 
@@ -218,29 +231,34 @@ def __init__(self, eval_metrics: dict = None, eval_datasets: dict = None):
             )
 
         # Default styles of different titles
-        self.suptitle_style = {
-            "size": "x-large",
-            "stretch": "expanded",
-            "style": "italic",
-            "variant": "small-caps",
-            "weight": "black",
-        }
-        self.row_col_title_style = {
-            "size": "large",
-            "weight": "bold",
-            "bbox": {
-                "boxstyle": "round",
-                "pad": 0.25,
-                "fc": "white",
-                "ec": "black",
-                "lw": 1,
-                "alpha": 0.5,
+        self._title_styles = {
+            'suptitle': {
+                "size": "x-large",
+                "stretch": "expanded",
+                "style": "italic",
+                "variant": "small-caps",
+                "weight": "black",
+            },
+            'row_col': {
+                "size": "large",
+                "weight": "bold",
+                "bbox": {
+                    "boxstyle": "round",
+                    "pad": 0.25,
+                    "fc": "white",
+                    "ec": "black",
+                    "lw": 1,
+                    "alpha": 0.5,
+                },
+            },
+            'axis': {
+                "size": "medium",
+                "weight": "bold",
             },
         }
-        self.axis_title_style = {
-            "size": "medium",
-            "weight": "bold",
-        }
+
+        # Container for data produced during fit
+        self._fit_data = _FitData()
 
     def fit(
         self,
@@ -284,13 +302,13 @@ def fit(
         eval_datasets_, eval_metrics_ = self.eval_datasets_, self.eval_metrics_
         # eval_datasets_ needs to be further validated in check_visualizer_ensembles
         (
-            self.ensembles_,
+            self._fit_data.ensembles_,
             self.eval_datasets_,
-            self.vis_format_,
+            self._fit_data.vis_format_,
         ) = check_visualizer_ensembles(ensembles, eval_datasets_, eval_metrics_)
         (
-            self.ensembles_n_training_samples_,
-            self.ensembles_has_n_training_samples_,
+            self._fit_data.ensembles_n_training_samples_,
+            self._fit_data.ensembles_has_n_training_samples_,
         ) = self._check_ensembles_stored_n_training_samples()
 
         # Check granularity
@@ -298,14 +316,14 @@ def fit(
             granularity_ = check_type(granularity, "granularity", numbers.Integral)
         else:
             max_n_estimators = max(
-                [len(ens.estimators_) for ens in self.ensembles_.values()]
+                [len(ens.estimators_) for ens in self._fit_data.ensembles_.values()]
             )
             granularity_ = max(int(max_n_estimators / 5), 1)
-        self.granularity_ = granularity_
+        self._fit_data.granularity_ = granularity_
 
         # Collect data for visualization
-        self.perf_dataframe_ = self._collect_all_ensemble_performance_data()
-        self.conf_matrices_ = self._collect_all_ensemble_confusion_matrix()
+        self._fit_data.perf_dataframe_ = self._collect_all_ensemble_performance_data()
+        self._fit_data.conf_matrices_ = self._collect_all_ensemble_confusion_matrix()
 
         self._fitted = True
 
@@ -330,7 +348,7 @@ def _check_ensembles_stored_n_training_samples(self):
 
         ensembles_has_n_training_samples = {}
         ensembles_n_training_samples = {}
-        for name, ensemble in self.ensembles_.items():
+        for name, ensemble in self._fit_data.ensembles_.items():
             (
                 estimators_n_training_samples_,
                 flag,
@@ -472,14 +490,14 @@ def _collect_ensemble_performance_data(
         """Private function for collecting performance data of a single ensemble model."""
 
         # Set local variables
-        granularity = self.granularity_
+        granularity = self._fit_data.granularity_
         eval_metrics = self.eval_metrics_
         classes_ = ensemble.classes_
         max_n_estimators = len(ensemble.estimators_)
-        estimators_n_training_samples_ = self.ensembles_n_training_samples_[
+        estimators_n_training_samples_ = self._fit_data.ensembles_n_training_samples_[
             ensemble_name
         ]
-        has_n_training_samples_ = self.ensembles_has_n_training_samples_[ensemble_name]
+        has_n_training_samples_ = self._fit_data.ensembles_has_n_training_samples_[ensemble_name]
 
         results = []
         # For each evaluation dataset
@@ -537,15 +555,15 @@ def _collect_all_ensemble_performance_data(self):
         """Private function for collecting performance data of all ensemble models."""
 
         # max_len_xxx_name was used to format the verbose
-        max_len_ensemble_name = max([len(_) for _ in self.ensembles_.keys()])
-        max_len_dataset_name = max([len(_) for _ in self.vis_format_["dataset_names"]])
+        max_len_ensemble_name = max([len(_) for _ in self._fit_data.ensembles_.keys()])
+        max_len_dataset_name = max([len(_) for _ in self._fit_data.vis_format_["dataset_names"]])
 
         return pd.concat(
             [
                 self._collect_ensemble_performance_data(
                     ensemble_name, ensemble, max_len_ensemble_name, max_len_dataset_name
                 )
-                for ensemble_name, ensemble in self.ensembles_.items()
+                for ensemble_name, ensemble in self._fit_data.ensembles_.items()
             ]
         )
 
@@ -554,7 +572,7 @@ def _split_ensembles_by_has_n_training_samples_(self, on_ensembles_):
         the array of numbers of training samples of each base estimator.
         """
 
-        ensembles_has_n_training_samples = self.ensembles_has_n_training_samples_
+        ensembles_has_n_training_samples = self._fit_data.ensembles_has_n_training_samples_
         positives, negatives = [], []
         for ensemble_name in on_ensembles_:
             if ensembles_has_n_training_samples[ensemble_name]:
@@ -564,71 +582,12 @@ def _split_ensembles_by_has_n_training_samples_(self, on_ensembles_):
 
         return tuple(positives), tuple(negatives)
 
-    def performance_lineplot(
-        self,
-        on_ensembles: list = None,
-        on_datasets: list = None,
-        on_metrics: list = None,
-        split_by: list = [],
-        n_samples_as_x_axis: bool = False,
-        sub_figsize: tuple = (4.0, 3.3),
-        sup_title: bool or str = True,
-        **lineplot_kwargs,
+    def _validate_performance_lineplot_params(
+        self, on_ensembles, on_datasets, on_metrics, split_by,
+        n_samples_as_x_axis, sub_figsize, sup_title, lineplot_kwargs
     ):
-        """Draw a performance line plot.
-
-        Parameters
-        ----------
-        on_ensembles : list of strings, default=None
-            The names of ensembles to include in the plot. It should be a
-            subset of ``self.ensembles_.keys()``. if ``None``, all ensembles
-            fitted by the visualizer will be included.
-
-        on_datasets : list of strings, default=None
-            The names of evaluation datasets to include in the plot. It
-            should be a subset of ``self.eval_datasets_.keys()``. if ``None``,
-            all evaluation datasets will be included.
-
-        on_metrics : list of strings, default=None
-            The names of evaluation metrics to include in the plot. It
-            should be a subset of ``self.eval_metrics_.keys()``. if ``None``,
-            all evaluation metrics will be included.
-
-        split_by : list of {'method', 'dataset'}, default=[]
-            How to group the results for visualization.
-
-            - if contains ``'method'``, the performance results of different
-              ensemble methods will be displayed in independent sub-figures.
-            - if contains ``'dataset'``, the performance results on different
-              evaluation datasets will be displayed in independent sub-figures.
-
-        n_samples_as_x_axis : bool, default=False
-            Whether to use the number of training samples as the x-axis.
-
-        sub_figsize: (float, float), default=(4.0, 3.3)
-            The size of an subfigure (width, height in inches).
-            The overall figure size will be automatically determined by
-            (sub_figsize[0] * num_columns, sub_figsize[1] * num_rows).
-
-        sup_title: bool or str, default=True
-            The super title of the figure.
-
-            - if ``True``, automatically determines the super title.
-            - if ``False``, no super title will be displayed.
-            - if ``string``, super title will be ``sup_title``.
-
-        **lineplot_kwargs : key, value mappings
-            Other keyword arguments are passed down to
-            :meth:`seaborn.lineplot`.
+        vis_perf_dataframe = self._fit_data.perf_dataframe_.copy()
 
-        Returns
-        -------
-        self : object
-        """
-
-        vis_perf_dataframe = self.perf_dataframe_.copy()
-
-        # Check if the visualizer is fitted
         if not self._fitted:
             raise NotFittedError(
                 f"This visualizer is not fitted yet."
@@ -636,21 +595,16 @@ def performance_lineplot(
                 f" 'performance_lineplot'."
             )
 
-        # Check parameters
         on_ensembles = self._check_is_subset(
-            on_ensembles, "on_ensembles", self.vis_format_["ensemble_names"]
+            on_ensembles, "on_ensembles", self._fit_data.vis_format_["ensemble_names"]
         )
         on_datasets = self._check_is_subset(
-            on_datasets, "on_datasets", self.vis_format_["dataset_names"]
+            on_datasets, "on_datasets", self._fit_data.vis_format_["dataset_names"]
         )
         on_metrics = self._check_is_subset(
-            on_metrics, "on_metrics", self.vis_format_["metric_names"]
-        )
-        n_ensembles, n_datasets, n_metrics = (
-            len(on_ensembles),
-            len(on_datasets),
-            len(on_metrics),
+            on_metrics, "on_metrics", self._fit_data.vis_format_["metric_names"]
         )
+        n_datasets, n_metrics = len(on_datasets), len(on_metrics)
 
         if not isinstance(split_by, list):
             raise TypeError(
@@ -681,50 +635,56 @@ def performance_lineplot(
                     f" performance_lineplot function."
                 )
 
-        # Select data for visualization
         on_ensembles_mask = vis_perf_dataframe["method"].map(
-            {k: k in on_ensembles for k in self.vis_format_["ensemble_names"]}
+            {k: k in on_ensembles for k in self._fit_data.vis_format_["ensemble_names"]}
         )
         on_datasets_mask = vis_perf_dataframe["dataset"].map(
-            {k: k in on_datasets for k in self.vis_format_["dataset_names"]}
+            {k: k in on_datasets for k in self._fit_data.vis_format_["dataset_names"]}
         )
         on_metrics_mask = vis_perf_dataframe["metric"].map(
-            {k: k in on_metrics for k in self.vis_format_["metric_names"]}
+            {k: k in on_metrics for k in self._fit_data.vis_format_["metric_names"]}
         )
         final_mask = on_ensembles_mask & on_datasets_mask & on_metrics_mask
         vis_perf_dataframe = vis_perf_dataframe.loc[final_mask]
 
-        # Set subfigure x-axis (# estimators or # training samples)
+        return (vis_perf_dataframe, on_ensembles, on_metrics, n_datasets,
+                n_metrics, n_samples_as_x_axis_, sub_fig_width,
+                sub_fig_height, sup_title_, lineplot_kwargs_)
+
+    def _configure_x_axis(
+        self, vis_perf_dataframe, on_ensembles, n_samples_as_x_axis_
+    ):
+        n_ensembles = len(on_ensembles)
+
         if n_samples_as_x_axis_:
             (
                 include_ensembles,
                 exclude_ensembles,
             ) = self._split_ensembles_by_has_n_training_samples_(on_ensembles)
-            # If any estimator does not have `estimators_n_training_samples_` attribute
             if len(exclude_ensembles) > 0:
-                # Raise a warning of excluding some methods from visualization
                 warn(
                     f"scikit-learn ensemble estimator(s) with name"
                     + f" {exclude_ensembles} does not record the number"
                     + f" training samples, they will be excluded from the"
                     + f" figure when `n_samples_as_x_axis` = `True`."
                 )
-                # Update the number of included ensemble methods
                 n_ensembles = len(include_ensembles)
-                # Select data
                 has_n_samples_mask = vis_perf_dataframe["method"].map(
-                    self.ensembles_has_n_training_samples_
+                    self._fit_data.ensembles_has_n_training_samples_
                 )
                 vis_perf_dataframe = vis_perf_dataframe.loc[has_n_samples_mask]
-            # Set x-axis to number of samples
             x_column = "n_samples"
             x_label = "# Training Samples"
         else:
-            # Set x-axis to number of base estimators
             x_column = "n_estimators"
             x_label = "# Base Estimators"
 
-        # Set figure size and layout
+        return x_column, x_label, vis_perf_dataframe, n_ensembles
+
+    def _setup_figure_layout(
+        self, n_ensembles, n_datasets, n_metrics, split_by,
+        sub_fig_width, sub_fig_height, sup_title_, on_metrics
+    ):
         n_rows_fig, n_columns_fig = 1, n_metrics
         if "method" in split_by:
             n_rows_fig *= n_ensembles
@@ -735,13 +695,16 @@ def performance_lineplot(
             sub_fig_width * n_columns_fig,
             sub_fig_height * n_rows_fig,
         )
-        # If has sup_title, add reserved space
         total_height += RESERVED_SUPTITLE_INCHES if sup_title_ else 0
         figsize = (total_width, total_height)
         fig, axes = plt.subplots(n_rows_fig, n_columns_fig, figsize=figsize)
         axes = np.array(axes).reshape(n_rows_fig, n_columns_fig)
 
-        # Set column titles
+        self._annotate_column_titles(axes, on_metrics)
+
+        return fig, axes, n_rows_fig, n_columns_fig, total_height
+
+    def _annotate_column_titles(self, axes, on_metrics):
         pad = 10
         col_titles = ["Metric: <{}>".format(metric) for metric in on_metrics]
         for ax, col_title in zip(axes[0], col_titles):
@@ -753,52 +716,152 @@ def performance_lineplot(
                 textcoords="offset points",
                 ha="center",
                 va="baseline",
-                **self.row_col_title_style,
+                **self._title_styles['row_col'],
             )
 
-        # Plot performances on each ax
-        vis_df_grp = vis_perf_dataframe.groupby(by=split_by + ["metric"])
-        for (key, _), ax in zip(vis_df_grp.groups.items(), axes.flatten()):
-            metric_name = key if len(split_by) == 0 else key[-1]
-            # Use seaborn.lineplot for visualization
-            kwargs = {
-                "data": vis_df_grp.get_group(key).reset_index(drop=True),
-                "x": x_column,
-                "y": "score",
-                "hue": "method",
-                "style": "dataset",
-                "ax": ax,
-            }
-            kwargs.update(lineplot_kwargs_)
-            ax = sns.lineplot(**kwargs)
-            # Set legend, border, x_label, y_label, and grid properties
-            ax.legend(
-                columnspacing=0.3,
-                borderaxespad=0.3,
-                handletextpad=0.3,
-                labelspacing=0.1,
-                handlelength=None,
-                borderpad=0.4,
+    @staticmethod
+    def _set_figure_suptitle(
+        fig, sup_title_, singular_text, plural_text, is_single, suptitle_style
+    ):
+        if sup_title_ is True:
+            fig.suptitle(
+                singular_text if is_single else plural_text, **suptitle_style
             )
-            ax = set_ax_border(ax, border_color="black", border_width=2)
-            ax.set_xlabel(x_label, **self.axis_title_style)
-            ax.set_ylabel(f"{metric_name}", **self.axis_title_style)
-            ax.grid(color="black", linestyle="-.", alpha=0.3)
+        elif isinstance(sup_title_, str):
+            fig.suptitle(sup_title_, **suptitle_style)
 
-        # Use tight layout
+    def _plot_single_performance_curve(
+        self, ax, data, x_column, x_label, metric_name, lineplot_kwargs
+    ):
+        kwargs = {
+            "data": data,
+            "x": x_column,
+            "y": "score",
+            "hue": "method",
+            "style": "dataset",
+            "ax": ax,
+        }
+        kwargs.update(lineplot_kwargs)
+        ax = sns.lineplot(**kwargs)
+        ax.legend(
+            columnspacing=0.3,
+            borderaxespad=0.3,
+            handletextpad=0.3,
+            labelspacing=0.1,
+            handlelength=None,
+            borderpad=0.4,
+        )
+        ax = set_ax_border(ax, border_color="black", border_width=2)
+        ax.set_xlabel(x_label, **self._title_styles['axis'])
+        ax.set_ylabel(f"{metric_name}", **self._title_styles['axis'])
+        ax.grid(color="black", linestyle="-.", alpha=0.3)
+
+    def _finalize_performance_figure(
+        self, fig, total_height, sup_title_, n_rows_fig, n_columns_fig
+    ):
         height_rect = (total_height - RESERVED_SUPTITLE_INCHES) / total_height
         plt.tight_layout(rect=(0, 0, 1, height_rect))
 
-        # Set super title
-        if sup_title_ == True:
-            if n_rows_fig == n_columns_fig == 1:
-                fig.suptitle("Performance Curve", **self.suptitle_style)
-            else:
-                fig.suptitle("Performance Curves", **self.suptitle_style)
-        elif sup_title_ == False:
-            pass
-        else:
-            fig.suptitle(sup_title_, **self.suptitle_style)
+        ImbalancedEnsembleVisualizer._set_figure_suptitle(
+            fig, sup_title_, "Performance Curve", "Performance Curves",
+            n_rows_fig == n_columns_fig == 1, self._title_styles['suptitle']
+        )
+
+    def performance_lineplot(
+        self,
+        on_ensembles: list = None,
+        on_datasets: list = None,
+        on_metrics: list = None,
+        split_by: list = [],
+        n_samples_as_x_axis: bool = False,
+        sub_figsize: tuple = (4.0, 3.3),
+        sup_title: bool or str = True,
+        **lineplot_kwargs,
+    ):
+        """Draw a performance line plot.
+
+        Parameters
+        ----------
+        on_ensembles : list of strings, default=None
+            The names of ensembles to include in the plot. It should be a
+            subset of ``self._fit_data.ensembles_.keys()``. if ``None``, all ensembles
+            fitted by the visualizer will be included.
+
+        on_datasets : list of strings, default=None
+            The names of evaluation datasets to include in the plot. It
+            should be a subset of ``self.eval_datasets_.keys()``. if ``None``,
+            all evaluation datasets will be included.
+
+        on_metrics : list of strings, default=None
+            The names of evaluation metrics to include in the plot. It
+            should be a subset of ``self.eval_metrics_.keys()``. if ``None``,
+            all evaluation metrics will be included.
+
+        split_by : list of {'method', 'dataset'}, default=[]
+            How to group the results for visualization.
+
+            - if contains ``'method'``, the performance results of different
+              ensemble methods will be displayed in independent sub-figures.
+            - if contains ``'dataset'``, the performance results on different
+              evaluation datasets will be displayed in independent sub-figures.
+
+        n_samples_as_x_axis : bool, default=False
+            Whether to use the number of training samples as the x-axis.
+
+        sub_figsize: (float, float), default=(4.0, 3.3)
+            The size of an subfigure (width, height in inches).
+            The overall figure size will be automatically determined by
+            (sub_figsize[0] * num_columns, sub_figsize[1] * num_rows).
+
+        sup_title: bool or str, default=True
+            The super title of the figure.
+
+            - if ``True``, automatically determines the super title.
+            - if ``False``, no super title will be displayed.
+            - if ``string``, super title will be ``sup_title``.
+
+        **lineplot_kwargs : key, value mappings
+            Other keyword arguments are passed down to
+            :meth:`seaborn.lineplot`.
+
+        Returns
+        -------
+        self : object
+        """
+
+        (vis_perf_dataframe, on_ensembles, on_metrics, n_datasets,
+         n_metrics, n_samples_as_x_axis_, sub_fig_width,
+         sub_fig_height, sup_title_, lineplot_kwargs_) = (
+            self._validate_performance_lineplot_params(
+                on_ensembles, on_datasets, on_metrics, split_by,
+                n_samples_as_x_axis, sub_figsize, sup_title, lineplot_kwargs
+            )
+        )
+
+        x_column, x_label, vis_perf_dataframe, n_ensembles = (
+            self._configure_x_axis(
+                vis_perf_dataframe, on_ensembles, n_samples_as_x_axis_
+            )
+        )
+
+        fig, axes, n_rows_fig, n_columns_fig, total_height = (
+            self._setup_figure_layout(
+                n_ensembles, n_datasets, n_metrics, split_by,
+                sub_fig_width, sub_fig_height, sup_title_, on_metrics
+            )
+        )
+
+        vis_df_grp = vis_perf_dataframe.groupby(by=split_by + ["metric"])
+        for (key, _), ax in zip(vis_df_grp.groups.items(), axes.flatten()):
+            metric_name = key if len(split_by) == 0 else key[-1]
+            data = vis_df_grp.get_group(key).reset_index(drop=True)
+            self._plot_single_performance_curve(
+                ax, data, x_column, x_label, metric_name, lineplot_kwargs_
+            )
+
+        self._finalize_performance_figure(
+            fig, total_height, sup_title_, n_rows_fig, n_columns_fig
+        )
 
         return fig, axes
 
@@ -809,7 +872,7 @@ def _collect_all_ensemble_confusion_matrix(
 
         print("Visualizer computing confusion matrices", end="")
         conf_matrices = {}
-        for ensemble_name, ensemble in self.ensembles_.items():
+        for ensemble_name, ensemble in self._fit_data.ensembles_.items():
             conf_matrices_ensemble = {}
             # Check whether the ensemble has verbose attribute
             has_verbose = hasattr(ensemble, "verbose")
@@ -854,6 +917,57 @@ def _check_is_subset(param: list, param_name: str, universal_set: list):
                 f" The possible values are {universal_set}."
             )
 
+    def _set_heatmap_titles(self, axes, on_datasets, on_ensembles):
+        pad = 10
+        row_titles = ["On dataset: <{}>".format(col) for col in on_datasets]
+        col_titles = ["Method: <{}>".format(row) for row in on_ensembles]
+        for ax, col_title in zip(axes[0], col_titles):
+            ax.annotate(
+                col_title,
+                xy=(0.5, 1),
+                xytext=(0, pad),
+                xycoords="axes fraction",
+                textcoords="offset points",
+                ha="center",
+                va="baseline",
+                **self._title_styles['row_col'],
+            )
+        for ax, row_title in zip(axes[:, 0], row_titles):
+            ax.annotate(
+                row_title,
+                xy=(0, 0.5),
+                xytext=(-ax.yaxis.labelpad - pad, 0),
+                xycoords=ax.yaxis.label,
+                textcoords="offset points",
+                ha="right",
+                va="center",
+                rotation=90,
+                **self._title_styles['row_col'],
+            )
+
+    def _plot_heatmap_subplots(
+        self, axes, on_datasets, on_ensembles, heatmap_kwargs, false_pred_only
+    ):
+        conf_matrices = self._fit_data.conf_matrices_
+        n_rows_fig, n_columns_fig = axes.shape
+        for dataset_name, i_row in zip(on_datasets, range(n_rows_fig)):
+            for ensemble_name, i_col in zip(on_ensembles, range(n_columns_fig)):
+                ax = axes[i_row, i_col]
+                conf_matrix_df = conf_matrices[ensemble_name][dataset_name]
+                kwargs = {
+                    "data": conf_matrix_df,
+                    "annot": True,
+                    "fmt": "d",
+                    "linewidths": 0.5,
+                    "ax": ax,
+                }
+                kwargs.update(heatmap_kwargs)
+                if false_pred_only:
+                    kwargs["mask"] = np.identity(conf_matrix_df.shape[0])
+                sns.heatmap(**kwargs)
+                ax.set_xlabel("Predicted Label", **self._title_styles['axis'])
+                ax.set_ylabel("Ground Truth", **self._title_styles['axis'])
+
     def confusion_matrix_heatmap(
         self,
         on_ensembles: list = None,
@@ -869,7 +983,7 @@ def confusion_matrix_heatmap(
         ----------
         on_ensembles : list of strings, default=None
             The names of ensembles to include in the plot. It should be a
-            subset of ``self.ensembles_.keys()``. if ``None``, all ensembles
+            subset of ``self._fit_data.ensembles_.keys()``. if ``None``, all ensembles
             fitted by the visualizer will be included.
 
         on_datasets : list of strings, default=None
@@ -911,10 +1025,10 @@ def confusion_matrix_heatmap(
             )
 
         on_ensembles = self._check_is_subset(
-            on_ensembles, "on_ensembles", self.vis_format_["ensemble_names"]
+            on_ensembles, "on_ensembles", self._fit_data.vis_format_["ensemble_names"]
         )
         on_datasets = self._check_is_subset(
-            on_datasets, "on_datasets", self.vis_format_["dataset_names"]
+            on_datasets, "on_datasets", self._fit_data.vis_format_["dataset_names"]
         )
         n_ensembles, n_datasets = len(on_ensembles), len(on_datasets)
 
@@ -945,73 +1059,18 @@ def confusion_matrix_heatmap(
         fig, axes = plt.subplots(n_rows_fig, n_columns_fig, figsize=figsize)
         axes = np.array(axes).reshape(n_rows_fig, n_columns_fig)
 
-        # Set titles for each column and row
-        pad = 10
-        row_titles = ["On dataset: <{}>".format(col) for col in on_datasets]
-        col_titles = ["Method: <{}>".format(row) for row in on_ensembles]
-        # Set column titles
-        for ax, col_title in zip(axes[0], col_titles):
-            ax.annotate(
-                col_title,
-                xy=(0.5, 1),
-                xytext=(0, pad),
-                xycoords="axes fraction",
-                textcoords="offset points",
-                ha="center",
-                va="baseline",
-                **self.row_col_title_style,
-            )
-
-        # Set row titles
-        for ax, row_title in zip(axes[:, 0], row_titles):
-            ax.annotate(
-                row_title,
-                xy=(0, 0.5),
-                xytext=(-ax.yaxis.labelpad - pad, 0),
-                xycoords=ax.yaxis.label,
-                textcoords="offset points",
-                ha="right",
-                va="center",
-                rotation=90,
-                **self.row_col_title_style,
-            )
-
-        # Plot confusion matrix heatmap on each ax
-        conf_matrices = self.conf_matrices_
-        for dataset_name, i_row in zip(on_datasets, range(n_rows_fig)):
-            for ensemble_name, i_col in zip(on_ensembles, range(n_columns_fig)):
-                ax = axes[i_row, i_col]
-                conf_matrix_df = conf_matrices[ensemble_name][dataset_name]
-                # Use seaborn.heatmap for visualization
-                kwargs = {
-                    "data": conf_matrix_df,
-                    "annot": True,
-                    "fmt": "d",
-                    "linewidths": 0.5,
-                    "ax": ax,
-                }
-                kwargs.update(heatmap_kwargs_)
-                # if false_pred_only, set mask to the heatmap
-                if false_pred_only:
-                    kwargs["mask"] = np.identity(conf_matrix_df.shape[0])
-                ax = sns.heatmap(**kwargs)
-                # Set x_label and y_label properties
-                ax.set_xlabel("Predicted Label", **self.axis_title_style)
-                ax.set_ylabel("Ground Truth", **self.axis_title_style)
+        self._set_heatmap_titles(axes, on_datasets, on_ensembles)
+        self._plot_heatmap_subplots(
+            axes, on_datasets, on_ensembles, heatmap_kwargs_, false_pred_only
+        )
 
         # Use tight layout
         height_rect = (total_height - RESERVED_SUPTITLE_INCHES * 1.3) / total_height
         plt.tight_layout(rect=(0, 0, 1, height_rect))
 
-        # Set super title
-        if sup_title_ == True:
-            if n_rows_fig == n_columns_fig == 1:
-                fig.suptitle("Confusion Matrix", **self.suptitle_style)
-            else:
-                fig.suptitle("Confusion Matrices", **self.suptitle_style)
-        elif sup_title_ == False:
-            pass
-        else:
-            fig.suptitle(sup_title_, **self.suptitle_style)
+        ImbalancedEnsembleVisualizer._set_figure_suptitle(
+            fig, sup_title_, "Confusion Matrix", "Confusion Matrices",
+            n_rows_fig == n_columns_fig == 1, self._title_styles['suptitle']
+        )
 
         return fig, axes