diff --git a/examples/classification_smoke.py b/examples/classification_smoke.py
new file mode 100644
index 00000000..66b5db6e
--- /dev/null
+++ b/examples/classification_smoke.py
@@ -0,0 +1,66 @@
+from ucimlrepo import fetch_ucirepo 
+import sklearn
+import numpy as np
+import pandas as pd
+
+from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder
+from sklearn.model_selection import cross_validate, train_test_split
+from sklearn.linear_model import Ridge, LogisticRegression
+from sklearn.utils import shuffle
+
+import suprb
+from sklearn.compose import make_column_transformer
+from suprb.utils import check_random_state
+from suprb.optimizer.rule.es import ES1xLambda
+from suprb.optimizer.rule.acceptance import MaxError
+from suprb.optimizer.solution.ga import GeneticAlgorithm
+from suprb.wrapper import SupRBWrapper
+
+from utils import log_scores
+
+
+if __name__ == '__main__':
+    random_state = 125
+    local_model = LogisticRegression(penalty='l1', C=0.1, random_state=random_state, solver='saga', tol=0.001, max_iter=1000)
+    #CLASSIFICATION
+    # fetch dataset 
+    iris = fetch_ucirepo(id=53)
+    X = iris.data.features.to_numpy()
+    y = iris.data.targets.to_numpy()
+    X, y = shuffle(X, y, random_state=random_state)
+    unique = np.unique(y)
+    toNum = dict(zip(unique, range(1, len(unique)+1)))
+    # Conversion of tragets to int required for mixing
+    # Similiar to sklearn.preprocessing.OrdinalEncoder
+    y = [toNum[x[0]] for x in y]
+    X = MinMaxScaler(feature_range=(-1, 1)).fit_transform(X)
+
+    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=random_state)
+
+    # Comparable with examples/example_2.py
+    model = SupRBWrapper(print_config=True,
+                        #n_iter=10,
+                         ## RULE GENERATION ##
+                         rule_generation=ES1xLambda(),
+                         rule_generation__n_iter=16,
+                         rule_generation__lmbda=16,
+                         rule_generation__operator='+',
+                         rule_generation__delay=10,
+                         rule_generation__random_state=random_state,
+                         rule_generation__n_jobs=4,
+                         rule_generation__init__model=local_model, 
+
+                         ## SOLUTION COMPOSITION ##
+                         solution_composition=GeneticAlgorithm(),
+                         solution_composition__init__mixing=suprb.solution.mixing_model.ErrorExperienceClassification(),
+                         solution_composition__n_iter=32,
+                         solution_composition__population_size=32,
+                         solution_composition__elitist_ratio=0.2,
+                         solution_composition__random_state=random_state,
+                         solution_composition__n_jobs=4)
+
+    scores = cross_validate(model, X, y, cv=4, n_jobs=4, verbose=10,
+                            scoring=['accuracy'],
+                            return_estimator=True, fit_params={'cleanup': True})
+
+    log_scores(scores)
diff --git a/examples/wrapper_example_1.py b/examples/wrapper_example_1.py
index 3a91854b..5c84dc4a 100644
--- a/examples/wrapper_example_1.py
+++ b/examples/wrapper_example_1.py
@@ -3,6 +3,7 @@
 
 from sklearn.preprocessing import StandardScaler, MinMaxScaler
 from sklearn.model_selection import cross_validate, train_test_split
+from sklearn.linear_model import Ridge
 
 
 from suprb.utils import check_random_state
@@ -49,16 +50,17 @@ def load_higdon_gramacy_lee(n_samples=1000, noise=0, random_state=None):
                          rule_generation__delay=150,
                          rule_generation__random_state=random_state,
                          rule_generation__n_jobs=1,
+                         rule_generation__init__model=Ridge(), 
 
                          ## SOLUTION COMPOSITION ##
                          solution_composition=GeneticAlgorithm(),
-                         solution_composition__n_iter=32,
+                         solution_composition__n_iter=10,
                          solution_composition__population_size=32,
                          solution_composition__elitist_ratio=0.2,
                          solution_composition__random_state=random_state,
                          solution_composition__n_jobs=1)
 
