Repository: NiklasMelton/AdaptiveResonanceLib (language: python)

Issue:
Match tracking update



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pypi-publish.yml b/.github/workflows/pypi-publish.yml
new file mode 100644
index 0000000..b7a704b
--- /dev/null
+++ b/.github/workflows/pypi-publish.yml
@@ -0,0 +1,39 @@
+# This workflow will upload a Python Package using Twine when a release is created
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python#publishing-to-package-registries
+
+# This workflow uses actions that are not certified by GitHub.
+# They are provided by a third-party and are governed by
+# separate terms of service, privacy policy, and support
+# documentation.
+
+name: Upload Python Package
+
+on:
+  release:
+    types: [published]
+
+permissions:
+  contents: read
+
+jobs:
+  deploy:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python
+      uses: actions/setup-python@v3
+      with:
+        python-version: '3.x'
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install build
+    - name: Build package
+      run: python -m build
+    - name: Publish package
+      uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}
diff --git a/README.md b/README.md
index 4a1d024..7a011be 100644
--- a/README.md
+++ b/README.md
@@ -17,18 +17,25 @@ AdaptiveResonanceLib includes implementations for the following ART models:
     - Fuzzy ART
     - Quadratic Neuron ART
     - Dual Vigilance ART
+    - Gram ART
+- #### Metric Informed 
+    - CVI ART  
+    - iCVI Fuzzy ART  
 - #### Topological
     - Topo ART
 - #### Classification
     - Simple ARTMAP
-      
 - #### Regression
     - ARTMAP
+    - Fusion ART
 - #### Hierarchical
     - DeepARTMAP
     - SMART
-- #### Data Fusion and Reinforcement Learning
+- #### Data Fusion
     - Fusion ART
+- #### Reinforcement Learning
+    - FALCON
+  
 - #### Biclustering
     - Biclustering ARTMAP
 
diff --git a/artlib/__init__.py b/artlib/__init__.py
index 6fb4b5b..2dfdd95 100644
--- a/artlib/__init__.py
+++ b/artlib/__init__.py
@@ -1,6 +1,7 @@
 from artlib.common.BaseART import BaseART
 from artlib.common.BaseARTMAP import BaseARTMAP
-from artlib.common.utils import normalize, compliment_code
+from artlib.common.utils import normalize, compliment_code, de_compliment_code, de_normalize
+from artlib.common.VAT import VAT
 
 from artlib.elementary.ART1 import ART1
 from artlib.elementary.ART2 import ART2A
@@ -12,6 +13,9 @@
 from artlib.elementary.HypersphereART import HypersphereART
 from artlib.elementary.QuadraticNeuronART import QuadraticNeuronART
 
+from artlib.cvi.iCVIFuzzyArt import iCVIFuzzyART
+from artlib.cvi.CVIART import CVIART
+
 from artlib.supervised.ARTMAP import ARTMAP, SimpleARTMAP
 
 from artlib.hierarchical.SMART import SMART
@@ -19,6 +23,8 @@
 
 from artlib.fusion.FusionART import FusionART
 
+from artlib.reinforcement.FALCON import FALCON, TD_FALCON
+
 from artlib.biclustering.BARTMAP import BARTMAP
 
 from artlib.topological.TopoART import TopoART
@@ -44,4 +50,6 @@
     "SMART",
     "FusionART",
     "BARTMAP",
+    "iCVIFuzzyART",
+    "CVIART"
 ]
\ No newline at end of file
diff --git a/artlib/biclustering/BARTMAP.py b/artlib/biclustering/BARTMAP.py
index 2b39fda..c8801fe 100644
--- a/artlib/biclustering/BARTMAP.py
+++ b/artlib/biclustering/BARTMAP.py
@@ -298,7 +298,6 @@ def fit(self, X: np.ndarray, max_iter=1):
 
         for _ in range(max_iter):
             for k in range(n):
-                print(k, self.module_a.n_clusters)
                 self.module_a.pre_step_fit(X)
                 c_a = self.step_fit(X_a, k)
                 self.module_a.labels_[k] = c_a
diff --git a/artlib/common/BaseART.py b/artlib/common/BaseART.py
index 115c025..3e11ce4 100644
--- a/artlib/common/BaseART.py
+++ b/artlib/common/BaseART.py
@@ -1,11 +1,13 @@
 import numpy as np
-from typing import Optional, Callable, Iterable
+from typing import Optional, Callable, Iterable, Literal, List, Tuple
+from copy import deepcopy
 from collections import defaultdict
 from matplotlib.axes import Axes
 from warnings import warn
 from sklearn.base import BaseEstimator, ClusterMixin
 from sklearn.utils.validation import check_is_fitted
-from artlib.common.utils import normalize
+from artlib.common.utils import normalize, de_normalize
+import operator
 
 
 class BaseART(BaseEstimator, ClusterMixin):
@@ -21,6 +23,8 @@ def __init__(self, params: dict):
         self.params = params
         self.sample_counter_ = 0
         self.weight_sample_counter_: list[int] = []
+        self.d_min_ = None
+        self.d_max_ = None
 
     def __getattr__(self, key):
         if key in self.params:
@@ -90,8 +94,7 @@ def set_params(self, **params):
         return self
 
 
-    @staticmethod
-    def prepare_data(X: np.ndarray) -> np.ndarray:
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
         """
         prepare data for clustering
 
@@ -101,7 +104,20 @@ def prepare_data(X: np.ndarray) -> np.ndarray:
         Returns:
             normalized data
         """
-        return normalize(X)
+        normalized, self.d_max_, self.d_min_ = normalize(X, self.d_max_, self.d_min_)
+        return normalized
+
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return de_normalize(X, d_max=self.d_max_, d_min=self.d_min_)
 
     @property
     def n_clusters(self) -> int:
@@ -183,7 +199,7 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         """
         raise NotImplementedError
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
+    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, op: Callable = operator.ge) -> tuple[bool, dict]:
         """
         get the binary match criterion of the cluster
 
@@ -197,7 +213,13 @@ def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache:
             cluster match criterion binary, cache used for later processing
 
         """
-        raise NotImplementedError
+        M, cache = self.match_criterion(i, w, params=params, cache=cache)
+        M_bin = op(M, params["rho"])
+        if cache is None:
+            cache = dict()
+        cache["match_criterion"] = M
+        cache["match_criterion_bin"] = M_bin
+        return M_bin, cache
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -253,7 +275,42 @@ def set_weight(self, idx: int, new_w: np.ndarray):
         self.weight_sample_counter_[idx] += 1
         self.W[idx] = new_w
 
-    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -> int:
+    def _match_tracking(self, cache: dict, epsilon: float, params: dict, method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        M = cache["match_criterion"]
+        if method == "MT+":
+            self.params["rho"] = M+epsilon
+            return True
+        elif method == "MT-":
+            self.params["rho"] = M - epsilon
+            return True
+        elif method == "MT0":
+            self.params["rho"] = M
+            return True
+        elif method == "MT1":
+            self.params["rho"] = np.inf
+            return False
+        elif method == "MT~":
+            return True
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
+
+    @staticmethod
+    def _match_tracking_operator(method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> Callable:
+        if method in ["MT+", "MT-","MT1"]:
+            return operator.ge
+        elif method in ["MT0", "MT~"]:
+            return operator.gt
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
+
+    def _set_params(self, new_params):
+        self.params = new_params
+
+    def _deep_copy_params(self) -> dict:
+        return deepcopy(self.params)
+
+
+    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0) -> int:
         """
         fit the model to a single sample
 
@@ -262,37 +319,64 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
         - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
             Permits external factors to influence cluster creation.
             Returns True if the cluster is valid for the sample, False otherwise
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
 
         Returns:
             cluster label of the input sample
 
         """
         self.sample_counter_ += 1
+        base_params = self._deep_copy_params()
+        mt_operator = self._match_tracking_operator(match_reset_method)
         if len(self.W) == 0:
             w_new = self.new_weight(x, self.params)
             self.add_weight(w_new)
             return 0
         else:
-            T_values, T_cache = zip(*[self.category_choice(x, w, params=self.params) for w in self.W])
+
+            if match_reset_method in ["MT~"] and match_reset_func is not None:
+                T_values, T_cache = zip(*[
+                    self.category_choice(x, w, params=self.params)
+                    if match_reset_func(x, w, c_, params=self.params, cache=None)
+                    else (np.nan, None)
+                    for c_, w in enumerate(self.W)
+                ])
+            else:
+                T_values, T_cache = zip(*[self.category_choice(x, w, params=self.params) for w in self.W])
             T = np.array(T_values)
-            while any(T > 0):
-                c_ = int(np.argmax(T))
+            while any(~np.isnan(T)):
+                c_ = int(np.nanargmax(T))
                 w = self.W[c_]
                 cache = T_cache[c_]
-                m, cache = self.match_criterion_bin(x, w, params=self.params, cache=cache)
-                no_match_reset = (
-                        match_reset_func is None or
-                        match_reset_func(x, w, c_, params=self.params, cache=cache)
-                )
+                m, cache = self.match_criterion_bin(x, w, params=self.params, cache=cache, op=mt_operator)
+                if match_reset_method in ["MT~"] and match_reset_func is not None:
+                    no_match_reset = True
+                else:
+                    no_match_reset = (
+                            match_reset_func is None or
+                            match_reset_func(x, w, c_, params=self.params, cache=cache)
+                    )
                 if m and no_match_reset:
                     self.set_weight(c_, self.update(x, w, self.params, cache=cache))
+                    self._set_params(base_params)
                     return c_
                 else:
-                    T[c_] = -1
+                    T[c_] = np.nan
+                    if m and not no_match_reset:
+                        keep_searching = self._match_tracking(cache, epsilon, self.params, match_reset_method)
+                        if not keep_searching:
+                            T[:] = np.nan
 
             c_new = len(self.W)
             w_new = self.new_weight(x, self.params)
             self.add_weight(w_new)
+            self._set_params(base_params)
             return c_new
 
     def step_pred(self, x) -> int:
@@ -334,8 +418,19 @@ def post_step_fit(self, X: np.ndarray):
         # this is where pruning steps can go
         pass
 
+    def post_fit(self, X: np.ndarray):
+        """
+        undefined function called after fit. Useful for cluster pruning
+
+        Parameters:
+        - X: data set
+
+        """
+        # this is where pruning steps can go
+        pass
+
 
-    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: Optional[Callable] = None, max_iter=1):
+    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: Optional[Callable] = None, max_iter=1, match_reset_method:Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0, verbose: bool = False):
         """
         Fit the model to the data
 
@@ -346,6 +441,12 @@ def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: O
             Permits external factors to influence cluster creation.
             Returns True if the cluster is valid for the sample, False otherwise
         - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         self.validate_data(X)
@@ -355,15 +456,21 @@ def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: O
         self.W: list[np.ndarray] = []
         self.labels_ = np.zeros((X.shape[0], ), dtype=int)
         for _ in range(max_iter):
-            for i, x in enumerate(X):
+            if verbose:
+                from tqdm import tqdm
+                x_iter = tqdm(enumerate(X), total=int(X.shape[0]))
+            else:
+                x_iter = enumerate(X)
+            for i, x in x_iter:
                 self.pre_step_fit(X)
