diff --git a/bob/bio/vein/algorithm/Correlate.py b/bob/bio/vein/algorithm/Correlate.py
index e35310617c82d465e96a244b07494f229f1125e1..416e26c2fa74fb9f7e9bd569f27383fa9fd80e8f 100644
--- a/bob/bio/vein/algorithm/Correlate.py
+++ b/bob/bio/vein/algorithm/Correlate.py
@@ -4,10 +4,10 @@
 import numpy
 import skimage.feature
 
-from bob.bio.base.algorithm import Algorithm
+from bob.bio.base.pipelines import BioAlgorithm
 
 
-class Correlate(Algorithm):
+class Correlate(BioAlgorithm):
     """Correlate probe and model without cropping
 
     The method is based on "cross-correlation" between a model and a probe image.
@@ -18,18 +18,34 @@ class Correlate(Algorithm):
 
     """
 
-    def __init__(self):
-
-        # call base class constructor
-        Algorithm.__init__(
-            self, multiple_model_scoring=None, multiple_probe_scoring=None
+    def __init__(
+        self,
+        probes_score_fusion="max",
+        enrolls_score_fusion="mean",
+        **kwargs,
+    ):
+        super().__init__(
+            probes_score_fusion=probes_score_fusion,
+            enrolls_score_fusion=enrolls_score_fusion,
+            **kwargs,
         )
 
-    def enroll(self, enroll_features):
-        """Enrolls the model by computing an average graph for each model"""
+    def create_templates(self, feature_sets, enroll):
+        return feature_sets
 
-        # return the generated model
-        return numpy.array(enroll_features)
+    def compare(self, enroll_templates, probe_templates):
+        # returns scores NxM where N is the number of enroll templates and M is the number of probe templates
+        # enroll_templates is Nx?1xD
+        # probe_templates is Mx?2xD
+        scores = []
+        for enroll in enroll_templates:
+            scores.append([])
+            for probe in probe_templates:
+                s = [[self.score(e, p) for p in probe] for e in enroll]
+                s = self.fuse_probe_scores(s, axis=1)
+                s = self.fuse_enroll_scores(s, axis=0)
+                scores[-1].append(s)
+        return numpy.array(scores)
 
     def score(self, model, probe):
         """Computes the score between the probe and the model.
@@ -49,21 +65,13 @@ class Correlate(Algorithm):
 
         image_ = probe.astype(numpy.float64)
 
-        if len(model.shape) == 2:
-            model = numpy.array([model])
-
-        scores = []
-
-        # iterate over all models for a given individual
-        for md in model:
-
-            R = md.astype(numpy.float64)
-            Nm = skimage.feature.match_template(image_, R)
+        R = model.astype(numpy.float64)
+        Nm = skimage.feature.match_template(image_, R)
 
-            # figures out where the maximum is on the resulting matrix
-            t0, s0 = numpy.unravel_index(Nm.argmax(), Nm.shape)
+        # figures out where the maximum is on the resulting matrix
+        t0, s0 = numpy.unravel_index(Nm.argmax(), Nm.shape)
 
-            # this is our output
-            scores.append(Nm[t0, s0])
+        # this is our output
+        score = Nm[t0, s0]
 
-        return numpy.mean(scores)
+        return score
diff --git a/bob/bio/vein/algorithm/HammingDistance.py b/bob/bio/vein/algorithm/HammingDistance.py
index ff38223a6d3eac13da316af2c69790c6b770e1ae..0ab27eee6cec98f4891ebd08da0c7a7ce73c4e0a 100644
--- a/bob/bio/vein/algorithm/HammingDistance.py
+++ b/bob/bio/vein/algorithm/HammingDistance.py
@@ -12,24 +12,22 @@ class HammingDistance(Distance):
     base :py:class:`bob.bio.base.algorithm.Distance`.
 
     The input to this function should be of binary nature (boolean arrays). Each
-    binary input is first flattened to form a one-dimensional vector. The `Hamming
-    distance <https://en.wikipedia.org/wiki/Hamming_distance>`_ is then
+    binary input is first flattened to form a one-dimensional vector. The
+    `Hamming distance <https://en.wikipedia.org/wiki/Hamming_distance>`_ is then
     calculated between these two binary vectors.
 
     The current implementation uses :py:func:`scipy.spatial.distance.hamming`,
     which returns a scalar 64-bit ``float`` to represent the proportion of
     mismatching corresponding bits between the two binary vectors.
 
-    The base class constructor parameter ``is_distance_function`` is set to
-    ``False`` on purpose to ensure that calculated distances are returned as
-    positive values rather than negative.
-
+    The base class constructor parameter ``factor`` is set to ``1`` on purpose
+    to ensure that calculated distances are returned as positive values rather
+    than negative.
     """
 
-    def __init__(self):
-        from scipy.spatial.distance import hamming
-
+    def __init__(self, **kwargs):
         super(HammingDistance, self).__init__(
-            distance_function=hamming,
-            is_distance_function=False,
+            distance_function="hamming",
+            factor=1,
+            **kwargs,
         )
diff --git a/bob/bio/vein/algorithm/MiuraMatch.py b/bob/bio/vein/algorithm/MiuraMatch.py
index 774c8f239127af0286e26e803b5ff5cc388f4482..4e9ebace9bfb17612fc45bb2e64ae94662e65865 100644
--- a/bob/bio/vein/algorithm/MiuraMatch.py
+++ b/bob/bio/vein/algorithm/MiuraMatch.py
@@ -47,25 +47,35 @@ class MiuraMatch(BioAlgorithm):
         self,
         ch=80,  # Maximum search displacement in y-direction
         cw=90,  # Maximum search displacement in x-direction
+        probes_score_fusion="max",
+        enrolls_score_fusion="mean",
+        **kwargs,
     ):