-    scores = cross_validate(model, X, y, cv=4, n_jobs=1, verbose=10,
+    scores = cross_validate(model, X, y, cv=4, n_jobs=4, verbose=10,
                             scoring=['r2', 'neg_mean_squared_error'],
                             return_estimator=True, fit_params={'cleanup': True})
 
diff --git a/requirements.txt b/requirements.txt
index a50b3f44..153e8220 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,3 +10,4 @@ joblib~=1.1.0
 tqdm~=4.62.3
 pytest~=6.2.5
 protobuf~=3.20.0
+ucimlrepo
\ No newline at end of file
diff --git a/suprb/base.py b/suprb/base.py
index 92118214..3ef60d7c 100644
--- a/suprb/base.py
+++ b/suprb/base.py
@@ -1,9 +1,9 @@
 from __future__ import annotations
 
-from abc import abstractmethod, ABCMeta
+from abc import abstractmethod, ABCMeta, ABC
 
 import numpy as np
-from sklearn.base import BaseEstimator, RegressorMixin
+from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin
 
 
 class BaseComponent(BaseEstimator, metaclass=ABCMeta):
@@ -47,13 +47,13 @@ def _more_str_attributes(self) -> dict:
         return {}
 
 
-class BaseRegressor(BaseEstimator, RegressorMixin, metaclass=ABCMeta):
-    """A base (composite) Regressor."""
+class BaseSupervised(BaseEstimator, metaclass=ABCMeta):
+    """A base (composite) Estimator for supervised learning."""
 
     is_fitted_: bool
 
     @abstractmethod
-    def fit(self, X: np.ndarray, y: np.ndarray) -> BaseRegressor:
+    def fit(self, X: np.ndarray, y: np.ndarray) -> BaseSupervised:
         """ A reference implementation of a fitting function.
 
         Parameters
@@ -68,7 +68,6 @@ def fit(self, X: np.ndarray, y: np.ndarray) -> BaseRegressor:
         self : BaseEstimator
             Returns self.
         """
-
         pass
 
     @abstractmethod
@@ -85,5 +84,43 @@ def predict(self, X: np.ndarray):
         y : np.ndarray
             Returns the estimation with shape (n_samples,).
         """
+        pass
+
+    @abstractmethod
+    def score(self, X, y, sample_weight=None):
+        pass
+
+class BaseRegressor(BaseSupervised, RegressorMixin, metaclass=ABCMeta):
+    """A base (composite) Regressor."""
+
+    is_fitted_: bool
+
+    @abstractmethod
+    def fit(self, X: np.ndarray, y: np.ndarray) -> BaseRegressor:
+        pass
+
+    @abstractmethod
+    def predict(self, X: np.ndarray):
+        pass
+
+
+class BaseClassifier(BaseSupervised, ClassifierMixin, metaclass=ABCMeta):
+    """A base (composite) Classifier."""
+
+    is_fitted_: bool
+
+    @abstractmethod
+    def fit(self, X: np.ndarray, y: np.ndarray) -> BaseClassifier:
+        pass
+
+    @abstractmethod
+    def predict(self, X: np.ndarray):
+        pass
 
+class SupervisedMixin(ABC):
+    @abstractmethod
+    def score(self, X, y, sample_weight=None):
         pass
+
+SupervisedMixin.register(RegressorMixin)
+SupervisedMixin.register(ClassifierMixin)
\ No newline at end of file
diff --git a/suprb/fitness.py b/suprb/fitness.py
index ae3d6948..d15fc56f 100644
--- a/suprb/fitness.py
+++ b/suprb/fitness.py
@@ -20,6 +20,14 @@ def pseudo_accuracy(error: float, beta=2) -> float:
     assert beta > 0
     return np.exp(-beta * error)
 
+def pseudo_error(accuracy: float) -> float:
+    return 1 - accuracy
+
+def actual_accuracy(error: float) -> float:
+    '''
+    for classification
+    '''
+    return 1 - error
 
 def emary(alpha: float, x1: float, x2: float) -> float:
     """
diff --git a/suprb/json.py b/suprb/json.py
index e0125d3a..a8b94fe0 100644
--- a/suprb/json.py
+++ b/suprb/json.py
@@ -86,7 +86,11 @@ def _save_config(suprb, json_config):
         elif isinstance(value, primitive):
             json_config["config"][key] = value
         else:
-            json_config["config"][key] = _get_full_class_name(value)
+            cname = _get_full_class_name(value)
+            if cname == "NoneType":
+                json_config["config"][key] = None
+            else:    
+                json_config["config"][key] = cname
 
 
 def _save_pool(pool, json_config):
@@ -113,6 +117,7 @@ def _convert_rule_to_json(rule):
     return {"error_":        rule.error_,
             "experience_":   rule.experience_,
             "match":         _convert_dict_to_json(vars(rule.match)),
+            "matching_type": _get_full_class_name(rule.match),
             "is_fitted_":    rule.is_fitted_,
             "model":         {"coef_":          _convert_to_json_format(getattr(rule.model, "coef_")),
                               "intercept_":     getattr(rule.model, "intercept_")}}
@@ -194,7 +199,7 @@ def _load_pool(json_dict, suprb):
 
 def _convert_json_to_rule(json_rule, json_dict):
 
-    rule = Rule(_convert_matching_type(json_rule["match"], json_dict["config"]["matching_type"]),
+    rule = Rule(_convert_matching_type(json_rule["match"], json_rule["matching_type"]),
                 _convert_from_json_to_array(json_dict["input_space"]),
                 _convert_model(json_rule["model"], json_dict["config"]["rule_generation__init__model"]),
                 _get_class(json_dict["config"]["rule_generation__init__fitness"]))
diff --git a/suprb/logging/base.py b/suprb/logging/base.py
index 4a8d5116..267e304a 100644
--- a/suprb/logging/base.py
+++ b/suprb/logging/base.py
@@ -2,28 +2,28 @@
 
 import numpy as np
 
-from suprb.base import BaseComponent, BaseRegressor
+from suprb.base import BaseComponent, BaseSupervised
 
 
 class BaseLogger(BaseComponent):
     """The base class for loggers."""
 
     @abstractmethod
-    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         """Logs initial parameters, before any iteration or fitting has taken place."""
         pass
 
     @abstractmethod
-    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor, iteration: int):
+    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised, iteration: int):
         """Logs an iteration of the estimator. May actually compute errors or scores, depending on the state."""
         pass
 
     @abstractmethod