-                c = self.step_fit(x, match_reset_func=match_reset_func)
+                c = self.step_fit(x, match_reset_func=match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
                 self.labels_[i] = c
                 self.post_step_fit(X)
+        self.post_fit(X)
         return self
 
 
-    def partial_fit(self, X: np.ndarray, match_reset_func: Optional[Callable] = None):
+    def partial_fit(self, X: np.ndarray, match_reset_func: Optional[Callable] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
         """
         iteratively fit the model to the data
 
@@ -372,6 +479,12 @@ def partial_fit(self, X: np.ndarray, match_reset_func: Optional[Callable] = None
         - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
             Permits external factors to influence cluster creation.
             Returns True if the cluster is valid for the sample, False otherwise
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
 
@@ -387,7 +500,7 @@ def partial_fit(self, X: np.ndarray, match_reset_func: Optional[Callable] = None
             j = len(self.labels_)
             self.labels_ = np.pad(self.labels_, [(0, X.shape[0])], mode='constant')
         for i, x in enumerate(X):
-            c = self.step_fit(x, match_reset_func=match_reset_func)
+            c = self.step_fit(x, match_reset_func=match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
             self.labels_[i+j] = c
         return self
 
@@ -414,6 +527,9 @@ def predict(self, X: np.ndarray) -> np.ndarray:
             y[i] = c
         return y
 
+    def shrink_clusters(self, shrink_ratio: float = 0.1):
+        return self
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         undefined function for visualizing the bounds of each cluster
@@ -426,6 +542,15 @@ def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         raise NotImplementedError
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        undefined function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        raise NotImplementedError
+
+
     def visualize(
             self,
             X: np.ndarray,
diff --git a/artlib/common/BaseARTMAP.py b/artlib/common/BaseARTMAP.py
index e92a52f..d54f4b1 100644
--- a/artlib/common/BaseARTMAP.py
+++ b/artlib/common/BaseARTMAP.py
@@ -1,5 +1,5 @@
 import numpy as np
-from typing import Union, Optional, Iterable
+from typing import Union, Optional, Iterable, Literal
 from collections import defaultdict
 from matplotlib.axes import Axes
 from sklearn.base import BaseEstimator, ClassifierMixin, ClusterMixin
@@ -75,7 +75,7 @@ def validate_data(self, X: np.ndarray, y: np.ndarray):
         """
         raise NotImplementedError
 
-    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
+    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10):
         """
         Fit the model to the data
 
@@ -87,7 +87,7 @@ def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
         """
         raise NotImplementedError
 
-    def partial_fit(self, X: np.ndarray, y: np.ndarray):
+    def partial_fit(self, X: np.ndarray, y: np.ndarray, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10):
         """
         Partial fit the model to the data
 
diff --git a/artlib/common/VAT.py b/artlib/common/VAT.py
new file mode 100644
index 0000000..37fc7ae
--- /dev/null
+++ b/artlib/common/VAT.py
@@ -0,0 +1,44 @@
+# Bezdek, J. C., & Hathaway, R. J. (2002). VAT: A tool for visual assessment of cluster tendency.
+# Proceedings of the 2002 International Joint Conference on Neural Networks. doi:10.1109/IJCNN.2002.1007487
+import numpy as np
+from typing import Optional, Tuple, Callable
+from scipy.spatial.distance import pdist, squareform
+
+
+def VAT(data: np.ndarray, distance_metric: Optional[Callable] = lambda X: pdist(X, "euclidean")) -> Tuple[np.ndarray, np.ndarray]:
+    """
+    Visual Assessment of cluster Tendency (VAT) algorithm.
+
+    Parameters:
+    - data: Input data set as a 2D numpy array where each row is a sample.
+    - distance_metric: Optional callable function to calculate pairwise distances.
+      Defaults to Euclidean distance using pdist.
+      If None, assume data is a pre-computed distance matrix
+
+    Returns:
+        Tuple containing:
+        - Reordered distance matrix reflecting cluster structure.
+        - Reordered list of indices indicating the optimal clustering order.
+    """
+    if distance_metric is None:
+        pairwise_dist = data
+    else:
+        pairwise_dist = squareform(distance_metric(data))
+
+    num_samples = data.shape[0]
+    indicies = []
+    remaining = list(range(num_samples))
+
+    ix, jx = np.unravel_index(pairwise_dist.argmax(), pairwise_dist.shape)
+    indicies.append(ix)
+    remaining.pop(ix)
+
+    while remaining:
+        sub_matrix = pairwise_dist[np.ix_(indicies, remaining)]
+        _, jx = np.unravel_index(sub_matrix.argmin(), sub_matrix.shape)
+        indicies.append(remaining[jx])
+        remaining.pop(jx)
+
+    return pairwise_dist[np.ix_(indicies, indicies)], np.array(indicies)
+
+
diff --git a/artlib/common/utils.py b/artlib/common/utils.py
index 106e029..285e804 100644
--- a/artlib/common/utils.py
+++ b/artlib/common/utils.py
@@ -1,17 +1,44 @@
 import numpy as np
+from typing import Tuple, Optional
 
-def normalize(data: np.ndarray) -> np.ndarray:
+
+def normalize(data: np.ndarray, d_max: Optional[np.ndarray] = None, d_min: Optional[np.ndarray] = None) -> Tuple[
+    np.ndarray, np.ndarray, np.ndarray]:
     """
-    normalize data betweeon 0 and 1
+    Normalize data column-wise between 0 and 1.
 
     Parameters:
-    - data: data set
+    - data: 2D array of data set (rows = samples, columns = features)
+    - d_max: Optional, maximum values for each column
+    - d_min: Optional, minimum values for each column
 
     Returns:
-        normalized data
+    - normalized: normalized data
+    - d_max: maximum values for each column
+    - d_min: minimum values for each column
     """
-    normalized = (data-np.min(data))/(np.max(data)-np.min(data))
-    return normalized
+    if d_min is None:
+        d_min = np.min(data, axis=0)
+    if d_max is None:
+        d_max = np.max(data, axis=0)
+
+    normalized = (data - d_min) / (d_max - d_min)
+    return normalized, d_max, d_min
+
+
+def de_normalize(data: np.ndarray, d_max: np.ndarray, d_min: np.ndarray) -> np.ndarray:
+    """
+    Restore column-wise normalized data to original scale.
+
+    Parameters:
+    - data: normalized data
+    - d_max: maximum values for each column
+    - d_min: minimum values for each column
+
+    Returns:
+    - De-normalized data
+    """
+    return data * (d_max - d_min) + d_min
 
 def compliment_code(data: np.ndarray) -> np.ndarray:
     """
@@ -26,6 +53,34 @@ def compliment_code(data: np.ndarray) -> np.ndarray:
     cc_data = np.hstack([data, 1.0-data])
     return cc_data
 
+def de_compliment_code(data: np.ndarray) -> np.ndarray:
+    """
+    finds centroid of compliment coded data
+
+    Parameters:
+    - data: data set
+
+    Returns:
+        compliment coded data
+    """
+    # Get the shape of the array
+    n, total_columns = data.shape
+
+    # Ensure the number of columns is even so that it can be split evenly
+    assert total_columns % 2 == 0, "The number of columns must be even"
+
+    # Calculate the number of columns in each resulting array
+    m = total_columns // 2
+
+    # Split the array into two arrays of shape n x m
+    arr1 = data[:, :m]
+    arr2 = 1-data[:, m:]
+
+    # Find the element-wise mean
+    mean_array = (arr1 + arr2) / 2
+
+    return mean_array
+
 def l1norm(x: np.ndarray) -> float:
     """
     get l1 norm of a vector
diff --git a/artlib/cvi/CVIART.py b/artlib/cvi/CVIART.py
new file mode 100644
index 0000000..62cfc81
--- /dev/null
+++ b/artlib/cvi/CVIART.py
@@ -0,0 +1,216 @@
+import numpy as np
+from copy import deepcopy
+import sklearn.metrics as metrics
+from typing import Optional, List, Callable, Literal, Iterable
+from matplotlib.axes import Axes
+from artlib.common.BaseART import BaseART
+
+
+class CVIART(BaseART):
+    """CVI Art Classification
+
+    Expanded version of Art that uses Cluster Validity Indicies to help with cluster selection.
+    PBM is not implemented, can be seen here.
+    https://git.mst.edu/acil-group/CVI-Fuzzy-ART/-/blob/master/PBM_index.m?ref_type=heads
+
+    Note, the default step_fit function in base ART evaluates the matching function even if
+    the other criteria has failed. This means it could run slower then it would otherwise.
+
+
+    Parameters:
+        rho: float [0,1] for the vigilance parameter.
+        alpha: float choice parameter. 1e-7 recommended value.
+        beta: float [0,1] learning parameters. beta = 1 is fast learning recommended value.
+        validity: int the cluster validity index being used.
+        W: list of weights, top down.
+        labels: class labels for data set.
+    """
+    CALINSKIHARABASZ = 1
+    DAVIESBOULDIN = 2
+    SILHOUETTE = 3
+    # PBM = 4
+
+    def validate_params(self, params: dict):
+        """
+        validate clustering parameters
+
+        Parameters:
+        - params: dict containing parameters for the algorithm
+
+        """
+        self.base_module.validate_params(params)
+        assert 'validity' in params
+        assert isinstance(params['validity'], int)
+        assert params["validity"] in [CVIART.CALINSKIHARABASZ, CVIART.DAVIESBOULDIN, CVIART.SILHOUETTE]
+
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            normalized data
+        """
+        return self.base_module.prepare_data(X)
+
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return self.base_module.restore_data(X)
+
+    @property
+    def W(self):
+        return self.base_module.W
+
+    @W.setter
+    def W(self, new_W):
+        self.base_module.W = new_W
+
+    @property
+    def labels_(self):
+        return self.base_module.labels_
+
+    @labels_.setter
+    def labels_(self, new_labels_):
+        self.base_module.labels_ = new_labels_
+
+
+    def __init__(self, base_module: BaseART, validity: int):
+        self.base_module = base_module
+        params = dict(base_module.params, **{"validity": validity})
+        super().__init__(params)
+        print(self.params)
+
+
+    def CVI_match(self, x, w, c_, params, extra, cache):
+        if len(self.W) < 2:
+            return True
+
+        if extra['validity'] == self.CALINSKIHARABASZ:
+            valid_func = metrics.calinski_harabasz_score
+        elif extra['validity'] == self.DAVIESBOULDIN:
+            valid_func = metrics.davies_bouldin_score
+        elif extra['validity'] == self.SILHOUETTE:
+            valid_func = metrics.silhouette_score
+        else:
+            raise ValueError(f"Invalid Validity Parameter: {extra['validity']}")
+
+        old_VI = valid_func(self.data, self.labels_)
+        new_labels = np.copy(self.labels_)
+        new_labels[extra["index"]] = c_
+        new_VI = valid_func(self.data, new_labels)
+        if extra['validity'] != self.DAVIESBOULDIN:
+            return new_VI > old_VI
+        else:
+            return new_VI < old_VI
+    def _match_tracking(self, cache: dict, epsilon: float, params: dict, method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        M = cache["match_criterion"]
+        if method == "MT+":
+            self.base_module.params["rho"] = M+epsilon
+            return True
+        elif method == "MT-":
+            self.base_module.params["rho"] = M - epsilon
+            return True
+        elif method == "MT0":
+            self.base_module.params["rho"] = M
+            return True
+        elif method == "MT1":
+            self.base_module.params["rho"] = np.inf
+            return False
+        elif method == "MT~":
+            return True
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
+
+    def _set_params(self, new_params):
+        self.base_module.params = new_params
+
+    def _deep_copy_params(self) -> dict:
+        return deepcopy(self.base_module.params)
+
+    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: Optional[Callable] = None,
+            max_iter=1, match_reset_method:Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
+        """
+        Fit the model to the data
+
+        Parameters:
+        - X: data set
+        - y: not used. For compatibility.
+        - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
+            Permits external factors to influence cluster creation.
+            Returns True if the cluster is valid for the sample, False otherwise
+        - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
+        """
+        self.data = X
+        self.base_module.validate_data(X)
+        self.base_module.check_dimensions(X)
+        self.is_fitted_ = True
+
+        self.W: list[np.ndarray] = []
+        self.labels_ = np.zeros((X.shape[0], ), dtype=int)
+        for _ in range(max_iter):
+            for index, x in enumerate(X):
+                self.pre_step_fit(X)
+                if match_reset_func is None:
+                    cvi_match_reset_func = lambda i, w, cluster_a, params, cache: self.CVI_match(i, w, cluster_a, params, {"index": index, "validity":self.params["validity"]}, cache)
+                else:
+                    cvi_match_reset_func = lambda i, w, cluster_a, params, cache: (match_reset_func(i,w,cluster_a,params,cache) and self.CVI_match(i, w, cluster_a, params, {"index": index, "validity":self.params["validity"]}, cache))
+                c = self.base_module.step_fit(x, match_reset_func=cvi_match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
+                self.labels_[index] = c
+                self.post_step_fit(X)
+
+
+    def pre_step_fit(self, X: np.ndarray):
+        return self.base_module.pre_step_fit(X)
+
+    def post_step_fit(self, X: np.ndarray):
+        return self.base_module.post_step_fit(X)
+
+    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0) -> int:
+        """
+        fit the model to a single sample
+
+        Parameters:
+        - x: data sample
+        - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
+            Permits external factors to influence cluster creation.
+            Returns True if the cluster is valid for the sample, False otherwise
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
+
+        Returns:
+            cluster label of the input sample
+
+        """
+        raise NotImplementedError
+
+    def step_pred(self, x) -> int:
+        return self.base_module.step_pred(x)
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        return self.base_module.get_cluster_centers()
+
+    def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
+        return self.base_module.plot_cluster_bounds(ax, colors,linewidth)
+
diff --git a/artlib/cvi/__init__.py b/artlib/cvi/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/artlib/cvi/iCVIFuzzyArt.py b/artlib/cvi/iCVIFuzzyArt.py
new file mode 100644
index 0000000..e16bf3e
--- /dev/null
+++ b/artlib/cvi/iCVIFuzzyArt.py
@@ -0,0 +1,94 @@
+"""
+Add Reference in correct format.
+The original matlab code can be found at https://github.com/ACIL-Group/iCVI-toolbox/tree/master
+The formulation is available at
+https://scholarsmine.mst.edu/cgi/viewcontent.cgi?article=3833&context=doctoral_dissertations Pages 314-316 and 319-320
+Extended icvi offline mode can be found at
+https://ieeexplore.ieee.org/document/9745260
+"""
+import numpy as np
+from typing import Optional, Literal, Callable
+from artlib.elementary.FuzzyART import FuzzyART
+from artlib.cvi.iCVIs.CalinkskiHarabasz import iCVI_CH
+
+
+class iCVIFuzzyART(FuzzyART):
+    """iCVI Fuzzy Art Classification
+
+    Parameters:
+        rho: float [0,1] for the vigilance parameter.
+        alpha: float choice parameter. 1e-7 recommended value.
+        beta: float [0,1] learning parameters. beta = 1 is fast learning recommended value.
+        validity: int the cluster validity index being used.
+        W: list of weights, top down.
+        labels: class labels for data set.
+    """
+    CALINSKIHARABASZ = 1
+
+    def __init__(self, rho: float, alpha: float, beta: float, validity: int, offline: bool = True):
+        super().__init__(rho, alpha, beta)
+        self.params['validity'] = validity  # Currently not used. Waiting for more algorithms.
+        self.offline = offline
+        assert 'validity' in self.params  # Because Fuzzy art doesn't accept validity, and makes the params the way it does, validations have to be done after init.
+        assert isinstance(self.params['validity'], int)
+
+    def iCVI_match(self, x, w, c_, params, cache):
+        if self.offline:
+            new = self.iCVI.switch_label(x, self.labels_[self.index], c_)
+        else:
+            new = self.iCVI.add_sample(x, c_)
+        # Eventually this should be an icvi function that you pass the params, and it handles if this is true or false.
+        return new['criterion_value'] > self.iCVI.criterion_value
+        # return self.iCVI.evalLabel(x, c_) This except pass params instead.
+
+    # Could add max epochs back in, but only if offline is true, or do something special...
+    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_func: Optional[Callable] = None, max_iter=1, match_reset_method:Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
+        """
+        Fit the model to the data
+
+        Parameters:
+        - X: data set
+        - y: not used. For compatibility.
+        - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
+            Permits external factors to influence cluster creation.
+            Returns True if the cluster is valid for the sample, False otherwise
+        - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
+        """
+        self.validate_data(X)
+        self.check_dimensions(X)
+        self.is_fitted_ = True
+
+        self.W: list[np.ndarray] = []
+        self.labels_ = np.zeros((X.shape[0], ), dtype=int)
+
+        self.iCVI = iCVI_CH(X[0])
+
+        if self.offline:
+            for x in X:
+                params = self.iCVI.add_sample(x, 0)
+                self.iCVI.update(params)
+
+        for i, x in enumerate(X):
+            self.pre_step_fit(X)
+            self.index = i
+            if match_reset_func is None:
+                c = self.step_fit(x, match_reset_func=self.iCVI_match, match_reset_method=match_reset_method, epsilon=epsilon)
+            else:
+                match_reset_func = lambda x, w, c_, params, cache: (match_reset_func(x, w, c_, params, cache) & self.iCVI_match(x, w, c_, params, cache))
+                c = self.step_fit(x, match_reset_func=match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
+
+            if self.offline:
+                params = self.iCVI.switch_label(x, self.labels_[i], c)
+            else:
+                params = self.iCVI.add_sample(x, c)
+            self.iCVI.update(params)
+
+            self.labels_[i] = c
+            self.post_step_fit(X)
diff --git a/artlib/cvi/iCVIs/CalinkskiHarabasz.py b/artlib/cvi/iCVIs/CalinkskiHarabasz.py
new file mode 100644
index 0000000..6f85989
--- /dev/null
+++ b/artlib/cvi/iCVIs/CalinkskiHarabasz.py
@@ -0,0 +1,289 @@
+"""
+Some things to consider in the future.
+
+Removing entire labels from the dataset should be possible.
+
+Right now I think its not possible for Fuzzy Art to delete a cluster, so no need to do so now.
+
+Creating functions to explictly add a sample, remove or switch, instead of requiring the update call would be nice
+Or at least change the function names for add, remove, and switch, since they imply that its done, not that it will
+update after update is called.
+"""
+import numpy as np
+
+
+def delta_add_sample_to_average(average, sample, total_samples):
+    """Calculate the new average if sample is added"""
+    return (sample - average) / total_samples
+
+
+def delta_remove_sample_from_average(average, sample, total_samples):
+    """Calculate the new average if sample is removed"""
+    return (average - sample) / (total_samples - 1)
+
+
+class iCVI_CH():
+    """Implementation of the Calinski Harabasz Validity Index in incremental form.
+
+    Expanded implementation of the incremental version of the Calinski Harabasz Cluster Validity Index.
+
+    The original matlab code can be found at https://github.com/ACIL-Group/iCVI-toolbox/blob/master/classes/CVI_CH.m
+    The formulation is available at
+    https://scholarsmine.mst.edu/cgi/viewcontent.cgi?article=3833&context=doctoral_dissertations Pages 314-316 and 319-320
+
+    This implementation returns a dictionary of updated parameters when calling functions, which can then be passed with
+    the update function to accept the changes. This allows for testing changes/additions to the categories without doing a
+    deep copy of the object.
+
+    In addition, the calculations for removing a sample, or switching the label of a sample from the dataset are included.
+    This allows for very efficient calculations on clustering algorithms that would like to prune or adjust the labels of
+    samples in the dataset.
+
+    For the Calinski Harabasz validity Index, larger values represent better clusters.
+
+    Parameters:
+        dim: an int storing the dimensionality of the input data.
+        n_samples: an int with the total number of samples.
+        mu: a numpy array representing the average of all data points.
+        CD: ClusterData... a dict holding all the data for each cluster, with the parameters
+            n: the number of samples belonging to each cluster/label.
+            v: the prototypes/centriods of each cluster.
+            CP: the compactness of each cluster. Not really needed since WGSS can be calculated incrementally.
+            G: the vector g of each cluster.
+        WGSS: Within Groupd sum of squares.
+        criterion_value: the calculated CH score.
+    """
+
+    def __init__(self, x: np.ndarray) -> None:
+        """Create the iCVI object.
+
+        Args:
+            x: a sample from the dataset used for recording data dimensionality.
+        """
+        self.dim = x.shape[0]  # Dimension of the input data
+        self.n_samples: int = 0  # number of samples encountered
+        self.mu = np.array([])  # geometric mean of the data
+        self.CD = {}            # Dict for each cluster label containing n,v,CP, and G
+        self.WGSS = 0           # within group sum of squares
+        self.criterion_value = 0  # calcualted CH index
+
+    def add_sample(self, x, label) -> dict:
+        """Calculate the result of adding a new sample with a given label.
+
+        Create a dictionary containing the updated values after assigning a label to a given sample.
+        To accept the outcome of the sample being added, pass the returned parameter dict to update.
+
+        In general, if newP['criterion_value'] > obj.criterion_value, the clustering has been improved.
+
+        Args:
+            x: The sample to add to the current validity index calculation
+            label: an int representing the sample category/cluster
+
+        Returns:
+            newP: a dictionary contained the values after the sample is added, to be passed to update call.
+        """
+        newP = {'x': x, 'label': label}  # New Parameters
+        newP['n_samples'] = self.n_samples + 1
+        if self.mu.size == 0:  # mu will be size 0 if no samples in dataset.
+            newP['mu'] = x
+        else:
+            newP['mu'] = self.mu + delta_add_sample_to_average(self.mu, x, newP['n_samples'])
+
+        CD = {}
+        newP['CD'] = CD
+        SEP = []  # separation between each cluster and the mean of the data
+        if label not in self.CD:
+            n_clusters = len(self.CD) + 1
+            CD['n'] = 1
+            CD['v'] = x
+            CD['CP'] = 0
+            CD['G'] = np.zeros(self.dim)
+            newP['CP_diff'] = 0
+            diff = x - newP['mu']
+            SEP.append(sum(diff ** 2))  # Handle SEP for new clusters now instead of later.
+        else:
+            n_clusters = len(self.CD)
+            Data = self.CD[label]
+            CD['n'] = Data['n'] + 1
+
+            # The paper defines deltaV = Vold - Vnew, so I need to switch this sign. Consider changing functions to do this.
+            deltaV = -1 * delta_add_sample_to_average(Data['v'], x, CD['n'])
+            CD['v'] = Data['v'] - deltaV  # Vnew = Vold - deltaV
+            diff_x_v = x - CD['v']
+
+            newP['CP_diff'] = (diff_x_v @ diff_x_v.T) + (CD['n'] - 1) * (deltaV @ deltaV.T) + 2 * (deltaV @ Data['G'])
+            CD['CP'] = Data['CP'] + newP['CP_diff']
+            CD['G'] = Data['G'] + diff_x_v + (CD['n'] - 1) * deltaV
+
+        if n_clusters < 2:
+            newP['criterion_value'] = 0
+        else:
+            for i in self.CD:  # A new label won't be processed, but handled earlier
+                if i == label:
+                    n = CD['n']
+                    v = CD['v']
+                else:
+                    n = self.CD[i]['n']
+                    v = self.CD[i]['v']
+                diff = v - newP['mu']
+                SEP.append(n * sum(diff ** 2))
+
+            WGSS = self.WGSS + newP['CP_diff']
+            BGSS = sum(SEP)  # between-group sum of squares
+            if WGSS == 0:  # this can be 0 if all samples in different clusters, which is a divide by 0 error.
+                newP['criterion_value'] = 0
+            else:
+                newP['criterion_value'] = (BGSS / WGSS) * (newP['n_samples'] - n_clusters) / (n_clusters - 1)
+        return newP
+
+    def update(self, params) -> None:
+        """Update the parameters of the object.
+
+        Takes the updated params from adding/removing a sample or switching its label, and updates the object.
+        Switching a label needs more updates, so those dicts have an extra set of things to update, signified with
+        the 'label2' key existing
+
+        Args:
+            params: dict containing the parameters to update.
+        """
+        self.n_samples = params['n_samples']
+        self.mu = params['mu']
+        self.criterion_value = params['criterion_value']
+        self.CD[params['label']] = params['CD']
+        self.WGSS += params['CP_diff']
+        if 'label2' in params:
+            self.CD[params['label2']] = params['CD2']
+            self.WGSS += params['CP_diff2']
+
+    def switch_label(self, x, label_old, label_new):
+        """Calculates the parameters if a sample has its label changed.
+
+        This essentially removes a sample with the old label from the clusters, then adds it back with the new sample.
+        There are a few optimizations, such as keeping mu the same since adding and removing it doesn't affect any calculations
+        that are needed.
+
+        Otherwise it should work the same as removing a sample and updating, then adding the sample back and updating, without
+        the need to create a deep copy of the object if just testing the operation.
+
+        Args:
+            x: The sample to switch the label of for the current validity index calculation
+            label_old: an int representing the sample category/cluster the sample belongs to
+            label_new: an int representing the sample category/cluster the sample will be assigned to
+
+        Returns:
+            newP: a dictionary contained the values after the sample is added, to be passed to update call."""
+        if label_new == label_old:
+            return {'n_samples': self.n_samples,
+                    'mu': self.mu,
+                    'criterion_value': self.criterion_value,
+                    'label': label_old,
+                    'CP_diff': 0,
+                    'CD': {
+                        'n': self.CD[label_old]['n'],
+                        'v': self.CD[label_old]['v'],
+                        'CP': self.CD[label_old]['CP'],
+                        'G': self.CD[label_old]['G']}
+                    }
+        if self.CD[label_old]['n'] <= 1:
+            raise Exception("Can't remove a value from a cluster of 1")
+
+        newP = {'x': x, 'label': label_old, 'label2': label_new}  # New Parameters
+        newP['mu'] = self.mu
+        newP['n_samples'] = self.n_samples
+        params_remove = self.remove_sample(x, label_old)
+        params_add = self.add_sample(x, label_new)
+        newP['CD'] = params_remove['CD']
+        newP['CP_diff'] = params_remove['CP_diff']
+        newP['CD2'] = params_add['CD']
+        newP['CP_diff2'] = params_add['CP_diff']
+
+        SEP = []  # separation between each cluster and the mean of the data
+        if label_new not in self.CD:
+            n_clusters = len(self.CD) + 1
+            diff = x - newP['mu']
+            SEP.append(sum(diff ** 2))  # Handle SEP for new clusters now instead of later.
+        else:
+            n_clusters = len(self.CD)
+
+        if n_clusters < 2:  # I don't think this can ever happen with remove label not deleting categories.
+            newP['criterion_value'] = 0
+        else:
+            for i in self.CD:
+                if i == label_old:
+                    n = params_remove['CD']['n']
+                    v = params_remove['CD']['v']
+                elif i == label_new:
+                    n = params_add['CD']['n']
+                    v = params_add['CD']['v']
+                else:
+                    n = self.CD[i]['n']
+                    v = self.CD[i]['v']
+                diff = v - newP['mu']
+                SEP.append(n * sum(diff ** 2))
+
+            WGSS = self.WGSS + newP['CP_diff']
+            WGSS += newP['CP_diff2']
+            BGSS = sum(SEP)  # between-group sum of squares
+            if WGSS == 0:  # this can be 0 if all samples in different clusters, which is a divide by 0 error.
+                newP['criterion_value'] = 0
+            else:
+                newP['criterion_value'] = (BGSS / WGSS) * (newP['n_samples'] - n_clusters) / (n_clusters - 1)
+        return newP
+
+    def remove_sample(self, x, label):  # This is left here mostly as an extra, and not really meant to be used.
+        """Remove a sample from the clusters
+
+         Calculates parameters after removing a sample from the clusters, or the opposite of an add operation.
+
+        Args:
+            x: The sample to remove from the current validity index calculation
+            label: an int representing the sample category/cluster
+
+        Returns:
+            newP: a dictionary contained the values after the sample is remove, to be passed to update call.
+        """
+        Data = self.CD[label]
+        if Data['n'] <= 1:
+            raise Exception("Can't remove a value from a cluster of 1")  # At least for now
+
+        newP = {'x': x, 'label': label}  # New Parameters after removal
+        newP['mu'] = self.mu - delta_remove_sample_from_average(self.mu, x, self.n_samples)
+        newP['n_samples'] = self.n_samples - 1
+
+        CD = {}
+        newP['CD'] = CD
+        CD['n'] = Data['n'] - 1
+
+        # We need the delta v from when the sample was added, but the paper defines deltaV = Vold - Vnew, so I need to keep these signs
+        deltaVPrior = delta_remove_sample_from_average(Data['v'], x, Data['n'])
+        CD['v'] = Data['v'] + deltaVPrior  # Vnew + deltaV = Vold
+        diff_x_vPrior = x - Data['v']
+
+        # We already subtracted 1 from n
+        CD['G'] = Data['G'] - (diff_x_vPrior + (CD['n']) * deltaVPrior)
+        # CD's G is the old G, which is what we added before.
+        newP['CP_diff'] = -1 * ((diff_x_vPrior @ diff_x_vPrior.T) + (CD['n']) * (deltaVPrior @ deltaVPrior.T) + 2 * (deltaVPrior @ CD['G']))
+        CD['CP'] = Data['CP'] + newP['CP_diff']
+
+        n_clusters = len(self.CD)  # Move this up if deleting clusters are allowed.
+        if n_clusters < 2:
+            newP['criterion_value'] = 0
+        else:
+            SEP = []  # separation between each cluster and the mean of the data
+            for i in self.CD:
+                if i == label:
+                    n = CD['n']
+                    v = CD['v']
+                else:
+                    n = self.CD[i]['n']
+                    v = self.CD[i]['v']
+                diff = v - newP['mu']
+                SEP.append(n * sum(diff ** 2))
+
+            WGSS = self.WGSS + newP['CP_diff']
+            BGSS = sum(SEP)  # between-group sum of squares
+            if WGSS == 0:  # this can be 0 if all samples in different clusters, which is a divide by 0 error.
+                newP['criterion_value'] = 0
+            else:
+                newP['criterion_value'] = (BGSS / WGSS) * (newP['n_samples'] - n_clusters) / (n_clusters - 1)
+        return newP
diff --git a/artlib/elementary/ART1.py b/artlib/elementary/ART1.py
index 79660fd..2a04404 100644
--- a/artlib/elementary/ART1.py
+++ b/artlib/elementary/ART1.py
@@ -5,13 +5,9 @@
 """
 
 import numpy as np
-from typing import Optional
+from typing import Optional, List
 from artlib.common.BaseART import BaseART
-from artlib.common.utils import normalize, l1norm
-
-def prepare_data(data: np.ndarray) -> np.ndarray:
-    normalized = normalize(data)
-    return normalized
+from artlib.common.utils import l1norm
 
 
 class ART1(BaseART):
@@ -94,22 +90,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         w_td = w[self.dim_:]
         return l1norm(np.logical_and(i, w_td)) / l1norm(i), cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params, cache)
-        return M >= params["rho"], cache
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -148,3 +128,11 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
         w_td_new = i
         w_bu_new = (params["L"] / (params["L"] - 1 + self.dim_))*w_td_new
         return np.concatenate([w_bu_new, w_td_new])
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return [w[self.dim_:] for w in self.W]
diff --git a/artlib/elementary/ART2.py b/artlib/elementary/ART2.py
index 1bc11e9..497e0a0 100644
--- a/artlib/elementary/ART2.py
+++ b/artlib/elementary/ART2.py
@@ -18,7 +18,7 @@
 """
 
 import numpy as np
-from typing import Optional
+from typing import Optional, List
 from warnings import warn
 from artlib.common.BaseART import BaseART
 
@@ -117,22 +117,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         else:
             return M, cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params, cache)
-        return M >= params["rho"], cache
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -164,3 +148,11 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
 
         """
         return i
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return self.W
\ No newline at end of file
diff --git a/artlib/elementary/BayesianART.py b/artlib/elementary/BayesianART.py
index a498c0b..d092384 100644
--- a/artlib/elementary/BayesianART.py
+++ b/artlib/elementary/BayesianART.py
@@ -4,16 +4,12 @@
 IEEE Transactions on Neural Networks, 18, 1628–1644. doi:10.1109/TNN.2007.900234.
 """
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, List, Callable, Literal, Tuple
+import operator
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
-from artlib.common.utils import normalize
 from artlib.common.visualization import plot_gaussian_contours_covariance
 
-def prepare_data(data: np.ndarray) -> np.ndarray:
-    normalized = normalize(data)
-    return normalized
-
 
 class BayesianART(BaseART):
     # implementation of Bayesian ART
@@ -118,7 +114,7 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         # return cache["det_cov"]
         return np.linalg.det(new_cov), cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
+    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, op: Callable = operator.ge) -> tuple[bool, dict]:
         """
         get the binary match criterion of the cluster
 
@@ -133,7 +129,33 @@ def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache:
 
         """
         M, cache = self.match_criterion(i, w, params=params, cache=cache)
-        return M <= params["rho"], cache
+        M_bin = op(params["rho"], M) # note that this is backwards from the base ART: rho >= M
+        if cache is None:
+            cache = dict()
+        cache["match_criterion"] = M
+        cache["match_criterion_bin"] = M_bin
+
+        return M_bin, cache
+
+    def _match_tracking(self, cache: dict, epsilon: float, params: dict, method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        M = cache["match_criterion"]
+        # we have to reverse some signs becayse bayesianART has an inverted vigilence check
+        if method == "MT+":
+            self.params["rho"] = M - epsilon
+            return True
+        elif method == "MT-":
+            self.params["rho"] = M + epsilon
+            return True
+        elif method == "MT0":
+            self.params["rho"] = M
+            return True
+        elif method == "MT1":
+            self.params["rho"] = -np.inf
+            return False
+        elif method == "MT~":
+            return True
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -184,6 +206,13 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
         """
         return np.concatenate([i, params["cov_init"].flatten(), [1]])
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return [w[:self.dim_] for w in self.W]
 
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
diff --git a/artlib/elementary/DualVigilanceART.py b/artlib/elementary/DualVigilanceART.py
index 60f3e36..cfb912a 100644
--- a/artlib/elementary/DualVigilanceART.py
+++ b/artlib/elementary/DualVigilanceART.py
@@ -4,8 +4,9 @@
 Neural Networks, 109, 1–5. doi:10.1016/j.neunet.2018.09.015.
 """
 import numpy as np
-from typing import Optional, Callable, Iterable
+from typing import Optional, Callable, Iterable, List, Literal
 from warnings import warn
+from copy import deepcopy
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 
@@ -42,6 +43,18 @@ def prepare_data(self, X: np.ndarray) -> np.ndarray:
         """
         return self.base_module.prepare_data(X)
 
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return self.base_module.restore_data(X)
+
     def get_params(self, deep: bool = True) -> dict:
         """
 
@@ -131,7 +144,32 @@ def validate_params(self, params: dict):
         assert params["rho_lower_bound"] >= 0
         assert isinstance(params["rho_lower_bound"], float)
 
-    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -> int:
+    def _match_tracking(self, cache: dict, epsilon: float, params: dict, method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        M = cache["match_criterion"]
+        if method == "MT+":
+            self.base_module.params["rho"] = M+epsilon
+            return True
+        elif method == "MT-":
+            self.base_module.params["rho"] = M - epsilon
+            return True
+        elif method == "MT0":
+            self.base_module.params["rho"] = M
+            return True
+        elif method == "MT1":
+            self.base_module.params["rho"] = np.inf
+            return False
+        elif method == "MT~":
+            return True
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
+
+    def _set_params(self, new_params):
+        self.base_module.params = new_params
+
+    def _deep_copy_params(self) -> dict:
+        return deepcopy(self.base_module.params)
+
+    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None,match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0) -> int:
         """
         fit the model to a single sample
 
@@ -145,6 +183,8 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
             cluster label of the input sample
 
         """
+        base_params = self._deep_copy_params()
+        mt_operator = self._match_tracking_operator(match_reset_method)
         self.sample_counter_ += 1
         if len(self.base_module.W) == 0:
             new_w = self.base_module.new_weight(x, self.base_module.params)
@@ -160,10 +200,10 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
             )
             T = np.array(T_values)
             while any(T > 0):
-                c_ = int(np.argmax(T))
+                c_ = int(np.nanargmax(T))
                 w = self.base_module.W[c_]
                 cache = T_cache[c_]
-                m1, cache = self.base_module.match_criterion_bin(x, w, params=self.base_module.params, cache=cache)
+                m1, cache = self.base_module.match_criterion_bin(x, w, params=self.base_module.params, cache=cache, op=mt_operator)
                 no_match_reset = (
                     match_reset_func is None or
                     match_reset_func(x, w, self.map[c_], params=self.base_module.params, cache=cache)
@@ -173,22 +213,29 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
                     if m1:
                         new_w = self.base_module.update(x, w, self.base_module.params, cache=cache)
                         self.base_module.set_weight(c_, new_w)
+                        self._set_params(base_params)
                         return self.map[c_]
                     else:
                         lb_params = dict(self.base_module.params, **{"rho": self.rho_lower_bound})
-                        m2, _ = self.base_module.match_criterion_bin(x, w, params=lb_params, cache=cache)
+                        m2, _ = self.base_module.match_criterion_bin(x, w, params=lb_params, cache=cache, op=mt_operator)
                         if m2:
                             c_new = len(self.base_module.W)
                             w_new = self.base_module.new_weight(x, self.base_module.params)
                             self.base_module.add_weight(w_new)
                             self.map[c_new] = self.map[c_]
+                            self._set_params(base_params)
                             return self.map[c_new]
-                T[c_] = -1
+                else:
+                    keep_searching = self._match_tracking(cache, epsilon, self.params, match_reset_method)
+                    if not keep_searching:
+                        T[:] = np.nan
+                T[c_] = np.nan
 
             c_new = len(self.base_module.W)
             w_new = self.base_module.new_weight(x, self.base_module.params)
             self.base_module.add_weight(w_new)
             self.map[c_new] = max(self.map.values()) + 1
+            self._set_params(base_params)
             return self.map[c_new]
 
     def step_pred(self, x) -> int:
@@ -212,9 +259,17 @@ def step_pred(self, x) -> int:
         c_ = int(np.argmax(T))
         return self.map[c_]
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return self.base_module.get_cluster_centers()
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
-        undefined function for visualizing the bounds of each cluster
+        function for visualizing the bounds of each cluster
 
         Parameters:
         - ax: figure axes
diff --git a/artlib/elementary/EllipsoidART.py b/artlib/elementary/EllipsoidART.py
index ecce7ed..fd0ccd2 100644
--- a/artlib/elementary/EllipsoidART.py
+++ b/artlib/elementary/EllipsoidART.py
@@ -10,7 +10,7 @@
 International Society for Optics and Photonics. doi:10.1117/12.421180.
 """
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, List
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 from artlib.common.utils import l2norm2
@@ -117,22 +117,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
 
         return 1 - (radius + max(radius, dist))/params["r_hat"], cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params=params, cache=cache)
-        return M >= params["rho"], cache
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -196,6 +180,14 @@ def get_2d_ellipsoids(self) -> list[tuple]:
 
         return ellipsoids
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return [w[:self.dim_] for w in self.W]
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         undefined function for visualizing the bounds of each cluster
diff --git a/artlib/elementary/FuzzyART.py b/artlib/elementary/FuzzyART.py
index 860f619..7ecef9b 100644
--- a/artlib/elementary/FuzzyART.py
+++ b/artlib/elementary/FuzzyART.py
@@ -4,10 +4,10 @@
 Neural Networks, 4, 759 – 771. doi:10.1016/0893-6080(91)90056-B.
 """
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, List
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
-from artlib.common.utils import normalize, compliment_code, l1norm, fuzzy_and
+from artlib.common.utils import normalize, compliment_code, l1norm, fuzzy_and, de_compliment_code
 
 
 def get_bounding_box(w: np.ndarray, n: Optional[int] = None) -> tuple[list[int], list[int]]:
@@ -56,8 +56,7 @@ def __init__(self, rho: float, alpha: float, beta: float):
         }
         super().__init__(params)
 
-    @staticmethod
-    def prepare_data(X: np.ndarray) -> np.ndarray:
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
         """
         prepare data for clustering
 
@@ -67,10 +66,23 @@ def prepare_data(X: np.ndarray) -> np.ndarray:
         Returns:
             normalized and compliment coded data
         """
-        normalized = normalize(X)
+        normalized, self.d_max_, self.d_min_ = normalize(X, self.d_max_, self.d_min_)
         cc_data = compliment_code(normalized)
         return cc_data
 
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        out = de_compliment_code(X)
+        return super(FuzzyART, self).restore_data(out)
+
     @staticmethod
     def validate_params(params: dict):
         """
@@ -149,22 +161,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         """
         return l1norm(fuzzy_and(i, w)) / self.dim_original, cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params)
