Repository: evotext/seqsim (language: python)

Issue:
Merge action
Attempt at triggering GitHub action


Discussion / hints:

Fix patch:
diff --git a/AUTHORS.md b/AUTHORS.md
index 33f8543..11b4ff4 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -4,7 +4,9 @@ Credits
 Development
 -----------
 
-* Tiago Tresoldi <tresoldi@gmail.com>
+* Tiago Tresoldi <tiago.tresoldi@lingfil.uu.se>
+* Luke Maurits
+* Michael Dunn
 
 Contributors
 ------------
diff --git a/README.md b/README.md
index 422ec1e..b95f8f7 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 # seqsim
 
 [![PyPI](https://img.shields.io/pypi/v/seqsim.svg)](https://pypi.org/project/seqsim)
-[![CI](https://github.com/tresoldi/seqsim/actions/workflows/main.yml/badge.svg)](https://github.com/tresoldi/seqsim/actions/workflows/main.yml)
+[![CI](https://github.com/evotext/seqsim/actions/workflows/main.yml/badge.svg)](https://github.com/evotext/seqsim/actions/workflows/main.yml)
 
 Python library for computing measures of similarity for sequences of hashable data types
 
@@ -28,11 +28,11 @@ The library offers different methods to compare sequences of hashable Python ele
 >> test3_seq_b = (1, 2, 4, 3, 6, 7)
 >> test4_seq_a = (1, 2, 3)
 >> test4_seq_b = ["a", "b", "c", "d"]
->> seqsim.distance.prev_edit_distance(test1_seq_a, test1_seq_b)
+>> seqsim.edit.prev_edit_distance(test1_seq_a, test1_seq_b)
 3.0
->> seqsim.distance.jaccard(test3_seq_a, test3_seq_b)
+>> seqsim.edit.jaccard_dist(test3_seq_a, test3_seq_b)
 0.4285714285714286
->> seqsim.distance.mmcwpa(test3_seq_a, test3_seq_b)
+>> seqsim.edit.mmcwpa_dist(test3_seq_a, test3_seq_b)
 0.5546382285848768
 ```
 
@@ -58,22 +58,24 @@ Contributing guidelines, including a code of conduct, can be found in the
 
 ## Author and citation
 
-The library is developed by Tiago Tresoldi (tiago.tresoldi@lingfil.uu.se).
+The library is developed in the context of "[Cultural Evolution of Text](https://www.evotext.se)",
+project, with funding from the Riksbankens Jubileumsfond (grant agreement ID:
+[MXM19-1087:1](https://www.rj.se/en/anslag/2019/cultural-evolution-of-texts/)).
 
 If you use `seqsim`, please cite it as:
 
-> Tresoldi, Tiago (2021). seqsim, a library for computing measures of similarity for
-> sequences of hashable data types. Version 0.2. Uppsala: Uppsala Universitet.
-> Available at: https://github.com/tresoldi/seqsim
+> Tresoldi, Tiago; Maurits, Luke; Dunn, Michael. (2021). seqsim, a library for computing
+> measures of similarity for sequences of hashable data types. Version 0.3.
+> Uppsala: Uppsala universitet. Available at: https://github.com/evotext/seqsim
 
 In BibTeX:
 
 ```
-@misc{Tresoldi2021titivillus,
-  author = {Tresoldi, Tiago},
-  title = {seqsim, a library for computing measures of similarity for sequences of hashable data types},
-  howpublished = {\url{https://github.com/tresoldi/seqsim}},
+@misc{Tresoldi2021seqsim,
+  author = {Tresoldi, Tiago; Maurits, Luke; Dunn, Michael},
+  title = {seqsim, a library for computing measures of similarity for sequences of hashable data types. Version 0.3},
+  howpublished = {\url{https://github.com/evotext/seqsim}},
   address = {Uppsala},
-  publisher = {Uppsala Universitet},
+  publisher = {Uppsala universitet},
   year = {2021},
 }
diff --git a/requirements.txt b/requirements.txt
index 296d654..20a9350 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,1 +1,2 @@
-numpy
\ No newline at end of file
+numpy
+textdistance
\ No newline at end of file
diff --git a/src/seqsim/__init__.py b/src/seqsim/__init__.py
index f17c18b..8c8b93b 100644
--- a/src/seqsim/__init__.py
+++ b/src/seqsim/__init__.py
@@ -1,15 +1,25 @@
 """
 seqsim __init__.py
+
+We follow the mathematical definitions for distinguishing between "similarity"
+and "distance", as the latter must have the following properties:
+
+  * positivity: d(x,y) >= 0
+  * symmetry: d(x,y) = d(y,x)
+  * identity-discerning: d(x,y) = 0 => x = y
+  * triangle inequality: d(x,z) <= d(x,y) + d(y,z)
 """
 
-# Version of the seqsim package
-__author__ = "Tiago Tresoldi"
+# Version of the `seqsim` package
+__author__ = "Tiago Tresoldi, Luke Maurits, Michael Dunn"
 __email__ = "tiago.tresoldi@lingfil.uu.se"
 __version__ = "0.3"
 
 # Import local modules
-from . import distance
-from . import similarity
+from . import edit
+from . import token
+from . import sequence
+from . import compression
 from .common import set_seeds, collect_subseqs
 from .ngrams import ngrams_iter, get_all_ngrams_by_order
 
diff --git a/src/seqsim/common.py b/src/seqsim/common.py
index 0c71198..c114d64 100644
--- a/src/seqsim/common.py
+++ b/src/seqsim/common.py
@@ -5,7 +5,7 @@
 # Import Python standard libraries
 import hashlib
 import random
-from typing import Union, List, Optional, Sequence
+from typing import Callable, Hashable, Union, List, Optional, Sequence
 
 # Import 3rd party libraries
 import numpy as np
@@ -13,7 +13,7 @@
 
 def set_seeds(seed: Union[str, float, int]) -> None:
     """
-    Set seeds globally from the user provided one.
+    Set seeds globally from the one provided by the user.
 
     The function takes care of reproducibility and allows to use strings and
     floats as seed for `numpy` as well.
@@ -63,6 +63,7 @@ def sequence_find(hay: Sequence, needle: Sequence) -> Optional[int]:
     return None
 
 
+# TODO: replace with the ngram collector module
 def collect_subseqs(sequence: Sequence, sort: bool = True) -> List[Sequence]:
     """
     Collects all possible sub-sequences in a given sequence.
@@ -113,6 +114,76 @@ def collect_subseqs(sequence: Sequence, sort: bool = True) -> List[Sequence]:
         length -= 1
 
     if sort:
-        ret = sorted(ret, key=lambda e: (len(e), e))
+        # We try to sort normally; if there is a TypeError, such as when the list has mixed
+        # ints and strings, we sort by the string representation of all elements
+        # TODO: do it in a better way
+        try:
+            ret = sorted(ret, key=lambda e: (len(e), e))
+        except TypeError:
+            ret = sorted(ret, key=lambda e: (len(str(e)), str(e)))
 
     return ret
+
+
+def _nwise(I, n):
+    """
+    Iterate through I n at a time, e.g.
+        _nwise("abcd", 2) -> ab, bc, cd
+    """
+    for i in range(len(I) + 1 - n):
+        yield I[i : i + n]
+
+
+def _indices(L: Sequence[Hashable], element: Hashable) -> list:
+    """
+    Find all _indices in `L` matching `element`, e.g.
+        _indices("abcab", "a") -> 0, 3
+    """
+
+    return [idx for idx, value in enumerate(L) if value == element]
+
+
+# TODO: continue `Callable` typing
+def _wagner_fischer(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    costs_fn: Callable,
+    d: Optional[List[List[float]]] = None,
+) -> float:
+    """
+    Computes arbitrary edit distances using the Wagner-Fischer algorithm.
+
+    See: https://en.wikipedia.org/wiki/Wagner%E2%80%93Fischer_algorithm
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param costs_fn: A cost function, specific to a particular edit distance.
+    @param d: An optional "starting matrix".
+    @return: The cost distance.
+    """
+
+    # Cache lengths
+    m = len(seq_x)
+    n = len(seq_y)
+
+    # If we haven't been provided a custom initial matrix specific to
+    # our particular distance measure, create a generic on here.  This
+    # fills out the first row and first column of the cost matrix,
+    # corresponding to inserting or deleting all characters.
+    if not d:
+        d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
+        for i in range(1, m + 1):
+            d[i][0] = i
+        for j in range(1, n + 1):
+            d[0][j] = j
+
+    # Flood fill the rest of the matrix with values computed using our
+    # cost function
+    for j in range(1, n + 1):
+        for i in range(1, m + 1):
+            costs = costs_fn(seq_x, seq_y, d, i, j)
+            d[i][j] = min(costs)
+
+    # Return lower-right element of matrix, which is the minimum cost
+    # for transforming s into t
+    return d[m][n]
diff --git a/src/seqsim/compression.py b/src/seqsim/compression.py
new file mode 100644
index 0000000..7da45d2
--- /dev/null
+++ b/src/seqsim/compression.py
@@ -0,0 +1,68 @@
+"""
+Module implementing various methods for similarity and distance from compression methods.
+
+Most of the methods are commonly used in string comparison from normalized
+compression distance, such as from Arithmetic Coding, but in this module we need
+to make sure we can operate on arbitrary iterable data structures.
+"""
+
+# Import Python standard libraries
+from typing import Hashable, Sequence
+import logging
+
+# Import 3rd-party libraries
+import textdistance
+
+# TODO: multiple sequences?
+# TODO: consider splitting in similarity/distance
+def arith_ncd(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a distance between two sequences based on Arithmetic Coding.
+
+    See: https://en.wikipedia.org/wiki/Arithmetic_coding
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Arithmetic Coding NCD between the two sequences.
+    """
+
+    # As the method uses .find, we need strings
+    # TODO: rewrite proper solution generalizing for all sequences
+    seq_x = "".join([str(elem) for elem in seq_x])
+    seq_y = "".join([str(elem) for elem in seq_y])
+
+    # TODO: have a normalization method based on seq_x*2 and seq_y*2
+    if normal:
+        logging.warning("Arithmetic Coding NCD cannot be normalized in range [0..1].")
+
+    return textdistance.arith_ncd(seq_x, seq_y)
+
+
+# TODO: multiple sequences?
+def entropy_ncd(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a distance between two sequences based on entropy.
+
+    See: https://en.wikipedia.org/wiki/Entropy_(information_theory)
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Entropy NCD between the two sequences.
+    """
+
+    # As the method uses .find, we need strings
+    # TODO: rewrite proper solution generalizing for all sequences
+    seq_x = "".join([str(elem) for elem in seq_x])
+    seq_y = "".join([str(elem) for elem in seq_y])
+
+    # TODO: have a normalization method based on seq_x*2 and seq_y*2
+    if normal:
+        logging.warning("Entropy NCD cannot be normalized in range [0..1].")
+
+    return textdistance.entropy_ncd(seq_x, seq_y)
diff --git a/src/seqsim/distance.py b/src/seqsim/distance.py
deleted file mode 100644
index 7c53860..0000000
--- a/src/seqsim/distance.py
+++ /dev/null
@@ -1,498 +0,0 @@
-"""
-Module implementing various methods for distance in sequences.
-
-Most of the methods are commonly used in string comparison, such as Levenshtein
-distance, but in this module we need to make sure we can operate on arbitrary
-iterable data structures.
-"""
-
-# TODO: redo distance/similarity
-# TODO: allow some normalization for all methods?
-
-# Import Python standard libraries
-from typing import Sequence, List, Tuple
-
-# Import 3rd party libraries
-import numpy as np
-
-# Import local modules
-from .common import sequence_find, collect_subseqs
-
-from . import similarity
-
-# TODO: review code, especially (N,N) (shoud probably delegate it to similarity)
-def birnbaum(M, N):
-    if len(N) > len(M):
-        M, N = N, M
-
-    distance = 1 - (similarity.birnbaum(M, N) / similarity.birnbaum(N, N))
-
-    # Distance can be less than zero if M contains a perfect N plus doublets
-    # matching part of N as well.
-    return max(0, distance)
-
-
-def levenshtein(seq_a, seq_b):
-    """
-    Compute the Levenshtein distance between two manuscripts, using the
-    Wagner-Fischer algorithm.
-    See: https://en.wikipedia.org/wiki/Levenshtein_distance
-    """
-    return _wagner_fischer(seq_a, seq_b, _levenshtein_costs)
-
-
-def _levenshtein_costs(seq_a, seq_b, d, i, j):
-    """
-    Compute the candidate costs for entry (i,j) of the Levenshtein distance
-    matrix, each cost corresponding to one of the available edit operations.
-    """
-    substitution_cost = 0 if seq_a[i - 1] == seq_b[j - 1] else 1
-    costs = (
-        d[i - 1][j] + 1,  # del
-        d[i][j - 1] + 1,  # ins
-        d[i - 1][j - 1] + substitution_cost,
-    )
-    return costs
-
-
-def norm_levenshtein(seq_a, seq_b):
-    """
-    Compute a normalised form of the Levenshtein distance between two
-    manuscripts, guaranteeing a result between 0 and 1 regardless of
-    manuscript length.
-    """
-    return levenshtein(seq_a, seq_b) / max([len(seq_a), len(seq_b)])
-
-
-def levdamerau(seq_a, seq_b):
-    """
-    Compute the Damerau-Levenshtein distance between two manuscripts,
-    using the Wagner-Fischer algorithm.
-    See: https://en.wikipedia.org/wiki/Damerau%E2%80%93Levenshtein_distance
-    """
-    return _wagner_fischer(seq_a, seq_b, _levdamerau_costs)
-
-
-def _levdamerau_costs(seq_a, seq_b, d, i, j):
-    """
-    Compute the candidate costs for entry (i,j) of the Damerau-Levenshtein
-    distance matrix, each cost corresponding to one of the available edit
-    operations.
-    """
-    # Start out as per Levenshtein
-    substitution_cost = 0 if seq_a[i - 1] == seq_b[j - 1] else 1
-    costs = [
-        d[i - 1][j] + 1,
-        d[i][j - 1] + 1,
-        d[i - 1][j - 1] + substitution_cost,
-    ]
-    # Add the transposition cost
-    if (
-        i > 1
-        and j > 1
-        and seq_a[i - 1] == seq_b[j - 2]
-        and seq_a[i - 2] == seq_b[j - 1]
-    ):
-        costs.append(d[i - 2][j - 2] + 1)
-    return costs
-
-
-def fragile_ends(seq_a, seq_b):
-    """
-    Compute the "fragile ends" distance between two manuscripts - like
-    Levenshtein but with deletions in the initial or final 10% of the
-    positions being discounted.
-    This distance measure is not used directly in the paper and was a
-    proof-of-concept while working toward the "Stemmatological distance".
-    """
-    d = _fragile_ends_initial_matrix(seq_a, seq_b)
-    return _wagner_fischer(seq_a, seq_b, _fragile_ends_costs, d)
-
-
-def _fragile_ends_initial_matrix(seq_a, seq_b):
-    """
-    Compute the starting matrix for the Wagner-Fischer algorithm when using
-    "fragile ends" costs.  The generic starting matrix isn't applicable here,
-    as the cost of deleting the entire starting manuscript is less than the
-    manuscript length.
-    """
-    m = len(seq_a)
-    n = len(seq_b)
-    d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
-    for i in range(1, m + 1):
-        if i <= round(0.1 * m) or i >= round(0.9 * m):
-            d[i][0] = d[i - 1][0] + 0.5
-        else:
-            d[i][0] = d[i - 1][0] + 1
-    # Insertions
-    for j in range(1, n + 1):
-        d[0][j] = j
-    return d
-
-
-def _fragile_ends_costs(seq_a, seq_b, d, i, j):
-    """
-    Compute the candidate costs for entry (i,j) of the "fragile ends"
-    distance matrix, each cost corresponding to one of the available edit
-    operations.
-    """
-    substitution_cost = 0 if seq_a[i - 1] == seq_b[j - 1] else 1
-    costs = [
-        d[i][j - 1] + 1,
-        d[i - 1][j - 1] + substitution_cost,
-    ]
-    # Discount deletion near ends
-    m = len(seq_a)
-    if i <= round(0.1 * m) or i >= round(0.9 * m):
-        costs.append(d[i - 1][j] + 0.5)
-    else:
-        costs.append(d[i - 1][j] + 1)
-    return costs
-
-
-def bulk_delete(seq_a, seq_b, max_del_len=5):
-    """
-    Compute the "bulk delete" distance between two manuscripts - like
-    Levenshtein but with deletions of consecutive stories up to the specified
-    maximum length counting as a single operation.
-    This distance measure is not used directly in the paper and was a
-    proof-of-concept while working toward the "Stemmatological distance".
-    """
-    d = _bulk_delete_initial_matrix(seq_a, seq_b, max_del_len)
-    _bulk_delete_costs = _bulk_delete_costs_factory(max_del_len)
-    return _wagner_fischer(seq_a, seq_b, _bulk_delete_costs, d)
-
-
-def _bulk_delete_initial_matrix(seq_a, seq_b, max_del_len):
-    """
-    Compute the starting matrix for the Wagner-Fischer algorithm when using the
-    "bulk delete distance" costs.  The generic starting matrix isn't applicable
-    here, as the cost of deleting the entire starting manuscript is less than
-    the manuscript length.
-    """
-    m = len(seq_a)
-    n = len(seq_b)
-    d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
-    # Bulk deletions
-    for i in range(1, m + 1):
-        bulk_deletions = int(i / max_del_len)
-        d[i][0] = bulk_deletions
-        remainder = i - bulk_deletions * max_del_len
-        if remainder:
-            d[i][0] += 1
-    # Insertions
-    for j in range(1, n + 1):
-        d[0][j] = j
-    return d
-
-
-def _bulk_delete_costs_factory(max_del_len=5):
-    """
-    Define and return a function which will compute the candidate costs for
-    entry (i,j) of the "bulk deletion" distance matrix with a particular maximum
-    deletion block length.
-    """
-
-    def _bulk_delete_costs(seq_a, seq_b, d, i, j):
-        substitution_cost = 0 if seq_a[i - 1] == seq_b[j - 1] else 1
-        costs = [
-            d[i][j - 1] + 1,  # ins
-            d[i - 1][j - 1] + substitution_cost,
-        ]
-        for n in range(1, min(max_del_len + 1, i)):
-            # Delete consecutive block of n
-            costs.append(d[i - n][j] + 1)
-        return costs
-
-    return _bulk_delete_costs
-
-
-def stemmatological(seq_a, seq_b, max_del_len=5, frag_start=10, frag_end=10):
-    """
-    Compute the "stemmatological distance" between two manuscripts.  This is
-    essentially a combination of the above "fragile ends" and "bulk delete"
-    methods, with the fragile ends method generalised a little to allow
-    specifying the size of both fragile regions.
-    """
-    d = _stemmatological_initial_matrix(seq_a, seq_b, max_del_len, frag_start, frag_end)
-    _stemmatology_costs = _stemmatological_costs_factory(
-        max_del_len, frag_start, frag_end
-    )
-    return _wagner_fischer(seq_a, seq_b, _stemmatology_costs, d)
-
-
-def _stemmatological_initial_matrix(
-    seq_a, seq_b, max_del_len=5, frag_start=10, frag_end=10
-):
-    """
-    Compute the starting matrix for the Wagner-Fischer algorithm when using the
-    "stemmatological distance" costs.  The generic starting matrix isn't
-    applicable here, as per the fragile ends and bulk delete methods above.
-    """
-    m = len(seq_a)
-    n = len(seq_b)
-    d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
-    lower = round(m * (frag_start) / 100.0)
-    upper = round(m * (100 - frag_end) / 100.0)
-    for i in range(1, m + 1):
-        if i <= lower or i >= upper:
-            d[i][0] = d[i - min(i, max_del_len)][0] + 0.5
-        else:
-            d[i][0] = d[i - min(i, max_del_len)][0] + 1.0
-    # Insertions
-    for j in range(1, n + 1):
-        d[0][j] = j
-    return d
-
-
-def _stemmatological_costs_factory(max_del_len=5, frag_start=10, frag_end=10):
-    """
-    Define and return a function which will compute the candidate costs for
-    entry (i,j) of the "stemmatological distance" matrix with specific parameter
-    values.
-    """
-
-    def _stemmatological_costs(seq_a, seq_b, d, i, j):
-        substitution_cost = 0 if seq_a[i - 1] == seq_b[j - 1] else 1
-        costs = [
-            d[i][j - 1] + 1,
-            d[i - 1][j - 1] + substitution_cost,
-        ]
-        m = len(seq_a)
-        lower = round(m * (frag_start) / 100.0)
-        upper = round(m * (100 - frag_end) / 100.0)
-        # Delete consecutive block of n
-        for n in range(1, min(max_del_len, i)):
-            # Discount bulk deletion near ends
-            if i <= lower or i >= upper:
-                costs.append(d[i - n][j] + 0.5)
-            else:
-                costs.append(d[i - n][j] + 1)
-        return costs
-
-    return _stemmatological_costs
-
-
-def norm_stemmatological(seq_a, seq_b):
-    """
-    Compute a normalised form of the "stemmatological distance"
-    between two manuscripts, guaranteeing a result between 0 and 1 regardless of
-    manuscript length.
-    """
-    return stemmatological(seq_a, seq_b) / max([len(seq_a), len(seq_b)])
-
-
-def norm_stemmatological_2030(seq_a, seq_b):
-    """
-    Compute a normalised form of the "stemmatological distance"
-    between two manuscripts, with the first 20% and final 30% of the manuscript
-    being "fragile".
-    """
-    return stemmatological(seq_a, seq_b, 5, 20, 30) / max([len(seq_a), len(seq_b)])
-
-
-def _wagner_fischer(seq_a, seq_b, costs_fn, d=None):
-    """
-    Implements the Wagner-Fischer algorithm to compute arbitrary edit
-    distances.
-    s is the "source" sequence (of length m)
-    t is the "target" sequence (of length n)
-    costs_fn is a cost function, specific to a particular edit distance
-    d is an optional "starting matrix"
-    See: https://en.wikipedia.org/wiki/Wagner%E2%80%93Fischer_algorithm
-    """
-    m = len(seq_a)
-    n = len(seq_b)
-
-    # If we haven't been provided a custom initial matrix specific to
-    # our particular distance measure, create a generic on here.  This
-    # fills out the first row and first column of the cost matrix,
-    # corresponding to inserting or deleting all characters.
-    if not d:
-        d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
-        for i in range(1, m + 1):
-            d[i][0] = i
-        for j in range(1, n + 1):
-            d[0][j] = j
-
-    # Flood fill the rest of the matrix with values computed using our
-    # cost function
-    for j in range(1, n + 1):
-        for i in range(1, m + 1):
-            costs = costs_fn(seq_a, seq_b, d, i, j)
-            d[i][j] = min(costs)
-
-    # Return lower-right element of matrix, which is the minimum cost
-    # for transforming s into t
-    return d[m][n]
-
-
-def jaccard(seq_a: Sequence, seq_b: Sequence) -> float:
-    """
-    Computes a Jaccard distance between two sequences.
-
-    :param seq_a: The first sequence of elements to be compared.
-    :param seq_b: The second sequence of elements to be compared.
-    :return: The Jaccard distance between the two sequences.
-    """
-
-    intersection = len(set(seq_a).intersection(seq_b))
-    union = len(seq_a) + len(seq_b) - intersection
-    return 1.0 - (float(intersection) / union)
-
-
-# TODO: rename from jaccard, as it is not really jaccard...
-def subseq_jaccard(seq_a: Sequence, seq_b: Sequence) -> float:
-    """
-    Computes a Jaccard distance between two sequences using sub-sequence occurrence.
-
-    :param seq_a: The first sequence of elements to be compared.
-    :param seq_b: The second sequence of elements to be compared.
-    :return: The Subseq-Jaccard distance between the two sequences.
-    """
-
-    subseqs1 = collect_subseqs(seq_a)
-    subseqs2 = collect_subseqs(seq_b)
-
-    # From the longest subseq, which is the length of the longest sequence,
-    # collect all subsequences of that given length in both sets, compute the
-    # Jaccard index, correct it by length of the subsequence (so that longer ones
-    # will score higher) and update the internal results before returning.
-    # TODO: collect beforehand in a list, so we don´t repeat the comprehension
-    # TODO: convert to tuple beforehand as well?
-    jaccard = []
-    max_length = max([len(seq_a), len(seq_b)])
-    for length in range(max_length, 0, -1):
-        l_subseq1 = [tuple(ss) for ss in subseqs1 if len(ss) == length]
-        l_subseq2 = [tuple(ss) for ss in subseqs2 if len(ss) == length]
-
-        intersection = len(set(l_subseq1).intersection(l_subseq2))
-        union = len(l_subseq1) + len(l_subseq2) - intersection
-        jaccard.append((float(intersection) / union) * length)
-
-    # Compute the denominator, as the highest possible value
-    den = (max_length * (max_length + 1)) / 2.0
-
-    return (1.0 - (sum(jaccard) / den)) ** max_length
-
-
-def _mmcwpa(
-    f_a: List[Sequence], f_b: List[Sequence], ssnc: float
-) -> Tuple[List[Sequence], List[Sequence], float]:
-    """
-    Internal function for MMCWPA implementation.
-
-    In this implementation of the Modified Moving Contracting Window Pattern Algorithm
-    (MMCWPA) to calculate sequence similarity, we return a list of non-overlapping,
-    non-contiguous fields Fx, a list of non-overlapping, non-contiguous fields Fy, and
-    the SSNC value (the Sum of the Square of the Number of the same characters). This
-    function separates the core method of the implementation and makes recursive calls
-    easier.
-
-    :param f_a: A list of sub-sequences, related to the first sequence.
-    :param f_b: A list of sub-sequences, related to the second sequence.
-    :param ssnc: The previous SSNC value.
-    :return: A tuple whose first element is a list of remaining sub-sequences from the
-             first sequence, the second element is a list of remaining sub-sequences
-             from the second sequence, and the third element is the updated SSNC.
-    """
-
-    # the boolean value indicating if a total or partial
-    # match was found between subfields Fx and Fy; when
-    # a match is found, the variable is used to cascade
-    # out of the loops of the function
-    match = False
-
-    # the variables where to store the new collections of
-    # subfields, if any match is found; if these values
-    # are not changed and the empty lists are returned,
-    # stringcomp() will break the loop of comparison,
-    # calculate the similarity ratio and return its value
-    new_f_a, new_f_b = [], []
-
-    # search patterns in all subfields of Fx; the index of
-    # the subfield in the list is used for upgrading the
-    # list, if a pattern is a found
-    for idx_a, sf_a in enumerate(f_a):
-        # 'length' stores the length of the sliding window,
-        # from full length to a single character
-        for length in range(len(sf_a), 0, -1):
-            # 'i' stores the starting index of the sliding
-            # window in Fx
-            for i in range(len(sf_a) - length + 1):
-                # extract the pattern for matching
-                pattern = sf_a[i : i + length]
-
-                # look for the pattern in Fy
-                for idx_b, sf_b in enumerate(f_b):
-                    # 'j' stores the starting index in Fy; the
-                    # Python find() function returns -1 if there
-                    # is no match
-                    j = sequence_find(sf_b, pattern)
-                    if j is not None:
-                        # the pattern was found; set 'new_fx' and
-                        # 'new_fy' to version of 'fx' and 'fy' with
-                        # the patterns removed, update the SSNC and
-                        # set 'match' as True, in order to cascade
-                        # out of the loops
-                        tmp_x = [sf_a[:i], sf_a[i + length :]]
-                        tmp_y = [sf_b[:j], sf_b[j + length :]]
-                        new_f_a = f_a[:idx_a] + tmp_x + f_a[idx_a + 1 :]
-                        new_f_b = f_b[:idx_b] + tmp_y + f_b[idx_b + 1 :]
-
-                        ssnc += (2 * length) ** 2
-
-                        match = True
-
-                        break
-
-                    # if the current match was found, end search
-                    if match:
-                        break
-
-                # if a match was found, end the sliding window
-                if match:
-                    break
-
-            # if a match was found, end Fx subfield enumeration
-            if match:
-                break
-
-        # remove any empty subfields due to pattern removal
-        new_f_a = [sf for sf in new_f_a if sf]
-        new_f_b = [sf for sf in new_f_b if sf]
-
-        return new_f_a, new_f_b, ssnc
-
-
-def mmcwpa(seq_a: Sequence, seq_b: Sequence) -> float:
-    """
-    Computes an MMCWPA distance between two sequences.
-
-    MMCWPA is the Modifier Moving Contracting Window Pattern Algorithm, modified by
-    Tiago Tresoldi from a method published by Yang et al. (2001). In order to simplify
-    the logic, the function uses the auxiliary internal function `_mmcwpa()`.
-
-    Reference for the original method:
-    Q. X. Yang, S. S. Yuan, L. Chun, L. Zhao, S. Peng. "Faster Algorithm of String
-    Comparison", eprint arXiv:cs/0112022, December 2001.
-
-    :param seq_a: The first sequence of elements to be compared.
-    :param seq_b: The second sequence of elements to be compared.
-    :return: The MMCWPA distance between the two sequences.
-
-    """
-
-    # Cache original lengths, before any modifications
-    # TODO: make sure it works in a function way, without touching `seq_a`
-    #       and `seq_b`, so we don't need to cache this
-    len_a, len_b = len(seq_a), len(seq_b)
-
-    # Initialize the `f_a` and `f_b` vectors with a single element each,
-    # the corresponding sequence; we also initialize the `ssnc` to zero.
-    f_a, f_b = [seq_a], [seq_b]
-    ssnc: float = 0.0
-    while f_a and f_b:
-        f_a, f_b, ssnc = _mmcwpa(f_a, f_b, ssnc)
-
-    return 1.0 - ((ssnc / ((len_a + len_b) ** 2.0)) ** 0.5)
diff --git a/src/seqsim/edit.py b/src/seqsim/edit.py
new file mode 100644
index 0000000..fe84b22
--- /dev/null
+++ b/src/seqsim/edit.py
@@ -0,0 +1,955 @@
+"""
+Module implementing various methods for similarity and distance from edit methods.
+
+Most of the methods are commonly used in string comparison, such as Levenshtein
+distance, but in this module we need to make sure we can operate on arbitrary
+iterable data structures.
+"""
+
+# Import Python standard libraries
+from typing import Callable, Hashable, List, Sequence, Tuple
+import difflib
+import logging
+
+# Import 3rd-party libraries
+import textdistance
+
+# Import local modules
+from .common import sequence_find, _nwise, _indices, _wagner_fischer
+
+# Methods for sequence similarity
+# -------------------------------
+#
+# Methods for "similarity" are those that do not guarantee all of the three
+# distance properties of positivity, symmetry, identity-discerning, and
+# triangle inequality
+
+
+def birnbaum_simil(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Birnbaum similarity score with the original method.
+
+    This implementation uses the method first developed for Göransson et al.,
+    following the description in Birnbaum (2003). Note that, on large
+    sequences, it is much slower and more memory-intensive than the one
+    implemented in the `fast_birnbaum_simil` method.  While in most cases
+    the results are comparable, particularly after scaling/normalization,
+    the values are *not* identical.
+
+    See: Birnbaum, David J. (2003). "Computer-Assisted Analysis and
+        Study of the Structure of Mixed-Content Miscellanies".
+        Scripta & Scripta 1:15-64.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The similarity score between the two sequences. The higher
+        the similarity score, the more similar the two sequences are;
+        a similarity score of zero is the theoretical maximum difference
+        between two sequences.
+    """
+
+    # Make sure `seq_x` is shorter or equal in length to `seq_y`
+    if len(seq_x) < len(seq_y):
+        seq_x, seq_y = seq_y, seq_x
+
+    # Cast to lists to make sure it works with arbitrary hashable elements
+    seq_x = [element for element in seq_x]
+    seq_y = [element for element in seq_y]
+
+    similarity: int = 0
+    subseq_len = len(seq_y)
+    while subseq_len:
+        # iterate though all subsequences of N
+        for Nsubseq in _nwise(seq_y, subseq_len):
+            # test whether each of these subsequences matches a section of M
+            # (go to each instance of Nsubseq[0] in M and test the
+            # appropriate length slice from there)
+            for i in _indices(seq_x, Nsubseq[0]):
+                if seq_x[i : i + subseq_len] == Nsubseq:
+                    similarity += 1
+        subseq_len -= 1
+
+    if normal:
+        return float(similarity) / max(
+            [
+                birnbaum_simil(seq_x, seq_x),
+                birnbaum_simil(seq_y, seq_y),
+            ]
+        )
+
+    return float(similarity)
+
+
+def fast_birnbaum_simil(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Birnbaum similarity score with the fast method.
+
+    This implementation uses the experimental method we developed following
+    the description in Birnbaum (2003), which is much faster and less
+    memory-intensive than the one implemented in the `birnbaum_simil()`
+    function. While in most cases the results are comparable, particularly
+    after scaling/normalization, and while the ones
+    provided by this method might be considered more adequate due to their
+    handling of duplicate information, the values are *not* identical.
+
+    See: Birnbaum, David J. (2003). "Computer-Assisted Analysis and
+        Study of the Structure of Mixed-Content Miscellanies".
+        Scripta & Scripta 1:15-64.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The similarity score between the two sequences. The higher
+        the similarity score, the more similar the two sequences are;
+        a similarity score of zero is the theoretical maximum difference
+        between two sequences.
+    """
+
+    # If the sequences are equal, we can just compute the score from length
+    if seq_x == seq_y:
+        length = len(seq_x)
+        return (length * (length + 1)) // 2
+
+    # Make sure `seq_x` is shorter or equal in length to `seq_y`
+    if len(seq_x) < len(seq_y):
+        seq_x, seq_y = seq_y, seq_x
+
+    # Get the opcodes from `SequenceMatcher` and drop the last one (a dummy)
+    sm = difflib.SequenceMatcher(None, seq_x, seq_y)
+    blocks = sm.get_matching_blocks()[:-1]
+
+    # Sum sizes and return
+    sizes = [match.size for match in blocks]
+    similarity = sum([(v * (v + 1)) // 2 for v in sizes])
+
+    if normal:
+        return float(similarity) / max(
+            [
+                fast_birnbaum_simil(seq_x, seq_x),
+                fast_birnbaum_simil(seq_y, seq_y),
+            ]
+        )
+
+    return float(similarity)
+
+
+def fragile_ends_simil(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the "fragile ends" similarity between two sequences.
+
+    The "fragile ends" similarity is defined as equal to the Levenshtein one,
+    but with deletions in the initial or final 10% of the positions being
+    discounted.
+
+    This function will use the standard Wagner-Fischer algorithm with the
+    default costs provided by the internal `_levdamerau_costs()`
+    function. This similarity measure is not used directly in the paper and was a
+    proof-of-concept while working toward the "Stemmatological distance".
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The computed "fragile ends" similarity.
+    """
+
+    d = _fragile_ends_initial_matrix(seq_x, seq_y)
+    similarity = _wagner_fischer(seq_x, seq_y, _fragile_ends_costs, d)
+
+    if normal:
+        seq_x_list = [elem for elem in seq_x]
+        seq_y_list = [elem for elem in seq_y]
+
+        return similarity / max(
+            [
+                similarity,
+                fragile_ends_simil(seq_x_list + seq_y_list, seq_x_list),
+                fragile_ends_simil(seq_x_list + seq_y_list, seq_y_list),
+            ]
+        )
+
+    return similarity
+
+
+# TODO: frag type and description
+# TODO: should `frag` be in range 0..1, as a float?
+def stemmatological_simil(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    frag_start: float = 10.0,
+    frag_end: float = 10.0,
+    max_del_len: int = 5,
+    normal: bool = False,
+) -> float:
+    """
+    Compute the "stemmatological" similarity between two sequences.
+
+    This function will use the standard Wagner-Fischer algorithm with the
+    default costs provided by the internal `_levdamerau_costs()`
+    function. This similarity measure is essentially a combination of the
+    "fragile ends" and "bulk delete" methods, with the first one
+    generalised a little to allow specifying the size of both fragile regions.
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param max_del_len: The maximum length of deletion block.
+    @param frag_start:
+    @param frag_end:
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The computed "stemmatological" similarity.
+    """
+
+    d = _stemmatological_initial_matrix(seq_x, seq_y, max_del_len, frag_start, frag_end)
+    _stemmatology_costs = _stemmatological_costs_factory(
+        max_del_len, frag_start, frag_end
+    )
+
+    dist = _wagner_fischer(seq_x, seq_y, _stemmatology_costs, d)
+
+    if normal:
+        return dist / max([len(seq_x), len(seq_y)])
+
+    return dist
+
+
+# Methods for sequence distance
+# -----------------------------
+#
+# Methods for "similarity" are those that guarantee all of the three
+# distance properties of positivity, symmetry, identity-discerning, and
+# triangle inequality
+
+
+def levenshtein_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Levenshtein distance between two sequences.
+
+    This function will use the standard Wagner-Fischer algorithm with the
+    default costs provided by the internal `_levenshtein_costs()`
+    function.
+
+    See: https://en.wikipedia.org/wiki/Levenshtein_distance
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The computed Levenshtein distance.
+    """
+
+    dist = _wagner_fischer(seq_x, seq_y, _levenshtein_costs)
+
+    if normal:
+        return dist / max([len(seq_x), len(seq_y)])
+
+    return dist
+
+
+def levdamerau_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Damerau-Levenshtein distance between two sequences.
+
+    This function will use the standard Wagner-Fischer algorithm with the
+    default costs provided by the internal `_levdamerau_costs()`
+    function.
+
+    See: https://en.wikipedia.org/wiki/Damerau%E2%80%93Levenshtein_distance
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The computed Levenshtein distance.
+    """
+
+    dist = _wagner_fischer(seq_x, seq_y, _levdamerau_costs)
+
+    if normal:
+        return dist / max([len(seq_x), len(seq_y)])
+
+    return dist
+
+
+# TODO: compute `max_del_len`, if not provided, from sequence length?
+def bulk_delete_dist(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    max_del_len: int = 5,
+    normal: bool = False,
+) -> float:
+    """
+    Compute the "bulk delete" distance between two sequences.
+
+    This function will use the standard Wagner-Fischer algorithm with the
+    default costs provided by the internal `_levdamerau_costs()`
+    function. This distance measure is not used directly in the paper and was a
+    proof-of-concept while working toward the "Stemmatological distance".
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param max_del_len: The maximum length of deletion block.
+    @param normal: Whether to normalize the similarity score in range
+        [0..1] using sequence lengths.
+    @return: The computed "bulk delete" distance.
+    """
+
+    d = _bulk_delete_initial_matrix(seq_x, seq_y, max_del_len)
+    _bulk_delete_costs = _bulk_delete_costs_factory(max_del_len)
+    dist = _wagner_fischer(seq_x, seq_y, _bulk_delete_costs, d)
+
+    if normal:
+        return dist / max([len(seq_x), len(seq_y)])
+
+    return dist
+
+
+# TODO: include `external` parameter?
+def jaro_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes the Jaro distance between two sequences.
+
+    This function returns the value from the implementation provided by
+    the `textdistance` library.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    See: https://en.wikipedia.org/wiki/Jaro%E2%80%93Winkler_distance
+
+    @param seq_x: The first sequence of elements to be compared.
+    @param seq_y: The second sequence of elements to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Jaro distance between the two sequences.
+    """
+
+    dist = textdistance.JaroWinkler(winklerize=False, external=False)(seq_x, seq_y)
+
+    if normal:
+        logging.warning(
+            "Jaro distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - dist
+
+
+# TODO: include `external` parameter?
+def jaro_winkler_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes the Jaro-Winkler distance between two sequences.
+
+    This function returns the value from the implementation provided by
+    the `textdistance` library.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    See: https://en.wikipedia.org/wiki/Jaro%E2%80%93Winkler_distance
+
+    @param seq_x: The first sequence of elements to be compared.
+    @param seq_y: The second sequence of elements to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Jaro-Winkler distance between the two sequences.
+    """
+
+    dist = textdistance.JaroWinkler(winklerize=True, external=False)(seq_x, seq_y)
+
+    if normal:
+        logging.warning(
+            "Jaro-Winkler distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - dist
+
+
+def mmcwpa_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes an MMCWPA distance between two sequences.
+
+    MMCWPA is the Modifier Moving Contracting Window Pattern Algorithm, modified by
+    Tiago Tresoldi from a method published by Yang et al. (2001). In order to simplify
+    the logic, the function uses the auxiliary internal function `_mmcwpa()`.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    See: Yang, Q. X.; Yuan, Sung S.; Chun, Lu; Zhao, Li; Peng Sun. "Faster Algorithm of String
+    Comparison", eprint arXiv:cs/0112022, December 2001.
+
+    @param seq_x: The first sequence of elements to be compared.
+    @param seq_y: The second sequence of elements to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The MMCWPA distance between the two sequences.
+    """
+
+    # Cache original lengths, before any modifications
+    len_a, len_b = len(seq_x), len(seq_y)
+
+    # Initialize the `f_x` and `f_y` vectors with a single element each,
+    # the corresponding sequence; we also initialize the `ssnc` to zero.
+    f_x, f_y = [seq_x], [seq_y]
+    ssnc: float = 0.0
+    while f_x and f_y:
+        f_x, f_y, ssnc = _mmcwpa(f_x, f_y, ssnc)
+
+    if normal:
+        logging.warning(
+            "MMCWPA distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - ((ssnc / ((len_a + len_b) ** 2.0)) ** 0.5)
+
+
+def birnbaum_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Birnbaum similarity distance with the original method.
+
+    This implementation uses the original method we developed following
+    the description in Birnbaum (2003). See comments for `birnbaum_simil()`.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    See: Birnbaum, David J. (2003). "Computer-Assisted Analysis and
+        Study of the Structure of Mixed-Content Miscellanies".
+        Scripta & Scripta 1:15-64.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The distance between the two sequences. A distance of 0.0 indicates
+        identical sequences, and a distance of 1.0 indicates the maximum
+        theoretical distance between two sequences.
+    """
+
+    # Make sure `seq_y` is longer or equal in length to `seq_x`
+    if len(seq_y) > len(seq_x):
+        seq_x, seq_y = seq_y, seq_x
+
+    distance = 1 - (birnbaum_simil(seq_x, seq_y) / birnbaum_simil(seq_y, seq_y))
+
+    # Compute the distance and correct if necessary (as the distance returned
+    # by the method can be less than zero when `seq_x` contains a perfect
+    # `seq_y` plus doublets matching parts of `seq_y` as well).
+    distance = max(0.0, distance)
+
+    if normal:
+        logging.warning(
+            "Birnbaum distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return distance
+
+
+def fast_birnbaum_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Compute the Birnbaum similarity distance with the fast method.
+
+    This implementation uses the experimental method we developed following
+    the description in Birnbaum (2003), which is much faster and less
+    memory-intensive than the one implemented in the `birnbaum()`
+    method. While in most cases the results are comparable, and the ones
+    provided by this method might be considered more adequate due to their
+    handling of duplicate information, the values are *not* identical.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    See: Birnbaum, David J. (2003). "Computer-Assisted Analysis and
+        Study of the Structure of Mixed-Content Miscellanies".
+        Scripta & Scripta 1:15-64.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The distance between the two sequences. A distance of 0.0 indicates
+        identical sequences, and a distance of 1.0 indicates the maximum
+        theoretical distance between two sequences.
+    """
+
+    # Make sure `seq_y` is longer or equal in length to `seq_x`
+    if len(seq_y) > len(seq_x):
+        seq_x, seq_y = seq_y, seq_x
+
+    # While the `.fast_birnbaum_simil()` already has an optimization for
+    # computing the similarity score of an identity, we can do it here as
+    # we know that `seq_y` is equal to itself
+    denom = (len(seq_y) * (len(seq_y) + 1)) / 2
+    distance = 1.0 - (fast_birnbaum_simil(seq_x, seq_y) / denom)
+
+    if normal:
+        logging.warning(
+            "Fast Birnbaum distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return distance
+
+
+# Supporting internal functions
+# -----------------------------
+#
+# These are internal functions supporting various different computations,
+# mostly related to the computation of matrices for edit distance
+# methods. Some of these could (should?) be moved to a new module in
+# the future for better organization and a shorter file.
+
+
+def _levenshtein_costs(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    d: List[List[float]],
+    i: int,
+    j: int,
+) -> Tuple[float, float, float]:
+    """
+    Computes candidate costs for an entry of an edit distance matrix.
+
+    This internal function will compute the candidate costs for an entry
+    (i, j) in terms of the Levenshtein distance matrix (seq_a, seq_b),
+    each cost corresponding to one of the available edit operations.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param d: The "starting matrix" for the cost computation.
+    @param i: The index of `seq_x` to be considered.
+    @param j: The index of `seq_y` to be considered.
+    @return: A tuple with the costs for deletion, insertion, and substitution.
+    """
+    substitution_cost = 0 if seq_x[i - 1] == seq_y[j - 1] else 1
+    costs = (
+        d[i - 1][j] + 1,  # del
+        d[i][j - 1] + 1,  # ins
+        d[i - 1][j - 1] + substitution_cost,
+    )
+
+    return costs
+
+
+def _levdamerau_costs(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    d: List[List[float]],
+    i: int,
+    j: int,
+) -> Tuple[float, ...]:
+    """
+    Computes candidate costs for an entry of a Damerau-Levenshtein distance matrix.
+
+    This internal function will compute the candidate costs for an entry
+    (i, j) in terms of the Damerau-Levenshtein distance matrix (seq_a, seq_b),
+    each cost corresponding to one of the available edit operations.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param d: The "starting matrix" for the cost computation.
+    @param i: The index of `seq_x` to be considered.
+    @param j: The index of `seq_y` to be considered.
+    @return: A tuple with the costs for deletion, insertion, and substitution,
+        and the transposition cost if necessary.
+    """
+
+    # Start out as per Levenshtein
+    substitution_cost = 0 if seq_x[i - 1] == seq_y[j - 1] else 1
+    costs = [
+        d[i - 1][j] + 1,
+        d[i][j - 1] + 1,
+        d[i - 1][j - 1] + substitution_cost,
+    ]
+
+    # Add the transposition cost
+    if (
+        i > 1
+        and j > 1
+        and seq_x[i - 1] == seq_y[j - 2]
+        and seq_x[i - 2] == seq_y[j - 1]
+    ):
+        costs.append(d[i - 2][j - 2] + 1)
+
+    return tuple(costs)
+
+
+# TODO: return type and description
+def _fragile_ends_initial_matrix(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable]
+) -> List[List[float]]:
+    """
+    Compute a starting matrix for the Wagner-Fischer algorithm with "fragile ends" costs.
+
+    The generic starting matrix isn't applicable here, as the cost of deleting the entire
+    starting hashable element ("manuscript", in the original application) is less
+    than then sequence length.
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @return:
+    """
+
+    m = len(seq_x)
+    n = len(seq_y)
+    d: List[List[float]] = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
+    for i in range(1, m + 1):
+        if i <= round(0.1 * m) or i >= round(0.9 * m):
+            d[i][0] = d[i - 1][0] + 0.5
+        else:
+            d[i][0] = d[i - 1][0] + 1
+
+    # Insertions
+    for j in range(1, n + 1):
+        d[0][j] = j
+
+    return d
+
+
+def _fragile_ends_costs(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    d: List[List[float]],
+    i: int,
+    j: int,
+) -> Tuple[float, ...]:
+    """
+    Computes candidate costs for an entry of a "fragile ends" distance matrix.
+
+    This internal function will compute the candidate costs for an entry
+    (i, j) in terms of the Levenshtein distance matrix (seq_a, seq_b),
+    each cost corresponding to one of the available edit operations.
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param d: The "starting matrix" for the cost computation.
+    @param i: The index of `seq_x` to be considered.
+    @param j: The index of `seq_y` to be considered.
+    @return: A tuple with the costs for deletion, insertion, and substitution.
+    """
+
+    substitution_cost = 0 if seq_x[i - 1] == seq_y[j - 1] else 1
+    costs = [
+        d[i][j - 1] + 1,
+        d[i - 1][j - 1] + substitution_cost,
+    ]
+
+    # Discount deletion near ends
+    m = len(seq_x)
+    if i <= round(0.1 * m) or i >= round(0.9 * m):
+        costs.append(d[i - 1][j] + 0.5)
+    else:
+        costs.append(d[i - 1][j] + 1)
+
+    return tuple(costs)
+
+
+# TODO: return type and description
+# TODO: allow max_del_len as a float with percentage length?
+def _bulk_delete_initial_matrix(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], max_del_len: int
+):
+    """
+    Compute a starting matrix for the Wagner-Fischer algorithm with "bulk delete" costs.
+
+    The generic starting matrix isn't applicable here, as the cost of deleting the entire
+    starting hashable element ("manuscript", in the original application) is less
+    than then sequence length.
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param max_del_len: The maximum length of deletion block.
+    @return:
+    """
+
+    m = len(seq_x)
+    n = len(seq_y)
+    d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
+
+    # Bulk deletions
+    for i in range(1, m + 1):
+        bulk_deletions = int(i / max_del_len)
+        d[i][0] = bulk_deletions
+        remainder = i - bulk_deletions * max_del_len
+        if remainder:
+            d[i][0] += 1
+
+    # Insertions
+    for j in range(1, n + 1):
+        d[0][j] = j
+
+    return d
+
+
+# TODO: return description
+def _bulk_delete_costs_factory(max_del_len: int = 5) -> Callable:
+    """
+    Define and return a function for computing candidate costs for a "bulk delete" distance matrix.
+
+    @param max_del_len: The maximum length of deletion block.
+    @return:
+    """
+
+    def _bulk_delete_costs(
+        seq_x: Sequence[Hashable],
+        seq_y: Sequence[Hashable],
+        d: List[List[float]],
+        i: int,
+        j: int,
+    ) -> Tuple[float, ...]:
+        """
+        Computes candidate costs for an entry of a "bulk delete" distance matrix.
+
+        This internal function will compute the candidate costs for an entry
+        (i, j) in terms of the Levenshtein distance matrix (seq_a, seq_b),
+        each cost corresponding to one of the available edit operations.
+
+        @param seq_x: The first sequence to be compared.
+        @param seq_y: The second sequence to be compared.
+        @param d: The "starting matrix" for the cost computation.
+        @param i: The index of `seq_x` to be considered.
+        @param j: The index of `seq_y` to be considered.
+        @return: A tuple with the costs for deletion, insertion, and substitution.
+        """
+
+        substitution_cost = 0 if seq_x[i - 1] == seq_y[j - 1] else 1
+        costs = [
+            d[i][j - 1] + 1,  # ins
+            d[i - 1][j - 1] + substitution_cost,
+        ]
+        for n in range(1, min(max_del_len + 1, i)):
+            # Delete consecutive block of n
+            costs.append(d[i - n][j] + 1)
+
+        return tuple(costs)
+
+    return _bulk_delete_costs
+
+
+# TODO: return type and description
+# TODO: frag type and description
+def _stemmatological_initial_matrix(
+    seq_x: Sequence[Hashable],
+    seq_y: Sequence[Hashable],
+    max_del_len: int = 5,
+    frag_start: float = 10.0,
+    frag_end: float = 10.0,
+):
+    """
+    Compute a starting matrix for the Wagner-Fischer algorithm with "stemmatological" costs.
+
+    The generic starting matrix isn't applicable here, as the cost of deleting the entire
+    starting hashable element ("manuscript", in the original application) is less
+    than then sequence length.
+
+    See: Göransson, Elisabet; Maurits, Luke; Dahlman, Britt; Sarkisian, Karine Å.;
+        Rubenson, Samuel; Dunn, Michael. "Improved distance measures for 'mixed-content
+        miscellania' (in prep.).
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param max_del_len: The maximum length of deletion block.
+    @param frag_start:
+    @param frag_end:
+    @return:
+    """
+    m = len(seq_x)
+    n = len(seq_y)
+    d = [[0 for i in range(0, n + 1)] for j in range(0, m + 1)]
+    lower = round(m * frag_start / 100.0)
+    upper = round(m * (100 - frag_end) / 100.0)
+    for i in range(1, m + 1):
+        if i <= lower or i >= upper:
+            d[i][0]: float = d[i - min(i, max_del_len)][0] + 0.5
+        else:
+            d[i][0]: float = d[i - min(i, max_del_len)][0] + 1.0
+
+    # Insertions
+    for j in range(1, n + 1):
+        d[0][j] = j
+
+    return d
+
+
+# TODO: frag type and description
+# TODO: return description
+def _stemmatological_costs_factory(
+    max_del_len: int = 5, frag_start: float = 10.0, frag_end: float = 10.0
+) -> Callable:
+    """
+    Define and return a function for computing candidate costs for a "stemmatological" distance matrix.
+
+    @param max_del_len: The maximum length of deletion block.
+    @param frag_start:
+    @param frag_end:
+    @return:
+    """
+
+    def _stemmatological_costs(
+        seq_x: Sequence[Hashable],
+        seq_y: Sequence[Hashable],
+        d: List[List[float]],
+        i: int,
+        j: int,
+    ):
+        """
+        Computes candidate costs for an entry of a "stemmatological" distance matrix.
+
+        This internal function will compute the candidate costs for an entry
+        (i, j) in terms of the Levenshtein distance matrix (seq_a, seq_b),
+        each cost corresponding to one of the available edit operations.
+
+        @param seq_x: The first sequence to be compared.
+        @param seq_y: The second sequence to be compared.
+        @param d: The "starting matrix" for the cost computation.
+        @param i: The index of `seq_x` to be considered.
+        @param j: The index of `seq_y` to be considered.
+        @return: A tuple with the costs for deletion, insertion, and substitution.
+        """
+        substitution_cost = 0 if seq_x[i - 1] == seq_y[j - 1] else 1
+        costs = [
+            d[i][j - 1] + 1,
+            d[i - 1][j - 1] + substitution_cost,
+        ]
+        m = len(seq_x)
+        lower = round(m * frag_start / 100.0)
+        upper = round(m * (100 - frag_end) / 100.0)
+
+        # Delete consecutive block of n
+        for n in range(1, min(max_del_len, i)):
+            # Discount bulk deletion near ends
+            if i <= lower or i >= upper:
+                costs.append(d[i - n][j] + 0.5)
+            else:
+                costs.append(d[i - n][j] + 1)
+
+        return costs
+
+    return _stemmatological_costs
+
+
+def _mmcwpa(
+    seq_x: List[Sequence[Hashable]], seq_y: List[Sequence[Hashable]], ssnc: float
+) -> Tuple[List[Sequence[Hashable]], List[Sequence[Hashable]], float]:
+    """
+    Internal function for MMCWPA implementation.
+
+    In this implementation of the Modified Moving Contracting Window Pattern Algorithm
+    (MMCWPA) to calculate sequence similarity, we return a list of non-overlapping,
+    non-contiguous fields Fx, a list of non-overlapping, non-contiguous fields Fy, and
+    the SSNC value (the Sum of the Square of the Number of the same characters). This
+    function separates the core method of the implementation and makes recursive calls
+    easier.
+
+    @param seq_x: A list of sub-sequences, related to the first sequence.
+    @param seq_y: A list of sub-sequences, related to the second sequence.
+    @param ssnc: The previous SSNC value.
+    @return: A tuple whose first element is a list of remaining sub-sequences from the
+             first sequence, the second element is a list of remaining sub-sequences
+             from the second sequence, and the third element is the updated SSNC.
+    """
+
+    # the boolean value indicating if a total or partial
+    # match was found between subfields Fx and Fy; when
+    # a match is found, the variable is used to cascade
+    # out of the loops of the function
+    match = False
+
+    # the variables where to store the new collections of
+    # subfields, if any match is found; if these values
+    # are not changed and the empty lists are returned,
+    # stringcomp() will break the loop of comparison,
+    # calculate the similarity ratio and return its value
+    new_f_x, new_f_y = [], []
+
+    # search patterns in all subfields of Fx; the index of
+    # the subfield in the list is used for upgrading the
+    # list, if a pattern is a found
+    for idx_x, sf_x in enumerate(seq_x):
+        # 'length' stores the length of the sliding window,
+        # from full length to a single character
+        for length in range(len(sf_x), 0, -1):
+            # 'i' stores the starting index of the sliding
+            # window in Fx
+            for i in range(len(sf_x) - length + 1):
+                # extract the pattern for matching
+                pattern = sf_x[i : i + length]
+
+                # look for the pattern in Fy
+                for idx_y, sf_y in enumerate(seq_y):
+                    # 'j' stores the starting index in Fy; the
+                    # Python find() function returns -1 if there
+                    # is no match
+                    j = sequence_find(sf_y, pattern)
+                    if j is not None and not match:
+                        # the pattern was found; set 'new_fx' and
+                        # 'new_fy' to version of 'fx' and 'fy' with
+                        # the patterns removed, update the SSNC and
+                        # set 'match' as True, in order to cascade
+                        # out of the loops
+                        tmp_x = [sf_x[:i], sf_x[i + length :]]
+                        tmp_y = [sf_y[:j], sf_y[j + length :]]
+                        new_f_x = seq_x[:idx_x] + tmp_x + seq_x[idx_x + 1 :]
+                        new_f_y = seq_y[:idx_y] + tmp_y + seq_y[idx_y + 1 :]
+
+                        ssnc += (2 * length) ** 2
+
+                        match = True
+
+                        break
+
+                # if a match was found, end the sliding window
+                if match:
+                    break
+
+            # if a match was found, end Fx subfield enumeration
+            if match:
+                break
+
+        # remove any empty subfields due to pattern removal
+        new_f_x = [sf for sf in new_f_x if sf]
+        new_f_y = [sf for sf in new_f_y if sf]
+
+        return new_f_x, new_f_y, ssnc
diff --git a/src/seqsim/ngrams.py b/src/seqsim/ngrams.py
index d2b07eb..2d8aac2 100644
--- a/src/seqsim/ngrams.py
+++ b/src/seqsim/ngrams.py
@@ -7,7 +7,7 @@
 
 # Import Python standard libraries
 from itertools import chain
-from typing import Sequence
+from typing import Hashable, Optional, Sequence
 
 _PAD_SYMBOL = "$$$"
 
@@ -17,8 +17,8 @@
 # primitive as fast as possible. This is intentionally not defaulting to any
 # value for the order, so that users won't confuse a given order to all
 # orders up to and including the given one.
-# TODO: typing for pad_symbol and return
-def ngrams_iter(sequence: Sequence, order: int, pad_symbol="$$$"):
+# TODO: typing for return
+def ngrams_iter(sequence: Sequence, order: int, pad_symbol: Optional[Hashable] = "$$$"):
     """
     Build an iterator for collecting all ngrams of a given order.
     The sequence can optionally be padded with boundary symbols which are
@@ -51,7 +51,8 @@ def ngrams_iter(sequence: Sequence, order: int, pad_symbol="$$$"):
     yield from zip(*[seq[i:] for i in range(order)])
 
 
-# TODO: rename to `collect` as above
+# TODO: rename to `collect` as in the comments above
+# TODO: typing and doc string (also changing example)
 def get_all_ngrams_by_order(sequence, orders=None, pad_symbol=_PAD_SYMBOL):
     """
     Build an iterator for collecting all ngrams of a given set of orders.
diff --git a/src/seqsim/sequence.py b/src/seqsim/sequence.py
new file mode 100644
index 0000000..83858e7
--- /dev/null
+++ b/src/seqsim/sequence.py
@@ -0,0 +1,46 @@
+"""
+Module implementing various methods for similarity and distance from sequence methods.
+
+Most of the methods are commonly used in string comparison, such as
+RatcliffObershelp, but in this module we need to make sure we can operate on arbitrary
+iterable data structures.
+"""
+
+# Import Python standard libraries
+from typing import Hashable, Sequence
+import logging
+
+# Import 3rd-party libraries
+import textdistance
+
+# TODO: add a pure ratcliff_obershelp similarity?
+# TODO: multiple sequences?
+def ratcliff_obershelp(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a distance between two sequences based on the Ratcliff-Obershelp similarity.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Ratcliff-Obershelp distance is already
+    in range [0..1].
+
+    See: https://en.wikipedia.org/wiki/Gestalt_Pattern_Matching
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Ratcliff-Obershelp distance between the two sequences.
+    """
+
+    # As the method uses .find, we need strings
+    # TODO: rewrite proper solution generalizing for all sequences
+    seq_x = "".join([str(elem) for elem in seq_x])
+    seq_y = "".join([str(elem) for elem in seq_y])
+
+    if normal:
+        logging.warning(
+            "Sørensen–Dice distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - textdistance.RatcliffObershelp(external=False)(seq_x, seq_y)
diff --git a/src/seqsim/similarity.py b/src/seqsim/similarity.py
deleted file mode 100644
index b1b46ca..0000000
--- a/src/seqsim/similarity.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import difflib
-
-
-def birnbaum(M, N):
-    sm = difflib.SequenceMatcher(None, M, N)
-    blocks = sm.get_matching_blocks()
-    sizes = [1 if match.size == 0 else match.size for match in blocks]
-    formula = [(v * (v + 1)) // 2 for v in sizes]
-
-    return sum(formula)
diff --git a/src/seqsim/token.py b/src/seqsim/token.py
new file mode 100644
index 0000000..1ab6e78
--- /dev/null
+++ b/src/seqsim/token.py
@@ -0,0 +1,118 @@
+"""
+Module implementing various methods for similarity and distance from token methods.
+
+Most of the methods are commonly used in string comparison, such as Jaccard
+index, but in this module we need to make sure we can operate on arbitrary
+iterable data structures.
+"""
+
+# Import Python standard libraries
+from typing import Hashable, Sequence
+import logging
+
+# Import 3rd-party libraries
+import textdistance
+
+# Import local modules
+from .common import collect_subseqs
+
+# TODO: have a jaccard similarity?
+def jaccard_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a Jaccard distance between two sequences.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Jaccard distance between the two sequences.
+    """
+
+    intersection = len(set(seq_x).intersection(seq_y))
+    union = len(seq_x) + len(seq_y) - intersection
+
+    if normal:
+        logging.warning(
+            "Jaccard distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - (float(intersection) / union)
+
+
+# TODO: rename appropriately with other methods, consider ngram usage
+# TODO: have a subseq_jaccard similarity?
+def subseq_jaccard_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a Jaccard distance between two sequences using sub-sequence occurrence.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Jaccard distance is already
+    in range [0..1].
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Subseq-Jaccard distance between the two sequences.
+    """
+
+    subseqs1 = collect_subseqs(seq_x)
+    subseqs2 = collect_subseqs(seq_y)
+
+    # From the longest subseq, which is the length of the longest sequence,
+    # collect all subsequences of that given length in both sets, compute the
+    # Jaccard index, correct it by length of the subsequence (so that longer ones
+    # will score higher) and update the internal results before returning.
+    # TODO: collect beforehand in a list, so we don´t repeat the comprehension
+    # TODO: convert to tuple beforehand as well?
+    jaccard_scores = []
+    max_length = max([len(seq_x), len(seq_y)])
+    for length in range(max_length, 0, -1):
+        l_subseq1 = [tuple(ss) for ss in subseqs1 if len(ss) == length]
+        l_subseq2 = [tuple(ss) for ss in subseqs2 if len(ss) == length]
+
+        intersection = len(set(l_subseq1).intersection(l_subseq2))
+        union = len(l_subseq1) + len(l_subseq2) - intersection
+        jaccard_scores.append((float(intersection) / union) * length)
+
+    # Compute the denominator, as the highest possible value
+    den = (max_length * (max_length + 1)) / 2.0
+
+    if normal:
+        logging.warning(
+            "Subseq-Jaccard distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return (1.0 - (sum(jaccard_scores) / den)) ** max_length
+
+
+def sorensen_dist(
+    seq_x: Sequence[Hashable], seq_y: Sequence[Hashable], normal: bool = False
+) -> float:
+    """
+    Computes a distance between two sequences based on the Sørensen–Dice coefficient.
+
+    The function accepts the `normal` parameter to have calls equivalent to those
+    of other methods, but it is redundant as the Sørensen–Dice distance is already
+    in range [0..1].
+
+    See: https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient
+
+    @param seq_x: The first sequence to be compared.
+    @param seq_y: The second sequence to be compared.
+    @param normal: Dummy parameter, see comment above.
+    @return: The Sørensen–Dice distance between the two sequences.
+    """
+
+    if normal:
+        logging.warning(
+            "Sørensen–Dice distance is always in [0..1] range, no need for `normal` parameter."
+        )
+
+    return 1.0 - textdistance.Sorensen(external=False)(seq_x, seq_y)


Test patch:
diff --git a/tests/test_compression.py b/tests/test_compression.py
new file mode 100644
index 0000000..57d850c
--- /dev/null
+++ b/tests/test_compression.py
@@ -0,0 +1,48 @@
+"""
+test_compression
+================
+
+Tests for the `compression` module of the `seqsim` package.
+"""
+
+# Import Python standard libraries
+import pytest
+
+# Import the library being tested
+from seqsim import compression
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 1.25, 0.0],
+        [(1, 2, 3), [1, 2, 3], 1.25, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 1.0, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.555555, 1e-6],
+    ],
+)
+def test_arith_ncd(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert compression.arith_ncd(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert compression.arith_ncd(seq_x, seq_y) == compression.arith_ncd(seq_y, seq_x)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.101340, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.114430, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 0.407464, 1e-6],
+    ],
+)
+def test_entropy_ncd(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert compression.entropy_ncd(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert compression.entropy_ncd(seq_x, seq_y) == compression.entropy_ncd(
+        seq_y, seq_x
+    )
diff --git a/tests/test_distance.py b/tests/test_distance.py
deleted file mode 100644
index 5ee91b7..0000000
--- a/tests/test_distance.py
+++ /dev/null
@@ -1,133 +0,0 @@
-"""
-test_distance
-=============
-
-Tests for the `distance` module of the `seqsim` package.
-"""
-
-# Import Python standard libraries
-import pytest
-
-# Import the library being tested
-from seqsim import distance
-
-test1_a = "kitten"
-test1_b = "sitting"
-test2_a = (1, 2, 3)
-test2_b = [1, 2, 3]
-test3_a = (1, 2, 3, 4, 5)
-test3_b = (1, 2, 4, 3, 6, 7)
-test4_a = (1, 2, 3)
-test4_b = ["a", "b", "c", "d"]
-
-
-def test_birnbaum_distance():
-    assert distance.birnbaum(test1_a, test1_b) == pytest.approx(
-        0.636363, abs=1e-6
-    )
-    assert distance.birnbaum(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.birnbaum(test3_a, test3_b) == pytest.approx(0.6875)
-    assert distance.birnbaum(test4_a, test4_b) == pytest.approx(
-        0.857142, abs=1e-6
-    )
-
-
-def test_levenshtein_distance():
-    assert distance.levenshtein(test1_a, test1_b) == pytest.approx(3.0)
-    assert distance.levenshtein(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.levenshtein(test3_a, test3_b) == pytest.approx(3.0)
-    assert distance.levenshtein(test4_a, test4_b) == pytest.approx(4.0)
-
-
-def test_norm_levenshtein_distance():
-    assert distance.norm_levenshtein(test1_a, test1_b) == pytest.approx(
-        0.428571, abs=1e-6
-    )
-    assert distance.norm_levenshtein(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.norm_levenshtein(test3_a, test3_b) == pytest.approx(0.5)
-    assert distance.norm_levenshtein(test4_a, test4_b) == pytest.approx(1.0)
-
-
-# TODO: add tests where levdamerau is different from levenshtein
-def test_levdamerau_distance():
-    assert distance.levdamerau(test1_a, test1_b) == pytest.approx(3.0)
-    assert distance.levdamerau(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.levdamerau(test3_a, test3_b) == pytest.approx(3.0)
-    assert distance.levdamerau(test4_a, test4_b) == pytest.approx(4.0)
-
-
-# TODO: add tests where fragile ends is different from levenshtein
-def test_fragile_ends_distance():
-    assert distance.fragile_ends(test1_a, test1_b) == pytest.approx(3.0)
-    assert distance.fragile_ends(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.fragile_ends(test3_a, test3_b) == pytest.approx(3.0)
-    assert distance.fragile_ends(test4_a, test4_b) == pytest.approx(4.0)
-
-
-# TODO: add tests where bulk delete is different from levenshtein
-def test_bulk_delete_distance():
-    assert distance.bulk_delete(test1_a, test1_b) == pytest.approx(3.0)
-    assert distance.bulk_delete(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.bulk_delete(test3_a, test3_b) == pytest.approx(3.0)
-    assert distance.bulk_delete(test4_a, test4_b) == pytest.approx(4.0)
-
-
-# TODO: add tests where stemmatology is different from levenshtein
-def test_stemmatology_distance():
-    assert distance.stemmatological(test1_a, test1_b) == pytest.approx(3.0)
-    assert distance.stemmatological(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.stemmatological(test3_a, test3_b) == pytest.approx(3.0)
-    assert distance.stemmatological(test4_a, test4_b) == pytest.approx(4.0)
-
-
-# TODO: add tests where stemmatology (norm) is different from levenshtein (norm)
-def test_norm_stemmatology_distance():
-    assert distance.norm_stemmatological(test1_a, test1_b) == pytest.approx(
-        0.428571, abs=1e-6
-    )
-    assert distance.norm_stemmatological(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.norm_stemmatological(test3_a, test3_b) == pytest.approx(0.5)
-    assert distance.norm_stemmatological(test4_a, test4_b) == pytest.approx(1.0)
-
-
-# TODO: add tests where stemmatology (norm/2030) is different from levenshtein (norm)
-def test_norm_stemmatology_2030_distance():
-    assert distance.norm_stemmatological_2030(
-        test1_a, test1_b
-    ) == pytest.approx(0.428571, abs=1e-6)
-    assert distance.norm_stemmatological_2030(
-        test2_a, test2_b
-    ) == pytest.approx(0.0)
-    assert distance.norm_stemmatological_2030(
-        test3_a, test3_b
-    ) == pytest.approx(0.5)
-    assert distance.norm_stemmatological_2030(
-        test4_a, test4_b
-    ) == pytest.approx(1.0)
-
-
-def test_jaccard_distance():
-    assert distance.jaccard(test1_a, test1_b) == pytest.approx(0.7)
-    assert distance.jaccard(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.jaccard(test3_a, test3_b) == pytest.approx(
-        0.428571, rel=1e-5
-    )
-    assert distance.jaccard(test4_a, test4_b) == pytest.approx(1.0)
-
-
-def test_subseq_jaccard_distance():
-    assert distance.subseq_jaccard(test1_a, test1_b) == pytest.approx(0.751556)
-    assert distance.subseq_jaccard(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.subseq_jaccard(test3_a, test3_b) == pytest.approx(0.787094)
-    assert distance.subseq_jaccard(test4_a, test4_b) == pytest.approx(1.0)
-
-
-def test_mmcwpa_distance():
-    assert distance.mmcwpa(test1_a, test1_b) == pytest.approx(
-        0.538461, rel=1e-5
-    )
-    assert distance.mmcwpa(test2_a, test2_b) == pytest.approx(0.0)
-    assert distance.mmcwpa(test3_a, test3_b) == pytest.approx(
-        0.554638, rel=1e-5
-    )
-    assert distance.mmcwpa(test4_a, test4_b) == pytest.approx(1.0)
diff --git a/tests/test_edit.py b/tests/test_edit.py
new file mode 100644
index 0000000..26e33b4
--- /dev/null
+++ b/tests/test_edit.py
@@ -0,0 +1,361 @@
+"""
+test_edit
+=========
+
+Tests for the `edit` module of the `seqsim` package.
+"""
+
+# TODO: add tests with empty sequences
+
+# Import Python standard libraries
+import pytest
+
+# Import the library being tested
+from seqsim import edit
+
+# TESTS FOR MEASURES OF SIMILARITY
+# --------------------------------
+#
+# Measures of similarity don't offer all the properties of distance measures (i.e.,
+# positivity, symmetry, identity-discerning, and triangle inequality), but may
+# offer some of these.
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 10.0, 0.3125, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 6.0, 1.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 5.0, 0.238095, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 0.0, 0.0, 0.0],
+    ],
+)
+def test_birnbaum_simil(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.birnbaum_simil(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value with normalization
+    assert edit.birnbaum_simil(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 7.0, 0.25, 0.0],
+        [(1, 2, 3), [1, 2, 3], 6.0, 1.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 4.0, 0.190476, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 0.0, 0.0, 0.0],
+    ],
+)
+def test_fast_birnbaum_simil(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.fast_birnbaum_simil(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value with normalization
+    assert edit.fast_birnbaum_simil(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+
+# TODO: add tests where fragile ends is different from levenshtein
+# TODO: move fragile_end to similarity as it is not (necessarily) symmetric
+# TODO: add a symmetric version with the mean?
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 3.0, 0.5, 0.0],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 3.0, 0.6, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 4.0, 1.0, 0.0],
+    ],
+)
+def test_fragile_ends_distance(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.fragile_ends_simil(seq_x, seq_y) == pytest.approx(expected)
+
+    # Test hard-coded expected value with normalization
+    assert edit.fragile_ends_simil(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.fragile_ends_simil(seq_x, seq_z)
+    dist_xy = edit.fragile_ends_simil(seq_x, seq_y)
+    dist_yz = edit.fragile_ends_simil(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+# TODO: add tests where stemmatology is different from levenshtein
+# TODO: move stemmatological to similarity as it is not (necessarily) symmetric
+# TODO: add a symmetric version with the mean?
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 3.0, 0.428571, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 3.0, 0.5, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 4.0, 1.0, 0.0],
+    ],
+)
+def test_stemmatology_distance(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.stemmatological_simil(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value without normalization
+    assert edit.stemmatological_simil(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.stemmatological_simil(seq_x, seq_z)
+    dist_xy = edit.stemmatological_simil(seq_x, seq_y)
+    dist_yz = edit.stemmatological_simil(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+# TODO: add tests where stemmatology (norm/2030) is different from levenshtein (norm)
+# TODO: move norm_stemmatological to similarity as it is not (necessarily) symmetric
+# TODO: add a symmetric version with the mean?
+# TODO: add non-normal test (if we keep this)
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.428571, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.5, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_stemmatology_2030_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.stemmatological_simil(
+        seq_x, seq_y, 20.0, 30.0, normal=True
+    ) == pytest.approx(expected, abs=tol)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.stemmatological_simil(seq_x, seq_z, 20.0, 30.0, normal=True)
+    dist_xy = edit.stemmatological_simil(seq_x, seq_y, 20.0, 30.0, normal=True)
+    dist_yz = edit.stemmatological_simil(seq_y, seq_z, 20.0, 30.0, normal=True)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+# TESTS FOR MEASURES OF DISTANCE
+# ------------------------------
+#
+# Measures of distance offer all the properties of distance measures (i.e.,
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 3.0, 0.428571, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 3.0, 0.5, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 4.0, 1.0, 0.0],
+    ],
+)
+def test_levenshtein_distance(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.levenshtein_dist(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value with normalization
+    assert edit.levenshtein_dist(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+    # Test symmetry
+    assert edit.levenshtein_dist(seq_x, seq_y) == edit.levenshtein_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.levenshtein_dist(seq_x, seq_z)
+    dist_xy = edit.levenshtein_dist(seq_x, seq_y)
+    dist_yz = edit.levenshtein_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+# TODO: add tests where levdamerau is different from levenshtein
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 3.0, 0.428571, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 3.0, 0.5, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 4.0, 1.0, 0.0],
+    ],
+)
+def test_levdamerau_distance(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.levdamerau_dist(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value with normalization
+    assert edit.levdamerau_dist(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+    # Test symmetry
+    assert edit.levdamerau_dist(seq_x, seq_y) == edit.levdamerau_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.levdamerau_dist(seq_x, seq_z)
+    dist_xy = edit.levdamerau_dist(seq_x, seq_y)
+    dist_yz = edit.levdamerau_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+# TODO: add tests where bulk delete is different from levenshtein
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,expected_norm,tol_norm",
+    [
+        ["kitten", "sitting", 3.0, 0.428571, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 3.0, 0.5, 0.0],
+        [(1, 2, 3), ["a", "b", "c", "d"], 4.0, 1.0, 0.0],
+    ],
+)
+def test_bulk_delete_distance(seq_x, seq_y, expected, expected_norm, tol_norm):
+    # Test hard-coded expected value without normalization
+    assert edit.bulk_delete_dist(seq_x, seq_y) == expected
+
+    # Test hard-coded expected value without normalization
+    assert edit.bulk_delete_dist(seq_x, seq_y, normal=True) == pytest.approx(
+        expected_norm, abs=tol_norm
+    )
+
+    # Test symmetry
+    assert edit.bulk_delete_dist(seq_x, seq_y) == edit.bulk_delete_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.bulk_delete_dist(seq_x, seq_z)
+    dist_xy = edit.bulk_delete_dist(seq_x, seq_y)
+    dist_yz = edit.bulk_delete_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.253968, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.261111, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_jaro_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.jaro_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert edit.jaro_dist(seq_x, seq_y) == edit.jaro_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.jaro_dist(seq_x, seq_z)
+    dist_xy = edit.jaro_dist(seq_x, seq_y)
+    dist_yz = edit.jaro_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.253968, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.208888, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_jarowinkler_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.jaro_winkler_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert edit.jaro_winkler_dist(seq_x, seq_y) == edit.jaro_winkler_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.jaro_winkler_dist(seq_x, seq_z)
+    dist_xy = edit.jaro_winkler_dist(seq_x, seq_y)
+    dist_yz = edit.jaro_winkler_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.538461, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.554638, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_mmcwpa_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.mmcwpa_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert edit.mmcwpa_dist(seq_x, seq_y) == edit.mmcwpa_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.mmcwpa_dist(seq_x, seq_z)
+    dist_xy = edit.mmcwpa_dist(seq_x, seq_y)
+    dist_yz = edit.mmcwpa_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.565217, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.666666, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_birnbaum_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.birnbaum_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert edit.birnbaum_dist(seq_x, seq_y) == edit.birnbaum_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.birnbaum_dist(seq_x, seq_z)
+    dist_xy = edit.birnbaum_dist(seq_x, seq_y)
+    dist_yz = edit.birnbaum_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.666666, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.733333, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_fast_birnbaum_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert edit.fast_birnbaum_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert edit.fast_birnbaum_dist(seq_x, seq_y) == edit.fast_birnbaum_dist(
+        seq_y, seq_x
+    )
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = edit.fast_birnbaum_dist(seq_x, seq_z)
+    dist_xy = edit.fast_birnbaum_dist(seq_x, seq_y)
+    dist_yz = edit.fast_birnbaum_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
diff --git a/tests/test_sequence.py b/tests/test_sequence.py
new file mode 100644
index 0000000..0df5611
--- /dev/null
+++ b/tests/test_sequence.py
@@ -0,0 +1,38 @@
+"""
+test_sequence
+=============
+
+Tests for the `sequence` module of the `seqsim` package.
+"""
+
+# Import Python standard libraries
+import pytest
+
+# Import the library being tested
+from seqsim import sequence
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting", 0.384615, 1e-6],
+        [(1, 2, 3), [1, 2, 3], 0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7), 0.454545, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"], 1.0, 0.0],
+    ],
+)
+def test_ratcliff_obershelp(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert sequence.ratcliff_obershelp(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert sequence.ratcliff_obershelp(seq_x, seq_y) == sequence.ratcliff_obershelp(
+        seq_y, seq_x
+    )
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = sequence.ratcliff_obershelp(seq_x, seq_z)
+    dist_xy = sequence.ratcliff_obershelp(seq_x, seq_y)
+    dist_yz = sequence.ratcliff_obershelp(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
diff --git a/tests/test_similarity.py b/tests/test_similarity.py
deleted file mode 100644
index 1600222..0000000
--- a/tests/test_similarity.py
+++ /dev/null
@@ -1,28 +0,0 @@
-"""
-test_similarity
-===============
-
-Tests for the `similarity` module of the `seqsim` package.
-"""
-
-# Import Python standard libraries
-import pytest
-
-# Import the library being tested
-from seqsim import similarity
-
-test1_a = "kitten"
-test1_b = "sitting"
-test2_a = (1, 2, 3)
-test2_b = [1, 2, 3]
-test3_a = (1, 2, 3, 4, 5)
-test3_b = (1, 2, 4, 3, 6, 7)
-test4_a = (1, 2, 3)
-test4_b = ["a", "b", "c", "d"]
-
-
-def test_birnbaum_similarity():
-    assert similarity.birnbaum(test1_a, test1_b) == 8
-    assert similarity.birnbaum(test2_a, test2_b) == 7
-    assert similarity.birnbaum(test3_a, test3_b) == 5
-    assert similarity.birnbaum(test4_a, test4_b) == 1
diff --git a/tests/test_token.py b/tests/test_token.py
new file mode 100644
index 0000000..352ffd7
--- /dev/null
+++ b/tests/test_token.py
@@ -0,0 +1,86 @@
+"""
+test_token
+==========
+
+Tests for the `token` module of the `seqsim` package.
+"""
+
+# Import Python standard libraries
+import pytest
+
+# Import the library being tested
+from seqsim import token
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting"] + [0.7, 0.0],
+        [(1, 2, 3), [1, 2, 3]] + [0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7)] + [0.428571, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"]] + [1.0, 0.0],
+    ],
+)
+def test_jaccard_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert token.jaccard_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert token.jaccard_dist(seq_x, seq_y) == token.jaccard_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = token.jaccard_dist(seq_x, seq_z)
+    dist_xy = token.jaccard_dist(seq_x, seq_y)
+    dist_yz = token.jaccard_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting"] + [0.751556, 1e-6],
+        [(1, 2, 3), [1, 2, 3]] + [0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7)] + [0.787094, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"]] + [1.0, 0.0],
+    ],
+)
+def test_subseq_jaccard_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert token.subseq_jaccard_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert token.subseq_jaccard_dist(seq_x, seq_y) == token.subseq_jaccard_dist(
+        seq_y, seq_x
+    )
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = token.subseq_jaccard_dist(seq_x, seq_z)
+    dist_xy = token.subseq_jaccard_dist(seq_x, seq_y)
+    dist_yz = token.subseq_jaccard_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)
+
+
+@pytest.mark.parametrize(
+    "seq_x,seq_y,expected,tol",
+    [
+        ["kitten", "sitting"] + [0.384615, 1e-6],
+        [(1, 2, 3), [1, 2, 3]] + [0.0, 0.0],
+        [(1, 2, 3, 4, 5), (1, 2, 4, 3, 6, 7)] + [0.272727, 1e-6],
+        [(1, 2, 3), ["a", "b", "c", "d"]] + [1.0, 0.0],
+    ],
+)
+def test_sorensen_distance(seq_x, seq_y, expected, tol):
+    # Test hard-coded expected value
+    assert token.sorensen_dist(seq_x, seq_y) == pytest.approx(expected, abs=tol)
+
+    # Test symmetry
+    assert token.sorensen_dist(seq_x, seq_y) == token.sorensen_dist(seq_y, seq_x)
+
+    # Test triangle-inequality
+    seq_z = [element for element in seq_x] + [element for element in seq_y]
+    dist_xz = token.sorensen_dist(seq_x, seq_z)
+    dist_xy = token.sorensen_dist(seq_x, seq_y)
+    dist_yz = token.sorensen_dist(seq_y, seq_z)
+    assert dist_xz <= (dist_xy + dist_yz)