-    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         """Log the final state of the estimator. It is assumed that the fitting process is already completed."""
         pass
 
     @abstractmethod
-    def get_elitist(self, estimator: BaseRegressor):
+    def get_elitist(self, estimator: BaseSupervised):
         """Log the final elitist"""
         pass
diff --git a/suprb/logging/combination.py b/suprb/logging/combination.py
index 7e3ed934..d6f1ab8a 100644
--- a/suprb/logging/combination.py
+++ b/suprb/logging/combination.py
@@ -3,7 +3,7 @@
 import numpy as np
 from sklearn import clone
 
-from suprb.base import BaseRegressor
+from suprb.base import BaseSupervised
 from . import BaseLogger
 
 
@@ -16,7 +16,7 @@ def __init__(self, loggers: list[tuple[str, BaseLogger]]):
         """An unique name for every logger must be supplied, such that the parameter get/set are well-defined."""
         self.loggers = loggers
 
-    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         if any(map(lambda logger: isinstance(logger, CombinedLogger), self.loggers)):
             warnings.warn("Nesting loggers is not recommended. Please add all loggers to this top-level logger.")
 
@@ -25,11 +25,11 @@ def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
         for _, logger in self.loggers_:
             logger.log_init(X=X, y=y, estimator=estimator)
 
-    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor, iteration: int):
+    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised, iteration: int):
         for _, logger in self.loggers_:
             logger.log_iteration(X=X, y=y, estimator=estimator, iteration=iteration)
 
-    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         for _, logger in self.loggers_:
             logger.log_final(X=X, y=y, estimator=estimator)
 
@@ -71,5 +71,5 @@ def _replace_logger(self, name: str, new_val: BaseLogger):
                 break
         self.loggers = new_loggers
 
-    def get_elitist(self, estimator: BaseRegressor):
+    def get_elitist(self, estimator: BaseSupervised):
         pass
diff --git a/suprb/logging/default.py b/suprb/logging/default.py
index 41d26c0e..d414231c 100644
--- a/suprb/logging/default.py
+++ b/suprb/logging/default.py
@@ -5,8 +5,8 @@
 
 from . import BaseLogger
 from .metrics import matched_training_samples, genome_diversity
-from .. import json as suprb_json
-from suprb.base import BaseRegressor
+from .. import json
+from suprb.base import BaseSupervised
 
 
 class DefaultLogger(BaseLogger):
@@ -25,13 +25,13 @@ def log_params(self, **kwargs):
         for key, value in kwargs.items():
             self.log_param(key=key, value=value)
 
-    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         self.params_ = {}
         self.metrics_ = defaultdict(dict)
 
         self.log_params(**estimator.get_params())
 
-    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor, iteration: int):
+    def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised, iteration: int):
         def log_metric(key, value):
             self.log_metric(key=key, value=value, step=estimator.step_)
 
@@ -61,17 +61,19 @@ def log_metric_min_max_mean(metric_name: str, attribute_name: str, lst: list):
         elitist = estimator.solution_composition_.elitist()
         log_metric("elitist_fitness", elitist.fitness_)
         log_metric("elitist_error", elitist.error_)
+        if elitist.isClassifier:
+            log_metric("elitist_accuracy", elitist.score(X, y))
         log_metric("elitist_complexity", elitist.complexity_)
         log_metric("elitist_matched", matched_training_samples(elitist.subpopulation))
-        # log_metric("elitist_rules", elitist.pool)
+        log_metric("elitist_rules", elitist.pool)
 
         # Log performance
         log_metric("training_score", elitist.score(X, y))
 
-    def get_elitist(self, estimator: BaseRegressor):
+    def get_elitist(self, estimator: BaseSupervised):
         json_data = {}
-        suprb_json._save_pool(estimator.solution_composition_.elitist().pool, json_data)
+        json._save_pool(estimator.solution_composition_.elitist().pool, json_data)
         return json_data
 
-    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor):
+    def log_final(self, X: np.ndarray, y: np.ndarray, estimator: BaseSupervised):
         pass
diff --git a/suprb/optimizer/rule/acceptance.py b/suprb/optimizer/rule/acceptance.py
index 86e078fc..222ba463 100644
--- a/suprb/optimizer/rule/acceptance.py
+++ b/suprb/optimizer/rule/acceptance.py
@@ -2,8 +2,11 @@
 
 import numpy as np
 