-        return M >= params["rho"], cache
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -203,6 +199,31 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
     def get_bounding_boxes(self, n: Optional[int] = None):
         return list(map(lambda w: get_bounding_box(w, n=n), self.W))
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        centers = []
+        for w in self.W:
+            ref_points, widths = get_bounding_box(w,None)
+            centers.append(np.array(ref_points)+0.5*np.array(widths))
+        return centers
+
+    def shrink_clusters(self, shrink_ratio: float = 0.1):
+        new_W = []
+        dim = len(self.W[0])//2
+        for w in self.W:
+            new_w = np.copy(w)
+            widths = (1-w[dim:]) - w[:dim]
+            new_w[:dim] += widths*shrink_ratio
+            new_w[dim:] += widths*shrink_ratio
+            new_W.append(new_w)
+        self.W = new_W
+        return self
+
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         undefined function for visualizing the bounds of each cluster
diff --git a/artlib/elementary/GaussianART.py b/artlib/elementary/GaussianART.py
index 8dd7fa1..ccf2b9f 100644
--- a/artlib/elementary/GaussianART.py
+++ b/artlib/elementary/GaussianART.py
@@ -5,7 +5,8 @@
 """
 
 import numpy as np
-from typing import Optional, Iterable
+from decimal import Decimal
+from typing import Optional, Iterable, List
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 from artlib.common.visualization import plot_gaussian_contours_fading
@@ -13,20 +14,22 @@
 
 class GaussianART(BaseART):
     # implementation of GaussianART
-    pi2 = np.pi*2
-    def __init__(self, rho: float, sigma_init: np.ndarray):
+    def __init__(self, rho: float, sigma_init: np.ndarray, alpha: float = 1e-10):
         """
         Parameters:
         - rho: vigilance parameter
         - sigma_init: initial estimate of the diagonal std