-
-        # call base class constructor
-        BioAlgorithm.__init__(
-            self,
-            ch=ch,
-            cw=cw,
-            # multiple_model_scoring = None,
-            # multiple_probe_scoring = None
+        super().__init__(
+            probes_score_fusion=probes_score_fusion,
+            enrolls_score_fusion=enrolls_score_fusion,
+            **kwargs,
         )
 
         self.ch = ch
         self.cw = cw
 
-    def enroll(self, enroll_features):
-        """Enrolls the model by computing an average graph for each model"""
+    def create_templates(self, feature_sets, enroll):
+        return feature_sets
 
-        # return the generated model
-        return numpy.array(enroll_features)
+    def compare(self, enroll_templates, probe_templates):
+        # returns scores NxM where N is the number of enroll templates and M is the number of probe templates
+        # enroll_templates is Nx?1xD
+        # probe_templates is Mx?2xD
+        scores = []
+        for enroll in enroll_templates:
+            scores.append([])
+            for probe in probe_templates:
+                s = [[self.score(e, p) for p in probe] for e in enroll]
+                s = self.fuse_probe_scores(s, axis=1)
+                s = self.fuse_enroll_scores(s, axis=0)
+                scores[-1].append(s)
+        return numpy.array(scores)
 
     def score(self, model, probe):
         """Computes the score between the probe and the model.
@@ -85,57 +95,36 @@ class MiuraMatch(BioAlgorithm):
 
         image_ = probe.astype(numpy.float64)
 
-        if len(model.shape) == 2:
-            model = numpy.array([model])
-
-        scores = []
+        md = model
+        # erode model by (ch, cw)
+        R = md.astype(numpy.float64)
+        h, w = R.shape  # same as I
+        crop_R = R[self.ch : h - self.ch, self.cw : w - self.cw]
+
+        # correlates using scipy - fastest option available iff the self.ch and
+        # self.cw are height (>30). In this case, the number of components
+        # returned by the convolution is high and using an FFT-based method
+        # yields best results. Otherwise, you may try  the other options bellow
+        # -> check our test_correlation() method on the test units for more
+        # details and benchmarks.
+        Nm = scipy.signal.fftconvolve(image_, numpy.rot90(crop_R, k=2), "valid")
+        # 2nd best: use convolve2d or correlate2d directly;
+        # Nm = scipy.signal.convolve2d(I, numpy.rot90(crop_R, k=2), 'valid')
+        # 3rd best: use correlate2d
+        # Nm = scipy.signal.correlate2d(I, crop_R, 'valid')
+
+        # figures out where the maximum is on the resulting matrix
+        t0, s0 = numpy.unravel_index(Nm.argmax(), Nm.shape)
+
+        # this is our output
+        Nmm = Nm[t0, s0]
+
+        # normalizes the output by the number of pixels lit on the input
+        # matrices, taking into consideration the surface that produced the
+        # result (i.e., the eroded model and part of the probe)
+        score = Nmm / (
+            crop_R.sum()
+            + image_[t0 : t0 + h - 2 * self.ch, s0 : s0 + w - 2 * self.cw].sum()
+        )
 
-        # iterate over all models for a given individual
-        for md in model:
-            # erode model by (ch, cw)
-            R = md.astype(numpy.float64)
-            h, w = R.shape  # same as I
-            crop_R = R[self.ch : h - self.ch, self.cw : w - self.cw]
-
-            # correlates using scipy - fastest option available iff the self.ch and
-            # self.cw are height (>30). In this case, the number of components
-            # returned by the convolution is high and using an FFT-based method
-            # yields best results. Otherwise, you may try  the other options bellow
-            # -> check our test_correlation() method on the test units for more
-            # details and benchmarks.
-            Nm = scipy.signal.fftconvolve(
-                image_, numpy.rot90(crop_R, k=2), "valid"
-            )
-            # 2nd best: use convolve2d or correlate2d directly;
-            # Nm = scipy.signal.convolve2d(I, numpy.rot90(crop_R, k=2), 'valid')
-            # 3rd best: use correlate2d
-            # Nm = scipy.signal.correlate2d(I, crop_R, 'valid')
-
-            # figures out where the maximum is on the resulting matrix
-            t0, s0 = numpy.unravel_index(Nm.argmax(), Nm.shape)
-
-            # this is our output
-            Nmm = Nm[t0, s0]
-
-            # normalizes the output by the number of pixels lit on the input
-            # matrices, taking into consideration the surface that produced the
-            # result (i.e., the eroded model and part of the probe)
-            scores.append(
-                Nmm
-                / (
-                    crop_R.sum()
-                    + image_[
-                        t0 : t0 + h - 2 * self.ch, s0 : s0 + w - 2 * self.cw
-                    ].sum()
-                )
-            )
-
-        return [numpy.mean(scores)]
-
-    def score_multiple_biometric_references(self, biometric_references, data):
-        if isinstance(biometric_references, list):
-            return [self.score(model, data) for model in biometric_references]
-        else:
-            raise ValueError(
-                "The model does not have the desired format (list, array, ...)"
-            )
+        return score
diff --git a/bob/bio/vein/database/roi_annotation.py b/bob/bio/vein/database/roi_annotation.py
index aa7b550a9f2bc052380c6d8db2bd8474bd0b5b29..f888aa9a96d8428da381d61e5861f1c1c05a2c0c 100644
--- a/bob/bio/vein/database/roi_annotation.py
+++ b/bob/bio/vein/database/roi_annotation.py
@@ -20,7 +20,6 @@ class ROIAnnotation(TransformerMixin, BaseEstimator):
 
     def _more_tags(self):
         return {
-            "stateless": True,
             "requires_fit": False,
         }
 
diff --git a/bob/bio/vein/database/utfvp.py b/bob/bio/vein/database/utfvp.py
index bc751026d7471197ef16fd7b4c92f1d2d85d37cd..a125724e30003b3d73fe151372695fdfb6919a7c 100644
--- a/bob/bio/vein/database/utfvp.py
+++ b/bob/bio/vein/database/utfvp.py
@@ -179,7 +179,7 @@ class UtfvpDatabase(CSVDataset):
                 ),
                 ROIAnnotation(roi_path=rc.get("bob.bio.vein.utfvp.roi", "")),
             ),
-            allow_scoring_with_all_biometric_references=True,
+            score_all_vs_all=True,
         )
 
     @staticmethod
diff --git a/bob/bio/vein/database/verafinger_contactless.py b/bob/bio/vein/database/verafinger_contactless.py
index 6730a5eac528b81c88a283085d23beab51109d6c..46e65851318bc3794a8b9fcd1806ca63230062d8 100644
--- a/bob/bio/vein/database/verafinger_contactless.py
+++ b/bob/bio/vein/database/verafinger_contactless.py
@@ -71,7 +71,7 @@ class VerafingerContactless(CSVDataset):
                 extension="",
                 reference_id_equal_subject_id=False,
             ),
-            allow_scoring_with_all_biometric_references=True,
+            score_all_vs_all=True,
         )
 
     @staticmethod
diff --git a/bob/bio/vein/tests/test_tools.py b/bob/bio/vein/tests/test_tools.py
index 0f0570cc2cc3aada6ded9da2ebd911344ff38951..8a7a257cfa1ee1fef01bffae03166bdfeae0db4c 100644
--- a/bob/bio/vein/tests/test_tools.py
+++ b/bob/bio/vein/tests/test_tools.py
@@ -674,16 +674,16 @@ def test_hamming_distance():
     # Tests on simple binary arrays:
     # 1.) Maximum HD (1.0):
     model_1 = numpy.array([0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0])
-    probe_1 = numpy.array([[1, 0, 0, 0, 1, 0], [0, 1, 1, 1, 0, 1]])
-    score_max = HD.score(model_1, probe_1)
+    probe_1 = numpy.array([[1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1]])
+    score_max = HD.compare([model_1], [probe_1])[0, 0]
     assert score_max == 1.0
     # 2.) Minimum HD (0.0):
     model_2 = numpy.array([0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1])
-    probe_2 = numpy.array([[0, 1, 1, 1, 0, 1], [0, 1, 1, 1, 0, 1]])
-    score_min = HD.score(model_2, probe_2)
+    probe_2 = numpy.array([[0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1]])
+    score_min = HD.compare([model_2], [probe_2])[0, 0]
     assert score_min == 0.0
     # 3.) HD of exactly half (0.5)
     model_3 = numpy.array([0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0])
-    probe_3 = numpy.array([[0, 1, 1, 1, 0, 1], [0, 1, 1, 1, 0, 1]])
-    score_half = HD.score(model_3, probe_3)
+    probe_3 = numpy.array([[0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1]])
+    score_half = HD.compare([model_3], [probe_3])[0, 0]
     assert score_half == 0.5