+from sklearn.metrics import precision_score, recall_score, f1_score
+
 from suprb.base import BaseComponent
 from suprb.rule import Rule
+from suprb.fitness import pseudo_error
 
 
 class RuleAcceptance(BaseComponent, metaclass=ABCMeta):
@@ -38,4 +41,34 @@ def __call__(self, rule: Rule, X: np.ndarray, y: np.ndarray) -> bool:
             return False
         local_y = y[rule.match_set_]
         default_error = np.sum(local_y ** 2) / (len(local_y) * self.beta)
+        if rule.isClassifier:
+            # default error is the trivial solution of always choosing the most common label
+            local_y = [round(y) for y in local_y]
+            default_accuracy = np.bincount(local_y).max() / (len(local_y) * self.beta)
+            default_error = pseudo_error(default_accuracy)
         return rule.error_ <= default_error
+
+
+class Precision(RuleAcceptance):
+    """Insert if the rule has a precision greater or equal to a threshold"""
+
+    def __init__(self, min_precission: float = 0.9):
+        self.min_precission = min_precission
+
+    def __call__(self, rule: Rule, X: np.ndarray, y: np.ndarray) -> bool:
+        if rule.experience_ < 1:
+            return False
+        local_y = y[rule.match_set_]
+        return precision_score(local_y, rule.predict(X[rule.match_set_]), average='macro', zero_division=np.nan) >= self.min_precission
+    
+class F1_Score(RuleAcceptance):
+    """Insert if the rule has a f1_score greater or equal to a threshold"""
+
+    def __init__(self, min_f1: float = 0.9):
+        self.min_precission = min_f1
+
+    def __call__(self, rule: Rule, X: np.ndarray, y: np.ndarray) -> bool:
+        if rule.experience_ < 1:
+            return False
+        local_y = y[rule.match_set_]
+        return f1_score(local_y, rule.predict(X[rule.match_set_])) >= self.min_f1
\ No newline at end of file
diff --git a/suprb/rule/base.py b/suprb/rule/base.py
index 252df0d1..3b4e5344 100644
--- a/suprb/rule/base.py
+++ b/suprb/rule/base.py
@@ -4,11 +4,11 @@
 from typing import Union
 
 import numpy as np
-from sklearn.base import RegressorMixin, clone
-from sklearn.metrics import mean_squared_error
+from sklearn.base import RegressorMixin, ClassifierMixin, clone
+from sklearn.metrics import mean_squared_error, accuracy_score
 
 from suprb.base import SolutionBase
-from suprb.fitness import BaseFitness
+from suprb.fitness import BaseFitness, pseudo_error
 from .matching import MatchingFunction
 
 
@@ -38,11 +38,15 @@ class Rule(SolutionBase):
     match_set_: np.ndarray
     pred_: Union[np.ndarray, None]  # only the prediction of matching points, so of x[match_]
 
-    def __init__(self, match: MatchingFunction, input_space: np.ndarray, model: RegressorMixin, fitness: RuleFitness):
+    def __init__(self, match: MatchingFunction, input_space: np.ndarray, model: RegressorMixin | ClassifierMixin, fitness: RuleFitness):
         self.match = match
         self.input_space = input_space
         self.model = model
         self.fitness = fitness
+        if isinstance(model, ClassifierMixin):
+            self.isClassifier = True
+        else:
+            self.isClassifier = False
 
     def fit(self, X: np.ndarray, y: np.ndarray) -> Rule:
 
@@ -70,11 +74,24 @@ def fit(self, X: np.ndarray, y: np.ndarray) -> Rule:
         # Get all data points which match the bounds.
         X, y = X[self.match_set_], y[self.match_set_]
 
+        # No reason to fit if only one label in match_set
+        if self.isClassifier:
+            if len(np.unique(y)) == 1:
+                self.pred_ = y[0]
+                self.error_ = 1
+                self.fitness_ = self.fitness(self)
+                self.experience_ = float(X.shape[0])
+                self.is_fitted_ = False
+                return self
+
         # Create and fit the model
         self.model.fit(X, y)
 
         self.pred_ = self.model.predict(X)
-        self.error_ = max(mean_squared_error(y, self.pred_), 1e-4)  # TODO: make min a parameter?
+        if not self.isClassifier:
+            self.error_ = max(mean_squared_error(y, self.pred_), 1e-4)  # TODO: make min a parameter?
+        elif self.isClassifier:
+            self.error_ = max(pseudo_error(accuracy_score(y, self.pred_)), 1e-4)
         self.fitness_ = self.fitness(self)
         self.experience_ = float(X.shape[0])
 
@@ -86,6 +103,9 @@ def volume_(self):
         return self.match.volume_
 
     def predict(self, X: np.ndarray):
+        # Make array 2dim in case of single sample
+        if np.ndim(X) == 1:
+            X = X.reshape(1, -1)
         return self.model.predict(X)
 
     def clone(self, **kwargs) -> Rule:
diff --git a/suprb/rule/fitness.py b/suprb/rule/fitness.py
index 59a2e724..93e14c8c 100644
--- a/suprb/rule/fitness.py
+++ b/suprb/rule/fitness.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from suprb.fitness import pseudo_accuracy, emary, wu
+from suprb.fitness import pseudo_accuracy, actual_accuracy, emary, wu
 from . import Rule, RuleFitness
 
 
@@ -16,6 +16,8 @@ class PseudoAccuracy(RuleFitness):
     def __call__(self, rule: Rule) -> float:
         # note that we multiply solely for readability reasons without
         # any performance impact
+        if rule.isClassifier:
+            return  actual_accuracy(rule.error_) * 100
         return pseudo_accuracy(rule.error_) * 100
 
 
@@ -29,10 +31,15 @@ def __init__(self, alpha: float = None):
 
     def __call__(self, rule: Rule) -> float:
         diff = rule.input_space[:, 1] - rule.input_space[:, 0]
+        diff[diff == 0] = 1.0  # avoid division by zero
         input_space_volume = np.prod(diff)
 
         volume_share = rule.volume_ / input_space_volume
-        return self.fitness_func_(alpha=self.alpha, x1=pseudo_accuracy(rule.error_, beta=2), x2=volume_share) * 100
+        if rule.isClassifier:
+            accuracy = actual_accuracy(rule.error_)
+        else:
+            accuracy = pseudo_accuracy(rule.error_, beta=2)
+        return self.fitness_func_(alpha=self.alpha, x1 = accuracy, x2=volume_share) * 100
 
 
 class VolumeEmary(VolumeRuleFitness):
diff --git a/suprb/rule/initialization.py b/suprb/rule/initialization.py
index a070053a..1cc763be 100644
--- a/suprb/rule/initialization.py
+++ b/suprb/rule/initialization.py
@@ -2,7 +2,7 @@
 from typing import Union
 import numpy as np
 from scipy.stats import halfnorm
-from sklearn.base import RegressorMixin, clone
+from sklearn.base import RegressorMixin, ClassifierMixin, clone
 from sklearn.linear_model import LinearRegression
 from sklearn.linear_model import Ridge
 
@@ -25,16 +25,19 @@ class RuleInit(BaseComponent, metaclass=ABCMeta):
     """
 
     def __init__(self, bounds: np.ndarray = None,
-                 model: RegressorMixin = None,
+                 model: BaseComponent = None,
                  fitness: RuleFitness = None,
                  matching_type: MatchingFunction = None):
         self.bounds = bounds
         self.model = model
         self.fitness = fitness
         self.matching_type = matching_type
-
         self._validate_components(model=Ridge(alpha=0.01),
                                   fitness=VolumeWu())
+        if isinstance(model, ClassifierMixin):
+            self.isClassifier = True
+        else:
+            self.isClassifier = False
 
     @property
     def matching_type(self):
diff --git a/suprb/rule/matching.py b/suprb/rule/matching.py
index 551e901d..abc51c8a 100644
--- a/suprb/rule/matching.py
+++ b/suprb/rule/matching.py
@@ -48,7 +48,7 @@ class OrderedBound(MatchingFunction):
     A standard interval-based matching function producing multi-dimensional
     hyperrectangular conditions. In effect, a lower (l) and upper bound (u) are
     specified for each dimension. Those bounds always fulfill l <= u
-    An example x is matched iff l_i <= x_i <= u_i for all dimensions i
+    An example x is matched if l_i <= x_i <= u_i for all dimensions i
     """
 
     def __init__(self, bounds: np.ndarray = np.array([])):
@@ -84,7 +84,7 @@ class UnorderedBound(MatchingFunction):
     A standard interval-based matching function producing multi-dimensional
     hyperrectangular conditions. Two bounds (p and q) exist which have
     no explicit ordering but are instead sorted during the matching process
-    An example x is matched iff q_i <= x_i <= p_i for all dimensions i
+    An example x is matched if q_i <= x_i <= p_i or p_i <= x_i <= q_i for all dimensions i
     """
 
     def __init__(self, bounds: np.ndarray):
@@ -173,7 +173,7 @@ class MinPercentage(MatchingFunction):
     A standard interval-based matching function producing multi-dimensional
     hyperrectangular conditions. In effect, a lower bound (l) and
     a distance proportion (p) are specified for each dimension.
-    An example x is matched iff
+    An example x is matched if
     l_i <= x_i <= l_i + p_i * (max_i - l_i) for all dimensions i
     """
 
diff --git a/suprb/solution/base.py b/suprb/solution/base.py
index 2078a43f..67278707 100644
--- a/suprb/solution/base.py
+++ b/suprb/solution/base.py
@@ -2,14 +2,15 @@
 
 import itertools
 from abc import ABCMeta, abstractmethod
+from typing import Union
 
 import numpy as np
-from sklearn.base import RegressorMixin
-from sklearn.metrics import mean_squared_error
+from sklearn.base import RegressorMixin, ClassifierMixin, BaseEstimator
+from sklearn.metrics import mean_squared_error, accuracy_score, r2_score
 
 from suprb.rule import Rule