+        - alpha: used to prevent division by zero errors
 
         """
         params = {
             "rho": rho,
             "sigma_init": sigma_init,
+            "alpha": alpha
         }
         super().__init__(params)
 
+
     @staticmethod
     def validate_params(params: dict):
         """
@@ -38,7 +41,9 @@ def validate_params(params: dict):
         """
         assert "rho" in params
         assert "sigma_init" in params
-        assert 1.0 >= params["rho"] > 0.
+        assert "alpha" in params
+        assert 1.0 >= params["rho"] >= 0.
+        assert params["alpha"] > 0.
         assert isinstance(params["rho"], float)
         assert isinstance(params["sigma_init"], np.ndarray)
 
@@ -57,15 +62,19 @@ def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict) -> tuple[f
 
         """
         mean = w[:self.dim_]
-        sigma = w[self.dim_:-1]
+        # sigma = w[self.dim_:2*self.dim]
+        inv_sig = w[2*self.dim_:3*self.dim_]
+        sqrt_det_sig = w[-2]
         n = w[-1]
-        sig = np.diag(np.multiply(sigma,sigma))
+
         dist = mean-i
-        exp_dist_sig_dist = np.exp(-0.5*np.matmul(dist.T, np.matmul(np.linalg.inv(sig), dist)))
+        exp_dist_sig_dist = np.exp(-0.5 * np.dot(dist, np.multiply(inv_sig, dist)))
+
         cache = {
             "exp_dist_sig_dist": exp_dist_sig_dist
         }
-        p_i_cj = exp_dist_sig_dist/np.sqrt((self.pi2**self.dim_)*np.linalg.det(sig))
+        # ignore the (2*pi)^d term as that is constant
+        p_i_cj = exp_dist_sig_dist/(params["alpha"]+sqrt_det_sig)
         p_cj = n/np.sum(w_[-1] for w_ in self.W)
 
         activation = p_i_cj*p_cj
@@ -93,23 +102,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         return exp_dist_sig_dist, cache
 
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params=params, cache=cache)
-        return M >= params["rho"], cache
-
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -126,14 +118,18 @@ def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dic
 
         """
         mean = w[:self.dim_]
-        sigma = w[self.dim_:-1]
+        sigma = w[self.dim_:2*self.dim_]
         n = w[-1]
 
         n_new = n+1
         mean_new = (1-(1/n_new))*mean + (1/n_new)*i
         sigma_new = np.sqrt((1-(1/n_new))*np.multiply(sigma, sigma) + (1/n_new)*((mean_new - i)**2))
 
-        return np.concatenate([mean_new, sigma_new, [n_new]])
+        sigma2 = np.multiply(sigma_new, sigma_new)
+        inv_sig = 1 / sigma2
+        det_sig = np.sqrt(np.prod(sigma2))
+
+        return np.concatenate([mean_new, sigma_new, inv_sig, [det_sig], [n_new]])
 
 
     def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
@@ -149,7 +145,18 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
             updated cluster weight
 
         """
-        return np.concatenate([i, params["sigma_init"], [1.]])
+        sigma2 = np.multiply(params["sigma_init"], params["sigma_init"])
+        inv_sig_init = 1 / sigma2
+        det_sig_init = np.sqrt(np.prod(sigma2))
+        return np.concatenate([i, params["sigma_init"], inv_sig_init, [det_sig_init], [1.]])
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return [w[:self.dim_] for w in self.W]
 
 
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
diff --git a/artlib/elementary/HypersphereART.py b/artlib/elementary/HypersphereART.py
index 30e95de..5716c38 100644
--- a/artlib/elementary/HypersphereART.py
+++ b/artlib/elementary/HypersphereART.py
@@ -5,7 +5,7 @@
 (pp. 59–64). volume 6. doi:10.1109/IJCNN.2000.859373.
 """
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, List
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 from artlib.common.utils import l2norm2
@@ -103,23 +103,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         return 1 - (max(radius, max_radius)/params["r_hat"]), cache
 
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params=params, cache=cache)
-        return M >= params["rho"], cache
-
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -163,6 +146,14 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
         """
         return np.concatenate([i, [0.]])
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return [w[:-1] for w in self.W]
+
 
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
diff --git a/artlib/elementary/QuadraticNeuronART.py b/artlib/elementary/QuadraticNeuronART.py
index c99a1ee..04b42ca 100644
--- a/artlib/elementary/QuadraticNeuronART.py
+++ b/artlib/elementary/QuadraticNeuronART.py
@@ -9,16 +9,12 @@
 """
 
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, List
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
-from artlib.common.utils import normalize, l2norm2
+from artlib.common.utils import l2norm2
 from artlib.common.visualization import plot_weight_matrix_as_ellipse
 
-def prepare_data(data: np.ndarray) -> np.ndarray:
-    normalized = normalize(data)
-    return normalized
-
 
 class QuadraticNeuronART(BaseART):
     # implementation of QuadraticNeuronART
@@ -114,23 +110,6 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
             raise ValueError("No cache provided")
         return cache["activation"], cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
-        """
-        get the binary match criterion of the cluster
-
-        Parameters:
-        - i: data sample
-        - w: cluster weight / info
-        - params: dict containing parameters for the algorithm
-        - cache: dict containing values cached from previous calculations
-
-        Returns:
-            cluster match criterion binary, cache used for later processing
-
-        """
-        M, cache = self.match_criterion(i, w, params, cache)
-        return M >= params["rho"], cache
-
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
         get the updated cluster weight
@@ -177,6 +156,15 @@ def new_weight(self, i: np.ndarray, params: dict) -> np.ndarray:
         w_new = np.identity(self.dim_)
         return np.concatenate([w_new.flatten(), i, [params["s_init"]]])
 
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        dim2 = self.dim_ * self.dim_
+        return [w[dim2:-1] for w in self.W]
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         undefined function for visualizing the bounds of each cluster
diff --git a/artlib/experimental/ConvexHullART.py b/artlib/experimental/ConvexHullART.py
new file mode 100644
index 0000000..84b2092
--- /dev/null
+++ b/artlib/experimental/ConvexHullART.py
@@ -0,0 +1,327 @@
+import numpy as np
+from matplotlib.axes import Axes
+from copy import deepcopy
+from typing import Optional, Iterable, List, Union
+from scipy.spatial import ConvexHull
+
+from artlib.common.BaseART import BaseART
+from artlib.experimental.merging import merge_objects
+
+
+def plot_convex_polygon(vertices: np.ndarray, ax: Axes, line_color: str = 'b', line_width: float = 1.0):
+    """
+    Plots a convex polygon given its vertices using Matplotlib.
+
+    Parameters:
+    - vertices: A list of vertices representing a convex polygon.
+    - ax: A matplotlib Axes object to plot on. If None, creates a new figure and axes.
+    - line_color: The color of the polygon lines.
+    - line_width: The width of the polygon lines.
+    """
+    vertices = np.array(vertices)
+    # Close the polygon by appending the first vertex at the end
+    vertices = np.vstack([vertices, vertices[0]])
+
+    ax.plot(vertices[:, 0], vertices[:, 1], linestyle='-', color=line_color, linewidth=line_width)
+
+
+def volume_of_simplex(vertices):
+    """
+    Calculates the n-dimensional volume of a simplex defined by its vertices.
+
+    Parameters:
+    - vertices: An (n+1) x n array representing the coordinates of the simplex vertices.
+
+    Returns:
+    - Volume of the simplex.
+    """
+    vertices = np.asarray(vertices)
+    # Subtract the first vertex from all vertices to form a matrix
+    matrix = vertices[1:] - vertices[0]
+    # Calculate the absolute value of the determinant divided by factorial(n) for volume
+    return np.abs(np.linalg.det(matrix)) / np.math.factorial(len(vertices) - 1)
+
+
+def minimum_distance(a1, a2):
+    def point_to_point_distance(P, Q):
+        """Calculate the Euclidean distance between two points P and Q."""
+        return np.linalg.norm(P - Q)
+
+    def point_to_line_segment_distance(P, Q1, Q2):
+        """Calculate the minimum distance between point P and line segment Q1-Q2."""
+        line_vec = Q2 - Q1
+        point_vec = P - Q1
+        line_len = np.dot(line_vec, line_vec)
+
+        if line_len == 0:  # Q1 and Q2 are the same point
+            return np.linalg.norm(P - Q1)
+
+        t = max(0, min(1, np.dot(point_vec, line_vec) / line_len))
+        projection = Q1 + t * line_vec
+        return np.linalg.norm(P - projection)
+
+    def line_segment_to_line_segment_distance(A1, A2, B1, B2):
+        """Calculate the minimum distance between two line segments A1-A2 and B1-B2."""
+        distances = [
+            point_to_line_segment_distance(A1, B1, B2),
+            point_to_line_segment_distance(A2, B1, B2),
+            point_to_line_segment_distance(B1, A1, A2),
+            point_to_line_segment_distance(B2, A1, A2),
+        ]
+        return min(distances)
+
+    # Determine the cases and compute the distance
+    if a1.shape[0] == 1 and a2.shape[0] == 1:
+        # Case 1: Point to point
+        return point_to_point_distance(a1[0], a2[0])
+
+    elif a1.shape[0] == 1 and a2.shape[0] == 2:
+        # Case 2: Point to line segment
+        return point_to_line_segment_distance(a1[0], a2[0], a2[1])
+
+    elif a1.shape[0] == 2 and a2.shape[0] == 1:
+        # Case 3: Line segment to point
+        return point_to_line_segment_distance(a2[0], a1[0], a1[1])
+
+    elif a1.shape[0] == 2 and a2.shape[0] == 2:
+        # Case 4: Line segment to line segment
+        return line_segment_to_line_segment_distance(a1[0], a1[1], a2[0], a2[1])
+    else:
+        raise RuntimeError("INVALID DISTANCE CALCULATION")
+
+
+
+class PseudoConvexHull:
+    def __init__(self, points: np.ndarray):
+        self.points = points
+
+    @property
+    def vertices(self):
+        return np.array([i for i in range(self.points.shape[0])])
+
+    def add_points(self, points):
+        self.points = np.vstack([self.points, points])
+
+
+HullTypes = Union[ConvexHull, PseudoConvexHull]
+
+
+def centroid_of_convex_hull(hull: HullTypes):
+    """
+    Finds the centroid of the volume of a convex hull in n-dimensional space.
+
+    Parameters:
+    - vertices: An array of shape (m, n), where m is the number of vertices and n is the dimension.
+
+    Returns:
+    - Centroid coordinates as a numpy array of length n.
+    """
+    hull_vertices = hull.points[hull.vertices]
+
+    centroid = np.zeros(hull_vertices.shape[1])
+    total_volume = 0
+
+    if isinstance(hull, PseudoConvexHull):
+        return np.mean(hull.points, axis=0)
+
+    for simplex in hull.simplices:
+        simplex_vertices = hull_vertices[simplex]
+        simplex_centroid = np.mean(simplex_vertices, axis=0)
+        simplex_volume = volume_of_simplex(simplex_vertices)
+
+        centroid += simplex_centroid * simplex_volume
+        total_volume += simplex_volume
+
+    centroid /= total_volume
+    return centroid
+
+
+class ConvexHullART(BaseART):
+    def __init__(self, rho: float, merge_rho: float):
+        """
+        Parameters:
+        - rho: vigilance parameter
+
+        """
+        params = {
+            "rho": rho,
+            "merge_rho": merge_rho
+        }
+        super().__init__(params)
+
+    @staticmethod
+    def validate_params(params: dict):
+        """
+        validate clustering parameters
+
+        Parameters:
+        - params: dict containing parameters for the algorithm
+
+        """
+        assert "rho" in params
+        assert params["rho"] >= 0.
+        assert isinstance(params["rho"], float)
+
+    def category_choice(self, i: np.ndarray, w: HullTypes, params: dict) -> tuple[float, Optional[dict]]:
+        """
+        get the activation of the cluster
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+
+        Returns:
+            cluster activation, cache used for later processing
+
+        """
+        if isinstance(w, PseudoConvexHull):
+            d = minimum_distance(i.reshape((1,-1)), w.points)
+            activation = 1. - d**len(i)
+            if w.points.shape[0] == 1:
+                new_w = deepcopy(w)
+                new_w.add_points(i.reshape((1,-1)))
+            else:
+                new_points = np.vstack([w.points[w.vertices,:], i.reshape((1,-1))])
+                new_w = ConvexHull(new_points, incremental=True)
+        else:
+            new_w = ConvexHull(w.points[w.vertices,:], incremental=True)
+            new_w.add_points(i.reshape((1,-1)))
+            activation = w.area / new_w.area
+
+        cache = {"new_w": new_w, "activation": activation}
+        return activation, cache
+
+    def match_criterion(self, i: np.ndarray, w: HullTypes, params: dict, cache: Optional[dict] = None) -> tuple[float, dict]:
+        """
+        get the match criterion of the cluster
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+        - cache: dict containing values cached from previous calculations
+
+        Returns:
+            cluster match criterion, cache used for later processing
+
+        """
+        return cache["activation"], cache
+
+
+    def update(self, i: np.ndarray, w: HullTypes, params: dict, cache: Optional[dict] = None) -> HullTypes:
+        """
+        get the updated cluster weight
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+        - cache: dict containing values cached from previous calculations
+
+        Returns:
+            updated cluster weight, cache used for later processing
+
+        """
+        return cache["new_w"]
+
+    def new_weight(self, i: np.ndarray, params: dict) -> HullTypes:
+        """
+        generate a new cluster weight
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+
+        Returns:
+            updated cluster weight
+
+        """
+        new_w = PseudoConvexHull(i.reshape((1,-1)))
+        return new_w
+
+    def merge_clusters(self):
+        def can_merge(w1, w2):
+            combined_points = np.vstack([w1.points[w1.vertices,:], w2.points[w2.vertices,:]])
+
+            if isinstance(w1, PseudoConvexHull) and isinstance(w2, PseudoConvexHull):
+                d = minimum_distance(w1.points, w2.points)
+                activation = 1.0 - d**w1.points.shape[1]
+            else:
+                new_w = ConvexHull(combined_points)
+                if isinstance(w1, ConvexHull):
+                    a1 = w1.area / new_w.area
+                else:
+                    a1 = np.nan
+                if isinstance(w2, ConvexHull):
+                    a2 = w2.area / new_w.area
+                else:
+                    a2 = np.nan
+                activation = np.max([a1,a2])
+
+            if activation > self.params["merge_rho"]:
+                return True
+            else:
+                return False
+
+        merges = merge_objects(self.W, can_merge)
+
+        new_W = []
+        new_sample_counter = []
+        new_labels = np.copy(self.labels_)
+        for i in range(len(merges)):
+            new_labels[np.isin(self.labels_, merges[i])] = i
+            new_sample_counter.append(sum(self.weight_sample_counter_[j] for j in merges[i]))
+            if len(merges[i]) > 1:
+                new_points = np.vstack([self.W[j].points for j in merges[i]])
+                if new_points.shape[0] > 2:
+                    new_W.append(ConvexHull(new_points, incremental=True))
+                else:
+                    new_W.append(PseudoConvexHull(new_points))
+            else:
+                new_W.append(self.W[merges[i][0]])
+        self.W = new_W
+        self.weight_sample_counter_ = new_sample_counter
+        self.labels_ = new_labels
+
+    def post_fit(self, X: np.ndarray):
+        """
+        function called after fit. Useful for cluster pruning
+
+        Parameters:
+        - X: data set
+
+        """
+        self.merge_clusters()
+
+
+
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        centers = []
+        for w in self.W:
+            centers.append(centroid_of_convex_hull(w))
+        return centers
+
+    def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
+        """
+        undefined function for visualizing the bounds of each cluster
+
+        Parameters:
+        - ax: figure axes
+        - colors: colors to use for each cluster
+        - linewidth: width of boundary line
+
+        """
+        for c, w in zip(colors, self.W):
+            if isinstance(w, ConvexHull):
+                vertices = w.points[w.vertices,:2]
+            else:
+                vertices = w.points[:, :2]
+            plot_convex_polygon(vertices, ax, line_width=linewidth, line_color=c)
diff --git a/artlib/experimental/SeqART.py b/artlib/experimental/SeqART.py
new file mode 100644
index 0000000..9900afe
--- /dev/null
+++ b/artlib/experimental/SeqART.py
@@ -0,0 +1,215 @@
+import numpy as np
+from typing import Optional, Callable
+from artlib import BaseART
+import operator
+import re
+
+def compress_dashes(input_string):
+    return re.sub('-+', '-', input_string)
+
+def arr2seq(x):
+    return "".join([str(i_) for i_ in x])
+
+def needleman_wunsch(seq1, seq2, match_score=1, gap_cost=-1, mismatch_cost=-1):
+    m, n = len(seq1), len(seq2)
+
+    # Initialize the scoring matrix
+    score_matrix = np.zeros((m + 1, n + 1), dtype=int)
+
+    # Initialize the gap penalties for the first row and column
+    for i in range(1, m + 1):
+        score_matrix[i][0] = score_matrix[i-1][0] + gap_cost
+    for j in range(1, n + 1):
+        score_matrix[0][j] = score_matrix[0][j-1] + gap_cost
+
+    # Fill in the scoring matrix
+    for i in range(1, m + 1):
+        for j in range(1, n + 1):
+            match = score_matrix[i-1][j-1] + (match_score if seq1[i-1] == seq2[j-1] else mismatch_cost)
+            delete = score_matrix[i-1][j] + gap_cost
+            insert = score_matrix[i][j-1] + gap_cost
+            score_matrix[i][j] = max(match, delete, insert)
+
+    # Traceback to get the alignment
+    align1, align2 = '', ''
+    i, j = m, n
+
+    while i > 0 and j > 0:
+        current_score = score_matrix[i][j]
+        diagonal_score = score_matrix[i-1][j-1]
+        up_score = score_matrix[i-1][j]
+        left_score = score_matrix[i][j-1]
+
+        if current_score == diagonal_score + (match_score if seq1[i-1] == seq2[j-1] else mismatch_cost):
+            align1 += seq1[i-1]
+            align2 += seq2[j-1]
+            i -= 1
+            j -= 1
+        elif current_score == up_score + gap_cost:
+            align1 += seq1[i-1]
+            align2 += '-'
+            i -= 1
+        elif current_score == left_score + gap_cost:
+            align1 += '-'
+            align2 += seq2[j-1]
+            j -= 1
+
+    # Fill the remaining sequence if the above while loop finishes
+    while i > 0:
+        align1 += seq1[i-1]
+        align2 += '-'
+        i -= 1
+
+    while j > 0:
+        align1 += '-'
+        align2 += seq2[j-1]
+        j -= 1
+
+    # The alignments are built from the end to the beginning, so we need to reverse them
+    align1 = align1[::-1]
+    align2 = align2[::-1]
+    alignment = ''.join([a if a == b else '-' for a, b in zip(align1, align2)])
+    l = max(len(seq1), len(seq2))
+    return alignment, float(score_matrix[m][n])/l
+
+
+def prepare_data(data: np.ndarray) -> np.ndarray:
+    return data
+
+
+class SeqART(BaseART):
+    # template for ART module
+    def __init__(self, rho: float, metric: Callable = needleman_wunsch):
+        """
+        Parameters:
+        - rho: vigilance parameter
+        - metric: allignment function. Should be in the format alignment, score = metric(seq_a, seq_b)
+
+        """
+        params = {
+            "rho": rho,
+        }
+        self.metric = metric
+        super().__init__(params)
+    @staticmethod
+    def validate_params(params: dict):
+        """
+        validate clustering parameters
+
+        Parameters:
+        - params: dict containing parameters for the algorithm
+
+        """
+        assert "rho" in params
+        assert isinstance(params["rho"], float)
+
+    def validate_data(self, X: np.ndarray):
+        """
+        validates the data prior to clustering
+
+        Parameters:
+        - X: data set
+
+        """
+        pass
+
+    def check_dimensions(self, X: np.ndarray):
+        """
+        check the data has the correct dimensions
+
+        Parameters:
+        - X: data set
+
+        """
+        pass
+
+
+    def category_choice(self, i: str, w: str, params: dict) -> tuple[float, Optional[dict]]:
+        """
+        get the activation of the cluster
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+
+        Returns:
+            cluster activation, cache used for later processing
+
+        """
+        alignment, score = self.metric(arr2seq(i), w)
+        cache = {'alignment': alignment, 'score': score}
+        return score, cache
+
+    def match_criterion(self, i: str, w: str, params: dict, cache: Optional[dict] = None) -> tuple[float, dict]:
+        """
+        get the match criterion of the cluster
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+        - cache: dict containing values cached from previous calculations
+
+        Returns:
+            cluster match criterion, cache used for later processing
+
+        """
+        # _, M = self.metric(cache['alignment'], w)
+
+        return cache['score'], cache
+
+    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, op: Callable = operator.ge) -> tuple[bool, dict]:
+        """
+        get the binary match criterion of the cluster
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+        - cache: dict containing values cached from previous calculations
+
+        Returns:
+            cluster match criterion binary, cache used for later processing
+
+        """
+        M, cache = self.match_criterion(arr2seq(i), w, params, cache)
+        M_bin = op(M, params["rho"])
+        if cache is None:
+            cache = dict()
+        cache["match_criterion"] = M
+        cache["match_criterion_bin"] = M_bin
+
+        return M_bin, cache
+
+    def update(self, i: str, w: str, params: dict, cache: Optional[dict] = None) -> str:
+        """
+        get the updated cluster weight
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+        - cache: dict containing values cached from previous calculations
+
+        Returns:
+            updated cluster weight, cache used for later processing
+
+        """
+        # print(cache['alignment'])
+        return compress_dashes(cache['alignment'])
+
+    def new_weight(self, i: str, params: dict) -> str:
+        """
+        generate a new cluster weight
+
+        Parameters:
+        - i: data sample
+        - w: cluster weight / info
+        - params: dict containing parameters for the algorithm
+
+        Returns:
+            updated cluster weight
+
+        """
+        return arr2seq(i)
\ No newline at end of file
diff --git a/artlib/experimental/__init__.py b/artlib/experimental/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/artlib/experimental/merging.py b/artlib/experimental/merging.py
new file mode 100644
index 0000000..615a9d0
--- /dev/null
+++ b/artlib/experimental/merging.py
@@ -0,0 +1,44 @@
+def find(parent, i):
+    if parent[i] == i:
+        return i
+    else:
+        parent[i] = find(parent, parent[i])
+        return parent[i]
+
+
+def union(parent, rank, x, y):
+    root_x = find(parent, x)
+    root_y = find(parent, y)
+
+    if root_x != root_y:
+        # Union by rank to keep the tree flat
+        if rank[root_x] > rank[root_y]:
+            parent[root_y] = root_x
+        elif rank[root_x] < rank[root_y]:
+            parent[root_x] = root_y
+        else:
+            parent[root_y] = root_x
+            rank[root_x] += 1
+
+
+def merge_objects(objects, can_merge):
+    # Initialize Union-Find structure
+    n = len(objects)
+    parent = list(range(n))
+    rank = [0] * n
+
+    # Check each pair of objects
+    for i in range(n):
+        for j in range(i + 1, n):
+            if can_merge(objects[i], objects[j]):
+                union(parent, rank, i, j)
+
+    # Collect merged groups by indices
+    groups = {}
+    for i in range(n):
+        root = find(parent, i)
+        if root not in groups:
+            groups[root] = []
+        groups[root].append(i)
+
+    return list(groups.values())
diff --git a/artlib/fusion/FusionART.py b/artlib/fusion/FusionART.py
index ebb3678..0728680 100644
--- a/artlib/fusion/FusionART.py
+++ b/artlib/fusion/FusionART.py
@@ -7,8 +7,11 @@
 doi:10.1007/ 978-3-540-72383-7_128.
 """
 import numpy as np
-from typing import Optional, Union
+from typing import Optional, Union, Callable, List, Literal, Tuple
+from copy import deepcopy
 from artlib.common.BaseART import BaseART
+from sklearn.utils.validation import check_is_fitted
+import operator
 
 def get_channel_position_tuples(channel_dims: list[int]) -> list[tuple[int, int]]:
     positions = []
@@ -54,6 +57,9 @@ def get_params(self, deep: bool = True) -> dict:
             out[f"module_{i}"] = module
         return out
 
+    @property
+    def n_clusters(self) -> int:
+        return self.modules[0].n_clusters
 
     @property
     def W(self):
@@ -115,7 +121,7 @@ def check_dimensions(self, X: np.ndarray):
         """
         assert X.shape[1] == self.dim_, "Invalid data shape"
 
-    def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict) -> tuple[float, Optional[dict]]:
+    def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict, skip_channels: List[int] = []) -> tuple[float, Optional[dict]]:
         """
         get the activation of the cluster
 
@@ -135,6 +141,8 @@ def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict) -> tuple[f
                     w[self._channel_indices[k][0]:self._channel_indices[k][1]],
                     self.modules[k].params
                 )
+                if k not in skip_channels
+                else (1., dict())
                 for k in range(self.n)
             ]
         )
@@ -142,7 +150,7 @@ def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict) -> tuple[f
         activation = sum([a*self.params["gamma_values"][k] for k, a in enumerate(activations)])
         return activation, cache
 
-    def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[list[float], dict]:
+    def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, skip_channels: List[int] = []) -> tuple[list[float], dict]:
         if cache is None:
             raise ValueError("No cache provided")
         M, caches = zip(
@@ -153,13 +161,15 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
                     self.modules[k].params,
                     cache[k]
                 )
+                if k not in skip_channels
+                else (np.inf, {"match_criterion": np.inf})
                 for k in range(self.n)
             ]
         )
         cache = {k: cache_k for k, cache_k in enumerate(caches)}
         return M, cache
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
+    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, skip_channels: List[int] = [], op: Callable = operator.ge) -> tuple[bool, dict]:
         """
         get the binary match criterion of the cluster
 
@@ -181,14 +191,110 @@ def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache:
                     i[self._channel_indices[k][0]:self._channel_indices[k][1]],
                     w[self._channel_indices[k][0]:self._channel_indices[k][1]],
                     self.modules[k].params,
-                    cache[k]
+                    cache[k],
+                    op
                 )
+                if k not in skip_channels
+                else (True, {"match_criterion": np.inf})
                 for k in range(self.n)
             ]
         )
         cache = {k: cache_k for k, cache_k in enumerate(caches)}
         return all(M_bin), cache
 
+
+    def _match_tracking(self, cache: List[dict], epsilon: float, params: List[dict], method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        keep_searching = []
+        for i in range(len(cache)):
+            if cache[i]["match_criterion_bin"]:
+                keep_searching_i = self.modules[i]._match_tracking(cache[i], epsilon, params[i], method)
+                keep_searching.append(keep_searching_i)
+            else:
+                keep_searching.append(True)
+        return all(keep_searching)
+
+
+    def _set_params(self, new_params):
+        for i in range(self.n):
+            self.modules[i].params = new_params[i]
+
+    def _deep_copy_params(self):
+        return {i: deepcopy(module.params) for i, module in enumerate(self.modules)}
+
+
+    def partial_fit(self, X: np.ndarray, match_reset_func: Optional[Callable] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
+        """
+        iteratively fit the model to the data
+
+        Parameters:
+        - X: data set
+        - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
+            Permits external factors to influence cluster creation.
+            Returns True if the cluster is valid for the sample, False otherwise
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
+        """
+
+        self.validate_data(X)
+        self.check_dimensions(X)
+        self.is_fitted_ =  True
+
+        if not hasattr(self.modules[0], 'W'):
+            self.W: list[np.ndarray] = []
+            self.labels_ = np.zeros((X.shape[0], ), dtype=int)
+            j = 0
+        else:
+            j = len(self.labels_)
+            self.labels_ = np.pad(self.labels_, [(0, X.shape[0])], mode='constant')
+        for i, x in enumerate(X):
+            c = self.step_fit(x, match_reset_func=match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
+            self.labels_[i+j] = c
+        return self
+
+    def step_pred(self, x, skip_channels: List[int] = []) -> int:
+        """
+        predict the label for a single sample
+
+        Parameters:
+        - x: data sample
+
+        Returns:
+            cluster label of the input sample
+
+        """
+        assert len(self.W) >= 0, "ART module is not fit."
+
+        T, _ = zip(*[self.category_choice(x, w, params=self.params, skip_channels=skip_channels) for w in self.W])
+        c_ = int(np.argmax(T))
+        return c_
+
+    def predict(self, X: np.ndarray, skip_channels: List[int] = []) -> np.ndarray:
+        """
+        predict labels for the data
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            labels for the data
+
+        """
+
+        check_is_fitted(self)
+        self.validate_data(X)
+        self.check_dimensions(X)
+
+        y = np.zeros((X.shape[0],), dtype=int)
+        for i, x in enumerate(X):
+            c = self.step_pred(x, skip_channels=skip_channels)
+            y[i] = c
+        return y
+
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
         get the updated cluster weight
@@ -260,3 +366,50 @@ def set_weight(self, idx: int, new_w: np.ndarray):
         for k in range(self.n):
             new_w_k = new_w[self._channel_indices[k][0]:self._channel_indices[k][1]]
             self.modules[k].set_weight(idx, new_w_k)
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        centers_ = [module.get_cluster_centers() for module in self.modules]
+        centers = [
+            np.concatenate(
+                [
+                    centers_[k][i]
+                    for k in range(self.n)
+                ]
+            )
+            for i
+            in range(self.n_clusters)
+        ]
+        return centers
+
+    def get_channel_centers(self, channel: int):
+        return self.modules[channel].get_cluster_centers()
+
+    def predict_regression(self, X: np.ndarray, target_channels: List[int] = [-1]) -> Union[np.ndarray, List[np.ndarray]]:
+        target_channels = [self.n+k if k < 0 else k for k in target_channels]
+        C = self.predict(X, skip_channels=target_channels)
+        centers = [self.get_channel_centers(k) for k in target_channels]
+        if len(target_channels) == 1:
+            return np.array([centers[0][c] for c in C])
+        else:
+            return [np.array([centers[k][c] for c in C]) for k in target_channels]
+
+    def join_channel_data(self, channel_data: List[np.ndarray], skip_channels: List[int] = []) -> np.ndarray:
+        skip_channels = [self.n+k if k < 0 else k for k in skip_channels]
+        n_samples = channel_data[0].shape[0]
+
+        formatted_channel_data = []
+        i = 0
+        for k in range(self.n):
+            if k not in skip_channels:
+                formatted_channel_data.append(channel_data[i])
+                i += 1
+            else:
+                formatted_channel_data.append(0.5*np.ones((n_samples, self._channel_indices[k][1]-self._channel_indices[k][0])))
+
+        X = np.hstack(formatted_channel_data)
+        return X
\ No newline at end of file
diff --git a/artlib/hierarchical/DeepARTMAP.py b/artlib/hierarchical/DeepARTMAP.py
index 0c5fcb7..61694c8 100644
--- a/artlib/hierarchical/DeepARTMAP.py
+++ b/artlib/hierarchical/DeepARTMAP.py
@@ -5,7 +5,7 @@
 """
 import numpy as np
 from sklearn.base import BaseEstimator, ClassifierMixin, ClusterMixin
-from typing import Optional, cast, Union
+from typing import Optional, cast, Union, Literal, Tuple
 from collections import defaultdict
 from artlib.common.BaseART import BaseART
 from artlib.common.BaseARTMAP import BaseARTMAP
@@ -149,8 +149,32 @@ def validate_data(
             n = X[0].shape[0]
         assert all(x.shape[0] == n for x in X), "Inconsistent sample number in input matrices"
 
+    def prepare_data(self,  X: list[np.ndarray], y: Optional[np.ndarray] = None) ->Tuple[list[np.ndarray], Optional[np.ndarray]]:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            prepared data
+        """
+        return [self.modules[i].prepare_data(X[i]) for i in range(self.n_modules)], y
+
+    def restore_data(self,  X: list[np.ndarray], y: Optional[np.ndarray] = None) ->Tuple[list[np.ndarray], Optional[np.ndarray]]:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            prepared data
+        """
+        return [self.modules[i].restore_data(X[i]) for i in range(self.n_modules)], y
+
 
-    def fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None, max_iter=1):
+    def fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None, max_iter=1, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
         """
         Fit the model to the data
 
@@ -158,34 +182,46 @@ def fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None, max_iter=1):
         - X: list of deep datasets
         - y: optional labels
         - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         self.validate_data(X, y)
         if y is not None:
             self.is_supervised = True
             self.layers = [SimpleARTMAP(self.modules[i]) for i in range(self.n_modules)]
-            self.layers[0] = self.layers[0].fit(X[0], y, max_iter=max_iter)
+            self.layers[0] = self.layers[0].fit(X[0], y, max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
         else:
             self.is_supervised = False
             assert self.n_modules >= 2, "Must provide at least two ART modules when providing cluster labels"
             self.layers = cast(list[BaseARTMAP], [ARTMAP(self.modules[1], self.modules[0])]) + \
                            cast(list[BaseARTMAP], [SimpleARTMAP(self.modules[i]) for i in range(2, self.n_modules)])
-            self.layers[0] = self.layers[0].fit(X[1], X[0], max_iter=max_iter)
+            self.layers[0] = self.layers[0].fit(X[1], X[0], max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
 
         for art_i in range(1, self.n_layers):
             y_i = self.layers[art_i-1].labels_a
-            self.layers[art_i] = self.layers[art_i].fit(X[art_i], y_i, max_iter=max_iter)
+            self.layers[art_i] = self.layers[art_i].fit(X[art_i], y_i, max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
 
         return self
 
 
-    def partial_fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None):
+    def partial_fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
         """
         Partial fit the model to the data
 
         Parameters:
         - X: list of deep datasets
         - y: optional labels
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         self.validate_data(X, y)
@@ -194,7 +230,7 @@ def partial_fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None):
                 self.is_supervised = True
                 self.layers = [SimpleARTMAP(self.modules[i]) for i in range(self.n_modules)]
             assert self.is_supervised, "Labels were previously provided. Must continue to provide labels for partial fit."
-            self.layers[0] = self.layers[0].partial_fit(X[0], y)
+            self.layers[0] = self.layers[0].partial_fit(X[0], y, match_reset_method=match_reset_method, epsilon=epsilon)
             x_i = 1
         else:
             if len(self.layers) == 0:
@@ -203,14 +239,14 @@ def partial_fit(self, X: list[np.ndarray], y: Optional[np.ndarray] = None):
                 self.layers = cast(list[BaseARTMAP], [ARTMAP(self.modules[1], self.modules[0])]) + \
                                cast(list[BaseARTMAP], [SimpleARTMAP(self.modules[i]) for i in range(2, self.n_modules)])
             assert not self.is_supervised, "Labels were not previously provided. Do not provide labels to continue partial fit."