-from suprb.base import BaseComponent, SolutionBase
-from suprb.fitness import BaseFitness
+from suprb.base import BaseComponent, SolutionBase, SupervisedMixin
+from suprb.fitness import BaseFitness, pseudo_error
 
 
 class MixingModel(BaseComponent, metaclass=ABCMeta):
@@ -33,7 +34,7 @@ def __call__(self, solution: Solution) -> float:
         pass
 
 
-class Solution(SolutionBase, RegressorMixin):
+class Solution(SolutionBase, SupervisedMixin):
     """Solution that mixes a subpopulation of rules."""
 
     input_size_: int
@@ -44,10 +45,38 @@ def __init__(self, genome: np.ndarray, pool: list[Rule], mixing: MixingModel, fi
         self.pool = pool
         self.mixing = mixing
         self.fitness = fitness
-
-    def fit(self, X: np.ndarray, y: np.ndarray) -> Solution:
-        pred = self.predict(X, cache=True)
-        self.error_ = max(mean_squared_error(y, pred), 1e-4)
+        self.isClassifier = None
+
+    def score(self, X, y, sample_weight=None):
+        if not self.pool:
+            return 0.0
+        if self.isClassifier is None:
+            if isinstance(self.pool[0].model, ClassifierMixin):
+                self.isClassifier = True
+            else:
+                self.isClassifier = False
+        if self.isClassifier:
+            return accuracy_score(y, self.predict(X), sample_weight=sample_weight)
+        else:
+            y_pred = self.predict(X)
+            return r2_score(y, y_pred, sample_weight=sample_weight)
+            
+
+    def fit(self, X: np.ndarray, y: np.ndarray, cache = True) -> Solution:
+        pred = self.predict(X, cache=cache)
+        if not self.pool:
+            self.error_ = 9999
+        else:
+            if self.isClassifier is None:
+                if isinstance(self.pool[0].model, ClassifierMixin):
+                    self.isClassifier = True
+                else:
+                    self.isClassifier = False
+            if self.isClassifier:
+                self.error_ = max(pseudo_error(accuracy_score(y, pred)), 1e-4)
+            else:
+                self.error_ = max(mean_squared_error(y, pred), 1e-4)
+            
         self.input_size_ = self.genome.shape[0]
         self.complexity_ = np.sum(self.genome).item()  # equivalent to np.count_nonzero, but possibly faster
         self.fitness_ = self.fitness(self)
diff --git a/suprb/solution/fitness.py b/suprb/solution/fitness.py
index 5fd4fbc9..7fbc2292 100644
--- a/suprb/solution/fitness.py
+++ b/suprb/solution/fitness.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from suprb.fitness import emary, pseudo_accuracy, wu
+from suprb.fitness import emary, pseudo_accuracy, wu, actual_accuracy
 from . import Solution, SolutionFitness
 
 
@@ -29,7 +29,8 @@ def __init__(self, alpha: float):
         self.alpha = alpha
 
     def __call__(self, solution: Solution) -> float:
-        return self.fitness_func_(self.alpha, pseudo_accuracy(solution.error_),
+        acc = actual_accuracy(solution.error_) if solution.isClassifier else pseudo_accuracy(solution.error_)
+        return self.fitness_func_(self.alpha, acc,
                                   c_norm(solution.complexity_, self.max_genome_length_)) * 100
 
 
diff --git a/suprb/solution/mixing_model.py b/suprb/solution/mixing_model.py
index 1c404632..ad0627fc 100644
--- a/suprb/solution/mixing_model.py
+++ b/suprb/solution/mixing_model.py
@@ -1,5 +1,7 @@
 import numpy as np
 
+from sklearn.base import RegressorMixin, ClassifierMixin
+
 from suprb.rule import Rule
 from suprb.utils import check_random_state, RandomState
 from . import MixingModel
@@ -67,6 +69,104 @@ def __call__(self, subpopulation: list[Rule], dim: int = None) -> list[Rule]:
         return np.clip(experiences, self.lower_bound * dim, self.upper_bound * dim)
 
 
+def get_local_pred(X: np.ndarray, subpopulation: list[Rule], cache: bool):
+        input_size = X.shape[0]
+        if subpopulation:
+            local_pred = np.zeros((len(subpopulation), input_size))
+
+        if cache:
+            # Use the precalculated matches and predictions from fit()
+            matches = [rule.match_set_ for rule in subpopulation]
+            for i, rule in enumerate(subpopulation):
+                local_pred[i][matches[i]] = rule.pred_
+        else:
+            # Generate all data new
+            matches = [rule.match(X) for rule in subpopulation]
+            for i, rule in enumerate(subpopulation):
+                if not matches[i].any():
+                    continue
+                local_pred[i][matches[i]] = rule.predict(X[matches[i]])
+
+        return local_pred, matches
+
+
+class MostVoted(MixingModel):
+
+    def __init__(self, filter_subpopulation: FilterSubpopulation = FilterSubpopulation()):
+        self.input_size = None
+        self.filter_subpopulation = filter_subpopulation
+
+    def __call__(self, X: np.ndarray, subpopulation: list[Rule], cache=False) -> np.ndarray:
+        self.input_size = X.shape[0]
+
+        # No need to perform any calculation if no rule was selected.
+        if not subpopulation:
+            return np.zeros(self.input_size)
+
+        subpopulation = self.filter_subpopulation(subpopulation)
+
+        local_pred, matches = get_local_pred(X, subpopulation, cache)
+        pred_per_sample = local_pred.transpose()
+        out = np.zeros(self.input_size)
+        for x, pred in enumerate(pred_per_sample):
+            # strip to ensure only valid predictions are counted
+            stripped_pred = [i for i in pred if i != 0]
+            if not stripped_pred:
+                out[x] = int(np.random.choice(local_pred.flatten()))
+                continue
+            out[x] = np.bincount(stripped_pred).argmax()
+        out = [int(label) for label in out]
+        return out
+
+
+class ErrorExperienceClassification(MixingModel):
+    """
+    Performs mixing similar to the Inverse Variance Heuristic from
+    https://researchportal.bath.ac.uk/en/studentTheses/learning-classifier-systems-from-first-principles-a-probabilistic,
+    but using (error / experience) as a mixing function.
+    """
+
+    def __init__(self, filter_subpopulation: FilterSubpopulation = FilterSubpopulation(),
+                 experience_calculation: ExperienceCalculation = ExperienceCalculation(),
+                 experience_weight: float = 1):
+        self.input_size = None
+        self.filter_subpopulation = filter_subpopulation
+        self.experience_calculation = experience_calculation
+        self.experience_weight = experience_weight
+
+    def __call__(self, X: np.ndarray, subpopulation: list[Rule], cache=False) -> np.ndarray:
+        self.input_size = X.shape[0]
+
+        # No need to perform any calculation if no rule was selected.
+        if not subpopulation:
+            return np.zeros(self.input_size)
+
+        subpopulation = self.filter_subpopulation(subpopulation)
+
+        local_pred, matches = get_local_pred(X, subpopulation, cache)
+        pred_per_sample = local_pred.transpose()
+        match_per_sample = np.array(matches).transpose()
+        taus = self._get_taus(subpopulation, X.shape[1])
+        out = np.zeros(self.input_size)
+        for x, pred in enumerate(pred_per_sample):
+            # strip to ensure only valid predictions are counted
+            stripped_pred = pred[match_per_sample[x]]
+            stripped_pred = [int(label) for label in stripped_pred]
+            if not stripped_pred:
+                out[x] = int(np.random.choice([i for i in local_pred.flatten() if i != 0]))
+                continue
+            local_taus = taus[match_per_sample[x]]
+            out[x] = np.bincount(stripped_pred, weights=local_taus).argmax()
+        return out
+
+    def _get_taus(self, subpopulation: list[Rule], dim: int):
+        # Get errors and experience of all rules in subpopulation
+        experiences = self.experience_calculation(subpopulation, dim)
+        errors = np.array([rule.error_ for rule in subpopulation])
+
+        return (1 / errors) * (experiences * self.experience_weight)
+
+    
 class ErrorExperienceHeuristic(MixingModel):
     """
     Performs mixing similar to the Inverse Variance Heuristic from
@@ -91,7 +191,7 @@ def __call__(self, X: np.ndarray, subpopulation: list[Rule], cache=False) -> np.
 
         subpopulation = self.filter_subpopulation(subpopulation)
 
-        local_pred, matches = self._get_local_pred(X, subpopulation, cache)
+        local_pred, matches = get_local_pred(X, subpopulation, cache)
         taus = self._get_taus(subpopulation, X.shape[1])
 
         # Stack all local predictions and sum them weighted with tau
@@ -100,26 +200,12 @@ def __call__(self, X: np.ndarray, subpopulation: list[Rule], cache=False) -> np.
         # Normalize
         tau_sum = self._get_tau_sum(subpopulation, matches, taus)
         out = pred / tau_sum
+        
+        # Round predicitions to int for Classification (intended for ordered labels only)
+        if isinstance(subpopulation[0].model, ClassifierMixin):
+            out = [round(x) for x in out]
         return out
 
-    def _get_local_pred(self, X: np.ndarray, subpopulation: list[Rule], cache: bool):
-        local_pred = np.zeros((len(subpopulation), self.input_size))
-
-        if cache:
-            # Use the precalculated matches and predictions from fit()
-            matches = [rule.match_set_ for rule in subpopulation]
-            for i, rule in enumerate(subpopulation):
-                local_pred[i][matches[i]] = rule.pred_
-        else:
-            # Generate all data new
-            matches = [rule.match(X) for rule in subpopulation]
-            for i, rule in enumerate(subpopulation):
-                if not matches[i].any():
-                    continue
-                local_pred[i][matches[i]] = rule.predict(X[matches[i]])
-
-        return local_pred, matches
-
     def _get_taus(self, subpopulation: list[Rule], dim: int):
         # Get errors and experience of all rules in subpopulation
         experiences = self.experience_calculation(subpopulation, dim)
diff --git a/suprb/suprb.py b/suprb/suprb.py
index ca4ebb7e..d54963c7 100644
--- a/suprb/suprb.py
+++ b/suprb/suprb.py
@@ -1,12 +1,15 @@
 import warnings
 import traceback
+import copy
 
 import numpy as np
 from sklearn import clone
 from sklearn.utils import check_X_y
 from sklearn.utils.validation import check_is_fitted, check_array
+from sklearn.base import RegressorMixin, ClassifierMixin, BaseEstimator
+from sklearn.metrics import mean_squared_error, accuracy_score, r2_score
 
-from .base import BaseRegressor
+from .base import BaseSupervised, BaseRegressor, BaseClassifier
 from .exceptions import PopulationEmptyWarning
 from .solution import Solution
 from .logging import BaseLogger
@@ -21,7 +24,7 @@
 from .solution.fitness import PseudoBIC
 
 
-class SupRB(BaseRegressor):
+class SupRB(BaseSupervised):
     """ The multi-solution batch learning LCS developed by the Organic Computing group at Universität Augsburg.
 
     Parameters
@@ -75,6 +78,7 @@ class SupRB(BaseRegressor):
     n_features_in_: int
 
     logger_: BaseLogger
+    is_fitted_: bool = False
 
     def __init__(self,
                  rule_generation: RuleGeneration = None,
@@ -102,6 +106,7 @@ def __init__(self,
         self.n_jobs = n_jobs
         self.early_stopping_patience = early_stopping_patience
         self.early_stopping_delta = early_stopping_delta
+        self.isClassifier = None
 
     def check_early_stopping(self):
         if self.early_stopping_patience > 0:
@@ -117,7 +122,23 @@ def check_early_stopping(self):
                     return True
         return False
 
-    def fit(self, X: np.ndarray, y: np.ndarray, cleanup=False):
+    
+    def score(self, X, y, sample_weight=None):
+        if not self.pool_:
+            return 0.0
+        if self.isClassifier is None:
+            if isinstance(self.pool_[0].model, ClassifierMixin):
+                self.isClassifier = True
+            else:
+                self.isClassifier = False
+        if self.isClassifier:
+            return accuracy_score(y, self.predict(X), sample_weight=sample_weight)
+        else:
+            y_pred = self.predict(X)
+            return r2_score(y, y_pred, sample_weight=sample_weight)
+        
+        
+    def fit(self, X: np.ndarray, y: np.ndarray, cleanup=False, force=False):
         """ Fit SupRB.2.
 
             Parameters
@@ -136,6 +157,9 @@ def fit(self, X: np.ndarray, y: np.ndarray, cleanup=False):
                 Returns self.
         """
 
+        if self.is_fitted_ and not force:
+            return self
+
         # Set these values so we gracefully exit on error
         self.early_stopping_counter_ = 0
         self.previous_fitness_ = 0
@@ -146,8 +170,12 @@ def fit(self, X: np.ndarray, y: np.ndarray, cleanup=False):
         self.elitist_.complexity_ = 99999
 
         # Check that x and y have correct shape
-        X, y = check_X_y(X, y, dtype='float64', y_numeric=True)
-        y = check_array(y, ensure_2d=False, dtype='float64')
+        if isinstance(self, BaseRegressor):
+            X, y = check_X_y(X, y, dtype='float64', y_numeric=True)
+            y = check_array(y, ensure_2d=False, dtype='float64')
+        else:
+            X, y = check_X_y(X, y, dtype=None, y_numeric=False)
+            y = check_array(y, ensure_2d=False, dtype=None)
 
         # Init sklearn interface
         self.n_features_in_ = X.shape[1]
@@ -271,7 +299,7 @@ def _compose_solution(self, X: np.ndarray, y: np.ndarray):
         self.solution_composition_.optimize(X, y)
 
     def predict(self, X: np.ndarray):
-        # Check is fit had been called
+        # Check if fit had been called
         check_is_fitted(self, ['is_fitted_'])
         # Input validation
         X = check_array(X)
@@ -339,3 +367,20 @@ def _more_tags(self):
         return {
             'poor_score': True,
         }
+
+    def model_swap(self, local_model):
+        solution = self.elitist_
+        pool = []
+        for old_rule in solution.pool:
+            rule = old_rule.clone(model = copy.deepcopy(local_model))
+            rule.isClassifier = old_rule.isClassifier
+            pool.append(rule)
+        solution.pool = pool
+        self.is_fitted_ = True
+        return solution
+    
+    def model_swap_fit(self, local_model, X, y):
+        solution = self.model_swap(local_model)
+        for rule in solution.pool:
+            rule.fit(X, y)
+        solution.fit(X, y, cache=False)
\ No newline at end of file