-            self.layers[0] = self.layers[0].partial_fit(X[1], X[0])
+            self.layers[0] = self.layers[0].partial_fit(X[1], X[0], match_reset_method=match_reset_method, epsilon=epsilon)
             x_i = 2
 
         n_samples = X[0].shape[0]
 
         for art_i in range(1, self.n_modules):
             y_i = self.layers[art_i-1].labels_a[-n_samples:]
-            self.layers[art_i] = self.layers[art_i].partial_fit(X[x_i], y_i)
+            self.layers[art_i] = self.layers[art_i].partial_fit(X[x_i], y_i, match_reset_method=match_reset_method, epsilon=epsilon)
             x_i += 1
         return self
 
diff --git a/artlib/hierarchical/SMART.py b/artlib/hierarchical/SMART.py
index 29777d4..a820618 100644
--- a/artlib/hierarchical/SMART.py
+++ b/artlib/hierarchical/SMART.py
@@ -6,7 +6,7 @@
 """
 
 import numpy as np
-from typing import Union, Type, Optional, Iterable
+from typing import Union, Type, Optional, Iterable, Literal
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 from artlib.hierarchical.DeepARTMAP import DeepARTMAP
@@ -32,7 +32,33 @@ def __init__(self, base_ART_class: Type, rho_values: Union[list[float], np.ndarr
             assert isinstance(layer, BaseART), "Only elementary ART-like objects are supported"
         super().__init__(layers)
 
-    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, max_iter=1):
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            prepared data
+        """
+        X_, _ = super(SMART, self).prepare_data([X]*self.n_modules)
+        return X_[0]
+
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        X_, _ = super(SMART, self).restore_data([X] * self.n_modules)
+        return X_[0]
+
+    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, max_iter=1, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
         """
         Fit the model to the data
 
@@ -40,14 +66,20 @@ def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, max_iter=1):
         - X: data set A
         - y: not used
         - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         X_list = [X]*self.n_modules
-        return super().fit(X_list, max_iter=max_iter)
+        return super().fit(X_list, max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
 
-    def partial_fit(self, X: np.ndarray, y: Optional[np.ndarray] = None):
+    def partial_fit(self, X: np.ndarray, y: Optional[np.ndarray] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0):
         X_list = [X] * self.n_modules
-        return self.partial_fit(X_list)
+        return self.partial_fit(X_list, match_reset_method=match_reset_method, epsilon=epsilon)
 
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
diff --git a/artlib/reinforcement/FALCON.py b/artlib/reinforcement/FALCON.py
new file mode 100644
index 0000000..389db87
--- /dev/null
+++ b/artlib/reinforcement/FALCON.py
@@ -0,0 +1,141 @@
+import numpy as np
+from typing import Optional, Literal
+from artlib import FusionART, BaseART, compliment_code, de_compliment_code
+
+
+class FALCON:
+    def __init__(
+            self,
+            state_art: BaseART,
+            action_art: BaseART,
+            reward_art: BaseART,
+            gamma_values: list[float] = np.array([0.33, 0.33, 0.34]),
+            channel_dims = list[int]
+    ):
+        self.fusion_art = FusionART(
+            modules=[state_art, action_art, reward_art],
+            gamma_values=gamma_values,
+            channel_dims=channel_dims
+        )
+
+    def fit(self, states: np.ndarray, actions: np.ndarray, rewards: np.ndarray):
+        data = self.fusion_art.join_channel_data([states, actions, rewards])
+        self.fusion_art = self.fusion_art.fit(data)
+        return self
+
+    def partial_fit(self, states: np.ndarray, actions: np.ndarray, rewards: np.ndarray):
+        data = self.fusion_art.join_channel_data([states, actions, rewards])
+        self.fusion_art = self.fusion_art.partial_fit(data)
+        return self
+
+    def get_actions_and_rewards(self, state: np.ndarray, action_space: Optional[np.ndarray] = None) -> np.ndarray:
+        reward_centers = self.fusion_art.get_channel_centers(2)
+        if action_space is None:
+            action_space = self.fusion_art.get_channel_centers(1)
+            action_space = np.array(action_space)
+        action_space_prepared = self.fusion_art.modules[0].prepare_data(action_space)
+        viable_clusters = []
+        for action in action_space_prepared:
+            data = self.fusion_art.join_channel_data([state.reshape(1, -1), action.reshape(1, -1)], skip_channels=[2])
+            c = self.fusion_art.predict(data, skip_channels=[2])
+            viable_clusters.append(c[0])
+
+        rewards = [reward_centers[c] for c in viable_clusters]
+
+        return action_space, np.array(rewards)
+    def get_action(self, state: np.ndarray, action_space: Optional[np.ndarray] = None, optimality: Literal["min", "max"] = "max") -> np.ndarray:
+        action_space, rewards = self.get_actions_and_rewards(state, action_space)
+        if optimality == "max":
+            c_winner = np.argmax(rewards)
+        else:
+            c_winner = np.argmin(rewards)
+        return action_space[c_winner]
+
+    def get_probabilistic_action(self, state: np.ndarray, action_space: Optional[np.ndarray] = None, offset: float = 0.1, optimality: Literal["min", "max"] = "max") -> np.ndarray:
+        action_space, rewards = self.get_actions_and_rewards(state, action_space)
+        action_indices = np.array(range(len(action_space)))
+
+
+        reward_dist = rewards
+        reward_dist /= np.sum(reward_dist)
+        reward_dist = reward_dist.reshape((-1,))
+
+        if optimality == "min":
+            reward_dist = 1.-reward_dist
+
+        reward_dist = np.maximum(np.minimum(reward_dist, offset), 0.0001)
+        reward_dist /= np.sum(reward_dist)
+
+        a_i = np.random.choice(action_indices, size=1, p=reward_dist)
+        return action_space[a_i[0]][0]
+
+    def get_rewards(self, states: np.ndarray, actions: np.ndarray) -> np.ndarray:
+        reward_centers = self.fusion_art.get_channel_centers(2)
+        data = self.fusion_art.join_channel_data([states, actions], skip_channels=[2])
+        C = self.fusion_art.predict(data, skip_channels=[2])
+        return np.array([reward_centers[c] for c in C])
+
+
+
+
+class TD_FALCON(FALCON):
+    # implements SARSA learning
+
+    def __init__(
+            self,
+            state_art: BaseART,
+            action_art: BaseART,
+            reward_art: BaseART,
+            gamma_values: list[float] = np.array([0.33, 0.33, 0.34]),
+            channel_dims = list[int],
+            td_alpha: float = 1.0,
+            td_lambda: float = 1.0,
+    ):
+        self.td_alpha = td_alpha
+        self.td_lambda = td_lambda
+        super(TD_FALCON, self).__init__(state_art, action_art, reward_art, gamma_values, channel_dims)
+
+    def fit(self, states: np.ndarray, actions: np.ndarray, rewards: np.ndarray):
+        raise NotImplementedError("TD-FALCON can only be trained with partial fit")
+
+    def calculate_SARSA(self, states: np.ndarray, actions: np.ndarray, rewards: np.ndarray, single_sample_reward: Optional[float] = None):
+        # calculate SARSA values
+        rewards_dcc = de_compliment_code(rewards)
+        if len(states) > 1:
+
+            if hasattr(self.fusion_art.modules[0], "W"):
+                # if FALCON has been trained get predicted rewards
+                Q = self.get_rewards(states, actions)
+            else:
+                # otherwise set predicted rewards to 0
+                Q = np.zeros_like(rewards_dcc)
+            # SARSA equation
+            sarsa_rewards = Q[:-1] + self.td_alpha * (rewards_dcc[:-1] + self.td_lambda * Q[1:] - Q[:-1])
+            # ensure SARSA values are between 0 and 1
+            sarsa_rewards = np.maximum(np.minimum(sarsa_rewards, 1.0), 0.0)
+            # compliment code rewards
+            sarsa_rewards_fit = compliment_code(sarsa_rewards)
+            # we cant train on the final state because no rewards are generated after it
+            states_fit = states[:-1, :]
+            actions_fit = actions[:-1, :]
+        else:
+            # if we only have a single sample, we cant learn from future samples
+            if single_sample_reward is None:
+                sarsa_rewards_fit = rewards
+            else:
+                sarsa_rewards_fit = compliment_code(np.array([[single_sample_reward]]))
+            states_fit = states
+            actions_fit = actions
+
+        return states_fit, actions_fit, sarsa_rewards_fit
+
+    def partial_fit(self, states: np.ndarray, actions: np.ndarray, rewards: np.ndarray, single_sample_reward: Optional[float] = None):
+
+        states_fit, actions_fit, sarsa_rewards_fit = self.calculate_SARSA(states, actions, rewards, single_sample_reward)
+        data = self.fusion_art.join_channel_data([states_fit, actions_fit, sarsa_rewards_fit])
+        self.fusion_art = self.fusion_art.partial_fit(data)
+        return self
+
+
+
+
diff --git a/artlib/supervised/ARTMAP.py b/artlib/supervised/ARTMAP.py
index ca2c743..0e61810 100644
--- a/artlib/supervised/ARTMAP.py
+++ b/artlib/supervised/ARTMAP.py
@@ -4,6 +4,7 @@
 Neural Networks, 4, 565 – 588. doi:10.1016/0893-6080(91)90012-T.
 """
 import numpy as np
+from typing import Literal, Tuple
 from artlib.common.BaseART import BaseART
 from artlib.supervised.SimpleARTMAP import SimpleARTMAP
 from sklearn.utils.validation import check_is_fitted
@@ -69,7 +70,31 @@ def validate_data(self, X: np.ndarray, y: np.ndarray):
         self.module_a.validate_data(X)
         self.module_b.validate_data(y)
 
-    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
+    def prepare_data(self, X: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            normalized data
+        """
+        return self.module_a.prepare_data(X), self.module_b.prepare_data(y)
+
+    def restore_data(self, X: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return self.module_a.restore_data(X), self.module_b.restore_data(y)
+
+    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10):
         """
         Fit the model to the data
 
@@ -77,32 +102,44 @@ def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
         - X: data set A
         - y: data set B
         - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         # Check that X and y have correct shape
         self.validate_data(X, y)
 
-        self.module_b.fit(y, max_iter=max_iter)
+        self.module_b.fit(y, max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
 
         y_c = self.module_b.labels_
 
-        super(ARTMAP, self).fit(X, y_c, max_iter=max_iter)
+        super(ARTMAP, self).fit(X, y_c, max_iter=max_iter, match_reset_method=match_reset_method, epsilon=epsilon)
 
         return self
 
 
-    def partial_fit(self, X: np.ndarray, y: np.ndarray):
+    def partial_fit(self, X: np.ndarray, y: np.ndarray, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10):
         """
         Partial fit the model to the data
 
         Parameters:
         - X: data set A
         - y: data set B
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         self.validate_data(X, y)
-        self.module_b.partial_fit(y)
-        super(ARTMAP, self).partial_fit(X, self.labels_b)
+        self.module_b.partial_fit(y, match_reset_method=match_reset_method, epsilon=epsilon)
+        super(ARTMAP, self).partial_fit(X, self.labels_b, match_reset_method=match_reset_method, epsilon=epsilon)
         return self
 
 
@@ -133,3 +170,20 @@ def predict_ab(self, X: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
         """
         check_is_fitted(self)
         return super(ARTMAP, self).predict_ab(X)
+
+    def predict_regression(self, X: np.ndarray) -> np.ndarray:
+        """
+        predict values for the data
+        ARTMAP is not recommended for regression. Use FusionART instead.
+
+        Parameters:
+        - X: data set A
+
+        Returns:
+            predicted values using cluster centers
+
+        """
+        check_is_fitted(self)
+        C = self.predict(X)
+        centers = self.module_b.get_cluster_centers()
+        return np.array([centers[c] for c in C])
diff --git a/artlib/supervised/SimpleARTMAP.py b/artlib/supervised/SimpleARTMAP.py
index 70d89a0..b0bf4d9 100644
--- a/artlib/supervised/SimpleARTMAP.py
+++ b/artlib/supervised/SimpleARTMAP.py
@@ -4,7 +4,7 @@
 Norwell, MA, USA: Kluwer Academic Publishers.
 """
 import numpy as np
-from typing import Optional, Iterable
+from typing import Optional, Iterable, Literal
 from matplotlib.axes import Axes
 from artlib.common.BaseART import BaseART
 from artlib.common.BaseARTMAP import BaseARTMAP
@@ -73,17 +73,47 @@ def validate_data(self, X: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.nd
         - y: data set B
 
         """
-        X, y = check_X_y(X, y)
+        X, y = check_X_y(X, y, dtype=None)
         self.module_a.validate_data(X)
         return X, y
 
-    def step_fit(self, x: np.ndarray, c_b: int) -> int:
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            prepared data
+        """
+        return self.module_a.prepare_data(X)
+
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return self.module_a.restore_data(X)
+
+    def step_fit(self, x: np.ndarray, c_b: int, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10) -> int:
         """
         Fit the model to a single sample
 
         Parameters:
         - x: data sample for side A
         - c_b: side b label
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         Returns:
             side A cluster label
@@ -92,14 +122,14 @@ def step_fit(self, x: np.ndarray, c_b: int) -> int:
         match_reset_func = lambda i, w, cluster, params, cache: self.match_reset_func(
             i, w, cluster, params=params, extra={"cluster_b": c_b}, cache=cache
         )
-        c_a = self.module_a.step_fit(x, match_reset_func=match_reset_func)
+        c_a = self.module_a.step_fit(x, match_reset_func=match_reset_func, match_reset_method=match_reset_method, epsilon=epsilon)
         if c_a not in self.map:
             self.map[c_a] = c_b
         else:
             assert self.map[c_a] == c_b
         return c_a
 
-    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
+    def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10, verbose: bool = False):
         """
         Fit the model to the data
 
@@ -107,6 +137,12 @@ def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
         - X: data set A
         - y: data set B
         - max_iter: number of iterations to fit the model on the same data set
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         # Check that X and y have correct shape
@@ -119,20 +155,31 @@ def fit(self, X: np.ndarray, y: np.ndarray, max_iter=1):
         self.module_a.labels_ = np.zeros((X.shape[0],), dtype=int)
 
         for _ in range(max_iter):
-            for i, (x, c_b) in enumerate(zip(X, y)):
-                c_a = self.step_fit(x, c_b)
+            if verbose:
+                from tqdm import tqdm
+                x_y_iter = tqdm(enumerate(zip(X, y)), total=int(X.shape[0]))
+            else:
+                x_y_iter = enumerate(zip(X, y))
+            for i, (x, c_b) in x_y_iter:
                 self.module_a.pre_step_fit(X)
+                c_a = self.step_fit(x, c_b, match_reset_method=match_reset_method, epsilon=epsilon)
                 self.module_a.labels_[i] = c_a
                 self.module_a.post_step_fit(X)
         return self
 
-    def partial_fit(self, X: np.ndarray, y: np.ndarray):
+    def partial_fit(self, X: np.ndarray, y: np.ndarray, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 1e-10):
         """
         Partial fit the model to the data
 
         Parameters:
         - X: data set A
         - y: data set B
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
 
         """
         SimpleARTMAP.validate_data(self, X, y)
@@ -148,7 +195,7 @@ def partial_fit(self, X: np.ndarray, y: np.ndarray):
             self.module_a.labels_ = np.pad(self.module_a.labels_, [(0, X.shape[0])], mode='constant')
         for i, (x, c_b) in enumerate(zip(X, y)):
             self.module_a.pre_step_fit(X)
-            c_a = self.step_fit(x, c_b)
+            c_a = self.step_fit(x, c_b, match_reset_method=match_reset_method, epsilon=epsilon)
             self.module_a.labels_[i+j] = c_a
             self.module_a.post_step_fit(X)
         return self
diff --git a/artlib/topological/TopoART.py b/artlib/topological/TopoART.py
index c89042e..c54994a 100644
--- a/artlib/topological/TopoART.py
+++ b/artlib/topological/TopoART.py
@@ -9,10 +9,12 @@
 """
 
 import numpy as np
-from typing import Optional, Callable, Iterable
+from typing import Optional, Callable, Iterable, List, Literal
 from matplotlib.axes import Axes
 from warnings import warn
+from copy import deepcopy
 from artlib.common.BaseART import BaseART
+import operator
 
 
 class TopoART(BaseART):
@@ -68,6 +70,30 @@ def validate_data(self, X: np.ndarray):
         """
         self.base_module.validate_data(X)
 
+    def prepare_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        prepare data for clustering
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            normalized data
+        """
+        return self.base_module.prepare_data(X)
+
+    def restore_data(self, X: np.ndarray) -> np.ndarray:
+        """
+        restore data to state prior to preparation
+
+        Parameters:
+        - X: data set
+
+        Returns:
+            restored data
+        """
+        return self.base_module.restore_data(X)
+
     def category_choice(self, i: np.ndarray, w: np.ndarray, params: dict) -> tuple[float, Optional[dict]]:
         """
         get the activation of the cluster
@@ -99,7 +125,7 @@ def match_criterion(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Opt
         """
         return self.base_module.match_criterion(i, w, params, cache)
 
-    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> tuple[bool, dict]:
+    def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None, op: Callable = operator.ge) -> tuple[bool, dict]:
         """
         get the binary match criterion of the cluster
 
@@ -113,7 +139,7 @@ def match_criterion_bin(self, i: np.ndarray, w: np.ndarray, params: dict, cache:
             cluster match criterion binary, cache used for later processing
 
         """
-        return self.base_module.match_criterion_bin(i, w, params, cache)
+        return self.base_module.match_criterion_bin(i, w, params, cache, op)
 
     def update(self, i: np.ndarray, w: np.ndarray, params: dict, cache: Optional[dict] = None) -> np.ndarray:
         """
@@ -202,7 +228,32 @@ def post_step_fit(self, X: np.ndarray):
         if self.sample_counter_ > 0 and self.sample_counter_ % self.tau == 0:
             self.prune(X)
 
-    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -> int:
+    def _match_tracking(self, cache: dict, epsilon: float, params: dict, method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"]) -> bool:
+        M = cache["match_criterion"]
+        if method == "MT+":
+            self.base_module.params["rho"] = M+epsilon
+            return True
+        elif method == "MT-":
+            self.base_module.params["rho"] = M - epsilon
+            return True
+        elif method == "MT0":
+            self.base_module.params["rho"] = M
+            return True
+        elif method == "MT1":
+            self.base_module.params["rho"] = np.inf
+            return False
+        elif method == "MT~":
+            return True
+        else:
+            raise ValueError(f"Invalid Match Tracking Method: {method}")
+
+    def _set_params(self, new_params):
+        self.base_module.params = new_params
+
+    def _deep_copy_params(self) -> dict:
+        return deepcopy(self.base_module.params)
+
+    def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None, match_reset_method: Literal["MT+", "MT-", "MT0", "MT1", "MT~"] = "MT+", epsilon: float = 0.0) -> int:
         """
         fit the model to a single sample
 
@@ -211,11 +262,20 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
         - match_reset_func: a callable accepting the data sample, a cluster weight, the params dict, and the cache dict
             Permits external factors to influence cluster creation.
             Returns True if the cluster is valid for the sample, False otherwise
+        - match_reset_method:
+            "MT+": Original method, rho=M+epsilon
+             "MT-": rho=M-epsilon
+             "MT0": rho=M, using > operator
+             "MT1": rho=1.0,  Immediately create a new cluster on mismatch
+             "MT~": do not change rho
+
 
         Returns:
             cluster label of the input sample
 
         """
+        base_params = self._deep_copy_params()
+        mt_operator = self._match_tracking_operator(match_reset_method)
         self.sample_counter_ += 1
         resonant_c: int = -1
 
@@ -226,22 +286,22 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
             self._permanent_mask = np.zeros((1, ), dtype=bool)
             return 0
         else:
-            T_values, T_cache = zip(*[self.category_choice(x, w, params=self.params) for w in self.W])
+            T_values, T_cache = zip(*[self.category_choice(x, w, params=self.base_module.params) for w in self.W])
             T = np.array(T_values)
-            while any(T > 0):
-                c_ = int(np.argmax(T))
+            while any(~np.isnan(T)):
+                c_ = int(np.nanargmax(T))
                 w = self.W[c_]
                 cache = T_cache[c_]
-                m, cache = self.match_criterion_bin(x, w, params=self.params, cache=cache)
+                m, cache = self.match_criterion_bin(x, w, params=self.base_module.params, cache=cache, op=mt_operator)
                 no_match_reset = (
                         match_reset_func is None or
-                        match_reset_func(x, w, c_, params=self.params, cache=cache)
+                        match_reset_func(x, w, c_, params=self.base_module.params, cache=cache)
                 )
                 if m and no_match_reset:
                     if resonant_c < 0:
-                        params = self.params
+                        params = self.base_module.params
                     else:
-                        params = dict(self.params, **{"beta": self.params["beta_lower"]})
+                        params = dict(self.base_module.params, **{"beta": self.params["beta_lower"]})
                     #TODO: make compatible with DualVigilanceART
                     new_w = self.update(
                         x,
@@ -252,12 +312,18 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
                     self.set_weight(c_, new_w)
                     if resonant_c < 0:
                         resonant_c = c_
-                        T[c_] = -1
+                        T[c_] = np.nan
                     else:
+                        self._set_params(base_params)
                         return resonant_c
                 else:
-                    T[c_] = -1
+                    T[c_] = np.nan
+                    if not no_match_reset:
+                        keep_searching = self._match_tracking(cache, epsilon, self.params, match_reset_method)
+                        if not keep_searching:
+                            T[:] = np.nan
 
+            self._set_params(base_params)
             if resonant_c < 0:
                 c_new = len(self.W)
                 w_new = self.new_weight(x, self.params)
@@ -266,6 +332,15 @@ def step_fit(self, x: np.ndarray, match_reset_func: Optional[Callable] = None) -
 
             return resonant_c
 
+
+    def get_cluster_centers(self) -> List[np.ndarray]:
+        """
+        function for getting centers of each cluster. Used for regression
+        Returns:
+            cluster centroid
+        """
+        return self.base_module.get_cluster_centers()
+
     def plot_cluster_bounds(self, ax: Axes, colors: Iterable, linewidth: int = 1):
         """
         undefined function for visualizing the bounds of each cluster


Test patch:
diff --git a/examples/test_VAT.py b/examples/test_VAT.py
new file mode 100644
index 0000000..6caa12a
--- /dev/null
+++ b/examples/test_VAT.py
@@ -0,0 +1,27 @@
+from sklearn.datasets import make_blobs
+import matplotlib.pyplot as plt
+
+from artlib import VAT
+import time
+
+
+def visualize_blobs():
+    data, target = make_blobs(n_samples=2000, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
+    print("Data has shape:", data.shape)
+
+    t0 = time.time()
+    R, P = VAT(data)
+    dt = time.time() - t0
+    print("time:", dt)
+
+
+    plt.figure()
+    plt.scatter(data[:,0], data[:,1], c=P, cmap="jet", s=10)
+
+    plt.figure()
+    plt.imshow(R)
+    plt.show()
+
+if __name__ == "__main__":
+
+    visualize_blobs()
diff --git a/examples/test_art1.py b/examples/test_art1.py
index 11a77d4..ff5c333 100644
--- a/examples/test_art1.py
+++ b/examples/test_art1.py
@@ -1,7 +1,7 @@
 from sklearn.datasets import make_blobs
 import matplotlib.pyplot as plt
 
-from artlib import ART1, normalize
+from artlib import ART1
 
 
 def cluster_blobs():
@@ -9,8 +9,6 @@ def cluster_blobs():
     data =  ( data > 0.5).astype(int)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
 
     params = {
         "rho": 0.7,
@@ -18,6 +16,10 @@ def cluster_blobs():
         "L": 1.0
     }
     cls = ART1(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_art2.py b/examples/test_art2.py
index e84fa58..9149a2b 100644
--- a/examples/test_art2.py
+++ b/examples/test_art2.py
@@ -1,6 +1,6 @@
 from sklearn.datasets import make_blobs
 import matplotlib.pyplot as plt
-from artlib import ART2A, normalize
+from artlib import ART2A
 
 """
 ==================================================================
@@ -15,15 +15,16 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=2, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
         "rho": 0.2,
         "alpha": 0.0,
         "beta": 1.0,
     }
     cls = ART2A(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_artmap.py b/examples/test_artmap.py
index 542037b..f40b581 100644
--- a/examples/test_artmap.py
+++ b/examples/test_artmap.py
@@ -2,7 +2,7 @@
 from sklearn.datasets import load_iris, make_blobs
 import matplotlib.pyplot as plt
 
-from artlib import FuzzyART, SimpleARTMAP
+from artlib import FuzzyART, SimpleARTMAP, QuadraticNeuronART
 
 
 def cluster_iris():
@@ -12,8 +12,6 @@ def cluster_iris():
     data, target = load_iris(return_X_y=True)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
 
     X_train, X_test, y_train, y_test = train_test_split(X, target, test_size=0.33, random_state=0)
 
@@ -23,6 +21,17 @@ def cluster_iris():
         "beta": 1.0
     }
     art = FuzzyART(**params)
+    # params = {
+    #     "rho": 0.0,
+    #     "s_init": 0.9,
+    #     "lr_b": 0.1,
+    #     "lr_w": 0.1,
+    #     "lr_s": 0.1
+    # }
+    # art = QuadraticNeuronART(**params)
+
+    X = art.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
 
     cls = SimpleARTMAP(art)
 
@@ -44,9 +53,6 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
         "rho": 0.9,
         "alpha": 0.0,
@@ -54,6 +60,18 @@ def cluster_blobs():
     }
     art = FuzzyART(**params)
 
+    # params = {
+    #     "rho": 0.0,
+    #     "s_init": 0.9,
+    #     "lr_b": 0.1,
+    #     "lr_w": 0.1,
+    #     "lr_s": 0.1
+    # }
+    # art = QuadraticNeuronART(**params)
+
+    X = art.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     cls = SimpleARTMAP(art)
 
     cls = cls.fit(X, target)
diff --git a/examples/test_bartmap.py b/examples/test_bartmap.py
index 2661fe6..8da0e97 100644
--- a/examples/test_bartmap.py
+++ b/examples/test_bartmap.py
@@ -1,7 +1,7 @@
 from sklearn.datasets import make_checkerboard
 import matplotlib.pyplot as plt
 
-from artlib import BARTMAP, FuzzyART, normalize
+from artlib import BARTMAP, FuzzyART
 
 
 def cluster_checkerboard():
@@ -11,9 +11,6 @@ def cluster_checkerboard():
     )
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params_a = {
         "rho": 0.6,
         "alpha": 0.0,
@@ -27,6 +24,9 @@ def cluster_checkerboard():
     art_a = FuzzyART(**params_a)
     art_b = FuzzyART(**params_b)
     cls = BARTMAP(art_a, art_b, eta=-1.)
+
+    X = data
+
     cls.fit(X)
 
     print(f"{cls.n_row_clusters} row clusters found")
diff --git a/examples/test_bayesian_art.py b/examples/test_bayesian_art.py
index fe2e6a9..86df104 100644
--- a/examples/test_bayesian_art.py
+++ b/examples/test_bayesian_art.py
@@ -11,7 +11,6 @@
 sys.path.append(directory.parent.parent)
 
 from artlib import BayesianART
-from artlib.common.utils import normalize
 import numpy as np
 
 
@@ -19,14 +18,16 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 2e-5,
+        "rho": 7e-5,
         "cov_init": np.array([[1e-4, 0.0], [0.0, 1e-4]]),
     }
     cls = BayesianART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_convex_hull_art.py b/examples/test_convex_hull_art.py
new file mode 100644
index 0000000..88c2073
--- /dev/null
+++ b/examples/test_convex_hull_art.py
@@ -0,0 +1,32 @@
+from sklearn.datasets import make_blobs
+import matplotlib.pyplot as plt
+import numpy as np
+
+from artlib.experimental.ConvexHullART import ConvexHullART, plot_convex_polygon
+from scipy.spatial import ConvexHull
+
+
+def cluster_blobs():
+    data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
+    print("Data has shape:", data.shape)
+
+
+    params = {
+        "rho": 0.85,
+        "merge_rho": 0.8
+    }
+    cls = ConvexHullART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
+    y = cls.fit_predict(X)
+
+    print(f"{cls.n_clusters} clusters found")
+
+    cls.visualize(X, y)
+    plt.show()
+
+
+if __name__ == "__main__":
+    cluster_blobs()
diff --git a/examples/test_cvi_art.py b/examples/test_cvi_art.py
new file mode 100644
index 0000000..d308443
--- /dev/null
+++ b/examples/test_cvi_art.py
@@ -0,0 +1,31 @@
+from sklearn.datasets import make_blobs
+import matplotlib.pyplot as plt
+
+from artlib import CVIART, FuzzyART
+
+
+def cluster_blobs():
+    data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
+    print("Data has shape:", data.shape)
+
+    params = {
+        "rho": 0.7,
+        "alpha": 0.0,
+        "beta": 1.0,
+    }
+    art = FuzzyART(**params)
+    cls = CVIART(art, validity=CVIART.CALINSKIHARABASZ)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
+    y = cls.fit_predict(X)
+
+    print(f"{cls.n_clusters} clusters found")
+
+    cls.visualize(X, y)
+    plt.show()
+
+
+if __name__ == "__main__":
+    cluster_blobs()
\ No newline at end of file
diff --git a/examples/test_dual_vigilance_art.py b/examples/test_dual_vigilance_art.py
index f6e323e..578e121 100644
--- a/examples/test_dual_vigilance_art.py
+++ b/examples/test_dual_vigilance_art.py
@@ -8,9 +8,6 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
         "rho": 0.85,
         "alpha": 0.8,
@@ -18,6 +15,10 @@ def cluster_blobs():
     }
     base_art = FuzzyART(**params)
     cls = DualVigilanceART(base_art, 0.78)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_ellipsoid_art.py b/examples/test_ellipsoid_art.py
index 6821da4..bf1b2b0 100644
--- a/examples/test_ellipsoid_art.py
+++ b/examples/test_ellipsoid_art.py
@@ -1,24 +1,25 @@
 from sklearn.datasets import make_blobs
 import matplotlib.pyplot as plt
 
-from artlib import EllipsoidART, normalize
+from artlib import EllipsoidART
 
 
 def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 0.2,
+        "rho": 0.01,
         "alpha": 0.0,
         "beta": 1.0,
-        "r_hat": 0.6,
+        "r_hat": 0.65,
         "mu": 0.8
     }
     cls = EllipsoidART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_fusion_art.py b/examples/test_fusion_art.py
index 3492035..01c6747 100644
--- a/examples/test_fusion_art.py
+++ b/examples/test_fusion_art.py
@@ -12,12 +12,6 @@ def cluster_blobs():
     data_channel_a = data[:,0].reshape((-1,1))
     data_channel_b = data[:,1].reshape((-1,1))
 
-    X_channel_a = FuzzyART.prepare_data(data_channel_a)
-    X_channel_b = FuzzyART.prepare_data(data_channel_b)
-
-    X = np.hstack([X_channel_a, X_channel_b])
-    print("Prepared data has shape:", X.shape)
-
     params = {
         "rho": 0.5,
         "alpha": 0.0,
@@ -25,7 +19,14 @@ def cluster_blobs():
     }
     art_a = FuzzyART(**params)
     art_b = FuzzyART(**params)
-    cls = FusionART([art_a, art_b], gamma_values=[0.5, 0.5], channel_dims=[2,2])
+    cls = FusionART([art_a, art_b], gamma_values=np.array([0.5, 0.5]), channel_dims=[2,2])
+
+    X_channel_a = art_a.prepare_data(data_channel_a)
+    X_channel_b = art_b.prepare_data(data_channel_b)
+
+    X = cls.join_channel_data(channel_data=[X_channel_a, X_channel_b])
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_fuzzy_art.py b/examples/test_fuzzy_art.py
index b5729be..29be523 100644
--- a/examples/test_fuzzy_art.py
+++ b/examples/test_fuzzy_art.py
@@ -8,15 +8,16 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 0.7,
+        "rho": 0.5,
         "alpha": 0.0,
         "beta": 1.0
     }
     cls = FuzzyART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_gaussian_art.py b/examples/test_gaussian_art.py
index 9bba697..f384a2f 100644
--- a/examples/test_gaussian_art.py
+++ b/examples/test_gaussian_art.py
@@ -2,21 +2,22 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from artlib import GaussianART, normalize
+from artlib import GaussianART
 
 
 def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 0.15,
+        "rho": 0.05,
         "sigma_init": np.array([0.5, 0.5]),
     }
     cls = GaussianART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_grid_search_cv.py b/examples/test_grid_search_cv.py
index 370b1d2..a63fb58 100644
--- a/examples/test_grid_search_cv.py
+++ b/examples/test_grid_search_cv.py
@@ -9,11 +9,6 @@ def grid_search_blobs():
     data, target = load_iris(return_X_y=True)
     print("Data has shape:", data.shape)
 
-    X = HypersphereART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
-    X_train, X_test, y_train, y_test = train_test_split(X, target, test_size=0.33, random_state=0)
-
     base_params = {
         "rho": 0.3,
         "alpha": 0.0,
@@ -21,6 +16,12 @@ def grid_search_blobs():
         "r_hat": 0.8
     }
     art = HypersphereART(**base_params)
+
+    X = art.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
+    X_train, X_test, y_train, y_test = train_test_split(X, target, test_size=0.33, random_state=0)
+
     cls = SimpleARTMAP(art)
     # print(cls.get_params(deep=False))
     # raise ValueError
diff --git a/examples/test_grid_search_cv_elementary.py b/examples/test_grid_search_cv_elementary.py
index 57230a5..c6b1e13 100644
--- a/examples/test_grid_search_cv_elementary.py
+++ b/examples/test_grid_search_cv_elementary.py
@@ -18,9 +18,6 @@ def grid_search_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     base_params = {
         "rho": 0.7,
         "alpha": 0.0,
@@ -28,6 +25,9 @@ def grid_search_blobs():
     }
     cls = FuzzyART(**base_params)
 
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     param_grid = {
         "rho": [0.0, 0.25, 0.5, 0.7, 0.9],
         "alpha": [0.0],
diff --git a/examples/test_hypersphere_art.py b/examples/test_hypersphere_art.py
index e1d76c1..99c7a36 100644
--- a/examples/test_hypersphere_art.py
+++ b/examples/test_hypersphere_art.py
@@ -1,23 +1,24 @@
 from sklearn.datasets import make_blobs
 import matplotlib.pyplot as plt
 
-from artlib import HypersphereART, normalize
+from artlib import HypersphereART
 
 
 def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 0.7,
+        "rho": 0.5,
         "alpha": 0.0,
         "beta": 1.0,
         "r_hat": 0.8
     }
     cls = HypersphereART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_icvi_art.py b/examples/test_icvi_art.py
new file mode 100644
index 0000000..3919708
--- /dev/null
+++ b/examples/test_icvi_art.py
@@ -0,0 +1,31 @@
+from sklearn.datasets import make_blobs
+import matplotlib.pyplot as plt
+
+from artlib import iCVIFuzzyART
+
+
+def cluster_blobs():
+    data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
+    print("Data has shape:", data.shape)
+
+    params = {
+        "rho": 0.0,
+        "alpha": 0.0,
+        "beta": 1.0,
+        "validity": iCVIFuzzyART.CALINSKIHARABASZ
+    }
+    cls = iCVIFuzzyART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
+    y = cls.fit_predict(X)
+
+    print(f"{cls.n_clusters} clusters found")
+
+    cls.visualize(X, y)
+    plt.show()
+
+
+if __name__ == "__main__":
+    cluster_blobs()
\ No newline at end of file
diff --git a/examples/test_quadratic_neuron_art.py b/examples/test_quadratic_neuron_art.py
index 099f655..c90cefb 100644
--- a/examples/test_quadratic_neuron_art.py
+++ b/examples/test_quadratic_neuron_art.py
@@ -1,16 +1,13 @@
 from sklearn.datasets import make_blobs
 import matplotlib.pyplot as plt
 
-from artlib import QuadraticNeuronART, normalize
+from artlib import QuadraticNeuronART
 
 
 def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = normalize(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
         "rho": 0.9,
         "s_init": 0.9,
@@ -19,6 +16,10 @@ def cluster_blobs():
         "lr_s": 0.1
     }
     cls = QuadraticNeuronART(**params)
+
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     y = cls.fit_predict(X)
 
     print(f"{cls.n_clusters} clusters found")
diff --git a/examples/test_regression.py b/examples/test_regression.py
new file mode 100644
index 0000000..d4c1a64
--- /dev/null
+++ b/examples/test_regression.py
@@ -0,0 +1,86 @@
+import numpy as np
+from sklearn.datasets import load_iris, make_blobs
+import matplotlib.pyplot as plt
+
+from artlib import FuzzyART, ARTMAP, FusionART, normalize
+
+def make_data():
+    x_ = np.arange(0.0,4.0, 0.01)
+    x = x_+0.5
+    y = np.sin(x)+np.cos(x**2)+np.log(x**3)+np.cos(x)*np.log(x**2)
+
+    x_norm, _, _ = normalize(x)
+    y_norm, _, _ = normalize(y)
+    return x_norm.reshape((-1,1)), y_norm.reshape((-1,1))
+
+
+def fit_regression_ARTMAP():
+    X_, y_ = make_data()
+    print("Data has shape:", X_.shape)
+
+    params = {
+        "rho": 0.95,
+        "alpha": 0.0,
+        "beta": 1.0
+    }
+    art_a = FuzzyART(**params)
+    art_b = FuzzyART(**params)
+
+    X = art_a.prepare_data(X_)
+    y = art_b.prepare_data(y_)
+    print("Prepared data has shape:", X.shape)
+
+    cls = ARTMAP(art_a, art_b)
+
+    cls = cls.fit(X, y)
+
+    print(f"{len(np.unique(cls.labels_a))} clusters found")
+
+    y_pred = cls.predict_regression(X)
+
+    plt.plot(X_, y_, 'r-', label="original")
+    plt.plot(X_, y_pred, "b-", label="ARTMAP")
+    plt.legend()
+    plt.show()
+
+def fit_regression_FusionART():
+    X_, y_ = make_data()
+    print("Data has shape:", X_.shape)
+
+    params = {
+        "rho": 0.95,
+        "alpha": 0.0,
+        "beta": 1.0
+    }
+    art_a = FuzzyART(**params)
+    art_b = FuzzyART(**params)
+
+    X = art_a.prepare_data(X_)
+    y = art_b.prepare_data(y_)
+    print("Prepared data has shape:", X.shape)
+
+    cls = FusionART([art_a, art_b], gamma_values=np.array([0.5, 0.5]), channel_dims=[2, 2])
+
+    Xy = cls.join_channel_data(channel_data=[X, y])
+
+    cls = cls.fit(Xy)
+
+    print(f"{cls.n_clusters} clusters found")
+
+    X_no_y = cls.join_channel_data(channel_data=[X], skip_channels=[-1])
+
+    y_pred = cls.predict_regression(X_no_y, target_channels=[-1])
+
+    plt.plot(X_, y_, 'r-', label="original")
+    plt.plot(X_, y_pred, "b-", label="FusionART")
+    plt.legend()
+    plt.show()
+
+def main():
+    # fit_regression_ARTMAP()
+    fit_regression_FusionART()
+
+
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/examples/test_reinforcement_learning.py b/examples/test_reinforcement_learning.py
new file mode 100644
index 0000000..faa0172
--- /dev/null
+++ b/examples/test_reinforcement_learning.py
@@ -0,0 +1,270 @@
+import numpy as np
+from tqdm import tqdm
+from artlib import TD_FALCON, FuzzyART, compliment_code
+import gymnasium as gym
+from collections import defaultdict
+import matplotlib.pyplot as plt
+from copy import deepcopy
+
+# only works with Fuzzy ART based FALCON
+def prune_clusters(cls):
+    # get existing state and action weights
+    state = np.array(cls.fusion_art.modules[0].W)
+    action = np.array(cls.fusion_art.modules[1].W)
+    # get predicted rewards for each cluster
+    reward = np.array(cls.fusion_art.get_channel_centers(2)).reshape((-1, 1))
+
+    # combine state and actions
+    combined = np.round(np.hstack((state, action)), decimals=3)
+    # identify unique combinations
+    unique_combined, indices = np.unique(combined, axis=0, return_inverse=True)
+    # get mean reward prediction for each unique state-action pair
+    unique_rewards = np.array([reward[indices == i, :].mean() for i in range(len(unique_combined))])
+
+    # split unique state-action pairs
+    unique_states = unique_combined[:, :state.shape[1]]
+    unique_actions = unique_combined[:, state.shape[1]:]
+
+    # update model to only have unique state-action pairs and their average rewards
+    cls.fusion_art.modules[0].W = [row for row in unique_states]
+    cls.fusion_art.modules[1].W = [row for row in unique_actions]
+    cls.fusion_art.modules[2].W = [row for row in compliment_code(unique_rewards.reshape((-1, 1)))]
+
+    return cls
+
+def update_FALCON(records, cls, shrink_ratio):
+    # convert records into arrays
+    states = np.array(records["states"]).reshape((-1,1))
+    actions = np.array(records["actions"]).reshape((-1,1))
+    rewards = np.array(records["rewards"]).reshape((-1,1))
+
+    # compliment code states and actions
+    states_cc = compliment_code(states)
+    actions_cc = compliment_code(actions)
+    rewards_cc = compliment_code(rewards)
+
+    # states_fit, actions_fit, sarsa_rewards_fit = cls.calculate_SARSA(states_cc, actions_cc, rewards_cc, single_sample_reward=1.0)
+
+    # if FALCON has been previously trained
+    if hasattr(cls.fusion_art.modules[0], "W"):
+        # remove any duplicate clusters
+        cls = prune_clusters(cls)
+        # # shrink clusters to account for dynamic programming changes
+        # for m in range(3):
+        #     cls.fusion_art.modules[m] = cls.fusion_art.modules[m].shrink_clusters(shrink_ratio)
+
+    # fit FALCON to data
+    # data = cls.fusion_art.join_channel_data([states_fit, actions_fit, sarsa_rewards_fit])
+    # cls.fusion_art = cls.fusion_art.partial_fit(data)
+    cls = cls.partial_fit(states_cc, actions_cc, rewards_cc, single_sample_reward=1.0)
+    for m in range(3):
+        cls.fusion_art.modules[m] = cls.fusion_art.modules[m].shrink_clusters(shrink_ratio)
+
+    return cls
+
+def training_cycle(cls, epochs, steps, sarsa_alpha, sarsa_gamma, render_mode=None, shrink_ratio=0.1, explore_rate=0.0, checkpoint_frequency=50, early_stopping_reward=-20):
+    # create the environment
+    env = gym.make('CliffWalking-v0', render_mode=render_mode)
+    # define some constants
+    ACTION_SPACE = np.array([[0], [1.], [2.], [3.]])
+    STATE_MAX = 47
+    ACTION_MAX = 3
+    REWARD_MAX = 150
+    cls.td_alpha = sarsa_alpha
+    cls.td_lambda = sarsa_gamma
+
+    best_reward = -np.inf
+    best_cls = None
+
+    # track reward history
+    reward_history = []
+
+    pbar = tqdm(range(epochs))
+    for e in pbar:
+        observation, info = env.reset()
+        n_observation = observation / STATE_MAX
+        records = {"states": [], "actions": [], "rewards": []}
+        past_states = []
+        for _ in range(steps):
+            # get an action
+            observation_cc = compliment_code(np.array([n_observation]).reshape(1, -1))
+            if np.random.random() < explore_rate:
+                action = int(np.random.choice(ACTION_SPACE.flatten()))
+            else:
+                cls_action = cls.get_action(observation_cc, action_space=ACTION_SPACE, optimality="min")
+                action = int(cls_action[0])
+            # normalize state and action
+            n_observation = observation / STATE_MAX
+            n_action = action / ACTION_MAX
+
+            # record state and action for training
+            records["states"].append(n_observation)
+            records["actions"].append(n_action)
+            past_states.append(observation)
+
+            # take a step
+            observation, reward, terminated, truncated, info = env.step(action)
+
+            # check reward value
+            if reward > -100:
+                # punish circular paths
+                if observation in past_states:
+                    reward = -2
+            # normalize and record reward from step
+            n_reward = abs(reward) / REWARD_MAX
+            records["rewards"].append(n_reward)
+
+            # check if epoch is done
+            if terminated or truncated or reward == -100:
+                break
+
+        # train FALCON
+        cls = update_FALCON(records, cls, shrink_ratio)
+        # record sum of rewards generated during this epoch
+        reward_history.append(-sum(records["rewards"])*REWARD_MAX)
+
+        # if this isnt random exploration
+        if explore_rate < 1.0:
+            # see if we should save a checkpoint
+            if (e+1)%checkpoint_frequency == 0 or e == epochs-1:
+                # test model
+                cls, test_reward_history = demo_cycle(cls, 1, steps, render_mode=None)
+                # check if our current model is better than the best previous model
+                if test_reward_history[0] >= best_reward or best_cls is None:
+                    # same a checkpoint
+                    best_cls = deepcopy(cls)
+                    best_reward = test_reward_history[0]
+                    # check early stopping condition
+                    if best_reward > early_stopping_reward:
+                        # show current best reward on progress bar
+                        pbar.set_postfix({'Best Reward': best_reward})
+                        return cls, reward_history
+                else:
+                    # restore previous best model
+                    cls = deepcopy(best_cls)
+            # show current best reward on progress bar
+            pbar.set_postfix({'Best Reward': best_reward})
+
+    env.close()
+    return cls, reward_history
+
+
+def demo_cycle(cls, epochs, steps, render_mode=None):
+    # create the environment
+    env = gym.make('CliffWalking-v0', render_mode=render_mode)
+    # define some constants
+    ACTION_SPACE = np.array([[0], [1.], [2.], [3.]])
+    STATE_MAX = 47
+    ACTION_MAX = 3
+    REWARD_MAX = 150
+
+    # track reward history
+    reward_history = []
+
+    for _ in range(epochs):
+        observation, info = env.reset()
+        n_observation = observation / STATE_MAX
+        records = {"states": [], "actions": [], "rewards": []}
+        past_states = []
+        for _ in range(steps):
+
+            # get an action
+            observation_cc = compliment_code(np.array([n_observation]).reshape(1, -1))
+            cls_action = cls.get_action(observation_cc, action_space=ACTION_SPACE, optimality="min")
+            action = int(cls_action[0])
+
+            # normalize state and action
+            n_observation = observation / STATE_MAX
+            n_action = action / ACTION_MAX
+
+            # record state and action
+            records["states"].append(n_observation)
+            records["actions"].append(n_action)
+            past_states.append(observation)
+
+            # take a step
+            observation, reward, terminated, truncated, info = env.step(action)
+
+            # check reward value
+            if reward > -100:
+                # punish circular paths
+                if observation in past_states:
+                    reward = -2
+            # normalize and record reward from step
+            n_reward = abs(reward) / REWARD_MAX
+            records["rewards"].append(n_reward)
+
+            # check if epoch is done
+            if terminated or truncated or reward == -100:
+                break
+        # record sum of rewards generated during this epoch
+        reward_history.append(-sum(records["rewards"])*REWARD_MAX)
+
+    env.close()
+    return cls, reward_history
+
+
+def max_up_to_each_element(lst):
+    max_list = []
+    current_max = float('-inf')  # Start with the smallest possible value
+    for num in lst:
+        current_max = max(current_max, num)
+        max_list.append(current_max)
+    return max_list
+
+
+def train_FALCON():
+    # define training regimen
+    training_regimen = [
+        {"name": "random", "epochs": 1000, "shrink_ratio": 0.3, "gamma": 0.0, "explore_rate": 1.0, "render_mode": None},
+        {"name": "explore 33%", "epochs": 500, "shrink_ratio": 0.3, "gamma": 0.2, "explore_rate": 0.333, "render_mode": None},
+        {"name": "explore 20%", "epochs": 500, "shrink_ratio": 0.3, "gamma": 0.2, "explore_rate": 0.20, "render_mode": None},
+        {"name": "explore 5%", "epochs": 1000, "shrink_ratio": 0.3, "gamma": 0.2, "explore_rate": 0.05, "render_mode": None},
+    ]
+    MAX_STEPS = 25
+    SARSA_ALPHA = 1.0
+
+    # define parameters for state, action, and reward modules
+    art_state = FuzzyART(rho=0.99,alpha=0.01, beta=1.0)
+    art_action = FuzzyART(rho=0.99,alpha=0.01, beta=1.0)
+    art_reward = FuzzyART(rho=0.0,alpha=0.01, beta=1.0)
+    # instantiate FALCON
+    cls = TD_FALCON(art_state, art_action, art_reward, channel_dims=[2, 2, 2])
+    # record rewards for each epoch
+    all_rewards = []
+    # initialize FALCON with random exploration
+    print("Starting Training")
+    for regimen in training_regimen:
+        print(f"Starting Training Cycle: {regimen['name']}")
+        cls, reward_history = training_cycle(
+            cls,
+            epochs=regimen["epochs"],
+            steps=MAX_STEPS,
+            sarsa_alpha=SARSA_ALPHA,
+            sarsa_gamma=regimen["gamma"],
+            render_mode=regimen["render_mode"],
+            shrink_ratio=regimen["shrink_ratio"],
+            explore_rate=regimen["explore_rate"]
+        )
+        all_rewards.extend(reward_history)
+        print("MAX REWARD:",max(reward_history))
+
+    # demo learned policy
+    cls, reward_history = demo_cycle(cls, epochs=2, steps=25, render_mode="human")
+    print(reward_history)
+    all_rewards.extend(reward_history)
+    max_rewards = max_up_to_each_element(all_rewards)
+
+    # plot reward history
+    plt.figure()
+    plt.plot(list(range(len(all_rewards))), all_rewards, "r-")
+    plt.plot(list(range(len(all_rewards))), max_rewards, "b-")
+    plt.xlabel("Epoch")
+    plt.ylabel("Reward")
+    plt.title("Rewards over Time")
+    plt.show()
+
+if __name__ == "__main__":
+    # This takes approximately 3 minutes
+    np.random.seed(42)
+    train_FALCON()
\ No newline at end of file
diff --git a/examples/test_seq_art.py b/examples/test_seq_art.py
new file mode 100644
index 0000000..1828c77
--- /dev/null
+++ b/examples/test_seq_art.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+from artlib.experimental.SeqART import SeqART
+
+
+sequences = [
+    "hello",
+    "helsinki",
+    "hella",
+    "hell",
+    "arkansas",
+    "kansas",
+    "advent",
+    "adventure",
+    "ventilation",
+    "addition",
+    "subtraction",
+    "substitution",
+]
+
+def cluster_sequences():
+    X = np.array(sequences, dtype=object)
+
+    params = {
+        "rho": -0.5
+    }
+    cls = SeqART(**params)
+
+    y = cls.fit_predict(X)
+
+    print(f"{cls.n_clusters} clusters found")
+
+    for i in range(cls.n_clusters):
+        print(f"Cluster {i}")
+        for s_i in range(len(sequences)):
+            if i == y[s_i]:
+                print(sequences[s_i])
+        print(cls.W[i])
+        print("=" * 20)
+
+
+if __name__ == "__main__":
+    cluster_sequences()
+
+
diff --git a/examples/test_smart.py b/examples/test_smart.py
index 090f147..11f76cc 100644
--- a/examples/test_smart.py
+++ b/examples/test_smart.py
@@ -8,16 +8,16 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     base_params = {
         "alpha": 0.0,
         "beta": 1.0
     }
-    rho_values = [0.7, 0.85, 0.9]
+    rho_values = [0.5, 0.85, 0.9]
     cls = SMART(FuzzyART, rho_values, base_params)
 
+    X = cls.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     cls = cls.fit(X)
     y = cls.layers[0].labels_
 
diff --git a/examples/test_topo_art.py b/examples/test_topo_art.py
index d346d0e..3e7bda7 100644
--- a/examples/test_topo_art.py
+++ b/examples/test_topo_art.py
@@ -8,15 +8,16 @@ def cluster_blobs():
     data, target = make_blobs(n_samples=150, centers=3, cluster_std=0.50, random_state=0, shuffle=False)
     print("Data has shape:", data.shape)
 
-    X = FuzzyART.prepare_data(data)
-    print("Prepared data has shape:", X.shape)
-
     params = {
-        "rho": 0.6,
+        "rho": 0.4,
         "alpha": 0.8,
         "beta": 1.0
     }
     base_art = FuzzyART(**params)
+
+    X = base_art.prepare_data(data)
+    print("Prepared data has shape:", X.shape)
+
     cls = TopoART(base_art, betta_lower=0.3, tau=150, phi=35)
     cls = cls.fit(X, max_iter=5)
     y = cls.labels_